Select Git revision
linearSystemGmm.h
TetNodeBasis.cpp 4.86 KiB
#include "TetNodeBasis.h"
#include "TetReferenceSpace.h"
#include "Legendre.h"
using namespace std;
TetNodeBasis::TetNodeBasis(size_t order){
// Reference Space //
refSpace = new TetReferenceSpace;
nRefSpace = getReferenceSpace().getNReferenceSpace();
const vector<vector<vector<size_t> > >&
edgeIdx = refSpace->getEdgeNodeIndex();
const vector<vector<vector<size_t> > >&
faceIdx = refSpace->getFaceNodeIndex();
// Set Basis Type //
this->order = order;
type = 0;
dim = 3;
nVertex = 4;
nEdge = 6 * (order - 1);
nFace = 2 * (order - 1) * (order - 2);
nCell = (order - 1) * (order - 2) * (order - 3) / 6;
nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
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);
// Lagrange Polynomial //
const Polynomial lagrange[4] =
{
Polynomial(Polynomial(1, 0, 0, 0) -
Polynomial(1, 1, 0, 0) -
Polynomial(1, 0, 1, 0) -
Polynomial(1, 0, 0, 1)),
Polynomial(Polynomial(1, 1, 0, 0)),
Polynomial(Polynomial(1, 0, 1, 0)),
Polynomial(Polynomial(1, 0, 0, 1))
};
// Basis //
basis = new Polynomial**[nRefSpace];
for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new Polynomial*[nFunction];
// Vertex Based //
for(size_t s = 0; s < nRefSpace; s++){
basis[s][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]);
basis[s][3] = new Polynomial(lagrange[3]);
}
// Edge Based //
for(size_t s = 0; s < nRefSpace; s++){
size_t i = nVertex;
for(int e = 0; e < 6; e++){
for(size_t l = 1; l < order; l++){
basis[s][i] =
new Polynomial(intLegendre[l].compose
(lagrange[edgeIdx[s][e][0]] -
lagrange[edgeIdx[s][e][1]]
,
lagrange[edgeIdx[s][e][0]] +
lagrange[edgeIdx[s][e][1]]));
i++;
}
}
}
// Face Based //
// TO CHECK: Are Triangles face matching tets ?
for(size_t s = 0; s < nRefSpace; s++){
size_t i = nVertex + nEdge;
for(int f = 0; f < 4; f++){
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
Polynomial sum =
lagrange[faceIdx[s][f][0]] +
lagrange[faceIdx[s][f][1]] +
lagrange[faceIdx[s][f][2]];
basis[s][i] =
new Polynomial(intLegendre[l1].compose
(lagrange[faceIdx[s][f][0]] -
lagrange[faceIdx[s][f][1]]
,
lagrange[faceIdx[s][f][0]] +
lagrange[faceIdx[s][f][1]])
*
lagrange[faceIdx[s][f][2]]
*
sclLegendre[l2].compose
(lagrange[faceIdx[s][f][2]] * 2 - sum, sum));
i++;
}
}
}
}
// Cell Based //
const Polynomial oneMinusFour = Polynomial(1, 0, 0, 0) - lagrange[3];
const Polynomial twoThreeOneMinusFour = lagrange[2] * 2 - oneMinusFour;
const Polynomial twoFourMinusOne = lagrange[3] * 2 - Polynomial(1, 0, 0, 0);
const Polynomial sub = lagrange[0] - lagrange[1];
const Polynomial add = lagrange[0] + lagrange[1];
for(size_t s = 0; s < nRefSpace; s++){
size_t i = nVertex + nEdge + nFace;
for(int l1 = 1; l1 < orderMinusTwo; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
basis[s][i] =
new Polynomial(intLegendre[l1].compose(sub, add) *
lagrange[2] *
sclLegendre[l2].compose(twoThreeOneMinusFour,
oneMinusFour) *
lagrange[3] *
legendre[l3].compose(twoFourMinusOne));
i++;
}
}
}
}
// Free Temporary Sapce //
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
}
TetNodeBasis::~TetNodeBasis(void){
// ReferenceSpace //
delete refSpace;
// Basis //
for(size_t i = 0; i < nRefSpace; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
}