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

TeT: Closure OK -- Poisson OK -- Polynomial OK -- \o/ -- Still Need Conv Test

parent 7ffc2cf0
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
FunctionSpace/BasisGenerator.cpp
+ 15
− 0
View file @ a0020629
......@@ -9,6 +9,8 @@
#include "TriEdgeBasis.h"
#include "TriNedelecBasis.h"
#include "TetNodeBasis.h"
#include "HexNodeBasis.h"
#include "HexEdgeBasis.h"
......@@ -25,6 +27,7 @@ Basis* BasisGenerator::generate(int elementType,
switch(elementType){
case TYPE_TRI: return triGen(basisType, order);
case TYPE_QUA: return quaGen(basisType, order);
case TYPE_TET: return tetGen(basisType, order);
case TYPE_HEX: return hexGen(basisType, order);
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
......@@ -59,6 +62,18 @@ Basis* BasisGenerator::quaGen(int basisType,
}
}
Basis* BasisGenerator::tetGen(int basisType,
int order){
switch(basisType){
case 0: return new TetNodeBasis(order);
case 1: throw Exception("1-form not implemented on Tetrahedrons");
case 2: throw Exception("2-form not implemented on Tetrahedrons");
case 3: throw Exception("3-form not implemented on Tetrahedrons");
default: throw Exception("There is no %d-form", basisType);
}
}
Basis* BasisGenerator::hexGen(int basisType,
int order){
switch(basisType){
......
FunctionSpace/BasisGenerator.h
+ 18
− 0
View file @ a0020629
......@@ -26,6 +26,7 @@ class BasisGenerator{
static Basis* triGen(int basisType, int order);
static Basis* quaGen(int basisType, int order);
static Basis* tetGen(int basisType, int order);
static Basis* hexGen(int basisType, int order);
};
......@@ -100,6 +101,23 @@ class BasisGenerator{
@li @c 3 for 3-Form
**
@fn BasisGenerator::tetGen
@param basisType The Basis type
@param order The order or the requested Basis
This method will @em instanciate the requested Basis,
with a @em Tetrahedron for support
@return Returns a @em pointer to a newly
@em instantiated Basis
@note Basis types are:
@li @c 0 for 0-Form
@li @c 1 for 1-Form
@li @c 2 for 2-Form
@li @c 3 for 3-Form
**
@fn BasisGenerator::hexGen
@param basisType The Basis type
@param order The order or the requested Basis
......
FunctionSpace/FunctionSpace.cpp
+ 2
− 75
View file @ a0020629
......@@ -16,17 +16,6 @@ FunctionSpace::~FunctionSpace(void){
// Basis //
delete basis;
// Closure //
/*
map<const MElement*, vector<bool>*>::iterator cIt
= edgeClosure->begin();
map<const MElement*, vector<bool>*>::iterator cStop
= edgeClosure->end();
for(; cIt != cStop; cIt++)
delete cIt->second;
delete edgeClosure;
*/
// Dof //
if(dof){
set<const Dof*>::iterator dStop = dof->end();
......@@ -89,72 +78,10 @@ void FunctionSpace::build(const GroupOfElement& goe,
fPerCell = basis->getNCellBased(); // We always got 1 cell
// Build Closure //
//buildClosure();
// Build Dof //
buildDof();
}
void FunctionSpace::buildClosure(void){
// Get Elements //
const vector<const MElement*>& element = goe->getAll();
const unsigned int size = element.size();
// Init //
edgeClosure = new map<const MElement*, vector<bool>*>;
// Iterate on elements //
for(unsigned int i = 0; i < size; i++){
// Get Element data
MElement& myElement =
const_cast<MElement&>(*element[i]);
const unsigned int nVertex = myElement.getNumVertices();
const unsigned int nEdge = myElement.getNumEdges();
const unsigned int nFace = myElement.getNumFaces();
const unsigned int nTotVertex = nVertex * fPerVertex;
const unsigned int nTotEdge = nEdge * fPerEdge;
const unsigned int nTotFace = nFace * fPerFace;
const unsigned int nTot = nTotVertex + nTotEdge + nTotFace + fPerCell;
// Closure
vector<bool>* closure = new vector<bool>(nTot);
unsigned int it = 0;
// Closure for vertices
for(unsigned int j = 0; j < nTotVertex; j++, it++)
(*closure)[it] = true;
// Closure for edges
for(unsigned int j = 0; j < fPerEdge; j++){
for(unsigned int k = 0; k < nEdge; k++, it++){
// Orientation
int orientation = mesh->getOrientation(myElement.getEdge(k));
if(orientation == 1)
(*closure)[it] = true;
else
(*closure)[it] = false;
}
}
// Closure for faces
// TODO
// Closure for cells
for(unsigned int j = 0; j < fPerCell; j++, it++)
(*closure)[it] = true;
// Add Closure
edgeClosure->insert
(pair<const MElement*, vector<bool>*>(element[i], closure));
}
}
void FunctionSpace::buildDof(void){
// Get Elements //
unsigned int nElement = goe->getNumber();
......@@ -264,7 +191,7 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
it++;
}
}
/*
// Add Face Based Dof //
for(int i = 0; i < nFFace; i++){
for(int j = 0; j < nFace; j++){
......@@ -272,7 +199,7 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
it++;
}
}
*/
// Add Cell Based Dof //
for(int i = 0; i < nFCell; i++){
myDof[it].setDof(mesh->getGlobalId(element), i);
......
FunctionSpace/FunctionSpaceScalar.cpp
+ 66
− 6
View file @ a0020629
......@@ -13,24 +13,69 @@ getLocalFunctions(const MElement& element) const{
const BasisScalar& basis = getBasis(element);
const vector <const Polynomial*>& node = basis.getNodeFunctions();
const vector<vector<const Polynomial*>*>& edge = basis.getEdgeFunctions();
const vector<vector<const Polynomial*>*>& face = basis.getFaceFunctions();
const vector <const Polynomial*>& cell = basis.getCellFunctions();
const unsigned int nFNode = basis.getNVertexBased();
const unsigned int nFEdge = basis.getNEdgeBased();
const unsigned int nFFace = basis.getNFaceBased();
const unsigned int nFCell = basis.getNCellBased();
// Get not const Element (gmsh problem, not mine !) //
MElement& eelement = const_cast<MElement&>(element);
// 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();
edge.getVertex(0)->getNum() < edge.getVertex(1)->getNum();
}
// Get Face Closure //
const unsigned int nFace = eelement.getNumFaces();
vector<int> faceClosure(nFace);
// Look Faces
for(unsigned int i = 0; i < nFace; i++){
MFace face = eelement.getFace(i);
unsigned int v[3];
v[0] = face.getVertex(0)->getNum();
v[1] = face.getVertex(1)->getNum();
v[2] = face.getVertex(2)->getNum();
bool c[3];
c[0] = v[0] < v[1];
c[1] = v[1] < v[2];
c[2] = v[2] < v[0];
if(c[0])
if(c[1])
faceClosure[i] = 0;
else
if(c[2])
faceClosure[i] = 4;
else
faceClosure[i] = 5;
else
if(c[1])
if(c[2])
faceClosure[i] = 2;
else
faceClosure[i] = 1;
else
faceClosure[i] = 3;
}
......@@ -63,7 +108,22 @@ getLocalFunctions(const MElement& element) const{
}
}
// Vertex Based
// Face Based
// Number of basis function *per* face
// --> should always be an integer !
const unsigned int nFPerFace = nFFace / nFace;
unsigned int fFace = 0;
for(unsigned int j = 0; j < nFPerFace; j++){
for(unsigned int k = 0; k < nFace; k++){
fun[i] = (*face[faceClosure[k]])[fFace];
fFace++;
i++;
}
}
// Cell Based
for(unsigned int j = 0; j < nFCell; j++){
fun[i] = cell[j];
i++;
......
FunctionSpace/HexNodeBasis.cpp
+ 6
− 6
View file @ a0020629
......@@ -119,8 +119,8 @@ HexNodeBasis::HexNodeBasis(const int order){
int i = 8;
// Points definig Edges
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};
const int edge1[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3};
const 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
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};
const int face1[6] = {0, 3, 2, 1, 5, 4};
const int face2[6] = {1, 7, 6, 0, 6, 7};
const int face3[6] = {2, 6, 5, 4, 7, 3};
const int face4[6] = {3, 2, 1, 5, 4, 0};
// 'Xi' Functions
for(int f = 0; f < 6; f++)
......
FunctionSpace/TetNodeBasis.cpp
+ 59
− 50
View file @ a0020629
......@@ -23,25 +23,34 @@ TetNodeBasis::TetNodeBasis(unsigned int 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;
const int orderMinus = order - 1;
const int orderMinusTwo = order - 2;
const 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};
// 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 //
unsigned int i;
int i;
// Basis //
node = new vector<const Polynomial*>(nVertex);
......@@ -81,16 +90,16 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
i = 0;
for(unsigned int l = 1; l < order; l++){
for(unsigned int e = 0; e < 6; e++){
for(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]]));
(*(*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)[edge1[e]] - *(*node)[edge0[e]],
*(*node)[edge1[e]] + *(*node)[edge0[e]]));
(*(*node)[edgeV[e][1]] - *(*node)[edgeV[e][0]],
*(*node)[edgeV[e][1]] + *(*node)[edgeV[e][0]]));
i++;
}
......@@ -100,73 +109,73 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
// 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++){
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)[face0[m]] +
*(*node)[face1[m]] +
*(*node)[face2[m]];
*(*node)[faceV[m][0]] +
*(*node)[faceV[m][1]] +
*(*node)[faceV[m][2]];
(*(*face)[0])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face0[m]] - *(*node)[face1[m]],
*(*node)[face0[m]] + *(*node)[face1[m]]) *
(*(*node)[faceV[m][0]] - *(*node)[faceV[m][1]],
*(*node)[faceV[m][0]] + *(*node)[faceV[m][1]]) *
*(*node)[face2[m]] *
*(*node)[faceV[m][2]] *
sclLegendre[l2].compose
(*(*node)[face2[m]] * 2 - sum, sum));
(*(*node)[faceV[m][2]] * 2 - sum, sum));
(*(*face)[1])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face1[m]] - *(*node)[face0[m]],
*(*node)[face1[m]] + *(*node)[face0[m]]) *
(*(*node)[faceV[m][1]] - *(*node)[faceV[m][0]],
*(*node)[faceV[m][1]] + *(*node)[faceV[m][0]]) *
*(*node)[face2[m]] *
*(*node)[faceV[m][2]] *
sclLegendre[l2].compose
(*(*node)[face2[m]] * 2 - sum, sum));
(*(*node)[faceV[m][2]] * 2 - sum, sum));
(*(*face)[2])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face1[m]] - *(*node)[face2[m]],
*(*node)[face1[m]] + *(*node)[face2[m]]) *
(*(*node)[faceV[m][1]] - *(*node)[faceV[m][2]],
*(*node)[faceV[m][1]] + *(*node)[faceV[m][2]]) *
*(*node)[face0[m]] *
*(*node)[faceV[m][0]] *
sclLegendre[l2].compose
(*(*node)[face0[m]] * 2 - sum, sum));
(*(*node)[faceV[m][0]] * 2 - sum, sum));
(*(*face)[3])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face2[m]] - *(*node)[face1[m]],
*(*node)[face2[m]] + *(*node)[face1[m]]) *
(*(*node)[faceV[m][2]] - *(*node)[faceV[m][1]],
*(*node)[faceV[m][2]] + *(*node)[faceV[m][1]]) *
*(*node)[face0[m]] *
*(*node)[faceV[m][0]] *
sclLegendre[l2].compose
(*(*node)[face0[m]] * 2 - sum, sum));
(*(*node)[faceV[m][0]] * 2 - sum, sum));
(*(*face)[4])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face2[m]] - *(*node)[face0[m]],
*(*node)[face2[m]] + *(*node)[face0[m]]) *
(*(*node)[faceV[m][2]] - *(*node)[faceV[m][0]],
*(*node)[faceV[m][2]] + *(*node)[faceV[m][0]]) *
*(*node)[face1[m]] *
*(*node)[faceV[m][1]] *
sclLegendre[l2].compose
(*(*node)[face1[m]] * 2 - sum, sum));
(*(*node)[faceV[m][1]] * 2 - sum, sum));
(*(*face)[5])[i] =
new Polynomial(intLegendre[l1].compose
(*(*node)[face0[m]] - *(*node)[face2[m]],
*(*node)[face0[m]] + *(*node)[face2[m]]) *
(*(*node)[faceV[m][0]] - *(*node)[faceV[m][2]],
*(*node)[faceV[m][0]] + *(*node)[faceV[m][2]]) *
*(*node)[face1[m]] *
*(*node)[faceV[m][1]] *
sclLegendre[l2].compose
(*(*node)[face1[m]] * 2 - sum, sum));
(*(*node)[faceV[m][1]] * 2 - sum, sum));
i++;
}
......@@ -184,9 +193,9 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
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++){
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++){
(*cell)[i] =
new Polynomial(intLegendre[l1].compose(sub, add) *
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment