From a0020629f044ef97662f911047fb170960450049 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Fri, 19 Oct 2012 14:02:16 +0000
Subject: [PATCH] TeT: Closure OK -- Poisson OK -- Polynomial OK -- \o/ --
Still Need Conv Test
---
FunctionSpace/BasisGenerator.cpp | 15 ++++
FunctionSpace/BasisGenerator.h | 18 +++++
FunctionSpace/FunctionSpace.cpp | 77 +-----------------
FunctionSpace/FunctionSpaceScalar.cpp | 72 +++++++++++++++--
FunctionSpace/HexNodeBasis.cpp | 12 +--
FunctionSpace/TetNodeBasis.cpp | 109 ++++++++++++++------------
6 files changed, 166 insertions(+), 137 deletions(-)
diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 15645c4a5b..2bdbb06439 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -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){
diff --git a/FunctionSpace/BasisGenerator.h b/FunctionSpace/BasisGenerator.h
index 38401febed..2f9ce3cf0a 100644
--- a/FunctionSpace/BasisGenerator.h
+++ b/FunctionSpace/BasisGenerator.h
@@ -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
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index b94195a6ce..28647e5dba 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -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);
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 677267d5ac..afd63a4c9e 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -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();
-
+ 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++;
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index 922d33ba38..43ca97949e 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -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++)
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 05770e8280..efd411af37 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -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) *
--
GitLab