diff --git a/FunctionSpace/LineLagrangeBasis.cpp b/FunctionSpace/LineLagrangeBasis.cpp
index c224d580214321e0063f38a01767819638bebd8a..ec2eb76cfcfcc9a152ecd5be8bca6f92a0943b5b 100644
--- a/FunctionSpace/LineLagrangeBasis.cpp
+++ b/FunctionSpace/LineLagrangeBasis.cpp
@@ -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.;
diff --git a/FunctionSpace/LineReferenceSpace.cpp b/FunctionSpace/LineReferenceSpace.cpp
index d4d41069dd02ad999328ea232a3f37a2400363b6..7ba0eb2d4aad1b8c1999650fe137835fdce09971 100644
--- a/FunctionSpace/LineReferenceSpace.cpp
+++ b/FunctionSpace/LineReferenceSpace.cpp
@@ -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;
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index afffb7e2fe982f0b8f3e7fd7ae5a6c468c6d61f5..73795d0341dae4168b0994a698efd5639c47ab39 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -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;
}
-
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 37971f323cde740a0898ed833b8bd1674a104789..27cb84d6ef6c631fc34fdf914332ead9a7430ead 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -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;
diff --git a/FunctionSpace/TetReferenceSpace.cpp b/FunctionSpace/TetReferenceSpace.cpp
index 247172bb05bf39fdf8b7d2f61e7aec9cc9f128d7..585744aa12d98b90f125ec7b5a98b748c46d2b65 100644
--- a/FunctionSpace/TetReferenceSpace.cpp
+++ b/FunctionSpace/TetReferenceSpace.cpp
@@ -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;
-}
-*/
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index f9c1fe7704ee626ab1a03fad95f226971e808e83..088215ee1211d87925a03499249bbcb6835fe874 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -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){
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index dfad8d928d268956994455bd76eebb0081b4050d..3a0cefd789f4ab94f295fe539027684fe1370bbf 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -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++;
}
}
diff --git a/FunctionSpace/TriReferenceSpace.cpp b/FunctionSpace/TriReferenceSpace.cpp
index e4a1d15850f16d43224d50b872964812253f13b5..92720d547cb1c07357c8a383e49f93ea7eff8af5 100644
--- a/FunctionSpace/TriReferenceSpace.cpp
+++ b/FunctionSpace/TriReferenceSpace.cpp
@@ -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;
-}
-*/