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

Notion of Reference Space Manager: cleaner code -- seems ok

parent 0802171b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 302 additions and 284 deletions
FunctionSpace/Basis.h
+ 43
− 31
View file @ 551681df
...@@ -3,7 +3,6 @@ ...@@ -3,7 +3,6 @@
#include <string> #include <string>
#include "MElement.h" #include "MElement.h"
#include "ReferenceSpace.h"
/** /**
@interface Basis @interface Basis
...@@ -20,38 +19,32 @@ ...@@ -20,38 +19,32 @@
The i-th row of these matrices is always refering to the The i-th row of these matrices is always refering to the
i-th function of the basis. i-th function of the basis.
Depending on the nature of the returned value Depending on the nature of the returned value (scalar or vector),
(scalar or vector), the columns are organized the columns are organized diferently.
diferently.
For scalar values, we have: For scalar values, we have:
@li The j-th column of the i-th row is @li The j-th column of the i-th row is
the evaluation of the i-th function at the j-th point the evaluation of the i-th function at the j-th point
For vectorial values, we have: For vectorial values, we have:
@li The j-th column of the i-th row is @li The j-th column of the i-th row is the first coordinate of
the first coordinate of
the evaluation of the i-th function at the 3 x j-th point 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 @li The (j-th + 1) column of the i-th row is the second coordinate of
the second coordinate of
the evaluation of the i-th function at the 3 x j-th point 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 @li The (j-th + 2) column of the i-th row is the third coordinate of
the third coordinate of
the evaluation of the i-th function at the 3 x j-th point the evaluation of the i-th function at the 3 x j-th point
*/ */
class Basis{ class Basis{
protected: protected:
ReferenceSpace* refSpace;
size_t nRefSpace;
bool scalar; bool scalar;
bool local; bool local;
size_t order; size_t order;
size_t type; size_t type;
size_t form;
size_t dim; size_t dim;
size_t nVertex; size_t nVertex;
...@@ -72,6 +65,7 @@ class Basis{ ...@@ -72,6 +65,7 @@ class Basis{
// Type of Basis // // Type of Basis //
size_t getOrder(void) const; size_t getOrder(void) const;
size_t getType(void) const; size_t getType(void) const;
size_t getForm(void) const;
size_t getDim(void) const; size_t getDim(void) const;
// Number of Functions // // Number of Functions //
...@@ -81,9 +75,6 @@ class Basis{ ...@@ -81,9 +75,6 @@ class Basis{
size_t getNCellBased(void) const; size_t getNCellBased(void) const;
size_t getNFunction(void) const; size_t getNFunction(void) const;
// Reference Element //
const ReferenceSpace& getReferenceSpace(void) const;
// Direct Access to Evaluated Functions // // Direct Access to Evaluated Functions //
virtual void getFunctions(fullMatrix<double>& retValues, virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element, const MElement& element,
...@@ -98,8 +89,10 @@ class Basis{ ...@@ -98,8 +89,10 @@ class Basis{
double u, double v, double w) const = 0; double u, double v, double w) const = 0;
// Precompute Functions // // Precompute Functions //
virtual void preEvaluateFunctions(const fullMatrix<double>& point) const = 0; virtual
virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const = 0; void preEvaluateFunctions(const fullMatrix<double>& point) const = 0;
virtual
void preEvaluateDerivatives(const fullMatrix<double>& point) const = 0;
// Access to Precomputed Functions // // Access to Precomputed Functions //
virtual const fullMatrix<double>& virtual const fullMatrix<double>&
...@@ -140,8 +133,8 @@ class Basis{ ...@@ -140,8 +133,8 @@ class Basis{
@li true, if this is a scalar Basis @li true, if this is a scalar Basis
@li false, if this is a vectorial Basis @li false, if this is a vectorial Basis
Scalar basis are sets of Polynomial%s, and Scalar basis are sets of Polynomial%s,
Vectorial basis are sets of Vector%s of Polynomial%s and Vectorial basis are sets of Vector%s of Polynomial%s
** **
@fn Basis::isLocal @fn Basis::isLocal
...@@ -155,7 +148,19 @@ class Basis{ ...@@ -155,7 +148,19 @@ class Basis{
** **
@fn Basis::getType @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 0 for 0-form
@li 1 for 1-form @li 1 for 1-form
@li 2 for 2-form @li 2 for 2-form
...@@ -187,10 +192,6 @@ class Basis{ ...@@ -187,10 +192,6 @@ class Basis{
(or Vector%s of Polynomial%s) Functions in this 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 @fn Basis::getFunctions(fullMatrix<double>&, const MElement&, double, double, double) const
@param retValues An allocated matrix @param retValues An allocated matrix
@param element A MElement @param element A MElement
...@@ -213,6 +214,17 @@ class Basis{ ...@@ -213,6 +214,17 @@ class Basis{
and for the given orientation 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 @fn Basis::preEvaluateFunctions
@param point A Matrix with points coordinate @param point A Matrix with points coordinate
(each line is a point and got 3 coordinates, i.e. 3 rows) (each line is a point and got 3 coordinates, i.e. 3 rows)
...@@ -243,8 +255,8 @@ class Basis{ ...@@ -243,8 +255,8 @@ class Basis{
@return Returns a Matrix with the PreEvaluated basis functions derivatives @return Returns a Matrix with the PreEvaluated basis functions derivatives
(see Basis::preEvaluateDerivatives()), with the given orientation (see Basis::preEvaluateDerivatives()), with the given orientation
If no PreEvaluation of the gradient has been done before calling this function, If no PreEvaluation of the gradient has been done
an Exception is thrown before calling this function, an Exception is thrown
** **
@fn Basis::getPreEvaluatedFunctions(const MElement&) const @fn Basis::getPreEvaluatedFunctions(const MElement&) const
...@@ -279,6 +291,10 @@ inline size_t Basis::getType(void) const{ ...@@ -279,6 +291,10 @@ inline size_t Basis::getType(void) const{
return type; return type;
} }
inline size_t Basis::getForm(void) const{
return form;
}
inline size_t Basis::getDim(void) const{ inline size_t Basis::getDim(void) const{
return dim; return dim;
} }
...@@ -303,8 +319,4 @@ inline size_t Basis::getNFunction(void) const{ ...@@ -303,8 +319,4 @@ inline size_t Basis::getNFunction(void) const{
return nFunction; return nFunction;
} }
inline const ReferenceSpace& Basis::getReferenceSpace(void) const{
return *refSpace;
}
#endif #endif
FunctionSpace/BasisHierarchical0Form.cpp
+ 46
− 34
View file @ 551681df
#include <sstream> #include <sstream>
#include "Exception.h" #include "Exception.h"
#include "ReferenceSpaceManager.h"
#include "BasisHierarchical0Form.h" #include "BasisHierarchical0Form.h"
using namespace std; using namespace std;
...@@ -8,6 +10,9 @@ BasisHierarchical0Form::BasisHierarchical0Form(void){ ...@@ -8,6 +10,9 @@ BasisHierarchical0Form::BasisHierarchical0Form(void){
// Scalar Basis ? // // Scalar Basis ? //
scalar = true; scalar = true;
// 0-Form //
form = 0;
// Grad Basis // // Grad Basis //
hasGrad = false; hasGrad = false;
grad = NULL; grad = NULL;
...@@ -21,9 +26,11 @@ BasisHierarchical0Form::BasisHierarchical0Form(void){ ...@@ -21,9 +26,11 @@ BasisHierarchical0Form::BasisHierarchical0Form(void){
} }
BasisHierarchical0Form::~BasisHierarchical0Form(void){ BasisHierarchical0Form::~BasisHierarchical0Form(void){
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Grad Basis // // Grad Basis //
if(hasGrad){ 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++) for(size_t j = 0; j < nFunction; j++)
delete grad[i][j]; delete grad[i][j];
...@@ -35,14 +42,14 @@ BasisHierarchical0Form::~BasisHierarchical0Form(void){ ...@@ -35,14 +42,14 @@ BasisHierarchical0Form::~BasisHierarchical0Form(void){
// PreEvaluation // // PreEvaluation //
if(preEvaluated){ if(preEvaluated){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i]; delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction; delete[] preEvaluatedFunction;
} }
if(preEvaluatedGrad){ if(preEvaluatedGrad){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedGradFunction[i]; delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction; delete[] preEvaluatedGradFunction;
...@@ -55,7 +62,7 @@ getFunctions(fullMatrix<double>& retValues, ...@@ -55,7 +62,7 @@ getFunctions(fullMatrix<double>& retValues,
double u, double v, double w) const{ double u, double v, double w) const{
// Define Orientation // // Define Orientation //
const size_t orientation = refSpace->getReferenceSpace(element); const size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Matrix // // Fill Matrix //
for(size_t i = 0; i < nFunction; i++) for(size_t i = 0; i < nFunction; i++)
...@@ -80,7 +87,7 @@ void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues, ...@@ -80,7 +87,7 @@ void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues,
getGrad(); getGrad();
// Define Orientation // // Define Orientation //
const size_t orientation = refSpace->getReferenceSpace(element); const size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Matrix // // Fill Matrix //
for(size_t i = 0; i < nFunction; i++){ for(size_t i = 0; i < nFunction; i++){
...@@ -92,62 +99,62 @@ void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues, ...@@ -92,62 +99,62 @@ void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues,
void BasisHierarchical0Form:: void BasisHierarchical0Form::
preEvaluateFunctions(const fullMatrix<double>& point) const{ preEvaluateFunctions(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Delete if older // // Delete if older //
if(preEvaluated){ if(preEvaluated){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i]; delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction; delete[] preEvaluatedFunction;
} }
// Alloc // // Alloc //
const size_t nPoint = point.size1(); const size_t nPoint = point.size1();
preEvaluatedFunction = new fullMatrix<double>*[nRefSpace]; preEvaluatedFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
preEvaluatedFunction[i] = preEvaluatedFunction[i] = new fullMatrix<double>(nFunction, nPoint);
new fullMatrix<double>(nFunction, nPoint);
// Fill Matrix // // 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 j = 0; j < nFunction; j++)
for(size_t k = 0; k < nPoint; k++) for(size_t k = 0; k < nPoint; k++)
(*preEvaluatedFunction[i])(j, k) = (*preEvaluatedFunction[i])(j, k) = basis[i][j]->at(point(k, 0),
basis[i][j]->at(point(k, 0), point(k, 1),
point(k, 1), point(k, 2));
point(k, 2));
// PreEvaluated // // PreEvaluated //
preEvaluated = true; preEvaluated = true;
} }
void BasisHierarchical0Form:: void BasisHierarchical0Form::
preEvaluateDerivatives(const fullMatrix<double>& point) const{ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Build Grad // // Build Grad //
if(!hasGrad) if(!hasGrad)
getGrad(); getGrad();
// Delete if older // // Delete if older //
if(preEvaluatedGrad){ if(preEvaluatedGrad){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedGradFunction[i]; delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction; delete[] preEvaluatedGradFunction;
} }
// Alloc // // Alloc //
const size_t nPoint = point.size1(); const size_t nPoint = point.size1();
const size_t nPoint3 = nPoint * 3; const size_t nPoint3 = nPoint * 3;
preEvaluatedGradFunction = new fullMatrix<double>*[nRefSpace]; preEvaluatedGradFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
preEvaluatedGradFunction[i] = preEvaluatedGradFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix // // Fill Matrix //
fullVector<double> tmp(3); 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 j = 0; j < nFunction; j++){
for(size_t k = 0; k < nPoint; k++){ for(size_t k = 0; k < nPoint; k++){
tmp = Polynomial::at(*grad[i][j], tmp = Polynomial::at(*grad[i][j],
...@@ -168,18 +175,21 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{ ...@@ -168,18 +175,21 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const fullMatrix<double>& BasisHierarchical0Form:: const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedFunctions(const MElement& element) const{ getPreEvaluatedFunctions(const MElement& element) const{
return getPreEvaluatedFunctions(refSpace->getReferenceSpace(element)); return
getPreEvaluatedFunctions(ReferenceSpaceManager::getOrientation(element));
} }
const fullMatrix<double>& BasisHierarchical0Form:: const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedDerivatives(const MElement& element) const{ getPreEvaluatedDerivatives(const MElement& element) const{
return getPreEvaluatedDerivatives(refSpace->getReferenceSpace(element)); return
getPreEvaluatedDerivatives(ReferenceSpaceManager::getOrientation(element));
} }
const fullMatrix<double>& BasisHierarchical0Form:: const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedFunctions(size_t orientation) const{ getPreEvaluatedFunctions(size_t orientation) const{
if(!preEvaluated) if(!preEvaluated)
throw Exception("getPreEvaluatedFunction: function has not been preEvaluated"); throw Exception
("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[orientation]; return *preEvaluatedFunction[orientation];
} }
...@@ -187,23 +197,25 @@ getPreEvaluatedFunctions(size_t orientation) const{ ...@@ -187,23 +197,25 @@ getPreEvaluatedFunctions(size_t orientation) const{
const fullMatrix<double>& BasisHierarchical0Form:: const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedDerivatives(size_t orientation) const{ getPreEvaluatedDerivatives(size_t orientation) const{
if(!preEvaluatedGrad) if(!preEvaluatedGrad)
throw Exception("getPreEvaluatedDerivative: gradient has not been preEvaluated"); throw Exception
("getPreEvaluatedDerivative: gradient has not been preEvaluated");
return *preEvaluatedGradFunction[orientation]; return *preEvaluatedGradFunction[orientation];
} }
void BasisHierarchical0Form::getGrad(void) const{ void BasisHierarchical0Form::getGrad(void) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Alloc // // 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[s] = new vector<Polynomial>*[nFunction];
// Grad // // 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++) for(size_t f = 0 ; f < nFunction; f++)
grad[s][f] = grad[s][f] = new vector<Polynomial>(basis[s][f]->gradient());
new vector<Polynomial>(basis[s][f]->gradient());
// Has Grad // // Has Grad //
hasGrad = true; hasGrad = true;
......
FunctionSpace/BasisHierarchical1Form.cpp
+ 43
− 31
View file @ 551681df
#include <sstream> #include <sstream>
#include "Exception.h" #include "Exception.h"
#include "ReferenceSpaceManager.h"
#include "BasisHierarchical1Form.h" #include "BasisHierarchical1Form.h"
using namespace std; using namespace std;
...@@ -8,6 +10,9 @@ BasisHierarchical1Form::BasisHierarchical1Form(void){ ...@@ -8,6 +10,9 @@ BasisHierarchical1Form::BasisHierarchical1Form(void){
// Scalar Basis ?// // Scalar Basis ?//
scalar = false; scalar = false;
// 1-Form //
form = 1;
// Curl Basis // // Curl Basis //
hasCurl = false; hasCurl = false;
curl = NULL; curl = NULL;
...@@ -21,9 +26,11 @@ BasisHierarchical1Form::BasisHierarchical1Form(void){ ...@@ -21,9 +26,11 @@ BasisHierarchical1Form::BasisHierarchical1Form(void){
} }
BasisHierarchical1Form::~BasisHierarchical1Form(void){ BasisHierarchical1Form::~BasisHierarchical1Form(void){
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Curl Basis // // Curl Basis //
if(hasCurl){ 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++) for(size_t j = 0; j < nFunction; j++)
delete curl[i][j]; delete curl[i][j];
...@@ -35,14 +42,14 @@ BasisHierarchical1Form::~BasisHierarchical1Form(void){ ...@@ -35,14 +42,14 @@ BasisHierarchical1Form::~BasisHierarchical1Form(void){
// PreEvaluation // // PreEvaluation //
if(preEvaluated){ if(preEvaluated){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i]; delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction; delete[] preEvaluatedFunction;
} }
if(preEvaluatedCurl){ if(preEvaluatedCurl){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedCurlFunction[i]; delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction; delete[] preEvaluatedCurlFunction;
...@@ -55,12 +62,11 @@ getFunctions(fullMatrix<double>& retValues, ...@@ -55,12 +62,11 @@ getFunctions(fullMatrix<double>& retValues,
double u, double v, double w) const{ double u, double v, double w) const{
// Define Orientation // // Define Orientation //
size_t orientation = refSpace->getReferenceSpace(element); size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Vector // // Fill Vector //
for(size_t i = 0; i < nFunction; i++){ for(size_t i = 0; i < nFunction; i++){
fullVector<double> eval = fullVector<double> eval = Polynomial::at(*basis[orientation][i], u, v, w);
Polynomial::at(*basis[orientation][i], u, v, w);
retValues(i, 0) = eval(0); retValues(i, 0) = eval(0);
retValues(i, 1) = eval(1); retValues(i, 1) = eval(1);
...@@ -75,8 +81,7 @@ getFunctions(fullMatrix<double>& retValues, ...@@ -75,8 +81,7 @@ getFunctions(fullMatrix<double>& retValues,
// Fill Vector // // Fill Vector //
for(size_t i = 0; i < nFunction; i++){ for(size_t i = 0; i < nFunction; i++){
fullVector<double> eval = fullVector<double> eval = Polynomial::at(*basis[orientation][i], u, v, w);
Polynomial::at(*basis[orientation][i], u, v, w);
retValues(i, 0) = eval(0); retValues(i, 0) = eval(0);
retValues(i, 1) = eval(1); retValues(i, 1) = eval(1);
...@@ -92,9 +97,11 @@ void BasisHierarchical1Form::getDerivative(fullMatrix<double>& retValues, ...@@ -92,9 +97,11 @@ void BasisHierarchical1Form::getDerivative(fullMatrix<double>& retValues,
void BasisHierarchical1Form:: void BasisHierarchical1Form::
preEvaluateFunctions(const fullMatrix<double>& point) const{ preEvaluateFunctions(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Delete if older // // Delete if older //
if(preEvaluated){ if(preEvaluated){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i]; delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction; delete[] preEvaluatedFunction;
...@@ -103,16 +110,15 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{ ...@@ -103,16 +110,15 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{
// Alloc // // Alloc //
const size_t nPoint = point.size1(); const size_t nPoint = point.size1();
const size_t nPoint3 = nPoint * 3; const size_t nPoint3 = nPoint * 3;
preEvaluatedFunction = new fullMatrix<double>*[nRefSpace]; preEvaluatedFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
preEvaluatedFunction[i] = preEvaluatedFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix // // Fill Matrix //
fullVector<double> tmp(3); 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 j = 0; j < nFunction; j++){
for(size_t k = 0; k < nPoint; k++){ for(size_t k = 0; k < nPoint; k++){
tmp = Polynomial::at(*basis[i][j], tmp = Polynomial::at(*basis[i][j],
...@@ -133,31 +139,32 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{ ...@@ -133,31 +139,32 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{
void BasisHierarchical1Form:: void BasisHierarchical1Form::
preEvaluateDerivatives(const fullMatrix<double>& point) const{ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Build Curl // // Build Curl //
if(!hasCurl) if(!hasCurl)
getCurl(); getCurl();
// Delete if older // // Delete if older //
if(preEvaluatedCurl){ if(preEvaluatedCurl){
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedCurlFunction[i]; delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction; delete[] preEvaluatedCurlFunction;
} }
// Alloc // // Alloc //
const size_t nPoint = point.size1(); const size_t nPoint = point.size1();
const size_t nPoint3 = nPoint * 3; const size_t nPoint3 = nPoint * 3;
preEvaluatedCurlFunction = new fullMatrix<double>*[nRefSpace]; preEvaluatedCurlFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nRefSpace; i++) for(size_t i = 0; i < nOrientation; i++)
preEvaluatedCurlFunction[i] = preEvaluatedCurlFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix // // Fill Matrix //
fullVector<double> tmp(3); 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 j = 0; j < nFunction; j++){
for(size_t k = 0; k < nPoint; k++){ for(size_t k = 0; k < nPoint; k++){
tmp = Polynomial::at(*curl[i][j], tmp = Polynomial::at(*curl[i][j],
...@@ -178,18 +185,21 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{ ...@@ -178,18 +185,21 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
const fullMatrix<double>& BasisHierarchical1Form:: const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedFunctions(const MElement& element) const{ getPreEvaluatedFunctions(const MElement& element) const{
return getPreEvaluatedFunctions(refSpace->getReferenceSpace(element)); return
getPreEvaluatedFunctions(ReferenceSpaceManager::getOrientation(element));
} }
const fullMatrix<double>& BasisHierarchical1Form:: const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedDerivatives(const MElement& element) const{ getPreEvaluatedDerivatives(const MElement& element) const{
return getPreEvaluatedDerivatives(refSpace->getReferenceSpace(element)); return
getPreEvaluatedDerivatives(ReferenceSpaceManager::getOrientation(element));
} }
const fullMatrix<double>& BasisHierarchical1Form:: const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedFunctions(size_t orientation) const{ getPreEvaluatedFunctions(size_t orientation) const{
if(!preEvaluated) if(!preEvaluated)
throw Exception("getPreEvaluatedFunction: function has not been preEvaluated"); throw
Exception("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[orientation]; return *preEvaluatedFunction[orientation];
} }
...@@ -197,23 +207,25 @@ getPreEvaluatedFunctions(size_t orientation) const{ ...@@ -197,23 +207,25 @@ getPreEvaluatedFunctions(size_t orientation) const{
const fullMatrix<double>& BasisHierarchical1Form:: const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedDerivatives(size_t orientation) const{ getPreEvaluatedDerivatives(size_t orientation) const{
if(!preEvaluatedCurl) if(!preEvaluatedCurl)
throw Exception("getPreEvaluatedDerivative: curl has not been preEvaluated"); throw
Exception("getPreEvaluatedDerivative: curl has not been preEvaluated");
return *preEvaluatedCurlFunction[orientation]; return *preEvaluatedCurlFunction[orientation];
} }
void BasisHierarchical1Form::getCurl(void) const{ void BasisHierarchical1Form::getCurl(void) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Alloc // // 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[s] = new vector<Polynomial>*[nFunction];
// Curl // // 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++) for(size_t f = 0 ; f < nFunction; f++)
curl[s][f] = curl[s][f] = new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
// Has Curl // // Has Curl //
hasCurl = true; hasCurl = true;
......
FunctionSpace/BasisLagrange.cpp
+ 3
− 1
View file @ 551681df
#include "ReferenceSpaceManager.h"
#include "BasisLagrange.h" #include "BasisLagrange.h"
using namespace std; using namespace std;
BasisLagrange::BasisLagrange(void){ BasisLagrange::BasisLagrange(void){
scalar = true; scalar = true;
form = 0;
preEvaluated = false; preEvaluated = false;
preEvaluatedGrad = false; preEvaluatedGrad = false;
...@@ -38,7 +40,7 @@ getFunctions(fullMatrix<double>& retValues, ...@@ -38,7 +40,7 @@ getFunctions(fullMatrix<double>& retValues,
retValues = tmp.transpose(); retValues = tmp.transpose();
// Permute retValues, accordingly to ReferenceSpace // Permute retValues, accordingly to ReferenceSpace
permutation(refSpace->getReferenceSpace(element), retValues); permutation(ReferenceSpaceManager::getOrientation(element), retValues);
} }
void BasisLagrange:: void BasisLagrange::
......
FunctionSpace/CMakeLists.txt
+ 2
− 0
View file @ 551681df
...@@ -18,6 +18,8 @@ set(SRC ...@@ -18,6 +18,8 @@ set(SRC
PyrReferenceSpace.cpp PyrReferenceSpace.cpp
PriReferenceSpace.cpp PriReferenceSpace.cpp
ReferenceSpaceManager.cpp
Basis.cpp Basis.cpp
BasisLocal.cpp BasisLocal.cpp
BasisGenerator.cpp BasisGenerator.cpp
......
FunctionSpace/FunctionSpace.cpp
+ 5
− 8
View file @ 551681df
#include <sstream> #include <sstream>
#include "FunctionSpace.h"
#include "ReferenceSpaceManager.h"
#include "BasisGenerator.h" #include "BasisGenerator.h"
#include "ElementType.h" #include "ElementType.h"
#include "Exception.h" #include "Exception.h"
#include "FunctionSpace.h"
using namespace std; using namespace std;
FunctionSpace::FunctionSpace(void){ FunctionSpace::FunctionSpace(void){
...@@ -35,9 +38,6 @@ FunctionSpace::~FunctionSpace(void){ ...@@ -35,9 +38,6 @@ FunctionSpace::~FunctionSpace(void){
// Element To GoD // // Element To GoD //
if(eToGod) if(eToGod)
delete eToGod; delete eToGod;
// Unorient GroupOfElement //
goe->unoriented();
} }
void FunctionSpace::build(GroupOfElement& goe, void FunctionSpace::build(GroupOfElement& goe,
...@@ -47,9 +47,6 @@ void FunctionSpace::build(GroupOfElement& goe, ...@@ -47,9 +47,6 @@ void FunctionSpace::build(GroupOfElement& goe,
this->goe = &goe; this->goe = &goe;
this->mesh = &(goe.getMesh()); this->mesh = &(goe.getMesh());
// Orient All Elements //
goe.orientAllElements(basis); // NOT SEXY: TO BE REMOVED
// Get Geo Data (WARNING HOMOGENE MESH REQUIRED)// // Get Geo Data (WARNING HOMOGENE MESH REQUIRED)//
const MElement& element = goe.get(0); const MElement& element = goe.get(0);
MElement& myElement = MElement& myElement =
...@@ -227,7 +224,7 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{ ...@@ -227,7 +224,7 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
// Create New Element With Permuted Vertices // // Create New Element With Permuted Vertices //
// Permutation // Permutation
const vector<size_t>& vPerm = const vector<size_t>& vPerm =
(*this->basis)[0]->getReferenceSpace().getNodeIndexFromABCtoUVW(elem); ReferenceSpaceManager::getNodeIndexFromABCtoUVW(elem);
// Permuted Vertices // Permuted Vertices
const size_t nVertex = element.getNumPrimaryVertices(); const size_t nVertex = element.getNumPrimaryVertices();
......
FunctionSpace/FunctionSpaceScalar.cpp
+ 1
− 3
View file @ 551681df
...@@ -38,9 +38,7 @@ interpolateDerivativeInABC(const MElement& element, ...@@ -38,9 +38,7 @@ interpolateDerivativeInABC(const MElement& element,
double abc[3]) const{ double abc[3]) const{
// Get Jacobian // // Get Jacobian //
fullMatrix<double> invJac(3, 3); fullMatrix<double> invJac(3, 3);
(*basis)[0]->getReferenceSpace().getJacobian(element, ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
abc[0], abc[1], abc[2],
invJac);
invJac.invertInPlace(); invJac.invertInPlace();
// Get Basis Functions // // Get Basis Functions //
......
FunctionSpace/FunctionSpaceScalar.h
+ 7
− 9
View file @ 551681df
#ifndef _FUNCTIONSPACESCALAR_H_ #ifndef _FUNCTIONSPACESCALAR_H_
#define _FUNCTIONSPACESCALAR_H_ #define _FUNCTIONSPACESCALAR_H_
#include "ReferenceSpaceManager.h"
#include "FunctionSpace.h" #include "FunctionSpace.h"
/** /**
...@@ -110,9 +111,8 @@ interpolate(const MElement& element, ...@@ -110,9 +111,8 @@ interpolate(const MElement& element,
// Get ABC Space coordinate // // Get ABC Space coordinate //
double abc[3]; double abc[3];
(*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element, ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
xyz(0), xyz(1), xyz(2),
abc);
// Interpolate in ABC // // Interpolate in ABC //
return interpolateInABC(element, coef, abc); return interpolateInABC(element, coef, abc);
} }
...@@ -124,9 +124,8 @@ interpolateInRefSpace(const MElement& element, ...@@ -124,9 +124,8 @@ interpolateInRefSpace(const MElement& element,
// Get ABC Space coordinate // // Get ABC Space coordinate //
double abc[3]; double abc[3];
(*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element, ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
uvw(0), uvw(1), uvw(2),
abc);
// Interpolate in ABC // // Interpolate in ABC //
return interpolateInABC(element, coef, abc); return interpolateInABC(element, coef, abc);
} }
...@@ -138,9 +137,8 @@ interpolateDerivative(const MElement& element, ...@@ -138,9 +137,8 @@ interpolateDerivative(const MElement& element,
// Get ABC Space coordinate // // Get ABC Space coordinate //
double abc[3]; double abc[3];
(*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element, ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
xyz(0), xyz(1), xyz(2),
abc);
// Interpolate in ABC // // Interpolate in ABC //
return interpolateDerivativeInABC(element, coef, abc); return interpolateDerivativeInABC(element, coef, abc);
} }
......
FunctionSpace/FunctionSpaceVector.cpp
+ 1
− 3
View file @ 551681df
...@@ -18,9 +18,7 @@ interpolateInABC(const MElement& element, ...@@ -18,9 +18,7 @@ interpolateInABC(const MElement& element,
// Get Jacobian // // Get Jacobian //
fullMatrix<double> invJac(3, 3); fullMatrix<double> invJac(3, 3);
(*basis)[0]->getReferenceSpace().getJacobian(element, ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
abc[0], abc[1], abc[2],
invJac);
invJac.invertInPlace(); invJac.invertInPlace();
// Get Basis Functions // // Get Basis Functions //
......
FunctionSpace/FunctionSpaceVector.h
+ 5
− 6
View file @ 551681df
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#define _FUNCTIONSPACEVECTOR_H_ #define _FUNCTIONSPACEVECTOR_H_
#include "fullMatrix.h" #include "fullMatrix.h"
#include "FunctionSpaceScalar.h"
#include "FunctionSpace.h" #include "FunctionSpace.h"
/** /**
...@@ -103,9 +104,8 @@ interpolate(const MElement& element, ...@@ -103,9 +104,8 @@ interpolate(const MElement& element,
// Get ABC Space coordinate // // Get ABC Space coordinate //
double abc[3]; double abc[3];
(*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element, ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
xyz(0), xyz(1), xyz(2),
abc);
// Interpolate in ABC // // Interpolate in ABC //
return interpolateInABC(element, coef, abc); return interpolateInABC(element, coef, abc);
} }
...@@ -117,9 +117,8 @@ interpolateInRefSpace(const MElement& element, ...@@ -117,9 +117,8 @@ interpolateInRefSpace(const MElement& element,
// Get ABC Space coordinate // // Get ABC Space coordinate //
double abc[3]; double abc[3];
(*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element, ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
uvw(0), uvw(1), uvw(2),
abc);
// Interpolate in ABC // // Interpolate in ABC //
return interpolateInABC(element, coef, abc); return interpolateInABC(element, coef, abc);
} }
......
FunctionSpace/HexLagrangeBasis.cpp
+ 5
− 9
View file @ 551681df
#include "HexLagrangeBasis.h"
#include "HexReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h" #include "ElementType.h"
#include "GmshDefines.h"
#include "pointsGenerators.h"
#include "HexLagrangeBasis.h"
HexLagrangeBasis::HexLagrangeBasis(size_t order){ HexLagrangeBasis::HexLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1 // If order 0 (Nedelec): use order 1
...@@ -11,7 +12,7 @@ HexLagrangeBasis::HexLagrangeBasis(size_t order){ ...@@ -11,7 +12,7 @@ HexLagrangeBasis::HexLagrangeBasis(size_t order){
// Set Basis Type // // Set Basis Type //
this->order = order; this->order = order;
type = 0; type = TYPE_HEX;
dim = 3; dim = 3;
nVertex = 8; nVertex = 8;
...@@ -25,14 +26,9 @@ HexLagrangeBasis::HexLagrangeBasis(size_t order){ ...@@ -25,14 +26,9 @@ HexLagrangeBasis::HexLagrangeBasis(size_t order){
// Init Lagrange Point // // Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsHexahedron(order, false)); lPoint = new fullMatrix<double>(gmshGeneratePointsHexahedron(order, false));
// Reference Space //
refSpace = new HexReferenceSpace;
nRefSpace = getReferenceSpace().getNReferenceSpace();
} }
HexLagrangeBasis::~HexLagrangeBasis(void){ HexLagrangeBasis::~HexLagrangeBasis(void){
delete lBasis; delete lBasis;
delete lPoint; delete lPoint;
delete refSpace;
} }
FunctionSpace/HexNodeBasis.cpp
+ 23
− 23
View file @ 551681df
#include "HexNodeBasis.h"
#include "HexReferenceSpace.h"
#include "Legendre.h" #include "Legendre.h"
#include "GmshDefines.h"
#include "ReferenceSpaceManager.h"
#include "HexNodeBasis.h"
using namespace std; using namespace std;
HexNodeBasis::HexNodeBasis(size_t order){ 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 // // Set Basis Type //
this->order = order; this->order = order;
type = 0; type = TYPE_HEX;
dim = 3; dim = 3;
nVertex = 8; nVertex = 8;
...@@ -27,6 +19,15 @@ HexNodeBasis::HexNodeBasis(size_t order){ ...@@ -27,6 +19,15 @@ HexNodeBasis::HexNodeBasis(size_t order){
nCell = (order - 1) * (order - 1) * (order - 1); nCell = (order - 1) * (order - 1) * (order - 1);
nFunction = nVertex + nEdge + nFace + nCell; 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 // // Legendre Polynomial //
Polynomial* legendre = new Polynomial[order]; Polynomial* legendre = new Polynomial[order];
Legendre::integrated(legendre, order); Legendre::integrated(legendre, order);
...@@ -103,13 +104,13 @@ HexNodeBasis::HexNodeBasis(size_t order){ ...@@ -103,13 +104,13 @@ HexNodeBasis::HexNodeBasis(size_t order){
}; };
// Basis // // 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]; basis[s] = new Polynomial*[nFunction];
// Vertex Based // // 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][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]); basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]); basis[s][2] = new Polynomial(lagrange[2]);
...@@ -121,7 +122,7 @@ HexNodeBasis::HexNodeBasis(size_t order){ ...@@ -121,7 +122,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
} }
// Edge Based // // Edge Based //
for(size_t s = 0; s < nRefSpace; s++){ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex; size_t i = nVertex;
for(size_t e = 0; e < 12; e++){ for(size_t e = 0; e < 12; e++){
...@@ -139,7 +140,7 @@ HexNodeBasis::HexNodeBasis(size_t order){ ...@@ -139,7 +140,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
} }
// Face Based // // Face Based //
for(size_t s = 0; s < nRefSpace; s++){ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex + nEdge; size_t i = nVertex + nEdge;
for(size_t f = 0; f < 6; f++){ for(size_t f = 0; f < 6; f++){
...@@ -173,7 +174,7 @@ HexNodeBasis::HexNodeBasis(size_t order){ ...@@ -173,7 +174,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
py = py - Polynomial(1, 0, 0, 0); py = py - Polynomial(1, 0, 0, 0);
pz = pz - 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; size_t i = nVertex + nEdge + nFace;
for(size_t l1 = 1; l1 < order; l1++){ for(size_t l1 = 1; l1 < order; l1++){
...@@ -206,7 +207,7 @@ HexNodeBasis::HexNodeBasis(size_t order){ ...@@ -206,7 +207,7 @@ HexNodeBasis::HexNodeBasis(size_t order){
Polynomial mapZ(Polynomial(0.5, 0, 0, 1) + Polynomial mapZ(Polynomial(0.5, 0, 0, 1) +
Polynomial(0.5, 0, 0, 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++){ for(size_t i = 0; i < nFunction; i++){
Polynomial* tmp; Polynomial* tmp;
tmp = basis[s][i]; tmp = basis[s][i];
...@@ -220,11 +221,10 @@ HexNodeBasis::HexNodeBasis(size_t order){ ...@@ -220,11 +221,10 @@ HexNodeBasis::HexNodeBasis(size_t order){
} }
HexNodeBasis::~HexNodeBasis(void){ HexNodeBasis::~HexNodeBasis(void){
// ReferenceSpace // const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
delete refSpace;
// Basis // // 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++) for(size_t j = 0; j < nFunction; j++)
delete basis[i][j]; delete basis[i][j];
......
FunctionSpace/LineEdgeBasis.cpp
+ 16
− 16
View file @ 551681df
#include "LineEdgeBasis.h"
#include "LineReferenceSpace.h"
#include "Legendre.h" #include "Legendre.h"
#include "GmshDefines.h"
#include "ReferenceSpaceManager.h"
#include "LineEdgeBasis.h"
using namespace std; using namespace std;
LineEdgeBasis::LineEdgeBasis(size_t order){ 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 // // Set Basis Type //
this->order = order; this->order = order;
type = 1; type = TYPE_LIN;
dim = 1; dim = 1;
nVertex = 0; nVertex = 0;
...@@ -24,6 +19,12 @@ LineEdgeBasis::LineEdgeBasis(size_t order){ ...@@ -24,6 +19,12 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
nCell = 0; nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell; 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 // // Legendre Polynomial //
const size_t orderPlus = order + 1; const size_t orderPlus = order + 1;
Polynomial* intLegendre = new Polynomial[orderPlus]; Polynomial* intLegendre = new Polynomial[orderPlus];
...@@ -41,13 +42,13 @@ LineEdgeBasis::LineEdgeBasis(size_t order){ ...@@ -41,13 +42,13 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
}; };
// Basis // // 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]; basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based // // Edge Based //
for(size_t s = 0; s < nRefSpace; s++){ for(size_t s = 0; s < nOrientation; s++){
size_t i = 0; size_t i = 0;
for(size_t l = 0; l < orderPlus; l++){ for(size_t l = 0; l < orderPlus; l++){
...@@ -87,11 +88,10 @@ LineEdgeBasis::LineEdgeBasis(size_t order){ ...@@ -87,11 +88,10 @@ LineEdgeBasis::LineEdgeBasis(size_t order){
} }
LineEdgeBasis::~LineEdgeBasis(void){ LineEdgeBasis::~LineEdgeBasis(void){
// ReferenceSpace // const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
delete refSpace;
// Basis // // 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++) for(size_t j = 0; j < nFunction; j++)
delete basis[i][j]; delete basis[i][j];
......
FunctionSpace/LineLagrangeBasis.cpp
+ 5
− 9
View file @ 551681df
#include "LineLagrangeBasis.h"
#include "LineReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h" #include "ElementType.h"
#include "GmshDefines.h"
#include "pointsGenerators.h"
#include "LineLagrangeBasis.h"
LineLagrangeBasis::LineLagrangeBasis(size_t order){ LineLagrangeBasis::LineLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1 // If order 0 (Nedelec): use order 1
...@@ -11,7 +12,7 @@ LineLagrangeBasis::LineLagrangeBasis(size_t order){ ...@@ -11,7 +12,7 @@ LineLagrangeBasis::LineLagrangeBasis(size_t order){
// Set Basis Type // // Set Basis Type //
this->order = order; this->order = order;
type = 0; type = TYPE_LIN;
dim = 1; dim = 1;
nVertex = 2; nVertex = 2;
...@@ -25,14 +26,9 @@ LineLagrangeBasis::LineLagrangeBasis(size_t order){ ...@@ -25,14 +26,9 @@ LineLagrangeBasis::LineLagrangeBasis(size_t order){
// Init Lagrange Point // // Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsLine(order)); lPoint = new fullMatrix<double>(gmshGeneratePointsLine(order));
// Reference Space //
refSpace = new LineReferenceSpace;
nRefSpace = getReferenceSpace().getNReferenceSpace();
} }
LineLagrangeBasis::~LineLagrangeBasis(void){ LineLagrangeBasis::~LineLagrangeBasis(void){
delete lBasis; delete lBasis;
delete lPoint; delete lPoint;
delete refSpace;
} }
FunctionSpace/LineNedelecBasis.cpp
+ 15
− 15
View file @ 551681df
#include "GmshDefines.h"
#include "ReferenceSpaceManager.h"
#include "LineNedelecBasis.h" #include "LineNedelecBasis.h"
#include "LineReferenceSpace.h"
using namespace std; using namespace std;
LineNedelecBasis::LineNedelecBasis(void){ LineNedelecBasis::LineNedelecBasis(void){
// Reference Space //
refSpace = new LineReferenceSpace;
nRefSpace = getReferenceSpace().getNReferenceSpace();
const vector<vector<vector<size_t> > >&
edgeIdx = refSpace->getEdgeNodeIndex();
// Set Basis Type // // Set Basis Type //
order = 0; order = 0;
type = 1; type = TYPE_LIN;
dim = 1; dim = 1;
nVertex = 0; nVertex = 0;
...@@ -23,6 +18,12 @@ LineNedelecBasis::LineNedelecBasis(void){ ...@@ -23,6 +18,12 @@ LineNedelecBasis::LineNedelecBasis(void){
nCell = 0; nCell = 0;
nFunction = 1; nFunction = 1;
// Reference Space //
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
const vector<vector<vector<size_t> > >&
edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
// Lagrange Polynomial // // Lagrange Polynomial //
const Polynomial lagrange[2] = const Polynomial lagrange[2] =
{ {
...@@ -34,13 +35,13 @@ LineNedelecBasis::LineNedelecBasis(void){ ...@@ -34,13 +35,13 @@ LineNedelecBasis::LineNedelecBasis(void){
}; };
// Basis // // 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]; basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) // // 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> tmp1 = lagrange[edgeIdx[s][0][1]].gradient();
vector<Polynomial> tmp2 = lagrange[edgeIdx[s][0][0]].gradient(); vector<Polynomial> tmp2 = lagrange[edgeIdx[s][0][0]].gradient();
...@@ -61,11 +62,10 @@ LineNedelecBasis::LineNedelecBasis(void){ ...@@ -61,11 +62,10 @@ LineNedelecBasis::LineNedelecBasis(void){
} }
LineNedelecBasis::~LineNedelecBasis(void){ LineNedelecBasis::~LineNedelecBasis(void){
// ReferenceSpace // const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
delete refSpace;
// Basis // // 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++) for(size_t j = 0; j < nFunction; j++)
delete basis[i][j]; delete basis[i][j];
......
FunctionSpace/LineNodeBasis.cpp
+ 17
− 17
View file @ 551681df
#include "LineNodeBasis.h"
#include "LineReferenceSpace.h"
#include "Legendre.h" #include "Legendre.h"
#include "GmshDefines.h"
#include "ReferenceSpaceManager.h"
#include "LineNodeBasis.h"
using namespace std; using namespace std;
LineNodeBasis::LineNodeBasis(size_t order){ 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 // // Set Basis Type //
this->order = order; this->order = order;
type = 0; type = TYPE_LIN;
dim = 1; dim = 1;
nVertex = 2; nVertex = 2;
...@@ -24,6 +19,12 @@ LineNodeBasis::LineNodeBasis(size_t order){ ...@@ -24,6 +19,12 @@ LineNodeBasis::LineNodeBasis(size_t order){
nCell = 0; nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell; 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 // // Legendre Polynomial //
Polynomial* intLegendre = new Polynomial[order]; Polynomial* intLegendre = new Polynomial[order];
Legendre::integrated(intLegendre, order); Legendre::integrated(intLegendre, order);
...@@ -39,19 +40,19 @@ LineNodeBasis::LineNodeBasis(size_t order){ ...@@ -39,19 +40,19 @@ LineNodeBasis::LineNodeBasis(size_t order){
}; };
// Basis // // 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]; basis[s] = new Polynomial*[nFunction];
// Vertex Based // // 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][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]); basis[s][1] = new Polynomial(lagrange[1]);
} }
// Edge Based // // Edge Based //
for(size_t s = 0; s < nRefSpace; s++){ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex; size_t i = nVertex;
for(size_t l = 1; l < order; l++){ for(size_t l = 1; l < order; l++){
...@@ -68,11 +69,10 @@ LineNodeBasis::LineNodeBasis(size_t order){ ...@@ -68,11 +69,10 @@ LineNodeBasis::LineNodeBasis(size_t order){
} }
LineNodeBasis::~LineNodeBasis(void){ LineNodeBasis::~LineNodeBasis(void){
// ReferenceSpace // const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
delete refSpace;
// Basis // // 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++) for(size_t j = 0; j < nFunction; j++)
delete basis[i][j]; delete basis[i][j];
......
FunctionSpace/QuadEdgeBasis.cpp
+ 22
− 21
View file @ 551681df
#include "QuadEdgeBasis.h"
#include "QuadReferenceSpace.h"
#include "Legendre.h" #include "Legendre.h"
#include "GmshDefines.h"
#include "ReferenceSpaceManager.h"
#include "QuadEdgeBasis.h"
using namespace std; using namespace std;
QuadEdgeBasis::QuadEdgeBasis(size_t order){ 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 // // Set Basis Type //
this->order = order; this->order = order;
type = 1; type = TYPE_QUA;
dim = 2; dim = 2;
nVertex = 0; nVertex = 0;
...@@ -27,6 +19,16 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){ ...@@ -27,6 +19,16 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
nCell = 0; nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell; 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 // // Legendre Polynomial //
const size_t orderPlus = order + 1; const size_t orderPlus = order + 1;
...@@ -68,13 +70,13 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){ ...@@ -68,13 +70,13 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
}; };
// Basis // // 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]; basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based // // Edge Based //
for(size_t s = 0; s < nRefSpace; s++){ for(size_t s = 0; s < nOrientation; s++){
size_t i = 0; size_t i = 0;
for(size_t e = 0; e < 4; e++){ for(size_t e = 0; e < 4; e++){
...@@ -115,7 +117,7 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){ ...@@ -115,7 +117,7 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
// where f = 0, because triangles // where f = 0, because triangles
// have only ONE face: the face '0' // 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; size_t i = nEdge;
// Type 1 // Type 1
...@@ -220,7 +222,7 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){ ...@@ -220,7 +222,7 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) + Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 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++){ for(size_t i = 0; i < nFunction; i++){
vector<Polynomial>* old; vector<Polynomial>* old;
vector<Polynomial> nxt(3); vector<Polynomial> nxt(3);
...@@ -241,11 +243,10 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){ ...@@ -241,11 +243,10 @@ QuadEdgeBasis::QuadEdgeBasis(size_t order){
} }
QuadEdgeBasis::~QuadEdgeBasis(void){ QuadEdgeBasis::~QuadEdgeBasis(void){
// ReferenceSpace // const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
delete refSpace;
// Basis // // 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++) for(size_t j = 0; j < nFunction; j++)
delete basis[i][j]; delete basis[i][j];
......
FunctionSpace/QuadLagrangeBasis.cpp
+ 4
− 9
View file @ 551681df
#include "QuadLagrangeBasis.h"
#include "QuadReferenceSpace.h"
#include "pointsGenerators.h"
#include "ElementType.h" #include "ElementType.h"
#include "GmshDefines.h"
#include "pointsGenerators.h"
#include "QuadLagrangeBasis.h"
QuadLagrangeBasis::QuadLagrangeBasis(size_t order){ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){
// If order 0 (Nedelec): use order 1 // If order 0 (Nedelec): use order 1
...@@ -11,7 +11,7 @@ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){ ...@@ -11,7 +11,7 @@ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){
// Set Basis Type // // Set Basis Type //
this->order = order; this->order = order;
type = 0; type = TYPE_QUA;
dim = 2; dim = 2;
nVertex = 4; nVertex = 4;
...@@ -25,14 +25,9 @@ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){ ...@@ -25,14 +25,9 @@ QuadLagrangeBasis::QuadLagrangeBasis(size_t order){
// Init Lagrange Point // // Init Lagrange Point //
lPoint = new fullMatrix<double>(gmshGeneratePointsQuadrangle(order, false)); lPoint = new fullMatrix<double>(gmshGeneratePointsQuadrangle(order, false));
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = getReferenceSpace().getNReferenceSpace();
} }
QuadLagrangeBasis::~QuadLagrangeBasis(void){ QuadLagrangeBasis::~QuadLagrangeBasis(void){
delete lBasis; delete lBasis;
delete lPoint; delete lPoint;
delete refSpace;
} }
FunctionSpace/QuadNedelecBasis.cpp
+ 17
− 17
View file @ 551681df
#include "QuadNedelecBasis.h"
#include "QuadReferenceSpace.h"
#include "Legendre.h" #include "Legendre.h"
#include "GmshDefines.h"
#include "ReferenceSpaceManager.h"
#include "QuadNedelecBasis.h"
using namespace std; using namespace std;
QuadNedelecBasis::QuadNedelecBasis(void){ QuadNedelecBasis::QuadNedelecBasis(void){
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = getReferenceSpace().getNReferenceSpace();
const vector<vector<vector<size_t> > >&
edgeIdx = refSpace->getEdgeNodeIndex();
// Set Basis Type // // Set Basis Type //
order = 0; order = 0;
type = 1; type = TYPE_QUA;
dim = 2; dim = 2;
nVertex = 0; nVertex = 0;
...@@ -24,6 +19,12 @@ QuadNedelecBasis::QuadNedelecBasis(void){ ...@@ -24,6 +19,12 @@ QuadNedelecBasis::QuadNedelecBasis(void){
nCell = 0; nCell = 0;
nFunction = 4; nFunction = 4;
// Reference Space //
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
const vector<vector<vector<size_t> > >&
edgeIdx = ReferenceSpaceManager::getEdgeNodeIndex(type);
// Lagrange & Lifting // // Lagrange & Lifting //
const Polynomial lagrange[4] = const Polynomial lagrange[4] =
{ {
...@@ -56,13 +57,13 @@ QuadNedelecBasis::QuadNedelecBasis(void){ ...@@ -56,13 +57,13 @@ QuadNedelecBasis::QuadNedelecBasis(void){
}; };
// Basis // // 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]; basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) // // 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++){ for(size_t e = 0; e < 4; e++){
Polynomial lambda = (lagrange[edgeIdx[s][e][0]] + Polynomial lambda = (lagrange[edgeIdx[s][e][0]] +
...@@ -91,7 +92,7 @@ QuadNedelecBasis::QuadNedelecBasis(void){ ...@@ -91,7 +92,7 @@ QuadNedelecBasis::QuadNedelecBasis(void){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) + Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 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++){ for(size_t i = 0; i < nFunction; i++){
vector<Polynomial>* old; vector<Polynomial>* old;
vector<Polynomial> nxt(3); vector<Polynomial> nxt(3);
...@@ -108,11 +109,10 @@ QuadNedelecBasis::QuadNedelecBasis(void){ ...@@ -108,11 +109,10 @@ QuadNedelecBasis::QuadNedelecBasis(void){
} }
QuadNedelecBasis::~QuadNedelecBasis(void){ QuadNedelecBasis::~QuadNedelecBasis(void){
// ReferenceSpace // const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
delete refSpace;
// Basis // // 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++) for(size_t j = 0; j < nFunction; j++)
delete basis[i][j]; delete basis[i][j];
......
FunctionSpace/QuadNodeBasis.cpp
+ 22
− 22
View file @ 551681df
#include "QuadNodeBasis.h"
#include "QuadReferenceSpace.h"
#include "Legendre.h" #include "Legendre.h"
#include "GmshDefines.h"
#include "ReferenceSpaceManager.h"
#include "QuadNodeBasis.h"
using namespace std; using namespace std;
QuadNodeBasis::QuadNodeBasis(size_t order){ 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 // // Set Basis Type //
this->order = order; this->order = order;
type = 0; type = TYPE_QUA;
dim = 2; dim = 2;
nVertex = 4; nVertex = 4;
...@@ -27,6 +19,15 @@ QuadNodeBasis::QuadNodeBasis(size_t order){ ...@@ -27,6 +19,15 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
nCell = 0; nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell; 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 // // Legendre Polynomial //
Polynomial* legendre = new Polynomial[order]; Polynomial* legendre = new Polynomial[order];
Legendre::integrated(legendre, order); Legendre::integrated(legendre, order);
...@@ -63,13 +64,13 @@ QuadNodeBasis::QuadNodeBasis(size_t order){ ...@@ -63,13 +64,13 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
}; };
// Basis // // 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]; basis[s] = new Polynomial*[nFunction];
// Vertex Based // // 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][0] = new Polynomial(lagrange[0]);
basis[s][1] = new Polynomial(lagrange[1]); basis[s][1] = new Polynomial(lagrange[1]);
basis[s][2] = new Polynomial(lagrange[2]); basis[s][2] = new Polynomial(lagrange[2]);
...@@ -77,7 +78,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){ ...@@ -77,7 +78,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
} }
// Edge Based // // Edge Based //
for(size_t s = 0; s < nRefSpace; s++){ for(size_t s = 0; s < nOrientation; s++){
size_t i = nVertex; size_t i = nVertex;
for(size_t e = 0; e < 4; e++){ for(size_t e = 0; e < 4; e++){
...@@ -100,7 +101,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){ ...@@ -100,7 +101,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
// where f = 0, because triangles // where f = 0, because triangles
// have only ONE face: the face '0' // 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; size_t i = nVertex + nEdge;
for(size_t l1 = 1; l1 < order; l1++){ for(size_t l1 = 1; l1 < order; l1++){
...@@ -131,7 +132,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){ ...@@ -131,7 +132,7 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) + Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 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++){ for(size_t i = 0; i < nFunction; i++){
Polynomial* tmp; Polynomial* tmp;
tmp = basis[s][i]; tmp = basis[s][i];
...@@ -145,11 +146,10 @@ QuadNodeBasis::QuadNodeBasis(size_t order){ ...@@ -145,11 +146,10 @@ QuadNodeBasis::QuadNodeBasis(size_t order){
} }
QuadNodeBasis::~QuadNodeBasis(void){ QuadNodeBasis::~QuadNodeBasis(void){
// ReferenceSpace // const size_t nOrientation = ReferenceSpaceManager::getNOrientation(type);
delete refSpace;
// Basis // // 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++) for(size_t j = 0; j < nFunction; j++)
delete basis[i][j]; delete basis[i][j];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment