diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index aa5901dd625807c3c5148f71c092954beb1934ff..dbb21ddb423b5e04e8b85af50b3dde0fe51d80a2 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -20,121 +20,6 @@ BasisLagrange::~BasisLagrange(void){
delete preEvaluatedGradFunction;
}
-/*
-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;
-}
-
-void BasisLagrange::mapFromXYZtoABC(const MElement& element,
- const fullVector<double>& xyz,
- double abc[3]) const{
-}
-
-vector<size_t> BasisLagrange::
-getFunctionOrdering(const MElement& element) const{
-
- static bool once = true;
-
- if(once){
- for(size_t i = 0; i < lBasis->fullClosures.size(); i++){
- vector<int>& closure = lBasis->fullClosures[i];
-
- for(size_t j =0; j < closure.size(); j++)
- cout << closure[j] << "\t";
- cout << endl;
- }
-
- once = false;
- }
-
-
- 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;
-
- vector<size_t> c(10);
-
- if(closureId == 0){
- c[0] = 0;
- c[1] = 1;
- c[2] = 2;
- c[3] = 3;
- c[4] = 4;
- c[5] = 5;
- c[6] = 6;
- c[7] = 8; // 7;
- c[8] = 7; // 8;
- c[9] = 9;
- }
-
- else{
- c[0] = 2;
- c[1] = 0;
- c[2] = 1;
- c[3] = 8; // 7;
- c[4] = 7; // 8;
- c[5] = 3;
- c[6] = 4;
- c[7] = 5;
- c[8] = 6;
- c[9] = 9;
- }
-
- return c;
-}
-*/
void BasisLagrange::
getFunctions(fullMatrix<double>& retValues,
const MElement& element,
@@ -151,6 +36,9 @@ getFunctions(fullMatrix<double>& retValues,
// Transpose 'tmp': otherwise not coherent with df !!
retValues = tmp.transpose();
+
+ // Permute retValues, accordingly to ReferenceSpace
+ permutation(refSpace->getReferenceSpace(element), retValues);
}
void BasisLagrange::
@@ -169,6 +57,9 @@ getFunctions(fullMatrix<double>& retValues,
// Transpose 'tmp': otherwise not coherent with df !!
retValues = tmp.transpose();
+
+ // Permute retValues, accordingly to ReferenceSpace
+ permutation(orientation, retValues);
}
void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
@@ -299,6 +190,10 @@ project(const MElement& element,
return newCoef;
}
+void BasisLagrange::permutation(size_t orientation,
+ fullMatrix<double>& function) const{
+}
+
std::string BasisLagrange::toString(void) const{
return std::string("BasisLagrange::toString() not Implemented");
}
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index 8abd5afe3ba56229bf52b2f121847bb2a894f3ba..216460ab3f00a05c5b76b994c5afbeb9980b9a02 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/FunctionSpace/BasisLagrange.h
@@ -75,6 +75,10 @@ class BasisLagrange: public BasisLocal{
protected:
BasisLagrange(void);
+
+ private:
+ void permutation(size_t orientation,
+ fullMatrix<double>& function) const;
};
diff --git a/FunctionSpace/HexLagrangeBasis.cpp b/FunctionSpace/HexLagrangeBasis.cpp
index 29c23a9330bc21d40d952bb3e7b3011a6c175890..7fc92e2ceed2c1f2cc3cd550caa293888a307395 100644
--- a/FunctionSpace/HexLagrangeBasis.cpp
+++ b/FunctionSpace/HexLagrangeBasis.cpp
@@ -1,4 +1,5 @@
#include "HexLagrangeBasis.h"
+#include "HexReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h"
@@ -24,9 +25,14 @@ HexLagrangeBasis::HexLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsHexahedron(order, false));
+
+ // Reference Space //
+ refSpace = new HexReferenceSpace;
+ nRefSpace = getReferenceSpace().getNReferenceSpace();
}
HexLagrangeBasis::~HexLagrangeBasis(void){
delete lBasis;
delete lPoint;
+ delete refSpace;
}
diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp
index 9dc6085605ab3a0e45fb84b8ffa3e171c4cf7319..f3ff36015b8d56e4e327bed96f3c87581cbc645a 100644
--- a/FunctionSpace/LineEdgeBasis.cpp
+++ b/FunctionSpace/LineEdgeBasis.cpp
@@ -24,27 +24,21 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
- // Alloc Temporary Space //
+ // Legendre Polynomial //
const size_t orderPlus = order + 1;
Polynomial* intLegendre = new Polynomial[orderPlus];
- vector<Polynomial> first(3);
- first[0] = Polynomial(-0.5, 0, 0, 0);
- first[1] = Polynomial( 0 , 0, 0, 0);
- first[2] = Polynomial( 0 , 0, 0, 0);
-
- vector<Polynomial> second(3);
- second[0] = Polynomial(+0.5, 0, 0, 0);
- second[1] = Polynomial( 0 , 0, 0, 0);
- second[2] = Polynomial( 0 , 0, 0, 0);
+ Legendre::integrated(intLegendre, orderPlus);
- const Polynomial x[2] = {
- Polynomial(-1, 1, 0, 0),
- Polynomial(+1, 1, 0, 0)
- };
+ // Lagrange Polynomial //
+ const Polynomial lagrange[2] =
+ {
+ Polynomial(Polynomial(0.5, 0, 0, 0) -
+ Polynomial(0.5, 1, 0, 0)),
- // Legendre Polynomial //
- Legendre::integrated(intLegendre, orderPlus);
+ Polynomial(Polynomial(0.5, 0, 0, 0) +
+ Polynomial(0.5, 1, 0, 0)),
+ };
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
@@ -52,19 +46,38 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
- // Edge Based (Nedelec) //
- basis[0][0] = new vector<Polynomial>(first);
- basis[1][0] = new vector<Polynomial>(second);
-
- // Edge Based (High Order) //
+ // Edge Based //
for(size_t s = 0; s < nRefSpace; s++){
- size_t i = 1;
-
- for(size_t l = 1; l < orderPlus; l++){
- basis[s][i] =
- new vector<Polynomial>((intLegendre[l].compose
- (x[edgeIdx[s][0][0]])).gradient());
-
+ size_t i = 0;
+
+ for(size_t l = 0; l < orderPlus; l++){
+ // Nedelec
+ if(l == 0){
+ vector<Polynomial> tmp1 = lagrange[edgeIdx[s][0][1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[edgeIdx[s][0][0]].gradient();
+
+ tmp1[0].mul(lagrange[edgeIdx[s][0][0]]);
+ tmp1[1].mul(lagrange[edgeIdx[s][0][0]]);
+ tmp1[2].mul(lagrange[edgeIdx[s][0][0]]);
+
+ tmp2[0].mul(lagrange[edgeIdx[s][0][1]]);
+ tmp2[1].mul(lagrange[edgeIdx[s][0][1]]);
+ tmp2[2].mul(lagrange[edgeIdx[s][0][1]]);
+
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
+
+ basis[s][i] = new vector<Polynomial>(tmp2);
+ }
+
+ // High Order
+ else{
+ basis[s][i] =
+ new vector<Polynomial>
+ ((intLegendre[l].compose(lagrange[edgeIdx[s][0][0]] -
+ lagrange[edgeIdx[s][0][1]]).gradient()));
+ }
i++;
}
}
diff --git a/FunctionSpace/LineLagrangeBasis.cpp b/FunctionSpace/LineLagrangeBasis.cpp
index 9a17e5034923e447bc7db9bc287c448c45971c0b..6842133c20b99439eaf7c7f4758f70cc0494237b 100644
--- a/FunctionSpace/LineLagrangeBasis.cpp
+++ b/FunctionSpace/LineLagrangeBasis.cpp
@@ -1,4 +1,5 @@
#include "LineLagrangeBasis.h"
+#include "LineReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h"
@@ -24,9 +25,14 @@ LineLagrangeBasis::LineLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsLine(order));
+
+ // Reference Space //
+ refSpace = new LineReferenceSpace;
+ nRefSpace = getReferenceSpace().getNReferenceSpace();
}
LineLagrangeBasis::~LineLagrangeBasis(void){
delete lBasis;
delete lPoint;
+ delete refSpace;
}
diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp
index 423ff30cc9c478b45e12c7699d0f5832f83319fe..ee1850ff1ccafa0f31d68c21312f38f1c4543ed8 100644
--- a/FunctionSpace/LineNedelecBasis.cpp
+++ b/FunctionSpace/LineNedelecBasis.cpp
@@ -1,6 +1,5 @@
#include "LineNedelecBasis.h"
#include "LineReferenceSpace.h"
-#include "Legendre.h"
using namespace std;
@@ -9,6 +8,9 @@ LineNedelecBasis::LineNedelecBasis(void){
refSpace = new LineReferenceSpace;
nRefSpace = getReferenceSpace().getNReferenceSpace();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
+
// Set Basis Type //
order = 0;
@@ -21,16 +23,15 @@ LineNedelecBasis::LineNedelecBasis(void){
nCell = 0;
nFunction = 1;
- // Alloc Temporary Space //
- vector<Polynomial> first(3);
- first[0] = Polynomial(-0.5, 0, 0, 0);
- first[1] = Polynomial( 0 , 0, 0, 0);
- first[2] = Polynomial( 0 , 0, 0, 0);
+ // Lagrange Polynomial //
+ const Polynomial lagrange[2] =
+ {
+ Polynomial(Polynomial(0.5, 0, 0, 0) -
+ Polynomial(0.5, 1, 0, 0)),
- vector<Polynomial> second(3);
- second[0] = Polynomial(+0.5, 0, 0, 0);
- second[1] = Polynomial( 0 , 0, 0, 0);
- second[2] = Polynomial( 0 , 0, 0, 0);
+ Polynomial(Polynomial(0.5, 0, 0, 0) +
+ Polynomial(0.5, 1, 0, 0)),
+ };
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
@@ -38,9 +39,25 @@ LineNedelecBasis::LineNedelecBasis(void){
for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
- // Nedelec //
- basis[0][0] = new vector<Polynomial>(first);
- basis[1][0] = new vector<Polynomial>(second);
+ // Edge Based (Nedelec) //
+ for(size_t s = 0; s < nRefSpace; s++){
+ vector<Polynomial> tmp1 = lagrange[edgeIdx[s][0][1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[edgeIdx[s][0][0]].gradient();
+
+ tmp1[0].mul(lagrange[edgeIdx[s][0][0]]);
+ tmp1[1].mul(lagrange[edgeIdx[s][0][0]]);
+ tmp1[2].mul(lagrange[edgeIdx[s][0][0]]);
+
+ tmp2[0].mul(lagrange[edgeIdx[s][0][1]]);
+ tmp2[1].mul(lagrange[edgeIdx[s][0][1]]);
+ tmp2[2].mul(lagrange[edgeIdx[s][0][1]]);
+
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
+
+ basis[s][0] = new vector<Polynomial>(tmp2);
+ }
}
LineNedelecBasis::~LineNedelecBasis(void){
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index feab18293dc63ea3cc32185114701b3dfe57483e..d4e9d5cb9915b935dc45cc474e0989aeb0aaac3c 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -24,16 +24,19 @@ LineNodeBasis::LineNodeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
- // Alloc Temporary Space //
+ // Legendre Polynomial //
Polynomial* intLegendre = new Polynomial[order];
+ Legendre::integrated(intLegendre, order);
- const Polynomial x[2] = {
- Polynomial(+1, 1, 0, 0),
- Polynomial(-1, 1, 0, 0)
- };
+ // Lagrange Polynomial //
+ const Polynomial lagrange[2] =
+ {
+ Polynomial(Polynomial(0.5, 0, 0, 0) -
+ Polynomial(0.5, 1, 0, 0)),
- // Legendre Polynomial //
- Legendre::integrated(intLegendre, order);
+ Polynomial(Polynomial(0.5, 0, 0, 0) +
+ Polynomial(0.5, 1, 0, 0)),
+ };
// Basis //
basis = new Polynomial**[nRefSpace];
@@ -41,15 +44,10 @@ LineNodeBasis::LineNodeBasis(size_t order){
for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new Polynomial*[nFunction];
- // Vertex Based (Lagrange) //
+ // Vertex Based //
for(size_t s = 0; s < nRefSpace; s++){
- basis[s][0] =
- new Polynomial(Polynomial(0.5, 0, 0, 0) -
- Polynomial(0.5, 1, 0, 0));
-
- basis[s][1] =
- new Polynomial(Polynomial(0.5, 0, 0, 0) +
- Polynomial(0.5, 1, 0, 0));
+ basis[s][0] = new Polynomial(lagrange[0]);
+ basis[s][1] = new Polynomial(lagrange[1]);
}
// Edge Based //
@@ -58,7 +56,8 @@ LineNodeBasis::LineNodeBasis(size_t order){
for(size_t l = 1; l < order; l++){
basis[s][i] =
- new Polynomial(intLegendre[l].compose(x[edgeIdx[s][0][0]]));
+ new Polynomial(intLegendre[l].compose(lagrange[edgeIdx[s][0][1]] -
+ lagrange[edgeIdx[s][0][0]]));
i++;
}
diff --git a/FunctionSpace/QuadLagrangeBasis.cpp b/FunctionSpace/QuadLagrangeBasis.cpp
index 9367c4f7273a8e50e0c2d6cce2c7411fe2164587..31efcf5555f255b9e2b78642c7b26ac531bba467 100644
--- a/FunctionSpace/QuadLagrangeBasis.cpp
+++ b/FunctionSpace/QuadLagrangeBasis.cpp
@@ -1,4 +1,5 @@
#include "QuadLagrangeBasis.h"
+#include "QuadReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h"
@@ -24,9 +25,14 @@ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsQuadrangle(order, false));
+
+ // Reference Space //
+ refSpace = new QuadReferenceSpace;
+ nRefSpace = getReferenceSpace().getNReferenceSpace();
}
QuadLagrangeBasis::~QuadLagrangeBasis(void){
delete lBasis;
delete lPoint;
+ delete refSpace;
}
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index 7898a4caac4b15cf484c52819bf997e930050b97..07d77b10aeee2663d0011d92150484b3d186583c 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -202,7 +202,7 @@ isCyclicPermutation(const vector<size_t>& pTest,
// Triplet to return
triplet tri;
-
+ /*
// If no Face, we have no Cyclic Permutation
if(!refFaceNodeIdx.size()){
tri.first = false;
@@ -240,11 +240,11 @@ isCyclicPermutation(const vector<size_t>& pTest,
// Test if we have a face cyclic permutation
size_t isCyclic = isFacePermutation(refNode, testNode);
-
- // Test ifwe have the same connectivity
+ */
+ // Test if we have the same connectivity
bool isSameConnectivity = isSameEdge(pTest, pRef);
- if(isCyclic && isSameConnectivity){
+ if(isSameConnectivity){
tri.first = true;
tri.second = getRefIndexPermutation(pRef, pTest);
tri.third = getReverseIndexPermutation(pRef, pTest);
@@ -258,7 +258,7 @@ isCyclicPermutation(const vector<size_t>& pTest,
return tri;
}
-
+/*
size_t ReferenceSpace::findCorrespondingFace(const vector<size_t>& face,
const vector<size_t>& node) const{
// Init Stuff //
@@ -347,7 +347,7 @@ bool ReferenceSpace::isFacePermutation(const vector<size_t>& refNode,
return match;
}
-
+*/
bool ReferenceSpace::isSameEdge(const std::vector<size_t>& pTest,
const std::vector<size_t>& pRef) const{
// Set of Reference Edges
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index 2e31620d439b9f634162b6fe663ac9e5fd58fedb..42856ec90d10f817f73554df6f481de6ac124ccc 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -148,7 +148,7 @@ class ReferenceSpace{
triplet isCyclicPermutation(const std::vector<size_t>& pTest,
const std::vector<size_t>& pRef) const;
-
+ /*
size_t findCorrespondingFace(const std::vector<size_t>& face,
const std::vector<size_t>& node) const;
@@ -157,7 +157,7 @@ class ReferenceSpace{
static bool isFacePermutation(const std::vector<size_t>& refNode,
const std::vector<size_t>& testNode);
-
+ */
bool isSameEdge(const std::vector<size_t>& pTest,
const std::vector<size_t>& pRef) const;
diff --git a/FunctionSpace/TetLagrangeBasis.cpp b/FunctionSpace/TetLagrangeBasis.cpp
index 3aaf26ca793ca0414b079f949d1d71d559042c49..527c6c33db9df391db9582ba4594d0950dff47b7 100644
--- a/FunctionSpace/TetLagrangeBasis.cpp
+++ b/FunctionSpace/TetLagrangeBasis.cpp
@@ -1,4 +1,5 @@
#include "TetLagrangeBasis.h"
+#include "TetReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h"
@@ -24,9 +25,14 @@ TetLagrangeBasis::TetLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsTetrahedron(order, false));
+
+ // Reference Space //
+ refSpace = new TetReferenceSpace;
+ nRefSpace = getReferenceSpace().getNReferenceSpace();
}
TetLagrangeBasis::~TetLagrangeBasis(void){
delete lBasis;
delete lPoint;
+ delete refSpace;
}
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 6a4a052731cdff18510dcb8d59a09112fc3fda93..e95fcdd1ae489efe170ef28fd3b53a7520b502c4 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -1,4 +1,5 @@
#include "TriLagrangeBasis.h"
+#include "TriReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h"
@@ -24,9 +25,14 @@ TriLagrangeBasis::TriLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsTriangle(order, false));
+
+ // Reference Space //
+ refSpace = new TriReferenceSpace;
+ nRefSpace = getReferenceSpace().getNReferenceSpace();
}
TriLagrangeBasis::~TriLagrangeBasis(void){
delete lBasis;
delete lPoint;
+ delete refSpace;
}