diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 75f7f2c0747a02dad1ca7558a14568c3c01b1f6f..d73e2e07e457eb815bc5f571579407579105239f 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -145,9 +145,8 @@ BasisLocal* BasisGenerator::tetHierarchicalGen(size_t basisType,
BasisLocal* BasisGenerator::hexHierarchicalGen(size_t basisType,
size_t order){
switch(basisType){
- //case 0: return new HexNodeBasis(order);
+ case 0: return new HexNodeBasis(order);
//case 1: return new HexEdgeBasis(order);
- case 0: throw Exception("0-form not implemented on Hexs");
case 1: throw Exception("1-form not implemented on Hexs");
case 2: throw Exception("2-form not implemented on Hexs");
case 3: throw Exception("3-form not implemented on Hexs");
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index 5e418137dab6d07d1cd113b44a7c3c4d62664c2c..7f742391d5c09aa83a7642615737670c9899f96f 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -1,153 +1,126 @@
#include "HexNodeBasis.h"
+#include "HexReferenceSpace.h"
#include "Legendre.h"
using namespace std;
HexNodeBasis::HexNodeBasis(size_t order){
- /*
+ // Reference Space //
+ refSpace = new HexReferenceSpace;
+ nRefSpace = getReferenceSpace().getNReferenceSpace();
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = refSpace->getFaceNodeIndex();
+
// Set Basis Type //
this->order = order;
type = 0;
dim = 3;
- nVertex = 8;
- nEdge = 12 * (order - 1);
- nFace = 6 * (order - 1) * (order - 1);
- nCell = (order - 1) * (order - 1) * (order - 1);
-
- nEdgeClosure = 2;
- nFaceClosure = 8;
-
- size = nVertex + nEdge + nFace + nCell;
-
- // Alloc Temporary Space //
- Polynomial* legendre = new Polynomial[order];
- Polynomial* lifting = new Polynomial[8];
-
- Polynomial* xi = new Polynomial[6];
- Polynomial* eta = new Polynomial[6];
- Polynomial* lambda = new Polynomial[6];
+ nVertex = 8;
+ nEdge = 12 * (order - 1);
+ nFace = 6 * (order - 1) * (order - 1);
+ nCell = (order - 1) * (order - 1) * (order - 1);
+ nFunction = nVertex + nEdge + nFace + nCell;
// Legendre Polynomial //
+ Polynomial* legendre = new Polynomial[order];
Legendre::integrated(legendre, order);
- // Vertices definig Edges & Permutations //
- const int edgeV[2][12][2] =
- {
- { {0, 1}, {0, 3}, {0, 4}, {1, 2}, {1, 5}, {2, 3},
- {2, 6}, {3, 7}, {4, 5}, {4, 7}, {5, 6}, {6, 7} },
-
- { {1, 0}, {3, 0}, {4, 0}, {2, 1}, {5, 1}, {3, 2},
- {6, 2}, {7, 3}, {5, 4}, {7, 4}, {6, 5}, {7, 6} },
- };
-
- const int faceV[6][4] =
+ // Lagrange & Lifting //
+ const Polynomial lagrange[8] =
{
- {0, 3, 2, 1},
- {0, 1, 5, 4},
- {0, 4, 7, 3},
- {1, 2, 6, 5},
- {2, 3, 7, 6},
- {4, 5, 6, 7}
- };
-
-
- // Lifting //
- lifting[0] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- lifting[1] =
- (Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- lifting[2] =
- (Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 1, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- lifting[3] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 1, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- lifting[4] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
- Polynomial(1, 0, 0, 1);
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 1))),
- lifting[5] =
- (Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
- Polynomial(1, 0, 0, 1);
+ Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 1))),
- lifting[6] =
- (Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 1, 0)) +
- Polynomial(1, 0, 0, 1);
+ Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 1))),
- lifting[7] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 1, 0)) +
- Polynomial(1, 0, 0, 1);
-
-
- // Basis //
- basis = new std::vector<const Polynomial*>(size);
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 1)))
+ };
+ const Polynomial lifting[8] =
+ {
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- // Vertex Based (Lagrange) //
- (*basis)[0] =
- new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
+ Polynomial((Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- (*basis)[1] =
- new Polynomial((Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
+ Polynomial((Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- (*basis)[2] =
- new Polynomial((Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1))),
- (*basis)[3] =
- new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 1))),
- (*basis)[4] =
- new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1));
+ Polynomial((Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 1))),
- (*basis)[5] =
- new Polynomial((Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1));
+ Polynomial((Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 1))),
- (*basis)[6] =
- new Polynomial((Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1));
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 1)))
+ };
- (*basis)[7] =
- new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1));
+ // Basis //
+ basis = new Polynomial**[nRefSpace];
+ for(size_t s = 0; s < nRefSpace; s++)
+ basis[s] = new Polynomial*[nFunction];
+ // Vertex Based //
+ 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]);
+ basis[s][3] = new Polynomial(lagrange[3]);
+ basis[s][4] = new Polynomial(lagrange[4]);
+ basis[s][5] = new Polynomial(lagrange[5]);
+ basis[s][6] = new Polynomial(lagrange[6]);
+ basis[s][7] = new Polynomial(lagrange[7]);
+ }
+ /*
// Edge Based //
- // Keep counting
- int i = 8;
-
- // Points definig Edges
- const int edge1[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3};
- const int edge2[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7};
for(int l = 1; l < order; l++){
for(int e = 0; e < 12; e++){
@@ -233,40 +206,16 @@ HexNodeBasis::HexNodeBasis(size_t order){
}
HexNodeBasis::~HexNodeBasis(void){
- /*
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
-
-
- // Edge Based //
- for(int c = 0; c < 2; c++){
- for(int i = 0; i < nEdge; i++)
- delete (*(*edge)[c])[i];
-
- delete (*edge)[c];
- }
+ // ReferenceSpace //
+ delete refSpace;
- delete edge;
-
-
- // Face Based //
- for(int c = 0; c < 6; c++){
- for(int i = 0; i < nFace; i++)
- delete (*(*face)[c])[i];
+ // Basis //
+ for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t j = 0; j < nFunction; j++)
+ delete basis[i][j];
- delete (*face)[c];
+ delete[] basis[i];
}
- delete face;
-
-
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
-
- delete cell;
- */
+ delete[] basis;
}
diff --git a/FunctionSpace/PermutationTree.cpp b/FunctionSpace/PermutationTree.cpp
index 018ba7dd8b08381396dc0fe0cb4c5f416ab28961..651aba0e6a00a6d7eb7d05d000ac264e8117ad48 100644
--- a/FunctionSpace/PermutationTree.cpp
+++ b/FunctionSpace/PermutationTree.cpp
@@ -1,5 +1,4 @@
#include <map>
-#include <sstream>
#include "Exception.h"
#include "PermutationTree.h"
@@ -24,6 +23,7 @@ PermutationTree::PermutationTree(const vector<size_t>& refSequence){
root->tag = -1;
// Tree //
+ nextNodeId = 0;
populate(root, listOfLeaf);
// Leaf //
@@ -40,6 +40,10 @@ PermutationTree::~PermutationTree(void){
void PermutationTree::populate(node_t* node,
list<node_t*>& listOfLeaf){
+ // Node Id //
+ node->nodeId = nextNodeId;
+ nextNodeId++;
+
// Number of son //
const size_t nSon = node->nxtChoice.size();
@@ -93,11 +97,6 @@ void PermutationTree::destroy(node_t* node){
delete node;
}
-
-size_t PermutationTree::getPermutationId(const vector<size_t>& sequence) const{
- return getLeaf(root, sequence, 0)->leafId;
-}
-
PermutationTree::node_t*
PermutationTree::getLeaf(node_t* root,
const vector<size_t>& sequence,
@@ -142,6 +141,33 @@ void PermutationTree::fillWithPermutation(size_t permutationId,
}
}
+void PermutationTree::compressTag(void){
+ // Populate Tag map
+ const size_t nLeaf = leaf.size();
+ multimap<size_t, node_t*> tag;
+
+ for(size_t i = 0; i < nLeaf; i++)
+ tag.insert(pair<size_t, node_t*>(leaf[i]->tag, leaf[i]));
+
+ // Compress Tag
+ const multimap<size_t, node_t*>::iterator end = tag.end();
+ multimap<size_t, node_t*>::iterator it = tag.begin();
+
+ size_t currentTag = it->first;
+ size_t nextTag = 0;
+
+ for(; it != end; it++){
+ // If new old tag --> increase new tag by one
+ if(it->first != currentTag){
+ currentTag = it->first;
+ nextTag++;
+ }
+
+ // Apply same new tag
+ it->second->tag = nextTag;
+ }
+}
+
std::vector<std::pair<size_t, size_t> >
PermutationTree::getAllTagsCount(void) const{
// Init Temp map
@@ -166,6 +192,20 @@ PermutationTree::getAllTagsCount(void) const{
}
string PermutationTree::toString(void) const{
+ stringstream stream;
+
+ // Number of Node //
+ stream << "nNode " << nextNodeId
+ << endl;
+
+ // Nodes (depth first travel) //
+ serial_t* tmp = new serial_t;
+ appendString(stream, root, tmp);
+ delete tmp;
+
+ return stream.str();
+
+ /*
// Number of Permutation //
const size_t max_t = -1;
const size_t nPerm = leaf.size();
@@ -195,4 +235,69 @@ string PermutationTree::toString(void) const{
// Return //
return stream.str();
+ */
+}
+
+void PermutationTree::toSerial(node_t* node, serial_t* serial){
+ serial->myChoice = node->myChoice;
+ serial->nNxtChoice = node->nxtChoice.size();
+ serial->nxtChoice.resize(serial->nNxtChoice);
+
+ for(size_t i = 0; i < serial->nNxtChoice; i++)
+ serial->nxtChoice[i] = node->nxtChoice[i];
+
+ if(!node->father)
+ serial->fatherId = (size_t)(-1);
+ else
+ serial->fatherId = node->father->nodeId;
+
+ serial->nSon = node->son.size();
+ serial->sonId.resize(serial->nSon);
+
+ for(size_t i = 0; i < serial->nSon; i++)
+ serial->sonId[i] = node->son[i]->nodeId;
+
+ serial->nodeId = node->nodeId;
+ serial->leafId = node->leafId;
+ serial->tag = node->tag;
+}
+
+void PermutationTree::appendString(stringstream& stream,
+ node_t* node,
+ serial_t* tmp){
+ // Write this node //
+ toSerial(node, tmp);
+ stream << toString(tmp) << endl;
+
+ // Write Son //
+ const size_t nSon = node->son.size();
+
+ for(size_t i = 0; i < nSon; i++)
+ appendString(stream, node->son[i], tmp);
+}
+
+string PermutationTree::toString(serial_t* serial){
+ stringstream stream;
+
+ stream << "Node " << serial->nodeId
+ << " - " << serial->leafId
+ << " - " << serial->tag
+ << endl
+
+ << "MyChoice " << serial->myChoice
+ << endl
+ << "NxtChoice " << serial->nNxtChoice;
+
+ for(size_t i = 0; i < serial->nNxtChoice; i++)
+ stream << " - " << serial->nxtChoice[i];
+
+ stream << endl
+ << "FatherId " << serial->fatherId
+ << endl
+ << "SonId " << serial->nSon;
+
+ for(size_t i = 0; i < serial->nNxtChoice; i++)
+ stream << " - " << serial->sonId[i];
+
+ return stream.str();
}
diff --git a/FunctionSpace/PermutationTree.h b/FunctionSpace/PermutationTree.h
index f0cd7e99a6f90036afc5b679b01f007faced0714..e41b05d1d1625833bbbf2c2792354d7bfd345fae 100644
--- a/FunctionSpace/PermutationTree.h
+++ b/FunctionSpace/PermutationTree.h
@@ -5,6 +5,7 @@
#include <vector>
#include <list>
#include <string>
+#include <sstream>
/**
@class PermutationTree
@@ -29,11 +30,31 @@ class PermutationTree{
node_s* father;
std::vector<node_s*> son;
+ size_t nodeId;
size_t leafId;
size_t tag;
} node_t;
private:
+ // Tree node for serialization //
+ typedef struct serial_s{
+ size_t myChoice;
+ size_t nNxtChoice;
+ std::vector<size_t> nxtChoice;
+
+ size_t fatherId;
+ size_t nSon;
+ std::vector<size_t> sonId;
+
+ size_t nodeId;
+ size_t leafId;
+ size_t tag;
+ } serial_t;
+
+ private:
+ // Next Node Id & Total number of node //
+ size_t nextNodeId;
+
// Sequence size //
size_t sequenceSize;
@@ -56,21 +77,29 @@ class PermutationTree{
std::vector<size_t>& vectorToFill) const;
void addTagToPermutation(size_t permutationId, size_t tag);
- size_t getTagFromPermutation(size_t permutationId);
+ void compressTag(void);
+ size_t getTagFromPermutation(size_t permutationId) const;
std::vector<std::pair<size_t, size_t> > getAllTagsCount(void) const;
std::string toString(void) const;
private:
- static void populate(node_t* node,
- std::list<node_t*>& listOfLeaf);
+ void populate(node_t* node,
+ std::list<node_t*>& listOfLeaf);
static void destroy(node_t* node);
static node_t* getLeaf(node_t* root,
const std::vector<size_t>& sequence,
size_t offset);
+
+ static void toSerial(node_t* node, serial_t* serial);
+ static void appendString(std::stringstream& stream,
+ node_t* node,
+ serial_t* tmp);
+
+ static std::string toString(serial_t* serial);
};
/**
@@ -115,6 +144,18 @@ class PermutationTree{
Associates the given permutation with the given tag
**
+ @fn PermutationTree::compressTag
+
+ Takes all the tags of this tree and replace them
+ by their compressed version.
+
+ When tags are compressed, it means that their value
+ is mapped in the range 0, ..., (number of different tag - 1).
+
+ The mapping is such that the smaller/bigger relation
+ between the tags is preserved.
+ **
+
@fn PermutationTree::getTagFromPermutation
@param permutationId A permuted sequence ID
@return Returns the tag of the given sequence
@@ -144,13 +185,18 @@ inline size_t PermutationTree::getNPermutation(void) const{
return leaf.size();
}
+inline size_t PermutationTree::
+getPermutationId(const std::vector<size_t>& sequence) const{
+ return getLeaf(root, sequence, 0)->leafId;
+}
+
inline void PermutationTree::
addTagToPermutation(size_t permutationId, size_t tag){
leaf[permutationId]->tag = tag;
}
inline size_t PermutationTree::
-getTagFromPermutation(size_t permutationId){
+getTagFromPermutation(size_t permutationId) const{
return leaf[permutationId]->tag;
}
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index e3a3bd3062138ea49d8f3ed991376a4628e10995..7af6d330d24490cad4cd93dedb1224173cd553e5 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -48,6 +48,10 @@ void ReferenceSpace::init(void){
pTree->fillWithPermutation(*refSpaceIt, refSpaceNodeId[i]);
}
+ // Compress Tag, so that they refer to refSpaceNodeId indices
+ // insted of leafId of pTree.
+ pTree->compressTag();
+
// ABC space to UVW space shape function index mapping //
const size_t nPermutation = pTree->getNPermutation();
ABCtoUVWIndex.resize(nPermutation);
@@ -784,7 +788,7 @@ string ReferenceSpace::toString(void) const{
stream << "Reference Space Node IDs:" << endl;
for(size_t s = 0; s < nRefSpace; s++){
- stream << " ** Reference Space #" << s + 1 << ": [";
+ stream << " ** Reference Space #" << s << ": [";
for(size_t v = 0; v < nVertex - 1; v++)
stream << refSpaceNodeId[s][v] << ", ";
@@ -796,7 +800,7 @@ string ReferenceSpace::toString(void) const{
stream << endl << "Ordered Edge Node Index:" << endl;
for(size_t s = 0; s < nRefSpace; s++){
- stream << " ** Reference Space #" << s + 1 << ":"
+ stream << " ** Reference Space #" << s << ":"
<< endl;
for(size_t e = 0; e < nEdge; e++)
@@ -810,7 +814,7 @@ string ReferenceSpace::toString(void) const{
stream << endl << "Ordered Face Node Index:" << endl;
for(size_t s = 0; s < nRefSpace; s++){
- stream << " ** Reference Space #" << s + 1 << ":"
+ stream << " ** Reference Space #" << s << ":"
<< endl;
for(size_t f = 0; f < nFace; f++){