Working on Dof Reordering -- New GroupOfDof Interface -- Cleaning some tabs...

Working on Dof Reordering -- New GroupOfDof Interface -- Cleaning some tabs (replaced by spaces) -- BasisLagrange Implementation of Function and Grad

FunctionSpace/Basis.h

@@ -77,13 +77,13 @@ class Basis{
   unsigned int getNCellBased(void) const;
   unsigned int getNFunction(void) const;
 
-  // Orientations //
+  // Reference Element //
   virtual unsigned int getNOrientation(void) const = 0;
   virtual unsigned int getOrientation(const MElement& element) const = 0;
 
-  // Functions Permutation //
-  virtual void getFunctionPermutation(const MElement& element,
-                                      unsigned int* indexPermutation) const = 0;
+  // Functions Ordering //
+  virtual std::vector<size_t>
+    getFunctionOrdering(const MElement& element) const = 0;
 
   // Direct Access to Evaluated Functions //
   virtual void getFunctions(fullMatrix<double>& retValues,

FunctionSpace/BasisHierarchical0From.cpp

@@ -59,9 +59,14 @@ getOrientation(const MElement& element) const{
   return refSpace->getPermutation(element);
 }
 
-void BasisHierarchical0From::
-getFunctionPermutation(const MElement& element,
-                       unsigned int* indexPermutation) const{
+vector<size_t> BasisHierarchical0From::
+getFunctionOrdering(const MElement& element) const{
+  vector<size_t> ordering(nFunction);
+
+  for(size_t i = 0; i < nFunction; i++)
+    ordering[i] = i;
+
+  return ordering;
 }
 
 void BasisHierarchical0From::

FunctionSpace/BasisHierarchical0From.h

@@ -39,8 +39,8 @@ class BasisHierarchical0From: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
-  virtual void getFunctionPermutation(const MElement& element,
-                                      unsigned int* indexPermutation) const;
+  virtual std::vector<size_t>
+    getFunctionOrdering(const MElement& element) const;
 
   virtual void getFunctions(fullMatrix<double>& retValues,
                             const MElement& element,

FunctionSpace/BasisHierarchical1From.cpp

@@ -59,9 +59,14 @@ getOrientation(const MElement& element) const{
   return refSpace->getPermutation(element);
 }
 
-void BasisHierarchical1From::
-getFunctionPermutation(const MElement& element,
-                       unsigned int* indexPermutation) const{
+vector<size_t> BasisHierarchical1From::
+getFunctionOrdering(const MElement& element) const{
+  vector<size_t> ordering(nFunction);
+
+  for(size_t i = 0; i < nFunction; i++)
+    ordering[i] = i;
+
+  return ordering;
 }
 
 void BasisHierarchical1From::

FunctionSpace/BasisHierarchical1From.h

@@ -39,8 +39,8 @@ class BasisHierarchical1From: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
-  virtual void getFunctionPermutation(const MElement& element,
-                                      unsigned int* indexPermutation) const;
+  virtual std::vector<size_t>
+    getFunctionOrdering(const MElement& element) const;
 
   virtual void getFunctions(fullMatrix<double>& retValues,
                             const MElement& element,

FunctionSpace/BasisLagrange.cpp

-#include "Exception.h"
 #include "BasisLagrange.h"
 
 using namespace std;
 
 BasisLagrange::BasisLagrange(void){
   scalar = true;
+
+  preEvaluated     = false;
+  preEvaluatedGrad = false;
+
+  preEvaluatedFunction     = NULL;
+  preEvaluatedGradFunction = NULL;
 }
 
 BasisLagrange::~BasisLagrange(void){
+  if(preEvaluated)
+    delete preEvaluatedFunction;
+
+  if(preEvaluatedGrad)
+    delete preEvaluatedGradFunction;
 }
 
 unsigned int BasisLagrange::
@@ -20,9 +30,67 @@ getOrientation(const MElement& element) const{
   return 0;
 }
 
-void BasisLagrange::
-getFunctionPermutation(const MElement& element,
-                       unsigned int* indexPermutation) const{
+static bool
+sortPredicate(const std::pair<size_t, size_t>& a,
+              const std::pair<size_t, size_t>& b){
+  return a.second < b.second;
+}
+
+static vector<int> reducedNodeId(const MElement& element){
+  const size_t nVertex = element.getNumPrimaryVertices();
+  vector<pair<size_t, size_t> > vertexGlobalId(nVertex);
+
+  for(size_t i = 0; i < nVertex; i++){
+    vertexGlobalId[i].first  = i;
+    vertexGlobalId[i].second = element.getVertex(i)->getNum();
+  }
+
+  std::sort(vertexGlobalId.begin(), vertexGlobalId.end(), sortPredicate);
+
+  vector<int> vertexReducedId(nVertex);
+
+  for(size_t i = 0; i < nVertex; i++)
+    vertexReducedId[vertexGlobalId[i].first] = i;
+
+  return vertexReducedId;
+}
+
+static size_t matchClosure(vector<int>& reduced,
+                           nodalBasis::clCont& closures){
+
+  const size_t nNode = reduced.size();
+  const size_t nPerm = closures.size();
+
+  size_t i = 0;
+  bool   match = false;
+
+  while(i < nPerm && !match){
+    match = true;
+
+    for(size_t j = 0; j < nNode && match; j++)
+      if(reduced[j] != closures[i][j])
+         match = false;
+
+    if(!match)
+      i++;
+  }
+
+  return i;
+}
+
+vector<size_t> BasisLagrange::
+getFunctionOrdering(const MElement& element) const{
+  vector<int> rNodeId = reducedNodeId(element);
+  const size_t closureId = matchClosure(rNodeId, lBasis->fullClosures);
+
+  vector<int>& closure = lBasis->fullClosures[closureId];
+
+  vector<size_t> myClosure(closure.size());
+
+  for(size_t i = 0; i < closure.size(); i++)
+    myClosure[i] = closure[i];
+
+  return myClosure;
 }
 
 void BasisLagrange::
@@ -62,32 +130,55 @@ getFunctions(fullMatrix<double>& retValues,
 }
 
 void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
-  throw Exception("BasisLagrange::Not Implemented");
+  // Delete if older //
+  if(preEvaluated)
+    delete preEvaluatedFunction;
+
+  // Fill Matrix //
+  fullMatrix<double> tmp;
+  lBasis->f(point, tmp);
+
+  // Transpose 'tmp': otherwise not coherent with df !!
+  preEvaluatedFunction = new fullMatrix<double>(tmp.transpose());
+
+  // PreEvaluated //
+  preEvaluated = true;
 }
 
 void BasisLagrange::
 preEvaluateDerivatives(const fullMatrix<double>& point) const{
-  throw Exception("BasisLagrange::Not Implemented");
+  // Delete if older //
+  if(preEvaluatedGrad)
+    delete preEvaluatedGradFunction;
+
+  // Alloc //
+  preEvaluatedGradFunction = new fullMatrix<double>;
+
+  // Fill Matrix //
+  lBasis->df(point, *preEvaluatedGradFunction);
+
+  // PreEvaluated //
+  preEvaluatedGrad = true;
 }
 
 const fullMatrix<double>&
 BasisLagrange::getPreEvaluatedFunctions(const MElement& element) const{
-  throw Exception("BasisLagrange::Not Implemented");
+  return *preEvaluatedFunction;
 }
 
 const fullMatrix<double>&
 BasisLagrange::getPreEvaluatedDerivatives(const MElement& element) const{
-  throw Exception("BasisLagrange::Not Implemented");
+  return *preEvaluatedGradFunction;
 }
 
 const fullMatrix<double>&
 BasisLagrange::getPreEvaluatedFunctions(unsigned int orientation) const{
-  throw Exception("BasisLagrange::Not Implemented");
+  return *preEvaluatedFunction;
 }
 
 const fullMatrix<double>&
 BasisLagrange::getPreEvaluatedDerivatives(unsigned int orientation) const{
-  throw Exception("BasisLagrange::Not Implemented");
+  return *preEvaluatedGradFunction;
 }
 
 vector<double> BasisLagrange::

FunctionSpace/BasisLagrange.h

@@ -6,7 +6,7 @@
 #include "FunctionSpaceVector.h"
 #include "fullMatrix.h"
 #include "polynomialBasis.h"
-#include "ReferenceSpaceLagrange.h"
+// #include "ReferenceSpace.h"
 
 /**
    @interface BasisLagrange
@@ -28,7 +28,14 @@ class BasisLagrange: public BasisLocal{
  protected:
   polynomialBasis*        lBasis;   // Lagrange Basis
   fullMatrix<double>*     lPoint;   // Lagrange Points
-  ReferenceSpaceLagrange* refSpace; // Reference Space
+  // ReferenceSpace*         refSpace; // RefSpace
+
+  // PreEvaluation //
+  mutable bool preEvaluated;
+  mutable bool preEvaluatedGrad;
+
+  mutable fullMatrix<double>* preEvaluatedFunction;
+  mutable fullMatrix<double>* preEvaluatedGradFunction;
 
  public:
   virtual ~BasisLagrange(void);
@@ -36,8 +43,8 @@ class BasisLagrange: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
-  virtual void getFunctionPermutation(const MElement& element,
-                                      unsigned int* indexPermutation) const;
+  virtual std::vector<size_t>
+    getFunctionOrdering(const MElement& element) const;
 
   virtual void getFunctions(fullMatrix<double>& retValues,
                             const MElement& element,

FunctionSpace/FunctionSpace.cpp

@@ -108,13 +108,18 @@ void FunctionSpace::buildDof(void){
     vector<Dof> myDof = getKeys(*(element[i]));
     unsigned int nDof = myDof.size();
 
-    // Create new GroupOfDof
-    GroupOfDof* god = new GroupOfDof(nDof, *(element[i]));
-    (*group)[i]     = god;
-
     // Add Dof
+    vector<const Dof*> trueDof(nDof);
+
     for(unsigned int j = 0; j < nDof; j++)
-      insertDof(myDof[j], god);
+      insertDof(myDof[j], trueDof, j);
+
+    // Dof Orderning
+    vector<size_t> dofOrder = (*basis)[0]->getFunctionOrdering(*(element[i]));
+
+    // Create new GroupOfDof
+    GroupOfDof* god = new GroupOfDof(*(element[i]), trueDof, dofOrder);
+    (*group)[i]     = god;
 
     // Map GOD
     eToGod->insert(pair<const MElement*, const GroupOfDof*>
@@ -122,7 +127,9 @@ void FunctionSpace::buildDof(void){
   }
 }
 
-void FunctionSpace::insertDof(Dof& d, GroupOfDof* god){
+void FunctionSpace::insertDof(Dof& d,
+                              vector<const Dof*>& trueDof,
+                              size_t index){
   // Copy 'd'
   const Dof* tmp = new Dof(d);
 
@@ -133,13 +140,13 @@ void FunctionSpace::insertDof(Dof& d, GroupOfDof* god){
   // If insertion is OK (Dof 'd' didn't exist) //
   //   --> Add new Dof in GoD
   if(p.second)
-    god->add(*tmp);
+    trueDof[index] = tmp;
 
   // If insertion failed (Dof 'd' already exists) //
-  //   --> delete 'd' and add existing Dof in GoD
+  //   --> delete 'tmp' and add existing Dof in GoD
   else{
     delete tmp;
-    god->add(*(*(p.first)));
+    trueDof[index] = *(p.first);
   }
 }
 

FunctionSpace/FunctionSpace.h

@@ -93,11 +93,13 @@ class FunctionSpace{
   FunctionSpace(void);
 
   void build(GroupOfElement& goe,
-	     const Basis& basis);
+             const Basis& basis);
 
   // Dof
   void buildDof(void);
-  void insertDof(Dof& d, GroupOfDof* god);
+  void insertDof(Dof& d,
+                 std::vector<const Dof*>& trueDof,
+                 size_t index);
 };
 
 

FunctionSpace/TriLagrangeBasis.cpp

-#include "Exception.h"
 #include "TriLagrangeBasis.h"
 #include "pointsGenerators.h"
-#include "TriLagrangeReferenceSpace.h"
+//#include "TriReferenceSpace.h"
 #include "ElementType.h"
 
-#include <iostream>
-
 TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
   // If order 0 (Nedelec): use order 1
   if(order == 0)
@@ -24,13 +21,14 @@ TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
   nFunction = nVertex + nEdge + nFace + nCell;
 
   // Init polynomialBasis //
-  lBasis = new polynomialBasis(getTag(order));
+  lBasis = new polynomialBasis(ElementType::getTag(TYPE_TRI, order, false));
 
   // Init Lagrange Point //
   lPoint = new fullMatrix<double>
     (gmshGeneratePointsTriangle(order, false));
 
   // ReferenceSpace //
+  // refSpace = new TriReferenceSpace;
   //  refSpace = new TriLagrangeReferenceSpace(order);
   //  std::cout << refSpace->toString() << std::endl;
 }
@@ -40,15 +38,3 @@ TriLagrangeBasis::~TriLagrangeBasis(void){
   delete lPoint;
   // delete refSpace;
 }
-
-unsigned int TriLagrangeBasis::getTag(unsigned int order){
-  unsigned int tag = ElementType::getTag(TYPE_TRI, order, false);
-
-  if(tag)
-    return tag;
-
-  else
-    throw Exception
-      ("Can't instanciate an order %d Lagrangian Basis for a Triangle",
-       order);
-}