diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 71c9565c09aff1e6776aa2b11c850e78d1e7ee8f..ddcfefd9b9a2f080c8e05fca9cf94bdf7409d2e6 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -106,8 +106,6 @@ TetEdgeBasis::TetEdgeBasis(size_t order){
}
// Face Based //
- // TO CHECK: Are Triangles face matching tets ?
-
for(size_t s = 0; s < nOrientation; s++){
size_t i = nEdge;
@@ -121,8 +119,8 @@ TetEdgeBasis::TetEdgeBasis(size_t order){
lagrange[faceIdx[s][f][2]];
Polynomial u =
- intLegendre[l1].compose(lagrange[faceIdx[s][f][0]] -
- lagrange[faceIdx[s][f][1]]
+ intLegendre[l1].compose(lagrange[faceIdx[s][f][1]] -
+ lagrange[faceIdx[s][f][0]]
,
lagrange[faceIdx[s][f][0]] +
lagrange[faceIdx[s][f][1]]);
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index b3005310f80ba7417fef4197a9b7307aa5055697..516c055bc0719a9fc5b5c8d3858f0497ff4e9882 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -104,137 +104,91 @@ TriEdgeBasis::TriEdgeBasis(size_t order){
// where f = 0, because triangles
// have only ONE face: the face '0'
- // TO CHECK: Are Triangles face matching tets ?
-
- // Alloc Temp
- Polynomial** u = new Polynomial*[nOrientation];
- Polynomial** v = new Polynomial*[nOrientation];
- Polynomial* p = new Polynomial[nOrientation];
- vector<Polynomial>** subGrad = new vector<Polynomial>*[nOrientation];
-
- for(size_t s = 0; s < nOrientation; s++)
- u[s] = new Polynomial[orderPlus];
-
- for(size_t s = 0; s < nOrientation; s++)
- v[s] = new Polynomial[orderPlus];
-
- // Preliminaries
- for(size_t s = 0; s < nOrientation; 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[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[faceIdx[s][0][2]]);
- }
-
- // Differences of grad(u) and grad(v)
- 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[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[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]);
- (*subGrad[s])[1].sub(l1GradL2[1]);
- (*subGrad[s])[2].sub(l1GradL2[2]);
- }
-
- // Face Basis
for(size_t s = 0; s < nOrientation; s++){
size_t i = nEdge;
- // Type 1
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();
-
- tmp1[0].mul(u[s][l1]);
- tmp1[1].mul(u[s][l1]);
- tmp1[2].mul(u[s][l1]);
-
- tmp2[0].mul(v[s][l2]);
- tmp2[1].mul(v[s][l2]);
- tmp2[2].mul(v[s][l2]);
-
- tmp2[0].add(tmp1[0]);
- tmp2[1].add(tmp1[1]);
- tmp2[2].add(tmp1[2]);
-
- basis[s][i] = new vector<Polynomial>(tmp2);
+ // Preliminary Type 1
+ Polynomial u =
+ intLegendre[l1].compose(lagrange[faceIdx[s][0][1]] -
+ lagrange[faceIdx[s][0][0]]
+ ,
+ lagrange[faceIdx[s][0][0]] +
+ lagrange[faceIdx[s][0][1]]);
+ Polynomial v =
+ lagrange[faceIdx[s][0][2]] *
+ legendre[l2].compose(lagrange[faceIdx[s][0][2]] * 2);
+
+ // 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[faceIdx[s][0][0]].gradient();
+ vector<Polynomial> gradL2 = lagrange[faceIdx[s][0][1]].gradient();
+
+ vector<Polynomial> l2GradL1(gradL1);
+ 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[faceIdx[s][0][0]]);
+ l1GradL2[1].mul(lagrange[faceIdx[s][0][0]]);
+ l1GradL2[2].mul(lagrange[faceIdx[s][0][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] =
+ new vector<Polynomial>((u * v).gradient());
i++;
- }
- }
-
- // Type 2
- 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();
-
- tmp1[0].mul(u[s][l1]);
- tmp1[1].mul(u[s][l1]);
- tmp1[2].mul(u[s][l1]);
- tmp2[0].mul(v[s][l2]);
- tmp2[1].mul(v[s][l2]);
- tmp2[2].mul(v[s][l2]);
-
- tmp2[0].sub(tmp1[0]);
- tmp2[1].sub(tmp1[1]);
- tmp2[2].sub(tmp1[2]);
-
- basis[s][i] = new vector<Polynomial>(tmp2);
+ // Type 2
+ basis[s][i] =
+ new vector<Polynomial>(subGradUV);
i++;
- }
- }
-
- // Type 3
- for(int l = 0; l < orderMinus; l++){
- vector<Polynomial> subGradL1L2V(*subGrad[s]);
- subGradL1L2V[0].mul(v[s][l]);
- subGradL1L2V[1].mul(v[s][l]);
- subGradL1L2V[2].mul(v[s][l]);
-
- basis[s][i] = new vector<Polynomial>(subGradL1L2V);
+ // Type 3
+ if(l1 == 1){
+ basis[s][i] =
+ new vector<Polynomial>(subGradL1L2V);
- i++;
+ i++;
+ }
+ }
}
}
- // Free Temporary Sapce //
- for(size_t s = 0; s < nOrientation; s++)
- delete[] u[s];
-
- for(size_t s = 0; s < nOrientation; s++)
- delete[] v[s];
-
- for(size_t s = 0; s < nOrientation; s++)
- delete subGrad[s];
-
+ // Clear //
delete[] legendre;
delete[] intLegendre;
- delete[] p;
- delete[] u;
- delete[] v;
- delete[] subGrad;
}
TriEdgeBasis::~TriEdgeBasis(void){