diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 49de2e09c269bfbf8a61f56050ad5005b3452c7b..b9f920ba7fbd95f6e4584d64f43d0571f8d308e3 100644
--- a/FunctionSpace/Basis.h
+++ b/FunctionSpace/Basis.h
@@ -3,7 +3,6 @@
#include <string>
#include "MElement.h"
-#include "ReferenceSpace.h"
/**
@interface Basis
@@ -20,38 +19,32 @@
The i-th row of these matrices is always refering to the
i-th function of the basis.
- Depending on the nature of the returned value
- (scalar or vector), the columns are organized
- diferently.
+ Depending on the nature of the returned value (scalar or vector),
+ the columns are organized diferently.
For scalar values, we have:
@li The j-th column of the i-th row is
the evaluation of the i-th function at the j-th point
For vectorial values, we have:
- @li The j-th column of the i-th row is
- the first coordinate of
+ @li The j-th column of the i-th row is the first coordinate of
the evaluation of the i-th function at the 3 x j-th point
- @li The (j-th + 1) column of the i-th row is
- the second coordinate of
+ @li The (j-th + 1) column of the i-th row is the second coordinate of
the evaluation of the i-th function at the 3 x j-th point
- @li The (j-th + 2) column of the i-th row is
- the third coordinate of
+ @li The (j-th + 2) column of the i-th row is the third coordinate of
the evaluation of the i-th function at the 3 x j-th point
*/
class Basis{
protected:
- ReferenceSpace* refSpace;
- size_t nRefSpace;
-
bool scalar;
bool local;
size_t order;
size_t type;
+ size_t form;
size_t dim;
size_t nVertex;
@@ -72,6 +65,7 @@ class Basis{
// Type of Basis //
size_t getOrder(void) const;
size_t getType(void) const;
+ size_t getForm(void) const;
size_t getDim(void) const;
// Number of Functions //
@@ -81,9 +75,6 @@ class Basis{
size_t getNCellBased(void) const;
size_t getNFunction(void) const;
- // Reference Element //
- const ReferenceSpace& getReferenceSpace(void) const;
-
// Direct Access to Evaluated Functions //
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
@@ -98,8 +89,10 @@ class Basis{
double u, double v, double w) const = 0;
// Precompute Functions //
- virtual void preEvaluateFunctions(const fullMatrix<double>& point) const = 0;
- virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const = 0;
+ virtual
+ void preEvaluateFunctions(const fullMatrix<double>& point) const = 0;
+ virtual
+ void preEvaluateDerivatives(const fullMatrix<double>& point) const = 0;
// Access to Precomputed Functions //
virtual const fullMatrix<double>&
@@ -140,8 +133,8 @@ class Basis{
@li true, if this is a scalar Basis
@li false, if this is a vectorial Basis
- Scalar basis are sets of Polynomial%s, and
- Vectorial basis are sets of Vector%s of Polynomial%s
+ Scalar basis are sets of Polynomial%s,
+ and Vectorial basis are sets of Vector%s of Polynomial%s
**
@fn Basis::isLocal
@@ -155,7 +148,19 @@ class Basis{
**
@fn Basis::getType
- @return Returns the type of the Basis:
+ @return Returns the type of the Basis (coherent with gmsh element types):
+ @li 1 for Points
+ @li 2 for Lines
+ @li 3 for Triangles
+ @li 4 for Quadrangles
+ @li 5 for Tetrahedra
+ @li 6 for Pyramids
+ @li 7 for Prisms
+ @li 8 for Hexahedra
+ **
+
+ @fn Basis::getForm
+ @return Returns the diferential form of the Basis:
@li 0 for 0-form
@li 1 for 1-form
@li 2 for 2-form
@@ -187,10 +192,6 @@ class Basis{
(or Vector%s of Polynomial%s) Functions in this Basis
**
- @fn Basis::getReferenceSpace
- @return Returns the ReferenceSpace associated to this basis
- **
-
@fn Basis::getFunctions(fullMatrix<double>&, const MElement&, double, double, double) const
@param retValues An allocated matrix
@param element A MElement
@@ -213,6 +214,17 @@ class Basis{
and for the given orientation
**
+ @fn Basis::getDerivative(fullMatrix<double>&, const MElement&, double, double, double) const
+ @param retValues An allocated matrix
+ @param element A MElement
+ @param u A u coordinate in the reference space of this Basis
+ @param v A v coordinate in the reference space of this Basis
+ @param w A w coordinate in the reference space of this Basis
+ @return The given matrix is populated with the evaluation
+ of the derivative of every basis function at the given coordinates,
+ and for the orientation of the given element
+ **
+
@fn Basis::preEvaluateFunctions
@param point A Matrix with points coordinate
(each line is a point and got 3 coordinates, i.e. 3 rows)
@@ -243,8 +255,8 @@ class Basis{
@return Returns a Matrix with the PreEvaluated basis functions derivatives
(see Basis::preEvaluateDerivatives()), with the given orientation
- If no PreEvaluation of the gradient has been done before calling this function,
- an Exception is thrown
+ If no PreEvaluation of the gradient has been done
+ before calling this function, an Exception is thrown
**
@fn Basis::getPreEvaluatedFunctions(const MElement&) const
@@ -279,6 +291,10 @@ inline size_t Basis::getType(void) const{
return type;
}
+inline size_t Basis::getForm(void) const{
+ return form;
+}
+
inline size_t Basis::getDim(void) const{
return dim;
}
@@ -303,8 +319,4 @@ inline size_t Basis::getNFunction(void) const{
return nFunction;
}
-inline const ReferenceSpace& Basis::getReferenceSpace(void) const{
- return *refSpace;
-}
-
#endif
diff --git a/FunctionSpace/BasisHierarchical0Form.cpp b/FunctionSpace/BasisHierarchical0Form.cpp
index 351001b9e25e6192ba9869077605074d7fee3251..df97eebade7b5e38655d354e7198356ddf5c84cd 100644
--- a/FunctionSpace/BasisHierarchical0Form.cpp
+++ b/FunctionSpace/BasisHierarchical0Form.cpp
@@ -1,5 +1,7 @@
#include <sstream>
#include "Exception.h"
+#include "ReferenceSpaceManager.h"
+
#include "BasisHierarchical0Form.h"
using namespace std;
@@ -8,6 +10,9 @@ BasisHierarchical0Form::BasisHierarchical0Form(void){
// Scalar Basis ? //
scalar = true;
+ // 0-Form //
+ form = 0;
+
// Grad Basis //
hasGrad = false;
grad = NULL;
@@ -21,9 +26,11 @@ BasisHierarchical0Form::BasisHierarchical0Form(void){
}
BasisHierarchical0Form::~BasisHierarchical0Form(void){
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Grad Basis //
if(hasGrad){
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete grad[i][j];
@@ -35,14 +42,14 @@ BasisHierarchical0Form::~BasisHierarchical0Form(void){
// PreEvaluation //
if(preEvaluated){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
if(preEvaluatedGrad){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction;
@@ -55,7 +62,7 @@ getFunctions(fullMatrix<double>& retValues,
double u, double v, double w) const{
// Define Orientation //
- const size_t orientation = refSpace->getReferenceSpace(element);
+ const size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Matrix //
for(size_t i = 0; i < nFunction; i++)
@@ -80,7 +87,7 @@ void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues,
getGrad();
// Define Orientation //
- const size_t orientation = refSpace->getReferenceSpace(element);
+ const size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Matrix //
for(size_t i = 0; i < nFunction; i++){
@@ -92,62 +99,62 @@ void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues,
void BasisHierarchical0Form::
preEvaluateFunctions(const fullMatrix<double>& point) const{
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Delete if older //
if(preEvaluated){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
// Alloc //
- const size_t nPoint = point.size1();
- preEvaluatedFunction = new fullMatrix<double>*[nRefSpace];
+ const size_t nPoint = point.size1();
+ preEvaluatedFunction = new fullMatrix<double>*[nOrientation];
- for(size_t i = 0; i < nRefSpace; i++)
- preEvaluatedFunction[i] =
- new fullMatrix<double>(nFunction, nPoint);
+ for(size_t i = 0; i < nOrientation; i++)
+ preEvaluatedFunction[i] = new fullMatrix<double>(nFunction, nPoint);
// Fill Matrix //
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
for(size_t j = 0; j < nFunction; j++)
for(size_t k = 0; k < nPoint; k++)
- (*preEvaluatedFunction[i])(j, k) =
- basis[i][j]->at(point(k, 0),
- point(k, 1),
- point(k, 2));
-
+ (*preEvaluatedFunction[i])(j, k) = basis[i][j]->at(point(k, 0),
+ point(k, 1),
+ point(k, 2));
// PreEvaluated //
preEvaluated = true;
}
void BasisHierarchical0Form::
preEvaluateDerivatives(const fullMatrix<double>& point) const{
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Build Grad //
if(!hasGrad)
getGrad();
// Delete if older //
if(preEvaluatedGrad){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction;
}
// Alloc //
- const size_t nPoint = point.size1();
- const size_t nPoint3 = nPoint * 3;
- preEvaluatedGradFunction = new fullMatrix<double>*[nRefSpace];
+ const size_t nPoint = point.size1();
+ const size_t nPoint3 = nPoint * 3;
+ preEvaluatedGradFunction = new fullMatrix<double>*[nOrientation];
- for(size_t i = 0; i < nRefSpace; i++)
- preEvaluatedGradFunction[i] =
- new fullMatrix<double>(nFunction, nPoint3);
+ for(size_t i = 0; i < nOrientation; i++)
+ preEvaluatedGradFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++){
for(size_t k = 0; k < nPoint; k++){
tmp = Polynomial::at(*grad[i][j],
@@ -168,18 +175,21 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedFunctions(const MElement& element) const{
- return getPreEvaluatedFunctions(refSpace->getReferenceSpace(element));
+ return
+ getPreEvaluatedFunctions(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedDerivatives(const MElement& element) const{
- return getPreEvaluatedDerivatives(refSpace->getReferenceSpace(element));
+ return
+ getPreEvaluatedDerivatives(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedFunctions(size_t orientation) const{
if(!preEvaluated)
- throw Exception("getPreEvaluatedFunction: function has not been preEvaluated");
+ throw Exception
+ ("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[orientation];
}
@@ -187,23 +197,25 @@ getPreEvaluatedFunctions(size_t orientation) const{
const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedDerivatives(size_t orientation) const{
if(!preEvaluatedGrad)
- throw Exception("getPreEvaluatedDerivative: gradient has not been preEvaluated");
+ throw Exception
+ ("getPreEvaluatedDerivative: gradient has not been preEvaluated");
return *preEvaluatedGradFunction[orientation];
}
void BasisHierarchical0Form::getGrad(void) const{
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Alloc //
- grad = new vector<Polynomial>**[nRefSpace];
+ grad = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
grad[s] = new vector<Polynomial>*[nFunction];
// Grad //
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
for(size_t f = 0 ; f < nFunction; f++)
- grad[s][f] =
- new vector<Polynomial>(basis[s][f]->gradient());
+ grad[s][f] = new vector<Polynomial>(basis[s][f]->gradient());
// Has Grad //
hasGrad = true;
diff --git a/FunctionSpace/BasisHierarchical1Form.cpp b/FunctionSpace/BasisHierarchical1Form.cpp
index 783bc19f30cf34ad71600084ecce4b99d5a0be3f..06b51302b2924f761ba944d1f82d5f741b665fea 100644
--- a/FunctionSpace/BasisHierarchical1Form.cpp
+++ b/FunctionSpace/BasisHierarchical1Form.cpp
@@ -1,5 +1,7 @@
#include <sstream>
#include "Exception.h"
+#include "ReferenceSpaceManager.h"
+
#include "BasisHierarchical1Form.h"
using namespace std;
@@ -8,6 +10,9 @@ BasisHierarchical1Form::BasisHierarchical1Form(void){
// Scalar Basis ?//
scalar = false;
+ // 1-Form //
+ form = 1;
+
// Curl Basis //
hasCurl = false;
curl = NULL;
@@ -21,9 +26,11 @@ BasisHierarchical1Form::BasisHierarchical1Form(void){
}
BasisHierarchical1Form::~BasisHierarchical1Form(void){
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Curl Basis //
if(hasCurl){
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete curl[i][j];
@@ -35,14 +42,14 @@ BasisHierarchical1Form::~BasisHierarchical1Form(void){
// PreEvaluation //
if(preEvaluated){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
if(preEvaluatedCurl){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction;
@@ -55,12 +62,11 @@ getFunctions(fullMatrix<double>& retValues,
double u, double v, double w) const{
// Define Orientation //
- size_t orientation = refSpace->getReferenceSpace(element);
+ size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Vector //
for(size_t i = 0; i < nFunction; i++){
- fullVector<double> eval =
- Polynomial::at(*basis[orientation][i], u, v, w);
+ fullVector<double> eval = Polynomial::at(*basis[orientation][i], u, v, w);
retValues(i, 0) = eval(0);
retValues(i, 1) = eval(1);
@@ -75,8 +81,7 @@ getFunctions(fullMatrix<double>& retValues,
// Fill Vector //
for(size_t i = 0; i < nFunction; i++){
- fullVector<double> eval =
- Polynomial::at(*basis[orientation][i], u, v, w);
+ fullVector<double> eval = Polynomial::at(*basis[orientation][i], u, v, w);
retValues(i, 0) = eval(0);
retValues(i, 1) = eval(1);
@@ -92,9 +97,11 @@ void BasisHierarchical1Form::getDerivative(fullMatrix<double>& retValues,
void BasisHierarchical1Form::
preEvaluateFunctions(const fullMatrix<double>& point) const{
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Delete if older //
if(preEvaluated){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
@@ -103,16 +110,15 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{
// Alloc //
const size_t nPoint = point.size1();
const size_t nPoint3 = nPoint * 3;
- preEvaluatedFunction = new fullMatrix<double>*[nRefSpace];
+ preEvaluatedFunction = new fullMatrix<double>*[nOrientation];
- for(size_t i = 0; i < nRefSpace; i++)
- preEvaluatedFunction[i] =
- new fullMatrix<double>(nFunction, nPoint3);
+ for(size_t i = 0; i < nOrientation; i++)
+ preEvaluatedFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++){
for(size_t k = 0; k < nPoint; k++){
tmp = Polynomial::at(*basis[i][j],
@@ -133,31 +139,32 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{
void BasisHierarchical1Form::
preEvaluateDerivatives(const fullMatrix<double>& point) const{
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Build Curl //
if(!hasCurl)
getCurl();
// Delete if older //
if(preEvaluatedCurl){
- for(size_t i = 0; i < nRefSpace; i++)
+ for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction;
}
// Alloc //
- const size_t nPoint = point.size1();
- const size_t nPoint3 = nPoint * 3;
- preEvaluatedCurlFunction = new fullMatrix<double>*[nRefSpace];
+ const size_t nPoint = point.size1();
+ const size_t nPoint3 = nPoint * 3;
+ preEvaluatedCurlFunction = new fullMatrix<double>*[nOrientation];
- for(size_t i = 0; i < nRefSpace; i++)
- preEvaluatedCurlFunction[i] =
- new fullMatrix<double>(nFunction, nPoint3);
+ for(size_t i = 0; i < nOrientation; i++)
+ preEvaluatedCurlFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++){
for(size_t k = 0; k < nPoint; k++){
tmp = Polynomial::at(*curl[i][j],
@@ -178,18 +185,21 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedFunctions(const MElement& element) const{
- return getPreEvaluatedFunctions(refSpace->getReferenceSpace(element));
+ return
+ getPreEvaluatedFunctions(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedDerivatives(const MElement& element) const{
- return getPreEvaluatedDerivatives(refSpace->getReferenceSpace(element));
+ return
+ getPreEvaluatedDerivatives(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedFunctions(size_t orientation) const{
if(!preEvaluated)
- throw Exception("getPreEvaluatedFunction: function has not been preEvaluated");
+ throw
+ Exception("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[orientation];
}
@@ -197,23 +207,25 @@ getPreEvaluatedFunctions(size_t orientation) const{
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedDerivatives(size_t orientation) const{
if(!preEvaluatedCurl)
- throw Exception("getPreEvaluatedDerivative: curl has not been preEvaluated");
+ throw
+ Exception("getPreEvaluatedDerivative: curl has not been preEvaluated");
return *preEvaluatedCurlFunction[orientation];
}
void BasisHierarchical1Form::getCurl(void) const{
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Alloc //
- curl = new vector<Polynomial>**[nRefSpace];
+ curl = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
curl[s] = new vector<Polynomial>*[nFunction];
// Curl //
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
for(size_t f = 0 ; f < nFunction; f++)
- curl[s][f] =
- new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
+ curl[s][f] = new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
// Has Curl //
hasCurl = true;
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index d5e9c2ff9c4981565f6183ea6ded6f66a394cc8d..ba4cc231b87ed6eae109f4547569489f167407a0 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -1,9 +1,11 @@
+#include "ReferenceSpaceManager.h"
#include "BasisLagrange.h"
using namespace std;
BasisLagrange::BasisLagrange(void){
scalar = true;
+ form = 0;
preEvaluated = false;
preEvaluatedGrad = false;
@@ -38,7 +40,7 @@ getFunctions(fullMatrix<double>& retValues,
retValues = tmp.transpose();
// Permute retValues, accordingly to ReferenceSpace
- permutation(refSpace->getReferenceSpace(element), retValues);
+ permutation(ReferenceSpaceManager::getOrientation(element), retValues);
}
void BasisLagrange::
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index a54b1450ec77c1500bf5b13b95c5afe3436f8fc3..a356168c42347bb235cd95f9aca1d99af863e1fd 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -18,6 +18,8 @@ set(SRC
PyrReferenceSpace.cpp
PriReferenceSpace.cpp
+ ReferenceSpaceManager.cpp
+
Basis.cpp
BasisLocal.cpp
BasisGenerator.cpp
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index c6aebe8293e0ed2ef776b2240bc851e5015ad5de..6b2550f5c513e335742985bc5f140e66ab82d6ef 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -1,9 +1,12 @@
#include <sstream>
-#include "FunctionSpace.h"
+
+#include "ReferenceSpaceManager.h"
#include "BasisGenerator.h"
#include "ElementType.h"
#include "Exception.h"
+#include "FunctionSpace.h"
+
using namespace std;
FunctionSpace::FunctionSpace(void){
@@ -35,9 +38,6 @@ FunctionSpace::~FunctionSpace(void){
// Element To GoD //
if(eToGod)
delete eToGod;
-
- // Unorient GroupOfElement //
- goe->unoriented();
}
void FunctionSpace::build(GroupOfElement& goe,
@@ -47,9 +47,6 @@ void FunctionSpace::build(GroupOfElement& goe,
this->goe = &goe;
this->mesh = &(goe.getMesh());
- // Orient All Elements //
- goe.orientAllElements(basis); // NOT SEXY: TO BE REMOVED
-
// Get Geo Data (WARNING HOMOGENE MESH REQUIRED)//
const MElement& element = goe.get(0);
MElement& myElement =
@@ -227,7 +224,7 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
// Create New Element With Permuted Vertices //
// Permutation
const vector<size_t>& vPerm =
- (*this->basis)[0]->getReferenceSpace().getNodeIndexFromABCtoUVW(elem);
+ ReferenceSpaceManager::getNodeIndexFromABCtoUVW(elem);
// Permuted Vertices
const size_t nVertex = element.getNumPrimaryVertices();
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index f04f030ebf1284a48471819fcf6b4c7ad6edbf9a..4808c53a72ee20dab8501d1733704c7a704e21ed 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -38,9 +38,7 @@ interpolateDerivativeInABC(const MElement& element,
double abc[3]) const{
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
- (*basis)[0]->getReferenceSpace().getJacobian(element,
- abc[0], abc[1], abc[2],
- invJac);
+ ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
invJac.invertInPlace();
// Get Basis Functions //
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 95dd1a4b91decdc56fc7a1752f0d64b5aa12822b..66fd1264ecdb55735ac0a79f9c009a22ed61d732 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -1,6 +1,7 @@
#ifndef _FUNCTIONSPACESCALAR_H_
#define _FUNCTIONSPACESCALAR_H_
+#include "ReferenceSpaceManager.h"
#include "FunctionSpace.h"
/**
@@ -110,9 +111,8 @@ interpolate(const MElement& element,
// Get ABC Space coordinate //
double abc[3];
- (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
- xyz(0), xyz(1), xyz(2),
- abc);
+ ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
+
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
@@ -124,9 +124,8 @@ interpolateInRefSpace(const MElement& element,
// Get ABC Space coordinate //
double abc[3];
- (*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element,
- uvw(0), uvw(1), uvw(2),
- abc);
+ ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
+
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
@@ -138,9 +137,8 @@ interpolateDerivative(const MElement& element,
// Get ABC Space coordinate //
double abc[3];
- (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
- xyz(0), xyz(1), xyz(2),
- abc);
+ ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
+
// Interpolate in ABC //
return interpolateDerivativeInABC(element, coef, abc);
}
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index 3f1e5b69ae1636642c74c8f7d77035b415346687..0cd0f96a954441fbcca4f923e24b1fc9d6b9ac0e 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -18,9 +18,7 @@ interpolateInABC(const MElement& element,
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
- (*basis)[0]->getReferenceSpace().getJacobian(element,
- abc[0], abc[1], abc[2],
- invJac);
+ ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
invJac.invertInPlace();
// Get Basis Functions //
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 6515f0045648f3c960a7648529bf792a5e14075f..2379d95cd84e5a8052e4037304b217db4a3fc1f1 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -2,6 +2,7 @@
#define _FUNCTIONSPACEVECTOR_H_
#include "fullMatrix.h"
+#include "FunctionSpaceScalar.h"
#include "FunctionSpace.h"
/**
@@ -103,9 +104,8 @@ interpolate(const MElement& element,
// Get ABC Space coordinate //
double abc[3];
- (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
- xyz(0), xyz(1), xyz(2),
- abc);
+ ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
+
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
@@ -117,9 +117,8 @@ interpolateInRefSpace(const MElement& element,
// Get ABC Space coordinate //
double abc[3];
- (*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element,
- uvw(0), uvw(1), uvw(2),
- abc);
+ ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
+
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
diff --git a/FunctionSpace/HexLagrangeBasis.cpp b/FunctionSpace/HexLagrangeBasis.cpp
index 7fc92e2ceed2c1f2cc3cd550caa293888a307395..71d2be0295b72247333c7a260756a2a76c28455b 100644
--- a/FunctionSpace/HexLagrangeBasis.cpp
+++ b/FunctionSpace/HexLagrangeBasis.cpp
@@ -1,7 +1,8 @@
-#include "HexLagrangeBasis.h"
-#include "HexReferenceSpace.h"
-#include "pointsGenerators.h"
#include "ElementType.h"
+#include "GmshDefines.h"
+#include "pointsGenerators.h"
+
+#include "HexLagrangeBasis.h"
HexLagrangeBasis::HexLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1
@@ -11,7 +12,7 @@ HexLagrangeBasis::HexLagrangeBasis(size_t order){
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_HEX;
dim = 3;
nVertex = 8;
@@ -25,14 +26,9 @@ HexLagrangeBasis::HexLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsHexahedron(order, false));
-
- // Reference Space //
- refSpace = new HexReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
}
HexLagrangeBasis::~HexLagrangeBasis(void){
delete lBasis;
delete lPoint;
- delete refSpace;
}
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index 6e197f4dd48538221aa8f04fa60a5f53c6f77580..ecbb3bebc100d4dbc02d91f898c9cdf0bf2a035f 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -1,24 +1,16 @@
-#include "HexNodeBasis.h"
-#include "HexReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "HexNodeBasis.h"
using namespace std;
HexNodeBasis::HexNodeBasis(size_t order){
- // Reference Space //
- refSpace = new HexReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
- const vector<vector<vector<size_t> > >&
- faceIdx = refSpace->getFaceNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_HEX;
dim = 3;
nVertex = 8;
@@ -27,6 +19,15 @@ HexNodeBasis::HexNodeBasis(size_t order){
nCell = (order - 1) * (order - 1) * (order - 1);
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = ReferenceSpaceManager::getFaceNodeIndex(type);
+
// Legendre Polynomial //
Polynomial* legendre = new Polynomial[order];
Legendre::integrated(legendre, order);
@@ -103,13 +104,13 @@ HexNodeBasis::HexNodeBasis(size_t order){
};
// Basis //
- basis = new Polynomial**[nRefSpace];
+ basis = new Polynomial**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new Polynomial*[nFunction];
// Vertex Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
basis[s][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]);
@@ -121,7 +122,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
}
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex;
for(size_t e = 0; e < 12; e++){
@@ -139,7 +140,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
}
// Face Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex + nEdge;
for(size_t f = 0; f < 6; f++){
@@ -173,7 +174,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
py = py - Polynomial(1, 0, 0, 0);
pz = pz - Polynomial(1, 0, 0, 0);
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex + nEdge + nFace;
for(size_t l1 = 1; l1 < order; l1++){
@@ -206,7 +207,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
Polynomial mapZ(Polynomial(0.5, 0, 0, 1) +
Polynomial(0.5, 0, 0, 0));
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
for(size_t i = 0; i < nFunction; i++){
Polynomial* tmp;
tmp = basis[s][i];
@@ -220,11 +221,10 @@ HexNodeBasis::HexNodeBasis(size_t order){
}
HexNodeBasis::~HexNodeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp
index f3ff36015b8d56e4e327bed96f3c87581cbc645a..1801458f47a67af3de74ceb2caa1d1dd0f397fce 100644
--- a/FunctionSpace/LineEdgeBasis.cpp
+++ b/FunctionSpace/LineEdgeBasis.cpp
@@ -1,21 +1,16 @@
-#include "LineEdgeBasis.h"
-#include "LineReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "LineEdgeBasis.h"
using namespace std;
LineEdgeBasis::LineEdgeBasis(size_t order){
- // Reference Space //
- refSpace = new LineReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 1;
+ type = TYPE_LIN;
dim = 1;
nVertex = 0;
@@ -24,6 +19,12 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
// Legendre Polynomial //
const size_t orderPlus = order + 1;
Polynomial* intLegendre = new Polynomial[orderPlus];
@@ -41,13 +42,13 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
};
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = 0;
for(size_t l = 0; l < orderPlus; l++){
@@ -87,11 +88,10 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
}
LineEdgeBasis::~LineEdgeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/LineLagrangeBasis.cpp b/FunctionSpace/LineLagrangeBasis.cpp
index 6842133c20b99439eaf7c7f4758f70cc0494237b..30792d36fc3ba4d871a08055cde167e01b393aee 100644
--- a/FunctionSpace/LineLagrangeBasis.cpp
+++ b/FunctionSpace/LineLagrangeBasis.cpp
@@ -1,7 +1,8 @@
-#include "LineLagrangeBasis.h"
-#include "LineReferenceSpace.h"
-#include "pointsGenerators.h"
#include "ElementType.h"
+#include "GmshDefines.h"
+#include "pointsGenerators.h"
+
+#include "LineLagrangeBasis.h"
LineLagrangeBasis::LineLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1
@@ -11,7 +12,7 @@ LineLagrangeBasis::LineLagrangeBasis(size_t order){
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_LIN;
dim = 1;
nVertex = 2;
@@ -25,14 +26,9 @@ LineLagrangeBasis::LineLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsLine(order));
-
- // Reference Space //
- refSpace = new LineReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
}
LineLagrangeBasis::~LineLagrangeBasis(void){
delete lBasis;
delete lPoint;
- delete refSpace;
}
diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp
index ee1850ff1ccafa0f31d68c21312f38f1c4543ed8..3f8e55bd4c04599e9a2e58ffc29edccfd1dccb58 100644
--- a/FunctionSpace/LineNedelecBasis.cpp
+++ b/FunctionSpace/LineNedelecBasis.cpp
@@ -1,20 +1,15 @@
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
#include "LineNedelecBasis.h"
-#include "LineReferenceSpace.h"
using namespace std;
LineNedelecBasis::LineNedelecBasis(void){
- // Reference Space //
- refSpace = new LineReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
// Set Basis Type //
order = 0;
- type = 1;
+ type = TYPE_LIN;
dim = 1;
nVertex = 0;
@@ -23,6 +18,12 @@ LineNedelecBasis::LineNedelecBasis(void){
nCell = 0;
nFunction = 1;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
// Lagrange Polynomial //
const Polynomial lagrange[2] =
{
@@ -34,13 +35,13 @@ LineNedelecBasis::LineNedelecBasis(void){
};
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
vector<Polynomial> tmp1 = lagrange[edgeIdx[s][0][1]].gradient();
vector<Polynomial> tmp2 = lagrange[edgeIdx[s][0][0]].gradient();
@@ -61,11 +62,10 @@ LineNedelecBasis::LineNedelecBasis(void){
}
LineNedelecBasis::~LineNedelecBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index d4e9d5cb9915b935dc45cc474e0989aeb0aaac3c..337801c7c8d66aaa89a8b6dbbcd3b2f0fbc85a76 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -1,21 +1,16 @@
-#include "LineNodeBasis.h"
-#include "LineReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "LineNodeBasis.h"
using namespace std;
LineNodeBasis::LineNodeBasis(size_t order){
- // Reference Space //
- refSpace = new LineReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_LIN;
dim = 1;
nVertex = 2;
@@ -24,6 +19,12 @@ LineNodeBasis::LineNodeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
// Legendre Polynomial //
Polynomial* intLegendre = new Polynomial[order];
Legendre::integrated(intLegendre, order);
@@ -39,19 +40,19 @@ LineNodeBasis::LineNodeBasis(size_t order){
};
// Basis //
- basis = new Polynomial**[nRefSpace];
+ basis = new Polynomial**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new Polynomial*[nFunction];
// Vertex Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
basis[s][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]);
}
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex;
for(size_t l = 1; l < order; l++){
@@ -68,11 +69,10 @@ LineNodeBasis::LineNodeBasis(size_t order){
}
LineNodeBasis::~LineNodeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 9fd2d3b753e82187ed6789c7503824fbbfc7ec17..337678818f3cc5afe80f4dbe1d0d2ce153ab6a44 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -1,24 +1,16 @@
-#include "QuadEdgeBasis.h"
-#include "QuadReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "QuadEdgeBasis.h"
using namespace std;
QuadEdgeBasis::QuadEdgeBasis(size_t order){
- // Reference Space //
- refSpace = new QuadReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
- const vector<vector<vector<size_t> > >&
- faceIdx = refSpace->getFaceNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 1;
+ type = TYPE_QUA;
dim = 2;
nVertex = 0;
@@ -27,6 +19,16 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = ReferenceSpaceManager::getFaceNodeIndex(type);
+
+
// Legendre Polynomial //
const size_t orderPlus = order + 1;
@@ -68,13 +70,13 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
};
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = 0;
for(size_t e = 0; e < 4; e++){
@@ -115,7 +117,7 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
// where f = 0, because triangles
// have only ONE face: the face '0'
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nEdge;
// Type 1
@@ -220,7 +222,7 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0));
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
for(size_t i = 0; i < nFunction; i++){
vector<Polynomial>* old;
vector<Polynomial> nxt(3);
@@ -241,11 +243,10 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
}
QuadEdgeBasis::~QuadEdgeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/QuadLagrangeBasis.cpp b/FunctionSpace/QuadLagrangeBasis.cpp
index 31efcf5555f255b9e2b78642c7b26ac531bba467..30974d5ec4d7af7d227f52d37de3686ff54a9ffc 100644
--- a/FunctionSpace/QuadLagrangeBasis.cpp
+++ b/FunctionSpace/QuadLagrangeBasis.cpp
@@ -1,7 +1,7 @@
-#include "QuadLagrangeBasis.h"
-#include "QuadReferenceSpace.h"
-#include "pointsGenerators.h"
#include "ElementType.h"
+#include "GmshDefines.h"
+#include "pointsGenerators.h"
+#include "QuadLagrangeBasis.h"
QuadLagrangeBasis::QuadLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1
@@ -11,7 +11,7 @@ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_QUA;
dim = 2;
nVertex = 4;
@@ -25,14 +25,9 @@ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsQuadrangle(order, false));
-
- // Reference Space //
- refSpace = new QuadReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
}
QuadLagrangeBasis::~QuadLagrangeBasis(void){
delete lBasis;
delete lPoint;
- delete refSpace;
}
diff --git a/FunctionSpace/QuadNedelecBasis.cpp b/FunctionSpace/QuadNedelecBasis.cpp
index b42175f3cce9085d23bc79ed9d44623d1bce1e82..2bf00eb8c8628667e7aa197e3ed2760889114c5d 100644
--- a/FunctionSpace/QuadNedelecBasis.cpp
+++ b/FunctionSpace/QuadNedelecBasis.cpp
@@ -1,21 +1,16 @@
-#include "QuadNedelecBasis.h"
-#include "QuadReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "QuadNedelecBasis.h"
using namespace std;
QuadNedelecBasis::QuadNedelecBasis(void){
- // Reference Space //
- refSpace = new QuadReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
// Set Basis Type //
order = 0;
- type = 1;
+ type = TYPE_QUA;
dim = 2;
nVertex = 0;
@@ -24,6 +19,12 @@ QuadNedelecBasis::QuadNedelecBasis(void){
nCell = 0;
nFunction = 4;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
// Lagrange & Lifting //
const Polynomial lagrange[4] =
{
@@ -56,13 +57,13 @@ QuadNedelecBasis::QuadNedelecBasis(void){
};
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
for(size_t e = 0; e < 4; e++){
Polynomial lambda = (lagrange[edgeIdx[s][e][0]] +
@@ -91,7 +92,7 @@ QuadNedelecBasis::QuadNedelecBasis(void){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0));
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
for(size_t i = 0; i < nFunction; i++){
vector<Polynomial>* old;
vector<Polynomial> nxt(3);
@@ -108,11 +109,10 @@ QuadNedelecBasis::QuadNedelecBasis(void){
}
QuadNedelecBasis::~QuadNedelecBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 68c01d21e92ced83b4f539d53a5fc23993fe4aa3..00b254051e38ffbcda69517df281ad104b169620 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -1,24 +1,16 @@
-#include "QuadNodeBasis.h"
-#include "QuadReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "QuadNodeBasis.h"
using namespace std;
QuadNodeBasis::QuadNodeBasis(size_t order){
- // Reference Space //
- refSpace = new QuadReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
- const vector<vector<vector<size_t> > >&
- faceIdx = refSpace->getFaceNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_QUA;
dim = 2;
nVertex = 4;
@@ -27,6 +19,15 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = ReferenceSpaceManager::getFaceNodeIndex(type);
+
// Legendre Polynomial //
Polynomial* legendre = new Polynomial[order];
Legendre::integrated(legendre, order);
@@ -63,13 +64,13 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
};
// Basis //
- basis = new Polynomial**[nRefSpace];
+ basis = new Polynomial**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new Polynomial*[nFunction];
// Vertex Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
basis[s][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]);
@@ -77,7 +78,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
}
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex;
for(size_t e = 0; e < 4; e++){
@@ -100,7 +101,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
// where f = 0, because triangles
// have only ONE face: the face '0'
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex + nEdge;
for(size_t l1 = 1; l1 < order; l1++){
@@ -131,7 +132,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0));
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
for(size_t i = 0; i < nFunction; i++){
Polynomial* tmp;
tmp = basis[s][i];
@@ -145,11 +146,10 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
}
QuadNodeBasis::~QuadNodeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index 2e31620d439b9f634162b6fe663ac9e5fd58fedb..b926f8f38f9eecd5249a10d7f63f51cf8ce1897a 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -92,8 +92,8 @@ class ReferenceSpace{
ReferenceSpace(const std::string& path);
virtual ~ReferenceSpace(void);
- size_t getNReferenceSpace(void) const;
- size_t getReferenceSpace(const MElement& element) const;
+ size_t getNOrientation(void) const;
+ size_t getOrientation(const MElement& element) const;
const std::vector<std::vector<std::vector<size_t> > >&
getEdgeNodeIndex(void) const;
@@ -348,10 +348,15 @@ class ReferenceSpace{
// Inline Functions //
//////////////////////
-inline size_t ReferenceSpace::getNReferenceSpace(void) const{
+inline size_t ReferenceSpace::getNOrientation(void) const{
return refSpaceNodeId.size();
}
+inline
+size_t ReferenceSpace::getOrientation(const MElement& element) const{
+ return pTree->getTagFromPermutation(getPermutationIdx(element));
+}
+
inline
const std::vector<std::vector<std::vector<size_t> > >&
ReferenceSpace::getEdgeNodeIndex(void) const{
@@ -370,11 +375,6 @@ ReferenceSpace::getNodeIndexFromABCtoUVW(const MElement& element) const{
return ABCtoUVWIndex[getPermutationIdx(element)];
}
-inline
-size_t ReferenceSpace::getReferenceSpace(const MElement& element) const{
- return pTree->getTagFromPermutation(getPermutationIdx(element));
-}
-
inline
bool
ReferenceSpace::sortPredicate(const std::pair<size_t, size_t>& a,
diff --git a/FunctionSpace/ReferenceSpaceManager.cpp b/FunctionSpace/ReferenceSpaceManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3023c2e36bae195cc5ec545fe798496a9ff9479d
--- /dev/null
+++ b/FunctionSpace/ReferenceSpaceManager.cpp
@@ -0,0 +1,43 @@
+#include "Exception.h"
+#include "ReferenceSpaceManager.h"
+
+using namespace std;
+
+const size_t ReferenceSpaceManager::nSpace = 9;
+vector<ReferenceSpace*> ReferenceSpaceManager::refSpace(nSpace, NULL);
+
+ReferenceSpaceManager::ReferenceSpaceManager(void){
+}
+
+ReferenceSpaceManager::~ReferenceSpaceManager(void){
+}
+
+void ReferenceSpaceManager::init(int elementType){
+ // Warning: no elementType == 0
+ // --> refSpace[0] will be alwys null --> not realy important
+
+ switch(elementType){
+ case 1: throw Exception("ReferenceSpace for Points not implemented");
+
+ case 2: refSpace[elementType] = new LineReferenceSpace; break;
+ case 3: refSpace[elementType] = new TriReferenceSpace; break;
+ case 4: refSpace[elementType] = new QuadReferenceSpace; break;
+ case 5: refSpace[elementType] = new TetReferenceSpace; break;
+ case 6: refSpace[elementType] = new PyrReferenceSpace; break;
+ case 7: refSpace[elementType] = new PriReferenceSpace; break;
+ case 8: refSpace[elementType] = new HexReferenceSpace; break;
+
+ case 9: throw Exception("ReferenceSpace for POLYG not implemented");
+ case 10: throw Exception("ReferenceSpace for POLYH not implemented");
+ case 11: throw Exception("ReferenceSpace for XFEM not implemented");
+ }
+}
+
+void ReferenceSpaceManager::clear(void){
+ for(size_t i = 0; i < nSpace; i++)
+ if(refSpace[i])
+ delete refSpace[i];
+
+ for(size_t i = 0; i < nSpace; i++)
+ refSpace[i] = NULL;
+}
diff --git a/FunctionSpace/ReferenceSpaceManager.h b/FunctionSpace/ReferenceSpaceManager.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ec773bd0be2b56c08bdba6883764039cd50833a
--- /dev/null
+++ b/FunctionSpace/ReferenceSpaceManager.h
@@ -0,0 +1,178 @@
+#ifndef _REFERENCESPACEMANAGER_H_
+#define _REFERENCESPACEMANAGER_H_
+
+#include "ReferenceSpace.h"
+
+#include "LineReferenceSpace.h"
+#include "TriReferenceSpace.h"
+#include "QuadReferenceSpace.h"
+#include "TetReferenceSpace.h"
+#include "PyrReferenceSpace.h"
+#include "PriReferenceSpace.h"
+#include "HexReferenceSpace.h"
+
+/**
+ @class ReferenceSpaceManager
+ @biref A way to handel ReferenceSpace%s
+
+ This class implements class method to handel and access ReferenceSpace%s
+ */
+
+class ReferenceSpaceManager{
+ private:
+ static const size_t nSpace;
+ static std::vector<ReferenceSpace*> refSpace;
+
+ public:
+ ReferenceSpaceManager(void);
+ ~ReferenceSpaceManager(void);
+
+ static void clear(void);
+ static const ReferenceSpace& getReferenceSpace(int elementType);
+
+ static size_t getNOrientation(int elementType);
+ static size_t getOrientation(const MElement& element);
+
+ static const std::vector<std::vector<std::vector<size_t> > >&
+ getEdgeNodeIndex(int elementType);
+
+ static const std::vector<std::vector<std::vector<size_t> > >&
+ getFaceNodeIndex(int elementType);
+
+ static const std::vector<size_t>&
+ getNodeIndexFromABCtoUVW(const MElement& element);
+
+ static void mapFromABCtoUVW(const MElement& element,
+ double a, double b, double c, double uvw[3]);
+
+ static void mapFromABCtoXYZ(const MElement& element,
+ double a, double b, double c, double xyz[3]);
+
+ static void mapFromUVWtoABC(const MElement& element,
+ double u, double v, double w, double abc[3]);
+
+ static void mapFromXYZtoABC(const MElement& element,
+ double x, double y, double z, double abc[3]);
+
+ static double getJacobian(const MElement& element,
+ double a, double b, double c,
+ fullMatrix<double>& jac);
+ private:
+ static void init(int elementType);
+};
+
+////////////////////
+// Inline Methods //
+////////////////////
+
+inline
+const ReferenceSpace& ReferenceSpaceManager::getReferenceSpace(int elementType){
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return *refSpace[elementType];
+}
+
+inline
+size_t ReferenceSpaceManager::getNOrientation(int elementType){
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->getNOrientation();
+}
+
+inline
+size_t ReferenceSpaceManager::getOrientation(const MElement& element){
+ const int elementType = element.getType();
+
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->getOrientation(element);
+}
+
+inline
+const std::vector<std::vector<std::vector<size_t> > >&
+ReferenceSpaceManager::getEdgeNodeIndex(int elementType){
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->getEdgeNodeIndex();
+}
+
+inline
+const std::vector<std::vector<std::vector<size_t> > >&
+ReferenceSpaceManager::getFaceNodeIndex(int elementType){
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->getFaceNodeIndex();
+}
+
+inline
+const std::vector<size_t>&
+ReferenceSpaceManager::getNodeIndexFromABCtoUVW(const MElement& element){
+ const int elementType = element.getType();
+
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->getNodeIndexFromABCtoUVW(element);
+}
+
+inline void ReferenceSpaceManager::mapFromABCtoUVW(const MElement& element,
+ double a, double b, double c,
+ double uvw[3]){
+ const int elementType = element.getType();
+
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->mapFromABCtoUVW(element, a, b, c, uvw);
+}
+
+inline void ReferenceSpaceManager::mapFromABCtoXYZ(const MElement& element,
+ double a, double b, double c,
+ double xyz[3]){
+ const int elementType = element.getType();
+
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->mapFromABCtoXYZ(element, a, b, c, xyz);
+}
+
+inline void ReferenceSpaceManager::mapFromUVWtoABC(const MElement& element,
+ double u, double v, double w,
+ double abc[3]){
+ const int elementType = element.getType();
+
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->mapFromUVWtoABC(element, u, v, w, abc);
+}
+
+inline void ReferenceSpaceManager::mapFromXYZtoABC(const MElement& element,
+ double x, double y, double z,
+ double abc[3]){
+ const int elementType = element.getType();
+
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->mapFromXYZtoABC(element, x, y, z, abc);
+}
+
+inline double ReferenceSpaceManager::getJacobian(const MElement& element,
+ double a, double b, double c,
+ fullMatrix<double>& jac){
+ const int elementType = element.getType();
+
+ if(!refSpace[elementType])
+ init(elementType);
+
+ return refSpace[elementType]->getJacobian(element, a, b, c, jac);
+}
+
+#endif
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 782593ee5ec190683338ba8e006f3f5c3dd6e043..71c9565c09aff1e6776aa2b11c850e78d1e7ee8f 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -1,24 +1,16 @@
-#include "TetEdgeBasis.h"
-#include "TetReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "TetEdgeBasis.h"
using namespace std;
TetEdgeBasis::TetEdgeBasis(size_t order){
- // Reference Space //
- refSpace = new TetReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
- const vector<vector<vector<size_t> > >&
- faceIdx = refSpace->getFaceNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 1;
+ type = TYPE_TET;
dim = 3;
nVertex = 0;
@@ -27,6 +19,15 @@ TetEdgeBasis::TetEdgeBasis(size_t order){
nCell = (order + 1) * (order - 1) * (order - 2) / 2;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = ReferenceSpaceManager::getFaceNodeIndex(type);
+
// Alloc Temporary Space //
const int orderPlus = order + 1;
const int orderMinus = order - 1;
@@ -58,13 +59,13 @@ TetEdgeBasis::TetEdgeBasis(size_t order){
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = 0;
for(int e = 0; e < 6; e++){
@@ -107,7 +108,7 @@ TetEdgeBasis::TetEdgeBasis(size_t order){
// Face Based //
// TO CHECK: Are Triangles face matching tets ?
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nEdge;
for(int f = 0; f < 4; f++){
@@ -199,7 +200,7 @@ TetEdgeBasis::TetEdgeBasis(size_t order){
// Cell Based //
const Polynomial one(1, 0, 0, 0);
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nEdge + nFace;
for(int l1 = 1; l1 < orderMinus; l1++){
@@ -324,11 +325,10 @@ TetEdgeBasis::TetEdgeBasis(size_t order){
}
TetEdgeBasis::~TetEdgeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/TetLagrangeBasis.cpp b/FunctionSpace/TetLagrangeBasis.cpp
index 527c6c33db9df391db9582ba4594d0950dff47b7..05ba166ed8ca03e54350a41d57414d0fec9d8b94 100644
--- a/FunctionSpace/TetLagrangeBasis.cpp
+++ b/FunctionSpace/TetLagrangeBasis.cpp
@@ -1,7 +1,8 @@
-#include "TetLagrangeBasis.h"
-#include "TetReferenceSpace.h"
-#include "pointsGenerators.h"
#include "ElementType.h"
+#include "GmshDefines.h"
+#include "pointsGenerators.h"
+
+#include "TetLagrangeBasis.h"
TetLagrangeBasis::TetLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1
@@ -11,7 +12,7 @@ TetLagrangeBasis::TetLagrangeBasis(size_t order){
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_TET;
dim = 3;
nVertex = 4;
@@ -25,14 +26,9 @@ TetLagrangeBasis::TetLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsTetrahedron(order, false));
-
- // Reference Space //
- refSpace = new TetReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
}
TetLagrangeBasis::~TetLagrangeBasis(void){
delete lBasis;
delete lPoint;
- delete refSpace;
}
diff --git a/FunctionSpace/TetNedelecBasis.cpp b/FunctionSpace/TetNedelecBasis.cpp
index a459231ecbdb39231b7496b3cad991f79d61a210..7a027ef689a4da06f8fc85ef774db93f8798dc0b 100644
--- a/FunctionSpace/TetNedelecBasis.cpp
+++ b/FunctionSpace/TetNedelecBasis.cpp
@@ -1,20 +1,15 @@
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
#include "TetNedelecBasis.h"
-#include "TetReferenceSpace.h"
using namespace std;
TetNedelecBasis::TetNedelecBasis(void){
- // Reference Space //
- refSpace = new TetReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
// Set Basis Type //
this->order = 0;
- type = 1;
+ type = TYPE_TET;
dim = 3;
nVertex = 0;
@@ -23,6 +18,12 @@ TetNedelecBasis::TetNedelecBasis(void){
nCell = 0;
nFunction = 6;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
// Lagrange Polynomial //
const Polynomial lagrange[4] =
{
@@ -40,13 +41,13 @@ TetNedelecBasis::TetNedelecBasis(void){
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
for(size_t e = 0; e < 6; e++){
vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
vector<Polynomial> tmp2 = lagrange[edgeIdx[s][e][0]].gradient();
@@ -69,11 +70,10 @@ TetNedelecBasis::TetNedelecBasis(void){
}
TetNedelecBasis::~TetNedelecBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 71fa5d79236cd7ac93cb787eadace310cc1c279c..9080b20e3944996e077d69420bf99bf5447c976d 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -1,24 +1,16 @@
-#include "TetNodeBasis.h"
-#include "TetReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "TetNodeBasis.h"
using namespace std;
TetNodeBasis::TetNodeBasis(size_t order){
- // Reference Space //
- refSpace = new TetReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
- const vector<vector<vector<size_t> > >&
- faceIdx = refSpace->getFaceNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_TET;
dim = 3;
nVertex = 4;
@@ -27,6 +19,15 @@ TetNodeBasis::TetNodeBasis(size_t order){
nCell = (order - 1) * (order - 2) * (order - 3) / 6;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = ReferenceSpaceManager::getFaceNodeIndex(type);
+
// Alloc Temporary Space //
const int orderMinus = order - 1;
const int orderMinusTwo = order - 2;
@@ -58,13 +59,13 @@ TetNodeBasis::TetNodeBasis(size_t order){
// Basis //
- basis = new Polynomial**[nRefSpace];
+ basis = new Polynomial**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new Polynomial*[nFunction];
// Vertex Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
basis[s][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]);
@@ -72,7 +73,7 @@ TetNodeBasis::TetNodeBasis(size_t order){
}
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex;
for(int e = 0; e < 6; e++){
@@ -89,7 +90,7 @@ TetNodeBasis::TetNodeBasis(size_t order){
}
}
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex + nEdge;
for(int f = 0; f < 4; f++){
@@ -123,12 +124,13 @@ TetNodeBasis::TetNodeBasis(size_t order){
// Cell Based //
const Polynomial oneMinusFour = Polynomial(1, 0, 0, 0) - lagrange[3];
const Polynomial twoThreeOneMinusFour = lagrange[2] * 2 - oneMinusFour;
- const Polynomial twoFourMinusOne = lagrange[3] * 2 - Polynomial(1, 0, 0, 0);
+ const Polynomial twoFourMinusOne =
+ lagrange[3] * 2 - Polynomial(1, 0, 0, 0);
const Polynomial sub = lagrange[0] - lagrange[1];
const Polynomial add = lagrange[0] + lagrange[1];
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex + nEdge + nFace;
for(int l1 = 1; l1 < orderMinusTwo; l1++){
@@ -154,11 +156,10 @@ TetNodeBasis::TetNodeBasis(size_t order){
}
TetNodeBasis::~TetNodeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 5a3a5883cc7f28201b18cbc6f803614e53d64f1e..b3005310f80ba7417fef4197a9b7307aa5055697 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -1,24 +1,16 @@
-#include "TriEdgeBasis.h"
-#include "TriReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "TriEdgeBasis.h"
using namespace std;
TriEdgeBasis::TriEdgeBasis(size_t order){
- // Reference Space //
- refSpace = new TriReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
- const vector<vector<vector<size_t> > >&
- faceIdx = refSpace->getFaceNodeIndex();
-
// Set Basis Type //
this->order = order;
- type = 1;
+ type = TYPE_TRI;
dim = 2;
nVertex = 0;
@@ -27,6 +19,15 @@ TriEdgeBasis::TriEdgeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = ReferenceSpaceManager::getFaceNodeIndex(type);
+
// Alloc Some Space //
const int orderPlus = order + 1;
const int orderMinus = order - 1;
@@ -51,13 +52,13 @@ TriEdgeBasis::TriEdgeBasis(size_t order){
};
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = 0;
for(int e = 0; e < 3; e++){
@@ -106,19 +107,19 @@ TriEdgeBasis::TriEdgeBasis(size_t order){
// TO CHECK: Are Triangles face matching tets ?
// Alloc Temp
- Polynomial** u = new Polynomial*[nRefSpace];
- Polynomial** v = new Polynomial*[nRefSpace];
- Polynomial* p = new Polynomial[nRefSpace];
- vector<Polynomial>** subGrad = new vector<Polynomial>*[nRefSpace];
+ Polynomial** u = new Polynomial*[nOrientation];
+ Polynomial** v = new Polynomial*[nOrientation];
+ Polynomial* p = new Polynomial[nOrientation];
+ vector<Polynomial>** subGrad = new vector<Polynomial>*[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
u[s] = new Polynomial[orderPlus];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
v[s] = new Polynomial[orderPlus];
// Preliminaries
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
p[s] = lagrange[faceIdx[s][0][2]] * 2 - Polynomial(1, 0, 0, 0);
// Polynomial u(x) & v(x)
@@ -153,7 +154,7 @@ TriEdgeBasis::TriEdgeBasis(size_t order){
}
// Face Basis
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nEdge;
// Type 1
@@ -219,13 +220,13 @@ TriEdgeBasis::TriEdgeBasis(size_t order){
}
// Free Temporary Sapce //
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
delete[] u[s];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
delete[] v[s];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
delete subGrad[s];
delete[] legendre;
@@ -237,11 +238,10 @@ TriEdgeBasis::TriEdgeBasis(size_t order){
}
TriEdgeBasis::~TriEdgeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index e95fcdd1ae489efe170ef28fd3b53a7520b502c4..96d16c7f1e294091a8ae8c954f820e1d0b46ac1b 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -1,7 +1,8 @@
-#include "TriLagrangeBasis.h"
-#include "TriReferenceSpace.h"
-#include "pointsGenerators.h"
#include "ElementType.h"
+#include "GmshDefines.h"
+#include "pointsGenerators.h"
+
+#include "TriLagrangeBasis.h"
TriLagrangeBasis::TriLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1
@@ -11,7 +12,7 @@ TriLagrangeBasis::TriLagrangeBasis(size_t order){
// Set Basis Type //
this->order = order;
- type = 0;
+ type = TYPE_TRI;
dim = 2;
nVertex = 3;
@@ -25,14 +26,9 @@ TriLagrangeBasis::TriLagrangeBasis(size_t order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsTriangle(order, false));
-
- // Reference Space //
- refSpace = new TriReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
}
TriLagrangeBasis::~TriLagrangeBasis(void){
delete lBasis;
delete lPoint;
- delete refSpace;
}
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index c8ac22bfb4a66d60eeadb681b65c30b1345f77ba..db87fb845d618c393b57afc34b8d7b47a8229fee 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -1,20 +1,15 @@
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
#include "TriNedelecBasis.h"
-#include "TriReferenceSpace.h"
using namespace std;
TriNedelecBasis::TriNedelecBasis(void){
- // Reference Space //
- refSpace = new TriReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
// Set Basis Type //
order = 0;
- type = 1;
+ type = TYPE_TRI;
dim = 2;
nVertex = 0;
@@ -23,6 +18,12 @@ TriNedelecBasis::TriNedelecBasis(void){
nCell = 0;
nFunction = 3;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
// Lagrange //
const Polynomial lagrange[3] =
{
@@ -36,13 +37,13 @@ TriNedelecBasis::TriNedelecBasis(void){
};
// Basis //
- basis = new vector<Polynomial>**[nRefSpace];
+ basis = new vector<Polynomial>**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
for(size_t e = 0; e < 3; e++){
vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
vector<Polynomial> tmp2 = lagrange[edgeIdx[s][e][0]].gradient();
@@ -65,11 +66,10 @@ TriNedelecBasis::TriNedelecBasis(void){
}
TriNedelecBasis::~TriNedelecBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index e6399f5a9a4db753427c91586bcffec425a70f2b..79e75c1708322730b367baa501f10bedf1caa6dc 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -1,24 +1,16 @@
-#include "TriNodeBasis.h"
-#include "TriReferenceSpace.h"
#include "Legendre.h"
+#include "GmshDefines.h"
+#include "ReferenceSpaceManager.h"
+
+#include "TriNodeBasis.h"
using namespace std;
TriNodeBasis::TriNodeBasis(size_t order){
- // Reference Space //
- refSpace = new TriReferenceSpace;
- nRefSpace = getReferenceSpace().getNReferenceSpace();
-
- const vector<vector<vector<size_t> > >&
- edgeIdx = refSpace->getEdgeNodeIndex();
-
- const vector<vector<vector<size_t> > >&
- faceIdx = refSpace->getFaceNodeIndex();
-
// Set BasisTwo Type //
this->order = order;
- type = 0;
+ type = TYPE_TRI;
dim = 2;
nVertex = 3;
@@ -27,6 +19,15 @@ TriNodeBasis::TriNodeBasis(size_t order){
nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
+ // Reference Space //
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
+
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
+
+ const vector<vector<vector<size_t> > >&
+ faceIdx = ReferenceSpaceManager::getFaceNodeIndex(type);
+
// Legendre Polynomial //
const int orderMinus = order - 1;
@@ -49,20 +50,20 @@ TriNodeBasis::TriNodeBasis(size_t order){
};
// Basis //
- basis = new Polynomial**[nRefSpace];
+ basis = new Polynomial**[nOrientation];
- for(size_t s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nOrientation; s++)
basis[s] = new Polynomial*[nFunction];
// Vertex Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
basis[s][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]);
}
// Edge Based //
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex;
for(size_t e = 0; e < 3; e++){
@@ -85,7 +86,7 @@ TriNodeBasis::TriNodeBasis(size_t order){
// have only ONE face: the face '0'
const int orderMinusTwo = order - 2;
- for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex + nEdge;
for(int l1 = 1; l1 < orderMinus; l1++){
@@ -115,11 +116,10 @@ TriNodeBasis::TriNodeBasis(size_t order){
}
TriNodeBasis::~TriNodeBasis(void){
- // ReferenceSpace //
- delete refSpace;
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
// Basis //
- for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];