diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 48297d445e317702b3b785a0b9a9a79e1ba538f3..f7109a9242d9a4a5b1d74e437c02d68a178ba11a 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -9,12 +9,16 @@
using namespace std;
+const size_t FunctionSpace::nGeoType = 9;
+
FunctionSpace::FunctionSpace(void){
}
FunctionSpace::~FunctionSpace(void){
- // Delete Vector of Basis
- // (FunctionSpace is not responsible for 'true' Basis Deletion)
+ if(self)
+ for(size_t i = 0; i < nGeoType; i++)
+ if(basis[i])
+ delete basis[i];
}
void FunctionSpace::build(GroupOfElement& goe,
@@ -24,40 +28,98 @@ void FunctionSpace::build(GroupOfElement& goe,
this->goe = &goe;
this->mesh = &(goe.getMesh());
- // Get Geo Data (WARNING HOMOGENE MESH REQUIRED)//
- const MElement& element = goe.get(0);
- MElement& myElement =
- const_cast<MElement&>(element);
+ // Check if homogene GoE and get geo type //
+ const vector<size_t>& gType = goe.getTypeStats();
+ const size_t nGType = gType.size();
+ size_t eType = (size_t)(-1);
+
+ for(size_t i = 0; i < nGType; i++)
+ if((eType == (size_t)(-1)) && (gType[i] != 0))
+ eType = i;
+ else if((eType != (size_t)(-1)) && (gType[i] != 0))
+ throw Exception("FunctionSpace needs a uniform mesh");
+
+ // Check if basis is matching type //
+ if(eType != basis.getType())
+ throw Exception("FunctionSpace: Basis is not matching type");
+
+ // Alloc Bases and save given Basis //
+ this->self = false;
+ this->basis.resize(nGeoType, NULL);
+ this->basis[eType] = &basis;
+
+ // Get Number of Function per Entity //
+ // Same for all basis since we have a uniform order
+ MElement& myElement = const_cast<MElement&>(goe.get(0));
int nVertex = myElement.getNumPrimaryVertices();
int nEdge = myElement.getNumEdges();
int nFace = myElement.getNumFaces();
- // Init Struct //
- this->basis.resize(1);
- this->basis[0] = &basis;
-
// Number of *Per* Entity functions //
- // Init
- fPerVertex.resize(1);
- fPerEdge.resize(1);
- fPerFace.resize(1);
- fPerCell.resize(1);
+ fPerVertex = this->basis[eType]->getNVertexBased() / nVertex;
+
+ if(nEdge)
+ fPerEdge = this->basis[eType]->getNEdgeBased() / nEdge;
+ else
+ fPerEdge = 0;
+
+ if(nFace)
+ fPerFace = this->basis[eType]->getNFaceBased() / nFace;
+ else
+ fPerFace = 0;
+
+ fPerCell = this->basis[eType]->getNCellBased();
+
+ // Build Dof //
+ buildDof();
+}
+
+void FunctionSpace::build(GroupOfElement& goe,
+ size_t order, size_t form, string family){
+
+ // Save GroupOfElement & Mesh //
+ this->goe = &goe;
+ this->mesh = &(goe.getMesh());
+
+ // Alloc Basis Vector for all possible geomtrical types //
+ self = true;
+ basis.resize(nGeoType, NULL);
+
+ // Generate Bases //
+ const vector<const MElement*>& element = goe.getAll();
+
+ // Get geomtrical type statistics
+ const vector<size_t>& geoTypeStat = goe.getTypeStats();
+ const size_t nGeoType = geoTypeStat.size();
+
+ // Buils basis for existing geomtrical type
+ for(size_t i = 0; i < nGeoType; i++)
+ if(geoTypeStat[i])
+ basis[i] = BasisGenerator::generate(i, form, order, family);
+
+ // Get Number of Function per Entity //
+ // Same for all basis since we have a uniform order
+ MElement& myElement = const_cast<MElement&>(*element[0]);
+
+ int nVertex = myElement.getNumPrimaryVertices();
+ int nEdge = myElement.getNumEdges();
+ int nFace = myElement.getNumFaces();
+ int type = myElement.getType();
- // Populate
- fPerVertex[0] = this->basis[0]->getNVertexBased() / nVertex;
+ fPerVertex = basis[type]->getNVertexBased() / nVertex;
if(nEdge)
- fPerEdge[0] = this->basis[0]->getNEdgeBased() / nEdge;
+ fPerEdge = this->basis[type]->getNEdgeBased() / nEdge;
else
- fPerEdge[0] = 0;
+ fPerEdge = 0;
if(nFace)
- fPerFace[0] = this->basis[0]->getNFaceBased() / nFace;
+ fPerFace = this->basis[type]->getNFaceBased() / nFace;
else
- fPerFace[0] = 0;
+ fPerFace = 0;
- fPerCell[0] = this->basis[0]->getNCellBased();
+ fPerCell = this->basis[type]->getNCellBased();
// Build Dof //
buildDof();
@@ -117,10 +179,10 @@ vector<Dof> FunctionSpace::getUnorderedKeys(const MElement& elem) const{
// Create Dof //
size_t nDof =
- fPerVertex[0] * nVertex +
- fPerEdge[0] * nEdge +
- fPerFace[0] * nFace +
- fPerCell[0] * nCell;
+ fPerVertex * nVertex +
+ fPerEdge * nEdge +
+ fPerFace * nFace +
+ fPerCell * nCell;
vector<Dof> myDof(nDof);
@@ -128,7 +190,7 @@ vector<Dof> FunctionSpace::getUnorderedKeys(const MElement& elem) const{
// Add Vertex Based Dof //
for(size_t i = 0; i < nVertex; i++){
- for(size_t j = 0; j < fPerVertex[0]; j++){
+ for(size_t j = 0; j < fPerVertex; j++){
myDof[it].setDof(mesh->getGlobalId(*vertex[i]), j);
it++;
}
@@ -136,7 +198,7 @@ vector<Dof> FunctionSpace::getUnorderedKeys(const MElement& elem) const{
// Add Edge Based Dof //
for(size_t i = 0; i < nEdge; i++){
- for(size_t j = 0; j < fPerEdge[0]; j++){
+ for(size_t j = 0; j < fPerEdge; j++){
myDof[it].setDof(mesh->getGlobalId(edge[i]), j);
it++;
}
@@ -144,7 +206,7 @@ vector<Dof> FunctionSpace::getUnorderedKeys(const MElement& elem) const{
// Add Face Based Dof //
for(size_t i = 0; i < nFace; i++){
- for(size_t j = 0; j < fPerFace[0]; j++){
+ for(size_t j = 0; j < fPerFace; j++){
myDof[it].setDof(mesh->getGlobalId(face[i]), j);
it++;
}
@@ -152,7 +214,7 @@ vector<Dof> FunctionSpace::getUnorderedKeys(const MElement& elem) const{
// Add Cell Based Dof //
for(size_t i = 0; i < nCell; i++){
- for(size_t j = 0; j < fPerCell[0]; j++){
+ for(size_t j = 0; j < fPerCell; j++){
myDof[it].setDof(mesh->getGlobalId(element), j);
it++;
}
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 6f27b6c38bbf26da6d5d0e48df13cc56f1afaf0a..0fb748d7ca109fb73f9e932f8c75886cf9afabcd 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/FunctionSpace/FunctionSpace.h
@@ -23,26 +23,29 @@
A FunctionSpace is also responsible for the generation of all
the Dof%s related to its geometrical Support.
-
- @todo
- Allow Hybrid Mesh
*/
class Mesh;
class GroupOfElement;
class FunctionSpace{
+ protected:
+ // Number of possible geomtrical topologies //
+ static const size_t nGeoType;
+
protected:
// Geometry //
const Mesh* mesh;
GroupOfElement* goe;
// Basis //
+ bool self;
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;
+
+ size_t fPerVertex;
+ size_t fPerEdge;
+ size_t fPerFace;
+ size_t fPerCell;
// Scalar Field ? //
bool scalar;
@@ -72,6 +75,8 @@ class FunctionSpace{
FunctionSpace(void);
void build(GroupOfElement& goe, const Basis& basis);
+ void build(GroupOfElement& goe,
+ size_t order, size_t form, std::string family);
void buildDof(void);
};
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 737def17af1d75f02ce3a180691d94b6ac88f662..1e4f3a9c4ce5a1d11c717770c94182a46c03ccb9 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -7,6 +7,11 @@ FunctionSpaceScalar::FunctionSpaceScalar(GroupOfElement& goe,
build(goe, basis);
}
+FunctionSpaceScalar::FunctionSpaceScalar(GroupOfElement& goe, size_t order){
+ scalar = true;
+ build(goe, order, 0, "hierarchical");
+}
+
FunctionSpaceScalar::~FunctionSpaceScalar(void){
// Done by FunctionSpace
}
@@ -15,12 +20,14 @@ double FunctionSpaceScalar::
interpolateInABC(const MElement& element,
const std::vector<double>& coef,
double abc[3]) const{
+ // Element Type //
+ const int eType = element.getType();
// Get Basis Functions //
- const size_t nFun = basis[0]->getNFunction();
+ const size_t nFun = basis[eType]->getNFunction();
fullMatrix<double> fun(nFun, 1);
- basis[0]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
+ basis[eType]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
// Interpolate (in Reference Place) //
double val = 0;
@@ -41,11 +48,14 @@ interpolateDerivativeInABC(const MElement& element,
ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
invJac.invertInPlace();
+ // Element Type //
+ const int eType = element.getType();
+
// Get Basis Functions //
- const size_t nFun = basis[0]->getNFunction();
+ const size_t nFun = basis[eType]->getNFunction();
fullMatrix<double> fun(nFun, 3);
- basis[0]->getDerivative(fun, element, abc[0], abc[1], abc[2]);
+ basis[eType]->getDerivative(fun, element, abc[0], abc[1], abc[2]);
// Interpolate (in Reference Place) //
fullMatrix<double> val(1, 3);
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 66fd1264ecdb55735ac0a79f9c009a22ed61d732..dec1f3938ff5ffa799f72a2ed32f0dbeb1725061 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -11,15 +11,14 @@
This class is a @em Scalar FunctionSpaces.@n
A FunctionSpaceScalar can be @em interpolated.
-
- @todo
- Allow interpolation with multiple basis
*/
class FunctionSpaceScalar : public FunctionSpace{
public:
FunctionSpaceScalar(GroupOfElement& goe, const Basis& basis);
+ FunctionSpaceScalar(GroupOfElement& goe, size_t order);
+
virtual ~FunctionSpaceScalar(void);
double
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index b2c2d6dabd69b93986a57b4cea9af9725a66b252..a1ffd51518acfc7a5d9a708e90722863716c770a 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -21,11 +21,14 @@ interpolateInABC(const MElement& element,
ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], invJac);
invJac.invertInPlace();
+ // Element Type //
+ const int eType = element.getType();
+
// Get Basis Functions //
- const size_t nFun = basis[0]->getNFunction();
+ const size_t nFun = basis[eType]->getNFunction();
fullMatrix<double> fun(nFun, 3);
- basis[0]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
+ basis[eType]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
// Interpolate (in Reference Place) //
fullMatrix<double> val(1, 3);
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 2379d95cd84e5a8052e4037304b217db4a3fc1f1..d230a37a47a8aaf4f8ad92b727b0aa81678ee77b 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -12,9 +12,6 @@
This class is a @em Vectorial FunctionSpaces.@n
A FunctionSpaceVector can be @em interpolated.
-
- @todo
- Allow interpolation with multiple basis
*/