From 8fb2313e85d174539029485ecb6429e770c81a0c Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 21 Jul 2004 22:19:56 +0000
Subject: [PATCH] - added Orientation() members for hexas, prisms and pyramids
- use det3x3() instead of ad-hoc code everywhere it makes sense
---
Mesh/Element.cpp | 45 ++++++++++++++++++++++++++++++++++++++++++++-
Mesh/Element.h | 3 +++
Mesh/Simplex.cpp | 6 ++----
Numeric/Numeric.cpp | 29 +++++++++++------------------
Numeric/Numeric.h | 2 +-
5 files changed, 61 insertions(+), 24 deletions(-)
diff --git a/Mesh/Element.cpp b/Mesh/Element.cpp
index 0fde9f97dc..11b7f96eb7 100644
--- a/Mesh/Element.cpp
+++ b/Mesh/Element.cpp
@@ -1,4 +1,4 @@
-// $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
//
@@ -21,6 +21,7 @@
#include "Gmsh.h"
#include "Mesh.h"
+#include "Numeric.h"
#include "Element.h"
int Element::TotalNumber = 0;
@@ -107,6 +108,20 @@ Hexahedron::Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4,
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()
{
if(VSUP) Free(VSUP);
@@ -160,6 +175,20 @@ Prism::Prism(Vertex *v1, Vertex *v2, Vertex *v3,
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()
{
if(VSUP) Free(VSUP);
@@ -200,6 +229,20 @@ Pyramid::Pyramid()
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)
{
iEnt = -1;
diff --git a/Mesh/Element.h b/Mesh/Element.h
index d50ccdaaf2..8bfe99f042 100644
--- a/Mesh/Element.h
+++ b/Mesh/Element.h
@@ -49,6 +49,7 @@ class Hexahedron : public Element{
Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4,
Vertex *v5, Vertex *v6, Vertex *v7, Vertex *v8);
~Hexahedron();
+ double Orientation();
};
class Prism : public Element{
@@ -58,6 +59,7 @@ class Prism : public Element{
Prism(Vertex *v1, Vertex *v2, Vertex *v3,
Vertex *v4, Vertex *v5, Vertex *v6);
~Prism();
+ double Orientation();
};
class Pyramid : public Element{
@@ -66,6 +68,7 @@ class Pyramid : public Element{
Pyramid();
Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5);
~Pyramid();
+ double Orientation();
};
// C interface
diff --git a/Mesh/Simplex.cpp b/Mesh/Simplex.cpp
index a34f667420..f819d0dfa0 100644
--- a/Mesh/Simplex.cpp
+++ b/Mesh/Simplex.cpp
@@ -1,4 +1,4 @@
-// $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
//
@@ -203,9 +203,7 @@ double Simplex::matsimpl(double mat[3][3])
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 (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
- 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]));
+ return det3x3(mat);
}
double Simplex::rhoin()
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index e707d2a382..472b6ee4ed 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1,4 +1,4 @@
-// $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
//
@@ -122,8 +122,7 @@ int sys2x2(double mat[2][2], double b[2], double res[2])
double det, ud, norm;
int i;
- norm =
- DSQR(mat[0][0]) + DSQR(mat[1][1]) + DSQR(mat[0][1]) + DSQR(mat[1][0]);
+ norm = 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];
@@ -146,14 +145,19 @@ int sys2x2(double mat[2][2], double b[2], double res[2])
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)
{
double ud;
int i;
- *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]);
+ *det = det3x3(mat);
if(*det == 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)
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)
{
double ud;
- *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]);
+ *det = det3x3(mat);
if(*det == 0.0)
return (0);
@@ -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][2] = ud * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
return (1);
-
}
double angle_02pi(double A3)
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index bddcf1eef4..6ada0a0bf1 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -65,7 +65,7 @@ double norme(double a[3]);
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_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);
double angle_02pi(double A3);
--
GitLab