From 2cdc4b49265a1c8a495a7d16168fa32cc687b5ce Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Mon, 29 Apr 2013 14:38:44 +0000
Subject: [PATCH] Minor Fix: TriNodeBasis face functions -- Add QuadFaces to
ReferenceSpace -- QuadNodeBasis: Edges OK -- QuadFace: Orientation is *WRONG*
---
FunctionSpace/BasisGenerator.cpp | 13 +-
FunctionSpace/CMakeLists.txt | 10 +-
FunctionSpace/LineReferenceSpace.cpp | 6 +-
FunctionSpace/QuadNodeBasis.cpp | 220 ++++++++++++---------------
FunctionSpace/QuadNodeBasis.h | 2 +-
FunctionSpace/QuadReferenceSpace.cpp | 4 +-
FunctionSpace/ReferenceSpace.cpp | 66 +++++++-
FunctionSpace/ReferenceSpace.h | 11 +-
FunctionSpace/TetReferenceSpace.cpp | 4 +-
FunctionSpace/TriNodeBasis.cpp | 14 +-
FunctionSpace/TriReferenceSpace.cpp | 4 +-
11 files changed, 202 insertions(+), 152 deletions(-)
diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 1d8b774cb5..82f64d5041 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -7,14 +7,14 @@
#include "LineNedelecBasis.h"
#include "LineLagrangeBasis.h"
-#include "QuadNodeBasis.h"
-#include "QuadEdgeBasis.h"
-
#include "TriNodeBasis.h"
#include "TriEdgeBasis.h"
#include "TriNedelecBasis.h"
#include "TriLagrangeBasis.h"
+#include "QuadNodeBasis.h"
+#include "QuadEdgeBasis.h"
+
#include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
#include "TetNedelecBasis.h"
@@ -113,8 +113,11 @@ BasisLocal* BasisGenerator::triHierarchicalGen(unsigned int basisType,
BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType,
unsigned int order){
switch(basisType){
- //case 0: return new QuadNodeBasis(order);
- //case 1: return new QuadEdgeBasis(order);
+ case 0: return new QuadNodeBasis(order);
+ case 1:
+ if (order == 0) throw Exception("Nedelec not implemented on Quads");
+ else throw Exception("1-form not implemented on Quads");
+
case 2: throw Exception("2-form not implemented on Quads");
case 3: throw Exception("3-form not implemented on Quads");
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index 193a314648..d93e55d851 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -26,22 +26,22 @@ set(SRC
LineNedelecBasis.cpp
LineLagrangeBasis.cpp
-# QuadNodeBasis.cpp
-# QuadEdgeBasis.cpp
-
TriNodeBasis.cpp
TriEdgeBasis.cpp
TriNedelecBasis.cpp
TriLagrangeBasis.cpp
-# HexNodeBasis.cpp
-# HexEdgeBasis.cpp
+ QuadNodeBasis.cpp
+# QuadEdgeBasis.cpp
TetNodeBasis.cpp
TetEdgeBasis.cpp
TetNedelecBasis.cpp
TetLagrangeBasis.cpp
+# HexNodeBasis.cpp
+# HexEdgeBasis.cpp
+
FunctionSpace.cpp
FunctionSpaceScalar.cpp
FunctionSpaceVector.cpp
diff --git a/FunctionSpace/LineReferenceSpace.cpp b/FunctionSpace/LineReferenceSpace.cpp
index 7ba0eb2d4a..b71ed0cda1 100644
--- a/FunctionSpace/LineReferenceSpace.cpp
+++ b/FunctionSpace/LineReferenceSpace.cpp
@@ -19,8 +19,10 @@ LineReferenceSpace::LineReferenceSpace(void){
nFace = 0;
refFace = NULL;
- // Init All //
- init();
+ // Init All (Use Tri Face -- //
+ // unused -- //
+ // just for interface) //
+ init(0);
}
LineReferenceSpace::~LineReferenceSpace(void){
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 2acd7dc202..4d4b4a6a9b 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -1,131 +1,130 @@
#include "QuadNodeBasis.h"
+#include "QuadReferenceSpace.h"
#include "Legendre.h"
using namespace std;
-QuadNodeBasis::QuadNodeBasis(int order){
+QuadNodeBasis::QuadNodeBasis(unsigned int order){
+ // Reference Space //
+ refSpace = new QuadReferenceSpace;
+ nRefSpace = refSpace->getNPermutation();
+
+ 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;
type = 0;
dim = 2;
- nVertex = 4;
- nEdge = 4 * (order - 1);
- nFace = (order - 1) * (order - 1);
- nCell = 0;
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
-
- size = nVertex + nEdge + nFace + nCell;
-
- // Alloc Temporary Space //
- Polynomial* legendre = new Polynomial[order];
- Polynomial lifting[4];
- Polynomial liftingSub[2][4];
+ nVertex = 4;
+ nEdge = 4 * (order - 1);
+ nFace = (order - 1) * (order - 1);
+ nCell = 0;
+ nFunction = nVertex + nEdge + nFace + nCell;
// Legendre Polynomial //
+ Polynomial* legendre = new Polynomial[order];
Legendre::integrated(legendre, order);
- // Vertices definig Edges & Permutations //
- const int edgeV[2][4][2] =
+ // Lagrange & Lifting //
+ const Polynomial lagrange[4] =
{
- { {0, 1}, {1, 2}, {2, 3}, {3, 0} },
- { {1, 0}, {2, 1}, {3, 2}, {0, 3} }
- };
-
- // Basis //
- node = new vector<Polynomial*>(nVertex);
- edge = new vector<vector<Polynomial*>*>(2);
- face = new vector<vector<Polynomial*>*>(0);
- cell = new vector<Polynomial*>(nCell);
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
- (*edge)[0] = new vector<Polynomial*>(nEdge);
- (*edge)[1] = new vector<Polynomial*>(nEdge);
+ Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
+ Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0))),
- // Lifting //
- lifting[0] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
-
- lifting[1] =
- (Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
-
- lifting[2] =
- (Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 1, 0));
-
- lifting[3] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
- (Polynomial(1, 0, 1, 0));
-
- // Lifting Sub //
- for(int e = 0; e < 4; e++){
- liftingSub[0][e] =
- lifting[edgeV[0][e][0]] -
- lifting[edgeV[0][e][1]];
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)))
+ };
- liftingSub[1][e] =
- lifting[edgeV[1][e][0]] -
- lifting[edgeV[1][e][1]];
- }
+ const Polynomial lifting[4] =
+ {
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
+ Polynomial((Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
- // Vertex Based (Lagrange) //
- (*node)[0] =
- new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
+ Polynomial((Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0))),
- (*node)[1] =
- new Polynomial((Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
+ Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)))
+ };
- (*node)[2] =
- new Polynomial((Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)));
+ // Basis //
+ basis = new Polynomial**[nRefSpace];
- (*node)[3] =
- new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)));
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ basis[s] = new Polynomial*[nFunction];
+ // Vertex Based //
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ basis[s][0] = new Polynomial(lagrange[0]);
+ basis[s][1] = new Polynomial(lagrange[1]);
+ basis[s][2] = new Polynomial(lagrange[2]);
+ basis[s][3] = new Polynomial(lagrange[3]);
+ }
// Edge Based //
- for(int c = 0; c < 2; c++){
- unsigned int i = 0;
-
- for(int l = 1; l < order; l++){
- for(int e = 0; e < 4; e++){
- (*(*edge)[c])[i] =
- new Polynomial(legendre[l].compose(liftingSub[c][e]) *
- (*(*node)[edgeV[c][e][0]] + *(*node)[edgeV[c][e][1]]));
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nVertex;
+
+ for(unsigned int e = 0; e < 4; e++){
+ for(unsigned int l = 1; l < order; l++){
+ basis[s][i] =
+ new Polynomial(legendre[l].compose(lifting[(*(*edgeV[s])[e])[1]] -
+ lifting[(*(*edgeV[s])[e])[0]])
+ *
+ (lagrange[(*(*edgeV[s])[e])[0]] +
+ lagrange[(*(*edgeV[s])[e])[1]]));
i++;
}
}
}
+ // Face Based //
- // Cell Based //
- Polynomial px = Polynomial(2, 1, 0, 0);
- Polynomial py = Polynomial(2, 0, 1, 0);
+ // NB: We use (*(*faceV[s])[f])[]
+ // where f = 0, because triangles
+ // have only ONE face: the face '0'
- px = px - Polynomial(1, 0, 0, 0);
- py = py - Polynomial(1, 0, 0, 0);
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nVertex + nEdge;
- unsigned int i = 0;
+ for(unsigned int l1 = 1; l1 < order; l1++){
+ for(unsigned int l2 = 1; l2 < order; l2++){
+ basis[s][i] =
+ new Polynomial(legendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
+ lifting[(*(*faceV[s])[0])[1]])
+ *
- for(int l1 = 1; l1 < order; l1++){
- for(int l2 = 1; l2 < order; l2++){
- (*cell)[i] =
- new Polynomial(legendre[l1].compose(px) * legendre[l2].compose(py));
+ legendre[l2].compose(lifting[(*(*faceV[s])[0])[0]] -
+ lifting[(*(*faceV[s])[0])[3]]));
- i++;
+ i++;
+ }
}
}
+ cout << (*(*faceV[1])[0])[0] << ", "
+ << (*(*faceV[1])[0])[1] << endl
+
+ << (*(*faceV[1])[0])[0] << ", "
+ << (*(*faceV[1])[0])[3] << endl
+
+ << (*(*faceV[1])[0])[2] << endl;
// Mapping to Gmsh Quad //
// x = (u + 1) / 2
@@ -140,47 +139,30 @@ QuadNodeBasis::QuadNodeBasis(int order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0));
- for(int i = 0; i < nVertex; i++)
- *(*node)[i] = (*node)[i]->compose(mapX, mapY);
-
- for(int i = 0; i < nEdgeClosure; i++)
- for(int j = 0; j < nEdge; j++)
- *(*(*edge)[i])[j] = (*(*edge)[i])[j]->compose(mapX, mapY);
-
- for(int i = 0; i < nCell; i++)
- *(*cell)[i] = (*cell)[i]->compose(mapX, mapY);
-
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(unsigned int i = 0; i < nFunction; i++){
+ Polynomial* tmp;
+ tmp = basis[s][i];
+ basis[s][i] = new Polynomial(tmp->compose(mapX, mapY));
+ delete tmp;
+ }
+ }
// Free Temporary Sapce //
delete[] legendre;
}
QuadNodeBasis::~QuadNodeBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
+ // ReferenceSpace //
+ delete refSpace;
- delete node;
-
-
- // 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 //
- delete face;
-
-
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
-
- delete cell;
+ delete[] basis;
}
diff --git a/FunctionSpace/QuadNodeBasis.h b/FunctionSpace/QuadNodeBasis.h
index 15f6205d20..7bc159aa1d 100644
--- a/FunctionSpace/QuadNodeBasis.h
+++ b/FunctionSpace/QuadNodeBasis.h
@@ -24,7 +24,7 @@ class QuadNodeBasis: public BasisHierarchical0From{
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Quads of the given order
- QuadNodeBasis(int order);
+ QuadNodeBasis(unsigned int order);
//! @return Deletes this Basis
//!
diff --git a/FunctionSpace/QuadReferenceSpace.cpp b/FunctionSpace/QuadReferenceSpace.cpp
index ecec74d447..5bc382730c 100644
--- a/FunctionSpace/QuadReferenceSpace.cpp
+++ b/FunctionSpace/QuadReferenceSpace.cpp
@@ -28,8 +28,8 @@ QuadReferenceSpace::QuadReferenceSpace(void){
refFace[0][2] = 2;
refFace[0][3] = 3;
- // Init All //
- init();
+ // Init All (Quad Face) //
+ init(1);
}
QuadReferenceSpace::~QuadReferenceSpace(void){
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index 7b9ef5967b..a2e89c8ade 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -36,7 +36,12 @@ ReferenceSpace::ReferenceSpace(void){
// Defining Ref Edge and Face in //
// Dervived Class //
- // And CALL INIT() //
+ // And CALL init(faceType) //
+ // //
+ // faceType is: //
+ // * 0 if triangular faces //
+ // * 1 if quadragular faces //
+ // * 3 if NOT IMPLEMENTED //
}
ReferenceSpace::~ReferenceSpace(void){
@@ -65,7 +70,13 @@ ReferenceSpace::~ReferenceSpace(void){
delete face;
}
-void ReferenceSpace::init(void){
+void ReferenceSpace::init(unsigned int faceType){
+ if(faceType > 3)
+ throw Exception("ReferenceSpace: unknown face type");
+
+ if(faceType == 3)
+ throw Exception("ReferenceSpace: prism not implemented");
+
// Init Root //
nPerm = 0;
nextLeafId = 0;
@@ -96,7 +107,13 @@ void ReferenceSpace::init(void){
// Get Edges & Faces //
getEdge();
- getFace();
+
+ switch(faceType){
+ case 0: getTriFace(); break;
+ case 1: getQuaFace(); break;
+
+ default: throw Exception("ReferenceSpace: impossible error");
+ }
}
void ReferenceSpace::populate(node* pTreeRoot){
@@ -188,7 +205,7 @@ void ReferenceSpace::getEdge(void){
}
}
-void ReferenceSpace::getFace(void){
+void ReferenceSpace::getTriFace(void){
vector<const vector<unsigned int>*>* tmp;
face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
@@ -204,6 +221,23 @@ void ReferenceSpace::getFace(void){
}
}
+void ReferenceSpace::getQuaFace(void){
+ vector<const vector<unsigned int>*>* tmp;
+ face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
+
+ for(unsigned int p = 0; p < nPerm; p++){
+ tmp = new vector<const vector<unsigned int>*>(nFace);
+
+ for(unsigned int i = 0; i < nFace; i++)
+ (*tmp)[i] = inOrder(p,
+ refFace[i][0],
+ refFace[i][1],
+ refFace[i][2],
+ refFace[i][3]);
+ (*face)[p] = tmp;
+ }
+}
+
const vector<unsigned int>* ReferenceSpace::
inOrder(unsigned int permutation,
unsigned int a,
@@ -249,6 +283,30 @@ inOrder(unsigned int permutation,
return inorder;
}
+const std::vector<unsigned int>* ReferenceSpace::
+inOrder(unsigned int permutation,
+ unsigned int a,
+ unsigned int b,
+ unsigned int c,
+ unsigned int d){
+
+ unsigned int v;
+ unsigned int k = 0;
+ vector<unsigned int>* inorder =
+ new vector<unsigned int>(4);
+
+ for(unsigned int i = 0; i < nVertex; i++){
+ v = perm[permutation][i];
+
+ if(v == a || v == b || v == c || v == d){
+ (*inorder)[k] = v;
+ k++;
+ }
+ }
+
+ return inorder;
+}
+
unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
// Const_Cast //
MElement& element = const_cast<MElement&>(elem);
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index 5fd8258b44..cf8f78a578 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -80,7 +80,7 @@ class ReferenceSpace{
protected:
ReferenceSpace(void);
- void init(void);
+ void init(unsigned int faceType);
void populate(node* pTreeRoot);
void destroy(node* node);
@@ -89,7 +89,8 @@ class ReferenceSpace{
void unconnect(void); // Unconnects leafs marked before
void getEdge(void);
- void getFace(void);
+ void getTriFace(void);
+ void getQuaFace(void);
const std::vector<unsigned int>* inOrder(unsigned int permutation,
unsigned int a,
@@ -100,6 +101,12 @@ class ReferenceSpace{
unsigned int b,
unsigned int c);
+ const std::vector<unsigned int>* inOrder(unsigned int permutation,
+ unsigned int a,
+ unsigned int b,
+ unsigned int c,
+ unsigned int d);
+
static bool sortPredicate(const std::pair<unsigned int, MVertex*>& a,
const std::pair<unsigned int, MVertex*>& b);
diff --git a/FunctionSpace/TetReferenceSpace.cpp b/FunctionSpace/TetReferenceSpace.cpp
index b591b39023..4fcc1f1db7 100644
--- a/FunctionSpace/TetReferenceSpace.cpp
+++ b/FunctionSpace/TetReferenceSpace.cpp
@@ -29,8 +29,8 @@ TetReferenceSpace::TetReferenceSpace(void){
refFace[i][2] = MTetrahedron::faces_tetra(i, 2);
}
- // Init All //
- init();
+ // Init All (Tri Face) //
+ init(0);
}
TetReferenceSpace::~TetReferenceSpace(void){
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 3a0cefd789..b6f1604b7d 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -27,13 +27,12 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
- // Alloc Some Space //
+ // Legendre Polynomial //
const int orderMinus = order - 1;
Polynomial* legendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order];
- // Legendre Polynomial //
Legendre::legendre(legendre, orderMinus);
Legendre::intScaled(intLegendre, order);
@@ -66,7 +65,7 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex;
- for(int e = 0; e < 3; e++){
+ for(unsigned int e = 0; e < 3; e++){
for(unsigned int l = 1; l < order; l++){
basis[s][i] =
new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
@@ -84,10 +83,6 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
// NB: We use (*(*faceV[s])[f])[]
// where f = 0, because triangles
// have only ONE face: the face '0'
-
- // TO CHECK: Are Triangles face matching tets ?
-
- const Polynomial p = (lagrange[2] * 2) - Polynomial(1, 0, 0, 0);
const int orderMinusTwo = order - 2;
for(unsigned int s = 0; s < nRefSpace; s++){
@@ -103,8 +98,11 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
lagrange[(*(*faceV[s])[0])[1]])
*
- legendre[l2].compose(lagrange[(*(*faceV[s])[0])[2]])
+ legendre[l2].compose((lagrange[(*(*faceV[s])[0])[2]] * 2)
+ -
+ Polynomial(1, 0, 0, 0))
*
+
lagrange[(*(*faceV[s])[0])[2]]);
i++;
}
diff --git a/FunctionSpace/TriReferenceSpace.cpp b/FunctionSpace/TriReferenceSpace.cpp
index a0c9b9f46d..bb9f03a24a 100644
--- a/FunctionSpace/TriReferenceSpace.cpp
+++ b/FunctionSpace/TriReferenceSpace.cpp
@@ -27,8 +27,8 @@ TriReferenceSpace::TriReferenceSpace(void){
refFace[0][1] = 1;
refFace[0][2] = 2;
- // Init All //
- init();
+ // Init All (Tri Face) //
+ init(0);
}
TriReferenceSpace::~TriReferenceSpace(void){
--
GitLab