Simplifying Implementation + Doc

FunctionSpace/FunctionSpace.cpp

@@ -19,14 +19,8 @@ FunctionSpace::~FunctionSpace(void){
   delete basis;
 
   // Dof //
-  if(dof){
-    set<const Dof*>::iterator dStop = dof->end();
-    set<const Dof*>::iterator dIt   = dof->begin();
-
-    for(; dIt != dStop; dIt++)
-      delete *dIt;
+  if(dof)
     delete dof;
-  }
 
   // Group //
   if(group){
@@ -96,7 +90,7 @@ void FunctionSpace::buildDof(void){
   const vector<const MElement*>& element = goe->getAll();
 
   // Init Struct //
-  dof      = new set<const Dof*, DofComparator>;
+  dof      = new set<Dof>;
   group    = new vector<GroupOfDof*>(nElement);
   eToGod   = new map<const MElement*,
                      const GroupOfDof*,
@@ -109,13 +103,15 @@ void FunctionSpace::buildDof(void){
     size_t nDof       = myDof.size();
 
     // Add Dof
-    vector<const Dof*> trueDof(nDof);
-
     for(size_t j = 0; j < nDof; j++)
-      insertDof(myDof[j], trueDof, j);
+      dof->insert(myDof[j]);
+    //vector<Dof> trueDof(nDof);
+
+    //for(size_t j = 0; j < nDof; j++)
+    //insertDof(myDof[j], trueDof, j);
 
     // Create new GroupOfDof
-    GroupOfDof* god = new GroupOfDof(*(element[i]), trueDof);
+    GroupOfDof* god = new GroupOfDof(*(element[i]), myDof);
     (*group)[i]     = god;
 
     // Map GOD
@@ -123,7 +119,7 @@ void FunctionSpace::buildDof(void){
                    (element[i], god));
   }
 }
-
+/*
 void FunctionSpace::insertDof(Dof& d,
                               vector<const Dof*>& trueDof,
                               size_t index){
@@ -146,7 +142,7 @@ void FunctionSpace::insertDof(Dof& d,
     trueDof[index] = *(p.first);
   }
 }
-
+*/
 vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
   // Const_Cast //
   MElement& element = const_cast<MElement&>(elem);
@@ -165,6 +161,7 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
 
   // New Element
   MElementFactory factory;
+  //MElement* permElement = factory.create(elem.getLowOrderTypeForMSH(), vertex);
   MElement* permElement = factory.create(elem.getTypeForMSH(), vertex);
 
   // Edge & Face from Permuted Element //

FunctionSpace/FunctionSpace.h

@@ -59,8 +59,8 @@ class FunctionSpace{
   bool scalar;
 
   // Dofs //
-  std::set<const Dof*, DofComparator>*     dof;
-  std::vector<GroupOfDof*>*                group;
+  std::set<Dof>*            dof;
+  std::vector<GroupOfDof*>* group;
   std::map<
     const MElement*,
     const GroupOfDof*, ElementComparator>* eToGod;
@@ -80,7 +80,7 @@ class FunctionSpace{
   std::vector<Dof> getKeys(const MEdge& edge) const;
   std::vector<Dof> getKeys(const MFace& face) const;
 
-  const std::vector<const Dof*>   getAllDofs(void) const;
+  const std::vector<Dof>          getAllDofs(void) const;
   const std::vector<GroupOfDof*>& getAllGroups(void) const;
 
   const GroupOfDof& getGoDFromElement(const MElement& element) const;
@@ -209,8 +209,8 @@ inline size_t FunctionSpace::groupNumber(void) const{
   return group->size();
 }
 
-inline const std::vector<const Dof*> FunctionSpace::getAllDofs(void) const{
-  return std::vector<const Dof*>(dof->begin(), dof->end());
+inline const std::vector<Dof> FunctionSpace::getAllDofs(void) const{
+  return std::vector<Dof>(dof->begin(), dof->end());
 }
 
 inline const std::vector<GroupOfDof*>& FunctionSpace::getAllGroups(void) const{

Geo/MElement.h

@@ -347,6 +347,7 @@ class MElement
   // info for specific IO formats (returning 0 means that the element
   // is not implemented in that format)
   virtual int getTypeForMSH() const { return 0; }
+  virtual int getLowOrderTypeForMSH() const { return 0; }
   virtual int getTypeForUNV() const { return 0; }
   virtual int getTypeForVTK() const { return 0; }
   virtual const char *getStringForPOS() const { return 0; }

Geo/MLine.h

@@ -65,6 +65,7 @@ class MLine : public MElement {
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){}
   virtual int getType() const { return TYPE_LIN; }
   virtual int getTypeForMSH() const { return MSH_LIN_2; }
+  virtual int getLowOrderTypeForMSH() const { return MSH_LIN_2; }
   virtual int getTypeForUNV() const { return 21; } // linear beam
   virtual int getTypeForVTK() const { return 3; }
   virtual const char *getStringForPOS() const { return "SL"; }

Geo/MQuadrangle.h

@@ -111,6 +111,7 @@ class MQuadrangle : public MElement {
   }
   virtual int getType() const { return TYPE_QUA; }
   virtual int getTypeForMSH() const { return MSH_QUA_4; }
+  virtual int getLowOrderTypeForMSH() const { return MSH_QUA_4; }
   virtual int getTypeForUNV() const { return 94; } // thin shell linear quadrilateral
   virtual int getTypeForVTK() const { return 9; }
   virtual const char *getStringForPOS() const { return "SQ"; }

Geo/MTetrahedron.h

@@ -100,6 +100,7 @@ class MTetrahedron : public MElement {
   }
   virtual int getType() const { return TYPE_TET; }
   virtual int getTypeForMSH() const { return MSH_TET_4; }
+  virtual int getLowOrderTypeForMSH() const { return MSH_TET_4; }
   virtual int getTypeForUNV() const { return 111; } // solid linear tetrahedron
   virtual int getTypeForVTK() const { return 10; }
   virtual const char *getStringForPOS() const { return "SS"; }

Geo/MTriangle.h

@@ -114,6 +114,7 @@ class MTriangle : public MElement {
   }
   virtual int getType() const { return TYPE_TRI; }
   virtual int getTypeForMSH() const { return MSH_TRI_3; }
+  virtual int getLowOrderTypeForMSH() const { return MSH_TRI_3; }
   virtual int getTypeForUNV() const { return 91; } // thin shell linear triangle
   virtual int getTypeForVTK() const { return 5; }
   virtual const char *getStringForPOS() const { return "ST"; }

Numeric/BasisFactory.cpp

@@ -46,9 +46,17 @@ const nodalBasis* BasisFactory::getNodalBasis(int tag)
   }
 
   // FIXME: check if already exists to deallocate if necessary
-  fs.insert(std::make_pair(tag, F));
+  std::pair<std::map<int, nodalBasis*>::const_iterator, bool> inserted;
 
-  return F;
+  #pragma omp critical
+    {
+      inserted = fs.insert(std::make_pair(tag, F));
+
+      if (!inserted.second)
+        delete F;
+    }
+
+  return inserted.first->second;
 }
 
 const JacobianBasis* BasisFactory::getJacobianBasis(int tag)

Numeric/fullMatrix.h

@@ -745,6 +745,25 @@ class fullMatrix
     printf("};\n");
   }
 
+
+  void printMatlab(const std::string name = "A") const
+  {
+    int ni = size1();
+    int nj = size2();
+
+    printf("%s = [", name.c_str());
+
+    for(int I = 0; I < ni; I++){
+      for(int J = 0; J < nj; J++){
+        printf("%+e,", (*this)(I, J));
+      }
+
+      printf(";");
+    }
+
+    printf("]\n");
+  }
+
   // specific functions for dgshell
   void mult_naiveBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol,
                        const int alpha, const int beta, fullVector<scalar> &c,