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

Triangle Cell is now a *FACE* with an *ORIENTATION* -- Seems to work

parent 695714d7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 124 additions and 78 deletions
FunctionSpace/LineLagrangeBasis.cpp
+ 4
− 2
View file @ 70e039d8
......@@ -44,8 +44,10 @@ unsigned int LineLagrangeBasis::getTag(unsigned int order){
fullMatrix<double> LineLagrangeBasis::
linePoint(unsigned int order){
fullMatrix<double> line(order + 1, 1);
line(0 ,0) = 0;
fullMatrix<double> line(order + 1, 3);
line(0, 0) = 0;
line(0, 1) = 0;
line(0, 2) = 0;
if(order > 0){
line(0, 0) = -1.;
......
FunctionSpace/LineReferenceSpace.cpp
+ 4
− 5
View file @ 70e039d8
......@@ -9,12 +9,11 @@ LineReferenceSpace::LineReferenceSpace(void){
// Edge Definition //
nEdge = 1;
refEdge = new unsigned int*[nEdge];
refEdge = new unsigned int*[nEdge];
refEdge[0] = new unsigned int[2];
for(unsigned int i = 0; i < nEdge; i++)
refEdge[i] = new unsigned int[2];
refEdge[0][0] = 0; refEdge[0][1] = 1;
refEdge[0][0] = 0;
refEdge[0][1] = 1;
// Face Definition //
nFace = 0;
......
FunctionSpace/TetEdgeBasis.cpp
+ 2
− 1
View file @ 70e039d8
......@@ -105,6 +105,8 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
}
// Face Based //
// TO CHECK: Are Triangles face matching tets ?
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge;
......@@ -335,4 +337,3 @@ TetEdgeBasis::~TetEdgeBasis(void){
delete[] basis;
}
FunctionSpace/TetNodeBasis.cpp
+ 2
− 0
View file @ 70e039d8
......@@ -90,6 +90,8 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
}
// Face Based //
// TO CHECK: Are Triangles face matching tets ?
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex + nEdge;
......
FunctionSpace/TetReferenceSpace.cpp
+ 0
− 11
View file @ 70e039d8
......@@ -89,14 +89,3 @@ string TetReferenceSpace::toLatex(void) const{
return stream.str();
}
/*
#include <iostream>
int main(void){
TetReferenceSpace p;
cout << p.toLatex() << endl;
return 0;
}
*/
FunctionSpace/TriEdgeBasis.cpp
+ 81
− 43
View file @ 70e039d8
......@@ -12,6 +12,9 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
const vector<const vector<const vector<unsigned int>*>*>&
faceV = refSpace->getAllFace();
// Set Basis Type //
this->order = order;
......@@ -30,8 +33,6 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus];
Polynomial* u = new Polynomial[orderPlus];
Polynomial* v = new Polynomial[orderPlus];
// Legendre Polynomial //
Legendre::legendre(legendre, order);
......@@ -99,50 +100,76 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Face Based //
// NB: We use (*(*faceV[s])[f])[]
// 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*[nRefSpace];
Polynomial** v = new Polynomial*[nRefSpace];
Polynomial* p = new Polynomial[nRefSpace];
vector<Polynomial>** subGrad = new vector<Polynomial>*[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
u[s] = new Polynomial[orderPlus];
for(unsigned int s = 0; s < nRefSpace; s++)
v[s] = new Polynomial[orderPlus];
// Preliminaries
const Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
for(unsigned int s = 0; s < nRefSpace; s++){
p[s] = lagrange[(*(*faceV[s])[0])[2]] * 2 - Polynomial(1, 0, 0, 0);
for(int l = 0; l < orderPlus; l++){
u[l] = intLegendre[l].compose(lagrange[1] - lagrange[0],
lagrange[0] + lagrange[1]);
v[l] = legendre[l].compose(p);
v[l].mul(lagrange[2]);
}
// Polynomial u(x) & v(x)
for(int l = 0; l < orderPlus; l++){
u[s][l] = intLegendre[l].compose(lagrange[(*(*faceV[s])[0])[1]] -
lagrange[(*(*faceV[s])[0])[0]],
lagrange[(*(*faceV[s])[0])[0]] +
lagrange[(*(*faceV[s])[0])[1]]);
vector<Polynomial> gradL1 = lagrange[0].gradient();
vector<Polynomial> gradL2 = lagrange[1].gradient();
v[s][l] = legendre[l].compose(p[s]);
v[s][l].mul(lagrange[(*(*faceV[s])[0])[2]]);
}
vector<Polynomial> l2GradL1(gradL1);
l2GradL1[0].mul(lagrange[1]);
l2GradL1[1].mul(lagrange[1]);
l2GradL1[2].mul(lagrange[1]);
// Differences of grad(u) and grad(v)
vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[0])[0]].gradient();
vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[0])[1]].gradient();
vector<Polynomial> l1GradL2(gradL2);
l1GradL2[0].mul(lagrange[0]);
l1GradL2[1].mul(lagrange[0]);
l1GradL2[2].mul(lagrange[0]);
vector<Polynomial> l2GradL1(gradL1);
l2GradL1[0].mul(lagrange[(*(*faceV[s])[0])[1]]);
l2GradL1[1].mul(lagrange[(*(*faceV[s])[0])[1]]);
l2GradL1[2].mul(lagrange[(*(*faceV[s])[0])[1]]);
vector<Polynomial> subGradL1L2(l2GradL1);
subGradL1L2[0].sub(l1GradL2[0]);
subGradL1L2[1].sub(l1GradL2[1]);
subGradL1L2[2].sub(l1GradL2[2]);
vector<Polynomial> l1GradL2(gradL2);
l1GradL2[0].mul(lagrange[(*(*faceV[s])[0])[0]]);
l1GradL2[1].mul(lagrange[(*(*faceV[s])[0])[0]]);
l1GradL2[2].mul(lagrange[(*(*faceV[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(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge;
// Type 1
for(unsigned int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[l2].gradient();
vector<Polynomial> tmp2 = u[l1].gradient();
vector<Polynomial> tmp1 = v[s][l2].gradient();
vector<Polynomial> tmp2 = u[s][l1].gradient();
tmp1[0].mul(u[l1]);
tmp1[1].mul(u[l1]);
tmp1[2].mul(u[l1]);
tmp1[0].mul(u[s][l1]);
tmp1[1].mul(u[s][l1]);
tmp1[2].mul(u[s][l1]);
tmp2[0].mul(v[l2]);
tmp2[1].mul(v[l2]);
tmp2[2].mul(v[l2]);
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]);
......@@ -157,16 +184,16 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Type 2
for(unsigned int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[l2].gradient();
vector<Polynomial> tmp2 = u[l1].gradient();
vector<Polynomial> tmp1 = v[s][l2].gradient();
vector<Polynomial> tmp2 = u[s][l1].gradient();
tmp1[0].mul(u[l1]);
tmp1[1].mul(u[l1]);
tmp1[2].mul(u[l1]);
tmp1[0].mul(u[s][l1]);
tmp1[1].mul(u[s][l1]);
tmp1[2].mul(u[s][l1]);
tmp2[0].mul(v[l2]);
tmp2[1].mul(v[l2]);
tmp2[2].mul(v[l2]);
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]);
......@@ -180,11 +207,11 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Type 3
for(int l = 0; l < orderMinus; l++){
vector<Polynomial> subGradL1L2V(subGradL1L2);
vector<Polynomial> subGradL1L2V(*subGrad[s]);
subGradL1L2V[0].mul(v[l]);
subGradL1L2V[1].mul(v[l]);
subGradL1L2V[2].mul(v[l]);
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);
......@@ -193,10 +220,21 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
}
// Free Temporary Sapce //
for(unsigned int s = 0; s < nRefSpace; s++)
delete[] u[s];
for(unsigned int s = 0; s < nRefSpace; s++)
delete[] v[s];
for(unsigned int s = 0; s < nRefSpace; s++)
delete subGrad[s];
delete[] legendre;
delete[] intLegendre;
delete[] p;
delete[] u;
delete[] v;
delete[] subGrad;
}
TriEdgeBasis::~TriEdgeBasis(void){
......
FunctionSpace/TriNodeBasis.cpp
+ 18
− 3
View file @ 70e039d8
......@@ -12,6 +12,9 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
const vector<const vector<const vector<unsigned int>*>*>&
faceV = refSpace->getAllFace();
// Set BasisTwo Type //
this->order = order;
......@@ -77,6 +80,13 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
}
// Face Based //
// NB: We use (*(*faceV[s])[f])[]
// where f = 0, because triangles
// have only ONE face: the face '0'
// TO CHECK: Are Triangles face matching tets ?
const Polynomial p = (lagrange[2] * 2) - Polynomial(1, 0, 0, 0);
const int orderMinusTwo = order - 2;
......@@ -86,11 +96,16 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
basis[s][i] =
new Polynomial(intLegendre[l1].compose(lagrange[1] - lagrange[0],
lagrange[1] + lagrange[0])
new Polynomial(intLegendre[l1].compose(lagrange[(*(*faceV[s])[0])[1]] -
lagrange[(*(*faceV[s])[0])[0]]
,
lagrange[(*(*faceV[s])[0])[0]] +
lagrange[(*(*faceV[s])[0])[1]])
*
legendre[l2].compose(p) * lagrange[2]);
legendre[l2].compose(lagrange[(*(*faceV[s])[0])[2]])
*
lagrange[(*(*faceV[s])[0])[2]]);
i++;
}
}
......
FunctionSpace/TriReferenceSpace.cpp
+ 13
− 13
View file @ 70e039d8
......@@ -19,8 +19,13 @@ TriReferenceSpace::TriReferenceSpace(void){
}
// Face Definition //
nFace = 0;
refFace = NULL;
nFace = 1;
refFace = new unsigned int*[nFace];
refFace[0] = new unsigned int[3];
refFace[0][0] = 0;
refFace[0][1] = 1;
refFace[0][2] = 2;
// Init All //
init();
......@@ -32,6 +37,12 @@ TriReferenceSpace::~TriReferenceSpace(void){
delete[] refEdge[i];
delete[] refEdge;
// Delete Ref Face //
for(unsigned int i = 0; i < nFace; i++)
delete[] refFace[i];
delete[] refFace;
}
string TriReferenceSpace::toLatex(void) const{
......@@ -75,14 +86,3 @@ string TriReferenceSpace::toLatex(void) const{
return stream.str();
}
/*
#include <iostream>
int main(void){
TriReferenceSpace p;
cout << p.toLatex() << endl;
return 0;
}
*/
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment