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

TetEdgeBasis: Problem with orientations -- RefBasis is wrong -- Node sorting... 

TetEdgeBasis: Problem with orientations -- RefBasis is wrong -- Node sorting is not the solution :-(
parent f53fb375
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 218 additions and 259 deletions
FunctionSpace/BasisGenerator.cpp
+ 1
− 1
View file @ 708a843e
......@@ -123,7 +123,7 @@ Basis* BasisGenerator::tetZaglmayrGen(unsigned int basisType,
unsigned int order){
switch(basisType){
case 0: return new TetNodeBasis(order);
//case 1: return new TetEdgeBasis(order);
case 1: return new TetEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Tetrahedrons");
case 3: throw Exception("3-form not implemented on Tetrahedrons");
......
FunctionSpace/CMakeLists.txt
+ 1
− 1
View file @ 708a843e
......@@ -40,7 +40,7 @@ set(SRC
# HexEdgeBasis.cpp
TetNodeBasis.cpp
# TetEdgeBasis.cpp
TetEdgeBasis.cpp
FunctionSpace.cpp
FunctionSpaceScalar.cpp
......
FunctionSpace/EvaluatedBasis.cpp
+ 1
− 1
View file @ 708a843e
FunctionSpace/EvaluatedBasis.h
+ 6
− 0
View file @ 708a843e
......@@ -100,6 +100,12 @@ EvaluatedBasis::getEvaluation(unsigned int refSpace) const{
inline const fullMatrix<double>&
EvaluatedBasis::getEvaluation(const MElement& element) const{
/*
std::cout << element.getNum()
<< ": "
<< refSpace->getReferenceSpace(element)
<< std::endl;
*/
return *(*eBasis)[refSpace->getReferenceSpace(element)];
}
......
FunctionSpace/FunctionSpace.cpp
+ 1
− 1
View file @ 708a843e
FunctionSpace/ReferenceSpace.cpp
+ 1
− 1
View file @ 708a843e
FunctionSpace/TetEdgeBasis.cpp
+ 206
− 253
View file @ 708a843e
#include "TetEdgeBasis.h"
#include "TetReferenceSpace.h"
#include "Legendre.h"
using namespace std;
TetEdgeBasis::TetEdgeBasis(int order){
TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Reference Space //
refSpace = new TetReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
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;
......@@ -14,19 +25,13 @@ TetEdgeBasis::TetEdgeBasis(int order){
nEdge = 6 * (order + 1);
nFace = 4 * (order + 1) * (order - 1);
nCell = (order + 1) * (order - 1) * (order - 2) / 2;
nEdgeClosure = 2;
nFaceClosure = 6;
size = nVertex + nEdge + nFace + nCell;
nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
const int orderPlus = order + 1;
const int orderMinus = order - 1;
const int orderMinusTwo = order - 2;
Polynomial lagrange[4];
Polynomial* legendre = new Polynomial[orderMinus];
Polynomial* sclLegendre = new Polynomial[orderMinus];
Polynomial* intLegendre = new Polynomial[orderPlus];
......@@ -36,121 +41,90 @@ TetEdgeBasis::TetEdgeBasis(int order){
Legendre::scaled(sclLegendre, orderMinusTwo);
Legendre::intScaled(intLegendre, orderPlus);
// Vertices definig Edges & Permutations //
const int edgeV[2][6][2] =
// Lagrange Polynomial //
const Polynomial lagrange[4] =
{
{ {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<vector<Polynomial>*>(nVertex);
edge = new vector<vector<vector<Polynomial>*>*>(2);
face = new vector<vector<vector<Polynomial>*>*>(6);
cell = new vector<vector<Polynomial>*>(nCell);
(*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
(*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
(*face)[0] = new vector<vector<Polynomial>*>(nFace);
(*face)[1] = new vector<vector<Polynomial>*>(nFace);
(*face)[2] = new vector<vector<Polynomial>*>(nFace);
(*face)[3] = new vector<vector<Polynomial>*>(nFace);
(*face)[4] = new vector<vector<Polynomial>*>(nFace);
(*face)[5] = new vector<vector<Polynomial>*>(nFace);
// Lagrange //
lagrange[0] =
Polynomial(1, 0, 0, 0) -
Polynomial(Polynomial(1, 0, 0, 0) -
Polynomial(1, 1, 0, 0) -
Polynomial(1, 0, 1, 0) -
Polynomial(1, 0, 0, 1);
Polynomial(1, 0, 0, 1)),
Polynomial(Polynomial(1, 1, 0, 0)),
lagrange[1] =
Polynomial(1, 1, 0, 0);
Polynomial(Polynomial(1, 0, 1, 0)),
lagrange[2] =
Polynomial(1, 0, 1, 0);
Polynomial(Polynomial(1, 0, 0, 1))
};
lagrange[3] =
Polynomial(1, 0, 0, 1);
// Basis //
basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
for(unsigned int s = 0; s < nRefSpace; s++)
(*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
// Edge Based (Nedelec) //
for(int c = 0; c < 2; c++){
for(unsigned int s = 0; s < nRefSpace; s++){
for(int e = 0; e < 6; e++){
vector<Polynomial> tmp = lagrange[edgeV[c][e][1]].gradient();
vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].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]]);
tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
(*(*edge)[c])[e] =
new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
(*(*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]]);
tmp1[0].sub(tmp2[0]);
tmp1[1].sub(tmp2[1]);
tmp1[2].sub(tmp2[2]);
(*(*edge)[c])[e]->at(0).sub(tmp[0]);
(*(*edge)[c])[e]->at(1).sub(tmp[1]);
(*(*edge)[c])[e]->at(2).sub(tmp[2]);
(*(*basis)[s])[e] = new vector<Polynomial>(tmp1);
}
}
// Edge Based (High Order) //
for(int c = 0; c < 2; c++){
unsigned int i = 0;
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = 6;
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 6; e++){
(*(*edge)[c])[i + 6] =
(*(*basis)[s])[i] =
new vector<Polynomial>
(intLegendre[l].compose
(lagrange[edgeV[c][e][0]] - lagrange[edgeV[c][e][1]],
lagrange[edgeV[c][e][0]] + lagrange[edgeV[c][e][1]]).gradient());
((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
lagrange[(*(*edgeV[s])[e])[1]]
,
lagrange[(*(*edgeV[s])[e])[0]] +
lagrange[(*(*edgeV[s])[e])[1]])).gradient());
i++;
}
}
}
// Face Based //
for(int c = 0; c < 6; c++){
unsigned int i = 0;
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge;
for(int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
for(unsigned int l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
for(int f = 0; f < 4; f++){
// Preliminary Type 1
Polynomial sum =
lagrange[faceV[c][f][0]] +
lagrange[faceV[c][f][1]] +
lagrange[faceV[c][f][2]];
Polynomial u = intLegendre[l1].compose
(lagrange[faceV[c][f][0]] - lagrange[faceV[c][f][1]],
lagrange[faceV[c][f][0]] + lagrange[faceV[c][f][1]]);
Polynomial v = lagrange[faceV[c][f][2]] * sclLegendre[l2].compose
(lagrange[faceV[c][f][2]] * 2 - sum, sum);
lagrange[(*(*faceV[s])[f])[0]] +
lagrange[(*(*faceV[s])[f])[1]] +
lagrange[(*(*faceV[s])[f])[2]];
Polynomial u =
intLegendre[l1].compose(lagrange[(*(*faceV[s])[f])[0]] -
lagrange[(*(*faceV[s])[f])[1]]
,
lagrange[(*(*faceV[s])[f])[0]] +
lagrange[(*(*faceV[s])[f])[1]]);
Polynomial v =
lagrange[(*(*faceV[s])[f])[2]] *
sclLegendre[l2].compose(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum);
// Preliminary Type 2
vector<Polynomial> gradU = u.gradient();
......@@ -172,18 +146,18 @@ TetEdgeBasis::TetEdgeBasis(int order){
subGradUV[2].sub(uGradV[2]);
// Preliminary Type 3
vector<Polynomial> gradL1 = lagrange[faceV[c][f][0]].gradient();
vector<Polynomial> gradL2 = lagrange[faceV[c][f][1]].gradient();
vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[f])[0]].gradient();
vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[f])[1]].gradient();
vector<Polynomial> l2GradL1(gradL1);
l2GradL1[0].mul(lagrange[faceV[c][f][1]]);
l2GradL1[1].mul(lagrange[faceV[c][f][1]]);
l2GradL1[2].mul(lagrange[faceV[c][f][1]]);
l2GradL1[0].mul(lagrange[(*(*faceV[s])[f])[1]]);
l2GradL1[1].mul(lagrange[(*(*faceV[s])[f])[1]]);
l2GradL1[2].mul(lagrange[(*(*faceV[s])[f])[1]]);
vector<Polynomial> l1GradL2(gradL2);
l1GradL2[0].mul(lagrange[faceV[c][f][0]]);
l1GradL2[1].mul(lagrange[faceV[c][f][0]]);
l1GradL2[2].mul(lagrange[faceV[c][f][0]]);
l1GradL2[0].mul(lagrange[(*(*faceV[s])[f])[0]]);
l1GradL2[1].mul(lagrange[(*(*faceV[s])[f])[0]]);
l1GradL2[2].mul(lagrange[(*(*faceV[s])[f])[0]]);
vector<Polynomial> subGradL1L2V(l2GradL1);
subGradL1L2V[0].sub(l1GradL2[0]);
......@@ -196,20 +170,20 @@ TetEdgeBasis::TetEdgeBasis(int order){
// Type 1
(*(*face)[c])[i] =
(*(*basis)[s])[i] =
new vector<Polynomial>((u * v).gradient());
i++;
// Type 2
(*(*face)[c])[i] =
(*(*basis)[s])[i] =
new vector<Polynomial>(subGradUV);
i++;
// Type 3
if(l1 == 1){
(*(*face)[c])[i] =
(*(*basis)[s])[i] =
new vector<Polynomial>(subGradL1L2V);
i++;
......@@ -219,11 +193,11 @@ TetEdgeBasis::TetEdgeBasis(int order){
}
}
// Cell Based //
Polynomial one(1, 0, 0, 0);
const Polynomial one(1, 0, 0, 0);
unsigned int i = 0;
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge + nFace;
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
......@@ -318,26 +292,26 @@ TetEdgeBasis::TetEdgeBasis(int order){
// Type 1
(*cell)[i] = new vector<Polynomial>((u * v * w).gradient());
(*(*basis)[s])[i] = new vector<Polynomial>((u * v * w).gradient());
i++;
// Type 2 -- Part 1
(*cell)[i] = new vector<Polynomial>(term1);
(*(*basis)[s])[i] = new vector<Polynomial>(term1);
i++;
// Type 2 -- Part 2
(*cell)[i] = new vector<Polynomial>(term2);
(*(*basis)[s])[i] = new vector<Polynomial>(term2);
i++;
// Type 3
if(l1 == 1){
(*cell)[i] = new vector<Polynomial>(subGradL1L2VW);
(*(*basis)[s])[i] = new vector<Polynomial>(subGradL1L2VW);
i++;
}
}
}
}
}
// Free Temporary Sapce //
delete[] legendre;
......@@ -346,38 +320,17 @@ TetEdgeBasis::TetEdgeBasis(int order){
}
TetEdgeBasis::~TetEdgeBasis(void){
// Vertex Based //
for(int i = 0; i < nVertex; i++)
delete (*node)[i];
delete node;
// ReferenceSpace //
delete refSpace;
// Edge Based //
for(int c = 0; c < 2; c++){
for(int i = 0; i < nEdge; i++)
delete (*(*edge)[c])[i];
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
delete (*(*basis)[i])[j];
delete (*edge)[c];
delete (*basis)[i];
}
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 basis;
}
delete face;
// Cell Based //
for(int i = 0; i < nCell; i++)
delete (*cell)[i];
delete cell;
}
FunctionSpace/TetEdgeBasis.h
+ 1
− 1
View file @ 708a843e
......@@ -20,7 +20,7 @@ class TetEdgeBasis: public BasisVector{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Tetrahedra of the given order
TetEdgeBasis(int order);
TetEdgeBasis(unsigned int order);
//! Deletes this Basis
//!
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment