Still a problem with SteadyWave in a Cube -- Try Cavity

FunctionSpace/BasisGenerator.cpp

@@ -16,6 +16,7 @@
 
 #include "TetNodeBasis.h"
 #include "TetEdgeBasis.h"
+#include "TetNedelecBasis.h"
 
 #include "HexNodeBasis.h"
 #include "HexEdgeBasis.h"
@@ -27,118 +28,121 @@ BasisGenerator::BasisGenerator(void){
 BasisGenerator::~BasisGenerator(void){
 }
 
-BasisLocal* BasisGenerator::generate(unsigned int elementType, 
-				     unsigned int basisType, 
+BasisLocal* BasisGenerator::generate(unsigned int elementType,
+				     unsigned int basisType,
 				     unsigned int order,
 				     std::string family){
-  
+
   if(!family.compare(std::string("hierarchical")))
     return generateHierarchical(elementType, basisType, order);
-  
+
   else if(!family.compare(std::string("lagrange")))
     return generateLagrange(elementType, basisType, order);
-  
+
   else
     throw Exception("Unknwown Basis Family: %s", family.c_str());
 }
 
-BasisLocal* BasisGenerator::generateHierarchical(unsigned int elementType, 
-						 unsigned int basisType, 
-						 unsigned int order){  
+BasisLocal* BasisGenerator::generateHierarchical(unsigned int elementType,
+						 unsigned int basisType,
+						 unsigned int order){
   switch(elementType){
   case TYPE_LIN: return linHierarchicalGen(basisType, order);
   case TYPE_TRI: return triHierarchicalGen(basisType, order);
   case TYPE_QUA: return quaHierarchicalGen(basisType, order);
   case TYPE_TET: return tetHierarchicalGen(basisType, order);
   case TYPE_HEX: return hexHierarchicalGen(basisType, order);
-    
-  default: throw Exception("Unknown Element Type (%d) for Basis Generation", 
+
+  default: throw Exception("Unknown Element Type (%d) for Basis Generation",
 			   elementType);
   }
 }
 
-BasisLocal* BasisGenerator::generateLagrange(unsigned int elementType, 
-					     unsigned int basisType, 
+BasisLocal* BasisGenerator::generateLagrange(unsigned int elementType,
+					     unsigned int basisType,
 					     unsigned int order){
   if(basisType != 0)
-    throw 
+    throw
       Exception("Cannot Have a %d-Form Lagrange Basis (0-Form only)",
 		basisType);
-  
+
   switch(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_HEX: throw Exception("Lagrange Basis on Hexs not Implemented");
-    
-  default: throw Exception("Unknown Element Type (%d) for Basis Generation", 
+
+  default: throw Exception("Unknown Element Type (%d) for Basis Generation",
 			   elementType);
   }
 }
 
-BasisLocal* BasisGenerator::linHierarchicalGen(unsigned int basisType, 
+BasisLocal* BasisGenerator::linHierarchicalGen(unsigned int basisType,
 					       unsigned int order){
-  switch(basisType){ 
+  switch(basisType){
   case  0: return new LineNodeBasis(order);
-  case  1: 
+  case  1:
     if (order == 0) return new LineNedelecBasis();
     else            return new LineEdgeBasis(order);
-    
+
   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);
-  }  
+  }
 }
 
-BasisLocal* BasisGenerator::triHierarchicalGen(unsigned int basisType, 
+BasisLocal* BasisGenerator::triHierarchicalGen(unsigned int basisType,
 					       unsigned int order){
   switch(basisType){
   case  0: return new TriNodeBasis(order);
-  case  1: 
+  case  1:
     if (order == 0) return new TriNedelecBasis();
     else            return new TriEdgeBasis(order);
-    
+
   case  2: throw Exception("2-form not implemented on Triangles");
   case  3: throw Exception("3-form not implemented on Triangles");
-    
+
   default: throw Exception("There is no %d-form", basisType);
-  }  
+  }
 }
 
-BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType, 
+BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType,
 					       unsigned int order){
   switch(basisType){
     //case  0: return new QuadNodeBasis(order);
     //case  1: return new QuadEdgeBasis(order);
   case  2: throw Exception("2-form not implemented on Quads");
   case  3: throw Exception("3-form not implemented on Quads");
-    
+
   default: throw Exception("There is no %d-form", basisType);
-  }  
+  }
 }
 
-BasisLocal* BasisGenerator::tetHierarchicalGen(unsigned int basisType, 
+BasisLocal* BasisGenerator::tetHierarchicalGen(unsigned int basisType,
 					       unsigned int order){
   switch(basisType){
   case  0: return new TetNodeBasis(order);
-  case  1: return new TetEdgeBasis(order);
+  case  1:
+    if (order == 0) return new TetNedelecBasis();
+    else            return new TetEdgeBasis(order);
+
   case  2: throw Exception("2-form not implemented on Tetrahedrons");
   case  3: throw Exception("3-form not implemented on Tetrahedrons");
-    
+
   default: throw Exception("There is no %d-form", basisType);
-  }  
+  }
 }
 
-BasisLocal* BasisGenerator::hexHierarchicalGen(unsigned int basisType, 
+BasisLocal* BasisGenerator::hexHierarchicalGen(unsigned int basisType,
 					       unsigned int order){
   switch(basisType){
     //case  0: return new HexNodeBasis(order);
     //case  1: return new HexEdgeBasis(order);
   case  2: throw Exception("2-form not implemented on Hexs");
   case  3: throw Exception("3-form not implemented on Hexs");
-    
+
   default: throw Exception("There is no %d-form", basisType);
-  }  
+  }
 }

FunctionSpace/CMakeLists.txt

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

FunctionSpace/TetNedelecBasis.cpp

+#include "TetNedelecBasis.h"
+#include "TetReferenceSpace.h"
+
+using namespace std;
+
+TetNedelecBasis::TetNedelecBasis(void){
+  // Reference Space //
+  refSpace  = new TetReferenceSpace;
+  nRefSpace = refSpace->getNPermutation();
+
+  const vector<const vector<const vector<unsigned int>*>*>&
+    edgeV = refSpace->getAllEdge();
+
+  // Set Basis Type //
+  this->order = 1;
+
+  type = 1;
+  dim  = 3;
+
+  nVertex   = 0;
+  nEdge     = 6;
+  nFace     = 0;
+  nCell     = 0;
+  nFunction = 6;
+
+  // Lagrange Polynomial //
+  const Polynomial lagrange[4] =
+    {
+      Polynomial(Polynomial(1, 0, 0, 0) -
+		 Polynomial(1, 1, 0, 0) -
+		 Polynomial(1, 0, 1, 0) -
+		 Polynomial(1, 0, 0, 1)),
+
+      Polynomial(Polynomial(1, 1, 0, 0)),
+
+      Polynomial(Polynomial(1, 0, 1, 0)),
+
+      Polynomial(Polynomial(1, 0, 0, 1))
+    };
+
+
+  // Basis //
+  basis = new vector<Polynomial>**[nRefSpace];
+
+  for(unsigned int s = 0; s < nRefSpace; s++)
+    basis[s] = new vector<Polynomial>*[nFunction];
+
+  // Edge Based (Nedelec) //
+  for(unsigned int s = 0; s < nRefSpace; s++){
+    for(unsigned int e = 0; e < 6; e++){
+      vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
+      vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
+
+      tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+      tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+      tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+
+      tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+      tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+      tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+
+      tmp2[0].sub(tmp1[0]);
+      tmp2[1].sub(tmp1[1]);
+      tmp2[2].sub(tmp1[2]);
+
+      basis[s][e] = new vector<Polynomial>(tmp2);
+    }
+  }
+}
+
+TetNedelecBasis::~TetNedelecBasis(void){
+  // ReferenceSpace //
+  delete refSpace;
+
+  // Basis //
+  for(unsigned int i = 0; i < nRefSpace; i++){
+    for(unsigned int j = 0; j < nFunction; j++)
+      delete basis[i][j];
+
+    delete[] basis[i];
+  }
+
+  delete[] basis;
+}

FunctionSpace/TetNedelecBasis.h

+#ifndef _TETNEDELECBASIS_H_
+#define _TETNEDELECBASIS_H_
+
+#include "BasisHierarchical1From.h"
+
+/**
+   @class TetNedelecBasis
+   @brief A Nedelec Basis for Tetrahedra
+
+   This class can instantiate
+   a Nedelec Basis for Tetrahedra.
+*/
+
+class TetNedelecBasis: public BasisHierarchical1From{
+ public:
+  //! @param order The order of the Basis
+  //!
+  //! Returns a new Nedelec Basis for Tetrahedra
+  TetNedelecBasis(void);
+
+  //! Deletes this Basis
+  //!
+  virtual ~TetNedelecBasis(void);
+};
+
+#endif