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

HUGE BUG FIX: TriEdgeBasis FACES were not confrom with TetEdgeBasis FACES!!

parent 6609db52
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FunctionSpace/TetEdgeBasis.cpp
+ 2
− 4
View file @ 6f76afc1
......@@ -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]]);
......
FunctionSpace/TriEdgeBasis.cpp
+ 68
− 114
View file @ 6f76afc1
......@@ -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){
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment