LineNodeBasis: Convergence Test -- Seems OK

FunctionSpace/BasisGenerator.cpp

@@ -2,6 +2,8 @@
 #include "GmshDefines.h"
 #include "Exception.h"
 
+#include "LineNodeBasis.h"
+
 #include "QuadNodeBasis.h"
 #include "QuadEdgeBasis.h"
 
@@ -26,6 +28,7 @@ Basis* BasisGenerator::generate(int elementType,
 				int basisType, 
 				int order){
   switch(elementType){
+  case TYPE_LIN: return linGen(basisType, order);
   case TYPE_TRI: return triGen(basisType, order);
   case TYPE_QUA: return quaGen(basisType, order);
   case TYPE_TET: return tetGen(basisType, order);
@@ -36,6 +39,18 @@ Basis* BasisGenerator::generate(int elementType,
   }
 }
 
+Basis* BasisGenerator::linGen(int basisType, 
+			      int order){
+  switch(basisType){ 
+  case  0: return new LineNodeBasis(order);
+  case  1: throw Exception("1-form not implemented on Lines");
+  case  2: throw Exception("2-form not implemented on Lines");
+  case  3: throw Exception("3-form not implemented on Lines");
+
+  default: throw Exception("There is no %d-form", basisType);
+  }  
+}
+
 Basis* BasisGenerator::triGen(int basisType, 
 			      int order){
   switch(basisType){

FunctionSpace/BasisGenerator.h

@@ -24,6 +24,7 @@ class BasisGenerator{
 			 int basisType, 
 			 int order);
 
+  static Basis* linGen(int basisType, int order);
   static Basis* triGen(int basisType, int order);
   static Basis* quaGen(int basisType, int order);
   static Basis* tetGen(int basisType, int order);
@@ -56,6 +57,7 @@ class BasisGenerator{
    @em instantiated Basis
 
    @note Element types are:
+   @li @c TYPE_LIN for Lines
    @li @c TYPE_TRI for Triangles
    @li @c TYPE_QUA for Quadrangles
    @li @c TYPE_TET for Tetrahedrons
@@ -68,6 +70,23 @@ class BasisGenerator{
    @li @c 3 for 3-Form
    **
 
+   @fn BasisGenerator::linGen
+   @param basisType The Basis type
+   @param order The order or the requested Basis
+
+   This method will @em instanciate the requested Basis,
+   with a @em Line for support
+   
+   @return Returns a @em pointer to a newly 
+   @em instantiated Basis
+
+   @note Basis types are:
+   @li @c 0 for 0-Form
+   @li @c 1 for 1-Form
+   @li @c 2 for 2-Form
+   @li @c 3 for 3-Form
+   **
+
    @fn BasisGenerator::triGen
    @param basisType The Basis type
    @param order The order or the requested Basis

FunctionSpace/FunctionSpace.cpp

@@ -321,34 +321,40 @@ FunctionSpace::locBasis(const MElement& element,
   // Edge Based
   // Number of basis function *per* edge
   //  --> should always be an integer !
-  const unsigned int nEdge     = edgeClosure.size();
-  const unsigned int nFPerEdge = nFEdge / nEdge;
-  unsigned int fEdge = 0;
+  const unsigned int nEdge = edgeClosure.size();
 
-  for(unsigned int j = 0; j < nFPerEdge; j++){
-    for(unsigned int k = 0; k < nEdge; k++){
-      fun[i] = 
-	&basis.getEdgeFunction(edgeClosure[k], fEdge);
-      
-      fEdge++;
-      i++;
+  if(nEdge){
+    const unsigned int nFPerEdge = nFEdge / nEdge;
+    unsigned int fEdge = 0;
+    
+    for(unsigned int j = 0; j < nFPerEdge; j++){
+      for(unsigned int k = 0; k < nEdge; k++){
+	fun[i] = 
+	  &basis.getEdgeFunction(edgeClosure[k], fEdge);
+	
+	fEdge++;
+	i++;
+      }
     }
   }
 
   // Face Based
   // Number of basis function *per* face
   //  --> should always be an integer !
-  const unsigned int nFace     = faceClosure.size();
-  const unsigned int nFPerFace = nFFace / nFace;
-  unsigned int fFace = 0;
-
-  for(unsigned int j = 0; j < nFPerFace; j++){
-    for(unsigned int k = 0; k < nFace; k++){
-      fun[i] = 
-	&basis.getFaceFunction(faceClosure[k], fFace);
-      
-      fFace++;
-      i++;
+  const unsigned int nFace = faceClosure.size();
+  
+  if(nFace){
+    const unsigned int nFPerFace = nFFace / nFace;
+    unsigned int fFace = 0;
+    
+    for(unsigned int j = 0; j < nFPerFace; j++){
+      for(unsigned int k = 0; k < nFace; k++){
+	fun[i] = 
+	  &basis.getFaceFunction(faceClosure[k], fFace);
+	
+	fFace++;
+	i++;
+      }
     }
   }
 
@@ -390,16 +396,19 @@ FunctionSpace::locBasis(const MElement& element,
   // Number of basis function *per* edge
   //  --> should always be an integer !
   const unsigned int nEdge     = edgeClosure.size();
-  const unsigned int nFPerEdge = nFEdge / nEdge;
-  unsigned int fEdge = 0;
 
-  for(unsigned int j = 0; j < nFPerEdge; j++){
-    for(unsigned int k = 0; k < nEdge; k++){
-      fun[i] = 
-	&basis.getEdgeFunction(edgeClosure[k], fEdge);
-      
-      fEdge++;
-      i++;
+  if(nEdge){
+    const unsigned int nFPerEdge = nFEdge / nEdge;
+    unsigned int fEdge = 0;
+    
+    for(unsigned int j = 0; j < nFPerEdge; j++){
+      for(unsigned int k = 0; k < nEdge; k++){
+	fun[i] = 
+	  &basis.getEdgeFunction(edgeClosure[k], fEdge);
+	
+	fEdge++;
+	i++;
+      }
     }
   }
 
@@ -407,16 +416,19 @@ FunctionSpace::locBasis(const MElement& element,
   // Number of basis function *per* face
   //  --> should always be an integer !
   const unsigned int nFace     = faceClosure.size();
-  const unsigned int nFPerFace = nFFace / nFace;
-  unsigned int fFace = 0;
 
-  for(unsigned int j = 0; j < nFPerFace; j++){
-    for(unsigned int k = 0; k < nFace; k++){
-      fun[i] = 
-	&basis.getFaceFunction(faceClosure[k], fFace);
-      
-      fFace++;
-      i++;
+  if(nFace){
+    const unsigned int nFPerFace = nFFace / nFace;
+    unsigned int fFace = 0;
+    
+    for(unsigned int j = 0; j < nFPerFace; j++){
+      for(unsigned int k = 0; k < nFace; k++){
+	fun[i] = 
+	  &basis.getFaceFunction(faceClosure[k], fFace);
+	
+	fFace++;
+	i++;
+      }
     }
   }
 

FunctionSpace/LineNodeBasis.cpp

@@ -36,16 +36,14 @@ LineNodeBasis::LineNodeBasis(int order){
 
   // Vertex Based (Lagrange) //
   (*node)[0] = 
-    new Polynomial((Polynomial(1, 0, 0, 0) - 
-		    Polynomial(1, 1, 0, 0)) *
-		   0.5);
+    new Polynomial(Polynomial(0.5, 0, 0, 0) - 
+		   Polynomial(0.5, 1, 0, 0));
   
   (*node)[1] = 
-    new Polynomial((Polynomial(1, 0, 0, 0) + 
-		    Polynomial(1, 1, 0, 0)) *
-		   0.5);
-  
-  
+    new Polynomial(Polynomial(0.5, 0, 0, 0) + 
+		   Polynomial(0.5, 1, 0, 0));
+
+
   // Edge Based //
   const int permutation[2] = {0, 1};
   const Polynomial    x[2] = {