Adding Lagrange Basis For Tetrahedron -- Adaptive Visu *ONLY* -- Seems Ok with Poisson on a Cube

FunctionSpace/BasisGenerator.cpp

@@ -17,6 +17,7 @@
 #include "TetNodeBasis.h"
 #include "TetEdgeBasis.h"
 #include "TetNedelecBasis.h"
+#include "TetLagrangeBasis.h"
 
 #include "HexNodeBasis.h"
 #include "HexEdgeBasis.h"
@@ -70,7 +71,7 @@ BasisLocal* BasisGenerator::generateLagrange(unsigned int elementType,
   case TYPE_LIN: throw Exception("Lagrange Basis on Lines not Implemented");
   case TYPE_TRI: return new TriLagrangeBasis(order);
   case TYPE_QUA: throw Exception("Lagrange Basis on Quads not Implemented");
-  case TYPE_TET: throw Exception("Lagrange Basis on Tets not Implemented");
+  case TYPE_TET: return new TetLagrangeBasis(order);
   case TYPE_HEX: throw Exception("Lagrange Basis on Hexs not Implemented");
 
   default: throw Exception("Unknown Element Type (%d) for Basis Generation",

FunctionSpace/CMakeLists.txt

@@ -38,6 +38,7 @@ set(SRC
   TetNodeBasis.cpp
   TetEdgeBasis.cpp
   TetNedelecBasis.cpp
+  TetLagrangeBasis.cpp
 
   FunctionSpace.cpp
   FunctionSpaceScalar.cpp

FunctionSpace/TetEdgeBasis.cpp

@@ -82,9 +82,9 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
 	  tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
 	  tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
 
-	  tmp1[0].sub(tmp2[0]);
-	  tmp1[1].sub(tmp2[1]);
-	  tmp1[2].sub(tmp2[2]);
+	  tmp2[0].sub(tmp1[0]);
+	  tmp2[1].sub(tmp1[1]);
+	  tmp2[2].sub(tmp1[2]);
 
 	  basis[s][i] = new vector<Polynomial>(tmp1);
 	}

FunctionSpace/TetLagrangeBasis.cpp

+#include "Exception.h"
+#include "TetLagrangeBasis.h"
+
+TetLagrangeBasis::TetLagrangeBasis(unsigned int order){
+  // Set Basis Type //
+  this->order = order;
+
+  type = 0;
+  dim  = 3;
+
+  nVertex   = 4;
+  nEdge     = 6 * (order - 1);
+  nFace     = 2 * (order - 1) * (order - 2);
+  nCell     =     (order - 1) * (order - 2) * (order - 3) / 6;
+  nFunction = nVertex + nEdge + nFace + nCell;
+
+  // Init polynomialBasis //
+  lBasis = new polynomialBasis(getTag(order));
+
+  // Init Lagrange Point //
+  lPoint = new fullMatrix<double>(tetPoint(order));
+}
+
+TetLagrangeBasis::~TetLagrangeBasis(void){
+  delete lBasis;
+  delete lPoint;
+}
+
+unsigned int TetLagrangeBasis::getTag(unsigned int order){
+  unsigned int tag = nodalBasis::getTag(TYPE_TET, order, false);
+
+  if(tag)
+    return tag;
+
+  else
+    throw Exception
+      ("Can't instanciate an order %d Lagrangian Basis for a Tetrahedron",
+       order);
+}
+
+fullMatrix<double> TetLagrangeBasis::
+tetPoint(unsigned int order){
+  unsigned int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
+  fullMatrix<double> point(nbPoints, 3);
+
+  double overOrder = (order == 0 ? 1. : 1. / order);
+
+  point(0, 0) = 0.;
+  point(0, 1) = 0.;
+  point(0, 2) = 0.;
+
+  if(order > 0){
+    point(1, 0) = order;
+    point(1, 1) = 0;
+    point(1, 2) = 0;
+
+    point(2, 0) = 0.;
+    point(2, 1) = order;
+    point(2, 2) = 0.;
+
+    point(3, 0) = 0.;
+    point(3, 1) = 0.;
+    point(3, 2) = order;
+
+    if(order > 1){
+      for(unsigned int k = 0; k < (order - 1); k++){
+        point(4 + k,                   0) = k + 1;
+        point(4 +      order - 1  + k, 0) = order - 1 - k;
+        point(4 + 2 * (order - 1) + k, 0) = 0.;
+        point(4 + 3 * (order - 1) + k, 0) = 0.;
+        point(4 + 4 * (order - 1) + k, 0) = 0.;
+        point(4 + 5 * (order - 1) + k, 0) = k+1;
+
+        point(4 + k,                   1) = 0.;
+        point(4 +      order - 1  + k, 1) = k + 1;
+        point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+        point(4 + 3 * (order - 1) + k, 1) = 0.;
+        point(4 + 4 * (order - 1) + k, 1) = k+1;
+        point(4 + 5 * (order - 1) + k, 1) = 0.;
+
+        point(4 + k,                   2) = 0.;
+        point(4 +      order - 1  + k, 2) = 0.;
+        point(4 + 2 * (order - 1) + k, 2) = 0.;
+        point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
+        point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+        point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
+      }
+
+      if(order > 2){
+        unsigned int ns = 4 + 6 * (order - 1);
+        unsigned int nbdofface = nbdoftriangle(order - 3);
+
+        double *u = new double[nbdofface];
+        double *v = new double[nbdofface];
+        double *w = new double[nbdofface];
+
+        nodepositionface0(order - 3, u, v, w);
+
+        // u-v plane
+
+        for(unsigned int i = 0; i < nbdofface; i++){
+          point(ns + i, 0) = u[i] * (order - 3) + 1.;
+          point(ns + i, 1) = v[i] * (order - 3) + 1.;
+          point(ns + i, 2) = w[i] * (order - 3);
+        }
+
+        ns = ns + nbdofface;
+
+        // u-w plane
+
+        nodepositionface1(order - 3, u, v, w);
+
+        for(unsigned int i=0; i < nbdofface; i++){
+          point(ns + i, 0) = u[i] * (order - 3) + 1.;
+          point(ns + i, 1) = v[i] * (order - 3) ;
+          point(ns + i, 2) = w[i] * (order - 3) + 1.;
+        }
+
+        // v-w plane
+
+        ns = ns + nbdofface;
+
+        nodepositionface2(order - 3, u, v, w);
+
+        for(unsigned int i = 0; i < nbdofface; i++){
+          point(ns + i, 0) = u[i] * (order - 3);
+          point(ns + i, 1) = v[i] * (order - 3) + 1.;
+          point(ns + i, 2) = w[i] * (order - 3) + 1.;
+        }
+
+        // u-v-w plane
+
+        ns = ns + nbdofface;
+
+        nodepositionface3(order - 3, u, v, w);
+
+        for(unsigned int i = 0; i < nbdofface; i++){
+          point(ns + i, 0) = u[i] * (order - 3) + 1.;
+          point(ns + i, 1) = v[i] * (order - 3) + 1.;
+          point(ns + i, 2) = w[i] * (order - 3) + 1.;
+        }
+
+        ns = ns + nbdofface;
+
+        delete [] u;
+        delete [] v;
+        delete [] w;
+
+        if(order > 3){
+          fullMatrix<double> interior = tetPoint(order - 4);
+
+          for(int k = 0; k < interior.size1(); k++){
+            point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
+            point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
+            point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
+          }
+        }
+      }
+    }
+  }
+
+  point.scale(overOrder);
+  return point;
+}
+
+unsigned int TetLagrangeBasis::nbdoftriangle(unsigned int order){
+  return (order + 1) * (order + 2) / 2;
+}
+
+void TetLagrangeBasis::
+nodepositionface0(unsigned int order, double *u, double *v, double *w){
+
+  unsigned int ndofT = nbdoftriangle(order);
+
+  if(order == 0){
+    u[0] = 0.;
+    v[0] = 0.;
+    w[0] = 0.;
+    return;
+  }
+
+  u[0]= 0.;    v[0]= 0.;    w[0] = 0.;
+  u[1]= 0.;    v[1]= order; w[1] = 0.;
+  u[2]= order; v[2]= 0.;    w[2] = 0.;
+
+  // edges
+  for(unsigned int k = 0; k < (order - 1); k++){
+    u[3 + k] = 0.;
+    v[3 + k] = k + 1;
+    w[3 + k] = 0.;
+
+    u[3 + order - 1 + k] = k + 1;
+    v[3 + order - 1 + k] = order - 1 - k ;
+    w[3 + order - 1 + k] = 0.;
+
+    u[3 + 2 * (order - 1) + k] = order - 1 - k;
+    v[3 + 2 * (order - 1) + k] = 0.;
+    w[3 + 2 * (order - 1) + k] = 0.;
+  }
+
+  if(order > 2){
+    unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+    nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+                      &w[3 + 3* (order - 1)]);
+
+    for(unsigned int k = 0; k < nbdoftemp; k++){
+      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
+    }
+  }
+
+  for(unsigned int k = 0; k < ndofT; k++){
+    u[k] = u[k] / order;
+    v[k] = v[k] / order;
+    w[k] = w[k] / order;
+  }
+}
+
+void TetLagrangeBasis::nodepositionface1
+(unsigned int order, double *u, double *v, double *w){
+
+  unsigned int ndofT = nbdoftriangle(order);
+
+  if(order == 0){
+    u[0] = 0.;
+    v[0] = 0.;
+    w[0] = 0.;
+    return;
+  }
+
+  u[0] = 0.;    v[0]= 0.;  w[0] = 0.;
+  u[1] = order; v[1]= 0.;  w[1] = 0.;
+  u[2] = 0.;    v[2]= 0.;  w[2] = order;
+
+  // edges
+  for(unsigned int k = 0; k < (order - 1); k++){
+    u[3 + k] = k + 1;
+    v[3 + k] = 0.;
+    w[3 + k] = 0.;
+
+    u[3 + order - 1 + k] = order - 1 - k;
+    v[3 + order - 1 + k] = 0.;
+    w[3 + order - 1+ k ] = k + 1;
+
+    u[3 + 2 * (order - 1) + k] = 0. ;
+    v[3 + 2 * (order - 1) + k] = 0.;
+    w[3 + 2 * (order - 1) + k] = order - 1 - k;
+  }
+
+  if(order > 2){
+    unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+    nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
+                      &w[3 + 3 * (order - 1)]);
+
+    for(unsigned int k = 0; k < nbdoftemp; k++){
+      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
+      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+    }
+  }
+
+  for (unsigned int k = 0; k < ndofT; k++){
+    u[k] = u[k] / order;
+    v[k] = v[k] / order;
+    w[k] = w[k] / order;
+  }
+}
+
+void TetLagrangeBasis::
+nodepositionface2(unsigned int order, double *u, double *v, double *w){
+
+  unsigned int ndofT = nbdoftriangle(order);
+
+  if(order == 0){
+    u[0] = 0.;
+    v[0] = 0.;
+    w[0] = 0.;
+    return;
+  }
+
+  u[0]= 0.; v[0]= 0.;    w[0] = 0.;
+  u[1]= 0.; v[1]= 0.;    w[1] = order;
+  u[2]= 0.; v[2]= order; w[2] = 0.;
+
+  // edges
+  for(unsigned int k = 0; k < (order - 1); k++){
+    u[3 + k] = 0.;
+    v[3 + k] = 0.;
+    w[3 + k] = k + 1;
+
+    u[3 + order - 1 + k] = 0.;
+    v[3 + order - 1 + k] = k + 1;
+    w[3 + order - 1 + k] = order - 1 - k;
+
+    u[3 + 2 * (order - 1) + k] = 0.;
+    v[3 + 2 * (order - 1) + k] = order - 1 - k;
+    w[3 + 2 * (order - 1) + k] = 0.;
+  }
+
+  if(order > 2){
+    unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+    nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+                      &w[3 + 3 * (order - 1)]);
+
+    for(unsigned int k = 0; k < nbdoftemp; k++){
+      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
+      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+    }
+  }
+
+  for(unsigned int k = 0; k < ndofT; k++){
+    u[k] = u[k] / order;
+    v[k] = v[k] / order;
+    w[k] = w[k] / order;
+  }
+}
+
+void TetLagrangeBasis::
+nodepositionface3(unsigned int order,  double *u,  double *v,  double *w){
+
+  unsigned int ndofT = nbdoftriangle(order);
+
+  if(order == 0){
+    u[0] = 0.;
+    v[0] = 0.;
+    w[0] = 0.;
+    return;
+  }
+
+  u[0]= 0.;    v[0]= 0.;    w[0] = order;
+  u[1]= order; v[1]= 0.;    w[1] = 0.;
+  u[2]= 0.;    v[2]= order; w[2] = 0.;
+
+  // edges
+  for(unsigned int k = 0; k < (order - 1); k++){
+    u[3 + k] = k + 1;
+    v[3 + k] = 0.;
+    w[3 + k] = order - 1 - k;
+
+    u[3 + order - 1 + k] = order - 1 - k;
+    v[3 + order - 1 + k] = k + 1;
+    w[3 + order - 1 + k] = 0.;
+
+    u[3 + 2 * (order - 1) + k] = 0.;
+    v[3 + 2 * (order - 1) + k] = order - 1 - k;
+    w[3 + 2 * (order - 1) + k] = k + 1;
+  }
+
+  if(order > 2){
+    unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+    nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+                       &w[3 + 3 * (order - 1)]);
+
+    for(unsigned int k = 0; k < nbdoftemp; k++){
+      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+    }
+  }
+
+  for (unsigned int k = 0; k < ndofT; k++){
+    u[k] = u[k] / order;
+    v[k] = v[k] / order;
+    w[k] = w[k] / order;
+  }
+}

FunctionSpace/TetLagrangeBasis.h

+#ifndef _TETLAGRANGEBASIS_H_
+#define _TETLAGRANGEBASIS_H_
+
+#include "BasisLagrange.h"
+
+/**
+   @class TetLagrangeBasis
+   @brief Lagrange Basis for Tetrahedra
+
+   This class can instantiate a @em Lagrange @em Basis
+   for a Tetrahedron and for a given Order.@n
+
+   It uses
+   <a href="http://geuz.org/gmsh/">gmsh</a> Basis.
+ */
+
+class TetLagrangeBasis: public BasisLagrange{
+ public:
+  //! @param order A natural number
+  //!
+  //! Returns a new TetLagrangeBasis
+  //! of the given Order
+  TetLagrangeBasis(unsigned int order);
+
+  //! Deletes this Basis
+  //!
+  virtual ~TetLagrangeBasis(void);
+
+ private:
+  //! @param order A natural number
+  //! @return Returns the @em tag of a @em Tetrahedron of
+  //! the given order
+  static unsigned int getTag(unsigned int order);
+
+  //! @param order A natural number
+  //! @return Returns Lagrangian Points on a Tetrahedron
+  //! for the given Order
+  static fullMatrix<double> tetPoint(unsigned int order);
+
+  //! Unknown gmsh function
+  static unsigned int nbdoftriangle(unsigned int order);
+
+  //! Unknown gmsh function
+  static void nodepositionface0(unsigned int order,
+                                double *u, double *v, double *w);
+
+  //! Unknown gmsh function
+  static void nodepositionface1(unsigned int order,
+                                double *u, double *v, double *w);
+
+  //! Unknown gmsh function
+  static void nodepositionface2(unsigned int order,
+                                double *u, double *v, double *w);
+
+  //! Unknown gmsh function
+  static void nodepositionface3(unsigned int order,
+                                double *u, double *v, double *w);
+};
+
+#endif

FunctionSpace/TriLagrangeBasis.cpp

@@ -27,24 +27,16 @@ TriLagrangeBasis::~TriLagrangeBasis(void){
 }
 
 unsigned int TriLagrangeBasis::getTag(unsigned int order){
-  switch(order){
-  case 1:  return MSH_TRI_3;
-  case 2:  return MSH_TRI_6;
-  case 3:  return MSH_TRI_10;
-  case 4:  return MSH_TRI_15;
-  case 5:  return MSH_TRI_21;
-  case 6:  return MSH_TRI_28;
-  case 7:  return MSH_TRI_36;
-  case 8:  return MSH_TRI_45;
-  case 9:  return MSH_TRI_55;
-  case 10: return MSH_TRI_66;
-
-  default:
+  unsigned int tag = nodalBasis::getTag(TYPE_TRI, order, false);
+
+  if(tag)
+    return tag;
+
+  else
     throw Exception
       ("Can't instanciate an order %d Lagrangian Basis for a Triangle",
        order);
 }
-}
 
 fullMatrix<double> TriLagrangeBasis::
 triPoint(unsigned int order){