From 29db91e265361b88f04c1283f4aca32e3d84ea72 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Wed, 17 Oct 2012 15:27:27 +0000
Subject: [PATCH] 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
---
FunctionSpace/CMakeLists.txt | 3 +-
FunctionSpace/FunctionSpace.cpp | 5 +-
FunctionSpace/FunctionSpaceScalar.cpp | 74 +++++---
FunctionSpace/HexNodeBasis.cpp | 12 +-
FunctionSpace/TetNodeBasis.cpp | 245 ++++++++++++++++++++++++++
FunctionSpace/TetNodeBasis.h | 30 ++++
6 files changed, 331 insertions(+), 38 deletions(-)
create mode 100644 FunctionSpace/TetNodeBasis.cpp
create mode 100644 FunctionSpace/TetNodeBasis.h
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index b31dd7a1b5..a1c28aa15b 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -19,7 +19,8 @@ set(SRC
TriNedelecBasis.cpp
HexNodeBasis.cpp
HexEdgeBasis.cpp
-
+ TetNodeBasis.cpp
+
FunctionSpace.cpp
FunctionSpaceScalar.cpp
FunctionSpaceVector.cpp
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 7cffb23b0d..b94195a6ce 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.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();
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index e5b7887f7c..677267d5ac 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -8,52 +8,68 @@ FunctionSpaceScalar::~FunctionSpaceScalar(void){
const vector<const Polynomial*> FunctionSpaceScalar::
getLocalFunctions(const MElement& element) const{
+
+ // 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();
- // Get Closure //
- map<const MElement*, vector<bool>*>::iterator it =
- edgeClosure->find(&element);
+ const unsigned int nFNode = basis.getNVertexBased();
+ const unsigned int nFEdge = basis.getNEdgeBased();
+ const unsigned int nFCell = basis.getNCellBased();
- if(it == edgeClosure->end())
- throw Exception("Element not found for closure");
- vector<bool>* closure = it->second;
+ // 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 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 Functions //
- 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());
+ vector<const Polynomial*> fun(basis.getSize());
+ unsigned int i = 0;
- // 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];
-
- else
- fun[i] = (*edge[1])[j];
+ // Edge Based
+ // Number of basis function *per* edge
+ // --> should always be an integer !
+ const unsigned int nFPerEdge = nFEdge / nEdge;
+ unsigned int fEdge = 0;
- i++;
+ 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])[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;
}
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index fef2ee1cf4..922d33ba38 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -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++)
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
new file mode 100644
index 0000000000..05770e8280
--- /dev/null
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -0,0 +1,245 @@
+#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;
+}
diff --git a/FunctionSpace/TetNodeBasis.h b/FunctionSpace/TetNodeBasis.h
new file mode 100644
index 0000000000..2083fddf59
--- /dev/null
+++ b/FunctionSpace/TetNodeBasis.h
@@ -0,0 +1,30 @@
+#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
--
GitLab