From 2c3ef322b22b9190acaa16e705d8e3608dd482a1 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Fri, 5 Jul 2013 13:35:38 +0000
Subject: [PATCH] Reworking ReferenceSpace -- Trying to have minimum reference
space -- find permutation and reference space for 2d and tet
---
FunctionSpace/BasisHierarchical0From.cpp | 10 +-
FunctionSpace/BasisHierarchical1From.cpp | 10 +-
FunctionSpace/BasisLagrange.cpp | 48 ++-
FunctionSpace/CMakeLists.txt | 2 +
FunctionSpace/FunctionSpace.cpp | 12 +-
FunctionSpace/FunctionSpaceScalar.cpp | 66 +++
FunctionSpace/FunctionSpaceScalar.h | 8 +-
FunctionSpace/FunctionSpaceVector.h | 8 +-
FunctionSpace/LineEdgeBasis.cpp | 6 +-
FunctionSpace/LineNedelecBasis.cpp | 2 +-
FunctionSpace/LineNodeBasis.cpp | 8 +-
FunctionSpace/LineReferenceSpace.cpp | 41 +-
FunctionSpace/PermutationTree.cpp | 174 ++++++++
FunctionSpace/PermutationTree.h | 83 ++++
FunctionSpace/QuadEdgeBasis.cpp | 8 +-
FunctionSpace/QuadNedelecBasis.cpp | 6 +-
FunctionSpace/QuadNodeBasis.cpp | 14 +-
FunctionSpace/QuadReferenceSpace.cpp | 34 +-
FunctionSpace/ReferenceSpace.cpp | 423 +++++++++++--------
FunctionSpace/ReferenceSpace.h | 87 ++--
FunctionSpace/ReferenceSpaceLagrange.cpp | 11 +-
FunctionSpace/ReferenceSpaceLagrange.h | 4 +-
FunctionSpace/TetEdgeBasis.cpp | 446 ++++++++++----------
FunctionSpace/TetNedelecBasis.cpp | 8 +-
FunctionSpace/TetNodeBasis.cpp | 88 ++--
FunctionSpace/TetReferenceSpace.cpp | 62 +--
FunctionSpace/TriEdgeBasis.cpp | 122 +++---
FunctionSpace/TriLagrangeReferenceSpace.cpp | 14 +-
FunctionSpace/TriNedelecBasis.cpp | 6 +-
FunctionSpace/TriNodeBasis.cpp | 32 +-
FunctionSpace/TriReferenceSpace.cpp | 32 +-
31 files changed, 1164 insertions(+), 711 deletions(-)
create mode 100644 FunctionSpace/PermutationTree.cpp
create mode 100644 FunctionSpace/PermutationTree.h
diff --git a/FunctionSpace/BasisHierarchical0From.cpp b/FunctionSpace/BasisHierarchical0From.cpp
index ec7848c270..26fd0ce1e0 100644
--- a/FunctionSpace/BasisHierarchical0From.cpp
+++ b/FunctionSpace/BasisHierarchical0From.cpp
@@ -51,12 +51,12 @@ BasisHierarchical0From::~BasisHierarchical0From(void){
unsigned int BasisHierarchical0From::
getNOrientation(void) const{
- return refSpace->getNPermutation();
+ return refSpace->getNReferenceSpace();
}
unsigned int BasisHierarchical0From::
getOrientation(const MElement& element) const{
- return refSpace->getPermutation(element);
+ return refSpace->getReferenceSpace(element);
}
vector<size_t> BasisHierarchical0From::
@@ -75,7 +75,7 @@ getFunctions(fullMatrix<double>& retValues,
double u, double v, double w) const{
// Define Orientation //
- unsigned int orientation = refSpace->getPermutation(element);
+ unsigned int orientation = refSpace->getReferenceSpace(element);
// Fill Matrix //
for(unsigned int i = 0; i < nFunction; i++)
@@ -170,12 +170,12 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const fullMatrix<double>& BasisHierarchical0From::
getPreEvaluatedFunctions(const MElement& element) const{
- return getPreEvaluatedFunctions(refSpace->getPermutation(element));
+ return getPreEvaluatedFunctions(refSpace->getReferenceSpace(element));
}
const fullMatrix<double>& BasisHierarchical0From::
getPreEvaluatedDerivatives(const MElement& element) const{
- return getPreEvaluatedDerivatives(refSpace->getPermutation(element));
+ return getPreEvaluatedDerivatives(refSpace->getReferenceSpace(element));
}
const fullMatrix<double>& BasisHierarchical0From::
diff --git a/FunctionSpace/BasisHierarchical1From.cpp b/FunctionSpace/BasisHierarchical1From.cpp
index e26d493118..7fe070a919 100644
--- a/FunctionSpace/BasisHierarchical1From.cpp
+++ b/FunctionSpace/BasisHierarchical1From.cpp
@@ -51,12 +51,12 @@ BasisHierarchical1From::~BasisHierarchical1From(void){
unsigned int BasisHierarchical1From::
getNOrientation(void) const{
- return refSpace->getNPermutation();
+ return refSpace->getNReferenceSpace();
}
unsigned int BasisHierarchical1From::
getOrientation(const MElement& element) const{
- return refSpace->getPermutation(element);
+ return refSpace->getReferenceSpace(element);
}
vector<size_t> BasisHierarchical1From::
@@ -75,7 +75,7 @@ getFunctions(fullMatrix<double>& retValues,
double u, double v, double w) const{
// Define Orientation //
- unsigned int orientation = refSpace->getPermutation(element);
+ unsigned int orientation = refSpace->getReferenceSpace(element);
// Fill Vector //
for(unsigned int i = 0; i < nFunction; i++){
@@ -192,12 +192,12 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const fullMatrix<double>& BasisHierarchical1From::
getPreEvaluatedFunctions(const MElement& element) const{
- return getPreEvaluatedFunctions(refSpace->getPermutation(element));
+ return getPreEvaluatedFunctions(refSpace->getReferenceSpace(element));
}
const fullMatrix<double>& BasisHierarchical1From::
getPreEvaluatedDerivatives(const MElement& element) const{
- return getPreEvaluatedDerivatives(refSpace->getPermutation(element));
+ return getPreEvaluatedDerivatives(refSpace->getReferenceSpace(element));
}
const fullMatrix<double>& BasisHierarchical1From::
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index 77c3077ebb..dcdb1e85e2 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -80,9 +80,25 @@ static size_t matchClosure(vector<int>& reduced,
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());
@@ -91,6 +107,36 @@ getFunctionOrdering(const MElement& element) const{
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::
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index 4b8ac0a8ed..e59333de60 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -7,6 +7,8 @@ set(SRC
Polynomial.cpp
Legendre.cpp
+ PermutationTree.cpp
+
ReferenceSpace.cpp
LineReferenceSpace.cpp
TriReferenceSpace.cpp
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index dd0b32a3eb..bc3488878c 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -116,9 +116,17 @@ void FunctionSpace::buildDof(void){
// Dof Orderning
vector<size_t> dofOrder = (*basis)[0]->getFunctionOrdering(*(element[i]));
-
+ vector<const Dof*> trueDofOrdered(nDof);
+
+ for(size_t j = 0; j < nDof; j++)
+ trueDofOrdered[dofOrder[j]] = trueDof[j];
+ /*
+ for(size_t j = 0; j < nDof; j++)
+ cout << dofOrder[j] << "\t";
+ cout << endl;
+ */
// Create new GroupOfDof
- GroupOfDof* god = new GroupOfDof(*(element[i]), trueDof, dofOrder);
+ GroupOfDof* god = new GroupOfDof(*(element[i]), trueDofOrdered);
(*group)[i] = god;
// Map GOD
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index fa675e8c97..1400403da7 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -22,9 +22,75 @@ interpolate(const MElement& element,
// Get Reference coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
double uvw[3];
+ double abc[3];
eelement.xyz2uvw(phys, uvw);
+ // Compute coordinate in abc Space //
+ std::vector<size_t> order = (*basis)[0]->getFunctionOrdering(element);
+ double abcMat[3][3];
+ double abcPnt[3][3];
+
+ element.getNode(0, abcPnt[0][0], abcPnt[0][1], abcPnt[0][2]);
+ element.getNode(1, abcPnt[1][0], abcPnt[1][1], abcPnt[1][2]);
+ element.getNode(2, abcPnt[2][0], abcPnt[2][1], abcPnt[2][2]);
+ /*
+ std::cout << order[0] << " | "
+ << order[1] << " | "
+ << order[2] << std::endl;
+ */
+ abcMat[0][0] = abcPnt[0][0];
+ abcMat[0][1] = abcPnt[1][0];
+ abcMat[0][2] = abcPnt[2][0];
+
+ abcMat[1][0] = abcPnt[0][1];
+ abcMat[1][1] = abcPnt[1][1];
+ abcMat[1][2] = abcPnt[2][1];
+
+ abcMat[2][0] = abcPnt[0][2];
+ abcMat[2][1] = abcPnt[1][2];
+ abcMat[2][2] = abcPnt[2][2];
+ /*
+ for(size_t i = 0; i < 3; i++){
+ for(size_t j = 0; j < 3; j++)
+ std::cout << abcMat[i][j] << "\t";
+ std::cout << std::endl;
+ }
+ std::cout << std::endl;
+ */
+ double phiUVW[3];
+ element.getShapeFunctions(uvw[0], uvw[1], uvw[2], phiUVW);
+
+ size_t trueOrder[3];
+
+ for(size_t i = 0; i < 3; i++){
+ size_t idx;
+
+ for(idx = 0; order[idx] != i; idx++)
+ ;
+
+ trueOrder[i] = idx;
+ }
+
+ abc[0] =
+ abcMat[0][0] * phiUVW[trueOrder[0]] +
+ abcMat[0][1] * phiUVW[trueOrder[1]] +
+ abcMat[0][2] * phiUVW[trueOrder[2]];
+
+ abc[1] =
+ abcMat[1][0] * phiUVW[trueOrder[0]] +
+ abcMat[1][1] * phiUVW[trueOrder[1]] +
+ abcMat[1][2] * phiUVW[trueOrder[2]];
+
+ abc[2] =
+ abcMat[2][0] * phiUVW[trueOrder[0]] +
+ abcMat[2][1] * phiUVW[trueOrder[1]] +
+ abcMat[2][2] * phiUVW[trueOrder[2]];
+
+ uvw[0] = abc[0];
+ uvw[1] = abc[1];
+ uvw[2] = abc[2];
+
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
fullMatrix<double> fun(nFun, 1);
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index ec9979570b..56350ae363 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -23,13 +23,13 @@ class FunctionSpaceScalar : public FunctionSpace{
double
interpolate(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& xyz) const;
+ const std::vector<double>& coef,
+ const fullVector<double>& xyz) const;
double
interpolateInRefSpace(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& uvw) const;
+ const std::vector<double>& coef,
+ const fullVector<double>& uvw) const;
};
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index df736855bb..9428e6ffe3 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -24,13 +24,13 @@ class FunctionSpaceVector : public FunctionSpace{
fullVector<double>
interpolate(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& xyz) const;
+ const std::vector<double>& coef,
+ const fullVector<double>& xyz) const;
fullVector<double>
interpolateInRefSpace(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& uvw) const;
+ const std::vector<double>& coef,
+ const fullVector<double>& uvw) const;
};
diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp
index 2fa03098b3..933d162ff1 100644
--- a/FunctionSpace/LineEdgeBasis.cpp
+++ b/FunctionSpace/LineEdgeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
LineEdgeBasis::LineEdgeBasis(unsigned int order){
// Reference Space //
refSpace = new LineReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -62,8 +62,8 @@ LineEdgeBasis::LineEdgeBasis(unsigned int order){
for(unsigned int l = 1; l < orderPlus; l++){
basis[s][i] =
- new vector<Polynomial>((intLegendre[l].compose
- (x[(*(*edgeV[s])[0])[0]])).gradient());
+ new vector<Polynomial>((intLegendre[l].compose
+ (x[(*(*edgeV[s])[0])[0]])).gradient());
i++;
}
diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp
index 3a9e9cff20..0e3727d7ef 100644
--- a/FunctionSpace/LineNedelecBasis.cpp
+++ b/FunctionSpace/LineNedelecBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
LineNedelecBasis::LineNedelecBasis(void){
// Reference Space //
refSpace = new LineReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
// Set Basis Type //
order = 0;
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index 3839f9b0fc..3c00f46ba1 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
LineNodeBasis::LineNodeBasis(unsigned int order){
// Reference Space //
refSpace = new LineReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -45,11 +45,11 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
for(unsigned int s = 0; s < nRefSpace; s++){
basis[s][0] =
new Polynomial(Polynomial(0.5, 0, 0, 0) -
- Polynomial(0.5, 1, 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));
+ Polynomial(0.5, 1, 0, 0));
}
// Edge Based //
@@ -58,7 +58,7 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
for(unsigned int l = 1; l < order; l++){
basis[s][i] =
- new Polynomial(intLegendre[l].compose(x[(*(*edgeV[s])[0])[0]]));
+ new Polynomial(intLegendre[l].compose(x[(*(*edgeV[s])[0])[0]]));
i++;
}
diff --git a/FunctionSpace/LineReferenceSpace.cpp b/FunctionSpace/LineReferenceSpace.cpp
index 7ba0eb2d4a..9251fdec6e 100644
--- a/FunctionSpace/LineReferenceSpace.cpp
+++ b/FunctionSpace/LineReferenceSpace.cpp
@@ -9,8 +9,8 @@ LineReferenceSpace::LineReferenceSpace(void){
// Edge Definition //
nEdge = 1;
- refEdge = new unsigned int*[nEdge];
- refEdge[0] = new unsigned int[2];
+ refEdge = new size_t*[nEdge];
+ refEdge[0] = new size_t[2];
refEdge[0][0] = 0;
refEdge[0][1] = 1;
@@ -25,44 +25,45 @@ LineReferenceSpace::LineReferenceSpace(void){
LineReferenceSpace::~LineReferenceSpace(void){
// Delete Ref Edge //
- for(unsigned int i = 0; i < nEdge; i++)
+ for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
delete[] refEdge;
}
string LineReferenceSpace::toLatex(void) const{
+ const size_t nPerm = pTree->getNPermutation();
stringstream stream;
stream << "\\documentclass{article}" << endl << endl
- << "\\usepackage{longtable}" << endl
- << "\\usepackage{tikz}" << endl
- << "\\usetikzlibrary{arrows}" << endl << endl
+ << "\\usepackage{longtable}" << endl
+ << "\\usepackage{tikz}" << endl
+ << "\\usetikzlibrary{arrows}" << endl << endl
- << "\\begin{document}" << endl
- << "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
- << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
+ << "\\begin{document}" << endl
+ << "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
+ << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
- << "\\begin{longtable}{c}" << endl << endl;
+ << "\\begin{longtable}{c}" << endl << endl;
- for(unsigned int p = 0; p < nPerm; p++){
+ for(size_t p = 0; p < nPerm; p++){
stream << "\\begin{tikzpicture}" << endl
- << "\\node[vertex] (n0) at(0, 0) {$0$};" << endl
- << "\\node[vertex] (n1) at(3, 0) {$1$};" << endl << endl
+ << "\\node[vertex] (n0) at(0, 0) {$0$};" << endl
+ << "\\node[vertex] (n1) at(3, 0) {$1$};" << endl << endl
- << "\\path[line]"
- << " (n" << (*(*(*edge)[p])[0])[0] << ")"
- << " -- "
- << " (n" << (*(*(*edge)[p])[0])[1] << ");"
- << endl
+ << "\\path[line]"
+ << " (n" << (*(*(*edge)[p])[0])[0] << ")"
+ << " -- "
+ << " (n" << (*(*(*edge)[p])[0])[1] << ");"
+ << endl
- << "\\end{tikzpicture} \\\\ \\\\" << endl << endl;
+ << "\\end{tikzpicture} \\\\ \\\\" << endl << endl;
}
stream << "\\end{longtable}" << endl
- << "\\end{document}" << endl;
+ << "\\end{document}" << endl;
return stream.str();
}
diff --git a/FunctionSpace/PermutationTree.cpp b/FunctionSpace/PermutationTree.cpp
new file mode 100644
index 0000000000..7f26080a18
--- /dev/null
+++ b/FunctionSpace/PermutationTree.cpp
@@ -0,0 +1,174 @@
+#include <sstream>
+
+#include "Exception.h"
+#include "PermutationTree.h"
+
+using namespace std;
+
+PermutationTree::PermutationTree(const vector<size_t>& refSequence){
+ // Sequence size //
+ sequenceSize = refSequence.size();
+
+ // List of leaf //
+ list<node_s*> listOfLeaf;
+
+ // Root //
+ root = new node_t;
+
+ root->myChoice = -1;
+ root->nxtChoice = refSequence;
+ root->father = NULL;
+
+ root->leafId = -1;
+ root->tag = -1;
+
+ // Tree //
+ populate(root, listOfLeaf);
+
+ // Leaf //
+ const size_t nLeaf = listOfLeaf.size();
+ leaf.assign(listOfLeaf.begin(), listOfLeaf.end());
+
+ for(size_t i = 0; i < nLeaf; i++)
+ leaf[i]->leafId = i;
+}
+
+PermutationTree::~PermutationTree(void){
+ destroy(root);
+}
+
+void PermutationTree::populate(node_t* node,
+ list<node_t*>& listOfLeaf){
+ // Number of son //
+ const size_t nSon = node->nxtChoice.size();
+
+ // Alloc space for son //
+ node->son.resize(nSon);
+
+ // Init each son //
+ for(size_t i = 0; i < nSon; i++){
+ // Alloc Son
+ node->son[i] = new node_t;
+
+ // Take father ith choice
+ node->son[i]->myChoice = node->nxtChoice[i];
+
+ // Next choices are all but ith
+ node->son[i]->nxtChoice.resize(nSon - 1);
+ for(size_t j = 0, k = 0; j < nSon; j++){
+ if(j != i){
+ node->son[i]->nxtChoice[k] = node->nxtChoice[j];
+ k++;
+ }
+
+ else{
+ }
+ }
+
+ // Dummy Stuff
+ node->son[i]->leafId = -1;
+ node->son[i]->tag = -1;
+
+ // Father
+ node->son[i]->father = node;
+
+ // Add Son
+ populate(node->son[i], listOfLeaf);
+ }
+
+ // If I am a leaf, add me to listOfLeaf
+ if(!nSon)
+ listOfLeaf.push_back(node);
+}
+
+void PermutationTree::destroy(node_t* node){
+ // Delete Son
+ const size_t nSon = node->son.size();
+
+ for(size_t i = 0; i < nSon; i++)
+ destroy(node->son[i]);
+
+ // Delete me
+ delete node;
+}
+
+
+size_t PermutationTree::getPermutationId(const vector<size_t>& sequence) const{
+ return getLeaf(root, sequence, 0)->leafId;
+}
+
+PermutationTree::node_t*
+PermutationTree::getLeaf(node_t* root,
+ const vector<size_t>& sequence,
+ size_t offset){
+ // Number of son
+ const size_t nSon = root->son.size();
+
+ // If a leaf return root
+ if(!nSon)
+ return root;
+
+ // Else, lookup
+ else{
+ // Find next son
+ size_t i = 0;
+ for(i = 0; i < nSon && root->nxtChoice[i] != sequence[offset]; i++)
+ ;
+
+ // If no math found, throw an Exception
+ if(i == nSon)
+ throw Exception("PermutationTree: cannot find given sequence");
+
+ // Get son leaf
+ return getLeaf(root->son[i], sequence, offset + 1);
+ }
+}
+
+void PermutationTree::fillWithPermutation(size_t permutationId,
+ vector<size_t>& vectorToFill) const{
+ // Go from ith leaf to root //
+ // And fill vectorToFill from end to start //
+
+ const size_t max = -1;
+
+ node_t* node = leaf[permutationId];
+ size_t i = sequenceSize - 1;
+
+ while(i != max && node){
+ vectorToFill[i] = node->myChoice;
+ i--;
+ node = node->father;
+ }
+}
+
+string PermutationTree::toString(void) const{
+ // Number of Permutation //
+ const size_t max_t = -1;
+ const size_t nPerm = leaf.size();
+
+ // Temporary //
+ stringstream stream;
+ vector<size_t> permutation(sequenceSize);
+
+ // Build String //
+ for(size_t i = 0; i < nPerm; i++){
+ fillWithPermutation(i, permutation);
+
+ stream << "Permutation #" << i << ":" << endl
+ << " -- " << permutation[0];
+
+ for(size_t j = 1; j < sequenceSize; j++)
+ stream << " " << permutation[j];
+
+ stream << endl << " -- Tag (";
+
+ if(leaf[i]-> tag != max_t)
+ stream << leaf[i]->tag;
+
+ stream << ")" << endl
+ << endl;
+ }
+
+ // Return //
+ return stream.str();
+}
diff --git a/FunctionSpace/PermutationTree.h b/FunctionSpace/PermutationTree.h
new file mode 100644
index 0000000000..7f13431d26
--- /dev/null
+++ b/FunctionSpace/PermutationTree.h
@@ -0,0 +1,83 @@
+#ifndef _PERMUTATIONTREE_H_
+#define _PERMUTATIONTREE_H_
+
+#include <cstdlib>
+#include <vector>
+#include <list>
+#include <string>
+
+class PermutationTree{
+ private:
+ // Tree node //
+ typedef struct node_s{
+ size_t myChoice;
+ std::vector<size_t> nxtChoice;
+
+ node_s* father;
+ std::vector<node_s*> son;
+
+ size_t leafId;
+ size_t tag;
+ } node_t;
+
+ private:
+ // Sequence size //
+ size_t sequenceSize;
+
+ // Root //
+ node_t* root;
+
+ // Leaf //
+ std::vector<node_t*> leaf;
+
+ public:
+ PermutationTree(const std::vector<size_t>& refSequence);
+ ~PermutationTree(void);
+
+ size_t getSequenceSize(void) const;
+
+ size_t getNPermutation(void) const;
+ size_t getPermutationId(const std::vector<size_t>& sequence) const;
+
+ void fillWithPermutation(size_t permutationId,
+ std::vector<size_t>& vectorToFill) const;
+
+ void addTagToPermutation(size_t permutationId, size_t tag);
+ size_t getTagFromPermutation(size_t permutationId);
+
+ std::string toString(void) const;
+
+ private:
+ static void populate(node_t* node,
+ std::list<node_t*>& listOfLeaf);
+
+ static void destroy(node_t* node);
+
+ static node_t* getLeaf(node_t* root,
+ const std::vector<size_t>& sequence,
+ size_t offset);
+};
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline size_t PermutationTree::getSequenceSize(void) const{
+ return sequenceSize;
+}
+
+inline size_t PermutationTree::getNPermutation(void) const{
+ return leaf.size();
+}
+
+inline void PermutationTree::
+addTagToPermutation(size_t permutationId, size_t tag){
+ leaf[permutationId]->tag = tag;
+}
+
+inline size_t PermutationTree::
+getTagFromPermutation(size_t permutationId){
+ return leaf[permutationId]->tag;
+}
+
+#endif
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 75e1962253..86c2a2c51f 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
// Reference Space //
refSpace = new QuadReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -104,7 +104,7 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
lagrange[(*(*edgeV[s])[e])[1]])).gradient());
}
- i++;
+ i++;
}
}
}
@@ -215,10 +215,10 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
// (u, v) = Gmsh Ref Quad
Polynomial mapX(Polynomial(0.5, 1, 0, 0) +
- Polynomial(0.5, 0, 0, 0));
+ Polynomial(0.5, 0, 0, 0));
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
- Polynomial(0.5, 0, 0, 0));
+ Polynomial(0.5, 0, 0, 0));
for(unsigned int s = 0; s < nRefSpace; s++){
for(unsigned int i = 0; i < nFunction; i++){
diff --git a/FunctionSpace/QuadNedelecBasis.cpp b/FunctionSpace/QuadNedelecBasis.cpp
index f8c6a2486b..7e4d245adb 100644
--- a/FunctionSpace/QuadNedelecBasis.cpp
+++ b/FunctionSpace/QuadNedelecBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
QuadNedelecBasis::QuadNedelecBasis(void){
// Reference Space //
refSpace = new QuadReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -86,10 +86,10 @@ QuadNedelecBasis::QuadNedelecBasis(void){
// (u, v) = Gmsh Ref Quad
Polynomial mapX(Polynomial(0.5, 1, 0, 0) +
- Polynomial(0.5, 0, 0, 0));
+ Polynomial(0.5, 0, 0, 0));
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
- Polynomial(0.5, 0, 0, 0));
+ Polynomial(0.5, 0, 0, 0));
for(unsigned int s = 0; s < nRefSpace; s++){
for(unsigned int i = 0; i < nFunction; i++){
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 200c9f1996..7034d5b9db 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
QuadNodeBasis::QuadNodeBasis(unsigned int order){
// Reference Space //
refSpace = new QuadReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -82,14 +82,14 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
for(unsigned int e = 0; e < 4; e++){
for(unsigned int l = 1; l < order; l++){
- basis[s][i] =
- new Polynomial(legendre[l].compose(lifting[(*(*edgeV[s])[e])[1]] -
+ basis[s][i] =
+ new Polynomial(legendre[l].compose(lifting[(*(*edgeV[s])[e])[1]] -
lifting[(*(*edgeV[s])[e])[0]])
*
- (lagrange[(*(*edgeV[s])[e])[0]] +
+ (lagrange[(*(*edgeV[s])[e])[0]] +
lagrange[(*(*edgeV[s])[e])[1]]));
- i++;
+ i++;
}
}
}
@@ -126,10 +126,10 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
// (u, v) = Gmsh Ref Quad
Polynomial mapX(Polynomial(0.5, 1, 0, 0) +
- Polynomial(0.5, 0, 0, 0));
+ Polynomial(0.5, 0, 0, 0));
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
- Polynomial(0.5, 0, 0, 0));
+ Polynomial(0.5, 0, 0, 0));
for(unsigned int s = 0; s < nRefSpace; s++){
for(unsigned int i = 0; i < nFunction; i++){
diff --git a/FunctionSpace/QuadReferenceSpace.cpp b/FunctionSpace/QuadReferenceSpace.cpp
index 5a55d5af8b..88e42b67b3 100644
--- a/FunctionSpace/QuadReferenceSpace.cpp
+++ b/FunctionSpace/QuadReferenceSpace.cpp
@@ -10,10 +10,10 @@ QuadReferenceSpace::QuadReferenceSpace(void){
// Edge Definition //
nEdge = 4;
- refEdge = new unsigned int*[nEdge];
+ refEdge = new size_t*[nEdge];
- for(unsigned int i = 0; i < nEdge; i++){
- refEdge[i] = new unsigned int[2];
+ for(size_t i = 0; i < nEdge; i++){
+ refEdge[i] = new size_t[2];
refEdge[i][0] = MQuadrangle::edges_quad(i, 0);
refEdge[i][1] = MQuadrangle::edges_quad(i, 1);
}
@@ -23,12 +23,12 @@ QuadReferenceSpace::QuadReferenceSpace(void){
nFace = 1;
// Number of node per face
- nNodeInFace = new unsigned int[nFace];
+ nNodeInFace = new size_t[nFace];
nNodeInFace[0] = 4;
// Reference Face
- refFace = new unsigned int*[nFace];
- refFace[0] = new unsigned int[4];
+ refFace = new size_t*[nFace];
+ refFace[0] = new size_t[4];
refFace[0][0] = 0;
refFace[0][1] = 1;
@@ -41,13 +41,13 @@ QuadReferenceSpace::QuadReferenceSpace(void){
QuadReferenceSpace::~QuadReferenceSpace(void){
// Delete Ref Edge //
- for(unsigned int i = 0; i < nEdge; i++)
+ for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
delete[] refEdge;
// Delete Ref Face //
- for(unsigned int i = 0; i < nFace; i++)
+ for(size_t i = 0; i < nFace; i++)
delete[] refFace[i];
delete[] refFace;
@@ -55,7 +55,9 @@ QuadReferenceSpace::~QuadReferenceSpace(void){
}
string QuadReferenceSpace::toLatex(void) const{
- stringstream stream;
+ const size_t nPerm = pTree->getNPermutation();
+ stringstream stream;
+ vector<size_t> perm(nVertex);
stream << "\\documentclass{article}" << endl << endl
@@ -69,16 +71,18 @@ string QuadReferenceSpace::toLatex(void) const{
<< "\\begin{longtable}{ccc}" << endl << endl;
- for(unsigned int p = 0; p < nPerm; p++){
+ for(size_t p = 0; p < nPerm; p++){
+ pTree->fillWithPermutation(p, perm);
+
stream << "\\begin{tikzpicture}" << endl
- << "\\node[vertex] (n0) at(0, 0) {$" << perm[p][0] << "$};" << endl
- << "\\node[vertex] (n1) at(3, 0) {$" << perm[p][1] << "$};" << endl
- << "\\node[vertex] (n2) at(3, 3) {$" << perm[p][2] << "$};" << endl
- << "\\node[vertex] (n3) at(0, 3) {$" << perm[p][3] << "$};" << endl
+ << "\\node[vertex] (n0) at(0, 0) {$" << perm[0] << "$};" << endl
+ << "\\node[vertex] (n1) at(3, 0) {$" << perm[1] << "$};" << endl
+ << "\\node[vertex] (n2) at(3, 3) {$" << perm[2] << "$};" << endl
+ << "\\node[vertex] (n3) at(0, 3) {$" << perm[3] << "$};" << endl
<< endl;
- for(unsigned int i = 0; i < 4; i++)
+ for(size_t i = 0; i < 4; i++)
stream << "\\path[line]"
<< " (n" << (*(*(*edge)[p])[i])[0] << ")"
<< " -- "
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index 241d87bd2d..110259b920 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -8,19 +8,6 @@ using namespace std;
ReferenceSpace::ReferenceSpace(void){
// Init to NULL //
- nVertex = 0;
- nextLeafId = 0;
-
- nPerm = 0;
- perm = NULL;
- lPerm = NULL;
-
- pTreeRoot.depth = 0;
- pTreeRoot.last = NULL;
- pTreeRoot.number = 0;
- pTreeRoot.possible = NULL;
- pTreeRoot.next = NULL;
-
nEdge = 0;
refEdge = NULL;
//permutedRefEdge = NULL;
@@ -33,17 +20,18 @@ ReferenceSpace::ReferenceSpace(void){
face = NULL;
// Defining Ref Edge and Face in //
- // Dervived Class //
+ // Derived Class //
}
ReferenceSpace::~ReferenceSpace(void){
+ const size_t nPerm = pTree->getNPermutation();
+
// Destroy Tree //
- destroy(&pTreeRoot);
- delete[] perm;
+ delete pTree;
// Delete Permutated Edge //
- for(unsigned int p = 0; p < nPerm; p++){
- for(unsigned int i = 0; i < nEdge; i++){
+ for(size_t p = 0; p < nPerm; p++){
+ for(size_t i = 0; i < nEdge; i++){
//delete[] permutedRefEdge[p][i];
delete (*(*edge)[p])[i];
}
@@ -56,8 +44,8 @@ ReferenceSpace::~ReferenceSpace(void){
delete edge;
// Delete Permutated Face //
- for(unsigned int p = 0; p < nPerm; p++){
- for(unsigned int i = 0; i < nFace; i++){
+ for(size_t p = 0; p < nPerm; p++){
+ for(size_t i = 0; i < nFace; i++){
//delete[] permutedRefFace[p][i];
delete (*(*face)[p])[i];
}
@@ -71,113 +59,233 @@ ReferenceSpace::~ReferenceSpace(void){
}
void ReferenceSpace::init(void){
- // Init Root //
- nPerm = 0;
- nextLeafId = 0;
-
- pTreeRoot.depth = 0;
- pTreeRoot.last = NULL;
- pTreeRoot.number = nVertex;
- pTreeRoot.possible = new unsigned int[pTreeRoot.number];
- pTreeRoot.next = NULL;
-
- for(unsigned int i = 0; i < pTreeRoot.number; i++)
- pTreeRoot.possible[i] = i;
-
- // Populate Tree //
- lPerm = new list<unsigned int*>;
- populate(&pTreeRoot);
-
- // Get Permutation //
- perm = new unsigned int*[nPerm];
- for(unsigned int i = 0; i < nPerm; i++){
- // Take Permutation for queue
- // (AND IN ORDER)
- perm[i] = lPerm->front();
- lPerm->pop_front();
- }
+ // Tree //
+ vector<size_t> vertexSeq(nVertex);
+
+ for(size_t i = 0; i < nVertex; i++)
+ vertexSeq[i] = i;
- delete lPerm;
+ pTree = new PermutationTree(vertexSeq);
- // Get Edges & Faces //
+ // Edges & Faces //
getEdge();
getFace();
+
+ // Cyclic Permutation //
+ findCyclicPermutation();
}
-void ReferenceSpace::populate(node* pTreeRoot){
- // Get Some Data on this Root //
- const unsigned int number = pTreeRoot->number;
- const unsigned int nextNumber = number - 1;
- const unsigned int depth = pTreeRoot->depth;
- const unsigned int nextDepth = depth + 1;
-
- // Temp Data //
- unsigned int offset;
-
- // If Leaf : a new permutation is found //
- if(!number){
- // Init Permutation
- pTreeRoot->next = NULL;
- pTreeRoot->leafId = nextLeafId;
-
- // Value for Next Permutation
- nextLeafId++;
- nPerm++;
-
- // Put this Permutation in queue
- // (AND IN ORDER)
- lPerm->push_back(pTreeRoot->last);
- }
+void ReferenceSpace::findCyclicPermutation(void){
+ // Alloc Some Data //
+ const size_t nPerm = pTree->getNPermutation();
+
+ list<size_t> listOfTrueReferenceSpace;
+ list<vector<size_t> > listOfRefNodeIndexPermutation;
+ list<vector<size_t> > listOfReverseNodeIndexPermutation;
+
+ vector<size_t> pTest(nVertex);
+ vector<size_t> pRef(nVertex);
- // Else: continue to build the tree //
- else{
- // We got 'number' child nodes
- pTreeRoot->next = new node[number];
+ list<size_t>::iterator it;
+ list<size_t>::iterator end;
+ triple match;
- // Init each child node
- for(unsigned int i = 0; i < number; i++){
- // Reset offset
- offset = 0;
+ vector<size_t> unPermutedIndex(nVertex);
- // Depth and Last Choices of child nodes
- pTreeRoot->next[i].depth = nextDepth;
- pTreeRoot->next[i].last = new unsigned int[nextDepth];
- pTreeRoot->next[i].last[depth] = pTreeRoot->possible[i];
+ for(size_t i = 0; i < nVertex; i++)
+ unPermutedIndex[i] = i;
- for(unsigned int j = 0; j < depth; j++)
- pTreeRoot->next[i].last[j] = pTreeRoot->last[j];
+ // Find Cyclic Permutation
+ for(size_t i = 0; i < nPerm; i++){
+ // No match
+ match.first = false;
- // Possibilities of child node
- pTreeRoot->next[i].number = nextNumber;
- pTreeRoot->next[i].possible = new unsigned int[nextNumber];
+ // Get Permutation 'i'
+ pTree->fillWithPermutation(i, pTest);
+
+ // Test it with already found Reference Space
+ it = listOfTrueReferenceSpace.begin();
+ end = listOfTrueReferenceSpace.end();
+
+ while(it != end && !match.first){
+ // Take Reference Space 'it'
+ pTree->fillWithPermutation(*it, pRef);
+
+ // Look if it matches
+ match = isCyclicPermutation(pTest, pRef);
+
+ // If not, go to next Reference Space
+ if(!match.first)
+ it++;
+ }
+
+ // If no Reference Space is found
+ // --> this Permutation is a new Reference Space
+ if(!match.first){
+ listOfTrueReferenceSpace.push_back(i);
+ listOfRefNodeIndexPermutation.push_back(unPermutedIndex);
+ listOfReverseNodeIndexPermutation.push_back(unPermutedIndex);
+
+ pTree->addTagToPermutation(i, i);
+ }
+
+ // If a ReferenceSpace is found, and the index permutations
+ else{
+ listOfRefNodeIndexPermutation.push_back(match.second);
+ listOfReverseNodeIndexPermutation.push_back(match.third);
+ pTree->addTagToPermutation(i, *it);
+ }
+ }
+}
- for(unsigned int j = 0; j < nextNumber; j++){
- if(pTreeRoot->possible[j] == pTreeRoot->possible[i])
- offset = 1;
+static bool isFacePermutation(vector<size_t>& refNode,
+ vector<size_t>& testNode){
+ const size_t size = refNode.size();
+ bool match = false;
- pTreeRoot->next[i].possible[j] = pTreeRoot->possible[j + offset];
+ if(size != testNode.size())
+ return false;
+
+ for(size_t i = 0; i < size && !match; i++){
+ bool submatch = true;
+
+ for(size_t j = 0; j < size && submatch; j++)
+ if(refNode[j] != testNode[j])
+ submatch = false;
+
+ if(!submatch){
+ size_t tmp0 = testNode[0];
+ size_t tmp1 = testNode[1];
+
+ for(size_t k = 1; k < size + 1; k++){
+ testNode[k % size] = tmp0;
+ tmp0 = tmp1;
+ tmp1 = testNode[(k + 1) % size];
}
+ }
+
+ else
+ match = true;
+ }
+
+ return match;
+}
- // Populate each child node (until a leaf is found)
- populate(&pTreeRoot->next[i]);
+static bool haveSameNode(vector<size_t>& face0,
+ vector<size_t>& face1){
+ const size_t size = face0.size();
+ bool matchIsPossible = true;
+
+ for(size_t i = 0; i < size && matchIsPossible; i++){
+ bool submatch = false;
+
+ for(size_t j = 0; j < size && !submatch; j++){
+ if(face0[i] == face1[j])
+ submatch = true;
}
+
+ matchIsPossible = submatch;
}
+
+ return matchIsPossible;
}
-void ReferenceSpace::destroy(node* node){
- const unsigned int number = node->number;
+size_t ReferenceSpace::findCorrespondingFace(vector<size_t>& face,
+ vector<size_t>& node){
+ // Init Stuff //
+ const size_t faceSize = face.size();
+ bool match = false;
+ size_t f = 0;
+
+ vector<size_t> testFace(faceSize);
+
+ // Test All Face until match
+ while(!match && f < nFace){
+
+ if(nNodeInFace[f] == faceSize){
+ // Get face f nodes
+ for(size_t i = 0; i < faceSize; i++)
+ testFace[i] = node[refFace[f][i]];
+
+ // Look if match
+ match = haveSameNode(testFace, face);
+ }
- for(unsigned int i = 0; i < number; i++){
- destroy(&node->next[i]);
- node->number--;
+ if(!match)
+ f++;
}
- delete[] node->possible;
- delete[] node->last;
- delete[] node->next;
+ return f;
+}
+
+static vector<size_t> getRefIndexPermutation(vector<size_t>& ref,
+ vector<size_t>& test){
+ const size_t size = ref.size();
+ vector<size_t> idxVec(ref.size());
+ size_t idx;
+
+ for(size_t i = 0; i < size; i++){
+ idx = 0;
+
+ while(test[i] != ref[idx])
+ idx++;
+
+ idxVec[i] = idx;
+ }
+
+ return idxVec;
+}
+
+static vector<size_t> getReverseIndexPermutation(vector<size_t>& ref,
+ vector<size_t>& test){
+ const size_t size = ref.size();
+ vector<size_t> idxVec(ref.size());
+ size_t idx;
+
+ for(size_t i = 0; i < size; i++){
+ idx = 0;
+
+ while(test[idx] != ref[i])
+ idx++;
+
+ idxVec[i] = idx;
+ }
+
+ return idxVec;
+}
+
+ReferenceSpace::triple ReferenceSpace::
+isCyclicPermutation(vector<size_t>& pTest,
+ vector<size_t>& pRef){
+
+ // Node IDs of Reference Space first Face
+ vector<size_t> refNode(nNodeInFace[0]);
+
+ for(size_t i = 0; i < nNodeInFace[0]; i++)
+ refNode[i] = pRef[refFace[0][i]];
+
+ // Corresponding Face in Test Permutation
+ size_t testFaceId = findCorrespondingFace(refNode, pTest);
+
+ // Node IDs of Test Permutation correspnding Face
+ vector<size_t> testNode(nNodeInFace[testFaceId]);
+
+ for(size_t i = 0; i < nNodeInFace[testFaceId]; i++)
+ testNode[i] = pTest[refFace[testFaceId][i]];
+
+ // Return Triple //
+ triple tri = {
+ isFacePermutation(refNode, testNode),
+ getRefIndexPermutation(pRef, pTest),
+ getReverseIndexPermutation(pRef, pTest)
+ };
+
+ return tri;
}
void ReferenceSpace::getEdge(void){
+ const size_t nPerm = pTree->getNPermutation();
+
// Alloc
vector<const vector<unsigned int>*>* tmp;
edge = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
@@ -199,10 +307,10 @@ void ReferenceSpace::getEdge(void){
delete[] nNodeInEdge;
*/
// Populate Edge
- for(unsigned int p = 0; p < nPerm; p++){
+ for(size_t p = 0; p < nPerm; p++){
tmp = new vector<const vector<unsigned int>*>(nEdge);
- for(unsigned int e = 0; e < nEdge; e++)
+ for(size_t e = 0; e < nEdge; e++)
(*tmp)[e] = getOrderedEdge(p, e);
(*edge)[p] = tmp;
@@ -210,6 +318,8 @@ void ReferenceSpace::getEdge(void){
}
void ReferenceSpace::getFace(void){
+ const size_t nPerm = pTree->getNPermutation();
+
// Alloc
vector<const vector<unsigned int>*>* tmp;
face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
@@ -223,10 +333,10 @@ void ReferenceSpace::getFace(void){
nFace);
*/
// Populate Face
- for(unsigned int p = 0; p < nPerm; p++){
+ for(size_t p = 0; p < nPerm; p++){
tmp = new vector<const vector<unsigned int>*>(nFace);
- for(unsigned int f = 0; f < nFace; f++){
+ for(size_t f = 0; f < nFace; f++){
// Dending on the number of node per face
// The ordering strategy is different
@@ -281,9 +391,13 @@ getOrderedEdge(unsigned int permutation,
// Alloc
vector<unsigned int>* ordered = new vector<unsigned int>(2);
+ // Permutation
+ vector<size_t> perm(nVertex);
+ pTree->fillWithPermutation(permutation, perm);
+
// Order refEdge
orderRefEntityForGivenPermutation(refEdge[edge],
- perm[permutation],
+ perm,
*ordered);
// Return ordered
@@ -301,9 +415,13 @@ getOrderedTriFace(unsigned int permutation,
// Alloc
vector<unsigned int>* ordered = new vector<unsigned int>(3);
+ // Permutation
+ vector<size_t> perm(nVertex);
+ pTree->fillWithPermutation(permutation, perm);
+
// Order refFace
orderRefEntityForGivenPermutation(refFace[face],
- perm[permutation],
+ perm,
*ordered);
// Return ordered
@@ -327,9 +445,13 @@ getOrderedQuadFace(unsigned int permutation,
// Alloc
vector<unsigned int>* ordered = new vector<unsigned int>(4);
+ // Permutation
+ vector<size_t> perm(nVertex);
+ pTree->fillWithPermutation(permutation, perm);
+
// Order refFace
orderRefEntityForGivenPermutation(refFace[face],
- perm[permutation],
+ perm,
*ordered);
// Get ordered[2] opposite to ordered[0]
@@ -353,8 +475,8 @@ getOrderedQuadFace(unsigned int permutation,
}
void ReferenceSpace::
-orderRefEntityForGivenPermutation(unsigned int* refEntity,
- unsigned int* permutation,
+orderRefEntityForGivenPermutation(size_t* refEntity,
+ vector<size_t>& permutation,
vector<unsigned int>& orderedEntity){
// Get Size
const unsigned int size = orderedEntity.size();
@@ -375,7 +497,7 @@ orderRefEntityForGivenPermutation(unsigned int* refEntity,
orderedEntity[i] = refEntity[sorted[i].first];
}
-unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
+size_t ReferenceSpace::getReferenceSpace(const MElement& elem) const{
// Const_Cast //
MElement& element = const_cast<MElement&>(elem);
@@ -393,14 +515,14 @@ unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
std::sort(vertexGlobalId.begin(), vertexGlobalId.end(), sortPredicate);
// Reduce Vertex Global ID //
- vector<unsigned int> vertexReducedId(nVertex);
+ vector<size_t> vertexReducedId(nVertex);
for(int i = 0; i < nVertex; i++)
vertexReducedId[vertexGlobalId[i].first] = i;
// Tree Lookup //
try{
- return treeLookup(&pTreeRoot, vertexReducedId);
+ return pTree->getPermutationId(vertexReducedId);
}
catch(...){
@@ -409,60 +531,20 @@ unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
}
}
-unsigned int ReferenceSpace::treeLookup(const node* node,
- vector<unsigned int>& vertexReducedId){
- // Temp Data //
- unsigned int choice;
- unsigned int i;
-
- // If Node is *not* a Leaf: Lookup //
- if(node->number){
- // Get This Choice
- choice = vertexReducedId[node->depth];
-
- // Look for next node corresponding to this Choice
- i = 0;
- while(node->possible[i] != choice)
- i++;
-
- // Look if a this Choice has been found
- if(i == node->number)
- throw Exception();
-
- // Go to next Node
- return treeLookup(&node->next[i], vertexReducedId);
- }
-
- // Else: Return Leaf ID //
- else
- return node->leafId;
-}
-
string ReferenceSpace::toString(void) const{
- stringstream stream;
-
- // Tree //
- stream << "Tree:" << endl;
- stream << toString(&pTreeRoot) << endl;
+ const size_t nPerm = pTree->getNPermutation();
+ stringstream stream;
// ReferenceSpaces //
- stream << "Reference Spaces:" << endl;
-
- for(unsigned int i = 0; i < nPerm; i++){
- stream << " * ";
-
- for(unsigned int j = 0; j < nVertex; j++)
- stream << perm[i][j] << " ";
-
- stream << " (# " << i + 1 << ")" << endl;
- }
+ stream << "Reference Spaces:" << endl
+ << pTree->toString() << endl;
stream << "Edges Permutations:" << endl;
- for(unsigned int i = 0; i < nPerm; i++){
+ for(size_t i = 0; i < nPerm; i++){
stream << " * RefSpace #" << i + 1 << ":" << endl;
- for(unsigned int j = 0; j < nEdge; j++)
+ for(size_t j = 0; j < nEdge; j++)
stream << " -- ["
<< edge->at(i)->at(j)->at(0) << ", "
<< edge->at(i)->at(j)->at(1) << "]" << endl;
@@ -470,10 +552,10 @@ string ReferenceSpace::toString(void) const{
stream << "Faces Permutations:" << endl;
- for(unsigned int i = 0; i < nPerm; i++){
+ for(size_t i = 0; i < nPerm; i++){
stream << " * RefSpace #" << i + 1 << ":" << endl;
- for(unsigned int j = 0; j < nFace; j++)
+ for(size_t j = 0; j < nFace; j++)
stream << " -- ["
<< face->at(i)->at(j)->at(0) << ", "
<< face->at(i)->at(j)->at(1) << ", "
@@ -483,23 +565,6 @@ string ReferenceSpace::toString(void) const{
return stream.str();
}
-string ReferenceSpace::toString(const node* node) const{
- const unsigned int number = node->number;
- stringstream stream;
-
- if(node->last)
- stream << "(" << node->last[node->depth - 1] << " ";
- else
- stream << "(root ";
-
- for(unsigned int i = 0; i < number; i++)
- stream << toString(&node->next[i]);
-
- stream << ")";
-
- return stream.str();
-}
-
string ReferenceSpace::toLatex(void) const{
stringstream stream;
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index e0dafe6441..aa04bf78b9 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -5,7 +5,9 @@
#include <list>
#include <stack>
#include <string>
+
#include "MElement.h"
+#include "PermutationTree.h"
/**
@interface ReferenceSpace
@@ -19,31 +21,15 @@
class ReferenceSpace{
private:
- // Permuation Tree Structure //
- struct node_s{
- unsigned int depth; // Depth
- unsigned int* last; // Last Choises
- unsigned int number; // Number of Next Choises
- unsigned int* possible; // Possible Next Choises
- node_s* next; // Next Choise
-
- unsigned int leafId; // If leaf: this leaf number
- // (from 0 to nLeaf - 1)
- // Else: no meaning
- };
-
- typedef node_s node;
+ typedef struct{
+ bool first;
+ std::vector<size_t> second;
+ std::vector<size_t> third;
+ } triple;
protected:
- // Permutation (Tree + Leaf) of Reduced IDs //
- unsigned int nVertex;
- unsigned int nextLeafId;
-
- unsigned int nPerm;
- unsigned int** perm;
- std::list<unsigned int*>* lPerm;
-
- node pTreeRoot;
+ // Permutation Tree //
+ PermutationTree* pTree;
//////////////////////////////////////////////////////////////////////////
// !!! CHANGE VARIABLES NAME !!! //
@@ -72,26 +58,32 @@ class ReferenceSpace{
// - - - - //
//////////////////////////////////////////////////////////////////////////
+ // Number of Vertices //
+ size_t nVertex;
+
// Edge Permutation //
- unsigned int nEdge;
- unsigned int** refEdge;
+ size_t nEdge;
+ size_t** refEdge;
//unsigned int*** permutedRefEdge;
std::vector<const std::vector<const std::vector<unsigned int>*>*>* edge;
// Face Permutation //
- unsigned int nFace;
- unsigned int* nNodeInFace;
- unsigned int** refFace;
+ size_t nFace;
+ size_t* nNodeInFace;
+ size_t** refFace;
//unsigned int*** permutedRefFace;
std::vector<const std::vector<const std::vector<unsigned int>*>*>* face;
public:
virtual ~ReferenceSpace(void);
- unsigned int getNPermutation(void) const;
-
- unsigned int getPermutation(const MElement& element) const;
+ size_t getNReferenceSpace(void) const;
+ size_t getReferenceSpace(const MElement& element) const;
+ /*
+ std::vector<std::vector<size_t> > getEdgeIndex(size_t referenceSpaceId) const;
+ std::vector<std::vector<size_t> > getFaceIndex(size_t referenceSpaceId) const;
+ */
const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
getAllEdge(void) const;
@@ -106,11 +98,17 @@ class ReferenceSpace{
ReferenceSpace(void);
void init(void);
- void populate(node* pTreeRoot);
- void destroy(node* node);
void getEdge(void);
void getFace(void);
+
+ void findCyclicPermutation(void);
+ triple isCyclicPermutation(std::vector<size_t>& pTest,
+ std::vector<size_t>& pRef);
+
+ size_t findCorrespondingFace(std::vector<size_t>& face,
+ std::vector<size_t>& node);
+
/*
void getPermutedRefEntity(unsigned int**** permutedRefEntity,
unsigned int** refEntity,
@@ -127,17 +125,12 @@ class ReferenceSpace{
unsigned int face);
static void
- orderRefEntityForGivenPermutation(unsigned int* refEntity,
- unsigned int* permutation,
+ orderRefEntityForGivenPermutation(size_t* refEntity,
+ std::vector<size_t>& permutation,
std::vector<unsigned int>& orderedEntity);
static bool sortPredicate(const std::pair<unsigned int, unsigned int>& a,
const std::pair<unsigned int, unsigned int>& b);
-
- static unsigned int treeLookup(const node* node,
- std::vector<unsigned int>& vertexReducedId);
-
- std::string toString(const node* node) const;
};
@@ -152,12 +145,12 @@ class ReferenceSpace{
Deletes this ReferenceSpace
**
- @fn ReferenceSpace::getNPermutation
+ @fn ReferenceSpace::getNReferenceSpace
@returns Returns the number of permutation of this
ReferenceSpace
**
- @fn ReferenceSpace::getPermutation
+ @fn ReferenceSpace::getReferenceSpace
@param element A MElement
@returns Returns the a natural number defining
the permutation of the given element
@@ -172,7 +165,7 @@ class ReferenceSpace{
@note
@li The fisrt vector represents a particular permutation
- (see ReferenceSpace::getPermutation())
+ (see ReferenceSpace::getReferenceSpace())
@li The second vector represents a particular edge
(for a given permutation)
@li The last vector represents the Vertex @c IDs of
@@ -184,7 +177,7 @@ class ReferenceSpace{
@note
@li The fisrt vector represents a particular permutation
- (see ReferenceSpace::getPermutation())
+ (see ReferenceSpace::getReferenceSpace())
@li The second vector represents a particular face
(for a given permutation)
@li The last vector represents the Vertex @c IDs of
@@ -204,10 +197,8 @@ class ReferenceSpace{
// Inline Functions //
//////////////////////
-inline
-unsigned int
-ReferenceSpace::getNPermutation(void) const{
- return nPerm;
+inline size_t ReferenceSpace::getNReferenceSpace(void) const{
+ return pTree->getNPermutation();
}
inline
diff --git a/FunctionSpace/ReferenceSpaceLagrange.cpp b/FunctionSpace/ReferenceSpaceLagrange.cpp
index 4a5ad4ec4f..0a618ea115 100644
--- a/FunctionSpace/ReferenceSpaceLagrange.cpp
+++ b/FunctionSpace/ReferenceSpaceLagrange.cpp
@@ -9,9 +9,11 @@ ReferenceSpaceLagrange::ReferenceSpaceLagrange(void){
}
ReferenceSpaceLagrange::~ReferenceSpaceLagrange(void){
+ const size_t nPerm = pTree->getNPermutation();
+
// If 'node' allocated //
if(node){
- for(unsigned int p = 0; p < nPerm; p++)
+ for(size_t p = 0; p < nPerm; p++)
delete (*node)[p];
delete node;
@@ -19,6 +21,8 @@ ReferenceSpaceLagrange::~ReferenceSpaceLagrange(void){
}
void ReferenceSpaceLagrange::getLagrangeNode(void){
+ const size_t nPerm = pTree->getNPermutation();
+
// Alloc //
vector<unsigned int>* tmp;
node = new vector<const vector<unsigned int>*>(nPerm);
@@ -60,7 +64,7 @@ void ReferenceSpaceLagrange::edgeSeq(vector<unsigned int>& vec,
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
- unsigned int* refEdge,
+ size_t* refEdge,
const vector<unsigned int>& edge){
// Is reverted ? //
@@ -93,7 +97,7 @@ void ReferenceSpaceLagrange::faceSeq(vector<unsigned int>& vec,
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
- unsigned int* refFace,
+ size_t* refFace,
const vector<unsigned int>& face,
unsigned int nNodePerEdge){
// Index //
@@ -108,6 +112,7 @@ void ReferenceSpaceLagrange::faceSeq(vector<unsigned int>& vec,
}
string ReferenceSpaceLagrange::toString(void) const{
+ const size_t nPerm = pTree->getNPermutation();
const unsigned int nNodeMinus = nNode - 1;
stringstream stream;
diff --git a/FunctionSpace/ReferenceSpaceLagrange.h b/FunctionSpace/ReferenceSpaceLagrange.h
index 0fd71a7912..92c088018b 100644
--- a/FunctionSpace/ReferenceSpaceLagrange.h
+++ b/FunctionSpace/ReferenceSpaceLagrange.h
@@ -50,14 +50,14 @@ class ReferenceSpaceLagrange: public ReferenceSpace{
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
- unsigned int* refEdge,
+ size_t* refEdge,
const std::vector<unsigned int>& edge);
static void faceSeq(std::vector<unsigned int>& vec,
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
- unsigned int* refFace,
+ size_t* refFace,
const std::vector<unsigned int>& face,
unsigned int nNodePerEdge);
};
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 73795d0341..405ce362d3 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Reference Space //
refSpace = new TetReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -45,9 +45,9 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
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(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0) -
+ Polynomial(1, 0, 0, 1)),
Polynomial(Polynomial(1, 1, 0, 0)),
@@ -69,37 +69,37 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
for(int e = 0; e < 6; e++){
for(int l = 0; l < orderPlus; l++){
- // Nedelec
- if(l == 0){
- 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][i] = new vector<Polynomial>(tmp2);
- }
-
- // High Order
- else{
- basis[s][i] =
- new vector<Polynomial>
- ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
- ,
- lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]])).gradient());
- }
- i++;
+ // Nedelec
+ if(l == 0){
+ 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][i] = new vector<Polynomial>(tmp2);
+ }
+
+ // High Order
+ else{
+ basis[s][i] =
+ new vector<Polynomial>
+ ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
+ lagrange[(*(*edgeV[s])[e])[1]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[0]] +
+ lagrange[(*(*edgeV[s])[e])[1]])).gradient());
+ }
+ i++;
}
}
}
@@ -112,86 +112,86 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
for(int f = 0; f < 4; f++){
for(unsigned int l1 = 1; l1 < order; l1++){
- for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
- // Preliminary Type 1
- Polynomial sum =
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]] +
- lagrange[(*(*faceV[s])[f])[2]];
-
- Polynomial u =
- intLegendre[l1].compose(lagrange[(*(*faceV[s])[f])[0]] -
- lagrange[(*(*faceV[s])[f])[1]]
- ,
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]]);
- Polynomial v =
- lagrange[(*(*faceV[s])[f])[2]] *
- sclLegendre[l2].compose(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum);
-
- // Preliminary Type 2
- vector<Polynomial> gradU = u.gradient();
- vector<Polynomial> gradV = v.gradient();
-
- vector<Polynomial> vGradU(gradU);
- vGradU[0].mul(v);
- vGradU[1].mul(v);
- vGradU[2].mul(v);
-
- vector<Polynomial> uGradV(gradV);
- uGradV[0].mul(u);
- uGradV[1].mul(u);
- uGradV[2].mul(u);
-
- vector<Polynomial> subGradUV(vGradU);
- subGradUV[0].sub(uGradV[0]);
- subGradUV[1].sub(uGradV[1]);
- subGradUV[2].sub(uGradV[2]);
-
- // Preliminary Type 3
- vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[f])[0]].gradient();
- vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[f])[1]].gradient();
-
- vector<Polynomial> l2GradL1(gradL1);
- l2GradL1[0].mul(lagrange[(*(*faceV[s])[f])[1]]);
- l2GradL1[1].mul(lagrange[(*(*faceV[s])[f])[1]]);
- l2GradL1[2].mul(lagrange[(*(*faceV[s])[f])[1]]);
-
- vector<Polynomial> l1GradL2(gradL2);
- l1GradL2[0].mul(lagrange[(*(*faceV[s])[f])[0]]);
- l1GradL2[1].mul(lagrange[(*(*faceV[s])[f])[0]]);
- l1GradL2[2].mul(lagrange[(*(*faceV[s])[f])[0]]);
-
- vector<Polynomial> subGradL1L2V(l2GradL1);
- subGradL1L2V[0].sub(l1GradL2[0]);
- subGradL1L2V[1].sub(l1GradL2[1]);
- subGradL1L2V[2].sub(l1GradL2[2]);
-
- subGradL1L2V[0].mul(v);
- subGradL1L2V[1].mul(v);
- subGradL1L2V[2].mul(v);
-
-
- // Type 1
- basis[s][i] =
- new vector<Polynomial>((u * v).gradient());
-
- i++;
-
- // Type 2
- basis[s][i] =
- new vector<Polynomial>(subGradUV);
-
- i++;
-
- // Type 3
- if(l1 == 1){
- basis[s][i] =
- new vector<Polynomial>(subGradL1L2V);
-
- i++;
- }
- }
+ for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
+ // Preliminary Type 1
+ Polynomial sum =
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]] +
+ lagrange[(*(*faceV[s])[f])[2]];
+
+ Polynomial u =
+ intLegendre[l1].compose(lagrange[(*(*faceV[s])[f])[0]] -
+ lagrange[(*(*faceV[s])[f])[1]]
+ ,
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]]);
+ Polynomial v =
+ lagrange[(*(*faceV[s])[f])[2]] *
+ sclLegendre[l2].compose(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum);
+
+ // Preliminary Type 2
+ vector<Polynomial> gradU = u.gradient();
+ vector<Polynomial> gradV = v.gradient();
+
+ vector<Polynomial> vGradU(gradU);
+ vGradU[0].mul(v);
+ vGradU[1].mul(v);
+ vGradU[2].mul(v);
+
+ vector<Polynomial> uGradV(gradV);
+ uGradV[0].mul(u);
+ uGradV[1].mul(u);
+ uGradV[2].mul(u);
+
+ vector<Polynomial> subGradUV(vGradU);
+ subGradUV[0].sub(uGradV[0]);
+ subGradUV[1].sub(uGradV[1]);
+ subGradUV[2].sub(uGradV[2]);
+
+ // Preliminary Type 3
+ vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[f])[0]].gradient();
+ vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[f])[1]].gradient();
+
+ vector<Polynomial> l2GradL1(gradL1);
+ l2GradL1[0].mul(lagrange[(*(*faceV[s])[f])[1]]);
+ l2GradL1[1].mul(lagrange[(*(*faceV[s])[f])[1]]);
+ l2GradL1[2].mul(lagrange[(*(*faceV[s])[f])[1]]);
+
+ vector<Polynomial> l1GradL2(gradL2);
+ l1GradL2[0].mul(lagrange[(*(*faceV[s])[f])[0]]);
+ l1GradL2[1].mul(lagrange[(*(*faceV[s])[f])[0]]);
+ l1GradL2[2].mul(lagrange[(*(*faceV[s])[f])[0]]);
+
+ vector<Polynomial> subGradL1L2V(l2GradL1);
+ subGradL1L2V[0].sub(l1GradL2[0]);
+ subGradL1L2V[1].sub(l1GradL2[1]);
+ subGradL1L2V[2].sub(l1GradL2[2]);
+
+ subGradL1L2V[0].mul(v);
+ subGradL1L2V[1].mul(v);
+ subGradL1L2V[2].mul(v);
+
+
+ // Type 1
+ basis[s][i] =
+ new vector<Polynomial>((u * v).gradient());
+
+ i++;
+
+ // Type 2
+ basis[s][i] =
+ new vector<Polynomial>(subGradUV);
+
+ i++;
+
+ // Type 3
+ if(l1 == 1){
+ basis[s][i] =
+ new vector<Polynomial>(subGradL1L2V);
+
+ i++;
+ }
+ }
}
}
}
@@ -204,114 +204,114 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
- for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusTwo; l3++){
- // Preliminary Type 1
- Polynomial u = intLegendre[l1].compose
- (lagrange[0] - lagrange[1],
- lagrange[0] + lagrange[1]);
-
- Polynomial v = lagrange[2] * sclLegendre[l2].compose
- (lagrange[2] * 2 - (one - lagrange[3]), one - lagrange[3]);
-
- Polynomial w = lagrange[3] * legendre[l3].compose
- (lagrange[3] * 2 - one);
-
- // Preliminary Type 2
- vector<Polynomial> gradU = u.gradient();
- vector<Polynomial> gradV = v.gradient();
- vector<Polynomial> gradW = w.gradient();
-
- vector<Polynomial> vwGradU(gradU);
- vwGradU[0].mul(v);
- vwGradU[1].mul(v);
- vwGradU[2].mul(v);
-
- vwGradU[0].mul(w);
- vwGradU[1].mul(w);
- vwGradU[2].mul(w);
-
- vector<Polynomial> uwGradV(gradV);
- uwGradV[0].mul(u);
- uwGradV[1].mul(u);
- uwGradV[2].mul(u);
-
- uwGradV[0].mul(w);
- uwGradV[1].mul(w);
- uwGradV[2].mul(w);
-
- vector<Polynomial> uvGradW(gradW);
- uvGradW[0].mul(u);
- uvGradW[1].mul(u);
- uvGradW[2].mul(u);
-
- uvGradW[0].mul(v);
- uvGradW[1].mul(v);
- uvGradW[2].mul(v);
-
- vector<Polynomial> term1(vwGradU);
- term1[0].sub(uwGradV[0]);
- term1[1].sub(uwGradV[1]);
- term1[2].sub(uwGradV[2]);
-
- term1[0].add(uvGradW[0]);
- term1[1].add(uvGradW[1]);
- term1[2].add(uvGradW[2]);
-
- vector<Polynomial> term2(vwGradU);
- term2[0].add(uwGradV[0]);
- term2[1].add(uwGradV[1]);
- term2[2].add(uwGradV[2]);
-
- term2[0].sub(uvGradW[0]);
- term2[1].sub(uvGradW[1]);
- term2[2].sub(uvGradW[2]);
-
- // Preliminary Type 3
- vector<Polynomial> gradL1 = lagrange[0].gradient();
- vector<Polynomial> gradL2 = lagrange[1].gradient();
-
- vector<Polynomial> l2GradL1(gradL1);
- l2GradL1[0].mul(lagrange[1]);
- l2GradL1[1].mul(lagrange[1]);
- l2GradL1[2].mul(lagrange[1]);
-
- vector<Polynomial> l1GradL2(gradL2);
- l1GradL2[0].mul(lagrange[0]);
- l1GradL2[1].mul(lagrange[0]);
- l1GradL2[2].mul(lagrange[0]);
-
- vector<Polynomial> subGradL1L2VW(l2GradL1);
- subGradL1L2VW[0].sub(l1GradL2[0]);
- subGradL1L2VW[1].sub(l1GradL2[1]);
- subGradL1L2VW[2].sub(l1GradL2[2]);
-
- subGradL1L2VW[0].mul(v);
- subGradL1L2VW[1].mul(v);
- subGradL1L2VW[2].mul(v);
-
- subGradL1L2VW[0].mul(w);
- subGradL1L2VW[1].mul(w);
- subGradL1L2VW[2].mul(w);
-
-
- // Type 1
- basis[s][i] = new vector<Polynomial>((u * v * w).gradient());
- i++;
-
- // Type 2 -- Part 1
- basis[s][i] = new vector<Polynomial>(term1);
- i++;
-
- // Type 2 -- Part 2
- basis[s][i] = new vector<Polynomial>(term2);
- i++;
-
- // Type 3
- if(l1 == 1){
- basis[s][i] = new vector<Polynomial>(subGradL1L2VW);
- i++;
- }
- }
+ for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusTwo; l3++){
+ // Preliminary Type 1
+ Polynomial u = intLegendre[l1].compose
+ (lagrange[0] - lagrange[1],
+ lagrange[0] + lagrange[1]);
+
+ Polynomial v = lagrange[2] * sclLegendre[l2].compose
+ (lagrange[2] * 2 - (one - lagrange[3]), one - lagrange[3]);
+
+ Polynomial w = lagrange[3] * legendre[l3].compose
+ (lagrange[3] * 2 - one);
+
+ // Preliminary Type 2
+ vector<Polynomial> gradU = u.gradient();
+ vector<Polynomial> gradV = v.gradient();
+ vector<Polynomial> gradW = w.gradient();
+
+ vector<Polynomial> vwGradU(gradU);
+ vwGradU[0].mul(v);
+ vwGradU[1].mul(v);
+ vwGradU[2].mul(v);
+
+ vwGradU[0].mul(w);
+ vwGradU[1].mul(w);
+ vwGradU[2].mul(w);
+
+ vector<Polynomial> uwGradV(gradV);
+ uwGradV[0].mul(u);
+ uwGradV[1].mul(u);
+ uwGradV[2].mul(u);
+
+ uwGradV[0].mul(w);
+ uwGradV[1].mul(w);
+ uwGradV[2].mul(w);
+
+ vector<Polynomial> uvGradW(gradW);
+ uvGradW[0].mul(u);
+ uvGradW[1].mul(u);
+ uvGradW[2].mul(u);
+
+ uvGradW[0].mul(v);
+ uvGradW[1].mul(v);
+ uvGradW[2].mul(v);
+
+ vector<Polynomial> term1(vwGradU);
+ term1[0].sub(uwGradV[0]);
+ term1[1].sub(uwGradV[1]);
+ term1[2].sub(uwGradV[2]);
+
+ term1[0].add(uvGradW[0]);
+ term1[1].add(uvGradW[1]);
+ term1[2].add(uvGradW[2]);
+
+ vector<Polynomial> term2(vwGradU);
+ term2[0].add(uwGradV[0]);
+ term2[1].add(uwGradV[1]);
+ term2[2].add(uwGradV[2]);
+
+ term2[0].sub(uvGradW[0]);
+ term2[1].sub(uvGradW[1]);
+ term2[2].sub(uvGradW[2]);
+
+ // Preliminary Type 3
+ vector<Polynomial> gradL1 = lagrange[0].gradient();
+ vector<Polynomial> gradL2 = lagrange[1].gradient();
+
+ vector<Polynomial> l2GradL1(gradL1);
+ l2GradL1[0].mul(lagrange[1]);
+ l2GradL1[1].mul(lagrange[1]);
+ l2GradL1[2].mul(lagrange[1]);
+
+ vector<Polynomial> l1GradL2(gradL2);
+ l1GradL2[0].mul(lagrange[0]);
+ l1GradL2[1].mul(lagrange[0]);
+ l1GradL2[2].mul(lagrange[0]);
+
+ vector<Polynomial> subGradL1L2VW(l2GradL1);
+ subGradL1L2VW[0].sub(l1GradL2[0]);
+ subGradL1L2VW[1].sub(l1GradL2[1]);
+ subGradL1L2VW[2].sub(l1GradL2[2]);
+
+ subGradL1L2VW[0].mul(v);
+ subGradL1L2VW[1].mul(v);
+ subGradL1L2VW[2].mul(v);
+
+ subGradL1L2VW[0].mul(w);
+ subGradL1L2VW[1].mul(w);
+ subGradL1L2VW[2].mul(w);
+
+
+ // Type 1
+ basis[s][i] = new vector<Polynomial>((u * v * w).gradient());
+ i++;
+
+ // Type 2 -- Part 1
+ basis[s][i] = new vector<Polynomial>(term1);
+ i++;
+
+ // Type 2 -- Part 2
+ basis[s][i] = new vector<Polynomial>(term2);
+ i++;
+
+ // Type 3
+ if(l1 == 1){
+ basis[s][i] = new vector<Polynomial>(subGradL1L2VW);
+ i++;
+ }
+ }
}
}
}
diff --git a/FunctionSpace/TetNedelecBasis.cpp b/FunctionSpace/TetNedelecBasis.cpp
index a2cf5ecc30..be15c1a6d9 100644
--- a/FunctionSpace/TetNedelecBasis.cpp
+++ b/FunctionSpace/TetNedelecBasis.cpp
@@ -6,7 +6,7 @@ using namespace std;
TetNedelecBasis::TetNedelecBasis(void){
// Reference Space //
refSpace = new TetReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -27,9 +27,9 @@ TetNedelecBasis::TetNedelecBasis(void){
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(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0) -
+ Polynomial(1, 0, 0, 1)),
Polynomial(Polynomial(1, 1, 0, 0)),
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 27cb84d6ef..cd221914e9 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
TetNodeBasis::TetNodeBasis(unsigned int order){
// Reference Space //
refSpace = new TetReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -45,9 +45,9 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
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(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0) -
+ Polynomial(1, 0, 0, 1)),
Polynomial(Polynomial(1, 1, 0, 0)),
@@ -77,14 +77,14 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
for(int e = 0; e < 6; e++){
for(unsigned int l = 1; l < order; l++){
- basis[s][i] =
- new Polynomial(intLegendre[l].compose
- (lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
- ,
- lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]]));
- i++;
+ basis[s][i] =
+ new Polynomial(intLegendre[l].compose
+ (lagrange[(*(*edgeV[s])[e])[0]] -
+ lagrange[(*(*edgeV[s])[e])[1]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[0]] +
+ lagrange[(*(*edgeV[s])[e])[1]]));
+ i++;
}
}
}
@@ -97,28 +97,28 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
for(int f = 0; f < 4; f++){
for(int l1 = 1; l1 < orderMinus; l1++){
- for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
- Polynomial sum =
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]] +
- lagrange[(*(*faceV[s])[f])[2]];
-
- basis[s][i] =
- new Polynomial(intLegendre[l1].compose
- (lagrange[(*(*faceV[s])[f])[0]] -
- lagrange[(*(*faceV[s])[f])[1]]
- ,
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]])
-
- *
-
- lagrange[(*(*faceV[s])[f])[2]]
- *
- sclLegendre[l2].compose
- (lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum));
- i++;
- }
+ for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
+ Polynomial sum =
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]] +
+ lagrange[(*(*faceV[s])[f])[2]];
+
+ basis[s][i] =
+ new Polynomial(intLegendre[l1].compose
+ (lagrange[(*(*faceV[s])[f])[0]] -
+ lagrange[(*(*faceV[s])[f])[1]]
+ ,
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]])
+
+ *
+
+ lagrange[(*(*faceV[s])[f])[2]]
+ *
+ sclLegendre[l2].compose
+ (lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum));
+ i++;
+ }
}
}
}
@@ -136,16 +136,16 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
for(int l1 = 1; l1 < orderMinusTwo; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
- for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
- basis[s][i] =
- new Polynomial(intLegendre[l1].compose(sub, add) *
- lagrange[2] *
- sclLegendre[l2].compose(twoThreeOneMinusFour,
- oneMinusFour) *
- lagrange[3] *
- legendre[l3].compose(twoFourMinusOne));
- i++;
- }
+ for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
+ basis[s][i] =
+ new Polynomial(intLegendre[l1].compose(sub, add) *
+ lagrange[2] *
+ sclLegendre[l2].compose(twoThreeOneMinusFour,
+ oneMinusFour) *
+ lagrange[3] *
+ legendre[l3].compose(twoFourMinusOne));
+ i++;
+ }
}
}
}
diff --git a/FunctionSpace/TetReferenceSpace.cpp b/FunctionSpace/TetReferenceSpace.cpp
index 0c214a408c..154b2a9e73 100644
--- a/FunctionSpace/TetReferenceSpace.cpp
+++ b/FunctionSpace/TetReferenceSpace.cpp
@@ -10,10 +10,10 @@ TetReferenceSpace::TetReferenceSpace(void){
// Edge Definition //
nEdge = 6;
- refEdge = new unsigned int*[nEdge];
+ refEdge = new size_t*[nEdge];
- for(unsigned int i = 0; i < nEdge; i++){
- refEdge[i] = new unsigned int[2];
+ for(size_t i = 0; i < nEdge; i++){
+ refEdge[i] = new size_t[2];
refEdge[i][0] = MTetrahedron::edges_tetra(i, 0);
refEdge[i][1] = MTetrahedron::edges_tetra(i, 1);
}
@@ -23,16 +23,16 @@ TetReferenceSpace::TetReferenceSpace(void){
nFace = 4;
// Number of node per face
- nNodeInFace = new unsigned int[nFace];
+ nNodeInFace = new size_t[nFace];
- for(unsigned int f = 0; f < nFace; f++)
+ for(size_t f = 0; f < nFace; f++)
nNodeInFace[f] = 3;
// Reference Face
- refFace = new unsigned int*[nFace];
+ refFace = new size_t*[nFace];
- for(unsigned int i = 0; i < nFace; i++){
- refFace[i] = new unsigned int[3];
+ for(size_t i = 0; i < nFace; i++){
+ refFace[i] = new size_t[3];
refFace[i][0] = MTetrahedron::faces_tetra(i, 0);
refFace[i][1] = MTetrahedron::faces_tetra(i, 1);
refFace[i][2] = MTetrahedron::faces_tetra(i, 2);
@@ -44,13 +44,13 @@ TetReferenceSpace::TetReferenceSpace(void){
TetReferenceSpace::~TetReferenceSpace(void){
// Delete Ref Edge //
- for(unsigned int i = 0; i < nEdge; i++)
+ for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
delete[] refEdge;
// Delete Ref Face //
- for(unsigned int i = 0; i < nFace; i++)
+ for(size_t i = 0; i < nFace; i++)
delete[] refFace[i];
delete[] refFace;
@@ -58,35 +58,39 @@ TetReferenceSpace::~TetReferenceSpace(void){
}
string TetReferenceSpace::toLatex(void) const{
- stringstream stream;
+ const size_t nPerm = pTree->getNPermutation();
+ stringstream stream;
+ vector<size_t> perm;
stream << "\\documentclass{article}" << endl << endl
- << "\\usepackage{longtable}" << endl
- << "\\usepackage{tikz}" << endl
- << "\\usetikzlibrary{arrows}" << endl << endl
+ << "\\usepackage{longtable}" << endl
+ << "\\usepackage{tikz}" << endl
+ << "\\usetikzlibrary{arrows}" << endl << endl
- << "\\begin{document}" << endl
- << "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
- << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
+ << "\\begin{document}" << endl
+ << "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
+ << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
- << "\\begin{longtable}{ccc}" << endl << endl;
+ << "\\begin{longtable}{ccc}" << endl << endl;
+
+ for(size_t p = 0; p < nPerm; p++){
+ pTree->fillWithPermutation(p, perm);
- for(unsigned int p = 0; p < nPerm; p++){
stream << "\\begin{tikzpicture}" << endl
- << "\\node[vertex] (n0) at(0, 0) {$" << perm[p][0] << "$};" << endl
- << "\\node[vertex] (n1) at(3, 0) {$" << perm[p][1] << "$};" << endl
- << "\\node[vertex] (n2) at(0, 3) {$" << perm[p][2] << "$};" << endl
- << "\\node[vertex] (n3) at(1, 1) {$" << perm[p][3] << "$};" << endl
+ << "\\node[vertex] (n0) at(0, 0) {$" << perm[0] << "$};" << endl
+ << "\\node[vertex] (n1) at(3, 0) {$" << perm[1] << "$};" << endl
+ << "\\node[vertex] (n2) at(0, 3) {$" << perm[2] << "$};" << endl
+ << "\\node[vertex] (n3) at(1, 1) {$" << perm[3] << "$};" << endl
<< endl;
- for(unsigned int i = 0; i < 6; i++)
+ for(size_t i = 0; i < 6; i++)
stream << "\\path[line]"
- << " (n" << (*(*(*edge)[p])[i])[0] << ")"
- << " -- "
- << " (n" << (*(*(*edge)[p])[i])[1] << ");"
- << endl;
+ << " (n" << (*(*(*edge)[p])[i])[0] << ")"
+ << " -- "
+ << " (n" << (*(*(*edge)[p])[i])[1] << ");"
+ << endl;
if((p + 1) % 3)
stream << "\\end{tikzpicture} & " << endl << endl;
@@ -96,7 +100,7 @@ string TetReferenceSpace::toLatex(void) const{
}
stream << "\\end{longtable}" << endl
- << "\\end{document}" << endl;
+ << "\\end{document}" << endl;
return stream.str();
}
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 088215ee12..251807a1c5 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Reference Space //
refSpace = new TriReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -42,8 +42,8 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
const Polynomial lagrange[3] =
{
Polynomial(Polynomial(1, 0, 0, 0) -
- Polynomial(1, 1, 0, 0) -
- Polynomial(1, 0, 1, 0)),
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0)),
Polynomial(Polynomial(1, 1, 0, 0)),
@@ -62,38 +62,38 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
for(int e = 0; e < 3; e++){
for(int l = 0; l < orderPlus; l++){
- // Nedelec
- if(l == 0){
- 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][i] = new vector<Polynomial>(tmp2);
- }
-
- // High Order
- else{
- basis[s][i] =
- new vector<Polynomial>
- ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
- ,
- lagrange[(*(*edgeV[s])[e])[1]] +
- lagrange[(*(*edgeV[s])[e])[0]])).gradient());
- }
- i++;
+ // Nedelec
+ if(l == 0){
+ 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][i] = new vector<Polynomial>(tmp2);
+ }
+
+ // High Order
+ else{
+ basis[s][i] =
+ new vector<Polynomial>
+ ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
+ lagrange[(*(*edgeV[s])[e])[1]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[1]] +
+ lagrange[(*(*edgeV[s])[e])[0]])).gradient());
+ }
+ i++;
}
}
}
@@ -160,48 +160,48 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Type 1
for(unsigned int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
- vector<Polynomial> tmp1 = v[s][l2].gradient();
- vector<Polynomial> tmp2 = u[s][l1].gradient();
+ vector<Polynomial> tmp1 = v[s][l2].gradient();
+ vector<Polynomial> tmp2 = u[s][l1].gradient();
- tmp1[0].mul(u[s][l1]);
- tmp1[1].mul(u[s][l1]);
- tmp1[2].mul(u[s][l1]);
+ tmp1[0].mul(u[s][l1]);
+ tmp1[1].mul(u[s][l1]);
+ tmp1[2].mul(u[s][l1]);
- tmp2[0].mul(v[s][l2]);
- tmp2[1].mul(v[s][l2]);
- tmp2[2].mul(v[s][l2]);
+ tmp2[0].mul(v[s][l2]);
+ tmp2[1].mul(v[s][l2]);
+ tmp2[2].mul(v[s][l2]);
- tmp2[0].add(tmp1[0]);
- tmp2[1].add(tmp1[1]);
- tmp2[2].add(tmp1[2]);
+ tmp2[0].add(tmp1[0]);
+ tmp2[1].add(tmp1[1]);
+ tmp2[2].add(tmp1[2]);
- basis[s][i] = new vector<Polynomial>(tmp2);
+ basis[s][i] = new vector<Polynomial>(tmp2);
- i++;
+ i++;
}
}
// Type 2
for(unsigned int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
- vector<Polynomial> tmp1 = v[s][l2].gradient();
- vector<Polynomial> tmp2 = u[s][l1].gradient();
+ vector<Polynomial> tmp1 = v[s][l2].gradient();
+ vector<Polynomial> tmp2 = u[s][l1].gradient();
- tmp1[0].mul(u[s][l1]);
- tmp1[1].mul(u[s][l1]);
- tmp1[2].mul(u[s][l1]);
+ tmp1[0].mul(u[s][l1]);
+ tmp1[1].mul(u[s][l1]);
+ tmp1[2].mul(u[s][l1]);
- tmp2[0].mul(v[s][l2]);
- tmp2[1].mul(v[s][l2]);
- tmp2[2].mul(v[s][l2]);
+ tmp2[0].mul(v[s][l2]);
+ tmp2[1].mul(v[s][l2]);
+ tmp2[2].mul(v[s][l2]);
- tmp2[0].sub(tmp1[0]);
- tmp2[1].sub(tmp1[1]);
- tmp2[2].sub(tmp1[2]);
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
- basis[s][i] = new vector<Polynomial>(tmp2);
+ basis[s][i] = new vector<Polynomial>(tmp2);
- i++;
+ i++;
}
}
diff --git a/FunctionSpace/TriLagrangeReferenceSpace.cpp b/FunctionSpace/TriLagrangeReferenceSpace.cpp
index d5e8772c69..b3c376f654 100644
--- a/FunctionSpace/TriLagrangeReferenceSpace.cpp
+++ b/FunctionSpace/TriLagrangeReferenceSpace.cpp
@@ -9,18 +9,18 @@ TriLagrangeReferenceSpace::TriLagrangeReferenceSpace(unsigned int order){
// Edge Definition //
nEdge = 3;
- refEdge = new unsigned int*[nEdge];
+ refEdge = new size_t*[nEdge];
- for(unsigned int i = 0; i < nEdge; i++){
- refEdge[i] = new unsigned int[2];
+ for(size_t i = 0; i < nEdge; i++){
+ refEdge[i] = new size_t[2];
refEdge[i][0] = MTriangle::edges_tri(i, 0);
refEdge[i][1] = MTriangle::edges_tri(i, 1);
}
// Face Definition //
nFace = 1;
- refFace = new unsigned int*[nFace];
- refFace[0] = new unsigned int[3];
+ refFace = new size_t*[nFace];
+ refFace[0] = new size_t[3];
refFace[0][0] = 0;
refFace[0][1] = 1;
@@ -45,13 +45,13 @@ TriLagrangeReferenceSpace::TriLagrangeReferenceSpace(unsigned int order){
TriLagrangeReferenceSpace::~TriLagrangeReferenceSpace(void){
// Delete Ref Edge //
- for(unsigned int i = 0; i < nEdge; i++)
+ for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
delete[] refEdge;
// Delete Ref Face //
- for(unsigned int i = 0; i < nFace; i++)
+ for(size_t i = 0; i < nFace; i++)
delete[] refFace[i];
delete[] refFace;
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index e16146f436..718ede16f2 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -6,7 +6,7 @@ using namespace std;
TriNedelecBasis::TriNedelecBasis(void){
// Reference Space //
refSpace = new TriReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -27,8 +27,8 @@ TriNedelecBasis::TriNedelecBasis(void){
const Polynomial lagrange[3] =
{
Polynomial(Polynomial(1, 0, 0, 0) -
- Polynomial(1, 1, 0, 0) -
- Polynomial(1, 0, 1, 0)),
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0)),
Polynomial(Polynomial(1, 1, 0, 0)),
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index b6f1604b7d..957722af4d 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -7,7 +7,7 @@ using namespace std;
TriNodeBasis::TriNodeBasis(unsigned int order){
// Reference Space //
refSpace = new TriReferenceSpace;
- nRefSpace = refSpace->getNPermutation();
+ nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
@@ -40,8 +40,8 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
const Polynomial lagrange[3] =
{
Polynomial(Polynomial(1, 0, 0, 0) -
- Polynomial(1, 1, 0, 0) -
- Polynomial(1, 0, 1, 0)),
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0)),
Polynomial(Polynomial(1, 1, 0, 0)),
@@ -67,13 +67,13 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
for(unsigned int e = 0; e < 3; e++){
for(unsigned int l = 1; l < order; l++){
- basis[s][i] =
- new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
- lagrange[(*(*edgeV[s])[e])[0]]
- ,
- lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]]));
- i++;
+ basis[s][i] =
+ new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
+ lagrange[(*(*edgeV[s])[e])[0]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[0]] +
+ lagrange[(*(*edgeV[s])[e])[1]]));
+ i++;
}
}
}
@@ -90,21 +90,21 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
- basis[s][i] =
- new Polynomial(intLegendre[l1].compose(lagrange[(*(*faceV[s])[0])[1]] -
+ basis[s][i] =
+ new Polynomial(intLegendre[l1].compose(lagrange[(*(*faceV[s])[0])[1]] -
lagrange[(*(*faceV[s])[0])[0]]
,
- lagrange[(*(*faceV[s])[0])[0]] +
+ lagrange[(*(*faceV[s])[0])[0]] +
lagrange[(*(*faceV[s])[0])[1]])
- *
+ *
- legendre[l2].compose((lagrange[(*(*faceV[s])[0])[2]] * 2)
+ legendre[l2].compose((lagrange[(*(*faceV[s])[0])[2]] * 2)
-
Polynomial(1, 0, 0, 0))
*
lagrange[(*(*faceV[s])[0])[2]]);
- i++;
+ i++;
}
}
}
diff --git a/FunctionSpace/TriReferenceSpace.cpp b/FunctionSpace/TriReferenceSpace.cpp
index 681f94f1e7..b976bc43ed 100644
--- a/FunctionSpace/TriReferenceSpace.cpp
+++ b/FunctionSpace/TriReferenceSpace.cpp
@@ -10,10 +10,10 @@ TriReferenceSpace::TriReferenceSpace(void){
// Edge Definition //
nEdge = 3;
- refEdge = new unsigned int*[nEdge];
+ refEdge = new size_t*[nEdge];
- for(unsigned int i = 0; i < nEdge; i++){
- refEdge[i] = new unsigned int[2];
+ for(size_t i = 0; i < nEdge; i++){
+ refEdge[i] = new size_t[2];
refEdge[i][0] = MTriangle::edges_tri(i, 0);
refEdge[i][1] = MTriangle::edges_tri(i, 1);
}
@@ -23,12 +23,12 @@ TriReferenceSpace::TriReferenceSpace(void){
nFace = 1;
// Number of node per face
- nNodeInFace = new unsigned int[nFace];
+ nNodeInFace = new size_t[nFace];
nNodeInFace[0] = 3;
// Reference Face
- refFace = new unsigned int*[nFace];
- refFace[0] = new unsigned int[3];
+ refFace = new size_t*[nFace];
+ refFace[0] = new size_t[3];
refFace[0][0] = 0;
refFace[0][1] = 1;
@@ -40,13 +40,13 @@ TriReferenceSpace::TriReferenceSpace(void){
TriReferenceSpace::~TriReferenceSpace(void){
// Delete Ref Edge //
- for(unsigned int i = 0; i < nEdge; i++)
+ for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
delete[] refEdge;
// Delete Ref Face //
- for(unsigned int i = 0; i < nFace; i++)
+ for(size_t i = 0; i < nFace; i++)
delete[] refFace[i];
delete[] refFace;
@@ -54,7 +54,9 @@ TriReferenceSpace::~TriReferenceSpace(void){
}
string TriReferenceSpace::toLatex(void) const{
- stringstream stream;
+ const size_t nPerm = pTree->getNPermutation();
+ stringstream stream;
+ vector<size_t> perm(nVertex);
stream << "\\documentclass{article}" << endl << endl
@@ -68,15 +70,17 @@ string TriReferenceSpace::toLatex(void) const{
<< "\\begin{longtable}{ccc}" << endl << endl;
- for(unsigned int p = 0; p < nPerm; p++){
+ for(size_t p = 0; p < nPerm; p++){
+ pTree->fillWithPermutation(p, perm);
+
stream << "\\begin{tikzpicture}" << endl
- << "\\node[vertex] (n0) at(0, 0) {$" << perm[p][0] << "$};" << endl
- << "\\node[vertex] (n1) at(3, 0) {$" << perm[p][1] << "$};" << endl
- << "\\node[vertex] (n2) at(0, 3) {$" << perm[p][2] << "$};" << endl
+ << "\\node[vertex] (n0) at(0, 0) {$" << perm[0] << "$};" << endl
+ << "\\node[vertex] (n1) at(3, 0) {$" << perm[1] << "$};" << endl
+ << "\\node[vertex] (n2) at(0, 3) {$" << perm[2] << "$};" << endl
<< endl;
- for(unsigned int i = 0; i < 3; i++)
+ for(size_t i = 0; i < 3; i++)
stream << "\\path[line]"
<< " (n" << (*(*(*edge)[p])[i])[0] << ")"
<< " -- "
--
GitLab