diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 6cd443357a881066c751f63b1deec58579ff3415..9c0ceb0e98e85f5c9b100725cc071336def3fafe 100644
--- a/FunctionSpace/Basis.h
+++ b/FunctionSpace/Basis.h
@@ -2,6 +2,7 @@
#define _BASIS_H_
#include "MElement.h"
+#include "ReferenceSpace.h"
/**
@interface Basis
@@ -43,6 +44,8 @@
class Basis{
protected:
+ ReferenceSpace* refSpace;
+
bool scalar;
bool local;
@@ -78,17 +81,11 @@ class Basis{
unsigned int getNFunction(void) const;
// Reference Element //
+ const ReferenceSpace& getReferenceSpace(void) const;
+
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;
-
// Direct Access to Evaluated Functions //
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
@@ -115,6 +112,8 @@ class Basis{
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(unsigned int orientation) const = 0;
+ virtual std::string toString(void) const = 0;
+
protected:
// 'Constructor' //
Basis(void);
@@ -327,4 +326,8 @@ inline unsigned int Basis::getNFunction(void) const{
return nFunction;
}
+inline const ReferenceSpace& Basis::getReferenceSpace(void) const{
+ return *refSpace;
+}
+
#endif
diff --git a/FunctionSpace/BasisHierarchical0From.cpp b/FunctionSpace/BasisHierarchical0From.cpp
index cc3a579c4c40fd2f5f7f2711945c4aeca7aba402..7f5a1597f5b5f5e17aa1890f13edd99cae3e178e 100644
--- a/FunctionSpace/BasisHierarchical0From.cpp
+++ b/FunctionSpace/BasisHierarchical0From.cpp
@@ -59,49 +59,6 @@ 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);
-
- 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;
-}
-
void BasisHierarchical0From::
getFunctions(fullMatrix<double>& retValues,
const MElement& element,
diff --git a/FunctionSpace/BasisHierarchical0From.h b/FunctionSpace/BasisHierarchical0From.h
index 96377520d72d392c672e0e6d9889b018f396593d..b2ee85c0edf9fbfdb6d05a87cfc33421ba41c763 100644
--- a/FunctionSpace/BasisHierarchical0From.h
+++ b/FunctionSpace/BasisHierarchical0From.h
@@ -16,8 +16,7 @@
class BasisHierarchical0From: public BasisLocal{
protected:
// Orientation //
- ReferenceSpace* refSpace;
- unsigned int nRefSpace;
+ unsigned int nRefSpace;
// Basis //
Polynomial*** basis;
@@ -39,13 +38,6 @@ 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;
-
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
diff --git a/FunctionSpace/BasisHierarchical1From.cpp b/FunctionSpace/BasisHierarchical1From.cpp
index fa94def81bbf88b53f34abd806d41e8cea0391c5..ea8e5442e1659fdb1cf4b9ab6bc710ddead33ffa 100644
--- a/FunctionSpace/BasisHierarchical1From.cpp
+++ b/FunctionSpace/BasisHierarchical1From.cpp
@@ -59,22 +59,6 @@ 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);
-
- for(size_t i = 0; i < nFunction; i++)
- ordering[i] = i;
-
- return ordering;
-}
-
void BasisHierarchical1From::
getFunctions(fullMatrix<double>& retValues,
const MElement& element,
diff --git a/FunctionSpace/BasisHierarchical1From.h b/FunctionSpace/BasisHierarchical1From.h
index f72d061ef5c28ea3bed87dbacd557237b2f29354..9c67fca2c4931723fac29193c5d11151361ef5c9 100644
--- a/FunctionSpace/BasisHierarchical1From.h
+++ b/FunctionSpace/BasisHierarchical1From.h
@@ -16,8 +16,7 @@
class BasisHierarchical1From: public BasisLocal{
protected:
// Orientation //
- ReferenceSpace* refSpace;
- unsigned int nRefSpace;
+ unsigned int nRefSpace;
// Basis //
std::vector<Polynomial>*** basis;
@@ -39,13 +38,6 @@ 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;
-
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index 592af30b3c9eaec9969623cda5b69011581745ea..7e5584cd16160cac62a7e3124e7a18e83d0801a9 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -29,7 +29,7 @@ unsigned int BasisLagrange::
getOrientation(const MElement& element) const{
return 0;
}
-
+/*
static bool
sortPredicate(const std::pair<size_t, size_t>& a,
const std::pair<size_t, size_t>& b){
@@ -85,7 +85,7 @@ void BasisLagrange::mapFromXYZtoABC(const MElement& element,
vector<size_t> BasisLagrange::
getFunctionOrdering(const MElement& element) const{
- /*
+
static bool once = true;
if(once){
@@ -99,11 +99,11 @@ getFunctionOrdering(const MElement& element) const{
once = false;
}
- */
+
vector<int> rNodeId = reducedNodeId(element);
const size_t closureId = matchClosure(rNodeId, lBasis->fullClosures);
- /*
+
vector<int>& closure = lBasis->fullClosures[closureId];
vector<size_t> myClosure(closure.size());
@@ -112,7 +112,7 @@ getFunctionOrdering(const MElement& element) const{
myClosure[i] = closure[i];
return myClosure;
- */
+
vector<size_t> c(10);
if(closureId == 0){
@@ -143,7 +143,7 @@ getFunctionOrdering(const MElement& element) const{
return c;
}
-
+*/
void BasisLagrange::
getFunctions(fullMatrix<double>& retValues,
const MElement& element,
@@ -307,3 +307,7 @@ project(const MElement& element,
// Return ;
return newCoef;
}
+
+std::string BasisLagrange::toString(void) const{
+ return std::string("BasisLagrange::toString() not Implemented");
+}
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index 23cd49da568d2d821051d0728e884be44b62003f..913be52da98aa8d617ccc5f834350eb781869c98 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/FunctionSpace/BasisLagrange.h
@@ -43,13 +43,6 @@ 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;
-
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
@@ -86,6 +79,8 @@ class BasisLagrange: public BasisLocal{
const std::vector<double>& coef,
const FunctionSpaceVector& fSpace);
+ virtual std::string toString(void) const;
+
protected:
BasisLagrange(void);
};
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index bc3488878c10072737750b2b2b7aca06f76be29c..ef98d9eae5c8fa0d40f98f71e3b4fc20c42dede9 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -114,19 +114,8 @@ void FunctionSpace::buildDof(void){
for(unsigned int j = 0; j < nDof; j++)
insertDof(myDof[j], trueDof, j);
- // Dof Orderning
- vector<size_t> dofOrder = (*basis)[0]->getFunctionOrdering(*(element[i]));
- vector<const Dof*> trueDofOrdered(nDof);
-
- for(size_t j = 0; j < nDof; j++)
- trueDofOrdered[dofOrder[j]] = trueDof[j];
- /*
- for(size_t j = 0; j < nDof; j++)
- cout << dofOrder[j] << "\t";
- cout << endl;
- */
// Create new GroupOfDof
- GroupOfDof* god = new GroupOfDof(*(element[i]), trueDofOrdered);
+ GroupOfDof* god = new GroupOfDof(*(element[i]), trueDof);
(*group)[i] = god;
// Map GOD
@@ -162,26 +151,39 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
// Const_Cast //
MElement& element = const_cast<MElement&>(elem);
- // Get Element Data //
- const unsigned int nVertex = element.getNumPrimaryVertices();
- const unsigned int nEdge = element.getNumEdges();
- const unsigned int nFace = element.getNumFaces();
+ // Create New Element With Permuted Vertices //
+ // Permutation
+ const vector<size_t>& vPerm =
+ (*this->basis)[0]->getReferenceSpace().getNodeIndexFromABCtoUVW(elem);
+ // Permuted Vertices
+ const size_t nVertex = element.getNumPrimaryVertices();
vector<MVertex*> vertex(nVertex);
+
+ for(size_t i = 0; i < nVertex; i++)
+ vertex[i] = element.getVertex(vPerm[i]);
+
+ // New Element
+ MElementFactory factory;
+ MElement* permElement = factory.create(elem.getTypeForMSH(), vertex);
+
+ // Edge & Face from Permuted Element //
+ const size_t nEdge = permElement->getNumEdges();
+ const size_t nFace = permElement->getNumFaces();
+
vector<MEdge> edge(nEdge);
vector<MFace> face(nFace);
- for(unsigned int i = 0; i < nVertex; i++)
- vertex[i] = element.getVertex(i);
+ for(size_t i = 0; i < nEdge; i++)
+ edge[i] = permElement->getEdge(i);
- for(unsigned int i = 0; i < nEdge; i++)
- edge[i] = element.getEdge(i);
+ for(size_t i = 0; i < nFace; i++)
+ face[i] = permElement->getFace(i);
- for(unsigned int i = 0; i < nFace; i++)
- face[i] = element.getFace(i);
+ delete permElement;
// Create Dof //
- unsigned int nDof =
+ size_t nDof =
fPerVertex * nVertex +
fPerEdge * nEdge +
fPerFace * nFace +
@@ -189,34 +191,34 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
vector<Dof> myDof(nDof);
- unsigned int it = 0;
+ size_t it = 0;
// Add Vertex Based Dof //
- for(unsigned int i = 0; i < nVertex; i++){
- for(unsigned int j = 0; j < fPerVertex; j++){
+ for(size_t i = 0; i < nVertex; i++){
+ for(size_t j = 0; j < fPerVertex; j++){
myDof[it].setDof(mesh->getGlobalId(*vertex[i]), j);
it++;
}
}
// Add Edge Based Dof //
- for(unsigned int i = 0; i < nEdge; i++){
- for(unsigned int j = 0; j < fPerEdge; j++){
+ for(size_t i = 0; i < nEdge; i++){
+ for(size_t j = 0; j < fPerEdge; j++){
myDof[it].setDof(mesh->getGlobalId(edge[i]), j);
it++;
}
}
// Add Face Based Dof //
- for(unsigned int i = 0; i < nFace; i++){
- for(unsigned int j = 0; j < fPerFace; j++){
+ for(size_t i = 0; i < nFace; i++){
+ for(size_t j = 0; j < fPerFace; j++){
myDof[it].setDof(mesh->getGlobalId(face[i]), j);
it++;
}
}
// Add Cell Based Dof //
- for(unsigned int j = 0; j < fPerCell; j++){
+ for(size_t j = 0; j < fPerCell; j++){
myDof[it].setDof(mesh->getGlobalId(element), j);
it++;
}
@@ -224,9 +226,11 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
return myDof;
}
-const GroupOfDof& FunctionSpace::getGoDFromElement(const MElement& element) const{
- const map<const MElement*, const GroupOfDof*, ElementComparator>::iterator it =
- eToGod->find(&element);
+const GroupOfDof& FunctionSpace::
+getGoDFromElement(const MElement& element) const{
+
+ const map<const MElement*, const GroupOfDof*, ElementComparator>::iterator
+ it = eToGod->find(&element);
if(it == eToGod->end())
throw
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 8e888e9c49a12cc6398bda1f81199f50329a9c1c..1c294d6a85add94ad1abb0b82f0298075c32beed 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -14,87 +14,10 @@ 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);
-
- // Get Reference coordinate //
- double phys[3] = {xyz(0), xyz(1), xyz(2)};
- double uvw[3];
- double abc[3];
-
- eelement.xyz2uvw(phys, uvw);
-
- // Compute coordinate in abc Space //
- std::vector<size_t> order = (*basis)[0]->getFunctionOrdering(element);
- double abcMat[3][3];
- double abcPnt[3][3];
-
- element.getNode(0, abcPnt[0][0], abcPnt[0][1], abcPnt[0][2]);
- element.getNode(1, abcPnt[1][0], abcPnt[1][1], abcPnt[1][2]);
- element.getNode(2, abcPnt[2][0], abcPnt[2][1], abcPnt[2][2]);
-
- std::cout << order[0] << " | "
- << order[1] << " | "
- << order[2] << std::endl;
-
- abcMat[0][0] = abcPnt[0][0];
- abcMat[0][1] = abcPnt[1][0];
- abcMat[0][2] = abcPnt[2][0];
-
- abcMat[1][0] = abcPnt[0][1];
- abcMat[1][1] = abcPnt[1][1];
- abcMat[1][2] = abcPnt[2][1];
-
- abcMat[2][0] = abcPnt[0][2];
- abcMat[2][1] = abcPnt[1][2];
- abcMat[2][2] = abcPnt[2][2];
-
- for(size_t i = 0; i < 3; i++){
- for(size_t j = 0; j < 3; j++)
- std::cout << abcMat[i][j] << "\t";
- std::cout << std::endl;
- }
- std::cout << std::endl;
-
- double phiUVW[3];
- element.getShapeFunctions(uvw[0], uvw[1], uvw[2], phiUVW);
-
- size_t trueOrder[3];
-
- for(size_t i = 0; i < 3; i++){
- size_t idx;
-
- for(idx = 0; order[idx] != i; idx++)
- ;
-
- trueOrder[i] = idx;
- }
-
- abc[0] =
- abcMat[0][0] * phiUVW[trueOrder[0]] +
- abcMat[0][1] * phiUVW[trueOrder[1]] +
- abcMat[0][2] * phiUVW[trueOrder[2]];
-
- abc[1] =
- abcMat[1][0] * phiUVW[trueOrder[0]] +
- abcMat[1][1] * phiUVW[trueOrder[1]] +
- abcMat[1][2] * phiUVW[trueOrder[2]];
-
- abc[2] =
- abcMat[2][0] * phiUVW[trueOrder[0]] +
- abcMat[2][1] * phiUVW[trueOrder[1]] +
- abcMat[2][2] * phiUVW[trueOrder[2]];
-
- uvw[0] = abc[0];
- uvw[1] = abc[1];
- uvw[2] = abc[2];
- */
// Get ABC Space coordinate //
double abc[3];
- (*basis)[0]->mapFromXYZtoABC(element, xyz, abc);
+ (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element, xyz, abc);
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
@@ -117,11 +40,15 @@ interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
+ // Get ABC Space coordinate //
+ double abc[3];
+ (*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element, uvw, 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/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 1e3c333d897c4647a5df1fc1c9097e0bb994fbe7..f75cdb87f6058ce8bff258299f3f250b59aa7d37 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -4,17 +4,16 @@
using namespace std;
-QuadNodeBasis::QuadNodeBasis(unsigned int order){
- /*
+QuadNodeBasis::QuadNodeBasis(size_t order){
// Reference Space //
refSpace = new QuadReferenceSpace;
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 Basis Type //
this->order = order;
@@ -66,11 +65,11 @@ QuadNodeBasis::QuadNodeBasis(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]);
@@ -78,17 +77,17 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
}
// 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 < 4; e++){
- for(unsigned int l = 1; l < order; l++){
+ for(size_t e = 0; e < 4; e++){
+ for(size_t l = 1; l < order; l++){
basis[s][i] =
- new Polynomial(legendre[l].compose(lifting[(*(*edgeV[s])[e])[1]] -
- lifting[(*(*edgeV[s])[e])[0]])
+ new Polynomial(legendre[l].compose(lifting[edgeIdx[s][e][1]] -
+ lifting[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++;
}
@@ -97,22 +96,22 @@ QuadNodeBasis::QuadNodeBasis(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'
- 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(unsigned int l1 = 1; l1 < order; l1++){
- for(unsigned int l2 = 1; l2 < order; l2++){
+ for(size_t l1 = 1; l1 < order; l1++){
+ for(size_t l2 = 1; l2 < order; l2++){
basis[s][i] =
- new Polynomial(legendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[1]])
+ new Polynomial(legendre[l1].compose(lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][1]])
*
- legendre[l2].compose(lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[3]]));
+ legendre[l2].compose(lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][3]]));
i++;
}
@@ -132,8 +131,8 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0));
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int i = 0; i < nFunction; i++){
+ for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t i = 0; i < nFunction; i++){
Polynomial* tmp;
tmp = basis[s][i];
basis[s][i] = new Polynomial(tmp->compose(mapX, mapY));
@@ -143,22 +142,19 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
// Free Temporary Sapce //
delete[] legendre;
- */
}
QuadNodeBasis::~QuadNodeBasis(void){
- /*
// ReferenceSpace //
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];
}
delete[] basis;
- */
}
diff --git a/FunctionSpace/QuadNodeBasis.h b/FunctionSpace/QuadNodeBasis.h
index 7bc159aa1d63c4cf9559fdf2a04f3cdb2ee4c85d..96c69c6d6e231a81d7b7093001d72bc0bea23a75 100644
--- a/FunctionSpace/QuadNodeBasis.h
+++ b/FunctionSpace/QuadNodeBasis.h
@@ -24,7 +24,7 @@ class QuadNodeBasis: public BasisHierarchical0From{
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Quads of the given order
- QuadNodeBasis(unsigned int order);
+ QuadNodeBasis(size_t order);
//! @return Deletes this Basis
//!
diff --git a/FunctionSpace/QuadReferenceSpace.cpp b/FunctionSpace/QuadReferenceSpace.cpp
index 111ed806e8cf9ca4a8db551decf64b6d9fb29a1b..47495dd8d718e063b1e65e4d0f3b0d504ab6cf53 100644
--- a/FunctionSpace/QuadReferenceSpace.cpp
+++ b/FunctionSpace/QuadReferenceSpace.cpp
@@ -5,64 +5,39 @@
using namespace std;
QuadReferenceSpace::QuadReferenceSpace(void){
- /*
// Vertex Definition //
nVertex = 4;
// Edge Definition //
- nEdge = 4;
- refEdge = new size_t*[nEdge];
+ const size_t nEdge = 4;
+ refEdgeNodeIdx.resize(nEdge);
for(size_t i = 0; i < nEdge; i++){
- refEdge[i] = new size_t[2];
- refEdge[i][0] = MQuadrangle::edges_quad(i, 0);
- refEdge[i][1] = MQuadrangle::edges_quad(i, 1);
+ refEdgeNodeIdx[i].resize(2); // Two Nodes per Edge
+ refEdgeNodeIdx[i][0] = MQuadrangle::edges_quad(i, 0);
+ refEdgeNodeIdx[i][1] = MQuadrangle::edges_quad(i, 1);
}
// Face Definition //
- // Number of face
- nFace = 1;
+ refFaceNodeIdx.resize(1); // One Face per Triangle
+ refFaceNodeIdx[0].resize(4); // Four Nodes per Face
- // Number of node per face
- nNodeInFace = new size_t[nFace];
- nNodeInFace[0] = 4;
-
- // Reference Face
- refFace = new size_t*[nFace];
- refFace[0] = new size_t[4];
-
- refFace[0][0] = 0;
- refFace[0][1] = 1;
- refFace[0][2] = 2;
- refFace[0][3] = 3;
+ refFaceNodeIdx[0][0] = 0;
+ refFaceNodeIdx[0][1] = 1;
+ refFaceNodeIdx[0][2] = 2;
+ refFaceNodeIdx[0][3] = 3;
// Init All //
init();
- */
}
QuadReferenceSpace::~QuadReferenceSpace(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 QuadReferenceSpace::toLatex(void) const{
- //const size_t nPerm = pTree->getNPermutation();
- stringstream stream;
- //vector<size_t> perm(nVertex);
- /*
+ const size_t nRefSpace = refSpaceNodeId.size();
+ stringstream stream;
+
stream << "\\documentclass{article}" << endl << endl
<< "\\usepackage{longtable}" << endl
@@ -71,29 +46,32 @@ string QuadReferenceSpace::toLatex(void) const{
<< "\\begin{document}" << endl
<< "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
- << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
+ << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl
+ << endl
<< "\\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(3, 3) {$" << perm[2] << "$};" << endl
- << "\\node[vertex] (n3) at(0, 3) {$" << perm[3] << "$};" << 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(3, 3) {$" << refSpaceNodeId[s][2] << "$};"
+ << endl
+ << "\\node[vertex] (n3) at(0, 3) {$" << refSpaceNodeId[s][3] << "$};"
+ << endl
<< endl;
- for(size_t i = 0; i < 4; i++)
+ for(size_t e = 0; e < 4; 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
@@ -102,6 +80,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 de56c4278e15b7df92106c676b442642d790a222..ed14cfcdd821204ebb4178170be5583a59376535 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -125,56 +125,35 @@ findCyclicPermutation(list<size_t>& listOfTrueReferenceSpace,
}
}
-static bool isFacePermutation(vector<size_t>& refNode,
- vector<size_t>& testNode){
- const size_t size = refNode.size();
- bool match = false;
-
- if(size != testNode.size())
- return false;
-
- for(size_t i = 0; i < size && !match; i++){
- bool submatch = true;
-
- for(size_t j = 0; j < size && submatch; j++)
- if(refNode[j] != testNode[j])
- submatch = false;
-
- if(!submatch){
- size_t tmp0 = testNode[0];
- size_t tmp1 = testNode[1];
-
- for(size_t k = 1; k < size + 1; k++){
- testNode[k % size] = tmp0;
- tmp0 = tmp1;
- tmp1 = testNode[(k + 1) % size];
- }
- }
+ReferenceSpace::triple ReferenceSpace::
+isCyclicPermutation(vector<size_t>& pTest,
+ vector<size_t>& pRef){
- else
- match = true;
- }
+ // Node IDs of Reference Space first Face
+ const size_t nNodeInFaceZero = refFaceNodeIdx[0].size();
+ vector<size_t> refNode(nNodeInFaceZero);
- return match;
-}
+ for(size_t i = 0; i < nNodeInFaceZero; i++)
+ refNode[i] = pRef[refFaceNodeIdx[0][i]];
-static bool haveSameNode(vector<size_t>& face0,
- vector<size_t>& face1){
- const size_t size = face0.size();
- bool matchIsPossible = true;
+ // Corresponding Face in Test Permutation
+ size_t testFaceId = findCorrespondingFace(refNode, pTest);
- for(size_t i = 0; i < size && matchIsPossible; i++){
- bool submatch = false;
+ // Node IDs of Test Permutation correspnding Face
+ const size_t nNodeInFaceTest = refFaceNodeIdx[testFaceId].size();
+ vector<size_t> testNode(nNodeInFaceTest);
- for(size_t j = 0; j < size && !submatch; j++){
- if(face0[i] == face1[j])
- submatch = true;
- }
+ for(size_t i = 0; i < nNodeInFaceTest; i++)
+ testNode[i] = pTest[refFaceNodeIdx[testFaceId][i]];
- matchIsPossible = submatch;
- }
+ // Return Triple //
+ triple tri = {
+ isFacePermutation(refNode, testNode),
+ getRefIndexPermutation(pRef, pTest),
+ getReverseIndexPermutation(pRef, pTest)
+ };
- return matchIsPossible;
+ return tri;
}
size_t ReferenceSpace::findCorrespondingFace(vector<size_t>& face,
@@ -208,8 +187,41 @@ size_t ReferenceSpace::findCorrespondingFace(vector<size_t>& face,
return f;
}
-static vector<size_t> getRefIndexPermutation(vector<size_t>& ref,
- vector<size_t>& test){
+bool ReferenceSpace::isFacePermutation(vector<size_t>& refNode,
+ vector<size_t>& testNode){
+ const size_t size = refNode.size();
+ bool match = false;
+
+ if(size != testNode.size())
+ return false;
+
+ for(size_t i = 0; i < size && !match; i++){
+ bool submatch = true;
+
+ for(size_t j = 0; j < size && submatch; j++)
+ if(refNode[j] != testNode[j])
+ submatch = false;
+
+ if(!submatch){
+ size_t tmp0 = testNode[0];
+ size_t tmp1 = testNode[1];
+
+ for(size_t k = 1; k < size + 1; k++){
+ testNode[k % size] = tmp0;
+ tmp0 = tmp1;
+ tmp1 = testNode[(k + 1) % size];
+ }
+ }
+
+ else
+ match = true;
+ }
+
+ return match;
+}
+
+vector<size_t> ReferenceSpace::getRefIndexPermutation(vector<size_t>& ref,
+ vector<size_t>& test){
const size_t size = ref.size();
vector<size_t> idxVec(ref.size());
size_t idx;
@@ -226,8 +238,8 @@ static vector<size_t> getRefIndexPermutation(vector<size_t>& ref,
return idxVec;
}
-static vector<size_t> getReverseIndexPermutation(vector<size_t>& ref,
- vector<size_t>& test){
+vector<size_t> ReferenceSpace::getReverseIndexPermutation(vector<size_t>& ref,
+ vector<size_t>& test){
const size_t size = ref.size();
vector<size_t> idxVec(ref.size());
size_t idx;
@@ -244,35 +256,23 @@ static vector<size_t> getReverseIndexPermutation(vector<size_t>& ref,
return idxVec;
}
-ReferenceSpace::triple ReferenceSpace::
-isCyclicPermutation(vector<size_t>& pTest,
- vector<size_t>& pRef){
-
- // Node IDs of Reference Space first Face
- const size_t nNodeInFaceZero = refFaceNodeIdx[0].size();
- vector<size_t> refNode(nNodeInFaceZero);
-
- 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);
+bool ReferenceSpace::haveSameNode(vector<size_t>& face0,
+ vector<size_t>& face1){
+ const size_t size = face0.size();
+ bool matchIsPossible = true;
- // Node IDs of Test Permutation correspnding Face
- const size_t nNodeInFaceTest = refFaceNodeIdx[testFaceId].size();
- vector<size_t> testNode(nNodeInFaceTest);
+ for(size_t i = 0; i < size && matchIsPossible; i++){
+ bool submatch = false;
- for(size_t i = 0; i < nNodeInFaceTest; i++)
- testNode[i] = pTest[refFaceNodeIdx[testFaceId][i]];
+ for(size_t j = 0; j < size && !submatch; j++){
+ if(face0[i] == face1[j])
+ submatch = true;
+ }
- // Return Triple //
- triple tri = {
- isFacePermutation(refNode, testNode),
- getRefIndexPermutation(pRef, pTest),
- getReverseIndexPermutation(pRef, pTest)
- };
+ matchIsPossible = submatch;
+ }
- return tri;
+ return matchIsPossible;
}
void ReferenceSpace::getOrderedEdge(void){
@@ -421,13 +421,25 @@ size_t ReferenceSpace::getPermutationIdx(const MElement& element) const{
void ReferenceSpace::mapFromXYZtoABC(const MElement& element,
const fullVector<double>& xyz,
- double abc[3]){
+ double abc[3]) const{
// Get UVW coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
double uvw[3];
element.xyz2uvw(phys, uvw);
+ // Get ABC coordinate //
+ fullVector<double> uuvw(3);
+ uuvw(0) = uvw[0];
+ uuvw(1) = uvw[1];
+ uuvw(2) = uvw[2];
+
+ mapFromUVWtoABC(element, uuvw, abc);
+}
+
+void ReferenceSpace::mapFromUVWtoABC(const MElement& element,
+ const fullVector<double>& uvw,
+ double abc[3]) const{
// Compute coordinate in ABC Space //
// ABC node coordinate
double** abcMat = new double*[nVertex];
@@ -440,7 +452,7 @@ void ReferenceSpace::mapFromXYZtoABC(const MElement& element,
// UVW (order 1) shape functions
double* phiUVW = new double[nVertex];
- element.getShapeFunctions(uvw[0], uvw[1], uvw[2], phiUVW, 1);
+ element.getShapeFunctions(uvw(0), uvw(1), uvw(2), phiUVW, 1);
// Element Permutation Index //
const size_t permutationIdx = getPermutationIdx(element);
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index ffdf867992f480b3999b604da35fd336122420f8..e951efb3b145fa902faef160301d16d908f30787 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -55,18 +55,22 @@ class ReferenceSpace{
size_t getNReferenceSpace(void) const;
size_t getReferenceSpace(const MElement& element) const;
- const std::vector<std::vector<size_t> >&
- getNodeId(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;
+ const std::vector<size_t>&
+ getNodeIndexFromABCtoUVW(const MElement& element) const;
+
+ void mapFromUVWtoABC(const MElement& element,
+ const fullVector<double>& xyz,
+ double abc[3]) const;
+
void mapFromXYZtoABC(const MElement& element,
const fullVector<double>& xyz,
- double abc[3]);
+ double abc[3]) const;
void getJacobian(const MElement& element,
const fullVector<double>& xyz,
@@ -80,6 +84,7 @@ class ReferenceSpace{
void init(void);
+ private:
void getOrderedEdge(void);
void getOrderedFace(void);
@@ -94,6 +99,20 @@ class ReferenceSpace{
size_t findCorrespondingFace(std::vector<size_t>& face,
std::vector<size_t>& node);
+ static bool isFacePermutation(std::vector<size_t>& refNode,
+ std::vector<size_t>& testNode);
+
+ std::vector<size_t> getRefIndexPermutation(std::vector<size_t>& ref,
+ std::vector<size_t>& test);
+
+ std::vector<size_t> getReverseIndexPermutation(std::vector<size_t>& ref,
+ std::vector<size_t>& test);
+
+ static bool haveSameNode(std::vector<size_t>& face0,
+ std::vector<size_t>& face1);
+
+ size_t getPermutationIdx(const MElement& element) const;
+
static void
orderRefEntityForGivenRefSpace(std::vector<size_t>& refEntityNodeIdx,
std::vector<size_t>& refSpaceNodeId,
@@ -102,8 +121,6 @@ class ReferenceSpace{
static void
correctQuadFaceNodeIdx(std::vector<size_t>& correctedQuadFaceNodeIdx);
- size_t getPermutationIdx(const MElement& element) const;
-
static bool sortPredicate(const std::pair<size_t, size_t>& a,
const std::pair<size_t, size_t>& b);
};
@@ -188,6 +205,12 @@ ReferenceSpace::getFaceNodeIndex(void) const{
return orderedFaceNodeIdx;
}
+inline
+const std::vector<size_t>&
+ReferenceSpace::getNodeIndexFromABCtoUVW(const MElement& element) const{
+ return ABCtoUVWIndex[getPermutationIdx(element)];
+}
+
inline
size_t ReferenceSpace::getReferenceSpace(const MElement& element) const{
return pTree->getTagFromPermutation(getPermutationIdx(element));
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 4af8626ad8b277d1acc17ec412c93fa710d8b939..c9854688a42b77c605fc85f5d13e87f95b69b88c 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -4,17 +4,16 @@
using namespace std;
-TetNodeBasis::TetNodeBasis(unsigned int order){
- /*
+TetNodeBasis::TetNodeBasis(size_t order){
// Reference Space //
refSpace = new TetReferenceSpace;
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 Basis Type //
this->order = order;
@@ -61,11 +60,11 @@ TetNodeBasis::TetNodeBasis(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]);
@@ -73,18 +72,18 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
}
// 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(int e = 0; e < 6; e++){
- for(unsigned int l = 1; l < order; l++){
+ for(size_t l = 1; l < order; l++){
basis[s][i] =
new Polynomial(intLegendre[l].compose
- (lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
+ (lagrange[edgeIdx[s][e][0]] -
+ lagrange[edgeIdx[s][e][1]]
,
- lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]]));
+ lagrange[edgeIdx[s][e][0]] +
+ lagrange[edgeIdx[s][e][1]]));
i++;
}
}
@@ -93,31 +92,31 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
// Face Based //
// TO CHECK: Are Triangles face matching tets ?
- 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 f = 0; f < 4; f++){
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
Polynomial sum =
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]] +
- lagrange[(*(*faceV[s])[f])[2]];
+ lagrange[faceIdx[s][f][0]] +
+ lagrange[faceIdx[s][f][1]] +
+ lagrange[faceIdx[s][f][2]];
basis[s][i] =
new Polynomial(intLegendre[l1].compose
- (lagrange[(*(*faceV[s])[f])[0]] -
- lagrange[(*(*faceV[s])[f])[1]]
+ (lagrange[faceIdx[s][f][0]] -
+ lagrange[faceIdx[s][f][1]]
,
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]])
+ lagrange[faceIdx[s][f][0]] +
+ lagrange[faceIdx[s][f][1]])
*
- lagrange[(*(*faceV[s])[f])[2]]
+ lagrange[faceIdx[s][f][2]]
*
sclLegendre[l2].compose
- (lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum));
+ (lagrange[faceIdx[s][f][2]] * 2 - sum, sum));
i++;
}
}
@@ -132,8 +131,8 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
const Polynomial sub = lagrange[0] - lagrange[1];
const Polynomial add = lagrange[0] + lagrange[1];
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nVertex + nEdge + nFace;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nVertex + nEdge + nFace;
for(int l1 = 1; l1 < orderMinusTwo; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
@@ -155,22 +154,19 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
- */
}
TetNodeBasis::~TetNodeBasis(void){
- /*
// ReferenceSpace //
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];
}
delete[] basis;
- */
}
diff --git a/FunctionSpace/TetNodeBasis.h b/FunctionSpace/TetNodeBasis.h
index b1ffb219306f65b57aeeab3e04076089efe62058..3eff8cf0754aeb284ef05d259484aa676096b458 100644
--- a/FunctionSpace/TetNodeBasis.h
+++ b/FunctionSpace/TetNodeBasis.h
@@ -20,7 +20,7 @@ class TetNodeBasis: public BasisHierarchical0From{
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Tetrahedra of the given order
- TetNodeBasis(unsigned int order);
+ TetNodeBasis(size_t order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TetReferenceSpace.cpp b/FunctionSpace/TetReferenceSpace.cpp
index 839f956229ad01e959fe4e4419faac98d38de60b..402d9104df9feed0af1cb70314c4e2a7c2e7c7df 100644
--- a/FunctionSpace/TetReferenceSpace.cpp
+++ b/FunctionSpace/TetReferenceSpace.cpp
@@ -5,67 +5,41 @@
using namespace std;
TetReferenceSpace::TetReferenceSpace(void){
- /*
// Vertex Definition //
nVertex = 4;
// Edge Definition //
- nEdge = 6;
- refEdge = new size_t*[nEdge];
+ const size_t nEdge = 6;
+ refEdgeNodeIdx.resize(nEdge);
for(size_t i = 0; i < nEdge; i++){
- refEdge[i] = new size_t[2];
- refEdge[i][0] = MTetrahedron::edges_tetra(i, 0);
- refEdge[i][1] = MTetrahedron::edges_tetra(i, 1);
+ refEdgeNodeIdx[i].resize(2); // Two Nodes per Edge
+ refEdgeNodeIdx[i][0] = MTetrahedron::edges_tetra(i, 0);
+ refEdgeNodeIdx[i][1] = MTetrahedron::edges_tetra(i, 1);
}
// Face Definition //
- // Number of face
- nFace = 4;
-
- // Number of node per face
- nNodeInFace = new size_t[nFace];
-
- for(size_t f = 0; f < nFace; f++)
- nNodeInFace[f] = 3;
-
- // Reference Face
- refFace = new size_t*[nFace];
+ size_t nFace = 4;
+ refFaceNodeIdx.resize(nFace);
for(size_t i = 0; i < nFace; i++){
- refFace[i] = new size_t[3];
- refFace[i][0] = MTetrahedron::faces_tetra(i, 0);
- refFace[i][1] = MTetrahedron::faces_tetra(i, 1);
- refFace[i][2] = MTetrahedron::faces_tetra(i, 2);
+ refFaceNodeIdx[i].resize(3); // Three Nodes per Face
+ refFaceNodeIdx[i][0] = MTetrahedron::faces_tetra(i, 0);
+ refFaceNodeIdx[i][1] = MTetrahedron::faces_tetra(i, 1);
+ refFaceNodeIdx[i][2] = MTetrahedron::faces_tetra(i, 2);
}
// Init All //
init();
- */
}
TetReferenceSpace::~TetReferenceSpace(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 TetReferenceSpace::toLatex(void) const{
- //const size_t nPerm = pTree->getNPermutation();
- stringstream stream;
- //vector<size_t> perm;
- /*
+ const size_t nRefSpace = refSpaceNodeId.size();
+ stringstream stream;
+
stream << "\\documentclass{article}" << endl << endl
<< "\\usepackage{longtable}" << endl
@@ -74,29 +48,32 @@ string TetReferenceSpace::toLatex(void) const{
<< "\\begin{document}" << endl
<< "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
- << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
+ << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl
+ << endl
<< "\\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] (n3) at(1, 1) {$" << perm[3] << "$};" << 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
+ << "\\node[vertex] (n3) at(1, 1) {$" << refSpaceNodeId[s][3] << "$};"
+ << endl
<< endl;
- for(size_t i = 0; i < 6; i++)
+ for(size_t e = 0; e < 6; 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
@@ -105,6 +82,6 @@ string TetReferenceSpace::toLatex(void) const{
stream << "\\end{longtable}" << endl
<< "\\end{document}" << endl;
- */
+
return stream.str();
}
diff --git a/FunctionSpace/TriReferenceSpace.cpp b/FunctionSpace/TriReferenceSpace.cpp
index 2b7acabe43cd6630796d7619774a36f270cbfb9e..e7b9399f6233fca7e7fe6fd624c97650226e5414 100644
--- a/FunctionSpace/TriReferenceSpace.cpp
+++ b/FunctionSpace/TriReferenceSpace.cpp
@@ -34,8 +34,8 @@ TriReferenceSpace::~TriReferenceSpace(void){
}
string TriReferenceSpace::toLatex(void) const{
- const size_t nRefSpace = refSpaceNodeId.size();
- stringstream stream;
+ const size_t nRefSpace = refSpaceNodeId.size();
+ stringstream stream;
stream << "\\documentclass{article}" << endl << endl
@@ -45,7 +45,8 @@ string TriReferenceSpace::toLatex(void) const{
<< "\\begin{document}" << endl
<< "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
- << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
+ << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl
+ << endl
<< "\\begin{longtable}{ccc}" << endl << endl;