From 1e704e72c6a598d4e85bc645362dcee3b28eaf74 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 16 Jul 2013 13:56:10 +0000
Subject: [PATCH] ** ReferenceSpace: -- Minimal number of reference space
(at least for tet, tri and quad) -- Automatic generation of these
reference space -- New interface -- Modification of TriNodeBasis
accordingly
** Need to compute Dof permutations automaticaly
** All Basis are commented (exept for TriNodeBasis)
** WARNING: Memory leak somewhere
---
FunctionSpace/Basis.h | 4 +
FunctionSpace/BasisHierarchical0From.cpp | 33 ++
FunctionSpace/BasisHierarchical0From.h | 4 +
FunctionSpace/BasisHierarchical1From.cpp | 6 +
FunctionSpace/BasisHierarchical1From.h | 4 +
FunctionSpace/BasisLagrange.cpp | 5 +
FunctionSpace/BasisLagrange.h | 4 +
FunctionSpace/FunctionSpaceScalar.cpp | 17 +-
FunctionSpace/LineEdgeBasis.cpp | 4 +
FunctionSpace/LineNedelecBasis.cpp | 4 +
FunctionSpace/LineNodeBasis.cpp | 4 +
FunctionSpace/LineReferenceSpace.cpp | 10 +-
FunctionSpace/QuadEdgeBasis.cpp | 4 +
FunctionSpace/QuadNedelecBasis.cpp | 4 +
FunctionSpace/QuadNodeBasis.cpp | 4 +
FunctionSpace/QuadReferenceSpace.cpp | 12 +-
FunctionSpace/ReferenceSpace.cpp | 499 +++++++++-----------
FunctionSpace/ReferenceSpace.h | 143 +++---
FunctionSpace/ReferenceSpaceLagrange.cpp | 2 +
FunctionSpace/TetEdgeBasis.cpp | 4 +
FunctionSpace/TetNedelecBasis.cpp | 4 +
FunctionSpace/TetNodeBasis.cpp | 4 +
FunctionSpace/TetReferenceSpace.cpp | 12 +-
FunctionSpace/TriEdgeBasis.cpp | 4 +
FunctionSpace/TriLagrangeReferenceSpace.cpp | 4 +
FunctionSpace/TriNedelecBasis.cpp | 4 +
FunctionSpace/TriNodeBasis.cpp | 52 +-
FunctionSpace/TriNodeBasis.h | 2 +-
FunctionSpace/TriReferenceSpace.cpp | 64 +--
29 files changed, 486 insertions(+), 435 deletions(-)
diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 7738a321f3..6cd443357a 100644
--- a/FunctionSpace/Basis.h
+++ b/FunctionSpace/Basis.h
@@ -81,6 +81,10 @@ class Basis{
virtual unsigned int getNOrientation(void) const = 0;
virtual unsigned int getOrientation(const MElement& element) const = 0;
+ virtual void mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]) const = 0;
+
// Functions Ordering //
virtual std::vector<size_t>
getFunctionOrdering(const MElement& element) const = 0;
diff --git a/FunctionSpace/BasisHierarchical0From.cpp b/FunctionSpace/BasisHierarchical0From.cpp
index 26fd0ce1e0..cc3a579c4c 100644
--- a/FunctionSpace/BasisHierarchical0From.cpp
+++ b/FunctionSpace/BasisHierarchical0From.cpp
@@ -59,6 +59,12 @@ getOrientation(const MElement& element) const{
return refSpace->getReferenceSpace(element);
}
+void BasisHierarchical0From::mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]) const{
+ return refSpace->mapFromXYZtoABC(element, xyz, abc);
+}
+
vector<size_t> BasisHierarchical0From::
getFunctionOrdering(const MElement& element) const{
vector<size_t> ordering(nFunction);
@@ -66,6 +72,33 @@ getFunctionOrdering(const MElement& element) const{
for(size_t i = 0; i < nFunction; i++)
ordering[i] = i;
+ /*
+ if(element.getNum() == 5){
+ ordering[0] = 0;
+ ordering[1] = 1;
+ ordering[2] = 2;
+ ordering[3] = 3;
+ ordering[4] = 4;
+ ordering[5] = 5;
+ ordering[6] = 6;
+ ordering[7] = 7;
+ ordering[8] = 8;
+ ordering[9] = 9;
+ }
+
+ else{
+ ordering[0] = 2;
+ ordering[1] = 0;
+ ordering[2] = 1;
+ ordering[3] = 7;
+ ordering[4] = 8;
+ ordering[5] = 3;
+ ordering[6] = 4;
+ ordering[7] = 5;
+ ordering[8] = 6;
+ ordering[9] = 9;
+ }
+ */
return ordering;
}
diff --git a/FunctionSpace/BasisHierarchical0From.h b/FunctionSpace/BasisHierarchical0From.h
index 97d0cba28c..96377520d7 100644
--- a/FunctionSpace/BasisHierarchical0From.h
+++ b/FunctionSpace/BasisHierarchical0From.h
@@ -39,6 +39,10 @@ class BasisHierarchical0From: public BasisLocal{
virtual unsigned int getNOrientation(void) const;
virtual unsigned int getOrientation(const MElement& element) const;
+ virtual void mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]) const;
+
virtual std::vector<size_t>
getFunctionOrdering(const MElement& element) const;
diff --git a/FunctionSpace/BasisHierarchical1From.cpp b/FunctionSpace/BasisHierarchical1From.cpp
index 7fe070a919..fa94def81b 100644
--- a/FunctionSpace/BasisHierarchical1From.cpp
+++ b/FunctionSpace/BasisHierarchical1From.cpp
@@ -59,6 +59,12 @@ getOrientation(const MElement& element) const{
return refSpace->getReferenceSpace(element);
}
+void BasisHierarchical1From::mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]) const{
+ return refSpace->mapFromXYZtoABC(element, xyz, abc);
+}
+
vector<size_t> BasisHierarchical1From::
getFunctionOrdering(const MElement& element) const{
vector<size_t> ordering(nFunction);
diff --git a/FunctionSpace/BasisHierarchical1From.h b/FunctionSpace/BasisHierarchical1From.h
index 3057837949..f72d061ef5 100644
--- a/FunctionSpace/BasisHierarchical1From.h
+++ b/FunctionSpace/BasisHierarchical1From.h
@@ -39,6 +39,10 @@ class BasisHierarchical1From: public BasisLocal{
virtual unsigned int getNOrientation(void) const;
virtual unsigned int getOrientation(const MElement& element) const;
+ virtual void mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]) const;
+
virtual std::vector<size_t>
getFunctionOrdering(const MElement& element) const;
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index dcdb1e85e2..592af30b3c 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -78,6 +78,11 @@ static size_t matchClosure(vector<int>& reduced,
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{
/*
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index b82487fc41..23cd49da56 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/FunctionSpace/BasisLagrange.h
@@ -43,6 +43,10 @@ class BasisLagrange: public BasisLocal{
virtual unsigned int getNOrientation(void) const;
virtual unsigned int getOrientation(const MElement& element) const;
+ virtual void mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]) const;
+
virtual std::vector<size_t>
getFunctionOrdering(const MElement& element) const;
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 1400403da7..8e888e9c49 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -14,7 +14,7 @@ double FunctionSpaceScalar::
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
-
+ /*
// Const Cast For MElement //
MElement& eelement =
const_cast<MElement&>(element);
@@ -34,11 +34,11 @@ interpolate(const MElement& element,
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];
@@ -50,14 +50,14 @@ interpolate(const MElement& element,
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);
@@ -90,12 +90,17 @@ interpolate(const MElement& element,
uvw[0] = abc[0];
uvw[1] = abc[1];
uvw[2] = abc[2];
+ */
+
+ // Get ABC Space coordinate //
+ double abc[3];
+ (*basis)[0]->mapFromXYZtoABC(element, xyz, abc);
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
fullMatrix<double> fun(nFun, 1);
- (*basis)[0]->getFunctions(fun, element, uvw[0], uvw[1], uvw[2]);
+ (*basis)[0]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
// Interpolate (in Reference Place) //
double val = 0;
diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp
index 933d162ff1..2371c893d3 100644
--- a/FunctionSpace/LineEdgeBasis.cpp
+++ b/FunctionSpace/LineEdgeBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
LineEdgeBasis::LineEdgeBasis(unsigned int order){
+ /*
// Reference Space //
refSpace = new LineReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -71,9 +72,11 @@ LineEdgeBasis::LineEdgeBasis(unsigned int order){
// Free Temporary Space //
delete[] intLegendre;
+ */
}
LineEdgeBasis::~LineEdgeBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -86,4 +89,5 @@ LineEdgeBasis::~LineEdgeBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp
index 0e3727d7ef..8de25c47bd 100644
--- a/FunctionSpace/LineNedelecBasis.cpp
+++ b/FunctionSpace/LineNedelecBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
LineNedelecBasis::LineNedelecBasis(void){
+ /*
// Reference Space //
refSpace = new LineReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -41,9 +42,11 @@ LineNedelecBasis::LineNedelecBasis(void){
// Nedelec //
basis[0][0] = new vector<Polynomial>(first);
basis[1][0] = new vector<Polynomial>(second);
+ */
}
LineNedelecBasis::~LineNedelecBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -56,4 +59,5 @@ LineNedelecBasis::~LineNedelecBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index 3c00f46ba1..1081cb263a 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
LineNodeBasis::LineNodeBasis(unsigned int order){
+ /*
// Reference Space //
refSpace = new LineReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -66,9 +67,11 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
// Free Temporary Sapce //
delete[] intLegendre;
+ */
}
LineNodeBasis::~LineNodeBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -81,4 +84,5 @@ LineNodeBasis::~LineNodeBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/LineReferenceSpace.cpp b/FunctionSpace/LineReferenceSpace.cpp
index 9251fdec6e..3e4b321556 100644
--- a/FunctionSpace/LineReferenceSpace.cpp
+++ b/FunctionSpace/LineReferenceSpace.cpp
@@ -4,6 +4,7 @@
using namespace std;
LineReferenceSpace::LineReferenceSpace(void){
+ /*
// Vertex Definition //
nVertex = 2;
@@ -21,20 +22,23 @@ LineReferenceSpace::LineReferenceSpace(void){
// Init All //
init();
+ */
}
LineReferenceSpace::~LineReferenceSpace(void){
+ /*
// Delete Ref Edge //
for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
delete[] refEdge;
+ */
}
string LineReferenceSpace::toLatex(void) const{
- const size_t nPerm = pTree->getNPermutation();
+ //const size_t nPerm = pTree->getNPermutation();
stringstream stream;
-
+ /*
stream << "\\documentclass{article}" << endl << endl
<< "\\usepackage{longtable}" << endl
@@ -64,6 +68,6 @@ string LineReferenceSpace::toLatex(void) const{
stream << "\\end{longtable}" << endl
<< "\\end{document}" << endl;
-
+ */
return stream.str();
}
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 86c2a2c51f..5fb906d080 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
+ /*
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -238,9 +239,11 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
+ */
}
QuadEdgeBasis::~QuadEdgeBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -253,4 +256,5 @@ QuadEdgeBasis::~QuadEdgeBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/QuadNedelecBasis.cpp b/FunctionSpace/QuadNedelecBasis.cpp
index 7e4d245adb..42eda0550c 100644
--- a/FunctionSpace/QuadNedelecBasis.cpp
+++ b/FunctionSpace/QuadNedelecBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
QuadNedelecBasis::QuadNedelecBasis(void){
+ /*
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -105,9 +106,11 @@ QuadNedelecBasis::QuadNedelecBasis(void){
delete old;
}
}
+ */
}
QuadNedelecBasis::~QuadNedelecBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -120,4 +123,5 @@ QuadNedelecBasis::~QuadNedelecBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 7034d5b9db..1e3c333d89 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
QuadNodeBasis::QuadNodeBasis(unsigned int order){
+ /*
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -142,9 +143,11 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
// Free Temporary Sapce //
delete[] legendre;
+ */
}
QuadNodeBasis::~QuadNodeBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -157,4 +160,5 @@ QuadNodeBasis::~QuadNodeBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/QuadReferenceSpace.cpp b/FunctionSpace/QuadReferenceSpace.cpp
index 88e42b67b3..111ed806e8 100644
--- a/FunctionSpace/QuadReferenceSpace.cpp
+++ b/FunctionSpace/QuadReferenceSpace.cpp
@@ -5,6 +5,7 @@
using namespace std;
QuadReferenceSpace::QuadReferenceSpace(void){
+ /*
// Vertex Definition //
nVertex = 4;
@@ -37,9 +38,11 @@ QuadReferenceSpace::QuadReferenceSpace(void){
// Init All //
init();
+ */
}
QuadReferenceSpace::~QuadReferenceSpace(void){
+ /*
// Delete Ref Edge //
for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
@@ -52,13 +55,14 @@ QuadReferenceSpace::~QuadReferenceSpace(void){
delete[] refFace;
delete[] nNodeInFace;
+ */
}
string QuadReferenceSpace::toLatex(void) const{
- const size_t nPerm = pTree->getNPermutation();
+ //const size_t nPerm = pTree->getNPermutation();
stringstream stream;
- vector<size_t> perm(nVertex);
-
+ //vector<size_t> perm(nVertex);
+ /*
stream << "\\documentclass{article}" << endl << endl
<< "\\usepackage{longtable}" << endl
@@ -98,6 +102,6 @@ string QuadReferenceSpace::toLatex(void) const{
stream << "\\end{longtable}" << endl
<< "\\end{document}" << endl;
-
+ */
return stream.str();
}
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index 110259b920..de56c4278e 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -7,55 +7,12 @@
using namespace std;
ReferenceSpace::ReferenceSpace(void){
- // Init to NULL //
- nEdge = 0;
- refEdge = NULL;
- //permutedRefEdge = NULL;
- edge = NULL;
-
- nFace = 0;
- nNodeInFace = NULL;
- refFace = NULL;
- //permutedRefFace = NULL;
- face = NULL;
-
// Defining Ref Edge and Face in //
// Derived Class //
}
ReferenceSpace::~ReferenceSpace(void){
- const size_t nPerm = pTree->getNPermutation();
-
- // Destroy Tree //
delete pTree;
-
- // Delete Permutated Edge //
- for(size_t p = 0; p < nPerm; p++){
- for(size_t i = 0; i < nEdge; i++){
- //delete[] permutedRefEdge[p][i];
- delete (*(*edge)[p])[i];
- }
-
- //delete[] permutedRefEdge[p];
- delete (*edge)[p];
- }
-
- //delete[] permutedRefEdge;
- delete edge;
-
- // Delete Permutated Face //
- for(size_t p = 0; p < nPerm; p++){
- for(size_t i = 0; i < nFace; i++){
- //delete[] permutedRefFace[p][i];
- delete (*(*face)[p])[i];
- }
-
- //delete[] permutedRefFace[p];
- delete (*face)[p];
- }
-
- //delete[] permutedRefFace;
- delete face;
}
void ReferenceSpace::init(void){
@@ -67,22 +24,51 @@ void ReferenceSpace::init(void){
pTree = new PermutationTree(vertexSeq);
- // Edges & Faces //
- getEdge();
- getFace();
-
// Cyclic Permutation //
- findCyclicPermutation();
+ list<size_t> listOfTrueReferenceSpace;
+ list<vector<size_t> > listOfRefNodeIndexPermutation;
+ list<vector<size_t> > listOfReverseNodeIndexPermutation;
+
+ findCyclicPermutation(listOfTrueReferenceSpace,
+ listOfRefNodeIndexPermutation,
+ listOfReverseNodeIndexPermutation);
+
+ // Iterators //
+ list<size_t>::iterator refSpaceIt = listOfTrueReferenceSpace.begin();
+
+ // Reference Spaces Node Id //
+ const size_t nRefSpace = listOfTrueReferenceSpace.size();
+ refSpaceNodeId.resize(nRefSpace);
+
+ for(size_t i = 0; i < nRefSpace; i++, refSpaceIt++){
+ refSpaceNodeId[i].resize(nVertex);
+
+ pTree->fillWithPermutation(*refSpaceIt, refSpaceNodeId[i]);
+ }
+
+ // ABC space to UVW space shape function index mapping //
+ const size_t nPermutation = pTree->getNPermutation();
+ ABCtoUVWIndex.resize(nPermutation);
+ ABCtoUVWIndex.assign(listOfReverseNodeIndexPermutation.begin(),
+ listOfReverseNodeIndexPermutation.end());
+
+ // UVW space to ABC space shape function index mapping //
+ UVWtoABCIndex.resize(nPermutation);
+ UVWtoABCIndex.assign(listOfRefNodeIndexPermutation.begin(),
+ listOfRefNodeIndexPermutation.end());
+
+ // Ordered Edges & Faces Node Index//
+ getOrderedEdge();
+ getOrderedFace();
}
-void ReferenceSpace::findCyclicPermutation(void){
+void ReferenceSpace::
+findCyclicPermutation(list<size_t>& listOfTrueReferenceSpace,
+ list<vector<size_t> >& listOfRefNodeIndexPermutation,
+ list<vector<size_t> >& listOfReverseNodeIndexPermutation){
// 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);
@@ -90,12 +76,13 @@ void ReferenceSpace::findCyclicPermutation(void){
list<size_t>::iterator end;
triple match;
+ // Values For unPermutedIndex //
vector<size_t> unPermutedIndex(nVertex);
for(size_t i = 0; i < nVertex; i++)
unPermutedIndex[i] = i;
- // Find Cyclic Permutation
+ // Find Cyclic Permutation //
for(size_t i = 0; i < nPerm; i++){
// No match
match.first = false;
@@ -193,19 +180,22 @@ static bool haveSameNode(vector<size_t>& face0,
size_t ReferenceSpace::findCorrespondingFace(vector<size_t>& face,
vector<size_t>& node){
// Init Stuff //
+ const size_t nFace = refFaceNodeIdx.size();
const size_t faceSize = face.size();
bool match = false;
size_t f = 0;
vector<size_t> testFace(faceSize);
+ size_t nNodeInFace;
// Test All Face until match
while(!match && f < nFace){
+ nNodeInFace = refFaceNodeIdx[f].size();
- if(nNodeInFace[f] == faceSize){
+ if(nNodeInFace == faceSize){
// Get face f nodes
for(size_t i = 0; i < faceSize; i++)
- testFace[i] = node[refFace[f][i]];
+ testFace[i] = node[refFaceNodeIdx[f][i]];
// Look if match
match = haveSameNode(testFace, face);
@@ -259,19 +249,21 @@ isCyclicPermutation(vector<size_t>& pTest,
vector<size_t>& pRef){
// Node IDs of Reference Space first Face
- vector<size_t> refNode(nNodeInFace[0]);
+ const size_t nNodeInFaceZero = refFaceNodeIdx[0].size();
+ vector<size_t> refNode(nNodeInFaceZero);
- for(size_t i = 0; i < nNodeInFace[0]; i++)
- refNode[i] = pRef[refFace[0][i]];
+ for(size_t i = 0; i < nNodeInFaceZero; i++)
+ refNode[i] = pRef[refFaceNodeIdx[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]);
+ const size_t nNodeInFaceTest = refFaceNodeIdx[testFaceId].size();
+ vector<size_t> testNode(nNodeInFaceTest);
- for(size_t i = 0; i < nNodeInFace[testFaceId]; i++)
- testNode[i] = pTest[refFace[testFaceId][i]];
+ for(size_t i = 0; i < nNodeInFaceTest; i++)
+ testNode[i] = pTest[refFaceNodeIdx[testFaceId][i]];
// Return Triple //
triple tri = {
@@ -283,229 +275,125 @@ isCyclicPermutation(vector<size_t>& pTest,
return tri;
}
-void ReferenceSpace::getEdge(void){
- const size_t nPerm = pTree->getNPermutation();
+void ReferenceSpace::getOrderedEdge(void){
+ // Fill orderedEdgeNodeIdx[s][e] //
+ // (for all reference space 's' and edge 'e') such that: //
+ // //
+ // refSpaceNodeId[s][orderedEdgeNodeIdx[s][e][0]] < //
+ // refSpaceNodeId[s][orderedEdgeNodeIdx[s][e][1]] //
// Alloc
- vector<const vector<unsigned int>*>* tmp;
- edge = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
- /*
- // All Edges have two nodes
- unsigned int* nNodeInEdge = new unsigned int[nEdge];
-
- for(unsigned int e = 0; e < nEdge; e++)
- nNodeInEdge[e] = 2;
-
- // Get Edge nodes depending on the orientation
- // (edge 1 = [1 2] for orientation 1)
- // ( = [4 1] for orientation 2)
- getPermutedRefEntity(&permutedRefEdge,
- refEdge,
- nNodeInEdge,
- nEdge);
-
- delete[] nNodeInEdge;
- */
+ const size_t nEdge = refEdgeNodeIdx.size();
+ const size_t nRefSpace = refSpaceNodeId.size();
+
+ orderedEdgeNodeIdx.resize(nRefSpace);
+
// Populate Edge
- for(size_t p = 0; p < nPerm; p++){
- tmp = new vector<const vector<unsigned int>*>(nEdge);
+ for(size_t s = 0; s < nRefSpace; s++){
+ orderedEdgeNodeIdx[s].resize(nEdge);
- for(size_t e = 0; e < nEdge; e++)
- (*tmp)[e] = getOrderedEdge(p, e);
+ for(size_t e = 0; e < nEdge; e++){
+ orderedEdgeNodeIdx[s][e].resize(2);
- (*edge)[p] = tmp;
+ orderRefEntityForGivenRefSpace(refEdgeNodeIdx[e],
+ refSpaceNodeId[s],
+ orderedEdgeNodeIdx[s][e]);
+ }
}
}
-void ReferenceSpace::getFace(void){
- const size_t nPerm = pTree->getNPermutation();
+void ReferenceSpace::getOrderedFace(void){
+ // Fill orderedEdgeNodeIdx[s][f] //
+ // (for all reference space 's' and face 'f') such that: //
+ // //
+ // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][0]] < //
+ // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][1]] < //
+ // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][2]] < //
+ // (refSpaceNodeId[s][orderedFaceNodeIdx[s][f][3]]) (Quad Face) //
+ // //
+ // If we have a Quad face, a correction is needed such that: //
+ // //
+ // orderedFaceNodeIdx[s][f][2] //
+ // is *opposite* to //
+ // orderedFaceNodeIdx[s][f][0] //
// Alloc
- vector<const vector<unsigned int>*>* tmp;
- face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
- /*
- // Get Face nodes depending on the orientation
- // (face 1 = [1 2 3 4] for orientation 1)
- // ( = [4 1 2 3] for orientation 2)
- getPermutedRefEntity(&permutedRefFace,
- refFace,
- nNodeInFace,
- nFace);
- */
- // Populate Face
- for(size_t p = 0; p < nPerm; p++){
- tmp = new vector<const vector<unsigned int>*>(nFace);
+ const size_t nFace = refFaceNodeIdx.size();
+ const size_t nRefSpace = refSpaceNodeId.size();
- for(size_t f = 0; f < nFace; f++){
- // Dending on the number of node per face
- // The ordering strategy is different
+ orderedFaceNodeIdx.resize(nRefSpace);
- switch(nNodeInFace[f]){
- case 3: (*tmp)[f] = getOrderedTriFace(p, f); break;
- case 4: (*tmp)[f] = getOrderedQuadFace(p, f); break;
+ // Populate Face
+ for(size_t s = 0; s < nRefSpace; s++){
+ orderedFaceNodeIdx[s].resize(nFace);
- default:
- throw Exception("I can handle face with %d nodes",
- nNodeInFace[f]);
- }
- }
+ for(size_t f = 0; f < nFace; f++){
+ size_t nNodeInFace = refFaceNodeIdx[f].size();
- (*face)[p] = tmp;
- }
-}
-/*
-void ReferenceSpace::getPermutedRefEntity(unsigned int**** permutedRefEntity,
- unsigned int** refEntity,
- unsigned int* nNodeInEntity,
- unsigned int nEntity){
- // Alloc PermutedRefEntity
- (*permutedRefEntity) = new unsigned int**[nPerm];
-
- for(unsigned int p = 0; p < nPerm; p++){
- (*permutedRefEntity)[p] = new unsigned int*[nEntity];
-
- for(unsigned int e = 0; e < nEntity; e++)
- (*permutedRefEntity)[p][e] = new unsigned int[nNodeInEntity[e]];
- }
+ orderedFaceNodeIdx[s][f].resize(nNodeInFace);
- // Generate Permuted Ref Entity
- for(unsigned int p = 0; p < nPerm; p++){
- for(unsigned int e = 0; e < nEntity; e++){
- for(unsigned int n = 0; n < nNodeInEntity[e]; n++)
- (*permutedRefEntity)[p][e][n] = perm[p][refEntity[e][n]];
+ orderRefEntityForGivenRefSpace(refFaceNodeIdx[f],
+ refSpaceNodeId[s],
+ orderedFaceNodeIdx[s][f]);
+ if(nNodeInFace == 4)
+ correctQuadFaceNodeIdx(orderedFaceNodeIdx[s][f]);
}
}
}
-*/
-const vector<unsigned int>* ReferenceSpace::
-getOrderedEdge(unsigned int permutation,
- unsigned int edge){
-
- /**************************************************************************
- * I want to know how I need to take the values of refEdge[edge] *
- * (i.e the node *index* of a given edge) so that the corresponding *
- * values in permutation (i.e perm[permutation][refEdge[edge]]) are going *
- * from the smallest value to the biggest. *
- **************************************************************************/
-
- // 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,
- *ordered);
-
- // Return ordered
- return ordered;
-}
-const std::vector<unsigned int>* ReferenceSpace::
-getOrderedTriFace(unsigned int permutation,
- unsigned int face){
-
- /***********************************************
- * Same as getOrderedEdge, but I with refFace. *
- ***********************************************/
+void ReferenceSpace::
+orderRefEntityForGivenRefSpace(vector<size_t>& refEntityNodeIdx,
+ vector<size_t>& refSpaceNodeId,
+ vector<size_t>& orderedEntityNodeIdx){
+ // Get Size
+ const size_t size = orderedEntityNodeIdx.size();
- // Alloc
- vector<unsigned int>* ordered = new vector<unsigned int>(3);
+ // Ref Entity Node *ID*
+ vector<pair<size_t, size_t> > refEntityNodeId(size);
- // Permutation
- vector<size_t> perm(nVertex);
- pTree->fillWithPermutation(permutation, perm);
+ for(size_t i = 0; i < size; i++){
+ refEntityNodeId[i].first = i; // Node Id INDEX
+ refEntityNodeId[i].second = refSpaceNodeId[refEntityNodeIdx[i]]; // Node Id
+ }
- // Order refFace
- orderRefEntityForGivenPermutation(refFace[face],
- perm,
- *ordered);
+ // Sort it with repsect to second entry (Node ID)
+ std::sort(refEntityNodeId.begin(), refEntityNodeId.end(), sortPredicate);
- // Return ordered
- return ordered;
+ // Populate orderedEntityNodeIdx (Usign sorted Node ID old Index)
+ for(size_t i = 0; i < size; i++)
+ orderedEntityNodeIdx[i] = refEntityNodeIdx[refEntityNodeId[i].first];
}
-const std::vector<unsigned int>* ReferenceSpace::
-getOrderedQuadFace(unsigned int permutation,
- unsigned int face){
-
- /*********************************************************************
- * Same as getOrderedFace, but refFace has 4 nodes. *
- * *
- * In that case, *
- * I need node at ordered[2] to be *opposite* to node at ordered[0]. *
- * --> If not make a permutation such *
- * that node opposite at ordered[2] *
- * is at ordered[0] *
- *********************************************************************/
-
- // Alloc
- vector<unsigned int>* ordered = new vector<unsigned int>(4);
-
- // Permutation
- vector<size_t> perm(nVertex);
- pTree->fillWithPermutation(permutation, perm);
+void ReferenceSpace::
+correctQuadFaceNodeIdx(vector<size_t>& correctedQuadFaceNodeIdx){
- // Order refFace
- orderRefEntityForGivenPermutation(refFace[face],
- perm,
- *ordered);
+ // Get : //
+ // correctedQuadFaceNodeIdx[2] //
+ // is *opposite* to //
+ // correctedQuadFaceNodeIdx[0] //
- // Get ordered[2] opposite to ordered[0]
- unsigned int opposite = ((*ordered)[0] + 2) % 4;
+ size_t opposite = (correctedQuadFaceNodeIdx[0] + 2) % 4;
- if((*ordered)[2] != opposite){
+ if(correctedQuadFaceNodeIdx[2] != opposite){
// Find opposite node index
- unsigned int tmp;
- unsigned int idx = 1;
- while((*ordered)[idx] != opposite)
- idx++;
-
- // Swap ordered[2] and ordered[idx]
- tmp = (*ordered)[2];
- (*ordered)[2] = opposite;
- (*ordered)[idx] = tmp;
- }
-
- // Return
- return ordered;
-}
-
-void ReferenceSpace::
-orderRefEntityForGivenPermutation(size_t* refEntity,
- vector<size_t>& permutation,
- vector<unsigned int>& orderedEntity){
- // Get Size
- const unsigned int size = orderedEntity.size();
+ size_t tmp;
+ size_t idx = 1;
- // Copy srcSort in sorted
- vector<pair<unsigned int, unsigned int> > sorted(size);
+ while(correctedQuadFaceNodeIdx[idx] != opposite)
+ idx++;
- for(unsigned int i = 0; i < size; i++){
- sorted[i].first = i;
- sorted[i].second = permutation[refEntity[i]];
+ // Swap correctedQuadFaceNodeIdx[2] and correctedQuadFaceNodeIdx[idx]
+ tmp = correctedQuadFaceNodeIdx[2];
+ correctedQuadFaceNodeIdx[2] = opposite;
+ correctedQuadFaceNodeIdx[idx] = tmp;
}
-
- // Sort it with repsect to second entry
- std::sort(sorted.begin(), sorted.end(), sortPredicate);
-
- // Swap with respect to srcSwap
- for(unsigned int i = 0; i < size; i++)
- orderedEntity[i] = refEntity[sorted[i].first];
}
-size_t ReferenceSpace::getReferenceSpace(const MElement& elem) const{
- // Const_Cast //
- MElement& element = const_cast<MElement&>(elem);
-
+size_t ReferenceSpace::getPermutationIdx(const MElement& element) const{
// Get Primary Vertices Global ID //
- const int nVertex = element.getNumPrimaryVertices();
- vector<pair<unsigned int, unsigned int> > vertexGlobalId(nVertex);
+ vector<pair<size_t, size_t> > vertexGlobalId(nVertex);
- for(int i = 0; i < nVertex; i++){
+ for(size_t i = 0; i < nVertex; i++){
vertexGlobalId[i].first = i;
vertexGlobalId[i].second = element.getVertex(i)->getNum();
}
@@ -517,7 +405,7 @@ size_t ReferenceSpace::getReferenceSpace(const MElement& elem) const{
// Reduce Vertex Global ID //
vector<size_t> vertexReducedId(nVertex);
- for(int i = 0; i < nVertex; i++)
+ for(size_t i = 0; i < nVertex; i++)
vertexReducedId[vertexGlobalId[i].first] = i;
// Tree Lookup //
@@ -531,37 +419,114 @@ size_t ReferenceSpace::getReferenceSpace(const MElement& elem) const{
}
}
+void ReferenceSpace::mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]){
+ // Get UVW coordinate //
+ double phys[3] = {xyz(0), xyz(1), xyz(2)};
+ double uvw[3];
+
+ element.xyz2uvw(phys, uvw);
+
+ // Compute coordinate in ABC Space //
+ // ABC node coordinate
+ double** abcMat = new double*[nVertex];
+
+ for(size_t i = 0; i < nVertex; i++)
+ abcMat[i] = new double[3];
+
+ for(size_t i = 0; i < nVertex; i++)
+ element.getNode(i, abcMat[i][0], abcMat[i][1], abcMat[i][2]);
+
+ // UVW (order 1) shape functions
+ double* phiUVW = new double[nVertex];
+ element.getShapeFunctions(uvw[0], uvw[1], uvw[2], phiUVW, 1);
+
+ // Element Permutation Index //
+ const size_t permutationIdx = getPermutationIdx(element);
+
+ // Map From UVW to ABC //
+ abc[0] = 0;
+ for(size_t i = 0; i < nVertex; i++)
+ abc[0] += abcMat[i][0] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+
+ abc[1] = 0;
+ for(size_t i = 0; i < nVertex; i++)
+ abc[1] += abcMat[i][1] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+
+ abc[2] = 0;
+ for(size_t i = 0; i < nVertex; i++)
+ abc[2] += abcMat[i][2] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+
+ // Free //
+ delete[] phiUVW;
+
+ for(size_t i = 0; i < nVertex; i++)
+ delete[] abcMat[i];
+
+ delete[] abcMat;
+}
+
+void ReferenceSpace::getJacobian(const MElement& element,
+ const fullVector<double>& xyz,
+ fullMatrix<double>& jac) const{
+}
+
string ReferenceSpace::toString(void) const{
- const size_t nPerm = pTree->getNPermutation();
+ const size_t nRefSpace = refSpaceNodeId.size();
+ const size_t nEdge = refEdgeNodeIdx.size();
+ const size_t nFace = refFaceNodeIdx.size();
stringstream stream;
- // ReferenceSpaces //
- stream << "Reference Spaces:" << endl
+ // Permutation Tree //
+ stream << "Permutation Tree:" << endl
<< pTree->toString() << endl;
- stream << "Edges Permutations:" << endl;
+ // Reference Space //
+ stream << "Reference Space Node IDs:" << endl;
- for(size_t i = 0; i < nPerm; i++){
- stream << " * RefSpace #" << i + 1 << ":" << endl;
+ for(size_t s = 0; s < nRefSpace; s++){
+ stream << " ** Reference Space #" << s + 1 << ": [";
- for(size_t j = 0; j < nEdge; j++)
- stream << " -- ["
- << edge->at(i)->at(j)->at(0) << ", "
- << edge->at(i)->at(j)->at(1) << "]" << endl;
+ for(size_t v = 0; v < nVertex - 1; v++)
+ stream << refSpaceNodeId[s][v] << ", ";
+
+ stream << refSpaceNodeId[s][nVertex - 1] << "]" << endl;
}
- stream << "Faces Permutations:" << endl;
+ stream << endl << "Ordered Edge Node Index:" << endl;
- for(size_t i = 0; i < nPerm; i++){
- stream << " * RefSpace #" << i + 1 << ":" << endl;
+ for(size_t s = 0; s < nRefSpace; s++){
+ stream << " ** Reference Space #" << s + 1 << ":"
+ << endl;
- for(size_t j = 0; j < nFace; j++)
+ for(size_t e = 0; e < nEdge; e++)
stream << " -- ["
- << face->at(i)->at(j)->at(0) << ", "
- << face->at(i)->at(j)->at(1) << ", "
- << face->at(i)->at(j)->at(2) << "]" << endl;
+ << orderedEdgeNodeIdx[s][e][0] << ", "
+ << orderedEdgeNodeIdx[s][e][1] << "]"
+ << endl;
+ }
+
+ stream << endl << "Ordered Face Node Index:" << endl;
+
+ for(size_t s = 0; s < nRefSpace; s++){
+ stream << " ** Reference Space #" << s + 1 << ":"
+ << endl;
+
+ for(size_t f = 0; f < nFace; f++){
+ const size_t nNodeInFace = orderedFaceNodeIdx[s][f].size();
+ stream << " -- [";
+
+ for(size_t n = 0; n < nNodeInFace - 1; n++)
+ stream << orderedFaceNodeIdx[s][f][n] << ", ";
+
+ stream << orderedFaceNodeIdx[s][f][nNodeInFace - 1] << "]"
+ << endl;
+ }
}
+ stream << endl;
+
return stream.str();
}
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index aa04bf78b9..ffdf867992 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -3,10 +3,10 @@
#include <vector>
#include <list>
-#include <stack>
#include <string>
#include "MElement.h"
+#include "fullMatrix.h"
#include "PermutationTree.h"
/**
@@ -28,70 +28,51 @@ class ReferenceSpace{
} triple;
protected:
+ // Element Definition //
+ size_t nVertex;
+ std::vector<std::vector<size_t> > refEdgeNodeIdx;
+ std::vector<std::vector<size_t> > refFaceNodeIdx;
+
// Permutation Tree //
PermutationTree* pTree;
- //////////////////////////////////////////////////////////////////////////
- // !!! CHANGE VARIABLES NAME !!! //
- // //
- // perm = nodeReducedGlobalId //
- // refEntity = entityNodeIndex //
- // entity = orientedEntityNodeIndex //
- // //
- // !!! I'm Working on *BOTH* indexes and global IDs !!! //
- // //
- // 16 + //
- // |\ ---> NodeIndex = [v0 v1 v2] //
- // | \ ---> GlobalID = [17 42 16] //
- // | \ ---> ReducedGlobalId = [1 2 0 ] ---> Permutation = 3 //
- // 17 +---+ 42 //
- // //
- // * refEdge[e][i] = ith INDEX of edge[e] //
- // For example refEdge[1] = [1 2] because edge[1] //
- // - - //
- // is made of nodes v1 and v2 //
- // - - //
- // //
- // * edge[p][e][i] = ith INDEX of edge[e] in orientation[p] //
- // such that GlobalID[edge[p][e][0]] > GlobalId[edge[p][e][0]] //
- // For example edge[3][1] = [2 1] because edge[1] goes from v2 to v1 //
- // - - - - //
- //////////////////////////////////////////////////////////////////////////
-
- // Number of Vertices //
- size_t nVertex;
+ // Reference Spaces Node Ids //
+ std::vector<std::vector<size_t> > refSpaceNodeId;
- // Edge Permutation //
- size_t nEdge;
- size_t** refEdge;
- //unsigned int*** permutedRefEdge;
- std::vector<const std::vector<const std::vector<unsigned int>*>*>* edge;
+ // ABC space to UVW space shape function index mapping //
+ std::vector<std::vector<size_t> > ABCtoUVWIndex;
- // Face Permutation //
- size_t nFace;
- size_t* nNodeInFace;
- size_t** refFace;
- //unsigned int*** permutedRefFace;
- std::vector<const std::vector<const std::vector<unsigned int>*>*>* face;
+ // UVW space to ABC space shape function index mapping //
+ std::vector<std::vector<size_t> > UVWtoABCIndex;
+
+ // Ordered Edge & Face Node Index //
+ std::vector<std::vector<std::vector<size_t> > > orderedEdgeNodeIdx;
+ std::vector<std::vector<std::vector<size_t> > > orderedFaceNodeIdx;
public:
virtual ~ReferenceSpace(void);
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;
- const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
- getAllFace(void) const;
+ const std::vector<std::vector<size_t> >&
+ getNodeId(void) const;
- virtual std::string toString(void) const;
+ const std::vector<std::vector<std::vector<size_t> > >&
+ getEdgeNodeIndex(void) const;
+
+ const std::vector<std::vector<std::vector<size_t> > >&
+ getFaceNodeIndex(void) const;
+
+ void mapFromXYZtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]);
+
+ void getJacobian(const MElement& element,
+ const fullVector<double>& xyz,
+ fullMatrix<double>& jac) const;
+ virtual std::string toString(void) const;
virtual std::string toLatex(void) const;
protected:
@@ -99,38 +80,32 @@ class ReferenceSpace{
void init(void);
- void getEdge(void);
- void getFace(void);
+ void getOrderedEdge(void);
+ void getOrderedFace(void);
+
+ void findCyclicPermutation
+ (std::list<size_t>& listOfTrueReferenceSpace,
+ std::list<std::vector<size_t> >& listOfRefNodeIndexPermutation,
+ std::list<std::vector<size_t> >& listOfReverseNodeIndexPermutation);
- 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,
- unsigned int* nNodeInEntity,
- unsigned int nEntity);
- */
- const std::vector<unsigned int>* getOrderedEdge(unsigned int permutation,
- unsigned int edge);
-
- const std::vector<unsigned int>* getOrderedTriFace(unsigned int permutation,
- unsigned int face);
-
- const std::vector<unsigned int>* getOrderedQuadFace(unsigned int permutation,
- unsigned int face);
+ static void
+ orderRefEntityForGivenRefSpace(std::vector<size_t>& refEntityNodeIdx,
+ std::vector<size_t>& refSpaceNodeId,
+ std::vector<size_t>& orderedEntityNodeIdx);
static void
- orderRefEntityForGivenPermutation(size_t* refEntity,
- std::vector<size_t>& permutation,
- std::vector<unsigned int>& orderedEntity);
+ correctQuadFaceNodeIdx(std::vector<size_t>& correctedQuadFaceNodeIdx);
+
+ size_t getPermutationIdx(const MElement& element) const;
- static bool sortPredicate(const std::pair<unsigned int, unsigned int>& a,
- const std::pair<unsigned int, unsigned int>& b);
+ static bool sortPredicate(const std::pair<size_t, size_t>& a,
+ const std::pair<size_t, size_t>& b);
};
@@ -198,26 +173,30 @@ class ReferenceSpace{
//////////////////////
inline size_t ReferenceSpace::getNReferenceSpace(void) const{
- return pTree->getNPermutation();
+ return refSpaceNodeId.size();
}
inline
-const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
-ReferenceSpace::getAllEdge(void) const{
- return *edge;
+const std::vector<std::vector<std::vector<size_t> > >&
+ReferenceSpace::getEdgeNodeIndex(void) const{
+ return orderedEdgeNodeIdx;
}
inline
-const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
-ReferenceSpace::getAllFace(void) const{
- return *face;
+const std::vector<std::vector<std::vector<size_t> > >&
+ReferenceSpace::getFaceNodeIndex(void) const{
+ return orderedFaceNodeIdx;
}
+inline
+size_t ReferenceSpace::getReferenceSpace(const MElement& element) const{
+ return pTree->getTagFromPermutation(getPermutationIdx(element));
+}
inline
bool
-ReferenceSpace::sortPredicate(const std::pair<unsigned int, unsigned int>& a,
- const std::pair<unsigned int, unsigned int>& b){
+ReferenceSpace::sortPredicate(const std::pair<size_t, size_t>& a,
+ const std::pair<size_t, size_t>& b){
return a.second < b.second;
}
diff --git a/FunctionSpace/ReferenceSpaceLagrange.cpp b/FunctionSpace/ReferenceSpaceLagrange.cpp
index 0a618ea115..932ec9ce4c 100644
--- a/FunctionSpace/ReferenceSpaceLagrange.cpp
+++ b/FunctionSpace/ReferenceSpaceLagrange.cpp
@@ -21,6 +21,7 @@ ReferenceSpaceLagrange::~ReferenceSpaceLagrange(void){
}
void ReferenceSpaceLagrange::getLagrangeNode(void){
+ /*
const size_t nPerm = pTree->getNPermutation();
// Alloc //
@@ -58,6 +59,7 @@ void ReferenceSpaceLagrange::getLagrangeNode(void){
// Insert in node
(*node)[p] = tmp;
}
+ */
}
void ReferenceSpaceLagrange::edgeSeq(vector<unsigned int>& vec,
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 405ce362d3..92bb3e6649 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
TetEdgeBasis::TetEdgeBasis(unsigned int order){
+ /*
// Reference Space //
refSpace = new TetReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -321,9 +322,11 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
+ */
}
TetEdgeBasis::~TetEdgeBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -336,4 +339,5 @@ TetEdgeBasis::~TetEdgeBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/TetNedelecBasis.cpp b/FunctionSpace/TetNedelecBasis.cpp
index be15c1a6d9..c803758152 100644
--- a/FunctionSpace/TetNedelecBasis.cpp
+++ b/FunctionSpace/TetNedelecBasis.cpp
@@ -4,6 +4,7 @@
using namespace std;
TetNedelecBasis::TetNedelecBasis(void){
+ /*
// Reference Space //
refSpace = new TetReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -66,9 +67,11 @@ TetNedelecBasis::TetNedelecBasis(void){
basis[s][e] = new vector<Polynomial>(tmp2);
}
}
+ */
}
TetNedelecBasis::~TetNedelecBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -81,4 +84,5 @@ TetNedelecBasis::~TetNedelecBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index cd221914e9..4af8626ad8 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
TetNodeBasis::TetNodeBasis(unsigned int order){
+ /*
// Reference Space //
refSpace = new TetReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -154,9 +155,11 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
+ */
}
TetNodeBasis::~TetNodeBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -169,4 +172,5 @@ TetNodeBasis::~TetNodeBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/TetReferenceSpace.cpp b/FunctionSpace/TetReferenceSpace.cpp
index 154b2a9e73..839f956229 100644
--- a/FunctionSpace/TetReferenceSpace.cpp
+++ b/FunctionSpace/TetReferenceSpace.cpp
@@ -5,6 +5,7 @@
using namespace std;
TetReferenceSpace::TetReferenceSpace(void){
+ /*
// Vertex Definition //
nVertex = 4;
@@ -40,9 +41,11 @@ TetReferenceSpace::TetReferenceSpace(void){
// Init All //
init();
+ */
}
TetReferenceSpace::~TetReferenceSpace(void){
+ /*
// Delete Ref Edge //
for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
@@ -55,13 +58,14 @@ TetReferenceSpace::~TetReferenceSpace(void){
delete[] refFace;
delete[] nNodeInFace;
+ */
}
string TetReferenceSpace::toLatex(void) const{
- const size_t nPerm = pTree->getNPermutation();
+ //const size_t nPerm = pTree->getNPermutation();
stringstream stream;
- vector<size_t> perm;
-
+ //vector<size_t> perm;
+ /*
stream << "\\documentclass{article}" << endl << endl
<< "\\usepackage{longtable}" << endl
@@ -101,6 +105,6 @@ string TetReferenceSpace::toLatex(void) const{
stream << "\\end{longtable}" << endl
<< "\\end{document}" << endl;
-
+ */
return stream.str();
}
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 251807a1c5..e39649f44e 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -5,6 +5,7 @@
using namespace std;
TriEdgeBasis::TriEdgeBasis(unsigned int order){
+ /*
// Reference Space //
refSpace = new TriReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -235,9 +236,11 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
delete[] u;
delete[] v;
delete[] subGrad;
+ */
}
TriEdgeBasis::~TriEdgeBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -250,4 +253,5 @@ TriEdgeBasis::~TriEdgeBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/TriLagrangeReferenceSpace.cpp b/FunctionSpace/TriLagrangeReferenceSpace.cpp
index b3c376f654..7ba864c88f 100644
--- a/FunctionSpace/TriLagrangeReferenceSpace.cpp
+++ b/FunctionSpace/TriLagrangeReferenceSpace.cpp
@@ -4,6 +4,7 @@
using namespace std;
TriLagrangeReferenceSpace::TriLagrangeReferenceSpace(unsigned int order){
+ /*
// Vertex Definition //
nVertex = 3;
@@ -41,9 +42,11 @@ TriLagrangeReferenceSpace::TriLagrangeReferenceSpace(unsigned int order){
nNodePerCell;
getLagrangeNode();
+ */
}
TriLagrangeReferenceSpace::~TriLagrangeReferenceSpace(void){
+ /*
// Delete Ref Edge //
for(size_t i = 0; i < nEdge; i++)
delete[] refEdge[i];
@@ -55,4 +58,5 @@ TriLagrangeReferenceSpace::~TriLagrangeReferenceSpace(void){
delete[] refFace[i];
delete[] refFace;
+ */
}
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index 718ede16f2..645aa03908 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -4,6 +4,7 @@
using namespace std;
TriNedelecBasis::TriNedelecBasis(void){
+ /*
// Reference Space //
refSpace = new TriReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
@@ -63,9 +64,11 @@ TriNedelecBasis::TriNedelecBasis(void){
basis[s][e] = new vector<Polynomial>(tmp2);
}
}
+ */
}
TriNedelecBasis::~TriNedelecBasis(void){
+ /*
// ReferenceSpace //
delete refSpace;
@@ -78,4 +81,5 @@ TriNedelecBasis::~TriNedelecBasis(void){
}
delete[] basis;
+ */
}
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 957722af4d..e5f5890c63 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -4,16 +4,16 @@
using namespace std;
-TriNodeBasis::TriNodeBasis(unsigned int order){
+TriNodeBasis::TriNodeBasis(size_t order){
// Reference Space //
refSpace = new TriReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
- const vector<const vector<const vector<unsigned int>*>*>&
- edgeV = refSpace->getAllEdge();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
- const vector<const vector<const vector<unsigned int>*>*>&
- faceV = refSpace->getAllFace();
+ const vector<vector<vector<size_t> > >&
+ faceIdx = refSpace->getFaceNodeIndex();
// Set BasisTwo Type //
this->order = order;
@@ -51,28 +51,28 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
// Basis //
basis = new Polynomial**[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new Polynomial*[nFunction];
// Vertex Based //
- for(unsigned int s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nRefSpace; s++){
basis[s][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]);
}
// Edge Based //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nVertex;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nVertex;
- for(unsigned int e = 0; e < 3; e++){
- for(unsigned int l = 1; l < order; l++){
+ for(size_t e = 0; e < 3; e++){
+ for(size_t l = 1; l < order; l++){
basis[s][i] =
- new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
- lagrange[(*(*edgeV[s])[e])[0]]
+ new Polynomial(intLegendre[l].compose(lagrange[edgeIdx[s][e][1]] -
+ lagrange[edgeIdx[s][e][0]]
,
- lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]]));
+ lagrange[edgeIdx[s][e][0]] +
+ lagrange[edgeIdx[s][e][1]]));
i++;
}
}
@@ -80,30 +80,30 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
// Face Based //
- // NB: We use (*(*faceV[s])[f])[]
+ // NB: We use (*(*faceIdx[s])[f])[]
// where f = 0, because triangles
// have only ONE face: the face '0'
const int orderMinusTwo = order - 2;
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nVertex + nEdge;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nVertex + nEdge;
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]] -
- lagrange[(*(*faceV[s])[0])[0]]
+ new Polynomial(intLegendre[l1].compose(lagrange[faceIdx[s][0][1]] -
+ lagrange[faceIdx[s][0][0]]
,
- lagrange[(*(*faceV[s])[0])[0]] +
- lagrange[(*(*faceV[s])[0])[1]])
+ lagrange[faceIdx[s][0][0]] +
+ lagrange[faceIdx[s][0][1]])
*
- legendre[l2].compose((lagrange[(*(*faceV[s])[0])[2]] * 2)
+ legendre[l2].compose((lagrange[faceIdx[s][0][2]] * 2)
-
Polynomial(1, 0, 0, 0))
*
- lagrange[(*(*faceV[s])[0])[2]]);
+ lagrange[faceIdx[s][0][2]]);
i++;
}
}
@@ -119,8 +119,8 @@ TriNodeBasis::~TriNodeBasis(void){
delete refSpace;
// Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
+ for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
diff --git a/FunctionSpace/TriNodeBasis.h b/FunctionSpace/TriNodeBasis.h
index 80d96f3e9d..2d1b98d304 100644
--- a/FunctionSpace/TriNodeBasis.h
+++ b/FunctionSpace/TriNodeBasis.h
@@ -20,7 +20,7 @@ class TriNodeBasis: public BasisHierarchical0From{
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Triangles of the given order
- TriNodeBasis(unsigned int order);
+ TriNodeBasis(size_t order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TriReferenceSpace.cpp b/FunctionSpace/TriReferenceSpace.cpp
index b976bc43ed..2b7acabe43 100644
--- a/FunctionSpace/TriReferenceSpace.cpp
+++ b/FunctionSpace/TriReferenceSpace.cpp
@@ -9,54 +9,33 @@ TriReferenceSpace::TriReferenceSpace(void){
nVertex = 3;
// Edge Definition //
- nEdge = 3;
- refEdge = new size_t*[nEdge];
+ const size_t nEdge = 3;
+ refEdgeNodeIdx.resize(nEdge);
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);
+ refEdgeNodeIdx[i].resize(2); // Two Nodes per Edge
+ refEdgeNodeIdx[i][0] = MTriangle::edges_tri(i, 0);
+ refEdgeNodeIdx[i][1] = MTriangle::edges_tri(i, 1);
}
// Face Definition //
- // Number of face
- nFace = 1;
+ refFaceNodeIdx.resize(1); // One Face per Triangle
+ refFaceNodeIdx[0].resize(3); // Three Nodes per Face
- // Number of node per face
- nNodeInFace = new size_t[nFace];
- nNodeInFace[0] = 3;
-
- // Reference Face
- refFace = new size_t*[nFace];
- refFace[0] = new size_t[3];
-
- refFace[0][0] = 0;
- refFace[0][1] = 1;
- refFace[0][2] = 2;
+ refFaceNodeIdx[0][0] = 0;
+ refFaceNodeIdx[0][1] = 1;
+ refFaceNodeIdx[0][2] = 2;
// Init All //
init();
}
TriReferenceSpace::~TriReferenceSpace(void){
- // Delete Ref Edge //
- for(size_t i = 0; i < nEdge; i++)
- delete[] refEdge[i];
-
- delete[] refEdge;
-
- // Delete Ref Face //
- for(size_t i = 0; i < nFace; i++)
- delete[] refFace[i];
-
- delete[] refFace;
- delete[] nNodeInFace;
}
string TriReferenceSpace::toLatex(void) const{
- const size_t nPerm = pTree->getNPermutation();
+ const size_t nRefSpace = refSpaceNodeId.size();
stringstream stream;
- vector<size_t> perm(nVertex);
stream << "\\documentclass{article}" << endl << endl
@@ -70,24 +49,25 @@ string TriReferenceSpace::toLatex(void) const{
<< "\\begin{longtable}{ccc}" << endl << endl;
- for(size_t p = 0; p < nPerm; p++){
- pTree->fillWithPermutation(p, perm);
-
+ for(size_t s = 0; s < nRefSpace; s++){
stream << "\\begin{tikzpicture}" << 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] (n0) at(0, 0) {$" << refSpaceNodeId[s][0] << "$};"
+ << endl
+ << "\\node[vertex] (n1) at(3, 0) {$" << refSpaceNodeId[s][1] << "$};"
+ << endl
+ << "\\node[vertex] (n2) at(0, 3) {$" << refSpaceNodeId[s][2] << "$};"
+ << endl
<< endl;
- for(size_t i = 0; i < 3; i++)
+ for(size_t e = 0; e < 3; e++)
stream << "\\path[line]"
- << " (n" << (*(*(*edge)[p])[i])[0] << ")"
+ << " (n" << orderedEdgeNodeIdx[s][e][0] << ")"
<< " -- "
- << " (n" << (*(*(*edge)[p])[i])[1] << ");"
+ << " (n" << orderedEdgeNodeIdx[s][e][1] << ");"
<< endl;
- if((p + 1) % 3)
+ if((s + 1) % 3)
stream << "\\end{tikzpicture} & " << endl << endl;
else
--
GitLab