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

Reorganisation of FunctionSpace -- Basis can be ONLY evaluated + Derivatives...

Reorganisation of FunctionSpace -- Basis can be ONLY evaluated + Derivatives -- New Dof ordering: Iterate on Entity and then on Function Order
parent 7beaa1cd
No related branches found
No related tags found
No related merge requests found
Showing
with 791 additions and 522 deletions
#ifndef _BASIS_H_ #ifndef _BASIS_H_
#define _BASIS_H_ #define _BASIS_H_
#include <string>
#include "ReferenceSpace.h"
/** /**
@interface Basis @interface Basis
@brief Common Interface of all Basis @brief Common Interface of all Basis
...@@ -19,8 +16,6 @@ ...@@ -19,8 +16,6 @@
class Basis{ class Basis{
protected: protected:
const ReferenceSpace* refSpace;
unsigned int nRefSpace;
bool scalar; bool scalar;
unsigned int order; unsigned int order;
...@@ -39,9 +34,6 @@ class Basis{ ...@@ -39,9 +34,6 @@ class Basis{
//! //!
virtual ~Basis(void); virtual ~Basis(void);
//! @return Returns the ReferenceSpace of this Basis
const ReferenceSpace& getReferenceSpace(void) const;
//! @return Returns: //! @return Returns:
//! @li @c true, if this is a //! @li @c true, if this is a
//! @em scalar Basis (BasisScalar) //! @em scalar Basis (BasisScalar)
...@@ -90,9 +82,6 @@ class Basis{ ...@@ -90,9 +82,6 @@ class Basis{
//! in this Basis //! in this Basis
unsigned int getNFunction(void) const; unsigned int getNFunction(void) const;
//! @return Returns the Basis String
virtual std::string toString(void) const = 0;
protected: protected:
//! @internal //! @internal
//! Instantiate a new Basis //! Instantiate a new Basis
...@@ -105,11 +94,6 @@ class Basis{ ...@@ -105,11 +94,6 @@ class Basis{
// Inline Functions // // Inline Functions //
////////////////////// //////////////////////
inline const ReferenceSpace&
Basis::getReferenceSpace(void) const{
return *refSpace;
}
inline bool Basis::isScalar(void) const{ inline bool Basis::isScalar(void) const{
return scalar; return scalar;
} }
......
...@@ -12,7 +12,7 @@ ...@@ -12,7 +12,7 @@
#include "TriNodeBasis.h" #include "TriNodeBasis.h"
#include "TriEdgeBasis.h" #include "TriEdgeBasis.h"
#include "TriNedelecBasis.h" #include "TriNedelecBasis.h"
#include "TriLagrangeBasis.h" //#include "TriLagrangeBasis.h"
#include "TetNodeBasis.h" #include "TetNodeBasis.h"
#include "TetEdgeBasis.h" #include "TetEdgeBasis.h"
...@@ -34,10 +34,10 @@ Basis* BasisGenerator::generate(unsigned int elementType, ...@@ -34,10 +34,10 @@ Basis* BasisGenerator::generate(unsigned int elementType,
if(!family.compare(std::string("zaglmayr"))) if(!family.compare(std::string("zaglmayr")))
return generateZaglmayr(elementType, basisType, order); return generateZaglmayr(elementType, basisType, order);
/*
else if(!family.compare(std::string("lagrange"))) else if(!family.compare(std::string("lagrange")))
return generateLagrange(elementType, basisType, order); return generateLagrange(elementType, basisType, order);
*/
else else
throw Exception("Unknwown Basis Family: %s", family.c_str()); throw Exception("Unknwown Basis Family: %s", family.c_str());
} }
...@@ -56,7 +56,7 @@ Basis* BasisGenerator::generateZaglmayr(unsigned int elementType, ...@@ -56,7 +56,7 @@ Basis* BasisGenerator::generateZaglmayr(unsigned int elementType,
elementType); elementType);
} }
} }
/*
Basis* BasisGenerator::generateLagrange(unsigned int elementType, Basis* BasisGenerator::generateLagrange(unsigned int elementType,
unsigned int basisType, unsigned int basisType,
unsigned int order){ unsigned int order){
...@@ -76,7 +76,7 @@ Basis* BasisGenerator::generateLagrange(unsigned int elementType, ...@@ -76,7 +76,7 @@ Basis* BasisGenerator::generateLagrange(unsigned int elementType,
elementType); elementType);
} }
} }
*/
Basis* BasisGenerator::linZaglmayrGen(unsigned int basisType, Basis* BasisGenerator::linZaglmayrGen(unsigned int basisType,
unsigned int order){ unsigned int order){
switch(basisType){ switch(basisType){
......
...@@ -40,10 +40,11 @@ class BasisGenerator{ ...@@ -40,10 +40,11 @@ class BasisGenerator{
static Basis* generateZaglmayr(unsigned int elementType, static Basis* generateZaglmayr(unsigned int elementType,
unsigned int basisType, unsigned int basisType,
unsigned int order); unsigned int order);
/*
static Basis* generateLagrange(unsigned int elementType, static Basis* generateLagrange(unsigned int elementType,
unsigned int basisType, unsigned int basisType,
unsigned int order); unsigned int order);
*/
}; };
......
#include <sstream>
#include "Exception.h"
#include "BasisHierarchicalScalar.h"
using namespace std;
BasisHierarchicalScalar::BasisHierarchicalScalar(void){
// Grad Basis //
hasGrad = false;
grad = NULL;
// PreEvaluation //
preEvaluated = false;
preEvaluatedGrad = false;
preEvaluatedFunction = NULL;
preEvaluatedGradFunction = NULL;
}
BasisHierarchicalScalar::~BasisHierarchicalScalar(void){
// Grad Basis //
if(hasGrad){
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
delete grad[i][j];
delete[] grad[i];
}
delete[] grad;
}
// PreEvaluation //
if(preEvaluated){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
if(preEvaluatedGrad){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction;
}
}
fullMatrix<double>* BasisHierarchicalScalar::
getFunctions(const MElement& element,
double u, double v, double w) const{
// Alloc Memory //
fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
// Define Permutation //
unsigned int permutation = refSpace->getPermutation(element);
// Fill Matrix //
for(unsigned int i = 0; i < nFunction; i++)
(*values)(i, 0) = basis[permutation][i]->at(u, v, w);
// Return //
return values;
}
fullMatrix<double>* BasisHierarchicalScalar::
getFunctions(unsigned int permutation,
double u, double v, double w) const{
// Alloc Memory //
fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
// Fill Matrix //
for(unsigned int i = 0; i < nFunction; i++)
(*values)(i, 0) = basis[permutation][i]->at(u, v, w);
// Return //
return values;
}
void BasisHierarchicalScalar::
preEvaluateFunctions(const fullMatrix<double>& point) const{
// Delete if older //
if(preEvaluated){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
// Alloc //
const unsigned int nPoint = point.size1();
preEvaluatedFunction = new fullMatrix<double>*[nRefSpace];
for(unsigned int i = 0; i < nRefSpace; i++)
preEvaluatedFunction[i] =
new fullMatrix<double>(nFunction, nPoint);
// Fill Matrix //
for(unsigned int i = 0; i < nRefSpace; i++)
for(unsigned int j = 0; j < nFunction; j++)
for(unsigned int k = 0; k < nPoint; k++)
(*preEvaluatedFunction[i])(j, k) =
basis[i][j]->at(point(k, 0),
point(k, 1),
point(k, 2));
// PreEvaluated //
preEvaluated = true;
}
void BasisHierarchicalScalar::
preEvaluateGradFunctions(const fullMatrix<double>& point) const{
// Build Grad //
if(!hasGrad)
getGrad();
// Delete if older //
if(preEvaluatedGrad){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction;
}
// Alloc //
const unsigned int nPoint = point.size1();
const unsigned int nPoint3 = nPoint * 3;
preEvaluatedGradFunction = new fullMatrix<double>*[nRefSpace];
for(unsigned int i = 0; i < nRefSpace; i++)
preEvaluatedGradFunction[i] =
new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++){
for(unsigned int k = 0; k < nPoint; k++){
tmp = Polynomial::at(*grad[i][j],
point(k, 0),
point(k, 1),
point(k, 2));
(*preEvaluatedGradFunction[i])(j, 3 * k) = tmp(0);
(*preEvaluatedGradFunction[i])(j, 3 * k + 1) = tmp(1);
(*preEvaluatedGradFunction[i])(j, 3 * k + 2) = tmp(2);
}
}
}
// PreEvaluated //
preEvaluatedGrad = true;
}
const fullMatrix<double>& BasisHierarchicalScalar::
getPreEvaluatedFunctions(const MElement& element) const{
if(!preEvaluated)
throw Exception("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[refSpace->getPermutation(element)];
}
const fullMatrix<double>& BasisHierarchicalScalar::
getPreEvaluatedGradFunctions(const MElement& element) const{
if(!preEvaluatedGrad)
throw Exception("getPreEvaluatedGradFunction: gradient has not been preEvaluated");
return *preEvaluatedGradFunction[refSpace->getPermutation(element)];
}
void BasisHierarchicalScalar::getGrad(void) const{
// Alloc //
grad = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
grad[s] = new vector<Polynomial>*[nFunction];
// Grad //
for(unsigned int s = 0; s < nRefSpace; s++)
for(unsigned int f = 0 ; f < nFunction; f++)
grad[s][f] =
new vector<Polynomial>(basis[s][f]->gradient());
// Has Grad //
hasGrad = true;
}
string BasisHierarchicalScalar::toString(void) const{
stringstream stream;
unsigned int i = 0;
const unsigned int refSpace = 0;
stream << "Vertex Based:" << endl;
for(; i < nVertex; i++)
stream << "f(" << i + 1 << ") = "
<< basis[refSpace][i]->toString() << endl;
stream << "Edge Based:" << endl;
for(; i < nVertex + nEdge; i++)
stream << "f(" << i + 1 << ") = "
<< basis[refSpace][i]->toString() << endl;
stream << "Face Based:" << endl;
for(; i < nVertex + nEdge + nFace; i++)
stream << "f(" << i + 1 << ") = "
<< basis[refSpace][i]->toString() << endl;
stream << "Cell Based:" << endl;
for(; i < nVertex + nEdge + nFace + nCell; i++)
stream << "f(" << i + 1 << ") = "
<< basis[refSpace][i]->toString() << endl;
return stream.str();
}
#ifndef _BASISHIERARCHICALSCALAR_H_
#define _BASISHIERARCHICALSCALAR_H_
#include <string>
#include "BasisScalar.h"
#include "Polynomial.h"
#include "ReferenceSpace.h"
class BasisHierarchicalScalar: public BasisScalar{
protected:
// Permutation //
ReferenceSpace* refSpace;
unsigned int nRefSpace;
// Basis //
Polynomial*** basis;
// Grad Basis //
mutable bool hasGrad;
mutable std::vector<Polynomial>*** grad;
// PreEvaluation //
mutable bool preEvaluated;
mutable bool preEvaluatedGrad;
mutable fullMatrix<double>** preEvaluatedFunction;
mutable fullMatrix<double>** preEvaluatedGradFunction;
public:
virtual ~BasisHierarchicalScalar(void);
virtual fullMatrix<double>* getFunctions(const MElement& element,
double u, double v, double w) const;
virtual fullMatrix<double>* getFunctions(unsigned int permutation,
double u, double v, double w) const;
virtual void preEvaluateFunctions (const fullMatrix<double>& point) const;
virtual void preEvaluateGradFunctions(const fullMatrix<double>& point) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedGradFunctions(const MElement& element) const;
std::string toString(void) const;
protected:
BasisHierarchicalScalar(void);
private:
void getGrad(void) const;
};
#endif
#include <sstream>
#include "Exception.h"
#include "BasisHierarchicalVector.h"
using namespace std;
BasisHierarchicalVector::BasisHierarchicalVector(void){
// Curl Basis //
hasCurl = false;
curl = NULL;
// Div Basis //
hasDiv = false;
div = NULL;
// PreEvaluation //
preEvaluated = false;
preEvaluatedCurl = false;
preEvaluatedDiv = false;
preEvaluatedFunction = NULL;
preEvaluatedCurlFunction = NULL;
preEvaluatedDivFunction = NULL;
}
BasisHierarchicalVector::~BasisHierarchicalVector(void){
// Curl Basis //
if(hasCurl){
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
delete curl[i][j];
delete[] curl[i];
}
delete[] curl;
}
// Div Basis //
if(hasDiv){
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
delete div[i][j];
delete[] div[i];
}
delete[] div;
}
// PreEvaluation //
if(preEvaluated){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
if(preEvaluatedCurl){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction;
}
}
fullMatrix<double>* BasisHierarchicalVector::
getFunctions(const MElement& element,
double u, double v, double w) const{
// Alloc Memory //
fullMatrix<double>* values = new fullMatrix<double>(nFunction, 3);
// Define Permutation //
unsigned int permutation = refSpace->getPermutation(element);
// Fill Vector //
for(unsigned int i = 0; i < nFunction; i++){
fullVector<double> eval =
Polynomial::at(*basis[permutation][i], u, v, w);
(*values)(i, 0) = eval(0);
(*values)(i, 1) = eval(1);
(*values)(i, 2) = eval(2);
}
// Return //
return values;
}
fullMatrix<double>* BasisHierarchicalVector::
getFunctions(unsigned int permutation,
double u, double v, double w) const{
// Alloc Memory //
fullMatrix<double>* values = new fullMatrix<double>(nFunction, 3);
// Fill Vector //
for(unsigned int i = 0; i < nFunction; i++){
fullVector<double> eval =
Polynomial::at(*basis[permutation][i], u, v, w);
(*values)(i, 0) = eval(0);
(*values)(i, 1) = eval(1);
(*values)(i, 2) = eval(2);
}
// Return //
return values;
}
void BasisHierarchicalVector::
preEvaluateFunctions(const fullMatrix<double>& point) const{
// Delete if older //
if(preEvaluated){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
// Alloc //
const unsigned int nPoint = point.size1();
const unsigned int nPoint3 = nPoint * 3;
preEvaluatedFunction = new fullMatrix<double>*[nRefSpace];
for(unsigned int i = 0; i < nRefSpace; i++)
preEvaluatedFunction[i] =
new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++){
for(unsigned int k = 0; k < nPoint; k++){
tmp = Polynomial::at(*basis[i][j],
point(k, 0),
point(k, 1),
point(k, 2));
(*preEvaluatedFunction[i])(j, 3 * k) = tmp(0);
(*preEvaluatedFunction[i])(j, 3 * k + 1) = tmp(1);
(*preEvaluatedFunction[i])(j, 3 * k + 2) = tmp(2);
}
}
}
// PreEvaluated //
preEvaluated = true;
}
void BasisHierarchicalVector::
preEvaluateCurlFunctions(const fullMatrix<double>& point) const{
// Build Curl //
if(!hasCurl)
getCurl();
// Delete if older //
if(preEvaluatedCurl){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction;
}
// Alloc //
const unsigned int nPoint = point.size1();
const unsigned int nPoint3 = nPoint * 3;
preEvaluatedCurlFunction = new fullMatrix<double>*[nRefSpace];
for(unsigned int i = 0; i < nRefSpace; i++)
preEvaluatedCurlFunction[i] =
new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++){
for(unsigned int k = 0; k < nPoint; k++){
tmp = Polynomial::at(*curl[i][j],
point(k, 0),
point(k, 1),
point(k, 2));
(*preEvaluatedCurlFunction[i])(j, 3 * k) = tmp(0);
(*preEvaluatedCurlFunction[i])(j, 3 * k + 1) = tmp(1);
(*preEvaluatedCurlFunction[i])(j, 3 * k + 2) = tmp(2);
}
}
}
// PreEvaluated //
preEvaluatedCurl = true;
}
void BasisHierarchicalVector::
preEvaluateDivFunctions(const fullMatrix<double>& point) const{
// Build Div //
if(!hasDiv)
getDiv();
// Delete if older //
if(preEvaluatedDiv){
for(unsigned int i = 0; i < nRefSpace; i++)
delete preEvaluatedDivFunction[i];
delete[] preEvaluatedDivFunction;
}
// Alloc //
const unsigned int nPoint = point.size1();
preEvaluatedDivFunction = new fullMatrix<double>*[nRefSpace];
for(unsigned int i = 0; i < nRefSpace; i++)
preEvaluatedDivFunction[i] =
new fullMatrix<double>(nFunction, nPoint);
// Fill Matrix //
for(unsigned int i = 0; i < nRefSpace; i++)
for(unsigned int j = 0; j < nFunction; j++)
for(unsigned int k = 0; k < nPoint; k++)
(*preEvaluatedDivFunction[i])(j, k) =
div[i][j]->at(point(k, 0),
point(k, 1),
point(k, 2));
// PreEvaluated //
preEvaluatedDiv = true;
}
const fullMatrix<double>& BasisHierarchicalVector::
getPreEvaluatedFunctions(const MElement& element) const{
if(!preEvaluated)
throw Exception("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[refSpace->getPermutation(element)];
}
const fullMatrix<double>& BasisHierarchicalVector::
getPreEvaluatedCurlFunctions(const MElement& element) const{
if(!preEvaluatedCurl)
throw Exception("getPreEvaluatedCurlFunction: curl has not been preEvaluated");
return *preEvaluatedCurlFunction[refSpace->getPermutation(element)];
}
const fullMatrix<double>& BasisHierarchicalVector::
getPreEvaluatedDivFunctions(const MElement& element) const{
if(!preEvaluatedDiv)
throw Exception("getPreEvaluatedDivFunction: divergence has not been preEvaluated");
return *preEvaluatedDivFunction[refSpace->getPermutation(element)];
}
void BasisHierarchicalVector::getCurl(void) const{
// Alloc //
curl = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
curl[s] = new vector<Polynomial>*[nFunction];
// Curl //
for(unsigned int s = 0; s < nRefSpace; s++)
for(unsigned int f = 0 ; f < nFunction; f++)
curl[s][f] =
new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
// Has Curl //
hasCurl = true;
}
void BasisHierarchicalVector::getDiv(void) const{
// Alloc //
div = new Polynomial**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
div[s] = new Polynomial*[nFunction];
// Div //
for(unsigned int s = 0; s < nRefSpace; s++)
for(unsigned int f = 0 ; f < nFunction; f++)
div[s][f] =
new Polynomial(Polynomial::divergence(*basis[s][f]));
// Has Div //
hasDiv = true;
}
string BasisHierarchicalVector::toString(void) const{
stringstream stream;
unsigned int i = 0;
const unsigned int refSpace = 0;
stream << "Vertex Based:" << endl;
for(; i < nVertex; i++)
stream << "f(" << i + 1 << ") = " << endl
<< "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
stream << "Edge Based:" << endl;
for(; i < nVertex + nEdge; i++)
stream << " f(" << i + 1 << ") = " << endl
<< "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
stream << "Face Based:" << endl;
for(; i < nVertex + nEdge + nFace; i++)
stream << " f(" << i + 1 << ") = " << endl
<< "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
stream << "Cell Based:" << endl;
for(; i < nVertex + nEdge + nFace + nCell; i++)
stream << " f(" << i + 1 << ") = " << endl
<< "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
<< "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
return stream.str();
}
#ifndef _BASISHIERARCHICALVECTOR_H_
#define _BASISHIERARCHICALVECTOR_H_
#include <string>
#include "BasisVector.h"
#include "Polynomial.h"
#include "ReferenceSpace.h"
class BasisHierarchicalVector: public BasisVector{
protected:
// Permutation //
ReferenceSpace* refSpace;
unsigned int nRefSpace;
// Basis //
std::vector<Polynomial>*** basis;
// Curl Basis //
mutable bool hasCurl;
mutable std::vector<Polynomial>*** curl;
// Div Basis //
mutable bool hasDiv;
mutable Polynomial*** div;
// PreEvaluation //
mutable bool preEvaluated;
mutable bool preEvaluatedCurl;
mutable bool preEvaluatedDiv;
mutable fullMatrix<double>** preEvaluatedFunction;
mutable fullMatrix<double>** preEvaluatedCurlFunction;
mutable fullMatrix<double>** preEvaluatedDivFunction;
public:
virtual ~BasisHierarchicalVector(void);
virtual fullMatrix<double>* getFunctions(const MElement& element,
double u, double v, double w) const;
virtual fullMatrix<double>* getFunctions(unsigned int permutation,
double u, double v, double w) const;
virtual void preEvaluateFunctions (const fullMatrix<double>& point) const;
virtual void preEvaluateCurlFunctions(const fullMatrix<double>& point) const;
virtual void preEvaluateDivFunctions (const fullMatrix<double>& point) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedCurlFunctions(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedDivFunctions (const MElement& element)const;
std::string toString(void) const;
protected:
BasisHierarchicalVector(void);
private:
void getCurl(void) const;
void getDiv(void) const;
};
#endif
#include "BasisLagrange.h"
BasisLagrange::BasisLagrange(void){
}
BasisLagrange::~BasisLagrange(void){
}
fullMatrix<double>* BasisLagrange::
getFunctions(const MElement& element,
double u, double v, double w) const{
return new fullMatrix<double>(nFunction, 1);
}
fullMatrix<double>* BasisLagrange::getFunctions(unsigned int permutation,
double u, double v, double w) const{
// Alloc Memory //
fullMatrix<double> tmp(nFunction, 1);
fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
// Fill Matrix //
fullMatrix<double> point(1, 3);
point(0, 0) = u;
point(0, 1) = v;
point(0, 2) = w;
basis->f(point, tmp);
*values = tmp.transpose(); // Otherwise not coherent with df !!
// Return //
return values;
}
void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
}
void BasisLagrange::preEvaluateGradFunctions(const fullMatrix<double>& point) const{
}
const fullMatrix<double>&
BasisLagrange::getPreEvaluatedFunctions(const MElement& element) const{
return fullMatrix<double>(nFunction, 1);
}
const fullMatrix<double>&
BasisLagrange::getPreEvaluatedGradFunctions(const MElement& element) const{
return fullMatrix<double>(nFunction, 1);
}
#ifndef _BASISLAGRANGE_H_
#define _BASISLAGRANGE_H_
#include "BasisScalar.h"
#include "polynomialBasis.h"
class BasisLagrange: public BasisScalar{
protected:
polynomialBasis* basis;
public:
virtual ~BasisLagrange(void);
virtual fullMatrix<double>* getFunctions(const MElement& element,
double u, double v, double w) const;
virtual fullMatrix<double>* getFunctions(unsigned int permutation,
double u, double v, double w) const;
virtual void preEvaluateFunctions (const fullMatrix<double>& point) const;
virtual void preEvaluateGradFunctions(const fullMatrix<double>& point) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedGradFunctions(const MElement& element) const;
protected:
BasisLagrange(void);
};
#endif
#include <sstream>
#include "BasisScalar.h" #include "BasisScalar.h"
using namespace std;
BasisScalar::BasisScalar(void){ BasisScalar::BasisScalar(void){
scalar = true; scalar = true;
} }
BasisScalar::~BasisScalar(void){ BasisScalar::~BasisScalar(void){
} }
string BasisScalar::toString(void) const{
stringstream stream;
unsigned int i = 0;
const unsigned int refSpace = 0;
stream << "Vertex Based:" << endl;
for(; i < nVertex; i++)
stream << "f(" << i + 1 << ") = "
<< (*(*basis)[refSpace])[i]->toString() << endl;
stream << "Edge Based:" << endl;
for(; i < nVertex + nEdge; i++)
stream << "f(" << i + 1 << ") = "
<< (*(*basis)[refSpace])[i]->toString() << endl;
stream << "Face Based:" << endl;
for(; i < nVertex + nEdge + nFace; i++)
stream << "f(" << i + 1 << ") = "
<< (*(*basis)[refSpace])[i]->toString() << endl;
stream << "Cell Based:" << endl;
for(; i < nVertex + nEdge + nFace + nCell; i++)
stream << "f(" << i + 1 << ") = "
<< (*(*basis)[refSpace])[i]->toString() << endl;
return stream.str();
}
#ifndef _BASISSCALAR_H_ #ifndef _BASISSCALAR_H_
#define _BASISSCALAR_H_ #define _BASISSCALAR_H_
#include <vector>
#include "Basis.h" #include "Basis.h"
#include "Polynomial.h" #include "MElement.h"
#include "fullMatrix.h"
/** /**
@interface BasisScalar @interface BasisScalar
...@@ -12,42 +12,28 @@ ...@@ -12,42 +12,28 @@
This class is the @em common @em interface for all This class is the @em common @em interface for all
@em scalar Basis.@n @em scalar Basis.@n
@note
A BasisScalar is an @em interface,
so it @em can't be instanciated
*/ */
class BasisScalar: public Basis{ class BasisScalar: public Basis{
protected:
std::vector<std::vector<const Polynomial*>*>* basis;
public: public:
//! Deletes this BasisScalar //! Deletes this BasisScalar
//! //!
virtual ~BasisScalar(void); virtual ~BasisScalar(void);
//! @param refSpace A natural number virtual fullMatrix<double>* getFunctions(const MElement& element,
//! @param i A natural number double u, double v, double w) const = 0;
//! @return Returns the @c i%th @em
//! Basis Function of the @c refSpace%th ReferenceSpace virtual fullMatrix<double>* getFunctions(unsigned int permutation,
const Polynomial& double u, double v, double w) const = 0;
getFunction(unsigned int refSpace, unsigned int i) const;
//! @param refSpace A natural number virtual void preEvaluateFunctions (const fullMatrix<double>& point) const = 0;
//! @return Returns the @em all virtual void preEvaluateGradFunctions(const fullMatrix<double>& point) const = 0;
//! Basis Function of the @c refSpace%th ReferenceSpace
const std::vector<const Polynomial*>&
getFunction(unsigned int refSpace) const;
//! @param element An Element virtual const fullMatrix<double>&
//! @return Returns the @em all getPreEvaluatedFunctions(const MElement& element) const = 0;
//! Basis Function in the @c given element
//! @em ReferenceSpace
const std::vector<const Polynomial*>&
getFunction(const MElement& element) const;
virtual std::string toString(void) const; virtual const fullMatrix<double>&
getPreEvaluatedGradFunctions(const MElement& element) const = 0;
protected: protected:
//! @internal //! @internal
...@@ -57,26 +43,4 @@ class BasisScalar: public Basis{ ...@@ -57,26 +43,4 @@ class BasisScalar: public Basis{
BasisScalar(void); BasisScalar(void);
}; };
//////////////////////
// Inline Function //
//////////////////////
inline
const Polynomial&
BasisScalar::getFunction(unsigned int refSpace, unsigned int i) const{
return *(*(*basis)[refSpace])[i];
}
inline
const std::vector<const Polynomial*>&
BasisScalar::getFunction(unsigned int refSpace) const{
return *(*basis)[refSpace];
}
inline
const std::vector<const Polynomial*>&
BasisScalar::getFunction(const MElement& element) const{
return *(*basis)[refSpace->getPermutation(element)];
}
#endif #endif
#include <sstream>
#include "BasisVector.h" #include "BasisVector.h"
using namespace std;
BasisVector::BasisVector(void){ BasisVector::BasisVector(void){
scalar = false; scalar = false;
} }
BasisVector::~BasisVector(void){ BasisVector::~BasisVector(void){
} }
string BasisVector::toString(void) const{
stringstream stream;
unsigned int i = 0;
const unsigned int refSpace = 0;
stream << "Vertex Based:" << endl;
for(; i < nVertex; i++)
stream << "f(" << i + 1 << ") = " << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
stream << "Edge Based:" << endl;
for(; i < nVertex + nEdge; i++)
stream << " f(" << i + 1 << ") = " << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
stream << "Face Based:" << endl;
for(; i < nVertex + nEdge + nFace; i++)
stream << " f(" << i + 1 << ") = " << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
stream << "Cell Based:" << endl;
for(; i < nVertex + nEdge + nFace + nCell; i++)
stream << " f(" << i + 1 << ") = " << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
<< "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
return stream.str();
}
#ifndef _BASISVECTOR_H_ #ifndef _BASISVECTOR_H_
#define _BASISVECTOR_H_ #define _BASISVECTOR_H_
#include <vector>
#include "Basis.h" #include "Basis.h"
#include "Polynomial.h" #include "MElement.h"
#include "fullMatrix.h"
/** /**
@interface BasisVector @interface BasisVector
...@@ -12,71 +12,39 @@ ...@@ -12,71 +12,39 @@
This class is the @em common @em interface for all This class is the @em common @em interface for all
@em vectorial Basis.@n @em vectorial Basis.@n
@note
A BasisVector is an @em interface,
so it @em can't be instanciated
*/ */
class BasisVector: public Basis{ class BasisVector: public Basis{
protected:
std::vector<std::vector<const std::vector<Polynomial>*>*>* basis;
public: public:
//! Deletes this BasisVector //! Deletes this BasisVector
//! //!
virtual ~BasisVector(void); virtual ~BasisVector(void);
//! @param refSpace A natural number virtual fullMatrix<double>* getFunctions(const MElement& element,
//! @param i A natural number double u, double v, double w) const = 0;
//! @return Returns the @c i%th @em
//! Basis Function of the @c refSpace%th ReferenceSpace virtual fullMatrix<double>* getFunctions(unsigned int permutation,
const std::vector<Polynomial>& double u, double v, double w) const = 0;
getFunction(unsigned int refSpace, unsigned int i) const;
virtual void preEvaluateFunctions (const fullMatrix<double>& point) const = 0;
virtual void preEvaluateCurlFunctions(const fullMatrix<double>& point) const = 0;
virtual void preEvaluateDivFunctions (const fullMatrix<double>& point) const = 0;
//! @param refSpace A natural number virtual const fullMatrix<double>&
//! @return Returns the @em all getPreEvaluatedFunctions(const MElement& element) const = 0;
//! Basis Function of the @c refSpace%th ReferenceSpace
const std::vector<const std::vector<Polynomial>*>&
getFunction(unsigned int refSpace) const;
//! @param element An Element virtual const fullMatrix<double>&
//! @return Returns the @em all getPreEvaluatedCurlFunctions(const MElement& element) const = 0;
//! Basis Function in the @c given element
//! @em ReferenceSpace
const std::vector<const std::vector<Polynomial>*>&
getFunction(const MElement& element) const;
virtual std::string toString(void) const; virtual const fullMatrix<double>&
getPreEvaluatedDivFunctions(const MElement& element) const = 0;
protected: protected:
//! @internal //! @internal
//! Instantiate a new BasisVector //! Instantiates a new BasisVector
//! //!
//! @endinternal //! @endinternal
BasisVector(void); BasisVector(void);
}; };
//////////////////////
// Inline Functions //
//////////////////////
inline
const std::vector<Polynomial>&
BasisVector::getFunction(unsigned int refSpace, unsigned int i) const{
return *(*(*basis)[refSpace])[i];
}
inline
const std::vector<const std::vector<Polynomial>*>&
BasisVector::getFunction(unsigned int refSpace) const{
return *(*basis)[refSpace];
}
inline
const std::vector<const std::vector<Polynomial>*>&
BasisVector::getFunction(const MElement& element) const{
return *(*basis)[refSpace->getPermutation(element)];
}
#endif #endif
...@@ -13,16 +13,13 @@ set(SRC ...@@ -13,16 +13,13 @@ set(SRC
TetReferenceSpace.cpp TetReferenceSpace.cpp
Basis.cpp Basis.cpp
LagrangeBasis.cpp
BasisScalar.cpp BasisScalar.cpp
BasisVector.cpp BasisVector.cpp
BasisGenerator.cpp BasisGenerator.cpp
GradBasis.cpp BasisLagrange.cpp
CurlBasis.cpp BasisHierarchicalScalar.cpp
DivBasis.cpp BasisHierarchicalVector.cpp
EvaluatedBasis.cpp
LineNodeBasis.cpp LineNodeBasis.cpp
LineEdgeBasis.cpp LineEdgeBasis.cpp
......
#include "CurlBasis.h"
using namespace std;
CurlBasis::CurlBasis(const BasisVector& other){
// Reference Space //
refSpace = &other.getReferenceSpace();
nRefSpace = other.getReferenceSpace().getNPermutation();
// Set Basis Type //
order = other.getOrder() - 1;
type = other.getType();
dim = other.getDim();
nVertex = other.getNVertexBased();
nEdge = other.getNEdgeBased();
nFace = other.getNFaceBased();
nCell = other.getNCellBased();
nFunction = other.getNFunction();
// Basis //
basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
for(unsigned int s = 0; s < nRefSpace; s++)
(*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
for(unsigned int s = 0; s < nRefSpace; s++)
for(unsigned int i = 0; i < nFunction; i++)
(*(*basis)[s])[i] =
new vector<Polynomial>
(Polynomial::curl(other.getFunction(s, i)));
}
CurlBasis::~CurlBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
delete (*(*basis)[i])[j];
delete (*basis)[i];
}
delete basis;
}
#ifndef _CURLBASIS_H_
#define _CURLBASIS_H_
#include "BasisVector.h"
/**
@class CurlBasis
@brief The Curl of an other Basis
This class can instantiate a Curl Basis.
*/
class CurlBasis: public BasisVector{
public:
//! @param other An Other Basis
//!
//! Returns a new Basis, which is the curl of
//! the given Basis
CurlBasis(const BasisVector& other);
//! Deletes this Basis
//!
virtual ~CurlBasis(void);
};
#endif
#include "DivBasis.h"
using namespace std;
DivBasis::DivBasis(const BasisVector& other){
// Reference Space //
refSpace = &other.getReferenceSpace();
nRefSpace = other.getReferenceSpace().getNPermutation();
// Set Basis Type //
order = other.getOrder() - 1;
type = other.getType();
dim = other.getDim();
nVertex = other.getNVertexBased();
nEdge = other.getNEdgeBased();
nFace = other.getNFaceBased();
nCell = other.getNCellBased();
nFunction = other.getNFunction();
// Basis //
basis = new vector<vector<const Polynomial*>*>(nRefSpace);
for(unsigned int s = 0; s < nRefSpace; s++)
(*basis)[s] = new vector<const Polynomial*>(nFunction);
for(unsigned int s = 0; s < nRefSpace; s++)
for(unsigned int i = 0; i < nFunction; i++)
(*(*basis)[s])[i] =
new Polynomial
(Polynomial::divergence(other.getFunction(s, i)));
}
DivBasis::~DivBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
delete (*(*basis)[i])[j];
delete (*basis)[i];
}
delete basis;
}
#ifndef _DIVBASIS_H_
#define _DIVBASIS_H_
#include "BasisScalar.h"
#include "BasisVector.h"
/**
@class DivBasis
@brief The Divergence of an other Basis
This class can instantiate a Divergence Basis.
*/
class DivBasis: public BasisScalar{
public:
//! @param other An Other Basis
//!
//! Returns a new Basis, which is the divergence of
//! the given Basis
DivBasis(const BasisVector& other);
//! Deletes this Basis
//!
virtual ~DivBasis(void);
};
#endif
#include "Polynomial.h"
#include "EvaluatedBasis.h"
using namespace std;
EvaluatedBasis::
EvaluatedBasis(const BasisScalar& basis,
const fullMatrix<double>& point){
// Data //
refSpace = &basis.getReferenceSpace();
nRefSpace = refSpace->getNPermutation();
scalar = true;
nFunction = basis.getNFunction();
nPoint = point.size1();
// Alloc //
eBasis = new vector<fullMatrix<double>*>(nRefSpace);
for(unsigned int i = 0; i < nRefSpace; i++)
(*eBasis)[i] = new fullMatrix<double>(nFunction, nPoint);
// Evaluate //
for(unsigned int i = 0; i < nRefSpace; i++)
for(unsigned int j = 0; j < nFunction; j++)
for(unsigned int k = 0; k < nPoint; k++)
(*(*eBasis)[i])(j, k) =
basis.getFunction(i, j).at(point(k, 0),
point(k, 1),
point(k, 2));
}
EvaluatedBasis::
EvaluatedBasis(const BasisVector& basis,
const fullMatrix<double>& point){
// Data //
refSpace = &basis.getReferenceSpace();
nRefSpace = refSpace->getNPermutation();
scalar = false;
nFunction = basis.getNFunction();
nPoint = point.size1();
// Alloc //
eBasis = new vector<fullMatrix<double>*>(nRefSpace);
// WARNING Each Evaluation Got *3* Component !
const unsigned int nPointThree = nPoint * 3;
for(unsigned int i = 0; i < nRefSpace; i++)
(*eBasis)[i] = new fullMatrix<double>(nFunction, nPointThree);
// Evaluate //
fullVector<double> tmp(3);
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++){
for(unsigned int k = 0; k < nPoint; k++){
tmp = Polynomial::at(basis.getFunction(i, j),
point(k, 0),
point(k, 1),
point(k, 2));
(*(*eBasis)[i])(j, 3 * k) = tmp(0);
(*(*eBasis)[i])(j, 3 * k + 1) = tmp(1);
(*(*eBasis)[i])(j, 3 * k + 2) = tmp(2);
}
}
}
}
EvaluatedBasis::~EvaluatedBasis(void){
for(unsigned int i = 0; i < nRefSpace; i++)
delete (*eBasis)[i];
delete eBasis;
}
#ifndef _EVALUATEDBASIS_H_
#define _EVALUATEDBASIS_H_
#include <vector>
#include "fullMatrix.h"
#include "BasisScalar.h"
#include "BasisVector.h"
#include "MElement.h"
#include "ReferenceSpace.h"
/**
@class EvaluatedBasis
@brief A Basis that has been Evaluated
This class is A basis that has been Evaluated at
some given points.@n
*/
class EvaluatedBasis{
private:
const ReferenceSpace* refSpace;
unsigned int nRefSpace;
bool scalar;
unsigned int nFunction;
unsigned int nPoint;
std::vector<fullMatrix<double>*>* eBasis;
public:
//! @param basis the Basis to Evaluate
//! @param point the Evaluation Points
//!
//! Instanciates the requested EvaluatedBasisScalar
//!
EvaluatedBasis(const BasisScalar& basis,
const fullMatrix<double>& point);
//! @param basis the Basis to Evaluate
//! @param point the Evaluation Points
//!
//! Instanciates the requested EvaluatedBasisVector
//!
EvaluatedBasis(const BasisVector& basis,
const fullMatrix<double>& point);
//! Deletes this EvaluatedBasis
//!
~EvaluatedBasis(void);
//! @return Returns:
//! @li @c true, if the evaluated
//! values are @em scalar
//! @li @c false, otherwise
bool isScalar(void) const;
//! @return Returns the number of
//! evaluated @em Functions
unsigned int getNFunction(void) const;
//! @return Returns the number of
//! evaluation @em Points
unsigned int getNPoint(void) const;
//! @param refSpace A natural number
//! @return Returns the evaluation of all Basis Functions
//! (at every evaluation Points)
//! for the @c refSpace%th @em ReferenceSpace
const fullMatrix<double>&
getEvaluation(unsigned int refSpace) const;
//! @param element A MElement
//! @return Returns the evaluation of all Basis Functions
//! (at every evaluation Points)
//! in the ReferenceSpace of the given element
const fullMatrix<double>&
getEvaluation(const MElement& element) const;
};
//////////////////////
// Inline Function //
//////////////////////
inline bool EvaluatedBasis::isScalar(void) const{
return scalar;
}
inline unsigned int EvaluatedBasis::getNFunction(void) const{
return nFunction;
}
inline unsigned int EvaluatedBasis::getNPoint(void) const{
return nPoint;
}
inline const fullMatrix<double>&
EvaluatedBasis::getEvaluation(unsigned int refSpace) const{
return *(*eBasis)[refSpace];
}
inline const fullMatrix<double>&
EvaluatedBasis::getEvaluation(const MElement& element) const{
return *(*eBasis)[refSpace->getPermutation(element)];
}
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment