diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 31d526f2c2f3a6e4674a9650d8c57f44b6fc469d..7eeb73f0e42055e7bc7962fd0e7e870b40ef7221 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -17,6 +17,7 @@
#include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
#include "TetNedelecBasis.h"
+#include "TetLagrangeBasis.h"
#include "HexNodeBasis.h"
#include "HexEdgeBasis.h"
@@ -70,7 +71,7 @@ BasisLocal* BasisGenerator::generateLagrange(unsigned int elementType,
case TYPE_LIN: throw Exception("Lagrange Basis on Lines not Implemented");
case TYPE_TRI: return new TriLagrangeBasis(order);
case TYPE_QUA: throw Exception("Lagrange Basis on Quads not Implemented");
- case TYPE_TET: throw Exception("Lagrange Basis on Tets not Implemented");
+ case TYPE_TET: return new TetLagrangeBasis(order);
case TYPE_HEX: throw Exception("Lagrange Basis on Hexs not Implemented");
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index a4a0b6894485c71846237dfcda847d331a5a9694..5e052f9b307e18e1fa017a4903c4b91d99341a78 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -38,6 +38,7 @@ set(SRC
TetNodeBasis.cpp
TetEdgeBasis.cpp
TetNedelecBasis.cpp
+ TetLagrangeBasis.cpp
FunctionSpace.cpp
FunctionSpaceScalar.cpp
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index e29bcb5599be0635cb27c930a595305cef3920f8..2b546f000fbda362cb8aa7ff39e98a802066a789 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -17,7 +17,7 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Set Basis Type //
this->order = order;
-
+
type = 1;
dim = 3;
@@ -40,19 +40,19 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
Legendre::legendre(legendre, orderMinusTwo);
Legendre::scaled(sclLegendre, orderMinusTwo);
Legendre::intScaled(intLegendre, orderPlus);
-
+
// Lagrange Polynomial //
- const Polynomial lagrange[4] =
+ const Polynomial lagrange[4] =
{
- Polynomial(Polynomial(1, 0, 0, 0) -
- Polynomial(1, 1, 0, 0) -
- Polynomial(1, 0, 1, 0) -
+ Polynomial(Polynomial(1, 0, 0, 0) -
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0) -
Polynomial(1, 0, 0, 1)),
-
+
Polynomial(Polynomial(1, 1, 0, 0)),
Polynomial(Polynomial(1, 0, 1, 0)),
-
+
Polynomial(Polynomial(1, 0, 0, 1))
};
@@ -62,7 +62,7 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
for(unsigned int 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;
@@ -73,130 +73,130 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
if(l == 0){
vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
vector<Polynomial> tmp2 = lagrange[(*(*edgeV[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]]);
-
+
tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
-
- tmp1[0].sub(tmp2[0]);
- tmp1[1].sub(tmp2[1]);
- tmp1[2].sub(tmp2[2]);
-
+ tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
+
basis[s][i] = new vector<Polynomial>(tmp1);
}
// High Order
else{
- basis[s][i] =
+ basis[s][i] =
new vector<Polynomial>
- ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
+ ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
lagrange[(*(*edgeV[s])[e])[1]]
,
- lagrange[(*(*edgeV[s])[e])[0]] +
+ lagrange[(*(*edgeV[s])[e])[0]] +
lagrange[(*(*edgeV[s])[e])[1]])).gradient());
}
i++;
}
}
}
-
+
// Face Based //
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge;
-
+
for(int f = 0; f < 4; f++){
for(unsigned int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
// Preliminary Type 1
- Polynomial sum =
+ Polynomial sum =
lagrange[(*(*faceV[s])[f])[0]] +
lagrange[(*(*faceV[s])[f])[1]] +
- lagrange[(*(*faceV[s])[f])[2]];
-
- Polynomial u =
- intLegendre[l1].compose(lagrange[(*(*faceV[s])[f])[0]] -
+ lagrange[(*(*faceV[s])[f])[2]];
+
+ Polynomial u =
+ intLegendre[l1].compose(lagrange[(*(*faceV[s])[f])[0]] -
lagrange[(*(*faceV[s])[f])[1]]
,
- lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[0]] +
lagrange[(*(*faceV[s])[f])[1]]);
- Polynomial v =
- lagrange[(*(*faceV[s])[f])[2]] *
+ Polynomial v =
+ lagrange[(*(*faceV[s])[f])[2]] *
sclLegendre[l2].compose(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum);
-
+
// Preliminary Type 2
vector<Polynomial> gradU = u.gradient();
vector<Polynomial> gradV = v.gradient();
-
+
vector<Polynomial> vGradU(gradU);
vGradU[0].mul(v);
vGradU[1].mul(v);
vGradU[2].mul(v);
-
+
vector<Polynomial> uGradV(gradV);
uGradV[0].mul(u);
uGradV[1].mul(u);
uGradV[2].mul(u);
-
+
vector<Polynomial> subGradUV(vGradU);
subGradUV[0].sub(uGradV[0]);
subGradUV[1].sub(uGradV[1]);
subGradUV[2].sub(uGradV[2]);
-
+
// Preliminary Type 3
vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[f])[0]].gradient();
vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[f])[1]].gradient();
-
+
vector<Polynomial> l2GradL1(gradL1);
l2GradL1[0].mul(lagrange[(*(*faceV[s])[f])[1]]);
l2GradL1[1].mul(lagrange[(*(*faceV[s])[f])[1]]);
l2GradL1[2].mul(lagrange[(*(*faceV[s])[f])[1]]);
-
+
vector<Polynomial> l1GradL2(gradL2);
l1GradL2[0].mul(lagrange[(*(*faceV[s])[f])[0]]);
l1GradL2[1].mul(lagrange[(*(*faceV[s])[f])[0]]);
l1GradL2[2].mul(lagrange[(*(*faceV[s])[f])[0]]);
-
+
vector<Polynomial> subGradL1L2V(l2GradL1);
subGradL1L2V[0].sub(l1GradL2[0]);
subGradL1L2V[1].sub(l1GradL2[1]);
subGradL1L2V[2].sub(l1GradL2[2]);
-
+
subGradL1L2V[0].mul(v);
subGradL1L2V[1].mul(v);
subGradL1L2V[2].mul(v);
-
-
+
+
// Type 1
- basis[s][i] =
+ basis[s][i] =
new vector<Polynomial>((u * v).gradient());
-
+
i++;
-
+
// Type 2
basis[s][i] =
new vector<Polynomial>(subGradUV);
-
+
i++;
-
+
// Type 3
if(l1 == 1){
basis[s][i] =
new vector<Polynomial>(subGradL1L2V);
-
+
i++;
}
}
}
}
- }
-
+ }
+
// Cell Based //
const Polynomial one(1, 0, 0, 0);
-
+
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge + nFace;
@@ -205,105 +205,105 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusTwo; l3++){
// Preliminary Type 1
Polynomial u = intLegendre[l1].compose
- (lagrange[0] - lagrange[1],
+ (lagrange[0] - lagrange[1],
lagrange[0] + lagrange[1]);
-
+
Polynomial v = lagrange[2] * sclLegendre[l2].compose
(lagrange[2] * 2 - (one - lagrange[3]), one - lagrange[3]);
-
+
Polynomial w = lagrange[3] * legendre[l3].compose
(lagrange[3] * 2 - one);
-
+
// Preliminary Type 2
vector<Polynomial> gradU = u.gradient();
vector<Polynomial> gradV = v.gradient();
vector<Polynomial> gradW = w.gradient();
-
+
vector<Polynomial> vwGradU(gradU);
vwGradU[0].mul(v);
vwGradU[1].mul(v);
vwGradU[2].mul(v);
-
+
vwGradU[0].mul(w);
vwGradU[1].mul(w);
vwGradU[2].mul(w);
-
+
vector<Polynomial> uwGradV(gradV);
uwGradV[0].mul(u);
uwGradV[1].mul(u);
uwGradV[2].mul(u);
-
+
uwGradV[0].mul(w);
uwGradV[1].mul(w);
uwGradV[2].mul(w);
-
+
vector<Polynomial> uvGradW(gradW);
uvGradW[0].mul(u);
uvGradW[1].mul(u);
uvGradW[2].mul(u);
-
+
uvGradW[0].mul(v);
uvGradW[1].mul(v);
uvGradW[2].mul(v);
-
+
vector<Polynomial> term1(vwGradU);
term1[0].sub(uwGradV[0]);
term1[1].sub(uwGradV[1]);
term1[2].sub(uwGradV[2]);
-
+
term1[0].add(uvGradW[0]);
term1[1].add(uvGradW[1]);
term1[2].add(uvGradW[2]);
-
+
vector<Polynomial> term2(vwGradU);
term2[0].add(uwGradV[0]);
term2[1].add(uwGradV[1]);
term2[2].add(uwGradV[2]);
-
+
term2[0].sub(uvGradW[0]);
term2[1].sub(uvGradW[1]);
term2[2].sub(uvGradW[2]);
-
+
// Preliminary Type 3
vector<Polynomial> gradL1 = lagrange[0].gradient();
vector<Polynomial> gradL2 = lagrange[1].gradient();
-
+
vector<Polynomial> l2GradL1(gradL1);
l2GradL1[0].mul(lagrange[1]);
l2GradL1[1].mul(lagrange[1]);
l2GradL1[2].mul(lagrange[1]);
-
+
vector<Polynomial> l1GradL2(gradL2);
l1GradL2[0].mul(lagrange[0]);
l1GradL2[1].mul(lagrange[0]);
l1GradL2[2].mul(lagrange[0]);
-
+
vector<Polynomial> subGradL1L2VW(l2GradL1);
subGradL1L2VW[0].sub(l1GradL2[0]);
subGradL1L2VW[1].sub(l1GradL2[1]);
subGradL1L2VW[2].sub(l1GradL2[2]);
-
+
subGradL1L2VW[0].mul(v);
subGradL1L2VW[1].mul(v);
subGradL1L2VW[2].mul(v);
-
+
subGradL1L2VW[0].mul(w);
subGradL1L2VW[1].mul(w);
subGradL1L2VW[2].mul(w);
-
-
+
+
// Type 1
basis[s][i] = new vector<Polynomial>((u * v * w).gradient());
i++;
-
+
// Type 2 -- Part 1
basis[s][i] = new vector<Polynomial>(term1);
i++;
-
+
// Type 2 -- Part 2
basis[s][i] = new vector<Polynomial>(term2);
i++;
-
+
// Type 3
if(l1 == 1){
basis[s][i] = new vector<Polynomial>(subGradL1L2VW);
@@ -313,7 +313,7 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
}
}
}
-
+
// Free Temporary Space //
// Legendre
delete[] legendre;
diff --git a/FunctionSpace/TetLagrangeBasis.cpp b/FunctionSpace/TetLagrangeBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4913e3e2c3bf6f58d27cf7a465cd203c68a33a9d
--- /dev/null
+++ b/FunctionSpace/TetLagrangeBasis.cpp
@@ -0,0 +1,372 @@
+#include "Exception.h"
+#include "TetLagrangeBasis.h"
+
+TetLagrangeBasis::TetLagrangeBasis(unsigned int order){
+ // Set Basis Type //
+ this->order = order;
+
+ type = 0;
+ dim = 3;
+
+ nVertex = 4;
+ nEdge = 6 * (order - 1);
+ nFace = 2 * (order - 1) * (order - 2);
+ nCell = (order - 1) * (order - 2) * (order - 3) / 6;
+ nFunction = nVertex + nEdge + nFace + nCell;
+
+ // Init polynomialBasis //
+ lBasis = new polynomialBasis(getTag(order));
+
+ // Init Lagrange Point //
+ lPoint = new fullMatrix<double>(tetPoint(order));
+}
+
+TetLagrangeBasis::~TetLagrangeBasis(void){
+ delete lBasis;
+ delete lPoint;
+}
+
+unsigned int TetLagrangeBasis::getTag(unsigned int order){
+ unsigned int tag = nodalBasis::getTag(TYPE_TET, order, false);
+
+ if(tag)
+ return tag;
+
+ else
+ throw Exception
+ ("Can't instanciate an order %d Lagrangian Basis for a Tetrahedron",
+ order);
+}
+
+fullMatrix<double> TetLagrangeBasis::
+tetPoint(unsigned int order){
+ unsigned int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
+ fullMatrix<double> point(nbPoints, 3);
+
+ double overOrder = (order == 0 ? 1. : 1. / order);
+
+ point(0, 0) = 0.;
+ point(0, 1) = 0.;
+ point(0, 2) = 0.;
+
+ if(order > 0){
+ point(1, 0) = order;
+ point(1, 1) = 0;
+ point(1, 2) = 0;
+
+ point(2, 0) = 0.;
+ point(2, 1) = order;
+ point(2, 2) = 0.;
+
+ point(3, 0) = 0.;
+ point(3, 1) = 0.;
+ point(3, 2) = order;
+
+ if(order > 1){
+ for(unsigned int k = 0; k < (order - 1); k++){
+ point(4 + k, 0) = k + 1;
+ point(4 + order - 1 + k, 0) = order - 1 - k;
+ point(4 + 2 * (order - 1) + k, 0) = 0.;
+ point(4 + 3 * (order - 1) + k, 0) = 0.;
+ point(4 + 4 * (order - 1) + k, 0) = 0.;
+ point(4 + 5 * (order - 1) + k, 0) = k+1;
+
+ point(4 + k, 1) = 0.;
+ point(4 + order - 1 + k, 1) = k + 1;
+ point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+ point(4 + 3 * (order - 1) + k, 1) = 0.;
+ point(4 + 4 * (order - 1) + k, 1) = k+1;
+ point(4 + 5 * (order - 1) + k, 1) = 0.;
+
+ point(4 + k, 2) = 0.;
+ point(4 + order - 1 + k, 2) = 0.;
+ point(4 + 2 * (order - 1) + k, 2) = 0.;
+ point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
+ point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+ point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
+ }
+
+ if(order > 2){
+ unsigned int ns = 4 + 6 * (order - 1);
+ unsigned int nbdofface = nbdoftriangle(order - 3);
+
+ double *u = new double[nbdofface];
+ double *v = new double[nbdofface];
+ double *w = new double[nbdofface];
+
+ nodepositionface0(order - 3, u, v, w);
+
+ // u-v plane
+
+ for(unsigned int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3);
+ }
+
+ ns = ns + nbdofface;
+
+ // u-w plane
+
+ nodepositionface1(order - 3, u, v, w);
+
+ for(unsigned int i=0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) ;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
+ }
+
+ // v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface2(order - 3, u, v, w);
+
+ for(unsigned int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3);
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
+ }
+
+ // u-v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface3(order - 3, u, v, w);
+
+ for(unsigned int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
+ }
+
+ ns = ns + nbdofface;
+
+ delete [] u;
+ delete [] v;
+ delete [] w;
+
+ if(order > 3){
+ fullMatrix<double> interior = tetPoint(order - 4);
+
+ for(int k = 0; k < interior.size1(); k++){
+ point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
+ point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
+ point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
+ }
+ }
+ }
+ }
+ }
+
+ point.scale(overOrder);
+ return point;
+}
+
+unsigned int TetLagrangeBasis::nbdoftriangle(unsigned int order){
+ return (order + 1) * (order + 2) / 2;
+}
+
+void TetLagrangeBasis::
+nodepositionface0(unsigned int order, double *u, double *v, double *w){
+
+ unsigned int ndofT = nbdoftriangle(order);
+
+ if(order == 0){
+ u[0] = 0.;
+ v[0] = 0.;
+ w[0] = 0.;
+ return;
+ }
+
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= order; w[1] = 0.;
+ u[2]= order; v[2]= 0.; w[2] = 0.;
+
+ // edges
+ for(unsigned int k = 0; k < (order - 1); k++){
+ u[3 + k] = 0.;
+ v[3 + k] = k + 1;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = k + 1;
+ v[3 + order - 1 + k] = order - 1 - k ;
+ w[3 + order - 1 + k] = 0.;
+
+ u[3 + 2 * (order - 1) + k] = order - 1 - k;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ w[3 + 2 * (order - 1) + k] = 0.;
+ }
+
+ if(order > 2){
+ unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+ nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3* (order - 1)]);
+
+ for(unsigned int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
+ }
+ }
+
+ for(unsigned int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+void TetLagrangeBasis::nodepositionface1
+(unsigned int order, double *u, double *v, double *w){
+
+ unsigned int ndofT = nbdoftriangle(order);
+
+ if(order == 0){
+ u[0] = 0.;
+ v[0] = 0.;
+ w[0] = 0.;
+ return;
+ }
+
+ u[0] = 0.; v[0]= 0.; w[0] = 0.;
+ u[1] = order; v[1]= 0.; w[1] = 0.;
+ u[2] = 0.; v[2]= 0.; w[2] = order;
+
+ // edges
+ for(unsigned int k = 0; k < (order - 1); k++){
+ u[3 + k] = k + 1;
+ v[3 + k] = 0.;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = order - 1 - k;
+ v[3 + order - 1 + k] = 0.;
+ w[3 + order - 1+ k ] = k + 1;
+
+ u[3 + 2 * (order - 1) + k] = 0. ;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ w[3 + 2 * (order - 1) + k] = order - 1 - k;
+ }
+
+ if(order > 2){
+ unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+ nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
+ &w[3 + 3 * (order - 1)]);
+
+ for(unsigned int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
+ }
+
+ for (unsigned int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+void TetLagrangeBasis::
+nodepositionface2(unsigned int order, double *u, double *v, double *w){
+
+ unsigned int ndofT = nbdoftriangle(order);
+
+ if(order == 0){
+ u[0] = 0.;
+ v[0] = 0.;
+ w[0] = 0.;
+ return;
+ }
+
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= 0.; w[1] = order;
+ u[2]= 0.; v[2]= order; w[2] = 0.;
+
+ // edges
+ for(unsigned int k = 0; k < (order - 1); k++){
+ u[3 + k] = 0.;
+ v[3 + k] = 0.;
+ w[3 + k] = k + 1;
+
+ u[3 + order - 1 + k] = 0.;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = order - 1 - k;
+
+ u[3 + 2 * (order - 1) + k] = 0.;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = 0.;
+ }
+
+ if(order > 2){
+ unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+ nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3 * (order - 1)]);
+
+ for(unsigned int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
+ }
+
+ for(unsigned int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+void TetLagrangeBasis::
+nodepositionface3(unsigned int order, double *u, double *v, double *w){
+
+ unsigned int ndofT = nbdoftriangle(order);
+
+ if(order == 0){
+ u[0] = 0.;
+ v[0] = 0.;
+ w[0] = 0.;
+ return;
+ }
+
+ u[0]= 0.; v[0]= 0.; w[0] = order;
+ u[1]= order; v[1]= 0.; w[1] = 0.;
+ u[2]= 0.; v[2]= order; w[2] = 0.;
+
+ // edges
+ for(unsigned int k = 0; k < (order - 1); k++){
+ u[3 + k] = k + 1;
+ v[3 + k] = 0.;
+ w[3 + k] = order - 1 - k;
+
+ u[3 + order - 1 + k] = order - 1 - k;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = 0.;
+
+ u[3 + 2 * (order - 1) + k] = 0.;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = k + 1;
+ }
+
+ if(order > 2){
+ unsigned int nbdoftemp = nbdoftriangle(order - 3);
+
+ nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3 * (order - 1)]);
+
+ for(unsigned int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
+ }
+
+ for (unsigned int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
diff --git a/FunctionSpace/TetLagrangeBasis.h b/FunctionSpace/TetLagrangeBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..27f7e5a6ecf181025c677e6d020f6d3a9bad7b5a
--- /dev/null
+++ b/FunctionSpace/TetLagrangeBasis.h
@@ -0,0 +1,60 @@
+#ifndef _TETLAGRANGEBASIS_H_
+#define _TETLAGRANGEBASIS_H_
+
+#include "BasisLagrange.h"
+
+/**
+ @class TetLagrangeBasis
+ @brief Lagrange Basis for Tetrahedra
+
+ This class can instantiate a @em Lagrange @em Basis
+ for a Tetrahedron and for a given Order.@n
+
+ It uses
+ <a href="http://geuz.org/gmsh/">gmsh</a> Basis.
+ */
+
+class TetLagrangeBasis: public BasisLagrange{
+ public:
+ //! @param order A natural number
+ //!
+ //! Returns a new TetLagrangeBasis
+ //! of the given Order
+ TetLagrangeBasis(unsigned int order);
+
+ //! Deletes this Basis
+ //!
+ virtual ~TetLagrangeBasis(void);
+
+ private:
+ //! @param order A natural number
+ //! @return Returns the @em tag of a @em Tetrahedron of
+ //! the given order
+ static unsigned int getTag(unsigned int order);
+
+ //! @param order A natural number
+ //! @return Returns Lagrangian Points on a Tetrahedron
+ //! for the given Order
+ static fullMatrix<double> tetPoint(unsigned int order);
+
+ //! Unknown gmsh function
+ static unsigned int nbdoftriangle(unsigned int order);
+
+ //! Unknown gmsh function
+ static void nodepositionface0(unsigned int order,
+ double *u, double *v, double *w);
+
+ //! Unknown gmsh function
+ static void nodepositionface1(unsigned int order,
+ double *u, double *v, double *w);
+
+ //! Unknown gmsh function
+ static void nodepositionface2(unsigned int order,
+ double *u, double *v, double *w);
+
+ //! Unknown gmsh function
+ static void nodepositionface3(unsigned int order,
+ double *u, double *v, double *w);
+};
+
+#endif
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 99ff0469da87fec71debcf68ee1923239fe7e8df..d9258206f8fee9580c6d31c436a727e96486cc0b 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -6,7 +6,7 @@ TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
this->order = order;
type = 0;
- dim = 2;
+ dim = 2;
nVertex = 3;
nEdge = 3 * (order - 1);
@@ -27,23 +27,15 @@ TriLagrangeBasis::~TriLagrangeBasis(void){
}
unsigned int TriLagrangeBasis::getTag(unsigned int order){
- switch(order){
- case 1: return MSH_TRI_3;
- case 2: return MSH_TRI_6;
- case 3: return MSH_TRI_10;
- case 4: return MSH_TRI_15;
- case 5: return MSH_TRI_21;
- case 6: return MSH_TRI_28;
- case 7: return MSH_TRI_36;
- case 8: return MSH_TRI_45;
- case 9: return MSH_TRI_55;
- case 10: return MSH_TRI_66;
-
- default:
+ unsigned int tag = nodalBasis::getTag(TYPE_TRI, order, false);
+
+ if(tag)
+ return tag;
+
+ else
throw Exception
("Can't instanciate an order %d Lagrangian Basis for a Triangle",
order);
- }
}
fullMatrix<double> TriLagrangeBasis::
@@ -57,17 +49,17 @@ triPoint(unsigned int order){
if(order > 0){
double dd = 1. / order;
-
+
point(1, 0) = 1;
point(1, 1) = 0;
point(1, 2) = 0;
-
+
point(2, 0) = 0;
point(2, 1) = 1;
point(2, 2) = 0;
-
+
unsigned int index = 3;
-
+
if(order > 1){
double ksi = 0;
double eta = 0;
@@ -104,7 +96,7 @@ triPoint(unsigned int order){
if(order > 2){
fullMatrix<double> inner = triPoint(order - 3);
-
+
inner.scale(1. - 3. * dd);
inner.add(dd);
point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);