From 9e26f41281805837d4df490e87c5d2b5df33c958 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Thu, 6 Jun 2013 15:49:28 +0000
Subject: [PATCH] ReferenceSpace reworking -- seems ok
---
FunctionSpace/FunctionSpace.cpp | 6 +-
FunctionSpace/FunctionSpaceScalar.cpp | 8 +-
FunctionSpace/FunctionSpaceVector.cpp | 8 +-
FunctionSpace/LineReferenceSpace.cpp | 6 +-
FunctionSpace/QuadReferenceSpace.cpp | 46 +--
FunctionSpace/ReferenceSpace.cpp | 358 ++++++++++++--------
FunctionSpace/ReferenceSpace.h | 94 +++--
FunctionSpace/TetReferenceSpace.cpp | 16 +-
FunctionSpace/TriLagrangeBasis.cpp | 6 +-
FunctionSpace/TriLagrangeReferenceSpace.cpp | 4 +-
FunctionSpace/TriReferenceSpace.cpp | 14 +-
11 files changed, 343 insertions(+), 223 deletions(-)
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 0fe11bea20..ae99c457bb 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -99,8 +99,8 @@ void FunctionSpace::buildDof(void){
dof = new set<const Dof*, DofComparator>;
group = new vector<GroupOfDof*>(nElement);
eToGod = new map<const MElement*,
- const GroupOfDof*,
- ElementComparator>;
+ const GroupOfDof*,
+ ElementComparator>;
// Create Dofs //
for(unsigned int i = 0; i < nElement; i++){
@@ -118,7 +118,7 @@ void FunctionSpace::buildDof(void){
// Map GOD
eToGod->insert(pair<const MElement*, const GroupOfDof*>
- (element[i], god));
+ (element[i], god));
}
}
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index f65cd8a383..42eb4218cb 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -12,8 +12,8 @@ FunctionSpaceScalar::~FunctionSpaceScalar(void){
double FunctionSpaceScalar::
interpolate(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& xyz) const{
+ const std::vector<double>& coef,
+ const fullVector<double>& xyz) const{
// Const Cast For MElement //
MElement& eelement =
@@ -44,8 +44,8 @@ interpolate(const MElement& element,
double FunctionSpaceScalar::
interpolateInRefSpace(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& uvw) const{
+ const std::vector<double>& coef,
+ const fullVector<double>& uvw) const{
// Get Basis Functions //
fullMatrix<double>* fun =
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index 8ca871ef95..ce4b89207b 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -13,8 +13,8 @@ FunctionSpaceVector::~FunctionSpaceVector(void){
fullVector<double> FunctionSpaceVector::
interpolate(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& xyz) const{
+ const std::vector<double>& coef,
+ const fullVector<double>& xyz) const{
// Const Cast For MElement //
MElement& eelement =
@@ -59,8 +59,8 @@ interpolate(const MElement& element,
fullVector<double> FunctionSpaceVector::
interpolateInRefSpace(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& uvw) const{
+ const std::vector<double>& coef,
+ const fullVector<double>& uvw) const{
// Const Cast For MElement //
MElement& eelement =
diff --git a/FunctionSpace/LineReferenceSpace.cpp b/FunctionSpace/LineReferenceSpace.cpp
index b71ed0cda1..7ba0eb2d4a 100644
--- a/FunctionSpace/LineReferenceSpace.cpp
+++ b/FunctionSpace/LineReferenceSpace.cpp
@@ -19,10 +19,8 @@ LineReferenceSpace::LineReferenceSpace(void){
nFace = 0;
refFace = NULL;
- // Init All (Use Tri Face -- //
- // unused -- //
- // just for interface) //
- init(0);
+ // Init All //
+ init();
}
LineReferenceSpace::~LineReferenceSpace(void){
diff --git a/FunctionSpace/QuadReferenceSpace.cpp b/FunctionSpace/QuadReferenceSpace.cpp
index 5bc382730c..5a55d5af8b 100644
--- a/FunctionSpace/QuadReferenceSpace.cpp
+++ b/FunctionSpace/QuadReferenceSpace.cpp
@@ -19,7 +19,14 @@ QuadReferenceSpace::QuadReferenceSpace(void){
}
// Face Definition //
- nFace = 1;
+ // Number of face
+ nFace = 1;
+
+ // Number of node per face
+ nNodeInFace = new unsigned int[nFace];
+ nNodeInFace[0] = 4;
+
+ // Reference Face
refFace = new unsigned int*[nFace];
refFace[0] = new unsigned int[4];
@@ -28,8 +35,8 @@ QuadReferenceSpace::QuadReferenceSpace(void){
refFace[0][2] = 2;
refFace[0][3] = 3;
- // Init All (Quad Face) //
- init(1);
+ // Init All //
+ init();
}
QuadReferenceSpace::~QuadReferenceSpace(void){
@@ -44,6 +51,7 @@ QuadReferenceSpace::~QuadReferenceSpace(void){
delete[] refFace[i];
delete[] refFace;
+ delete[] nNodeInFace;
}
string QuadReferenceSpace::toLatex(void) const{
@@ -51,31 +59,31 @@ string QuadReferenceSpace::toLatex(void) const{
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(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(3, 3) {$" << perm[p][2] << "$};" << endl
- << "\\node[vertex] (n3) at(0, 3) {$" << perm[p][3] << "$};" << 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
<< endl;
for(unsigned int i = 0; i < 4; 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;
@@ -85,7 +93,7 @@ string QuadReferenceSpace::toLatex(void) const{
}
stream << "\\end{longtable}" << endl
- << "\\end{document}" << endl;
+ << "\\end{document}" << endl;
return stream.str();
}
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index bc7c8c21f5..419ec28a79 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -15,33 +15,25 @@ ReferenceSpace::ReferenceSpace(void){
perm = NULL;
lPerm = NULL;
- nUnconnected = 0;
- unconnected = NULL;
- toBeUnconnected = NULL;
- reduceBy = 0;
-
pTreeRoot.depth = 0;
pTreeRoot.last = NULL;
pTreeRoot.number = 0;
pTreeRoot.possible = NULL;
pTreeRoot.next = NULL;
- nEdge = 0;
- refEdge = NULL;
- edge = NULL;
+ nEdge = 0;
+ refEdge = NULL;
+ permutedRefEdge = NULL;
+ edge = NULL;
- nFace = 0;
- refFace = NULL;
- face = NULL;
+ nFace = 0;
+ nNodeInFace = NULL;
+ refFace = NULL;
+ permutedRefFace = NULL;
+ face = NULL;
// Defining Ref Edge and Face in //
// Dervived Class //
- // And CALL init(faceType) //
- // //
- // faceType is: //
- // * 0 if triangular faces //
- // * 1 if quadragular faces //
- // * 3 if NOT IMPLEMENTED //
}
ReferenceSpace::~ReferenceSpace(void){
@@ -51,32 +43,34 @@ ReferenceSpace::~ReferenceSpace(void){
// Delete Permutated Edge //
for(unsigned int p = 0; p < nPerm; p++){
- for(unsigned int i = 0; i < nEdge; i++)
- delete (*(*edge)[p])[i];
+ for(unsigned int i = 0; i < nEdge; i++){
+ delete[] permutedRefEdge[p][i];
+ delete (*(*edge)[p])[i];
+ }
- delete (*edge)[p];
+ delete[] permutedRefEdge[p];
+ delete (*edge)[p];
}
- delete edge;
+ delete[] permutedRefEdge;
+ delete edge;
// Delete Permutated Face //
for(unsigned int p = 0; p < nPerm; p++){
- for(unsigned int i = 0; i < nFace; i++)
+ for(unsigned int 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(unsigned int faceType){
- if(faceType > 3)
- throw Exception("ReferenceSpace: unknown face type");
-
- if(faceType == 3)
- throw Exception("ReferenceSpace: prism not implemented");
-
+void ReferenceSpace::init(void){
// Init Root //
nPerm = 0;
nextLeafId = 0;
@@ -94,7 +88,7 @@ void ReferenceSpace::init(unsigned int faceType){
lPerm = new list<unsigned int*>;
populate(&pTreeRoot);
- // Get Permutations //
+ // Get Permutation //
perm = new unsigned int*[nPerm];
for(unsigned int i = 0; i < nPerm; i++){
// Take Permutation for queue
@@ -107,13 +101,7 @@ void ReferenceSpace::init(unsigned int faceType){
// Get Edges & Faces //
getEdge();
-
- switch(faceType){
- case 0: getTriFace(); break;
- case 1: getQuaFace(); break;
-
- default: throw Exception("ReferenceSpace: impossible error");
- }
+ getFace();
}
void ReferenceSpace::populate(node* pTreeRoot){
@@ -121,10 +109,9 @@ void ReferenceSpace::populate(node* pTreeRoot){
const unsigned int number = pTreeRoot->number;
const unsigned int nextNumber = number - 1;
const unsigned int depth = pTreeRoot->depth;
- const unsigned int nextDepth = pTreeRoot->depth + 1;
+ const unsigned int nextDepth = depth + 1;
// Temp Data //
- unsigned int nextLast;
unsigned int offset;
// If Leaf : a new permutation is found //
@@ -149,13 +136,13 @@ void ReferenceSpace::populate(node* pTreeRoot){
// Init each child node
for(unsigned int i = 0; i < number; i++){
- nextLast = pTreeRoot->possible[i];
- offset = 0;
+ // Reset offset
+ offset = 0;
// 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] = nextLast;
+ pTreeRoot->next[i].last[depth] = pTreeRoot->possible[i];
for(unsigned int j = 0; j < depth; j++)
pTreeRoot->next[i].last[j] = pTreeRoot->last[j];
@@ -165,7 +152,7 @@ void ReferenceSpace::populate(node* pTreeRoot){
pTreeRoot->next[i].possible = new unsigned int[nextNumber];
for(unsigned int j = 0; j < nextNumber; j++){
- if(pTreeRoot->possible[j] == nextLast)
+ if(pTreeRoot->possible[j] == pTreeRoot->possible[i])
offset = 1;
pTreeRoot->next[i].possible[j] = pTreeRoot->possible[j + offset];
@@ -191,165 +178,247 @@ void ReferenceSpace::destroy(node* node){
}
void ReferenceSpace::getEdge(void){
+ // 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;
+
+ // Populate Edge
for(unsigned int p = 0; p < nPerm; p++){
tmp = new vector<const vector<unsigned int>*>(nEdge);
- for(unsigned int i = 0; i < nEdge; i++)
- (*tmp)[i] = inOrder(p,
- refEdge[i][0],
- refEdge[i][1]);
+ for(unsigned int e = 0; e < nEdge; e++)
+ (*tmp)[e] = getOrderedEdge(p, e);
+
(*edge)[p] = tmp;
}
}
-void ReferenceSpace::getTriFace(void){
+void ReferenceSpace::getFace(void){
+ // 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(unsigned int p = 0; p < nPerm; p++){
tmp = new vector<const vector<unsigned int>*>(nFace);
- for(unsigned int i = 0; i < nFace; i++)
- (*tmp)[i] = inOrder(p,
- refFace[i][0],
- refFace[i][1],
- refFace[i][2]);
+ for(unsigned int f = 0; f < nFace; f++){
+ // Dending on the number of node per face
+ // The ordering strategy is different
+
+ switch(nNodeInFace[f]){
+ case 3: (*tmp)[f] = getOrderedTriFace(p, f); break;
+ case 4: (*tmp)[f] = getOrderedQuadFace(p, f); break;
+
+ default:
+ throw Exception("I can handle face with %d nodes",
+ nNodeInFace[f]);
+ }
+ }
+
(*face)[p] = tmp;
}
}
-void ReferenceSpace::getQuaFace(void){
- vector<const vector<unsigned int>*>* tmp;
- face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
+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++){
- tmp = new vector<const vector<unsigned int>*>(nFace);
+ (*permutedRefEntity)[p] = new unsigned int*[nEntity];
- for(unsigned int i = 0; i < nFace; i++)
- (*tmp)[i] = inOrder(p,
- refFace[i][0],
- refFace[i][1],
- refFace[i][2],
- refFace[i][3]);
- (*face)[p] = tmp;
+ for(unsigned int e = 0; e < nEntity; e++)
+ (*permutedRefEntity)[p][e] = new unsigned int[nNodeInEntity[e]];
}
-}
-const vector<unsigned int>* ReferenceSpace::
-inOrder(unsigned int permutation,
- unsigned int a,
- unsigned int b){
+ // Alloc Idx
+ unsigned int** idx = new unsigned int*[nEntity];
- unsigned int v;
- unsigned int k = 0;
- vector<unsigned int>* inorder =
- new vector<unsigned int>(2);
+ for(unsigned int e = 0; e < nEntity; e++)
+ idx[e] = new unsigned int[nNodeInEntity[e]];
- for(unsigned int i = 0; i < nVertex; i++){
- v = perm[permutation][i];
+ // Get Reference Index
+ // Use refEntity and perm[0] as reference element
+ for(unsigned int e = 0; e < nEntity; e++)
+ getIndex(nNodeInEntity[e], refEntity[e], nVertex, perm[0], &idx[e]);
- if(v == a || v == b){
- (*inorder)[k] = v;
- k++;
+ // 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][idx[e][n]];
}
}
- return inorder;
-}
+ // Delete Temp
+ for(unsigned int e = 0; e < nEntity; e++)
+ delete[] idx[e];
-const std::vector<unsigned int>* ReferenceSpace::
-inOrder(unsigned int permutation,
- unsigned int a,
- unsigned int b,
- unsigned int c){
-
- unsigned int v;
- unsigned int k = 0;
- vector<unsigned int>* inorder =
- new vector<unsigned int>(3);
-
- for(unsigned int i = 0; i < nVertex; i++){
- v = perm[permutation][i];
-
- if(v == a || v == b || v == c){
- (*inorder)[k] = v;
- k++;
- }
- }
+ delete[] idx;
+}
- return inorder;
+const vector<unsigned int>* ReferenceSpace::
+getOrderedEdge(unsigned int permutation,
+ unsigned int edge){
+
+ /*************************************************************************
+ * Let say we have we have permutedRefEdge[4][2] = [0 1]. *
+ * I want to get back to same node ID than refEdge[2] = [2 0]. *
+ * That is, I want to return [0 2]. *
+ * *
+ * I need to sort permutedRefEdge and keep the index swaping. *
+ * I can use this index swaping to get the right permutation of refEdge. *
+ *************************************************************************/
+
+ // Alloc
+ vector<unsigned int>* ordered = new vector<unsigned int>(2);
+
+ // Sort & Swap
+ sortAndSwap(permutedRefEdge[permutation][edge],
+ refEdge[edge],
+ *ordered);
+
+ // Return ordered
+ return ordered;
}
const std::vector<unsigned int>* ReferenceSpace::
-inOrder(unsigned int permutation,
- unsigned int a,
- unsigned int b,
- unsigned int c,
- unsigned int d){
-
- unsigned int opposite;
-
- unsigned int v;
- unsigned int k = 0;
- vector<unsigned int>* inorder =
- new vector<unsigned int>(4);
-
- // Take nodes in order //
- for(unsigned int i = 0; i < nVertex; i++){
- v = perm[permutation][i];
-
- if(v == a || v == b || v == c || v == d){
- (*inorder)[k] = v;
- k++;
- }
- }
+getOrderedTriFace(unsigned int permutation,
+ unsigned int face){
+
+ /******************************************************
+ * Same as getOrderedEdge, but I need to use 3 nodes. *
+ ******************************************************/
+
+ // Alloc
+ vector<unsigned int>* ordered = new vector<unsigned int>(3);
- // We need node at inoder[2] to be //
- // opposite to node at inorder[0] //
- // --> if not make a permutation such //
- // that node opposite at inorder[2] //
- // is at inorder[0] //
+ // Sort & Swap
+ sortAndSwap(permutedRefFace[permutation][face],
+ refFace[face],
+ *ordered);
- opposite = ((*inorder)[0] + 2) % 4;
- if((*inorder)[2] != opposite){
+ // Return ordered
+ return ordered;
+}
+
+const std::vector<unsigned int>* ReferenceSpace::
+getOrderedQuadFace(unsigned int permutation,
+ unsigned int face){
+
+ /*********************************************************************
+ * Same as getOrderedEdge, but I need to use 4 nodes. *
+ * *
+ * Moreover, *
+ * 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);
+
+ // Sort & Swap
+ sortAndSwap(permutedRefFace[permutation][face],
+ refFace[face],
+ *ordered);
+
+ // Get ordered[2] opposite to ordered[0]
+ unsigned int opposite = ((*ordered)[0] + 2) % 4;
+
+ if((*ordered)[2] != opposite){
// Find opposite node index
unsigned int tmp;
unsigned int idx = 1;
- while((*inorder)[idx] != opposite)
+ while((*ordered)[idx] != opposite)
idx++;
- // Swap inorder[2] and inorder[idx]
- tmp = (*inorder)[2];
- (*inorder)[2] = opposite;
- (*inorder)[idx] = tmp;
+ // Swap ordered[2] and ordered[idx]
+ tmp = (*ordered)[2];
+ (*ordered)[2] = opposite;
+ (*ordered)[idx] = tmp;
+ }
+
+ // Return
+ return ordered;
+}
+
+void ReferenceSpace::sortAndSwap(unsigned int* srcSort,
+ unsigned int* srcSwap,
+ vector<unsigned int>& dest){
+ // Get Size
+ const unsigned int size = dest.size();
+
+ // Copy srcSort in sorted
+ vector<pair<unsigned int, unsigned int> > sorted(size);
+
+ for(unsigned int i = 0; i < size; i++){
+ sorted[i].first = i;
+ sorted[i].second = srcSort[i];
}
- return inorder;
+ // 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++)
+ dest[i] = srcSwap[sorted[i].first];
}
unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
// Const_Cast //
MElement& element = const_cast<MElement&>(elem);
- // Get Primary Vertices //
+ // Get Primary Vertices Global ID //
const int nVertex = element.getNumPrimaryVertices();
- vector<pair<unsigned int, MVertex*> > vertex(nVertex);
+ vector<pair<unsigned int, unsigned int> > vertexGlobalId(nVertex);
for(int i = 0; i < nVertex; i++){
- vertex[i].first = i;
- vertex[i].second = element.getVertex(i);
+ vertexGlobalId[i].first = i;
+ vertexGlobalId[i].second = element.getVertex(i)->getNum();
}
// Sort Them with repsect to Vertex Global ID //
// (vertex[i].second->getNum) //
- std::sort(vertex.begin(), vertex.end(), sortPredicate);
+ std::sort(vertexGlobalId.begin(), vertexGlobalId.end(), sortPredicate);
+
+ // Reduce Vertex Global ID //
+ vector<unsigned int> vertexReducedId(nVertex);
+
+ for(int i = 0; i < nVertex; i++)
+ vertexReducedId[vertexGlobalId[i].first] = i;
// Tree Lookup //
try{
- return treeLookup(&pTreeRoot, vertex);
+ return treeLookup(&pTreeRoot, vertexReducedId);
}
catch(...){
@@ -359,8 +428,7 @@ unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
}
unsigned int ReferenceSpace::treeLookup(const node* root,
- vector<pair<unsigned int, MVertex*> >&
- sortedArray){
+ vector<unsigned int>& vertexReducedId){
// Temp Data //
unsigned int choice;
unsigned int i;
@@ -368,7 +436,7 @@ unsigned int ReferenceSpace::treeLookup(const node* root,
// If Root is *not* a Leaf: Lookup //
if(root->number){
// Get This Choice
- choice = sortedArray[root->depth].first;
+ choice = vertexReducedId[root->depth];
// Look for next node corresponding to this Choice
i = 0;
@@ -380,7 +448,7 @@ unsigned int ReferenceSpace::treeLookup(const node* root,
throw Exception();
// Go to next Node
- return treeLookup(&root->next[i], sortedArray);
+ return treeLookup(&root->next[i], vertexReducedId);
}
// Else: Return Leaf ID //
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index f376d9f2ee..d051686ff6 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -23,8 +23,8 @@ class ReferenceSpace{
struct node_s{
unsigned int depth; // Depth
unsigned int* last; // Last Choises
- unsigned int number; // Number of Next Choise
- unsigned int* possible; // Possible Next Choise
+ 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
@@ -43,21 +43,19 @@ class ReferenceSpace{
unsigned int** perm;
std::list<unsigned int*>* lPerm;
- unsigned int nUnconnected;
- std::pair<unsigned int, unsigned int>* unconnected;
- std::stack<node*>* toBeUnconnected;
- unsigned int reduceBy;
-
node pTreeRoot;
// Edge Permutation //
unsigned int nEdge;
unsigned int** 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;
+ unsigned int*** permutedRefFace;
std::vector<const std::vector<const std::vector<unsigned int>*>*>* face;
public:
@@ -80,39 +78,46 @@ class ReferenceSpace{
protected:
ReferenceSpace(void);
- void init(unsigned int faceType);
+ void init(void);
void populate(node* pTreeRoot);
void destroy(node* node);
- void unconnectWalk(node* pTreeRoot); // Find wrong permutations
- void markAsUnconnect(node* pTreeRoot); // Mark leafs, with pTreeRoot as root, to be 'unconnected'
- void unconnect(void); // Unconnects leafs marked before
-
void getEdge(void);
- void getTriFace(void);
- void getQuaFace(void);
+ void getFace(void);
+
+ 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>* inOrder(unsigned int permutation,
- unsigned int a,
- unsigned int b);
+ const std::vector<unsigned int>* getOrderedTriFace(unsigned int permutation,
+ unsigned int face);
- const std::vector<unsigned int>* inOrder(unsigned int permutation,
- unsigned int a,
- unsigned int b,
- unsigned int c);
+ const std::vector<unsigned int>* getOrderedQuadFace(unsigned int permutation,
+ unsigned int face);
- const std::vector<unsigned int>* inOrder(unsigned int permutation,
- unsigned int a,
- unsigned int b,
- unsigned int c,
- unsigned int d);
+ static void sortAndSwap(unsigned int* srcSort,
+ unsigned int* srcSwap,
+ std::vector<unsigned int>& dest);
- static bool sortPredicate(const std::pair<unsigned int, MVertex*>& a,
- const std::pair<unsigned int, MVertex*>& b);
+ static unsigned int whereIsIncluded(unsigned int elem,
+ unsigned int* vec,
+ unsigned int size);
+
+ static void getIndex(unsigned int sizeRef,
+ unsigned int* ref,
+ unsigned int sizeVec,
+ unsigned int* vec,
+ unsigned int** idx);
+
+ 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* root,
- std::vector<std::pair<unsigned int, MVertex*> >&
- sortedArray);
+ std::vector<unsigned int>& vertexReducedId);
std::string toString(const node* node) const;
};
@@ -199,11 +204,34 @@ ReferenceSpace::getAllFace(void) const{
return *face;
}
+inline
+void ReferenceSpace::getIndex(unsigned int sizeRef,
+ unsigned int* ref,
+ unsigned int sizeVec,
+ unsigned int* vec,
+ unsigned int** idx){
+
+ for(unsigned int i = 0; i < sizeRef; i++)
+ (*idx)[i] = whereIsIncluded(ref[i], vec, sizeVec);
+}
+
+inline
+unsigned int ReferenceSpace::whereIsIncluded(unsigned int elem,
+ unsigned int* vec,
+ unsigned int size){
+
+ for(unsigned int i = 0; i < size; i++)
+ if(vec[i] == elem)
+ return i;
+
+ return -1;
+}
+
inline
bool
-ReferenceSpace::sortPredicate(const std::pair<unsigned int, MVertex*>& a,
- const std::pair<unsigned int, MVertex*>& b){
- return a.second->getNum() < b.second->getNum();
+ReferenceSpace::sortPredicate(const std::pair<unsigned int, unsigned int>& a,
+ const std::pair<unsigned int, unsigned int>& b){
+ return a.second < b.second;
}
#endif
diff --git a/FunctionSpace/TetReferenceSpace.cpp b/FunctionSpace/TetReferenceSpace.cpp
index 4fcc1f1db7..0c214a408c 100644
--- a/FunctionSpace/TetReferenceSpace.cpp
+++ b/FunctionSpace/TetReferenceSpace.cpp
@@ -19,7 +19,16 @@ TetReferenceSpace::TetReferenceSpace(void){
}
// Face Definition //
- nFace = 4;
+ // Number of face
+ nFace = 4;
+
+ // Number of node per face
+ nNodeInFace = new unsigned int[nFace];
+
+ for(unsigned int f = 0; f < nFace; f++)
+ nNodeInFace[f] = 3;
+
+ // Reference Face
refFace = new unsigned int*[nFace];
for(unsigned int i = 0; i < nFace; i++){
@@ -29,8 +38,8 @@ TetReferenceSpace::TetReferenceSpace(void){
refFace[i][2] = MTetrahedron::faces_tetra(i, 2);
}
- // Init All (Tri Face) //
- init(0);
+ // Init All //
+ init();
}
TetReferenceSpace::~TetReferenceSpace(void){
@@ -45,6 +54,7 @@ TetReferenceSpace::~TetReferenceSpace(void){
delete[] refFace[i];
delete[] refFace;
+ delete[] nNodeInFace;
}
string TetReferenceSpace::toLatex(void) const{
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index ac249adab1..e99e07974c 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -30,14 +30,14 @@ TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
(gmshGeneratePointsTriangle(order, false));
// ReferenceSpace //
- refSpace = new TriLagrangeReferenceSpace(order);
- std::cout << refSpace->toString() << std::endl;
+ // refSpace = new TriLagrangeReferenceSpace(order);
+ // std::cout << refSpace->toString() << std::endl;
}
TriLagrangeBasis::~TriLagrangeBasis(void){
delete lBasis;
delete lPoint;
- delete refSpace;
+ // delete refSpace;
}
unsigned int TriLagrangeBasis::getTag(unsigned int order){
diff --git a/FunctionSpace/TriLagrangeReferenceSpace.cpp b/FunctionSpace/TriLagrangeReferenceSpace.cpp
index ede923049f..d5e8772c69 100644
--- a/FunctionSpace/TriLagrangeReferenceSpace.cpp
+++ b/FunctionSpace/TriLagrangeReferenceSpace.cpp
@@ -26,8 +26,8 @@ TriLagrangeReferenceSpace::TriLagrangeReferenceSpace(unsigned int order){
refFace[0][1] = 1;
refFace[0][2] = 2;
- // Init ReferenceSpace (Tri Face) //
- init(0);
+ // Init ReferenceSpace //
+ init();
// Get Lagrange Node //
nNodePerEdge = 3 * (order - 1) / nEdge;
diff --git a/FunctionSpace/TriReferenceSpace.cpp b/FunctionSpace/TriReferenceSpace.cpp
index c56567d901..681f94f1e7 100644
--- a/FunctionSpace/TriReferenceSpace.cpp
+++ b/FunctionSpace/TriReferenceSpace.cpp
@@ -19,7 +19,14 @@ TriReferenceSpace::TriReferenceSpace(void){
}
// Face Definition //
- nFace = 1;
+ // Number of face
+ nFace = 1;
+
+ // Number of node per face
+ nNodeInFace = new unsigned int[nFace];
+ nNodeInFace[0] = 3;
+
+ // Reference Face
refFace = new unsigned int*[nFace];
refFace[0] = new unsigned int[3];
@@ -27,8 +34,8 @@ TriReferenceSpace::TriReferenceSpace(void){
refFace[0][1] = 1;
refFace[0][2] = 2;
- // Init All (Tri Face) //
- init(0);
+ // Init All //
+ init();
}
TriReferenceSpace::~TriReferenceSpace(void){
@@ -43,6 +50,7 @@ TriReferenceSpace::~TriReferenceSpace(void){
delete[] refFace[i];
delete[] refFace;
+ delete[] nNodeInFace;
}
string TriReferenceSpace::toLatex(void) const{
--
GitLab