Skip to content
Snippets Groups Projects
Commit 8fb2313e authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files 

- added Orientation() members for hexas, prisms and pyramids

- use det3x3() instead of ad-hoc code everywhere it makes sense
parent c6affe0c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/Element.cpp
+ 44
− 1
View file @ 8fb2313e
// $Id: Element.cpp,v 1.1 2004-05-25 04:10:05 geuzaine Exp $ // $Id: Element.cpp,v 1.2 2004-07-21 22:19:56 geuzaine Exp $
// //
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
// //
...@@ -21,6 +21,7 @@ ...@@ -21,6 +21,7 @@
#include "Gmsh.h" #include "Gmsh.h"
#include "Mesh.h" #include "Mesh.h"
#include "Numeric.h"
#include "Element.h" #include "Element.h"
int Element::TotalNumber = 0; int Element::TotalNumber = 0;
...@@ -107,6 +108,20 @@ Hexahedron::Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, ...@@ -107,6 +108,20 @@ Hexahedron::Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4,
VSUP = NULL; VSUP = NULL;
} }
double Hexahedron::Orientation(){
double mat[3][3];
mat[0][0] = V[1]->Pos.X - V[0]->Pos.X;
mat[0][1] = V[3]->Pos.X - V[0]->Pos.X;
mat[0][2] = V[4]->Pos.X - V[0]->Pos.X;
mat[1][0] = V[1]->Pos.Y - V[0]->Pos.Y;
mat[1][1] = V[3]->Pos.Y - V[0]->Pos.Y;
mat[1][2] = V[4]->Pos.Y - V[0]->Pos.Y;
mat[2][0] = V[1]->Pos.Z - V[0]->Pos.Z;
mat[2][1] = V[3]->Pos.Z - V[0]->Pos.Z;
mat[2][2] = V[4]->Pos.Z - V[0]->Pos.Z;
return det3x3(mat);
}
Hexahedron::~Hexahedron() Hexahedron::~Hexahedron()
{ {
if(VSUP) Free(VSUP); if(VSUP) Free(VSUP);
...@@ -160,6 +175,20 @@ Prism::Prism(Vertex *v1, Vertex *v2, Vertex *v3, ...@@ -160,6 +175,20 @@ Prism::Prism(Vertex *v1, Vertex *v2, Vertex *v3,
VSUP = NULL; VSUP = NULL;
} }
double Prism::Orientation(){
double mat[3][3];
mat[0][0] = V[1]->Pos.X - V[0]->Pos.X;
mat[0][1] = V[2]->Pos.X - V[0]->Pos.X;
mat[0][2] = V[3]->Pos.X - V[0]->Pos.X;
mat[1][0] = V[1]->Pos.Y - V[0]->Pos.Y;
mat[1][1] = V[2]->Pos.Y - V[0]->Pos.Y;
mat[1][2] = V[3]->Pos.Y - V[0]->Pos.Y;
mat[2][0] = V[1]->Pos.Z - V[0]->Pos.Z;
mat[2][1] = V[2]->Pos.Z - V[0]->Pos.Z;
mat[2][2] = V[3]->Pos.Z - V[0]->Pos.Z;
return det3x3(mat);
}
Prism::~Prism() Prism::~Prism()
{ {
if(VSUP) Free(VSUP); if(VSUP) Free(VSUP);
...@@ -200,6 +229,20 @@ Pyramid::Pyramid() ...@@ -200,6 +229,20 @@ Pyramid::Pyramid()
VSUP = NULL; VSUP = NULL;
} }
double Pyramid::Orientation(){
double mat[3][3];
mat[0][0] = V[1]->Pos.X - V[0]->Pos.X;
mat[0][1] = V[3]->Pos.X - V[0]->Pos.X;
mat[0][2] = V[4]->Pos.X - V[0]->Pos.X;
mat[1][0] = V[1]->Pos.Y - V[0]->Pos.Y;
mat[1][1] = V[3]->Pos.Y - V[0]->Pos.Y;
mat[1][2] = V[4]->Pos.Y - V[0]->Pos.Y;
mat[2][0] = V[1]->Pos.Z - V[0]->Pos.Z;
mat[2][1] = V[3]->Pos.Z - V[0]->Pos.Z;
mat[2][2] = V[4]->Pos.Z - V[0]->Pos.Z;
return det3x3(mat);
}
Pyramid::Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5) Pyramid::Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5)
{ {
iEnt = -1; iEnt = -1;
......
Mesh/Element.h
+ 3
− 0
View file @ 8fb2313e
...@@ -49,6 +49,7 @@ class Hexahedron : public Element{ ...@@ -49,6 +49,7 @@ class Hexahedron : public Element{
Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4,
Vertex *v5, Vertex *v6, Vertex *v7, Vertex *v8); Vertex *v5, Vertex *v6, Vertex *v7, Vertex *v8);
~Hexahedron(); ~Hexahedron();
double Orientation();
}; };
class Prism : public Element{ class Prism : public Element{
...@@ -58,6 +59,7 @@ class Prism : public Element{ ...@@ -58,6 +59,7 @@ class Prism : public Element{
Prism(Vertex *v1, Vertex *v2, Vertex *v3, Prism(Vertex *v1, Vertex *v2, Vertex *v3,
Vertex *v4, Vertex *v5, Vertex *v6); Vertex *v4, Vertex *v5, Vertex *v6);
~Prism(); ~Prism();
double Orientation();
}; };
class Pyramid : public Element{ class Pyramid : public Element{
...@@ -66,6 +68,7 @@ class Pyramid : public Element{ ...@@ -66,6 +68,7 @@ class Pyramid : public Element{
Pyramid(); Pyramid();
Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5); Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5);
~Pyramid(); ~Pyramid();
double Orientation();
}; };
// C interface // C interface
......
Mesh/Simplex.cpp
+ 2
− 4
View file @ 8fb2313e
// $Id: Simplex.cpp,v 1.33 2004-05-25 23:16:27 geuzaine Exp $ // $Id: Simplex.cpp,v 1.34 2004-07-21 22:19:56 geuzaine Exp $
// //
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
// //
...@@ -203,9 +203,7 @@ double Simplex::matsimpl(double mat[3][3]) ...@@ -203,9 +203,7 @@ double Simplex::matsimpl(double mat[3][3])
mat[2][0] = V[1]->Pos.Z - V[0]->Pos.Z; mat[2][0] = V[1]->Pos.Z - V[0]->Pos.Z;
mat[2][1] = V[2]->Pos.Z - V[0]->Pos.Z; mat[2][1] = V[2]->Pos.Z - V[0]->Pos.Z;
mat[2][2] = V[3]->Pos.Z - V[0]->Pos.Z; mat[2][2] = V[3]->Pos.Z - V[0]->Pos.Z;
return (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) - return det3x3(mat);
mat[1][0] * (mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2]) +
mat[2][0] * (mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]));
} }
double Simplex::rhoin() double Simplex::rhoin()
......
Numeric/Numeric.cpp
+ 11
− 18
View file @ 8fb2313e
// $Id: Numeric.cpp,v 1.16 2004-05-08 00:36:00 geuzaine Exp $ // $Id: Numeric.cpp,v 1.17 2004-07-21 22:19:56 geuzaine Exp $
// //
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
// //
...@@ -122,8 +122,7 @@ int sys2x2(double mat[2][2], double b[2], double res[2]) ...@@ -122,8 +122,7 @@ int sys2x2(double mat[2][2], double b[2], double res[2])
double det, ud, norm; double det, ud, norm;
int i; int i;
norm = norm = DSQR(mat[0][0]) + DSQR(mat[1][1]) + DSQR(mat[0][1]) + DSQR(mat[1][0]);
DSQR(mat[0][0]) + DSQR(mat[1][1]) + DSQR(mat[0][1]) + DSQR(mat[1][0]);
det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]; det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
...@@ -146,14 +145,19 @@ int sys2x2(double mat[2][2], double b[2], double res[2]) ...@@ -146,14 +145,19 @@ int sys2x2(double mat[2][2], double b[2], double res[2])
return (1); return (1);
} }
double det3x3(double mat[3][3])
{
return (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]));
}
int sys3x3(double mat[3][3], double b[3], double res[3], double *det) int sys3x3(double mat[3][3], double b[3], double res[3], double *det)
{ {
double ud; double ud;
int i; int i;
*det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) - *det = det3x3(mat);
mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
if(*det == 0.0) { if(*det == 0.0) {
res[0] = res[1] = res[2] = 0.0; res[0] = res[1] = res[2] = 0.0;
...@@ -203,21 +207,11 @@ int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det) ...@@ -203,21 +207,11 @@ int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
return out; return out;
} }
int det3x3(double mat[3][3], double *det)
{
*det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
return 1;
}
int inv3x3(double mat[3][3], double inv[3][3], double *det) int inv3x3(double mat[3][3], double inv[3][3], double *det)
{ {
double ud; double ud;
*det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) - *det = det3x3(mat);
mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
if(*det == 0.0) if(*det == 0.0)
return (0); return (0);
...@@ -234,7 +228,6 @@ int inv3x3(double mat[3][3], double inv[3][3], double *det) ...@@ -234,7 +228,6 @@ int inv3x3(double mat[3][3], double inv[3][3], double *det)
inv[2][1] = -ud * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]); inv[2][1] = -ud * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
inv[2][2] = ud * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]); inv[2][2] = ud * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
return (1); return (1);
} }
double angle_02pi(double A3) double angle_02pi(double A3)
......
Numeric/Numeric.h
+ 1
− 1
View file @ 8fb2313e
...@@ -65,7 +65,7 @@ double norme(double a[3]); ...@@ -65,7 +65,7 @@ double norme(double a[3]);
int sys2x2(double mat[2][2], double b[2], double res[2]); int sys2x2(double mat[2][2], double b[2], double res[2]);
int sys3x3(double mat[3][3], double b[3], double res[3], double *det); int sys3x3(double mat[3][3], double b[3], double res[3], double *det);
int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det); int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det);
int det3x3(double mat[3][3], double *det); double det3x3(double mat[3][3]);
int inv3x3(double mat[3][3], double inv[3][3], double *det); int inv3x3(double mat[3][3], double inv[3][3], double *det);
double angle_02pi(double A3); double angle_02pi(double A3);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment