From aead963103198be8e42751bd75dd7f5f2e0a946d Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 8 Jan 2013 16:45:35 +0000
Subject: [PATCH] TetEdgeBasis: Seems OK now -- YEAAAAHHH \o/ -- Functions were
good but with bad number
---
FunctionSpace/TetEdgeBasis.cpp | 216 +++++++++++++++++++++++----------
1 file changed, 149 insertions(+), 67 deletions(-)
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 053641bc47..176ddef8d1 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -62,7 +62,7 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
for(unsigned int s = 0; s < nRefSpace; s++)
(*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
-
+
// Edge Based (Nedelec) //
for(unsigned int s = 0; s < nRefSpace; s++){
for(int e = 0; e < 6; e++){
@@ -102,97 +102,152 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
}
}
}
+
+ // Face Based (Preliminary) //
+ unsigned int faceOffset = 0;
+
+ // Sum
+ Polynomial** sum;
+ sum = new Polynomial*[nRefSpace];
+
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ sum[s] = new Polynomial[4];
+
+ for(unsigned int f = 0; f < 4; f++)
+ sum[s][f] =
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]] +
+ lagrange[(*(*faceV[s])[f])[2]];
+ }
+
+ // Polynomial u
+ Polynomial*** u;
+ u = new Polynomial**[nRefSpace];
+
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ u[s] = new Polynomial*[orderPlus];
+
+ for(int l = 0; l < orderPlus; l++){
+ u[s][l] = new Polynomial[4];
+
+ for(unsigned int f = 0; f < 4; f++)
+ u[s][l][f] =
+ intLegendre[l].compose(lagrange[(*(*faceV[s])[f])[0]] -
+ lagrange[(*(*faceV[s])[f])[1]]
+ ,
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]]);
+ }
+ }
+
+ // Polynomial v
+ Polynomial*** v;
+ v = new Polynomial**[nRefSpace];
+
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ v[s] = new Polynomial*[orderMinus];
- // Face Based //
+ for(int l = 0; l < orderMinus; l++){
+ v[s][l] = new Polynomial[4];
+
+ for(unsigned int f = 0; f < 4; f++)
+ v[s][l][f] =
+ lagrange[(*(*faceV[s])[f])[2]] *
+ sclLegendre[l].compose
+ (lagrange[(*(*faceV[s])[f])[2]] * 2 - sum[s][f], sum[s][f]);
+ }
+ }
+
+ // Face Based (Type 1) //
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge;
for(unsigned int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
for(int f = 0; f < 4; f++){
- // Preliminary Type 1
- 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])[1]]
- ,
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]]);
- Polynomial v =
- lagrange[(*(*faceV[s])[f])[2]] *
- sclLegendre[l2].compose(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum);
+ (*(*basis)[s])[i] =
+ new vector<Polynomial>((u[s][l1][f] * v[s][l2][f]).gradient());
- // Preliminary Type 2
- vector<Polynomial> gradU = u.gradient();
- vector<Polynomial> gradV = v.gradient();
+ i++;
+
+ if(s == 0)
+ faceOffset++;
+ }
+ }
+ }
+ }
+
+ // Face Based (Type 2) //
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nEdge + faceOffset;
+
+ for(unsigned int l1 = 1; l1 < order; l1++){
+ for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
+ for(int f = 0; f < 4; f++){
+ vector<Polynomial> gradU = u[s][l1][f].gradient();
+ vector<Polynomial> gradV = v[s][l2][f].gradient();
vector<Polynomial> vGradU(gradU);
- vGradU[0].mul(v);
- vGradU[1].mul(v);
- vGradU[2].mul(v);
+ vGradU[0].mul(v[s][l2][f]);
+ vGradU[1].mul(v[s][l2][f]);
+ vGradU[2].mul(v[s][l2][f]);
vector<Polynomial> uGradV(gradV);
- uGradV[0].mul(u);
- uGradV[1].mul(u);
- uGradV[2].mul(u);
+ uGradV[0].mul(u[s][l1][f]);
+ uGradV[1].mul(u[s][l1][f]);
+ uGradV[2].mul(u[s][l1][f]);
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] =
- 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++;
- }
}
}
}
}
-
+
+ // Face Based (Type 3) //
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nEdge + faceOffset * 2;
+
+ for(int l2 = 0; l2 < orderMinus; l2++){
+ for(int f = 0; f < 4; f++){
+ 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[s][l2][f]);
+ subGradL1L2V[1].mul(v[s][l2][f]);
+ subGradL1L2V[2].mul(v[s][l2][f]);
+
+ (*(*basis)[s])[i] =
+ new vector<Polynomial>(subGradL1L2V);
+
+ i++;
+ }
+ }
+ }
+
// Cell Based //
const Polynomial one(1, 0, 0, 0);
@@ -312,11 +367,38 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
}
}
}
-
- // Free Temporary Sapce //
+
+ // Free Temporary Space //
+ // Legendre
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
+
+ // Sum
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ delete[] sum[s];
+
+ delete[] sum;
+
+ // Polynomial u
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(int l = 0; l < orderPlus; l++)
+ delete[] u[s][l];
+
+ delete[] u[s];
+ }
+
+ delete[] u;
+
+ // Polynomial v
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(int l = 0; l < orderMinus; l++)
+ delete[] v[s][l];
+
+ delete[] v[s];
+ }
+
+ delete[] v;
}
TetEdgeBasis::~TetEdgeBasis(void){
--
GitLab