Skip to content
Snippets Groups Projects
Commit 95eaa84c authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

Tet Edge + Test Codes

parent 8fc764be
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
FunctionSpace/BasisGenerator.cpp
+ 2
− 1
View file @ 95eaa84c
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include "TriNedelecBasis.h" #include "TriNedelecBasis.h"
#include "TetNodeBasis.h" #include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
#include "HexNodeBasis.h" #include "HexNodeBasis.h"
#include "HexEdgeBasis.h" #include "HexEdgeBasis.h"
...@@ -66,7 +67,7 @@ Basis* BasisGenerator::tetGen(int basisType, ...@@ -66,7 +67,7 @@ Basis* BasisGenerator::tetGen(int basisType,
int order){ int order){
switch(basisType){ switch(basisType){
case 0: return new TetNodeBasis(order); case 0: return new TetNodeBasis(order);
case 1: throw Exception("1-form not implemented on Tetrahedrons"); case 1: return new TetEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Tetrahedrons"); case 2: throw Exception("2-form not implemented on Tetrahedrons");
case 3: throw Exception("3-form not implemented on Tetrahedrons"); case 3: throw Exception("3-form not implemented on Tetrahedrons");
......
FunctionSpace/TetEdgeBasis.cpp
+ 129
− 74
View file @ 95eaa84c
...@@ -73,7 +73,10 @@ TetEdgeBasis::TetEdgeBasis(int order){ ...@@ -73,7 +73,10 @@ TetEdgeBasis::TetEdgeBasis(int order){
// Lagrange // // Lagrange //
lagrange[0] = lagrange[0] =
Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0); Polynomial(1, 0, 0, 0) -
Polynomial(1, 1, 0, 0) -
Polynomial(1, 0, 1, 0) -
Polynomial(1, 0, 0, 1);
lagrange[1] = lagrange[1] =
Polynomial(1, 1, 0, 0); Polynomial(1, 1, 0, 0);
...@@ -154,7 +157,7 @@ TetEdgeBasis::TetEdgeBasis(int order){ ...@@ -154,7 +157,7 @@ TetEdgeBasis::TetEdgeBasis(int order){
vGradU[1].mul(v); vGradU[1].mul(v);
vGradU[2].mul(v); vGradU[2].mul(v);
vector<Polynomial> uGradV(gradU); vector<Polynomial> uGradV(gradV);
uGradV[0].mul(u); uGradV[0].mul(u);
uGradV[1].mul(u); uGradV[1].mul(u);
uGradV[2].mul(u); uGradV[2].mul(u);
...@@ -169,14 +172,14 @@ TetEdgeBasis::TetEdgeBasis(int order){ ...@@ -169,14 +172,14 @@ TetEdgeBasis::TetEdgeBasis(int order){
vector<Polynomial> gradL2 = lagrange[faceV[c][f][1]].gradient(); vector<Polynomial> gradL2 = lagrange[faceV[c][f][1]].gradient();
vector<Polynomial> l2GradL1(gradL1); vector<Polynomial> l2GradL1(gradL1);
l2GradL1[0].mul(l2); l2GradL1[0].mul(lagrange[faceV[c][f][1]]);
l2GradL1[1].mul(l2); l2GradL1[1].mul(lagrange[faceV[c][f][1]]);
l2GradL1[2].mul(l2); l2GradL1[2].mul(lagrange[faceV[c][f][1]]);
vector<Polynomial> l1GradL2(gradL2); vector<Polynomial> l1GradL2(gradL2);
l1GradL2[0].mul(l1); l1GradL2[0].mul(lagrange[faceV[c][f][0]]);
l1GradL2[1].mul(l1); l1GradL2[1].mul(lagrange[faceV[c][f][0]]);
l1GradL2[2].mul(l1); l1GradL2[2].mul(lagrange[faceV[c][f][0]]);
vector<Polynomial> subGradL1L2V(l2GradL1); vector<Polynomial> subGradL1L2V(l2GradL1);
subGradL1L2V[0].sub(l1GradL2[0]); subGradL1L2V[0].sub(l1GradL2[0]);
...@@ -201,6 +204,7 @@ TetEdgeBasis::TetEdgeBasis(int order){ ...@@ -201,6 +204,7 @@ TetEdgeBasis::TetEdgeBasis(int order){
i++; i++;
// Type 3 // Type 3
if(l1 == 1){
(*(*face)[c])[i] = (*(*face)[c])[i] =
new vector<Polynomial>(subGradL1L2V); new vector<Polynomial>(subGradL1L2V);
...@@ -209,77 +213,128 @@ TetEdgeBasis::TetEdgeBasis(int order){ ...@@ -209,77 +213,128 @@ TetEdgeBasis::TetEdgeBasis(int order){
} }
} }
} }
/*
// Cell Based (Preliminary) //
Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
for(int l = 0; l < orderPlus; l++){
u[l] = intLegendre[l].compose(lagrangeSub[0] * -1, lagrangeSum[0]);
v[l] = legendre[l].compose(p);
v[l].mul(lagrange[2]);
} }
// Cell Based (Type 1) //
for(int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
std::vector<Polynomial> tmp = v[l2].gradient();
tmp[0].mul(u[l1]);
tmp[1].mul(u[l1]);
tmp[2].mul(u[l1]);
basis[i] = new std::vector<Polynomial>(u[l1].gradient()); // Cell Based //
Polynomial one(1, 0, 0, 0);
basis[i]->at(0).mul(v[l2]); unsigned int i = 0;
basis[i]->at(1).mul(v[l2]);
basis[i]->at(2).mul(v[l2]);
basis[i]->at(0).add(tmp[0]); for(int l1 = 1; l1 < orderMinus; l1++){
basis[i]->at(1).add(tmp[1]); for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
basis[i]->at(2).add(tmp[2]); 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]);
revBasis[i] = basis[i]; Polynomial v = lagrange[2] * sclLegendre[l2].compose
(lagrange[2] * 2 - (one - lagrange[3]), one - lagrange[3]);
i++; Polynomial w = lagrange[3] * legendre[l3].compose
} (lagrange[3] * 2 - one);
}
// Cell Based (Type 2) // // Preliminary Type 2
for(int l1 = 1; l1 < order; l1++){ vector<Polynomial> gradU = u.gradient();
for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){ vector<Polynomial> gradV = v.gradient();
std::vector<Polynomial> tmp = v[l2].gradient(); vector<Polynomial> gradW = w.gradient();
tmp[0].mul(u[l1]);
tmp[1].mul(u[l1]);
tmp[2].mul(u[l1]);
basis[i] = new std::vector<Polynomial>(u[l1].gradient()); vector<Polynomial> vwGradU(gradU);
vwGradU[0].mul(v);
vwGradU[1].mul(v);
vwGradU[2].mul(v);
basis[i]->at(0).mul(v[l2]); vwGradU[0].mul(w);
basis[i]->at(1).mul(v[l2]); vwGradU[1].mul(w);
basis[i]->at(2).mul(v[l2]); vwGradU[2].mul(w);
basis[i]->at(0).sub(tmp[0]); vector<Polynomial> uwGradV(gradV);
basis[i]->at(1).sub(tmp[1]); uwGradV[0].mul(u);
basis[i]->at(2).sub(tmp[2]); uwGradV[1].mul(u);
uwGradV[2].mul(u);
revBasis[i] = basis[i]; uwGradV[0].mul(w);
uwGradV[1].mul(w);
uwGradV[2].mul(w);
i++; 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]);
// Cell Based (Type 3) // term1[0].add(uvGradW[0]);
for(int l = 0; l < orderMinus; l++){ term1[1].add(uvGradW[1]);
basis[i] = new std::vector<Polynomial>(*basis[0]); term1[2].add(uvGradW[2]);
basis[i]->at(0).mul(v[l]); vector<Polynomial> term2(vwGradU);
basis[i]->at(1).mul(v[l]); term2[0].add(uwGradV[0]);
basis[i]->at(2).mul(v[l]); term2[1].add(uwGradV[1]);
term2[2].add(uwGradV[2]);
revBasis[i] = basis[i]; 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
(*cell)[i] = new vector<Polynomial>((u * v * w).gradient());
i++;
// Type 2 -- Part 1
(*cell)[i] = new vector<Polynomial>(term1);
i++;
// Type 2 -- Part 2
(*cell)[i] = new vector<Polynomial>(term2);
i++;
// Type 3
if(l1 == 1){
(*cell)[i] = new vector<Polynomial>(subGradL1L2VW);
i++; i++;
} }
*/ }
}
}
// Free Temporary Sapce // // Free Temporary Sapce //
delete[] legendre; delete[] legendre;
delete[] sclLegendre; delete[] sclLegendre;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment