diff --git a/Mesh/Element.cpp b/Mesh/Element.cpp
index 0fde9f97dc8eadaf979e6ce1ae16a5260b60e0a3..11b7f96eb786f1af4a29e1413c1488d5b75e9ad0 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 d50ccdaaf219b44429ad72d86fb9c7c025decb69..8bfe99f042ff007f0dec0f7061f6c861b1c276cc 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 a34f6674205aa731056a3863d943ee78b674bf51..f819d0dfa046bac0d76b73587af72681132248d3 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 e707d2a38215749fb3614087e445bac4318d3bd3..472b6ee4ed56fb25423759b4e099e608c89e910b 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 bddcf1eef4e1b9f4f0579bfa4a34401044aea99f..6ada0a0bf1cb46b4b84ca0a5c32b7abf4806b112 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);