diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index e1e515091f34950e9f7e8e1de80799fa579c94d8..7738a321f36d6fdfb130ab7d3b935fc71834e44d 100644
--- a/FunctionSpace/Basis.h
+++ b/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,
diff --git a/FunctionSpace/BasisHierarchical0From.cpp b/FunctionSpace/BasisHierarchical0From.cpp
index 36778207b728b765f595b9a21e5931742732c55d..ec7848c2700973a071b4ac036d64fba093928e3f 100644
--- a/FunctionSpace/BasisHierarchical0From.cpp
+++ b/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::
diff --git a/FunctionSpace/BasisHierarchical0From.h b/FunctionSpace/BasisHierarchical0From.h
index 15c180f28b0690dd5ec0ba58b1de0a5c78e84ced..97d0cba28c972e1e42a21cac86d4deb4b6db8582 100644
--- a/FunctionSpace/BasisHierarchical0From.h
+++ b/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,
diff --git a/FunctionSpace/BasisHierarchical1From.cpp b/FunctionSpace/BasisHierarchical1From.cpp
index 826d4b64c9dac718f4f1e0425ffb4b59559cfccb..e26d493118c5996e7fd5f5a113b6e240095ce6cd 100644
--- a/FunctionSpace/BasisHierarchical1From.cpp
+++ b/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::
diff --git a/FunctionSpace/BasisHierarchical1From.h b/FunctionSpace/BasisHierarchical1From.h
index ad0f99e352086b8f614323cafcedd306d9de0c15..305783794902a474d1b86a6d89e581520e7021c5 100644
--- a/FunctionSpace/BasisHierarchical1From.h
+++ b/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,
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index 36a0d8732972fc0b7dcf49adc12115838c8dcb78..77c3077ebb65d4b8d6ac91d982d7cfe6d1eea1f6 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -1,13 +1,23 @@
-#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::
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index 105802265c1436a7f55fe4c922135257922c7adf..b82487fc4157c5d4acb5fc3edffd8a67d5910785 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/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,
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index ae99c457bb0b38dc4f0043e185763dbcee15c96e..dd0b32a3eb755833ea923ec2a62d4dda237dd9a3 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/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);
}
}
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 6f50b3f0f24e87a71c41a0b2050b286e10e3fa02..2c9881eda1bd1b6730b6552dd142468c5bb5459c 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/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);
};
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 263ffd397619272ee2979b7d7cf08474a857ffed..a6d4073b739984fba46efa51061dfc4ff47a2550 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -1,11 +1,8 @@
-#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);
-}