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

moving FunctionSpace back in small_fem: becomes functionspace

parent 311e5afa
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 0 additions and 3409 deletions
FunctionSpace/Basis.cpp deleted 100644 → 0
+ 0
− 7
View file @ 311e5afa
#include "Basis.h"
Basis::Basis(void){
}
Basis::~Basis(void){
}
FunctionSpace/Basis.h deleted 100644 → 0
+ 0
− 329
View file @ 311e5afa
#ifndef _BASIS_H_
#define _BASIS_H_
#include <string>
#include "MElement.h"
/**
@interface Basis
@brief Common Interface of all Basis
This class is the common interface for all Basis.
A Basis is set of linearly independent Polynomial%s
(or Vector%s of Polynomial%s).
The returned matrices are the result of the evaluation
of the basis functions (at N points).
The i-th row of these matrices is always refering to the
i-th function of the basis.
Depending on the nature of the returned value (scalar or vector),
the columns are organized diferently.
For scalar values, we have:
@li The j-th column of the i-th row is
the evaluation of the i-th function at the j-th point
For vectorial values, we have:
@li The j-th column of the i-th row is the first coordinate of
the evaluation of the i-th function at the 3 x j-th point
@li The (j-th + 1) column of the i-th row is the second coordinate of
the evaluation of the i-th function at the 3 x j-th point
@li The (j-th + 2) column of the i-th row is the third coordinate of
the evaluation of the i-th function at the 3 x j-th point
*/
class Basis{
protected:
bool scalar;
bool local;
size_t order;
size_t type;
size_t form;
size_t dim;
size_t nVertex;
size_t nEdge;
size_t nFace;
size_t nCell;
size_t nFunction;
public:
// Destructor //
virtual ~Basis(void);
// Scalar & Local //
bool isScalar(void) const;
bool isLocal(void) const;
// Type of Basis //
size_t getOrder(void) const;
size_t getType(void) const;
size_t getForm(void) const;
size_t getDim(void) const;
// Number of Functions //
size_t getNVertexBased(void) const;
size_t getNEdgeBased(void) const;
size_t getNFaceBased(void) const;
size_t getNCellBased(void) const;
size_t getNFunction(void) const;
// Direct Access to Evaluated Functions //
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const = 0;
virtual void getFunctions(fullMatrix<double>& retValues,
size_t orientation,
double u, double v, double w) const = 0;
virtual void getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const = 0;
// Precompute Functions //
virtual
void preEvaluateFunctions(const fullMatrix<double>& point) const = 0;
virtual
void preEvaluateDerivatives(const fullMatrix<double>& point) const = 0;
// Access to Precomputed Functions //
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(const MElement& element) const = 0;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(const MElement& element) const = 0;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(size_t orientation) const = 0;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(size_t orientation) const = 0;
virtual std::string toString(void) const = 0;
protected:
// 'Constructor' //
Basis(void);
};
/**
@internal
@fn Basis::Basis
Instantiate a new Basis
@endinternal
**
@fn Basis::~Basis
Deletes this Basis
**
@fn Basis::isScalar
@return Returns:
@li true, if this is a scalar Basis
@li false, if this is a vectorial Basis
Scalar basis are sets of Polynomial%s,
and Vectorial basis are sets of Vector%s of Polynomial%s
**
@fn Basis::isLocal
@return Returns:
@li true, if this is a Local Basis
@li false, if this is a Global Basis
**
@fn Basis::getOrder
@return Returns the polynomial order of the Basis
**
@fn Basis::getType
@return Returns the type of the Basis (coherent with gmsh element types):
@li 1 for Points
@li 2 for Lines
@li 3 for Triangles
@li 4 for Quadrangles
@li 5 for Tetrahedra
@li 6 for Pyramids
@li 7 for Prisms
@li 8 for Hexahedra
**
@fn Basis::getForm
@return Returns the diferential form of the Basis:
@li 0 for 0-form
@li 1 for 1-form
@li 2 for 2-form
@li 3 for 3-form
**
@fn Basis::getDim
@return Returns the dimension (1D, 2D or 3D) of the Basis
**
@fn Basis::getNVertexBased
@return Returns the number of Vertex based functions of this Basis
**
@fn Basis::getNEdgeBased
@return Returns the number of Edge based functions of this Basis
**
@fn Basis::getNFaceBased
@return Returns the number of Face based functions of this Basis
**
@fn Basis::getNCellBased
@return Returns the number of Cell based functions of this Basis
**
@fn Basis::getNFunction
@return Returns the number of Polynomial%s
(or Vector%s of Polynomial%s) Functions in this Basis
**
@fn Basis::getFunctions(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
The given matrix is populated with the evaluation of every basis function
at the given coordinates, and for the orientation of the given element
**
@fn Basis::getFunctions(fullMatrix<double>&, size_t, double, double, double) const
@param retValues An allocated matrix
@param orientation A integer
@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
The given matrix is populated with the evaluation of every basis function
at the given coordinates, 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
The given matrix is populated with the evaluation of the derivative
of every basis function at the given coordinates,
and for the orientation of the given element
**
@fn Basis::preEvaluateFunctions
@param point A Matrix with points coordinate
(each line is a point and got 3 coordinates, i.e. 3 rows)
Pre Evaluates every basis function at the given points
**
@fn Basis::preEvaluateDerivatives
@param point A Matrix with points coordinate
(each line is a point and got 3 coordinates, i.e. 3 rows)
Pre Evaluates every basis function derivative at the given points
@li For 0-Form it computes the gradient
@li For 1-Form it computes the curl
@li For 2-Form it computes the divergence
**
@fn Basis::getPreEvaluatedFunctions(size_t) const
@param orientation A natural number defining the reference space orientation
@return Returns a Matrix with the PreEvaluated basis functions
(see Basis::preEvaluateFunctions()), with the given orientation
If no PreEvaluation has been done before calling this function,
an Exception is thrown
**
@fn Basis::getPreEvaluatedDerivatives(size_t) const
@param orientation A natural number defining the reference space orientation
@return Returns a Matrix with the PreEvaluated basis functions derivatives
(see Basis::preEvaluateDerivatives()), with the given orientation
If no PreEvaluation of the gradient has been done
before calling this function, an Exception is thrown
**
@fn Basis::getPreEvaluatedFunctions(const MElement&) const
@param element A MElement
@return Same as Basis::getPreEvaluatedFunctions,
but the orientation is computed with the given element
**
@fn Basis::getPreEvaluatedDerivatives(const MElement&) const
@param element A MElement
@return Same as Basis::getPreEvaluatedFunctions,
but the orientation is computed with the given element
**
@fn Basis::toString
@return Returns a string describing this Basis
*/
//////////////////////
// Inline Functions //
//////////////////////
inline bool Basis::isScalar(void) const{
return scalar;
}
inline bool Basis::isLocal(void) const{
return local;
}
inline size_t Basis::getOrder(void) const{
return order;
}
inline size_t Basis::getType(void) const{
return type;
}
inline size_t Basis::getForm(void) const{
return form;
}
inline size_t Basis::getDim(void) const{
return dim;
}
inline size_t Basis::getNVertexBased(void) const{
return nVertex;
}
inline size_t Basis::getNEdgeBased(void) const{
return nEdge;
}
inline size_t Basis::getNFaceBased(void) const{
return nFace;
}
inline size_t Basis::getNCellBased(void) const{
return nCell;
}
inline size_t Basis::getNFunction(void) const{
return nFunction;
}
#endif
FunctionSpace/BasisGenerator.cpp deleted 100644 → 0
+ 0
− 156
View file @ 311e5afa
#include "BasisGenerator.h"
#include "GmshDefines.h"
#include "Exception.h"
#include "LineNodeBasis.h"
#include "LineEdgeBasis.h"
#include "LineNedelecBasis.h"
#include "LineLagrangeBasis.h"
#include "TriNodeBasis.h"
#include "TriEdgeBasis.h"
#include "TriNedelecBasis.h"
#include "TriLagrangeBasis.h"
#include "QuadNodeBasis.h"
#include "QuadEdgeBasis.h"
#include "QuadNedelecBasis.h"
#include "QuadLagrangeBasis.h"
#include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
#include "TetNedelecBasis.h"
#include "TetLagrangeBasis.h"
#include "HexNodeBasis.h"
#include "HexEdgeBasis.h"
#include "HexLagrangeBasis.h"
BasisGenerator::BasisGenerator(void){
}
BasisGenerator::~BasisGenerator(void){
}
BasisLocal* BasisGenerator::generate(size_t elementType,
size_t basisType,
size_t order,
std::string family){
if(!family.compare(std::string("hierarchical")))
return generateHierarchical(elementType, basisType, order);
else if(!family.compare(std::string("lagrange")))
return generateLagrange(elementType, basisType, order);
else
throw Exception("Unknwown Basis Family: %s", family.c_str());
}
BasisLocal* BasisGenerator::generateHierarchical(size_t elementType,
size_t basisType,
size_t order){
switch(elementType){
case TYPE_LIN: return linHierarchicalGen(basisType, order);
case TYPE_TRI: return triHierarchicalGen(basisType, order);
case TYPE_QUA: return quaHierarchicalGen(basisType, order);
case TYPE_TET: return tetHierarchicalGen(basisType, order);
case TYPE_HEX: return hexHierarchicalGen(basisType, order);
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
elementType);
}
}
BasisLocal* BasisGenerator::generateLagrange(size_t elementType,
size_t basisType,
size_t order){
if(basisType != 0)
throw
Exception("Cannot Have a %d-Form Lagrange Basis (0-Form only)",
basisType);
switch(elementType){
case TYPE_LIN: return new LineLagrangeBasis(order);
case TYPE_TRI: return new TriLagrangeBasis(order);
case TYPE_QUA: return new QuadLagrangeBasis(order);
case TYPE_TET: return new TetLagrangeBasis(order);
case TYPE_HEX: return new HexLagrangeBasis(order);
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
elementType);
}
}
BasisLocal* BasisGenerator::linHierarchicalGen(size_t basisType,
size_t order){
switch(basisType){
case 0: return new LineNodeBasis(order);
case 1:
if (order == 0) return new LineNedelecBasis();
else return new LineEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Lines");
case 3: throw Exception("3-form not implemented on Lines");
default: throw Exception("There is no %d-form", basisType);
}
}
BasisLocal* BasisGenerator::triHierarchicalGen(size_t basisType,
size_t order){
switch(basisType){
case 0: return new TriNodeBasis(order);
case 1:
if (order == 0) return new TriNedelecBasis();
else return new TriEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Triangles");
case 3: throw Exception("3-form not implemented on Triangles");
default: throw Exception("There is no %d-form", basisType);
}
}
BasisLocal* BasisGenerator::quaHierarchicalGen(size_t basisType,
size_t order){
switch(basisType){
case 0: return new QuadNodeBasis(order);
case 1:
if (order == 0) return new QuadNedelecBasis();
else return new QuadEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Quads");
case 3: throw Exception("3-form not implemented on Quads");
default: throw Exception("There is no %d-form", basisType);
}
}
BasisLocal* BasisGenerator::tetHierarchicalGen(size_t basisType,
size_t order){
switch(basisType){
case 0: return new TetNodeBasis(order);
case 1:
if (order == 0) return new TetNedelecBasis();
else return new TetEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Tetrahedrons");
case 3: throw Exception("3-form not implemented on Tetrahedrons");
default: throw Exception("There is no %d-form", basisType);
}
}
BasisLocal* BasisGenerator::hexHierarchicalGen(size_t basisType,
size_t order){
switch(basisType){
case 0: return new HexNodeBasis(order);
//case 1: return new HexEdgeBasis(order);
case 1: throw Exception("1-form not implemented on Hexs");
case 2: throw Exception("2-form not implemented on Hexs");
case 3: throw Exception("3-form not implemented on Hexs");
default: throw Exception("There is no %d-form", basisType);
}
}
FunctionSpace/BasisGenerator.h deleted 100644 → 0
+ 0
− 156
View file @ 311e5afa
#ifndef _BASISGENERATOR_H_
#define _BASISGENERATOR_H_
#include <string>
#include "BasisLocal.h"
/**
@class BasisGenerator
@brief A bunch of class method to generate a Local Basis
A BasisGenerator is a bunch of class methods to generate a local Basis.
*/
class BasisGenerator{
public:
BasisGenerator(void);
~BasisGenerator(void);
static BasisLocal* generate(size_t elementType,
size_t basisType,
size_t order,
std::string family);
static BasisLocal* generate(size_t elementType,
size_t basisType,
size_t order);
static BasisLocal* linHierarchicalGen(size_t basisType, size_t order);
static BasisLocal* triHierarchicalGen(size_t basisType, size_t order);
static BasisLocal* quaHierarchicalGen(size_t basisType, size_t order);
static BasisLocal* tetHierarchicalGen(size_t basisType, size_t order);
static BasisLocal* hexHierarchicalGen(size_t basisType, size_t order);
private:
static BasisLocal* generateHierarchical(size_t elementType,
size_t basisType,
size_t order);
static BasisLocal* generateLagrange(size_t elementType,
size_t basisType,
size_t order);
};
/**
@fn BasisGenerator::BasisGenerator
Instantiates a new BasisGenerator
This class got only class methods, so it is not required to instanciate it.
**
@fn BasisGenerator::~BasisGenerator
Deletes this BasisGenerator
**
@fn BasisGenerator::generate(size_t, size_t, size_t, std::string)
@param elementType The element type of the requested Basis
@param basisType The Basis type
@param order The order or the requested Basis
@param family A string
This method will instanciate the requested Basis of the requested family
@return Returns a pointer to a newly instantiated Basis
Element types are:
@li @c TYPE_LIN for Lines
@li @c TYPE_TRI for Triangles
@li @c TYPE_QUA for Quadrangles
@li @c TYPE_TET for Tetrahedrons
@li @c TYPE_HEX for Hexahedrons
Basis types are:
@li 0 for 0-Form
@li 1 for 1-Form
@li 2 for 2-Form
@li 3 for 3-Form
Families are:
@li hierarchical for
<a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>
Basis functions
@li lagrange for Lagrange's Basis functions
**
@fn BasisGenerator::generate(size_t, size_t, size_t)
@param elementType The element type of the requested Basis
@param basisType The Basis type
@param order The order or the requested Basis
Same as
BasisGenerator::generate(elementType, basisType, order, @em hierarchical)
@return Returns a pointer to a newly instantiated Basis
**
@fn BasisGenerator::linHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
This method will instanciate the requested Basis with a Line support
@return Returns a pointer to a newly instantiated Basis
Basis types are:
@li 0 for 0-Form
@li 1 for 1-Form
@li 2 for 2-Form
@li 3 for 3-Form
The Basis family will be hierarchical
**
@fn BasisGenerator::triHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
Same as BasisGenerator::linHierarchicalGen() but for Triangles
**
@fn BasisGenerator::quaHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
Same as BasisGenerator::linHierarchicalGen() but for Quadrangles
**
@fn BasisGenerator::tetHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
Same as BasisGenerator::linHierarchicalGen() but for Tetrahedra
**
@fn BasisGenerator::hexHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
Same as BasisGenerator::linHierarchicalGen() but for Hexahedra
*/
//////////////////////
// Inline Functions //
//////////////////////
inline BasisLocal* BasisGenerator::generate(size_t elementType,
size_t basisType,
size_t order){
return BasisGenerator::generate(elementType,
basisType,
order,
"hierarchical");
}
#endif
FunctionSpace/BasisHierarchical0Form.cpp deleted 100644 → 0
+ 0
− 250
View file @ 311e5afa
#include <sstream>
#include "Exception.h"
#include "ReferenceSpaceManager.h"
#include "BasisHierarchical0Form.h"
using namespace std;
BasisHierarchical0Form::BasisHierarchical0Form(void){
// Scalar Basis ? //
scalar = true;
// 0-Form //
form = 0;
// Grad Basis //
hasGrad = false;
grad = NULL;
// PreEvaluation //
preEvaluated = false;
preEvaluatedGrad = false;
preEvaluatedFunction = NULL;
preEvaluatedGradFunction = NULL;
}
BasisHierarchical0Form::~BasisHierarchical0Form(void){
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Grad Basis //
if(hasGrad){
for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete grad[i][j];
delete[] grad[i];
}
delete[] grad;
}
// PreEvaluation //
if(preEvaluated){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
if(preEvaluatedGrad){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction;
}
}
void BasisHierarchical0Form::
getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
// Define Orientation //
const size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Matrix //
for(size_t i = 0; i < nFunction; i++)
retValues(i, 0) = basis[orientation][i]->at(u, v, w);
}
void BasisHierarchical0Form::
getFunctions(fullMatrix<double>& retValues,
size_t orientation,
double u, double v, double w) const{
// Fill Matrix //
for(size_t i = 0; i < nFunction; i++)
retValues(i, 0) = basis[orientation][i]->at(u, v, w);
}
void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
// Get Grad //
if(!hasGrad)
getGrad();
// Define Orientation //
const size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Matrix //
for(size_t i = 0; i < nFunction; i++){
retValues(i, 0) = grad[orientation][i]->at(0).at(u, v, w);
retValues(i, 1) = grad[orientation][i]->at(1).at(u, v, w);
retValues(i, 2) = grad[orientation][i]->at(2).at(u, v, w);
}
}
void BasisHierarchical0Form::
preEvaluateFunctions(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Delete if older //
if(preEvaluated){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
// Alloc //
const size_t nPoint = point.size1();
preEvaluatedFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nOrientation; i++)
preEvaluatedFunction[i] = new fullMatrix<double>(nFunction, nPoint);
// Fill Matrix //
for(size_t i = 0; i < nOrientation; i++)
for(size_t j = 0; j < nFunction; j++)
for(size_t k = 0; k < nPoint; k++)
(*preEvaluatedFunction[i])(j, k) = basis[i][j]->at(point(k, 0),
point(k, 1),
point(k, 2));
// PreEvaluated //
preEvaluated = true;
}
void BasisHierarchical0Form::
preEvaluateDerivatives(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Build Grad //
if(!hasGrad)
getGrad();
// Delete if older //
if(preEvaluatedGrad){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedGradFunction[i];
delete[] preEvaluatedGradFunction;
}
// Alloc //
const size_t nPoint = point.size1();
const size_t nPoint3 = nPoint * 3;
preEvaluatedGradFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nOrientation; i++)
preEvaluatedGradFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++){
for(size_t 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>& BasisHierarchical0Form::
getPreEvaluatedFunctions(const MElement& element) const{
return
getPreEvaluatedFunctions(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedDerivatives(const MElement& element) const{
return
getPreEvaluatedDerivatives(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedFunctions(size_t orientation) const{
if(!preEvaluated)
throw Exception
("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[orientation];
}
const fullMatrix<double>& BasisHierarchical0Form::
getPreEvaluatedDerivatives(size_t orientation) const{
if(!preEvaluatedGrad)
throw Exception
("getPreEvaluatedDerivative: gradient has not been preEvaluated");
return *preEvaluatedGradFunction[orientation];
}
void BasisHierarchical0Form::getGrad(void) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Alloc //
grad = new vector<Polynomial>**[nOrientation];
for(size_t s = 0; s < nOrientation; s++)
grad[s] = new vector<Polynomial>*[nFunction];
// Grad //
for(size_t s = 0; s < nOrientation; s++)
for(size_t f = 0 ; f < nFunction; f++)
grad[s][f] = new vector<Polynomial>(basis[s][f]->gradient());
// Has Grad //
hasGrad = true;
}
string BasisHierarchical0Form::toString(void) const{
stringstream stream;
size_t i = 0;
const size_t 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();
}
FunctionSpace/BasisHierarchical0Form.h deleted 100644 → 0
+ 0
− 80
View file @ 311e5afa
#ifndef _BASISHIERARCHICAL0FORM_H_
#define _BASISHIERARCHICAL0FORM_H_
#include "BasisLocal.h"
#include "Polynomial.h"
/**
@interface BasisHierarchical0Form
@brief Interface for hierarchical 0-form local Basis
This is an interface for hierarchical 0-form local Basis.@n
*/
class BasisHierarchical0Form: public BasisLocal{
protected:
// 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 ~BasisHierarchical0Form(void);
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
virtual void getFunctions(fullMatrix<double>& retValues,
size_t orientation,
double u, double v, double w) const;
virtual void getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(size_t orientation) const;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(size_t orientation) const;
virtual std::string toString(void) const;
protected:
BasisHierarchical0Form(void);
private:
void getGrad(void) const;
};
/**
@internal
@fn BasisHierarchical0Form::BasisHierarchical0Form
Instanciates an new BasisHierarchical0Form
@endinternal
**
@fn BasisHierarchical0Form::~BasisHierarchical0Form
Deletes this BasisHierarchical0Form
*/
#endif
FunctionSpace/BasisHierarchical1Form.cpp deleted 100644 → 0
+ 0
− 280
View file @ 311e5afa
#include <sstream>
#include "Exception.h"
#include "ReferenceSpaceManager.h"
#include "BasisHierarchical1Form.h"
using namespace std;
BasisHierarchical1Form::BasisHierarchical1Form(void){
// Scalar Basis ?//
scalar = false;
// 1-Form //
form = 1;
// Curl Basis //
hasCurl = false;
curl = NULL;
// PreEvaluation //
preEvaluated = false;
preEvaluatedCurl = false;
preEvaluatedFunction = NULL;
preEvaluatedCurlFunction = NULL;
}
BasisHierarchical1Form::~BasisHierarchical1Form(void){
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Curl Basis //
if(hasCurl){
for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++)
delete curl[i][j];
delete[] curl[i];
}
delete[] curl;
}
// PreEvaluation //
if(preEvaluated){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
if(preEvaluatedCurl){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction;
}
}
void BasisHierarchical1Form::
getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
// Define Orientation //
size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Vector //
for(size_t i = 0; i < nFunction; i++){
fullVector<double> eval = Polynomial::at(*basis[orientation][i], u, v, w);
retValues(i, 0) = eval(0);
retValues(i, 1) = eval(1);
retValues(i, 2) = eval(2);
}
}
void BasisHierarchical1Form::
getFunctions(fullMatrix<double>& retValues,
size_t orientation,
double u, double v, double w) const{
// Fill Vector //
for(size_t i = 0; i < nFunction; i++){
fullVector<double> eval = Polynomial::at(*basis[orientation][i], u, v, w);
retValues(i, 0) = eval(0);
retValues(i, 1) = eval(1);
retValues(i, 2) = eval(2);
}
}
void BasisHierarchical1Form::getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
// Build Curl //
if(!hasCurl)
getCurl();
// Define Orientation //
const size_t orientation = ReferenceSpaceManager::getOrientation(element);
// Fill Matrix //
for(size_t i = 0; i < nFunction; i++){
retValues(i, 0) = curl[orientation][i]->at(0).at(u, v, w);
retValues(i, 1) = curl[orientation][i]->at(1).at(u, v, w);
retValues(i, 2) = curl[orientation][i]->at(2).at(u, v, w);
}
}
void BasisHierarchical1Form::
preEvaluateFunctions(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Delete if older //
if(preEvaluated){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedFunction[i];
delete[] preEvaluatedFunction;
}
// Alloc //
const size_t nPoint = point.size1();
const size_t nPoint3 = nPoint * 3;
preEvaluatedFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nOrientation; i++)
preEvaluatedFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++){
for(size_t 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 BasisHierarchical1Form::
preEvaluateDerivatives(const fullMatrix<double>& point) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Build Curl //
if(!hasCurl)
getCurl();
// Delete if older //
if(preEvaluatedCurl){
for(size_t i = 0; i < nOrientation; i++)
delete preEvaluatedCurlFunction[i];
delete[] preEvaluatedCurlFunction;
}
// Alloc //
const size_t nPoint = point.size1();
const size_t nPoint3 = nPoint * 3;
preEvaluatedCurlFunction = new fullMatrix<double>*[nOrientation];
for(size_t i = 0; i < nOrientation; i++)
preEvaluatedCurlFunction[i] = new fullMatrix<double>(nFunction, nPoint3);
// Fill Matrix //
fullVector<double> tmp(3);
for(size_t i = 0; i < nOrientation; i++){
for(size_t j = 0; j < nFunction; j++){
for(size_t 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;
}
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedFunctions(const MElement& element) const{
return
getPreEvaluatedFunctions(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedDerivatives(const MElement& element) const{
return
getPreEvaluatedDerivatives(ReferenceSpaceManager::getOrientation(element));
}
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedFunctions(size_t orientation) const{
if(!preEvaluated)
throw
Exception("getPreEvaluatedFunction: function has not been preEvaluated");
return *preEvaluatedFunction[orientation];
}
const fullMatrix<double>& BasisHierarchical1Form::
getPreEvaluatedDerivatives(size_t orientation) const{
if(!preEvaluatedCurl)
throw
Exception("getPreEvaluatedDerivative: curl has not been preEvaluated");
return *preEvaluatedCurlFunction[orientation];
}
void BasisHierarchical1Form::getCurl(void) const{
const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
// Alloc //
curl = new vector<Polynomial>**[nOrientation];
for(size_t s = 0; s < nOrientation; s++)
curl[s] = new vector<Polynomial>*[nFunction];
// Curl //
for(size_t s = 0; s < nOrientation; s++)
for(size_t f = 0 ; f < nFunction; f++)
curl[s][f] = new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
// Has Curl //
hasCurl = true;
}
string BasisHierarchical1Form::toString(void) const{
stringstream stream;
size_t i = 0;
const size_t 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();
}
FunctionSpace/BasisHierarchical1Form.h deleted 100644 → 0
+ 0
− 80
View file @ 311e5afa
#ifndef _BASISHIERARCHICAL1FORM_H_
#define _BASISHIERARCHICAL1FORM_H_
#include "BasisLocal.h"
#include "Polynomial.h"
/**
@interface BasisHierarchical1Form
@brief Interface for hierarchical 1-form local Basis
This is an interface for hierarchical 1-form local Basis.@n
*/
class BasisHierarchical1Form: public BasisLocal{
protected:
// Basis //
std::vector<Polynomial>*** basis;
// Curl Basis //
mutable bool hasCurl;
mutable std::vector<Polynomial>*** curl;
// PreEvaluation //
mutable bool preEvaluated;
mutable bool preEvaluatedCurl;
mutable fullMatrix<double>** preEvaluatedFunction;
mutable fullMatrix<double>** preEvaluatedCurlFunction;
public:
virtual ~BasisHierarchical1Form(void);
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
virtual void getFunctions(fullMatrix<double>& retValues,
size_t orientation,
double u, double v, double w) const;
virtual void getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(size_t orientation) const;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(size_t orientation) const;
virtual std::string toString(void) const;
protected:
BasisHierarchical1Form(void);
private:
void getCurl(void) const;
};
/**
@internal
@fn BasisHierarchical1Form::BasisHierarchical1Form
Instanciates an new BasisHierarchical1Form
@endinternal
**
@fn BasisHierarchical1Form::~BasisHierarchical1Form
Deletes this BasisHierarchical1Form
*/
#endif
FunctionSpace/BasisLagrange.cpp deleted 100644 → 0
+ 0
− 214
View file @ 311e5afa
#include "ReferenceSpaceManager.h"
#include "BasisLagrange.h"
#include "Exception.h"
using namespace std;
BasisLagrange::BasisLagrange(void){
scalar = true;
form = 0;
preEvaluated = false;
preEvaluatedGrad = false;
preEvaluatedFunction = NULL;
preEvaluatedGradFunction = NULL;
}
BasisLagrange::~BasisLagrange(void){
if(preEvaluated)
delete preEvaluatedFunction;
if(preEvaluatedGrad)
delete preEvaluatedGradFunction;
}
void BasisLagrange::
getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
// Fill Matrix //
fullMatrix<double> tmp;
fullMatrix<double> point(1, 3);
point(0, 0) = u;
point(0, 1) = v;
point(0, 2) = w;
lBasis->f(point, tmp);
// Transpose 'tmp': otherwise not coherent with df !!
retValues = tmp.transpose();
// Permute retValues, accordingly to ReferenceSpace
permutation(ReferenceSpaceManager::getOrientation(element), retValues);
}
void BasisLagrange::
getFunctions(fullMatrix<double>& retValues,
size_t orientation,
double u, double v, double w) const{
// Fill Matrix //
fullMatrix<double> tmp;
fullMatrix<double> point(1, 3);
point(0, 0) = u;
point(0, 1) = v;
point(0, 2) = w;
lBasis->f(point, tmp);
// Transpose 'tmp': otherwise not coherent with df !!
retValues = tmp.transpose();
// Permute retValues, accordingly to ReferenceSpace
permutation(orientation, retValues);
}
void BasisLagrange::getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
throw Exception("Not Implemented");
}
void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
// Delete if older //
if(preEvaluated)
delete preEvaluatedFunction;
// Fill Matrix //
fullMatrix<double> tmp;
lBasis->f(point, tmp);
// Transpose 'tmp': otherwise not coherent with df !!
preEvaluatedFunction = new fullMatrix<double>(tmp.transpose());
// PreEvaluated //
preEvaluated = true;
}
void BasisLagrange::
preEvaluateDerivatives(const fullMatrix<double>& point) const{
// Delete if older //
if(preEvaluatedGrad)
delete preEvaluatedGradFunction;
// Alloc //
preEvaluatedGradFunction = new fullMatrix<double>;
// Fill Matrix //
lBasis->df(point, *preEvaluatedGradFunction);
// PreEvaluated //
preEvaluatedGrad = true;
}
const fullMatrix<double>&
BasisLagrange::getPreEvaluatedFunctions(const MElement& element) const{
return *preEvaluatedFunction;
}
const fullMatrix<double>&
BasisLagrange::getPreEvaluatedDerivatives(const MElement& element) const{
return *preEvaluatedGradFunction;
}
const fullMatrix<double>&
BasisLagrange::getPreEvaluatedFunctions(size_t orientation) const{
return *preEvaluatedFunction;
}
const fullMatrix<double>&
BasisLagrange::getPreEvaluatedDerivatives(size_t orientation) const{
return *preEvaluatedGradFunction;
}
vector<double> BasisLagrange::
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceScalar& fSpace){
// Init New Coefs //
const size_t size = lPoint->size1();
const size_t dim = lPoint->size2();
vector<double> newCoef(size);
// Interpolation at Lagrange Points //
for(size_t i = 0; i < size; i++){
fullVector<double> uvw(3);
if(dim > 0)
uvw(0) = (*lPoint)(i, 0);
else
uvw(0) = 0;
if(dim > 1)
uvw(1) = (*lPoint)(i, 1);
else
uvw(1) = 0;
if(dim > 2)
uvw(2) = (*lPoint)(i, 2);
else
uvw(2) = 0;
newCoef[i] = fSpace.interpolateInRefSpace(element,
coef,
uvw);
}
// Return ;
return newCoef;
}
std::vector<double> BasisLagrange::
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceVector& fSpace){
// Init New Coefs //
const size_t size = lPoint->size1();
const size_t dim = lPoint->size2();
vector<double> newCoef(size * 3);
fullVector<double> tmp(3);
// Interpolation at Lagrange Points //
for(size_t i = 0; i < size; i++){
fullVector<double> uvw(3);
if(dim > 0)
uvw(0) = (*lPoint)(i, 0);
else
uvw(0) = 0;
if(dim > 1)
uvw(1) = (*lPoint)(i, 1);
else
uvw(1) = 0;
if(dim > 2)
uvw(2) = (*lPoint)(i, 2);
else
uvw(2) = 0;
tmp = fSpace.interpolateInRefSpace(element,
coef,
uvw);
newCoef[i * 3 + 0] = tmp(0);
newCoef[i * 3 + 1] = tmp(1);
newCoef[i * 3 + 2] = tmp(2);
}
// Return ;
return newCoef;
}
void BasisLagrange::permutation(size_t orientation,
fullMatrix<double>& function) const{
}
std::string BasisLagrange::toString(void) const{
return std::string("BasisLagrange::toString() not Implemented");
}
FunctionSpace/BasisLagrange.h deleted 100644 → 0
+ 0
− 148
View file @ 311e5afa
#ifndef _BASISLAGRANGE_H_
#define _BASISLAGRANGE_H_
#include "BasisLocal.h"
#include "FunctionSpaceScalar.h"
#include "FunctionSpaceVector.h"
#include "fullMatrix.h"
#include "polynomialBasis.h"
/**
@interface BasisLagrange
@brief Interface for Lagrange Basis
This is an interface for Lagrange Basis.
These local scalar Basis allow a coefficient matrix and a monomial matrix
to be consulted.
Coefficients from an other Basis can be projected into a Lagrange Basis.
*/
class BasisLagrange: public BasisLocal{
protected:
polynomialBasis* lBasis; // Lagrange Basis
fullMatrix<double>* lPoint; // Lagrange Points
// PreEvaluation //
mutable bool preEvaluated;
mutable bool preEvaluatedGrad;
mutable fullMatrix<double>* preEvaluatedFunction;
mutable fullMatrix<double>* preEvaluatedGradFunction;
public:
virtual ~BasisLagrange(void);
virtual void getFunctions(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
virtual void getFunctions(fullMatrix<double>& retValues,
size_t orientation,
double u, double v, double w) const;
virtual void getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const;
virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(const MElement& element) const;
virtual const fullMatrix<double>&
getPreEvaluatedFunctions(size_t orientation) const;
virtual const fullMatrix<double>&
getPreEvaluatedDerivatives(size_t orientation) const;
const fullMatrix<double>& getCoefficient(void) const;
const fullMatrix<double>& getMonomial(void) const;
std::vector<double>
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceScalar& fSpace);
std::vector<double>
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceVector& fSpace);
virtual std::string toString(void) const;
protected:
BasisLagrange(void);
private:
void permutation(size_t orientation,
fullMatrix<double>& function) const;
};
/**
@internal
@fn BasisLagrange::BasisLagrange
Instanciates an new BasisLagrange
@endinternal
**
@fn BasisLagrange::~BasisLagrange
Deletes this BasisLagrange
**
@fn BasisLagrange::getCoefficient
@return Returns the Coefficient Matrix
**
@fn BasisLagrange::getMonomial
@return Returns the Monomial Matrix
**
@fn BasisLagrange::project(const MElement&, const std::vector<double>&, const FunctionSpaceScalar&)
@param element A MElement
@param coef A vector of coefficient associated to the given element
@param fSpace The (scalar) FunctionSpace of the given coefficients
@return Returns a vector with the projection of the given coefficients
in this BasisLagrange
**
@fn BasisLagrange::project(const MElement&, const std::vector<double>&, const FunctionSpaceVector&)
@param element A MElement
@param coef A vector of coefficient associated to the given element
@param fSpace The (vectorial) FunctionSpace of the given coefficients
@return Returns a vector with the projection of the given coefficients
in this BasisLagrange
The returned vector has a size three times bigger than coef,
since we need three coefficients with a Lagrange basis,
when we project a vectorial basis on it (one per direction).
The Lagranges coefficients corresponding to the ith entry of coef
are located, in the returned vector, at the index i * 3,
with the following padding:
@li index i * 3 + 0 is the projected coefficient for direction x
@li index i * 3 + 1 is the projected coefficient for direction y
@li index i * 3 + 2 is the projected coefficient for direction z
*/
//////////////////////
// Inline Functions //
//////////////////////
inline const fullMatrix<double>& BasisLagrange::
getCoefficient(void) const{
return lBasis->coefficients;
}
inline const fullMatrix<double>& BasisLagrange::
getMonomial(void) const{
return lBasis->monomials;
}
#endif
FunctionSpace/BasisLocal.cpp deleted 100644 → 0
+ 0
− 8
View file @ 311e5afa
#include "BasisLocal.h"
BasisLocal::BasisLocal(void){
local = true;
}
BasisLocal::~BasisLocal(void){
}
FunctionSpace/BasisLocal.h deleted 100644 → 0
+ 0
− 27
View file @ 311e5afa
#ifndef _BASISLOCAL_H_
#define _BASISLOCAL_H_
#include "Basis.h"
/**
@interface BasisLocal
@brief Common interface of all local Basis
This class is the common interface for all local Basis.
*/
class BasisLocal: public Basis{
public:
//! Deletes this BasisLocal
//!
virtual ~BasisLocal(void);
protected:
//! @internal
//! Instantiate a new BasisLocal
//!
//! @endinternal
BasisLocal(void);
};
#endif
FunctionSpace/CMakeLists.txt deleted 100644 → 0
+ 0
− 64
View file @ 311e5afa
# Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
#
# See the LICENSE.txt file for license information. Please report all
# bugs and problems to the public mailing list <gmsh@onelab.info>.
set(SRC
Polynomial.cpp
Legendre.cpp
PermutationTree.cpp
ReferenceSpace.cpp
LineReferenceSpace.cpp
TriReferenceSpace.cpp
QuadReferenceSpace.cpp
TetReferenceSpace.cpp
HexReferenceSpace.cpp
PyrReferenceSpace.cpp
PriReferenceSpace.cpp
ReferenceSpaceManager.cpp
Basis.cpp
BasisLocal.cpp
BasisGenerator.cpp
BasisLagrange.cpp
BasisHierarchical0Form.cpp
BasisHierarchical1Form.cpp
LineNodeBasis.cpp
LineEdgeBasis.cpp
LineNedelecBasis.cpp
LineLagrangeBasis.cpp
TriNodeBasis.cpp
TriEdgeBasis.cpp
TriNedelecBasis.cpp
TriLagrangeBasis.cpp
QuadNodeBasis.cpp
QuadEdgeBasis.cpp
QuadNedelecBasis.cpp
QuadLagrangeBasis.cpp
TetNodeBasis.cpp
TetEdgeBasis.cpp
TetNedelecBasis.cpp
TetLagrangeBasis.cpp
HexNodeBasis.cpp
HexEdgeBasis.cpp
HexLagrangeBasis.cpp
FunctionSpace.cpp
FunctionSpaceScalar.cpp
FunctionSpaceVector.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
append_gmsh_src(FunctionSpace "${SRC};${HDR}")
## Compatibility with SmallFEM (TO BE REMOVED !!!)
add_sources_in_gmsh(FunctionSpace "${SRC}")
FunctionSpace/FunctionSpace.cpp deleted 100644 → 0
+ 0
− 322
View file @ 311e5afa
#include <sstream>
#include "ReferenceSpaceManager.h"
#include "BasisGenerator.h"
#include "ElementType.h"
#include "FunctionSpace.h"
using namespace std;
const size_t FunctionSpace::nGeoType = 9;
size_t FunctionSpace::nxtOffset = 0;
FunctionSpace::FunctionSpace(void){
// Clear //
dof.clear();
rejected.clear();
// Alloc Basis Vector for all possible geomtrical types //
basis.resize(nGeoType, NULL);
// Alloc Function per Entity //
fPerVertex.resize(nGeoType, 0);
fPerEdge.resize(nGeoType , 0);
fPerFace.resize(nGeoType , 0);
fPerCell.resize(nGeoType , 0);
}
FunctionSpace::~FunctionSpace(void){
for(size_t i = 0; i < nGeoType; i++)
if(basis[i])
delete basis[i];
}
void FunctionSpace::build(const vector<const GroupOfElement*>& goe,
const vector<const GroupOfElement*>& exl,
string family){
// Save Dof type offset //
offset = nxtOffset;
// Save Mesh & Get number of GoE//
const size_t nGoe = goe.size();
const size_t nExl = exl.size();
this->mesh = &(goe[0]->getMesh());
// Build Bases //
for(size_t i = 0; i < nGoe; i++)
getBases(*goe[i], family);
// Build Dof to reject //
for(size_t i = 0; i < nExl; i++)
getRejec(*exl[i]);
// Build Dof //
for(size_t i = 0; i < nGoe; i++)
getMyDof(*goe[i]);
// Next Offset for next FunctionSpace
nxtOffset = findMaxType() + 1;
}
void FunctionSpace::getBases(const GroupOfElement& goe, string family){
// Generate Bases //
const vector<size_t>& geoTypeStat = goe.getTypeStats();
const size_t nGeoType = geoTypeStat.size();
for(size_t i = 0; i < nGeoType; i++)
if(geoTypeStat[i] != 0 && basis[i] == NULL)
basis[i] = BasisGenerator::generate(i, form, order, family);
// Get Number of Function per Entity //
for(size_t i = 0; i < nGeoType; i++){
if(geoTypeStat[i] != 0 && fPerVertex[i] == 0){
int nVertex = ReferenceSpaceManager::getNVertex(i);
int nEdge = ReferenceSpaceManager::getNEdge(i);
int nFace = ReferenceSpaceManager::getNFace(i);
fPerVertex[i] = basis[i]->getNVertexBased() / nVertex;
if(nEdge)
fPerEdge[i] = this->basis[i]->getNEdgeBased() / nEdge;
else
fPerEdge[i] = 0;
if(nFace)
fPerFace[i] = this->basis[i]->getNFaceBased() / nFace;
else
fPerFace[i] = 0;
fPerCell[i] = this->basis[i]->getNCellBased();
}
}
}
void FunctionSpace::getMyDof(const GroupOfElement& goe){
// Get Elements //
const size_t nElement = goe.getNumber();
const vector<const MElement*>& element = goe.getAll();
// Push GroupOfElement into map //
pair<size_t, vector<vector<Dof> > > toInsert;
pair<map<size_t, vector<vector<Dof> > >::iterator, bool> isInserted;
toInsert.first = goe.getId();
toInsert.second = vector<vector<Dof> >(0);
isInserted = dof.insert(toInsert);
if(!isInserted.second)
throw Exception("FunctionSpace: cannot computed Dofs for GroupOfElement %d",
goe.getId());
// Reference & Allocate //
vector<vector<Dof> >& myDof = isInserted.first->second;
myDof.resize(nElement);
// Create Dofs //
for(size_t i = 0; i < nElement; i++)
getKeys(*(element[i]), myDof[i]);
}
void FunctionSpace::getRejec(const GroupOfElement& goe){
// Get Elements //
const size_t nElement = goe.getNumber();
const vector<const MElement*>& element = goe.getAll();
// Allocate //
vector<vector<Dof> > myDof(nElement);
// Create Dofs //
for(size_t i = 0; i < nElement; i++)
getKeys(*(element[i]), myDof[i]);
// Push in rejection map //
for(size_t i = 0; i < nElement; i++){
size_t nDof = myDof[i].size();
for(size_t j = 0; j < nDof; j++)
rejected.insert(myDof[i][j]);
}
}
size_t FunctionSpace::findMaxType(void){
// Maximum type //
size_t maxType = 0;
// Iterate on GroupOfElement Id //
map<size_t, vector<vector<Dof> > >::iterator it = dof.begin();
map<size_t, vector<vector<Dof> > >::iterator end = dof.end();
size_t nElement;
size_t nDof;
size_t type;
for(; it != end; it++){
// Iterate on Elements of this GroupOfElement //
nElement = it->second.size();
for(size_t e = 0; e < nElement; e++){
// Iterate on Dofs of this Element //
nDof = it->second[e].size();
for(size_t d = 0; d < nDof; d++){
// Check if not a RejectedDof
if(it->second[e][d] != Dof::RejectedDof()){
// This Dof Type
type = it->second[e][d].getType();
// If this Dof type is bigger, it becomes the new 'maxType'
if(type > maxType)
maxType = type;
}
}
}
}
// Return maxType //
return maxType;
}
void FunctionSpace::getUnorderedKeys(const MElement& elem,
std::vector<Dof>& dof) const{
// Const_Cast //
MElement& element = const_cast<MElement&>(elem);
// Vertex, Edge & Face //
const size_t nVertex = element.getNumPrimaryVertices();
const size_t nEdge = element.getNumEdges();
const size_t nFace = element.getNumFaces();
vector<MVertex*> vertex(nVertex);
vector<MEdge> edge(nEdge);
vector<MFace> face(nFace);
for(size_t i = 0; i < nVertex; i++)
vertex[i] = element.getVertex(i);
for(size_t i = 0; i < nEdge; i++)
edge[i] = element.getEdge(i);
for(size_t i = 0; i < nFace; i++)
face[i] = element.getFace(i);
// Create Dof //
const size_t type = element.getType();
const size_t fPerVertex = this->fPerVertex[type];
const size_t fPerEdge = this->fPerEdge[type];
const size_t fPerFace = this->fPerFace[type];
const size_t fPerCell = this->fPerCell[type];
size_t it = 0;
size_t nDof =
fPerVertex * nVertex +
fPerEdge * nEdge +
fPerFace * nFace +
fPerCell;
dof.resize(nDof);
// Add Vertex Based Dof //
for(size_t i = 0; i < nVertex; i++){
for(size_t j = 0; j < fPerVertex; j++){
dof[it].setDof(mesh->getGlobalId(*vertex[i]), j + offset);
it++;
}
}
// Add Edge Based Dof //
for(size_t i = 0; i < nEdge; i++){
for(size_t j = 0; j < fPerEdge; j++){
dof[it].setDof(mesh->getGlobalId(edge[i]), j + offset);
it++;
}
}
// Add Face Based Dof //
for(size_t i = 0; i < nFace; i++){
for(size_t j = 0; j < fPerFace; j++){
dof[it].setDof(mesh->getGlobalId(face[i]), j + offset);
it++;
}
}
// Add Cell Based Dof //
for(size_t j = 0; j < fPerCell; j++){
dof[it].setDof(mesh->getGlobalId(element), j + offset);
it++;
}
// Mark rejected keys //
markKeys(dof);
}
void FunctionSpace::markKeys(vector<Dof>& dof) const{
if(rejected.size() != 0){
const size_t nDof = dof.size();
for(size_t i = 0; i < nDof; i++)
if(rejected.count(dof[i]) == 1)
dof[i] = Dof::RejectedDof();
}
}
void FunctionSpace::getKeys(const MElement& elem, std::vector<Dof>& dof) const{
// Const_Cast //
MElement& element = const_cast<MElement&>(elem);
// Create New Element With Permuted Vertices //
// Permutation
const vector<size_t>& vPerm =
ReferenceSpaceManager::getNodeIndexFromABCtoUVW(elem);
// Permuted Vertices
const size_t nVertex = element.getNumPrimaryVertices();
vector<MVertex*> vertex(nVertex);
for(size_t i = 0; i < nVertex; i++)
vertex[i] = element.getVertex(vPerm[i]);
// New Element
MElementFactory factory;
int parentTag = ElementType::ParentTypeFromTag(elem.getTypeForMSH());
int lowOrderTag = ElementType::getTag(parentTag, 1, false);
MElement* permElement = factory.create(lowOrderTag, vertex, element.getNum());
// Get Dofs from permuted Element //
getUnorderedKeys(*permElement, dof);
// Free and Return //
delete permElement;
}
void FunctionSpace::getKeys(const GroupOfElement& goe,
std::set<Dof>& dof) const{
// Get Dofs //
const vector<vector<Dof> >& allDofs = getKeys(goe);
// Add them into map //
const size_t size = allDofs.size();
size_t nDof;
for(size_t i = 0; i < size; i++){
nDof = allDofs[i].size();
for(size_t j = 0; j < nDof; j++)
if(allDofs[i][j] != Dof::RejectedDof())
dof.insert(allDofs[i][j]);
}
}
const std::vector<std::vector<Dof> >&
FunctionSpace::getKeys(const GroupOfElement& goe) const{
// Find vector of Dof from map //
map<size_t, vector<vector<Dof> > >::const_iterator it = dof.find(goe.getId());
if(it == dof.end())
throw Exception("FunctionSpace: cannot find Dofs of GroupOfElement %d",
goe.getId());
// Return vector //
return it->second;
}
FunctionSpace/FunctionSpace.h deleted 100644 → 0
+ 0
− 179
View file @ 311e5afa
#ifndef _FUNCTIONSPACE_H_
#define _FUNCTIONSPACE_H_
#include <map>
#include <vector>
#include "Dof.h"
#include "Mesh.h"
#include "Basis.h"
#include "MElement.h"
#include "Exception.h"
#include "GroupOfElement.h"
/**
@interface FunctionSpace
@brief Common Interface of all Function Spaces
This is the common interface of all Function Spaces.
A FunctionSpace is defined on a support,
which is a collection of MElement%s (GroupOfElement).
Those MElement%s must belong to the same Mesh.
A FunctionSpace is also responsible for the generation of all
the Dof%s related to its geometrical Support.
*/
class Mesh;
class GroupOfElement;
class FunctionSpace{
protected:
// Number of possible geomtrical topologies & Dof Type offset //
static const size_t nGeoType;
static size_t nxtOffset;
protected:
// Offset //
size_t offset;
// Mesh //
const Mesh* mesh;
// Basis //
std::vector<const Basis*> basis;
std::vector<size_t> fPerVertex;
std::vector<size_t> fPerEdge;
std::vector<size_t> fPerFace;
std::vector<size_t> fPerCell;
// Differential From & Order //
bool scalar;
size_t form;
size_t order;
// Rejected Dofs //
std::set<Dof> rejected;
// Dofs //
std::map<size_t, std::vector<std::vector<Dof> > > dof;
public:
virtual ~FunctionSpace(void);
bool isScalar(void) const;
size_t getForm(void) const;
size_t getOrder(void) const;
const Basis& getBasis(const MElement& element) const;
const Basis& getBasis(size_t eType) const;
void getKeys(const MElement& element, std::vector<Dof>& dof) const;
void getKeys(const GroupOfElement& goe, std::set<Dof>& dof) const;
const std::vector<std::vector<Dof> >& getKeys(const GroupOfElement& goe)const;
protected:
FunctionSpace(void);
void build(const std::vector<const GroupOfElement*>& goe,
const std::vector<const GroupOfElement*>& exl,
std::string family);
void getBases(const GroupOfElement& goe, std::string family);
void getMyDof(const GroupOfElement& goe);
void getRejec(const GroupOfElement& goe);
size_t findMaxType(void);
void getUnorderedKeys(const MElement& element, std::vector<Dof>& dof) const;
void markKeys(std::vector<Dof>& dof) const;
};
/**
@internal
@fn FunctionSpace::FunctionSpace
Instatiate a new FunctionSpace
@endinternal
**
@fn FunctionSpace::~FunctionSpace
Deletes this FunctionSpace
**
@fn FunctionSpace::isScalar
@return Returns:
@li true, if this FunstionSpace is scalar
@li flase, otherwise
**
@fn FunctionSpace::getForm
@return Returns this FunctionSpace differential form (0, 1, 2 or 3)
**
@fn FunctionSpace::getOrder
@return Returns this FunctionSpace order
**
@fn FunctionSpace::getBasis(const MElement& element) const
@param element A MElement
@return Returns the Basis associated to the given element
**
@fn FunctionSpace::getBasis(size_t eType) const
@param eType A geomtrical element type tag
@return Returns the Basis associated to the given geomtrical element type tag
**
@fn void FunctionSpace::getKeys(const MElement&,std::vector<Dof>&) const
@param element A MElement
@param dof A vector of Dof%s
Populates the given vector with the Dof%s associated to the given MElement
**
@fn void FunctionSpace::getKeys(const GroupOfElement&, std::set<Dof>&) const
@param goe A GroupOfElement
@param dof A set of Dof%s
Populates the given set with the Dof%s associated to the MElement%s
of the given GroupOfElement
**
@fn const std::vector<std::vector<Dof> >& FunctionSpace::getKeys(const GroupOfElement&) const
@param goe A GroupOfElement
@return Returns a vector of vector of Dof such that:
dof[i][j] is the jth Dof of the ith element of the given GroupOfElement
*/
//////////////////////
// Inline Functions //
//////////////////////
inline bool FunctionSpace::isScalar(void) const{
return scalar;
}
inline size_t FunctionSpace::getForm(void) const{
return form;
}
inline size_t FunctionSpace::getOrder(void) const{
return order;
}
inline const Basis& FunctionSpace::getBasis(const MElement& element) const{
return *basis[element.getType()];
}
inline const Basis& FunctionSpace::getBasis(size_t eType) const{
if(eType >= basis.size())
throw Exception("FunctionSpace::getBasis() -- unknown geometrical type %u",
eType);
return *basis[eType];
}
#endif
FunctionSpace/FunctionSpaceScalar.cpp deleted 100644 → 0
+ 0
− 122
View file @ 311e5afa
#include "Mapper.h"
#include "Exception.h"
#include "FunctionSpaceScalar.h"
using namespace std;
FunctionSpaceScalar::
FunctionSpaceScalar(const vector<const GroupOfElement*>& goe,
const vector<const GroupOfElement*>& exclude,
size_t order, string family){
// Init
init(goe, exclude, order, family);
}
FunctionSpaceScalar::
FunctionSpaceScalar(const vector<const GroupOfElement*>& goe,
size_t order, string family){
// Dummy exclude
vector<const GroupOfElement*> dummy;
// Init
init(goe, dummy, order, family);
}
FunctionSpaceScalar::
FunctionSpaceScalar(const GroupOfElement& goe, size_t order, string family){
// Temp vector
vector<const GroupOfElement*> tmp(1);
tmp[0] = &goe;
// Dummy exclude
vector<const GroupOfElement*> dummy;
// Init
init(tmp, dummy, order, family);
}
FunctionSpaceScalar::~FunctionSpaceScalar(void){
// Done by FunctionSpace
}
void FunctionSpaceScalar::init(const vector<const GroupOfElement*>& goe,
const vector<const GroupOfElement*>& exclude,
size_t order, string family){
// Check
if(order == 0)
throw Exception("FunctionSpaceScalar: "
"Cannot have a order 0 scalar function space");
// Init
this->scalar = true;
this->form = 0;
this->order = order;
// Build FunctionSpace
build(goe, exclude, family);
}
double FunctionSpaceScalar::interpolateInABC(const MElement& element,
const vector<double>& coef,
double abc[3]) const{
// Get Basis Functions //
const Basis& basis = getBasis(element);
const size_t nFun = basis.getNFunction();
fullMatrix<double> fun(nFun, 1);
basis.getFunctions(fun, element, abc[0], abc[1], abc[2]);
// Get All Dofs
vector<Dof> myDof;
getKeys(element, myDof);
// Interpolate (in Reference Place) //
double val = 0;
for(size_t i = 0; i < nFun; i++){
if(myDof[i] != Dof::RejectedDof())
val += fun(i, 0) * coef[i];
}
// Return Interpolated Value //
return val;
}
fullVector<double> FunctionSpaceScalar::
interpolateDerivativeInABC(const MElement& element,
const vector<double>& coef,
double abc[3]) const{
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
invJac.invertInPlace();
// Get Basis Functions //
const Basis& basis = getBasis(element);
const size_t nFun = basis.getNFunction();
fullMatrix<double> fun(nFun, 3);
basis.getDerivative(fun, element, abc[0], abc[1], abc[2]);
// Get All Dofs
vector<Dof> myDof;
getKeys(element, myDof);
// Interpolate (in Reference Place) //
fullMatrix<double> val(1, 3);
val(0, 0) = 0;
val(0, 1) = 0;
val(0, 2) = 0;
for(size_t i = 0; i < nFun; i++){
if(myDof[i] != Dof::RejectedDof()){
val(0, 0) += fun(i, 0) * coef[i];
val(0, 1) += fun(i, 1) * coef[i];
val(0, 2) += fun(i, 2) * coef[i];
}
}
// Return Interpolated Value //
fullVector<double> map(3);
Mapper::hCurl(val, 0, 0, invJac, map);
return map;
}
FunctionSpace/FunctionSpaceScalar.h deleted 100644 → 0
+ 0
− 178
View file @ 311e5afa
#ifndef _FUNCTIONSPACESCALAR_H_
#define _FUNCTIONSPACESCALAR_H_
#include "ReferenceSpaceManager.h"
#include "FunctionSpace.h"
/**
@class FunctionSpaceScalar
@brief A scalar FunctionSpace
This class is a scalar FunctionSpaces.
A FunctionSpaceScalar can be interpolated.
*/
class FunctionSpaceScalar : public FunctionSpace{
public:
FunctionSpaceScalar(const std::vector<const GroupOfElement*>& goe,
const std::vector<const GroupOfElement*>& exclude,
size_t order, std::string family = "hierarchical");
FunctionSpaceScalar(const std::vector<const GroupOfElement*>& goe,
size_t order, std::string family = "hierarchical");
FunctionSpaceScalar(const GroupOfElement& goe,
size_t order, std::string family = "hierarchical");
virtual ~FunctionSpaceScalar(void);
double
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const;
double
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const;
fullVector<double>
interpolateDerivative(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const;
private:
void init(const std::vector<const GroupOfElement*>& goe,
const std::vector<const GroupOfElement*>& exclude,
size_t order, std::string family);
double interpolateInABC(const MElement& element,
const std::vector<double>& coef,
double abc[3]) const;
fullVector<double>
interpolateDerivativeInABC(const MElement& element,
const std::vector<double>& coef,
double abc[3]) const;
};
/**
@fn FunctionSpaceScalar::FunctionSpaceScalar(const std::vector<const GroupOfElement*>&,const std::vector<const GroupOfElement*>&,size_t,std::string)
@param goe A vector of GroupOfElement
@param exclude An other of GroupOfElement
@param order A natural number
@param family A string (defaulted to 'hierarchical')
Instanciates a new FunctionSpaceScalar
on the GroupOfElement%s of 'goe',
with the exception of the GroupOfElement%s of 'exclude',
and with the given order
The instanciated FunctionSpace will use the requested Basis family:
@li If family is equal to 'lagrange' a Lagrange Basis will be used
@li If family is equal to 'hierarchical' a hierarchical Basis will be used
@see See BasisGenerator::generate()
**
@fn FunctionSpaceScalar::FunctionSpaceScalar(const std::vector<const GroupOfElement*>&,size_t,std::string)
@param goe A vector of GroupOfElement
@param order A natural number
@param family A string (defaulted to 'hierarchical')
Same as FunctionSpaceScalar::FunctionSpaceScalar(goe, [], order, family)
**
@fn FunctionSpaceScalar::FunctionSpaceScalar(const GroupOfElement&,size_t,std::string)
@param goe A GroupOfElement
@param order A natural number
@param family A string (defaulted to 'hierarchical')
Same as FunctionSpaceScalar::FunctionSpaceScalar([goe], [], order, family)
**
@fn FunctionSpaceScalar::~FunctionSpaceScalar
Deletes this FunctionSpaceScalar
**
@fn FunctionSpaceScalar::interpolate
@param element The MElement to interpolate on
@param coef The coefficients of the interpolation
@param xyz The coordinate (of a point inside the given element)
of the interpolation in the @em physical space
@return Returns the (scalar) interpolated value
If the given coordinate are not in the given
element @em Bad @em Things may happend
**
@fn FunctionSpaceScalar::interpolateInRefSpace
@param element The MElement to interpolate on
@param coef The coefficients of the interpolation
@param uvw The coordinate (of a point inside the given element)
of the interpolation in the @em reference space
@return Returns the (scalar) interpolated value
If the given coordinate are not in the given
element @em Bad @em Things may happend
**
@fn FunctionSpaceScalar::interpolateDerivative
@param element The MElement to interpolate on
@param coef The coefficients of the interpolation
@param xyz The coordinate (of a point inside the given element)
of the interpolation in the @em physical space
Same as FunctionSpaceScalar::interpolate(element, coef, xyz),
but this method iterpolates the derivative.
*/
//////////////////////
// Inline Functions //
//////////////////////
inline double FunctionSpaceScalar::
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Get ABC Space coordinate //
double abc[3];
ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
inline double FunctionSpaceScalar::
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
// Get ABC Space coordinate //
double abc[3];
ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
inline fullVector<double> FunctionSpaceScalar::
interpolateDerivative(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Get ABC Space coordinate //
double abc[3];
ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
// Interpolate in ABC //
return interpolateDerivativeInABC(element, coef, abc);
}
#endif
FunctionSpace/FunctionSpaceVector.cpp deleted 100644 → 0
+ 0
− 133
View file @ 311e5afa
#include "Mapper.h"
#include "FunctionSpaceVector.h"
using namespace std;
FunctionSpaceVector::
FunctionSpaceVector(const vector<const GroupOfElement*>& goe,
const vector<const GroupOfElement*>& exclude,
size_t order, string family){
// Init
init(goe, exclude, order, family);
}
FunctionSpaceVector::
FunctionSpaceVector(const vector<const GroupOfElement*>& goe,
size_t order, string family){
// Dummy Exclude
vector<const GroupOfElement*> dummy;
// Init
init(goe, dummy, order, family);
}
FunctionSpaceVector::
FunctionSpaceVector(const GroupOfElement& goe,
size_t order, string family){
// Temp vector
vector<const GroupOfElement*> tmp(1);
tmp[0] = &goe;
// Dummy Exclude
vector<const GroupOfElement*> dummy;
// Init
init(tmp, dummy, order, family);
}
FunctionSpaceVector::~FunctionSpaceVector(void){
// Done by FunctionSpace
}
void FunctionSpaceVector::init(const vector<const GroupOfElement*>& goe,
const vector<const GroupOfElement*>& exclude,
size_t order, string family){
// Init
this->scalar = false;
this->form = 1;
this->order = order;
// Build FunctionSpace
build(goe, exclude, family);
}
fullVector<double> FunctionSpaceVector::
interpolateInABC(const MElement& element,
const vector<double>& coef,
double abc[3]) const{
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
invJac.invertInPlace();
// Get Basis Functions //
const Basis& basis = getBasis(element);
const size_t nFun = basis.getNFunction();
fullMatrix<double> fun(nFun, 3);
basis.getFunctions(fun, element, abc[0], abc[1], abc[2]);
// Get All Dofs
vector<Dof> myDof;
getKeys(element, myDof);
// Interpolate (in Reference Place) //
fullMatrix<double> val(1, 3);
val(0, 0) = 0;
val(0, 1) = 0;
val(0, 2) = 0;
for(size_t i = 0; i < nFun; i++){
if(myDof[i] != Dof::RejectedDof()){
val(0, 0) += fun(i, 0) * coef[i];
val(0, 1) += fun(i, 1) * coef[i];
val(0, 2) += fun(i, 2) * coef[i];
}
}
// Return Interpolated Value //
fullVector<double> map(3);
Mapper::hCurl(val, 0, 0, invJac, map);
return map;
}
fullVector<double> FunctionSpaceVector::
interpolateDerivativeInABC(const MElement& element,
const vector<double>& coef,
double abc[3]) const{
// Get Jacobian //
fullMatrix<double> jac(3, 3);
double det =
ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], jac);
// Get Basis Functions //
const Basis& basis = getBasis(element);
const size_t nFun = basis.getNFunction();
fullMatrix<double> fun(nFun, 3);
basis.getDerivative(fun, element, abc[0], abc[1], abc[2]);
// Get All Dofs
vector<Dof> myDof;
getKeys(element, myDof);
// Interpolate (in Reference Place) //
fullMatrix<double> val(1, 3);
val(0, 0) = 0;
val(0, 1) = 0;
val(0, 2) = 0;
for(size_t i = 0; i < nFun; i++){
if(myDof[i] != Dof::RejectedDof()){
val(0, 0) += fun(i, 0) * coef[i];
val(0, 1) += fun(i, 1) * coef[i];
val(0, 2) += fun(i, 2) * coef[i];
}
}
// Return Interpolated Value //
fullVector<double> map(3);
Mapper::hDiv(val, 0, 0, jac, det, map);
return map;
}
FunctionSpace/FunctionSpaceVector.h deleted 100644 → 0
+ 0
− 208
View file @ 311e5afa
#ifndef _FUNCTIONSPACEVECTOR_H_
#define _FUNCTIONSPACEVECTOR_H_
#include "fullMatrix.h"
#include "FunctionSpaceScalar.h"
#include "FunctionSpace.h"
/**
@class FunctionSpaceVector
@brief A vectorial FunctionSpaces
This class is a vectorial FunctionSpaces.
A FunctionSpaceVector can be interpolated.
*/
class FunctionSpaceVector : public FunctionSpace{
public:
FunctionSpaceVector(const std::vector<const GroupOfElement*>& goe,
const std::vector<const GroupOfElement*>& exclude,
size_t order, std::string family = "hierarchical");
FunctionSpaceVector(const std::vector<const GroupOfElement*>& goe,
size_t order, std::string family = "hierarchical");
FunctionSpaceVector(const GroupOfElement& goe,
size_t order, std::string family = "hierarchical");
virtual ~FunctionSpaceVector(void);
fullVector<double>
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const;
fullVector<double>
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const;
fullVector<double>
interpolateDerivative(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const;
fullVector<double>
interpolateDerivativeInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const;
private:
void init(const std::vector<const GroupOfElement*>& goe,
const std::vector<const GroupOfElement*>& exclude,
size_t order, std::string family);
fullVector<double>
interpolateInABC(const MElement& element,
const std::vector<double>& coef,
double abc[3]) const;
fullVector<double>
interpolateDerivativeInABC(const MElement& element,
const std::vector<double>& coef,
double abc[3]) const;
};
/**
@fn FunctionSpaceVector::FunctionSpaceVector(const std::vector<const GroupOfElement*>&,const std::vector<const GroupOfElement*>&,size_t,std::string)
@param goe A vector of GroupOfElement
@param exclude An other of GroupOfElement
@param order A natural number
@param family A string (defaulted to 'hierarchical')
Instanciates a new FunctionSpaceVector
on the GroupOfElement%s of 'goe',
with the exception of the GroupOfElement%s of 'exclude',
and with the given order
The instanciated FunctionSpace will use the requested Basis family:
@li If family is equal to 'lagrange' a Lagrange Basis will be used
@li If family is equal to 'hierarchical' a hierarchical Basis will be used
@see See BasisGenerator::generate()
**
@fn FunctionSpaceVector::FunctionSpaceVector(const std::vector<const GroupOfElement*>&,size_t,std::string)
@param goe A vector of GroupOfElement
@param order A natural number
@param family A string (defaulted to 'hierarchical')
Same as FunctionSpaceVector::FunctionSpaceVector(goe, [], order, family)
**
@fn FunctionSpaceVector::FunctionSpaceVector(const GroupOfElement&,size_t,std::string)
@param goe A GroupOfElement
@param order A natural number
@param family A string (defaulted to 'hierarchical')
Same as FunctionSpaceVector::FunctionSpaceVector([goe], [], order, family)
**
@fn FunctionSpaceVector::~FunctionSpaceVector
Deletes this FunctionSpaceVector
**
@fn FunctionSpaceVector::interpolate
@param element The MElement to interpolate on
@param coef The coefficients of the interpolation
@param xyz The coordinate (of a point inside the given element)
of the interpolation in the @em physical space
@return Returns the (vectorial) interpolated value
If the given coordinate are not in the given
element @em Bad @em Things may happend
**
@fn FunctionSpaceVector::interpolateInRefSpace
@param element The MElement to interpolate on
@param coef The coefficients of the interpolation
@param uvw The coordinate (of a point inside the given element)
of the interpolation in the @em reference space
@return Returns the (vectorial) interpolated value
If the given coordinate are not in the given
element @em Bad @em Things may happend
**
@fn FunctionSpaceVector::interpolateDerivative
@param element The MElement to interpolate on
@param coef The coefficients of the interpolation
@param xyz The coordinate (of a point inside the given element)
of the interpolation in the @em physical space
Same as FunctionSpaceVector::interpolate(element, coef, xyz),
but this method iterpolates the derivative.
**
@fn FunctionSpaceVector::interpolateDerivativeInRefSpace
@param element The MElement to interpolate on
@param coef The coefficients of the interpolation
@param uvw The coordinate (of a point inside the given element)
of the interpolation in the @em reference space
Same as FunctionSpaceVector::interpolateInRefSpace(element, coef, uvw),
but this method iterpolates the derivative.
*/
//////////////////////
// Inline Functions //
//////////////////////
inline fullVector<double> FunctionSpaceVector::
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Get ABC Space coordinate //
double abc[3];
ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
inline fullVector<double> FunctionSpaceVector::
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
// Get ABC Space coordinate //
double abc[3];
ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
inline fullVector<double> FunctionSpaceVector::
interpolateDerivative(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Get ABC Space coordinate //
double abc[3];
ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
// Interpolate in ABC //
return interpolateDerivativeInABC(element, coef, abc);
}
inline fullVector<double> FunctionSpaceVector::
interpolateDerivativeInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
// Get ABC Space coordinate //
double abc[3];
ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
// Interpolate in ABC //
return interpolateDerivativeInABC(element, coef, abc);
}
#endif
FunctionSpace/HexEdgeBasis.cpp deleted 100644 → 0
+ 0
− 468
View file @ 311e5afa
#include <vector>
#include "HexEdgeBasis.h"
#include "Legendre.h"
using namespace std;
HexEdgeBasis::HexEdgeBasis(size_t order){
/*
// Set Basis Type //
this->order = order;
type = 1;
dim = 3;
nVertex = 0;
nEdge = 12 * (order + 1);
nFace = 12 * order * (order + 1);
nCell = 3 * order * order * (order + 1);
nEdgeClosure = 2;
nFaceClosure = 8;
size = 3 * (order + 2) * (order + 2) * (order + 1);
// Alloc Temporary Space //
const int orderPlus = order + 1;
Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus];
Polynomial* xi = new Polynomial[6];
Polynomial* eta = new Polynomial[6];
Polynomial* lambda = new Polynomial[6];
vector<Polynomial>* grXi = new vector<Polynomial>[6];
vector<Polynomial>* grEta = new vector<Polynomial>[6];
Polynomial** iLegendreXi = new Polynomial*[orderPlus];
Polynomial** iLegendreEta = new Polynomial*[orderPlus];
Polynomial** legendreXi = new Polynomial*[orderPlus];
Polynomial** legendreEta = new Polynomial*[orderPlus];
Polynomial* iLegendreX = new Polynomial[orderPlus];
Polynomial* iLegendreY = new Polynomial[orderPlus];
Polynomial* iLegendreZ = new Polynomial[orderPlus];
Polynomial* lagrange = new Polynomial[8];
Polynomial* lagrangeSum = new Polynomial[12];
Polynomial* lifting = new Polynomial[8];
Polynomial* liftingSub = new Polynomial[12];
// Legendre Polynomial //
Legendre::integrated(intLegendre, orderPlus);
Legendre::legendre(legendre, order);
// Points definig Edges //
int edge1[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3};
int edge2[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7};
// Lagrange //
lagrange[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lagrange[1] =
(Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lagrange[2] =
(Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lagrange[3] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lagrange[4] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
Polynomial(1, 0, 0, 1);
lagrange[5] =
(Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
Polynomial(1, 0, 0, 1);
lagrange[6] =
(Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) *
Polynomial(1, 0, 0, 1);
lagrange[7] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) *
Polynomial(1, 0, 0, 1);
// Lagrange Sum //
for(int i = 0; i < 12; i++)
lagrangeSum[i] = lagrange[edge1[i]] + lagrange[edge2[i]];
// Lifting //
lifting[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lifting[1] =
(Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lifting[2] =
(Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lifting[3] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
lifting[4] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
Polynomial(1, 0, 0, 1);
lifting[5] =
(Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
Polynomial(1, 0, 0, 1);
lifting[6] =
(Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0)) +
Polynomial(1, 0, 0, 1);
lifting[7] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0)) +
Polynomial(1, 0, 0, 1);
// Lifting Sub //
for(int i = 0; i < 12; i++)
liftingSub[i] = lifting[edge1[i]] - lifting[edge2[i]];
// Basis (temporary --- *no* const) //
std::vector<std::vector<Polynomial>*> basis(size);
// Edge Based (Nedelec) //
int i = 0; // Function Counter
Polynomial oneHalf(0.5, 0, 0, 0);
for(int e = 0; e < 12; e++){
basis[i] =
new std::vector<Polynomial>((liftingSub[e]).gradient());
basis[i]->at(0).mul(lagrangeSum[e]);
basis[i]->at(1).mul(lagrangeSum[e]);
basis[i]->at(2).mul(lagrangeSum[e]);
basis[i]->at(0).mul(oneHalf);
basis[i]->at(1).mul(oneHalf);
basis[i]->at(2).mul(oneHalf);
i++;
}
// Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 12; e++){
basis[i] =
new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient());
i++;
}
}
// Face Based (Preliminary) //
// Points definig Faces
int face1[6] = {0, 3, 2, 1, 5, 4};
int face2[6] = {1, 7, 6, 0, 6, 7};
int face3[6] = {2, 6, 5, 4, 7, 3};
int face4[6] = {3, 2, 1, 5, 4, 0};
// 'Xi' Functions
for(int f = 0; f < 6; f++)
xi[f] = lifting[face1[f]] - lifting[face2[f]];
// 'Eta' Functions
for(int f = 0; f < 6; f++)
eta[f] = lifting[face1[f]] - lifting[face4[f]];
// 'Lambda' Functions
for(int f = 0; f < 6; f++)
lambda[f] =
lagrange[face1[f]] +
lagrange[face2[f]] +
lagrange[face3[f]] +
lagrange[face4[f]];
// Gradients
for(int f = 0; f < 6; f++){
grXi[f] = xi[f].gradient();
grEta[f] = eta[f].gradient();
}
// Compositions
for(int l = 0; l < orderPlus; l++){
iLegendreXi[l] = new Polynomial[6];
iLegendreEta[l] = new Polynomial[6];
legendreXi[l] = new Polynomial[6];
legendreEta[l] = new Polynomial[6];
for(int f = 0; f < 6; f++){
iLegendreXi[l][f] = intLegendre[l].compose(xi[f]);
iLegendreEta[l][f] = intLegendre[l].compose(eta[f]);
legendreXi[l][f] = legendre[l].compose(xi[f]);
legendreEta[l][f] = legendre[l].compose(eta[f]);
}
}
// Face Based (Type 1) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
for(int f = 0; f < 6; f++){
basis[i] = new std::vector<Polynomial>((iLegendreXi[l1][f] *
iLegendreEta[l2][f] *
lambda[f]).gradient());
i++;
}
}
}
// Face Based (Type 2) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
for(int f = 0; f < 6; f++){
basis[i] = new std::vector<Polynomial>(3);
Polynomial tmp1 =
legendreXi[l1][f] *
iLegendreEta[l2][f];
Polynomial tmp2 =
iLegendreXi[l1][f] *
legendreEta[l2][f];
vector<Polynomial> gr1 = grXi[f];
gr1[0].mul(tmp1);
gr1[1].mul(tmp1);
gr1[2].mul(tmp1);
vector<Polynomial> gr2 = grEta[f];
gr2[0].mul(tmp2);
gr2[1].mul(tmp2);
gr2[2].mul(tmp2);
basis[i]->at(0) = (gr1[0] - gr2[0]) * lambda[f];
basis[i]->at(1) = (gr1[1] - gr2[1]) * lambda[f];
basis[i]->at(2) = (gr1[2] - gr2[2]) * lambda[f];
i++;
}
}
}
// Face Based (Type 3 -- Xi) //
for(int l = 1; l < orderPlus; l++){
for(int f = 0; f < 6; f++){
Polynomial tmp = iLegendreEta[l][f] * lambda[f];
basis[i] = new std::vector<Polynomial>(grXi[f]);
basis[i]->at(0).mul(tmp);
basis[i]->at(1).mul(tmp);
basis[i]->at(2).mul(tmp);
i++;
}
}
// Face Based (Type 3 -- Eta) //
for(int l = 1; l < orderPlus; l++){
for(int f = 0; f < 6; f++){
Polynomial tmp = iLegendreXi[l][f] * lambda[f];
basis[i] = new std::vector<Polynomial>(grEta[f]);
basis[i]->at(0).mul(tmp);
basis[i]->at(1).mul(tmp);
basis[i]->at(2).mul(tmp);
i++;
}
}
// Cell Based (Preliminary) //
Polynomial px = Polynomial(2, 1, 0, 0);
Polynomial py = Polynomial(2, 0, 1, 0);
Polynomial pz = Polynomial(2, 0, 0, 1);
Polynomial zero = Polynomial(0, 0, 0, 0);
px = px - Polynomial(1, 0, 0, 0);
py = py - Polynomial(1, 0, 0, 0);
pz = pz - Polynomial(1, 0, 0, 0);
for(int l = 0; l < orderPlus; l++){
iLegendreX[l] = intLegendre[l].compose(px);
iLegendreY[l] = intLegendre[l].compose(py);
iLegendreZ[l] = intLegendre[l].compose(pz);
}
// Cell Based (Type 1) //
int cellStart = i; // Type one cell base counter
int cellNumber = 0;
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
for(int l3 = 1; l3 < orderPlus; l3++){
basis[i] = new std::vector<Polynomial>((iLegendreX[l1] *
iLegendreY[l2] *
iLegendreZ[l3]).gradient());
i++;
cellNumber++;
}
}
}
// Cell Based (Type 2 -- First Part) //
for(int j = 0; j < cellNumber; j++){
basis[i] = new std::vector<Polynomial>(3);
int off = j + cellStart;
basis[i]->at(0) = basis[off]->at(0);
basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
basis[i]->at(2) = basis[off]->at(2);
i++;
}
// Cell Based (Type 2 -- Second Part) //
for(int j = 0; j < cellNumber; j++){
basis[i] = new std::vector<Polynomial>(3);
int off = j + cellStart;
basis[i]->at(0) = basis[off]->at(0);
basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
basis[i]->at(2) = basis[off]->at(2) * Polynomial(-1, 0, 0, 0);
i++;
}
// Cell Based (Type 3 -- First Part) //
for(int l2 = 1; l2 < orderPlus; l2++){
for(int l3 = 1; l3 < orderPlus; l3++){
basis[i] = new std::vector<Polynomial>(3);
basis[i]->at(0) = iLegendreY[l2] * iLegendreZ[l3];
basis[i]->at(1) = zero;
basis[i]->at(2) = zero;
i++;
}
}
// Cell Based (Type 3 -- Second Part) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l3 = 1; l3 < orderPlus; l3++){
basis[i] = new std::vector<Polynomial>(3);
basis[i]->at(0) = zero;
basis[i]->at(1) = iLegendreX[l1] * iLegendreZ[l3];
basis[i]->at(2) = zero;
i++;
}
}
// Cell Based (Type 3 -- Thrid Part) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
basis[i] = new std::vector<Polynomial>(3);
basis[i]->at(0) = zero;
basis[i]->at(1) = zero;
basis[i]->at(2) = iLegendreX[l1] * iLegendreY[l2];
i++;
}
}
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
delete[] lagrange;
delete[] lagrangeSum;
delete[] lifting;
delete[] liftingSub;
delete[] xi;
delete[] eta;
delete[] lambda;
delete[] grXi;
delete[] grEta;
for(int l = 0; l < orderPlus; l++){
delete[] iLegendreXi[l];
delete[] iLegendreEta[l];
delete[] legendreXi[l];
delete[] legendreEta[l];
}
delete[] iLegendreXi;
delete[] iLegendreEta;
delete[] legendreXi;
delete[] legendreEta;
delete[] iLegendreX;
delete[] iLegendreY;
delete[] iLegendreZ;
// Set Basis //
this->basis = new std::vector<const std::vector<Polynomial>*>
(basis.begin(), basis.end());
*/
}
HexEdgeBasis::~HexEdgeBasis(void){
/*
for(int i = 0; i < size; i++)
delete (*basis)[i];
delete basis;
*/
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment