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

Nodal Element on Tet (with closure for Edges and Faces) -- Local Functions...

Nodal Element on Tet (with closure for Edges and Faces) -- Local Functions closure for Edges on the fly (no more orientation map) -- Still need Local Functions closure for Faces on the fly AND Face based Dof
parent ace2f4fe
No related branches found
No related tags found
No related merge requests found
......@@ -19,6 +19,7 @@ set(SRC
TriNedelecBasis.cpp
HexNodeBasis.cpp
HexEdgeBasis.cpp
TetNodeBasis.cpp
FunctionSpace.cpp
FunctionSpaceScalar.cpp
......
......@@ -17,6 +17,7 @@ FunctionSpace::~FunctionSpace(void){
delete basis;
// Closure //
/*
map<const MElement*, vector<bool>*>::iterator cIt
= edgeClosure->begin();
map<const MElement*, vector<bool>*>::iterator cStop
......@@ -25,7 +26,7 @@ FunctionSpace::~FunctionSpace(void){
for(; cIt != cStop; cIt++)
delete cIt->second;
delete edgeClosure;
*/
// Dof //
if(dof){
set<const Dof*>::iterator dStop = dof->end();
......@@ -89,7 +90,7 @@ void FunctionSpace::build(const GroupOfElement& goe,
fPerCell = basis->getNCellBased(); // We always got 1 cell
// Build Closure //
buildClosure();
//buildClosure();
// Build Dof //
buildDof();
......
......@@ -9,51 +9,67 @@ FunctionSpaceScalar::~FunctionSpaceScalar(void){
const vector<const Polynomial*> FunctionSpaceScalar::
getLocalFunctions(const MElement& element) const{
// Get Closure //
map<const MElement*, vector<bool>*>::iterator it =
edgeClosure->find(&element);
// Get Basis //
const BasisScalar& basis = getBasis(element);
const vector <const Polynomial*>& node = basis.getNodeFunctions();
const vector<vector<const Polynomial*>*>& edge = basis.getEdgeFunctions();
const vector <const Polynomial*>& cell = basis.getCellFunctions();
if(it == edgeClosure->end())
throw Exception("Element not found for closure");
const unsigned int nFNode = basis.getNVertexBased();
const unsigned int nFEdge = basis.getNEdgeBased();
const unsigned int nFCell = basis.getNCellBased();
vector<bool>* closure = it->second;
// Get Basis //
const vector <const Polynomial*>& node = getBasis(element).getNodeFunctions();
const vector<vector<const Polynomial*>*>& edge = getBasis(element).getEdgeFunctions();
const vector <const Polynomial*>& cell = getBasis(element).getCellFunctions();
// Get Edge Closure //
MElement& eelement = const_cast<MElement&>(element);
const unsigned int nEdge = eelement.getNumEdges();
vector<bool> edgeClosure(nEdge);
// Look Edges
for(unsigned int i = 0; i < nEdge; i++){
MEdge edge = eelement.getEdge(i);
edgeClosure[i] =
edge.getVertex(0)->getNum() > edge.getVertex(1)->getNum();
}
// Get Functions //
vector<const Polynomial*> fun(basis.getSize());
unsigned int i = 0;
const unsigned int nNode = node.size();
const unsigned int nEdge = edge[0]->size();
const unsigned int nCell = cell.size();
vector<const Polynomial*> fun(getBasis(element).getSize());
// Vertex Based //
for(unsigned int j = 0; j < nNode; j++){
// Vertex Based
for(unsigned int j = 0; j < nFNode; j++){
fun[i] = node[j];
i++;
}
// Edge Based //
for(unsigned int j = 0; j < nEdge; j++){
if((*closure)[i])
fun[i] = (*edge[0])[j];
// Edge Based
// Number of basis function *per* edge
// --> should always be an integer !
const unsigned int nFPerEdge = nFEdge / nEdge;
unsigned int fEdge = 0;
for(unsigned int j = 0; j < nFPerEdge; j++){
for(unsigned int k = 0; k < nEdge; k++){
if(edgeClosure[k])
fun[i] = (*edge[0])[fEdge];
else
fun[i] = (*edge[1])[j];
fun[i] = (*edge[1])[fEdge];
fEdge++;
i++;
}
}
// Vertex Based //
for(unsigned int j = 0; j < nCell; j++){
// Vertex Based
for(unsigned int j = 0; j < nFCell; j++){
fun[i] = cell[j];
i++;
}
// Return
// Return //
return fun;
}
......@@ -119,8 +119,8 @@ HexNodeBasis::HexNodeBasis(const int order){
int i = 8;
// Points definig Edges
int edge1[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3};
int edge2[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7};
const unsigned int edge1[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3};
const unsigned int edge2[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7};
for(int l = 1; l < order; l++){
for(int e = 0; e < 12; e++){
......@@ -135,10 +135,10 @@ HexNodeBasis::HexNodeBasis(const int order){
// Face Based (Preliminary) //
// Points definig Faces
int face1[6] = {0, 3, 2, 1, 5, 4};
int face2[6] = {1, 7, 6, 0, 6, 7};
int face3[6] = {2, 6, 5, 4, 7, 3};
int face4[6] = {3, 2, 1, 5, 4, 0};
const unsigned int face1[6] = {0, 3, 2, 1, 5, 4};
const unsigned int face2[6] = {1, 7, 6, 0, 6, 7};
const unsigned int face3[6] = {2, 6, 5, 4, 7, 3};
const unsigned int face4[6] = {3, 2, 1, 5, 4, 0};
// 'Xi' Functions
for(int f = 0; f < 6; f++)
......
#include "TetNodeBasis.h"
#include "Legendre.h"
using namespace std;
TetNodeBasis::TetNodeBasis(unsigned int order){
// 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;
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 unsigned int orderMinus = order - 1;
const unsigned int orderMinusTwo = order - 2;
const unsigned int orderMinusThree = order - 3;
Legendre::legendre(legendre, orderMinus);
Legendre::scaled(sclLegendre, orderMinus);
Legendre::intScaled(intLegendre, order);
// Points definig Edges //
const unsigned int edge0[6] = {0, 1, 2, 2, 3, 3};
const unsigned int edge1[6] = {1, 2, 0, 3, 0, 1};
// Points definig Faces //
const unsigned int face0[4] = {0, 0, 0, 1};
const unsigned int face1[4] = {1, 1, 2, 2};
const unsigned int face2[4] = {3, 2, 3, 3};
// Counter //
unsigned int i;
// Basis //
node = new vector<const Polynomial*>(nVertex);
edge = new vector<vector<const Polynomial*>*>(2);
face = new vector<vector<const Polynomial*>*>(6);
cell = new vector<const Polynomial*>(nCell);
(*edge)[0] = new vector<const Polynomial*>(nEdge);
(*edge)[1] = new vector<const Polynomial*>(nEdge);
(*face)[0] = new vector<const Polynomial*>(nFace);
(*face)[1] = new vector<const Polynomial*>(nFace);
(*face)[2] = new vector<const Polynomial*>(nFace);
(*face)[3] = new vector<const Polynomial*>(nFace);
(*face)[4] = new vector<const Polynomial*>(nFace);
(*face)[5] = new vector<const Polynomial*>(nFace);
// Vertex Based (Lagrange) //
(*node)[0] =
new Polynomial(Polynomial(1, 0, 0, 0) -
Polynomial(1, 1, 0, 0) -
Polynomial(1, 0, 1, 0) -
Polynomial(1, 0, 0, 1));
(*node)[1] =
new Polynomial(Polynomial(1, 1, 0, 0));
(*node)[2] =
new Polynomial(Polynomial(1, 0, 1, 0));
(*node)[3] =
new Polynomial(Polynomial(1, 0, 0, 1));
// Edge Based //
i = 0;
for(unsigned int l = 1; l < order; l++){
for(unsigned int e = 0; e < 6; e++){
(*(*edge)[0])[i] =
new Polynomial(intLegendre[l].compose
(*(*node)[edge0[e]] - *(*node)[edge1[e]],
*(*node)[edge0[e]] + *(*node)[edge1[e]]));
(*(*edge)[1])[i] =
new Polynomial(intLegendre[l].compose
(*(*node)[edge1[e]] - *(*node)[edge0[e]],
*(*node)[edge1[e]] + *(*node)[edge0[e]]));
i++;
}
}
// Face Based //
i = 0;
for(unsigned int l1 = 1; l1 < orderMinus; l1++){
for(unsigned int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
for(unsigned int m = 0; m < 4; m++){
Polynomial sum =
*(*node)[face0[m]] +
*(*node)[face1[m]] +
*(*node)[face2[m]];
(*(*face)[0])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face0[m]] - *(*node)[face1[m]],
*(*node)[face0[m]] + *(*node)[face1[m]]) *
*(*node)[face2[m]] *
sclLegendre[l2].compose
(*(*node)[face2[m]] * 2 - sum, sum));
(*(*face)[1])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face1[m]] - *(*node)[face0[m]],
*(*node)[face1[m]] + *(*node)[face0[m]]) *
*(*node)[face2[m]] *
sclLegendre[l2].compose
(*(*node)[face2[m]] * 2 - sum, sum));
(*(*face)[2])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face1[m]] - *(*node)[face2[m]],
*(*node)[face1[m]] + *(*node)[face2[m]]) *
*(*node)[face0[m]] *
sclLegendre[l2].compose
(*(*node)[face0[m]] * 2 - sum, sum));
(*(*face)[3])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face2[m]] - *(*node)[face1[m]],
*(*node)[face2[m]] + *(*node)[face1[m]]) *
*(*node)[face0[m]] *
sclLegendre[l2].compose
(*(*node)[face0[m]] * 2 - sum, sum));
(*(*face)[4])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face2[m]] - *(*node)[face0[m]],
*(*node)[face2[m]] + *(*node)[face0[m]]) *
*(*node)[face1[m]] *
sclLegendre[l2].compose
(*(*node)[face1[m]] * 2 - sum, sum));
(*(*face)[5])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face0[m]] - *(*node)[face2[m]],
*(*node)[face0[m]] + *(*node)[face2[m]]) *
*(*node)[face1[m]] *
sclLegendre[l2].compose
(*(*node)[face1[m]] * 2 - sum, sum));
i++;
}
}
}
// Cell Based //
Polynomial oneMinusFour = Polynomial(1, 0, 0, 0) - *(*node)[3];
Polynomial twoThreeOneMinusFour = *(*node)[2] * 2 - oneMinusFour;
Polynomial twoFourMinusOne = *(*node)[3] * 2 - Polynomial(1, 0, 0, 0);
Polynomial sub = *(*node)[0] - *(*node)[1];
Polynomial add = *(*node)[0] + *(*node)[1];
i = 0;
for(unsigned int l1 = 1; l1 < orderMinusTwo; l1++){
for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
for(unsigned int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
(*cell)[i] =
new Polynomial(intLegendre[l1].compose(sub, add) *
*(*node)[2] *
sclLegendre[l2].compose(twoThreeOneMinusFour,
oneMinusFour) *
*(*node)[3] *
legendre[l3].compose(twoFourMinusOne));
i++;
}
}
}
// Free Temporary Sapce //
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
}
TetNodeBasis::~TetNodeBasis(void){
// 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 cell;
}
#ifndef _TETNODEBASIS_H_
#define _TETNODEBASIS_H_
#include "BasisScalar.h"
/**
@class TetNodeBasis
@brief A Node Basis for Tetrahedrons
This class can instantiate a Node-Based Basis
(high or low order) for Tetrahedrons.@n
It uses
<a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>
Basis for @em high @em order Polynomial%s generation.@n
*/
class TetNodeBasis: public BasisScalar{
public:
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Tetrahedrons of the given order
TetNodeBasis(unsigned int order);
//! Deletes this Basis
//!
virtual ~TetNodeBasis(void);
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment