diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index f7175273e8e7e8c8ca02b0539bcbb4d9b59ec32d..90eda2048f1742c42811cc0cb3f8243c9da78ccf 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -33,9 +33,9 @@ BasisGenerator::~BasisGenerator(void){
}
BasisLocal* BasisGenerator::generate(unsigned int elementType,
- unsigned int basisType,
- unsigned int order,
- std::string family){
+ unsigned int basisType,
+ unsigned int order,
+ std::string family){
if(!family.compare(std::string("hierarchical")))
return generateHierarchical(elementType, basisType, order);
@@ -48,8 +48,8 @@ BasisLocal* BasisGenerator::generate(unsigned int elementType,
}
BasisLocal* BasisGenerator::generateHierarchical(unsigned int elementType,
- unsigned int basisType,
- unsigned int order){
+ unsigned int basisType,
+ unsigned int order){
switch(elementType){
case TYPE_LIN: return linHierarchicalGen(basisType, order);
case TYPE_TRI: return triHierarchicalGen(basisType, order);
@@ -58,17 +58,17 @@ BasisLocal* BasisGenerator::generateHierarchical(unsigned int elementType,
case TYPE_HEX: return hexHierarchicalGen(basisType, order);
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
- elementType);
+ elementType);
}
}
BasisLocal* BasisGenerator::generateLagrange(unsigned int elementType,
- unsigned int basisType,
- unsigned int order){
+ unsigned int basisType,
+ unsigned int order){
if(basisType != 0)
throw
Exception("Cannot Have a %d-Form Lagrange Basis (0-Form only)",
- basisType);
+ basisType);
switch(elementType){
case TYPE_LIN: return new LineLagrangeBasis(order);
@@ -78,12 +78,12 @@ BasisLocal* BasisGenerator::generateLagrange(unsigned int elementType,
case TYPE_HEX: throw Exception("Lagrange Basis on Hexs not Implemented");
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
- elementType);
+ elementType);
}
}
BasisLocal* BasisGenerator::linHierarchicalGen(unsigned int basisType,
- unsigned int order){
+ unsigned int order){
switch(basisType){
case 0: return new LineNodeBasis(order);
case 1:
@@ -98,7 +98,7 @@ BasisLocal* BasisGenerator::linHierarchicalGen(unsigned int basisType,
}
BasisLocal* BasisGenerator::triHierarchicalGen(unsigned int basisType,
- unsigned int order){
+ unsigned int order){
switch(basisType){
case 0: return new TriNodeBasis(order);
case 1:
@@ -113,7 +113,7 @@ BasisLocal* BasisGenerator::triHierarchicalGen(unsigned int basisType,
}
BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType,
- unsigned int order){
+ unsigned int order){
switch(basisType){
case 0: return new QuadNodeBasis(order);
case 1:
@@ -128,7 +128,7 @@ BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType,
}
BasisLocal* BasisGenerator::tetHierarchicalGen(unsigned int basisType,
- unsigned int order){
+ unsigned int order){
switch(basisType){
case 0: return new TetNodeBasis(order);
case 1:
@@ -143,10 +143,12 @@ BasisLocal* BasisGenerator::tetHierarchicalGen(unsigned int basisType,
}
BasisLocal* BasisGenerator::hexHierarchicalGen(unsigned int basisType,
- unsigned int order){
+ unsigned int order){
switch(basisType){
//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/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 1c294d6a85add94ad1abb0b82f0298075c32beed..b1e9e1eaa1f35accf8ad458379bc568f4a671819 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -11,38 +11,9 @@ FunctionSpaceScalar::~FunctionSpaceScalar(void){
}
double FunctionSpaceScalar::
-interpolate(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& xyz) const{
-
- // Get ABC Space coordinate //
- double abc[3];
- (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element, xyz, abc);
-
- // Get Basis Functions //
- const unsigned int nFun = (*basis)[0]->getNFunction();
- fullMatrix<double> fun(nFun, 1);
-
- (*basis)[0]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
-
- // Interpolate (in Reference Place) //
- double val = 0;
-
- for(unsigned int i = 0; i < nFun; i++)
- val += fun(i, 0) * coef[i];
-
- // Return Interpolated Value //
- return val;
-}
-
-double FunctionSpaceScalar::
-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);
+interpolateInABC(const MElement& element,
+ const std::vector<double>& coef,
+ double abc[3]) const{
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 56350ae3635a70707811b04f7894e637fb9d6645..33bb161870d91120cef2561ee7a52b276d9fa80d 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -30,6 +30,11 @@ class FunctionSpaceScalar : public FunctionSpace{
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const;
+
+ private:
+ double interpolateInABC(const MElement& element,
+ const std::vector<double>& coef,
+ double abc[3]) const;
};
@@ -84,4 +89,37 @@ class FunctionSpaceScalar : public FunctionSpace{
---> check
*/
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline double FunctionSpaceScalar::
+interpolate(const MElement& element,
+ const std::vector<double>& coef,
+ const fullVector<double>& xyz) const{
+
+ // Get ABC Space coordinate //
+ double abc[3];
+ (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
+ xyz(0), xyz(1), xyz(2),
+ abc);
+ // Interpolate in ABC //
+ return interpolateInABC(element, coef, abc);
+}
+
+inline double FunctionSpaceScalar::
+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(0), uvw(1), uvw(2),
+ abc);
+ // Interpolate in ABC //
+ return interpolateInABC(element, coef, abc);
+}
+
#endif
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index cc7fab61a605713e48dab5087e2067fdb242c739..412d206e3bd53469d57f61a0c7634a58845d826b 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -12,68 +12,22 @@ FunctionSpaceVector::~FunctionSpaceVector(void){
}
fullVector<double> FunctionSpaceVector::
-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];
-
- eelement.xyz2uvw(phys, uvw);
+interpolateInABC(const MElement& element,
+ const std::vector<double>& coef,
+ double abc[3]) const{
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
- eelement.getJacobian(uvw[0], uvw[1], uvw[2], invJac);
- invJac.invertInPlace();
-
- // Get Basis Functions //
- const unsigned int nFun = (*basis)[0]->getNFunction();
- fullMatrix<double> fun(nFun, 3);
-
- (*basis)[0]->getFunctions(fun, element, uvw[0], uvw[1], uvw[2]);
-
- // Interpolate (in Reference Place) //
- fullMatrix<double> val(1, 3);
- val(0, 0) = 0;
- val(0, 1) = 0;
- val(0, 2) = 0;
-
- for(unsigned int i = 0; i < nFun; i++){
- val(0, 0) += fun(i, 0) * coef[i];
- val(0, 1) += fun(i, 1) * coef[i];
- val(0, 2) += fun(i, 2) * coef[i];
- }
-
- // Return Interpolated Value //
- fullVector<double> map(3);
- Mapper::hCurl(val, 0, 0, invJac, map);
- return map;
-}
-
-fullVector<double> FunctionSpaceVector::
-interpolateInRefSpace(const MElement& element,
- const std::vector<double>& coef,
- const fullVector<double>& uvw) const{
-
- // Const Cast For MElement //
- MElement& eelement =
- const_cast<MElement&>(element);
-
- // Get Jacobian //
- fullMatrix<double> invJac(3, 3);
- eelement.getJacobian(uvw(0), uvw(1), uvw(2), invJac);
+ (*basis)[0]->getReferenceSpace().getJacobian(element,
+ abc[0], abc[1], abc[2],
+ invJac);
invJac.invertInPlace();
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
fullMatrix<double> fun(nFun, 3);
- (*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) //
fullMatrix<double> val(1, 3);
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 9428e6ffe3a6c9a73a4ffec48a4057f434c811be..6515f0045648f3c960a7648529bf792a5e14075f 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -31,6 +31,12 @@ class FunctionSpaceVector : public FunctionSpace{
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const;
+
+ private:
+ fullVector<double>
+ interpolateInABC(const MElement& element,
+ const std::vector<double>& coef,
+ double abc[3]) const;
};
@@ -85,4 +91,37 @@ class FunctionSpaceVector : public FunctionSpace{
---> check
*/
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline fullVector<double> FunctionSpaceVector::
+interpolate(const MElement& element,
+ const std::vector<double>& coef,
+ const fullVector<double>& xyz) const{
+
+ // Get ABC Space coordinate //
+ double abc[3];
+ (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
+ xyz(0), xyz(1), xyz(2),
+ abc);
+ // Interpolate in ABC //
+ return interpolateInABC(element, coef, abc);
+}
+
+inline fullVector<double> FunctionSpaceVector::
+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(0), uvw(1), uvw(2),
+ abc);
+ // Interpolate in ABC //
+ return interpolateInABC(element, coef, abc);
+}
+
#endif
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index ed14cfcdd821204ebb4178170be5583a59376535..3c98ebec033fd483d5c71fc5de071bb9b31ad90c 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -2,6 +2,7 @@
#include <sstream>
#include "Exception.h"
+#include "Numeric.h"
#include "ReferenceSpace.h"
using namespace std;
@@ -419,69 +420,231 @@ size_t ReferenceSpace::getPermutationIdx(const MElement& element) const{
}
}
+void ReferenceSpace::mapFromABCtoUVW(const MElement& element,
+ double a, double b, double c,
+ double uvw[3]) const{
+ // Get Index Permutation
+ const size_t permutationIdx = getPermutationIdx(element);
+
+ // UVW node coordinate
+ double** uvwNode = new double*[nVertex];
+
+ for(size_t i = 0; i < nVertex; i++)
+ uvwNode[i] = new double[3];
+
+ for(size_t i = 0; i < nVertex; i++)
+ element.getNode(i,
+ uvwNode[UVWtoABCIndex[permutationIdx][i]][0],
+ uvwNode[UVWtoABCIndex[permutationIdx][i]][1],
+ uvwNode[UVWtoABCIndex[permutationIdx][0]][2]);
+
+ // ABC (order 1) grad shape functions
+ double* phiABC = new double[nVertex];
+ element.getShapeFunctions(a, b, c, phiABC, 1);
+
+ // Map From UVW to UVW //
+ uvw[0] = 0;
+ for(size_t i = 0; i < nVertex; i++)
+ uvw[0] += uvwNode[i][0] * phiABC[i];
+
+ uvw[1] = 0;
+ for(size_t i = 0; i < nVertex; i++)
+ uvw[1] += uvwNode[i][1] * phiABC[i];
+
+ uvw[2] = 0;
+ for(size_t i = 0; i < nVertex; i++)
+ uvw[2] += uvwNode[i][2] * phiABC[i];
+
+ // Free //
+ delete[] phiABC;
+
+ for(size_t i = 0; i < nVertex; i++)
+ delete[] uvwNode[i];
+
+ delete[] uvwNode;
+}
+
void ReferenceSpace::mapFromXYZtoABC(const MElement& element,
- const fullVector<double>& xyz,
+ double x, double y, double z,
double abc[3]) const{
// Get UVW coordinate //
- double phys[3] = {xyz(0), xyz(1), xyz(2)};
+ double xyz[3] = {x, y, z};
double uvw[3];
- element.xyz2uvw(phys, uvw);
+ element.xyz2uvw(xyz, uvw);
// Get ABC coordinate //
- fullVector<double> uuvw(3);
- uuvw(0) = uvw[0];
- uuvw(1) = uvw[1];
- uuvw(2) = uvw[2];
-
- mapFromUVWtoABC(element, uuvw, abc);
+ mapFromUVWtoABC(element, uvw[0], uvw[1], uvw[2], abc);
}
void ReferenceSpace::mapFromUVWtoABC(const MElement& element,
- const fullVector<double>& uvw,
+ double u, double v, double w,
double abc[3]) const{
- // Compute coordinate in ABC Space //
- // ABC node coordinate
- double** abcMat = new double*[nVertex];
+ // ABC node coordinate
+ double** abcNode = new double*[nVertex];
for(size_t i = 0; i < nVertex; i++)
- abcMat[i] = new double[3];
+ abcNode[i] = new double[3];
for(size_t i = 0; i < nVertex; i++)
- element.getNode(i, abcMat[i][0], abcMat[i][1], abcMat[i][2]);
+ element.getNode(i, abcNode[i][0], abcNode[i][1], abcNode[i][2]);
// UVW (order 1) shape functions
double* phiUVW = new double[nVertex];
- element.getShapeFunctions(uvw(0), uvw(1), uvw(2), phiUVW, 1);
+ element.getShapeFunctions(u, v, w, phiUVW, 1);
- // Element Permutation Index //
+ // Element Permutation Index
const size_t permutationIdx = getPermutationIdx(element);
- // Map From UVW to ABC //
+ // Map From UVW to ABC
abc[0] = 0;
for(size_t i = 0; i < nVertex; i++)
- abc[0] += abcMat[i][0] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+ abc[0] += abcNode[i][0] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
abc[1] = 0;
for(size_t i = 0; i < nVertex; i++)
- abc[1] += abcMat[i][1] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+ abc[1] += abcNode[i][1] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
abc[2] = 0;
for(size_t i = 0; i < nVertex; i++)
- abc[2] += abcMat[i][2] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+ abc[2] += abcNode[i][2] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
- // Free //
+ // Free
delete[] phiUVW;
for(size_t i = 0; i < nVertex; i++)
- delete[] abcMat[i];
+ delete[] abcNode[i];
+
+ delete[] abcNode;
+}
+
+double ReferenceSpace::getJacobian(const MElement& element,
+ double a, double b, double c,
+ fullMatrix<double>& jac) const{
+ // ABC to UVW Jacobian //
+ // Get Index Permutation
+ const size_t permutationIdx = getPermutationIdx(element);
+
+ // UVW node coordinate
+ double** uvwNode = new double*[nVertex];
+
+ for(size_t i = 0; i < nVertex; i++)
+ uvwNode[i] = new double[3];
+
+ for(size_t i = 0; i < nVertex; i++)
+ element.getNode(i,
+ uvwNode[UVWtoABCIndex[permutationIdx][i]][0],
+ uvwNode[UVWtoABCIndex[permutationIdx][i]][1],
+ uvwNode[UVWtoABCIndex[permutationIdx][i]][2]);
- delete[] abcMat;
+ // ABC (order 1) grad shape functions
+ double phiABC[1256][3]; // Cannot be dynamicaly allocated since
+ // GMSH wants a double[][3] for phiABC !!!
+
+ element.getGradShapeFunctions(a, b, c, phiABC, 1);
+
+ // Jacobian
+ fullMatrix<double> jacABCtoUVW(3, 3);
+ jacABCtoUVW.setAll(0);
+
+ for(size_t i = 0; i < nVertex; i++){
+ jacABCtoUVW(0, 0) += uvwNode[i][0] * phiABC[i][0];
+ jacABCtoUVW(0, 1) += uvwNode[i][1] * phiABC[i][0];
+ jacABCtoUVW(0, 2) += uvwNode[i][2] * phiABC[i][0];
+ }
+
+ for(size_t i = 0; i < nVertex; i++){
+ jacABCtoUVW(1, 0) += uvwNode[i][0] * phiABC[i][1];
+ jacABCtoUVW(1, 1) += uvwNode[i][1] * phiABC[i][1];
+ jacABCtoUVW(1, 2) += uvwNode[i][2] * phiABC[i][1];
+ }
+
+ for(size_t i = 0; i < nVertex; i++){
+ jacABCtoUVW(2, 0) += uvwNode[i][0] * phiABC[i][2];
+ jacABCtoUVW(2, 1) += uvwNode[i][1] * phiABC[i][2];
+ jacABCtoUVW(2, 2) += uvwNode[i][2] * phiABC[i][2];
+ }
+
+ // Regularize Jacobian
+ regularize(element.getDim(), jacABCtoUVW);
+
+ // Free
+ for(size_t i = 0; i < nVertex; i++)
+ delete[] uvwNode[i];
+
+ delete[] uvwNode;
+
+ // Map ABC coordinate to UVW point //
+ double uvw[3];
+ mapFromABCtoUVW(element, a, b, c, uvw);
+
+ // UVW to XYZ Jacobian + Determinant //
+ // NB: Volume is not modified when we go from ABC to UVW //
+ // --> Determinant unchanged //
+ fullMatrix<double> jacUVWtoXYZ(3, 3);
+ double det = element.getJacobian(uvw[0], uvw[1], uvw[2], jacUVWtoXYZ);
+
+ // Product of the two Jacobians & Return //
+ jac.gemm(jacABCtoUVW, jacUVWtoXYZ, 1, 0);
+ return det;
}
-void ReferenceSpace::getJacobian(const MElement& element,
- const fullVector<double>& xyz,
- fullMatrix<double>& jac) const{
+void ReferenceSpace::regularize(size_t dim, fullMatrix<double>& jac){
+ double a[3];
+ double b[3];
+ double c[3];
+
+ switch(dim){
+ case 0:
+ jac(0, 0) = jac(1, 1) = jac(2, 2) = 1.0;
+ jac(0, 1) = jac(1, 0) = jac(2, 0) = 0.0;
+ jac(0, 2) = jac(1, 2) = jac(2, 1) = 0.0;
+ break;
+
+ case 1:
+ a[0] = jac(0, 0);
+ a[1] = jac(0, 1);
+ a[2] = jac(0, 2);
+
+ if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
+ (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))){
+ b[0] = a[1];
+ b[1] = -a[0];
+ b[2] = 0.;
+ }
+
+ else{
+ b[0] = 0.;
+ b[1] = a[2];
+ b[2] = -a[1];
+ }
+
+ norme(b);
+ prodve(a, b, c);
+ norme(c);
+
+ jac(1, 0) = b[0]; jac(1, 1) = b[1]; jac(1, 2) = b[2];
+ jac(2, 0) = c[0]; jac(2, 1) = c[1]; jac(2, 2) = c[2];
+ break;
+
+ case 2:
+ a[0] = jac(0, 0);
+ a[1] = jac(0, 1);
+ a[2] = jac(0, 2);
+
+ b[0] = jac(1, 0);
+ b[1] = jac(1, 1);
+ b[2] = jac(1, 2);
+
+ prodve(a, b, c);
+ norme(c);
+
+ jac(2, 0) = c[0]; jac(2, 1) = c[1]; jac(2, 2) = c[2];
+ break;
+
+ case 3:
+ break;
+ }
}
string ReferenceSpace::toString(void) const{
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index e951efb3b145fa902faef160301d16d908f30787..ba92877fb98b6fcc3b16ac1a92c819ffdcb39c37 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -64,17 +64,21 @@ class ReferenceSpace{
const std::vector<size_t>&
getNodeIndexFromABCtoUVW(const MElement& element) const;
+ void mapFromABCtoUVW(const MElement& element,
+ double a, double b, double c,
+ double uvw[3]) const;
+
void mapFromUVWtoABC(const MElement& element,
- const fullVector<double>& xyz,
+ double u, double v, double w,
double abc[3]) const;
void mapFromXYZtoABC(const MElement& element,
- const fullVector<double>& xyz,
+ double x, double y, double z,
double abc[3]) const;
- void getJacobian(const MElement& element,
- const fullVector<double>& xyz,
- fullMatrix<double>& jac) const;
+ double getJacobian(const MElement& element,
+ double a, double b, double c,
+ fullMatrix<double>& jac) const;
virtual std::string toString(void) const;
virtual std::string toLatex(void) const;
@@ -123,6 +127,8 @@ class ReferenceSpace{
static bool sortPredicate(const std::pair<size_t, size_t>& a,
const std::pair<size_t, size_t>& b);
+
+ static void regularize(size_t dim, fullMatrix<double>& jac);
};
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index e39649f44e8ff7e9777e2d0a24a1bd125f633fc2..c1faa4c2493231bf219d07bf0b6544a6d75de0e2 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -4,17 +4,16 @@
using namespace std;
-TriEdgeBasis::TriEdgeBasis(unsigned int order){
- /*
+TriEdgeBasis::TriEdgeBasis(size_t order){
// Reference Space //
refSpace = new TriReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
- const vector<const vector<const vector<unsigned int>*>*>&
- edgeV = refSpace->getAllEdge();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
- const vector<const vector<const vector<unsigned int>*>*>&
- faceV = refSpace->getAllFace();
+ const vector<vector<vector<size_t> > >&
+ faceIdx = refSpace->getFaceNodeIndex();
// Set Basis Type //
this->order = order;
@@ -54,28 +53,27 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = 0;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = 0;
for(int e = 0; e < 3; e++){
for(int l = 0; l < orderPlus; l++){
// Nedelec
if(l == 0){
- vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
- vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
+ vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[edgeIdx[s][e][0]].gradient();
- tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[0].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[1].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[2].mul(lagrange[edgeIdx[s][e][0]]);
-
- tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[0].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[1].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[2].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[0].sub(tmp1[0]);
tmp2[1].sub(tmp1[1]);
@@ -88,11 +86,11 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
else{
basis[s][i] =
new vector<Polynomial>
- ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
+ ((intLegendre[l].compose(lagrange[edgeIdx[s][e][0]] -
+ lagrange[edgeIdx[s][e][1]]
,
- lagrange[(*(*edgeV[s])[e])[1]] +
- lagrange[(*(*edgeV[s])[e])[0]])).gradient());
+ lagrange[edgeIdx[s][e][1]] +
+ lagrange[edgeIdx[s][e][0]])).gradient());
}
i++;
}
@@ -101,7 +99,7 @@ TriEdgeBasis::TriEdgeBasis(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'
@@ -113,40 +111,40 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
Polynomial* p = new Polynomial[nRefSpace];
vector<Polynomial>** subGrad = new vector<Polynomial>*[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
u[s] = new Polynomial[orderPlus];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
v[s] = new Polynomial[orderPlus];
// Preliminaries
- for(unsigned int s = 0; s < nRefSpace; s++){
- p[s] = lagrange[(*(*faceV[s])[0])[2]] * 2 - Polynomial(1, 0, 0, 0);
+ for(size_t s = 0; s < nRefSpace; s++){
+ p[s] = lagrange[faceIdx[s][0][2]] * 2 - Polynomial(1, 0, 0, 0);
// Polynomial u(x) & v(x)
for(int l = 0; l < orderPlus; l++){
- u[s][l] = intLegendre[l].compose(lagrange[(*(*faceV[s])[0])[1]] -
- lagrange[(*(*faceV[s])[0])[0]],
- lagrange[(*(*faceV[s])[0])[0]] +
- lagrange[(*(*faceV[s])[0])[1]]);
+ u[s][l] = intLegendre[l].compose(lagrange[faceIdx[s][0][1]] -
+ lagrange[faceIdx[s][0][0]],
+ lagrange[faceIdx[s][0][0]] +
+ lagrange[faceIdx[s][0][1]]);
v[s][l] = legendre[l].compose(p[s]);
- v[s][l].mul(lagrange[(*(*faceV[s])[0])[2]]);
+ v[s][l].mul(lagrange[faceIdx[s][0][2]]);
}
// Differences of grad(u) and grad(v)
- vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[0])[0]].gradient();
- vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[0])[1]].gradient();
+ vector<Polynomial> gradL1 = lagrange[faceIdx[s][0][0]].gradient();
+ vector<Polynomial> gradL2 = lagrange[faceIdx[s][0][1]].gradient();
vector<Polynomial> l2GradL1(gradL1);
- l2GradL1[0].mul(lagrange[(*(*faceV[s])[0])[1]]);
- l2GradL1[1].mul(lagrange[(*(*faceV[s])[0])[1]]);
- l2GradL1[2].mul(lagrange[(*(*faceV[s])[0])[1]]);
+ l2GradL1[0].mul(lagrange[faceIdx[s][0][1]]);
+ l2GradL1[1].mul(lagrange[faceIdx[s][0][1]]);
+ l2GradL1[2].mul(lagrange[faceIdx[s][0][1]]);
vector<Polynomial> l1GradL2(gradL2);
- l1GradL2[0].mul(lagrange[(*(*faceV[s])[0])[0]]);
- l1GradL2[1].mul(lagrange[(*(*faceV[s])[0])[0]]);
- l1GradL2[2].mul(lagrange[(*(*faceV[s])[0])[0]]);
+ l1GradL2[0].mul(lagrange[faceIdx[s][0][0]]);
+ l1GradL2[1].mul(lagrange[faceIdx[s][0][0]]);
+ l1GradL2[2].mul(lagrange[faceIdx[s][0][0]]);
subGrad[s] = new vector<Polynomial>(3);
(*subGrad[s])[0].sub(l1GradL2[0]);
@@ -155,11 +153,11 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
}
// Face Basis
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nEdge;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nEdge;
// Type 1
- for(unsigned int l1 = 1; l1 < order; l1++){
+ for(size_t l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[s][l2].gradient();
vector<Polynomial> tmp2 = u[s][l1].gradient();
@@ -183,7 +181,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
}
// Type 2
- for(unsigned int l1 = 1; l1 < order; l1++){
+ for(size_t l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[s][l2].gradient();
vector<Polynomial> tmp2 = u[s][l1].gradient();
@@ -221,13 +219,13 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
}
// Free Temporary Sapce //
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
delete[] u[s];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
delete[] v[s];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
delete subGrad[s];
delete[] legendre;
@@ -236,22 +234,19 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
delete[] u;
delete[] v;
delete[] subGrad;
- */
}
TriEdgeBasis::~TriEdgeBasis(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/TriEdgeBasis.h b/FunctionSpace/TriEdgeBasis.h
index 856b6e0dae07b9f3fb62abbc21eab506fc2334de..99c7a4d08a69555120a9263ef7467beed172e4a4 100644
--- a/FunctionSpace/TriEdgeBasis.h
+++ b/FunctionSpace/TriEdgeBasis.h
@@ -20,7 +20,7 @@ class TriEdgeBasis: public BasisHierarchical1From{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Triangles of the given order
- TriEdgeBasis(unsigned int order);
+ TriEdgeBasis(size_t order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index 645aa039089fceb5a051c69f23c89fc0c2f2436c..7cabff1df297b7e04efa942ea393bb86078e8619 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -4,13 +4,12 @@
using namespace std;
TriNedelecBasis::TriNedelecBasis(void){
- /*
// Reference Space //
refSpace = new TriReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
- const vector<const vector<const vector<unsigned int>*>*>&
- edgeV = refSpace->getAllEdge();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
// Set Basis Type //
order = 0;
@@ -39,23 +38,22 @@ TriNedelecBasis::TriNedelecBasis(void){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int e = 0; e < 3; e++){
- vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
- vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
+ for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t e = 0; e < 3; e++){
+ vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[edgeIdx[s][e][0]].gradient();
- tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[0].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[1].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[2].mul(lagrange[edgeIdx[s][e][0]]);
-
- tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[0].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[1].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[2].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[0].sub(tmp1[0]);
tmp2[1].sub(tmp1[1]);
@@ -64,22 +62,19 @@ TriNedelecBasis::TriNedelecBasis(void){
basis[s][e] = new vector<Polynomial>(tmp2);
}
}
- */
}
TriNedelecBasis::~TriNedelecBasis(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;
- */
}