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

Same Canva for all Basis -- Quads OK -- Hex ToBeDone

parent ce437a9f
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 426 additions and 366 deletions
FunctionSpace/HexNodeBasis.cpp
+ 59
− 6
View file @ 58bedb0a
#include "HexNodeBasis.h"
#include "Legendre.h"
using namespace std;
HexNodeBasis::HexNodeBasis(const int order){
// Set Basis Type //
this->order = order;
......@@ -13,7 +15,7 @@ HexNodeBasis::HexNodeBasis(const int order){
nFace = 6 * (order - 1) * (order - 1);
nCell = (order - 1) * (order - 1) * (order - 1);
size = (order + 1) * (order + 1) * (order + 1);
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
......@@ -25,6 +27,27 @@ HexNodeBasis::HexNodeBasis(const int order){
// Legendre Polynomial //
Legendre::integrated(legendre, order);
// Vertices definig Edges & Permutations //
const int edgeV[2][12][2] =
{
{ {0, 1}, {0, 3}, {0, 4}, {1, 2}, {1, 5}, {2, 3},
{2, 6}, {3, 7}, {4, 5}, {4, 7}, {5, 6}, {6, 7} },
{ {1, 0}, {3, 0}, {4, 0}, {2, 1}, {5, 1}, {3, 2},
{6, 2}, {7, 3}, {5, 4}, {7, 4}, {6, 5}, {7, 6} },
};
const int faceV[6][4] =
{
{0, 3, 2, 1},
{0, 1, 5, 4},
{0, 4, 7, 3},
{1, 2, 6, 5},
{2, 3, 7, 6},
{4, 5, 6, 7}
};
// Lifting //
lifting[0] =
......@@ -67,7 +90,7 @@ HexNodeBasis::HexNodeBasis(const int order){
(Polynomial(1, 0, 1, 0)) +
Polynomial(1, 0, 0, 1);
/*
// Basis //
basis = new std::vector<const Polynomial*>(size);
......@@ -194,7 +217,7 @@ HexNodeBasis::HexNodeBasis(const int order){
}
}
*/
// Free Temporary Sapce //
delete[] legendre;
delete[] lifting;
......@@ -205,8 +228,38 @@ HexNodeBasis::HexNodeBasis(const int order){
}
HexNodeBasis::~HexNodeBasis(void){
for(int i = 0; i < size; i++)
delete (*basis)[i];
// Vertex Based //
for(int i = 0; i < nVertex; i++)
delete (*node)[i];
delete node;
// Edge Based //
for(int c = 0; c < 2; c++){
for(int i = 0; i < nEdge; i++)
delete (*(*edge)[c])[i];
delete (*edge)[c];
}
delete edge;
// Face Based //
for(int c = 0; c < 6; c++){
for(int i = 0; i < nFace; i++)
delete (*(*face)[c])[i];
delete (*face)[c];
}
delete face;
// Cell Based //
for(int i = 0; i < nCell; i++)
delete (*cell)[i];
delete basis;
delete cell;
}
FunctionSpace/QuadEdgeBasis.cpp
+ 101
− 84
View file @ 58bedb0a
#include "QuadEdgeBasis.h"
#include "Legendre.h"
using namespace std;
QuadEdgeBasis::QuadEdgeBasis(const int order){
// Set Basis Type //
this->order = order;
......@@ -8,15 +10,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
type = 1;
dim = 2;
nVertex = 0 ;
nEdge = 4 * (order + 1) ;
nFace = 0 ;
nVertex = 0;
nEdge = 4 * (order + 1);
nFace = 0;
nCell = 2 * (order + 1) * order;
size = 2 * (order + 2) * (order + 1);
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
const int orderPlus = order + 1;
Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus];
......@@ -28,13 +31,29 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
Polynomial lagrange[4];
Polynomial lifting[4];
Polynomial lagrangeSum[4];
Polynomial liftingSub[4];
Polynomial rLiftingSub[4];
Polynomial liftingSub[2][4];
// Integrated and classical Legendre Polynomial //
// Legendre Polynomial //
Legendre::integrated(intLegendre, orderPlus);
Legendre::legendre(legendre, order);
// Vertices definig Edges & Permutations //
const int edgeV[2][4][2] =
{
{ {0, 1}, {1, 2}, {2, 3}, {3, 0} },
{ {1, 0}, {2, 1}, {3, 2}, {0, 3} }
};
// Basis //
node = new vector<vector<Polynomial>*>(nVertex);
edge = new vector<vector<vector<Polynomial>*>*>(2);
face = new vector<vector<vector<Polynomial>*>*>(0);
cell = new vector<vector<Polynomial>*>(nCell);
(*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
(*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
// Lagrange //
lagrange[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
......@@ -53,8 +72,10 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
(Polynomial(1, 0, 1, 0));
// Lagrange Sum //
for(int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4)
lagrangeSum[i] = lagrange[i] + lagrange[j];
for(int e = 0; e < 4; e++)
lagrangeSum[e] =
lagrange[edgeV[0][e][0]] +
lagrange[edgeV[0][e][1]];
// Lifting //
lifting[0] =
......@@ -74,59 +95,48 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
(Polynomial(1, 0, 1, 0));
// Lifting Sub //
for(int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4){
liftingSub[i] = lifting[j] - lifting[i];
rLiftingSub[i] = lifting[i] - lifting[j];
}
for(int e = 0; e < 4; e++){
liftingSub[0][e] =
lifting[edgeV[0][e][0]] -
lifting[edgeV[0][e][1]];
// Basis (temporary --- *no* const) //
std::vector<std::vector<Polynomial>*> basis(size);
std::vector<std::vector<Polynomial>*> revBasis(size);
liftingSub[1][e] =
lifting[edgeV[1][e][0]] -
lifting[edgeV[1][e][1]];
}
// Edge Based (Nedelec) //
int i = 0;
Polynomial oneHalf(0.5, 0, 0, 0);
for(int e = 0; e < 4; e++){
// Direct
basis[i] = new std::vector<Polynomial>(liftingSub[e].gradient());
basis[i]->at(0).mul(lagrangeSum[e]);
basis[i]->at(1).mul(lagrangeSum[e]);
basis[i]->at(2).mul(lagrangeSum[e]);
basis[i]->at(0).mul(oneHalf);
basis[i]->at(1).mul(oneHalf);
basis[i]->at(2).mul(oneHalf);
// Revert
revBasis[i] = new std::vector<Polynomial>(rLiftingSub[e].gradient());
for(int c = 0; c < 2; c++){
for(int e = 0; e < 4; e++){
(*(*edge)[c])[e] =
new vector<Polynomial>(liftingSub[c][e].gradient());
revBasis[i]->at(0).mul(lagrangeSum[e]);
revBasis[i]->at(1).mul(lagrangeSum[e]);
revBasis[i]->at(2).mul(lagrangeSum[e]);
(*(*edge)[c])[e]->at(0).mul(lagrangeSum[e]);
(*(*edge)[c])[e]->at(1).mul(lagrangeSum[e]);
(*(*edge)[c])[e]->at(2).mul(lagrangeSum[e]);
revBasis[i]->at(0).mul(oneHalf);
revBasis[i]->at(1).mul(oneHalf);
revBasis[i]->at(2).mul(oneHalf);
// Next
i++;
(*(*edge)[c])[e]->at(0).mul(oneHalf);
(*(*edge)[c])[e]->at(1).mul(oneHalf);
(*(*edge)[c])[e]->at(2).mul(oneHalf);
}
}
// Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 4; e++){
basis[i] =
new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) *
lagrangeSum[e]).gradient());
revBasis[i] =
new std::vector<Polynomial>((intLegendre[l].compose(rLiftingSub[e]) *
lagrangeSum[e]).gradient());
i++;
for(int c = 0; c < 2; c++){
unsigned int i = 0;
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 4; e++){
(*(*edge)[c])[i + 4] =
new vector<Polynomial>((intLegendre[l].compose(liftingSub[c][e]) *
lagrangeSum[e]).gradient());
i++;
}
}
}
......@@ -139,6 +149,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
px = px - Polynomial(1, 0, 0, 0);
py = py - Polynomial(1, 0, 0, 0);
unsigned int i = 0;
for(int l = 0; l < orderPlus; l++){
iLegendreX[l] = intLegendre[l].compose(px);
iLegendreY[l] = intLegendre[l].compose(py);
......@@ -149,9 +161,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 1) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
basis[i] = new std::vector<Polynomial>((iLegendreX[l1] * iLegendreY[l2]).gradient());
revBasis[i] = basis[i];
(*cell)[i] =
new vector<Polynomial>((iLegendreX[l1] * iLegendreY[l2]).gradient());
i++;
}
......@@ -160,13 +171,11 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 2) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
basis[i] = new std::vector<Polynomial>(3);
basis[i]->at(0) = legendreX[l1] * iLegendreY[l2];
basis[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1;
basis[i]->at(2) = zero;
(*cell)[i] = new vector<Polynomial>(3);
revBasis[i] = basis[i];
(*cell)[i]->at(0) = legendreX[l1] * iLegendreY[l2];
(*cell)[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1;
(*cell)[i]->at(2) = zero;
i++;
}
......@@ -174,18 +183,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 3) //
for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){
basis[i] = new std::vector<Polynomial>(3);
basis[iPlus] = new std::vector<Polynomial>(3);
(*cell)[i] = new vector<Polynomial>(3);
(*cell)[iPlus] = new vector<Polynomial>(3);
basis[i]->at(0) = iLegendreY[l];
basis[i]->at(1) = zero;
basis[i]->at(2) = zero;
(*cell)[i]->at(0) = iLegendreY[l];
(*cell)[i]->at(1) = zero;
(*cell)[i]->at(2) = zero;
basis[iPlus]->at(0) = zero;
basis[iPlus]->at(1) = iLegendreX[l];
basis[iPlus]->at(2) = zero;
revBasis[i] = basis[i];
(*cell)[iPlus]->at(0) = zero;
(*cell)[iPlus]->at(1) = iLegendreX[l];
(*cell)[iPlus]->at(2) = zero;
i++;
}
......@@ -199,24 +206,34 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
delete[] iLegendreY;
delete[] legendreX;
delete[] legendreY;
// Set Basis //
this->basis = new std::vector<const std::vector<Polynomial>*>
(basis.begin(), basis.end());
this->revBasis = new std::vector<const std::vector<Polynomial>*>
(revBasis.begin(), revBasis.end());
}
QuadEdgeBasis::~QuadEdgeBasis(void){
for(int i = 0; i < size; i++){
delete (*basis)[i];
// Vertex Based //
for(int i = 0; i < nVertex; i++)
delete (*node)[i];
delete node;
if(i >= nVertex && i < nVertex + nEdge)
delete (*revBasis)[i];
// Edge Based //
for(int c = 0; c < 2; c++){
for(int i = 0; i < nEdge; i++)
delete (*(*edge)[c])[i];
delete (*edge)[c];
}
delete edge;
// Face Based //
delete face;
// Cell Based //
for(int i = 0; i < nCell; i++)
delete (*cell)[i];
delete basis;
delete revBasis;
delete cell;
}
FunctionSpace/QuadNodeBasis.cpp
+ 74
− 35
View file @ 58bedb0a
#include "QuadNodeBasis.h"
#include "Legendre.h"
using namespace std;
QuadNodeBasis::QuadNodeBasis(const int order){
// Set Basis Type //
this->order = order;
......@@ -8,20 +10,38 @@ QuadNodeBasis::QuadNodeBasis(const int order){
type = 0;
dim = 2;
nVertex = 4 ;
nEdge = 4 * (order - 1) ;
nFace = 0 ;
nVertex = 4;
nEdge = 4 * (order - 1);
nFace = 0;
nCell = (order - 1) * (order - 1);
size = (order + 1) * (order + 1);
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
Polynomial lifting[4];
Polynomial liftingSub[2][4];
// Legendre Polynomial //
Legendre::integrated(legendre, order);
// Vertices definig Edges & Permutations //
const int edgeV[2][4][2] =
{
{ {0, 1}, {1, 2}, {2, 3}, {3, 0} },
{ {1, 0}, {2, 1}, {3, 2}, {0, 3} }
};
// Basis //
node = new vector<Polynomial*>(nVertex);
edge = new vector<vector<Polynomial*>*>(2);
face = new vector<vector<Polynomial*>*>(0);
cell = new vector<Polynomial*>(nCell);
(*edge)[0] = new vector<Polynomial*>(nEdge);
(*edge)[1] = new vector<Polynomial*>(nEdge);
// Lifting //
lifting[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
......@@ -39,47 +59,48 @@ QuadNodeBasis::QuadNodeBasis(const int order){
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0));
// Lifting Sub //
for(int e = 0; e < 4; e++){
liftingSub[0][e] =
lifting[edgeV[0][e][0]] -
lifting[edgeV[0][e][1]];
liftingSub[1][e] =
lifting[edgeV[1][e][0]] -
lifting[edgeV[1][e][1]];
}
// Basis //
basis = new std::vector<const Polynomial*>(size);
revBasis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) //
(*basis)[0] =
(*node)[0] =
new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
(*basis)[1] =
(*node)[1] =
new Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
(*basis)[2] =
(*node)[2] =
new Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)));
(*basis)[3] =
(*node)[3] =
new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)));
// Vertex Based (Revert) //
for(int i = 0; i < 3; i++)
(*revBasis)[i] = (*basis)[i];
// Edge Based //
int i = 4;
for(int c = 0; c < 2; c++){
unsigned int i = 0;
for(int l = 1; l < order; l++){
for(int e1 = 0, e2 = 1; e1 < 4; e1++, e2 = (e2 + 1) % 4){
(*basis)[i] =
new Polynomial(legendre[l].compose(lifting[e2] - lifting[e1]) *
(*(*basis)[e1] + *(*basis)[e2]));
for(int l = 1; l < order; l++){
for(int e = 0; e < 4; e++){
(*(*edge)[c])[i] =
new Polynomial(legendre[l].compose(liftingSub[c][e]) *
(*(*node)[edgeV[c][e][0]] + *(*node)[edgeV[c][e][1]]));
(*revBasis)[i] =
new Polynomial(legendre[l].compose(lifting[e1] - lifting[e2]) *
(*(*basis)[e1] + *(*basis)[e2]));
i++;
i++;
}
}
}
......@@ -91,13 +112,13 @@ QuadNodeBasis::QuadNodeBasis(const int order){
px = px - Polynomial(1, 0, 0, 0);
py = py - Polynomial(1, 0, 0, 0);
unsigned int i = 0;
for(int l1 = 1; l1 < order; l1++){
for(int l2 = 1; l2 < order; l2++){
(*basis)[i] =
(*cell)[i] =
new Polynomial(legendre[l1].compose(px) * legendre[l2].compose(py));
(*revBasis)[i] = (*basis)[i];
i++;
}
}
......@@ -108,13 +129,31 @@ QuadNodeBasis::QuadNodeBasis(const int order){
}
QuadNodeBasis::~QuadNodeBasis(void){
for(int i = 0; i < size; i++){
delete (*basis)[i];
// Vertex Based //
for(int i = 0; i < nVertex; i++)
delete (*node)[i];
delete node;
if(i >= nVertex && i < nVertex + nEdge)
delete (*revBasis)[i];
// Edge Based //
for(int c = 0; c < 2; c++){
for(int i = 0; i < nEdge; i++)
delete (*(*edge)[c])[i];
delete (*edge)[c];
}
delete edge;
// Face Based //
delete face;
// Cell Based //
for(int i = 0; i < nCell; i++)
delete (*cell)[i];
delete basis;
delete revBasis;
delete cell;
}
FunctionSpace/TetEdgeBasis.cpp
+ 3
− 2
View file @ 58bedb0a
......@@ -28,7 +28,7 @@ TetEdgeBasis::TetEdgeBasis(int order){
Polynomial* sclLegendre = new Polynomial[orderMinus];
Polynomial* intLegendre = new Polynomial[orderPlus];
// Classical and Intrated-Scaled Legendre Polynomial //
// Legendre Polynomial //
Legendre::legendre(legendre, orderMinusTwo);
Legendre::scaled(sclLegendre, orderMinusTwo);
Legendre::intScaled(intLegendre, orderPlus);
......@@ -97,7 +97,8 @@ TetEdgeBasis::TetEdgeBasis(int order){
tmp[1].mul(lagrange[edgeV[c][e][0]]);
tmp[2].mul(lagrange[edgeV[c][e][0]]);
(*(*edge)[c])[e] = new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
(*(*edge)[c])[e] =
new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
(*(*edge)[c])[e]->at(0).mul(lagrange[edgeV[c][e][1]]);
(*(*edge)[c])[e]->at(1).mul(lagrange[edgeV[c][e][1]]);
......
FunctionSpace/TetNodeBasis.cpp
+ 59
− 110
View file @ 58bedb0a
......@@ -18,39 +18,39 @@ TetNodeBasis::TetNodeBasis(int order){
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
Polynomial* sclLegendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order];
// Classical, Scaled and Intrated-Scaled Legendre Polynomial //
const int orderMinus = order - 1;
const int orderMinusTwo = order - 2;
const int orderMinusThree = order - 3;
Polynomial* legendre = new Polynomial[order];
Polynomial* sclLegendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order];
// Legendre Polynomial //
Legendre::legendre(legendre, orderMinus);
Legendre::scaled(sclLegendre, orderMinus);
Legendre::intScaled(intLegendre, order);
// Vertices definig Edges //
const int edgeV[6][2] = {
{0, 1},
{1, 2},
{2, 0},
{3, 0},
{3, 2},
{3, 1}
};
// Vertices definig Faces //
const int faceV[4][3] = {
{0, 2, 1},
{0, 1, 3},
{0, 3, 2},
{3, 1, 2}
};
// Counter //
int i;
// Vertices definig Edges & Permutations //
const int edgeV[2][6][2] =
{
{ {0, 1}, {1, 2}, {2, 0},
{3, 0}, {3, 2}, {3, 1} },
{ {1, 0}, {2, 1}, {0, 2},
{0, 3}, {2, 3}, {1, 3} }
};
// Vertices definig Faces & Permutations //
const int faceV[6][4][3] =
{
{ {0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2} },
{ {2, 0, 1}, {1, 0, 3}, {3, 0, 2}, {1, 3, 2} },
{ {2, 1, 0}, {1, 3, 0}, {3, 2, 0}, {1, 2, 3} },
{ {1, 2, 0}, {3, 1, 0}, {2, 3, 0}, {2, 1, 3} },
{ {1, 0, 2}, {3, 0, 1}, {2, 0, 3}, {2, 3, 1} },
{ {0, 1, 2}, {0, 3, 1}, {0, 2, 3}, {3, 2, 1} },
};
// Basis //
node = new vector<Polynomial*>(nVertex);
......@@ -87,97 +87,46 @@ TetNodeBasis::TetNodeBasis(int order){
// Edge Based //
i = 0;
for(int l = 1; l < order; l++){
for(int e = 0; e < 6; e++){
(*(*edge)[0])[i] =
new Polynomial(intLegendre[l].compose
(*(*node)[edgeV[e][0]] - *(*node)[edgeV[e][1]],
*(*node)[edgeV[e][0]] + *(*node)[edgeV[e][1]]));
(*(*edge)[1])[i] =
new Polynomial(intLegendre[l].compose
(*(*node)[edgeV[e][1]] - *(*node)[edgeV[e][0]],
*(*node)[edgeV[e][1]] + *(*node)[edgeV[e][0]]));
i++;
for(int c = 0; c < 2; c++){
unsigned int i = 0;
for(int l = 1; l < order; l++){
for(int e = 0; e < 6; e++){
(*(*edge)[c])[i] =
new Polynomial(intLegendre[l].compose
(*(*node)[edgeV[c][e][0]] - *(*node)[edgeV[c][e][1]],
*(*node)[edgeV[c][e][0]] + *(*node)[edgeV[c][e][1]]));
i++;
}
}
}
// Face Based //
i = 0;
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
for(int m = 0; m < 4; m++){
Polynomial sum =
*(*node)[faceV[m][0]] +
*(*node)[faceV[m][1]] +
*(*node)[faceV[m][2]];
(*(*face)[0])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[faceV[m][0]] - *(*node)[faceV[m][1]],
*(*node)[faceV[m][0]] + *(*node)[faceV[m][1]]) *
*(*node)[faceV[m][2]] *
sclLegendre[l2].compose
(*(*node)[faceV[m][2]] * 2 - sum, sum));
(*(*face)[1])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[faceV[m][1]] - *(*node)[faceV[m][0]],
*(*node)[faceV[m][1]] + *(*node)[faceV[m][0]]) *
*(*node)[faceV[m][2]] *
sclLegendre[l2].compose
(*(*node)[faceV[m][2]] * 2 - sum, sum));
(*(*face)[2])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[faceV[m][1]] - *(*node)[faceV[m][2]],
*(*node)[faceV[m][1]] + *(*node)[faceV[m][2]]) *
*(*node)[faceV[m][0]] *
sclLegendre[l2].compose
(*(*node)[faceV[m][0]] * 2 - sum, sum));
(*(*face)[3])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[faceV[m][2]] - *(*node)[faceV[m][1]],
*(*node)[faceV[m][2]] + *(*node)[faceV[m][1]]) *
*(*node)[faceV[m][0]] *
sclLegendre[l2].compose
(*(*node)[faceV[m][0]] * 2 - sum, sum));
(*(*face)[4])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[faceV[m][2]] - *(*node)[faceV[m][0]],
*(*node)[faceV[m][2]] + *(*node)[faceV[m][0]]) *
*(*node)[faceV[m][1]] *
sclLegendre[l2].compose
(*(*node)[faceV[m][1]] * 2 - sum, sum));
(*(*face)[5])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[faceV[m][0]] - *(*node)[faceV[m][2]],
*(*node)[faceV[m][0]] + *(*node)[faceV[m][2]]) *
for(int c = 0; c < 6; c++){
unsigned int i = 0;
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
for(int f = 0; f < 4; f++){
Polynomial sum =
*(*node)[faceV[c][f][0]] +
*(*node)[faceV[c][f][1]] +
*(*node)[faceV[c][f][2]];
(*(*face)[c])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[faceV[c][f][0]] - *(*node)[faceV[c][f][1]],
*(*node)[faceV[c][f][0]] + *(*node)[faceV[c][f][1]]) *
*(*node)[faceV[m][1]] *
*(*node)[faceV[c][f][2]] *
sclLegendre[l2].compose
(*(*node)[faceV[m][1]] * 2 - sum, sum));
i++;
sclLegendre[l2].compose
(*(*node)[faceV[c][f][2]] * 2 - sum, sum));
i++;
}
}
}
}
......@@ -191,7 +140,7 @@ TetNodeBasis::TetNodeBasis(int order){
Polynomial sub = *(*node)[0] - *(*node)[1];
Polynomial add = *(*node)[0] + *(*node)[1];
i = 0;
unsigned int i = 0;
for(int l1 = 1; l1 < orderMinusTwo; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
......
FunctionSpace/TriEdgeBasis.cpp
+ 66
− 70
View file @ 58bedb0a
......@@ -15,7 +15,7 @@ TriEdgeBasis::TriEdgeBasis(int order){
nFace = 0;
nCell = ((order - 1) * order + order - 1);
size = (order + 1) * (order + 2);
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
const int orderPlus = order + 1;
......@@ -30,30 +30,16 @@ TriEdgeBasis::TriEdgeBasis(int order){
Polynomial lagrangeSum [3];
Polynomial lagrangeSub [2][3];
// Classical and Intrated-Scaled Legendre Polynomial //
// Legendre Polynomial //
Legendre::legendre(legendre, order);
Legendre::intScaled(intLegendre, orderPlus);
// Lagrange //
lagrange[0] =
Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
lagrange[1] =
Polynomial(1, 1, 0, 0);
lagrange[2] =
Polynomial(1, 0, 1, 0);
// Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
lagrangeSum[i] = lagrange[j] + lagrange[i];
// Lagrange Sub (& revert) //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
lagrangeSub[0][i] = lagrange[i] - lagrange[j];
lagrangeSub[1][i] = lagrange[j] - lagrange[i];
}
// Vertices definig Edges & Permutations //
const int edgeV[2][3][2] =
{
{ {0, 1}, {1, 2}, {2, 0} },
{ {1, 0}, {2, 1}, {0, 2} }
};
// Basis //
node = new vector<vector<Polynomial>*>(nVertex);
......@@ -63,66 +49,76 @@ TriEdgeBasis::TriEdgeBasis(int order){
(*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
(*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
// Counter
unsigned int i = 0;
// Edge Based (Nedelec) //
for(int j = 1; i < 3; j = (j + 1) % 3){
// Temp
vector<Polynomial> tmp = lagrange[j].gradient();
// Edge[0]
tmp[0].mul(lagrange[i]);
tmp[1].mul(lagrange[i]);
tmp[2].mul(lagrange[i]);
(*(*edge)[0])[i] = new vector<Polynomial>(lagrange[i].gradient());
// Lagrange //
lagrange[0] =
Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
(*(*edge)[0])[i]->at(0).mul(lagrange[j]);
(*(*edge)[0])[i]->at(1).mul(lagrange[j]);
(*(*edge)[0])[i]->at(2).mul(lagrange[j]);
lagrange[1] =
Polynomial(1, 1, 0, 0);
(*(*edge)[0])[i]->at(0).sub(tmp[0]);
(*(*edge)[0])[i]->at(1).sub(tmp[1]);
(*(*edge)[0])[i]->at(2).sub(tmp[2]);
lagrange[2] =
Polynomial(1, 0, 1, 0);
//Edge[1]
tmp = lagrange[i].gradient();
// Lagrange Sum //
for(int e = 0; e < 3; e++)
lagrangeSum[e] =
lagrange[edgeV[0][e][0]] +
lagrange[edgeV[0][e][1]];
// Lagrange Sub //
for(int e = 0; e < 3; e++){
lagrangeSub[0][e] =
lagrange[edgeV[0][e][0]] -
lagrange[edgeV[0][e][1]];
lagrangeSub[1][e] =
lagrange[edgeV[1][e][0]] -
lagrange[edgeV[1][e][1]];
}
tmp[0].mul(lagrange[j]);
tmp[1].mul(lagrange[j]);
tmp[2].mul(lagrange[j]);
// Edge Based (Nedelec) //
for(int c = 0; c < 2; c++){
for(int e = 0; e < 3; e++){
vector<Polynomial> tmp = lagrange[edgeV[c][e][1]].gradient();
tmp[0].mul(lagrange[edgeV[c][e][0]]);
tmp[1].mul(lagrange[edgeV[c][e][0]]);
tmp[2].mul(lagrange[edgeV[c][e][0]]);
(*(*edge)[c])[e] =
new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
(*(*edge)[1])[i] = new vector<Polynomial>(lagrange[j].gradient());
(*(*edge)[c])[e]->at(0).mul(lagrange[edgeV[c][e][1]]);
(*(*edge)[c])[e]->at(1).mul(lagrange[edgeV[c][e][1]]);
(*(*edge)[c])[e]->at(2).mul(lagrange[edgeV[c][e][1]]);
(*(*edge)[c])[e]->at(0).sub(tmp[0]);
(*(*edge)[c])[e]->at(1).sub(tmp[1]);
(*(*edge)[c])[e]->at(2).sub(tmp[2]);
}
}
(*(*edge)[1])[i]->at(0).mul(lagrange[i]);
(*(*edge)[1])[i]->at(1).mul(lagrange[i]);
(*(*edge)[1])[i]->at(2).mul(lagrange[i]);
(*(*edge)[1])[i]->at(0).sub(tmp[0]);
(*(*edge)[1])[i]->at(1).sub(tmp[1]);
(*(*edge)[1])[i]->at(2).sub(tmp[2]);
// Next
i++;
}
// Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 3; e++){
(*(*edge)[0])[i] =
new vector<Polynomial>((intLegendre[l].compose(lagrangeSub[0][e],
lagrangeSum[e])).gradient());
for(int c = 0; c < 2; c++){
unsigned int i = 0;
(*(*edge)[1])[i] =
new vector<Polynomial>((intLegendre[l].compose(lagrangeSub[1][e],
lagrangeSum[e])).gradient());
i++;
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 3; e++){
(*(*edge)[c])[i + 3] =
new vector<Polynomial>
((intLegendre[l].compose(lagrangeSub[c][e],
lagrangeSum[e])).gradient());
i++;
}
}
}
// Cell Based (Preliminary) //
Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
......@@ -132,7 +128,7 @@ TriEdgeBasis::TriEdgeBasis(int order){
v[l].mul(lagrange[2]);
}
i = 0;
unsigned int i = 0;
// Cell Based (Type 1) //
for(int l1 = 1; l1 < order; l1++){
......
FunctionSpace/TriNedelecBasis.cpp
+ 37
− 46
View file @ 58bedb0a
......@@ -16,6 +16,23 @@ TriNedelecBasis::TriNedelecBasis(void){
size = 3;
// Vertices definig Edges & Permutations //
const int edgeV[2][3][2] =
{
{ {0, 1}, {1, 2}, {2, 0} },
{ {1, 0}, {2, 1}, {0, 2} }
};
// Basis //
node = new vector<vector<Polynomial>*>(nVertex);
edge = new vector<vector<vector<Polynomial>*>*>(2);
face = new vector<vector<vector<Polynomial>*>*>(0);
cell = new vector<vector<Polynomial>*>(nCell);
(*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
(*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
// Lagrange //
Polynomial lagrange[3];
......@@ -29,53 +46,27 @@ TriNedelecBasis::TriNedelecBasis(void){
Polynomial(1, 0, 1, 0);
// Basis //
node = new vector<vector<Polynomial>*>(nVertex);
edge = new vector<vector<vector<Polynomial>*>*>(2);
face = new vector<vector<vector<Polynomial>*>*>(0);
cell = new vector<vector<Polynomial>*>(nCell);
(*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
(*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
// Nedelec //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
// Temp
vector<Polynomial> tmp = lagrange[j].gradient();
// Edge[0]
tmp[0].mul(lagrange[i]);
tmp[1].mul(lagrange[i]);
tmp[2].mul(lagrange[i]);
(*(*edge)[0])[i] = new vector<Polynomial>(lagrange[i].gradient());
(*(*edge)[0])[i]->at(0).mul(lagrange[j]);
(*(*edge)[0])[i]->at(1).mul(lagrange[j]);
(*(*edge)[0])[i]->at(2).mul(lagrange[j]);
(*(*edge)[0])[i]->at(0).sub(tmp[0]);
(*(*edge)[0])[i]->at(1).sub(tmp[1]);
(*(*edge)[0])[i]->at(2).sub(tmp[2]);
//Edge[1]
tmp = lagrange[i].gradient();
tmp[0].mul(lagrange[j]);
tmp[1].mul(lagrange[j]);
tmp[2].mul(lagrange[j]);
(*(*edge)[1])[i] = new vector<Polynomial>(lagrange[j].gradient());
(*(*edge)[1])[i]->at(0).mul(lagrange[i]);
(*(*edge)[1])[i]->at(1).mul(lagrange[i]);
(*(*edge)[1])[i]->at(2).mul(lagrange[i]);
(*(*edge)[1])[i]->at(0).sub(tmp[0]);
(*(*edge)[1])[i]->at(1).sub(tmp[1]);
(*(*edge)[1])[i]->at(2).sub(tmp[2]);
}
for(int c = 0; c < 2; c++){
for(int e = 0; e < 3; e++){
vector<Polynomial> tmp = lagrange[edgeV[c][e][1]].gradient();
tmp[0].mul(lagrange[edgeV[c][e][0]]);
tmp[1].mul(lagrange[edgeV[c][e][0]]);
tmp[2].mul(lagrange[edgeV[c][e][0]]);
(*(*edge)[c])[e] =
new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
(*(*edge)[c])[e]->at(0).mul(lagrange[edgeV[c][e][1]]);
(*(*edge)[c])[e]->at(1).mul(lagrange[edgeV[c][e][1]]);
(*(*edge)[c])[e]->at(2).mul(lagrange[edgeV[c][e][1]]);
(*(*edge)[c])[e]->at(0).sub(tmp[0]);
(*(*edge)[c])[e]->at(1).sub(tmp[1]);
(*(*edge)[c])[e]->at(2).sub(tmp[2]);
}
}
}
TriNedelecBasis::~TriNedelecBasis(void){
......
FunctionSpace/TriNodeBasis.cpp
+ 27
− 13
View file @ 58bedb0a
......@@ -15,21 +15,27 @@ TriNodeBasis::TriNodeBasis(int order){
nFace = 0;
nCell = (order - 1) * (order - 2) / 2;
size = (order + 1) * (order + 2) / 2;
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order];
const int orderMinus = order - 1;
Polynomial lagrangeSum[3];
Polynomial lagrangeSub[2][3];
Polynomial* legendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order];
// Classical and Intrated-Scaled Legendre Polynomial //
const int orderMinus = order - 1;
Polynomial lagrangeSum[3];
Polynomial lagrangeSub[2][3];
// Legendre Polynomial //
Legendre::legendre(legendre, orderMinus);
Legendre::intScaled(intLegendre, order);
// Vertices definig Edges & Permutations //
const int edgeV[2][3][2] =
{
{ {0, 1}, {1, 2}, {2, 0} },
{ {1, 0}, {2, 1}, {0, 2} }
};
// Basis //
node = new vector<Polynomial*>(nVertex);
......@@ -55,13 +61,20 @@ TriNodeBasis::TriNodeBasis(int order){
// Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
lagrangeSum[i] = *(*node)[i] + *(*node)[j];
for(int e = 0; e < 3; e++)
lagrangeSum[e] =
*(*node)[edgeV[0][e][0]] +
*(*node)[edgeV[0][e][1]];
// Lagrange Sub //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
lagrangeSub[0][i] = *(*node)[j] - *(*node)[i];
lagrangeSub[1][i] = *(*node)[i] - *(*node)[j];
for(int e = 0; e < 3; e++){
lagrangeSub[0][e] =
*(*node)[edgeV[0][e][0]] -
*(*node)[edgeV[0][e][1]];
lagrangeSub[1][e] =
*(*node)[edgeV[1][e][0]] -
*(*node)[edgeV[1][e][1]];
}
......@@ -96,6 +109,7 @@ TriNodeBasis::TriNodeBasis(int order){
}
}
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment