Skip to content
Snippets Groups Projects
Commit 9e26f412 authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

ReferenceSpace reworking -- seems ok

parent cd1427a6
No related branches found
No related tags found
No related merge requests found
......@@ -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){
......
......@@ -19,7 +19,14 @@ QuadReferenceSpace::QuadReferenceSpace(void){
}
// Face Definition //
// 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{
......
......@@ -15,11 +15,6 @@ ReferenceSpace::ReferenceSpace(void){
perm = NULL;
lPerm = NULL;
nUnconnected = 0;
unconnected = NULL;
toBeUnconnected = NULL;
reduceBy = 0;
pTreeRoot.depth = 0;
pTreeRoot.last = NULL;
pTreeRoot.number = 0;
......@@ -28,20 +23,17 @@ ReferenceSpace::ReferenceSpace(void){
nEdge = 0;
refEdge = NULL;
permutedRefEdge = NULL;
edge = 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++)
for(unsigned int 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(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];
// 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];
delete[] idx;
}
const std::vector<unsigned int>* ReferenceSpace::
inOrder(unsigned int permutation,
unsigned int a,
unsigned int b,
unsigned int c){
const vector<unsigned int>* ReferenceSpace::
getOrderedEdge(unsigned int permutation,
unsigned int edge){
unsigned int v;
unsigned int k = 0;
vector<unsigned int>* inorder =
new vector<unsigned int>(3);
/*************************************************************************
* 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. *
*************************************************************************/
for(unsigned int i = 0; i < nVertex; i++){
v = perm[permutation][i];
// Alloc
vector<unsigned int>* ordered = new vector<unsigned int>(2);
if(v == a || v == b || v == c){
(*inorder)[k] = v;
k++;
}
}
// Sort & Swap
sortAndSwap(permutedRefEdge[permutation][edge],
refEdge[edge],
*ordered);
return inorder;
// 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){
getOrderedTriFace(unsigned int permutation,
unsigned int face){
unsigned int opposite;
/******************************************************
* Same as getOrderedEdge, but I need to use 3 nodes. *
******************************************************/
unsigned int v;
unsigned int k = 0;
vector<unsigned int>* inorder =
new vector<unsigned int>(4);
// Alloc
vector<unsigned int>* ordered = new vector<unsigned int>(3);
// Take nodes in order //
for(unsigned int i = 0; i < nVertex; i++){
v = perm[permutation][i];
// Sort & Swap
sortAndSwap(permutedRefFace[permutation][face],
refFace[face],
*ordered);
if(v == a || v == b || v == c || v == d){
(*inorder)[k] = v;
k++;
// Return ordered
return ordered;
}
}
// 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] //
opposite = ((*inorder)[0] + 2) % 4;
if((*inorder)[2] != opposite){
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 //
......
......@@ -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
......@@ -19,7 +19,16 @@ TetReferenceSpace::TetReferenceSpace(void){
}
// Face Definition //
// 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{
......
......@@ -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){
......
......@@ -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;
......
......@@ -19,7 +19,14 @@ TriReferenceSpace::TriReferenceSpace(void){
}
// Face Definition //
// 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{
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment