From d651b7cd4d79c2ba98f59e1af611eee9800fc244 Mon Sep 17 00:00:00 2001
From: Eric Bechet <eric.bechet@ulg.ac.be>
Date: Thu, 10 Dec 2009 10:38:33 +0000
Subject: [PATCH] Messed up with function spaces / working example in linear
elasiticity
---
Solver/CMakeLists.txt | 1 +
Solver/SElement.cpp | 4 +-
Solver/dofManager.h | 14 +--
Solver/elasticitySolver.cpp | 170 +++++++++++++++++++++++++++++
Solver/elasticitySolver.h | 14 ++-
Solver/femTerm.h | 38 +++++++
Solver/functionSpace.cpp | 14 +++
Solver/functionSpace.h | 206 +++++++++++++++++++++++++++++++-----
8 files changed, 427 insertions(+), 34 deletions(-)
create mode 100644 Solver/functionSpace.cpp
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index b230ce4275..add58f8e9a 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -21,6 +21,7 @@ set(SRC
function.cpp
functionDerivator.cpp
luaFunction.cpp
+ functionSpace.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h)
diff --git a/Solver/SElement.cpp b/Solver/SElement.cpp
index 8c9edcfd34..1c7a02354f 100644
--- a/Solver/SElement.cpp
+++ b/Solver/SElement.cpp
@@ -42,7 +42,7 @@ void SElement::gradNodalFunctions (double u, double v, double w, double invjac[3
_e->getShapeFunctions(u, v, w, sf);
// FIXME : re-use sf for computing coordinates
_e->pnt(u, v, w, p);
- double E = (*_enrichement_s)(p.x(), p.y(), p.z());
+ double E = (*_enrichement)(p.x(), p.y(), p.z());
double dEdx, dEdy, dEdz;
_enrichement_s->gradient(p.x(), p.y(), p.z(), dEdx, dEdy, dEdz);
for (int i = 0; i < N; i++){
@@ -62,7 +62,7 @@ void SElement::nodalFunctions (double u, double v, double w, double s[],
SPoint3 p;
// FIXME : re-use s for computing coordinates
_e->pnt(u, v, w, p);
- double E = (*_enrichement_s)(p.x(), p.y(), p.z());
+ double E = (*_enrichement)(p.x(), p.y(), p.z());
for (int i = 0; i < N; i++){
s[i] *= E;
}
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 45802227b5..5ffcbf3b69 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -23,11 +23,11 @@ class Dof{
Dof(long int entity, int type) : _entity(entity), _type(type) {}
inline long int getEntity() const { return _entity; }
inline int getType() const { return _type; }
- static int createTypeWithTwoInts(int i1, int i2)
+ inline static int createTypeWithTwoInts(int i1, int i2)
{
return i1 + 10000 * i2;
}
- static void getTwoIntsFromType(int t, int &i1, int &i2)
+ inline static void getTwoIntsFromType(int t, int &i1, int &i2)
{
i1 = t % 10000;
i2 = t / 10000;
@@ -50,18 +50,18 @@ template<class T> struct dofTraits
typedef T MatType;
};
-template<> struct dofTraits<fullVector<double> >
+template<class T> struct dofTraits<fullVector<T> >
{
- typedef fullVector<double> VecType;
- typedef fullMatrix<double> MatType;
+ typedef fullVector<T> VecType;
+ typedef fullMatrix<T> MatType;
};
-
+/*
template<> struct dofTraits<fullVector<std::complex<double> > >
{
typedef fullVector<std::complex<double> > VecType;
typedef fullMatrix<std::complex<double> > MatType;
};
-
+*/
template<class T>
class DofAffineConstraint{
public:
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 4358067e65..b2f6a33203 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -13,6 +13,8 @@
#include "linearSystemPETSc.h"
#include "linearSystemGMM.h"
#include "Numeric.h"
+#include "functionSpace.h"
+#include "terms.h"
#if defined(HAVE_POST)
#include "PView.h"
@@ -321,3 +323,171 @@ void elasticitySolver::solve()
// solving
lsys->systemSolve();
}
+
+
+
+
+void MyelasticitySolver::solve()
+{
+#if defined(HAVE_TAUCS)
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+ linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#else
+ linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+ lsys->setNoisy(2);
+#endif
+ pAssembler = new dofManager<double>(lsys);
+
+ std::map<int, std::vector<GEntity*> > groups[4];
+ pModel->getPhysicalGroups(groups);
+
+ // we first do all fixations. the behavior of the dofManager is to
+ // give priority to fixations : when a dof is fixed, it cannot be
+ // numbered afterwards
+
+ for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
+ it != nodalDisplacements.end(); ++it){
+ MVertex *v = pModel->getMeshVertexByTag(it->first.first);
+ pAssembler->fixVertex(v, it->first.second, _tag, it->second);
+ printf("-- Fixing node %3d comp %1d to %8.5f\n",
+ it->first.first, it->first.second, it->second);
+ }
+
+ for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
+ it != edgeDisplacements.end(); ++it){
+ DummyfemTerm El(pModel);
+ El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag,
+ simpleFunction<double>(it->second), *pAssembler);
+ printf("-- Fixing edge %3d comp %1d to %8.5f\n",
+ it->first.first, it->first.second, it->second);
+ }
+
+ for (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
+ it != faceDisplacements.end(); ++it){
+ DummyfemTerm El(pModel);
+ El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag,
+ simpleFunction<double>(it->second), *pAssembler);
+ printf("-- Fixing face %3d comp %1d to %8.5f\n",
+ it->first.first, it->first.second, it->second);
+ }
+
+ // we number the nodes : when a node is numbered, it cannot be numbered
+ // again with another number. so we loop over elements
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
+ groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+ for ( ; it != elasticFields[i].g->end(); ++it){
+ SElement se(*it);
+ for (int j = 0; j < se.getNumNodalShapeFunctions(); ++j){
+ pAssembler->numberVertex(se.getVertex(j), 0, elasticFields[i]._tag);
+ pAssembler->numberVertex(se.getVertex(j), 1, elasticFields[i]._tag);
+ pAssembler->numberVertex(se.getVertex(j), 2, elasticFields[i]._tag);
+ }
+ }
+ }
+
+ // Now we start the assembly process
+ // First build the force vector
+ // Nodes at geometrical vertices
+ for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
+ it != nodalForces.end(); ++it){
+ int iVertex = it->first;
+ SVector3 f = it->second;
+ std::vector<GEntity*> ent = groups[0][iVertex];
+ for (unsigned int i = 0; i < ent.size(); i++){
+ pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
+ pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
+ pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
+ printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
+ ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
+ }
+ }
+
+ // line forces
+ for (std::map<int, SVector3 >::iterator it = lineForces.begin();
+ it != lineForces.end(); ++it){
+ DummyfemTerm El(pModel);
+ int iEdge = it->first;
+ SVector3 f = it->second;
+ El.neumannNodalBC(iEdge, 1, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
+ El.neumannNodalBC(iEdge, 1, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
+ El.neumannNodalBC(iEdge, 1, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+ printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
+ }
+
+ // face forces
+ for (std::map<int, SVector3 >::iterator it = faceForces.begin();
+ it != faceForces.end(); ++it){
+ DummyfemTerm El(pModel);
+ int iFace = it->first;
+ SVector3 f = it->second;
+ El.neumannNodalBC(iFace, 2, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
+ El.neumannNodalBC(iFace, 2, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
+ El.neumannNodalBC(iFace, 2, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+ printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
+ }
+
+ // volume forces
+ for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
+ it != volumeForces.end(); ++it){
+ DummyfemTerm El(pModel);
+ int iVolume = it->first;
+ SVector3 f = it->second;
+// El.setVector(f);
+ printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
+ std::vector<GEntity*> ent = groups[_dim][iVolume];
+ for (unsigned int i = 0; i < ent.size(); i++){
+ // to do
+ // El.addToRightHandSide(*pAssembler, ent[i]);
+ }
+ }
+
+ // assembling the stifness matrix
+ for (unsigned int i = 0; i < elasticFields.size(); i++){
+ SElement::setShapeEnrichement(elasticFields[i]._enrichment);
+ for (unsigned int j = 0; j < elasticFields.size(); j++){
+ elasticityTerm El(0, elasticFields[i]._E, elasticFields[i]._nu,
+ elasticFields[i]._tag, elasticFields[j]._tag);
+ SElement::setTestEnrichement(elasticFields[j]._enrichment);
+ // printf("%d %d\n",elasticFields[i]._tag,elasticFields[j]._tag);
+ El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);
+ }
+ }
+
+/* VectorLagrangeFunctionSpace L1(Dof::createTypeWithTwoInts(0,_tag),VectorLagrangeFunctionSpace::VECTOR_X);
+ VectorLagrangeFunctionSpace L2(Dof::createTypeWithTwoInts(1,_tag),VectorLagrangeFunctionSpace::VECTOR_Y);
+ VectorLagrangeFunctionSpace L3(Dof::createTypeWithTwoInts(2,_tag),VectorLagrangeFunctionSpace::VECTOR_Z);
+ CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType> P123(L1,L2,L3);*/
+ VectorLagrangeFunctionSpace P123(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+
+ for (unsigned int i = 0; i < elasticFields.size(); i++)
+ {
+ DummyfemTerm El(pModel);
+// ElasticTerm<CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType>,CompositeFunctionSpace<VectorLagrangeFunctionSpace::ValType> > Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
+ ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
+ fullMatrix<double> localMatrix(12,12);
+ std::vector<Dof> R;R.reserve(100);
+ for ( groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end() ; ++it)
+ {
+ MElement *e = *it;
+ R.clear();
+ int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
+ int npts=0;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ Eterm.get(e,npts,GP,localMatrix);
+ El.addToMatrix(*pAssembler,localMatrix,R);
+ }
+ }
+
+ printf("-- done assembling!\n");
+ // for (int i=0;i<pAssembler->sizeOfR();i++){
+ // printf("%g ",lsys->getFromRightHandSide(i));
+ // }
+ // printf("\n");
+
+ // solving
+ lsys->systemSolve();
+}
+
+
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index c0fedc0bb2..c3bb8026d5 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -63,7 +63,7 @@ class elasticitySolver{
void setMesh(const std::string &meshFileName);
virtual void solve();
void readInputFile(const std::string &meshFileName);
- PView *buildDisplacementView(const std::string &postFileName);
+ virtual PView *buildDisplacementView(const std::string &postFileName);
// PView *buildVonMisesView(const std::string &postFileName);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int);
@@ -71,4 +71,16 @@ class elasticitySolver{
// (const std::string &errorFileName, const elasticityData &ref, double, int);
};
+
+
+// another elastic solver ...
+class MyelasticitySolver : public elasticitySolver
+{
+ public:
+ MyelasticitySolver(int tag) : elasticitySolver(tag) {}
+ virtual void solve();
+};
+
+
+
#endif
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 3536058089..ee759aa75c 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -179,4 +179,42 @@ class femTerm {
}
};
+
+
+class DummyfemTerm : public femTerm<double>
+{
+ public:
+ typedef dofTraits<double>::VecType dataVec;
+ typedef dofTraits<double>::MatType dataMat;
+ DummyfemTerm(GModel *gm) : femTerm<double>(gm) {}
+ virtual ~DummyfemTerm() {}
+ private : // i dont want to mess with this anymore
+ virtual int sizeOfC(SElement *se) const {return 0;}
+ virtual int sizeOfR(SElement *se) const {return 0;}
+ virtual Dof getLocalDofR(SElement *se, int iRow) const {return Dof(0,0);}
+ virtual Dof getLocalDofC(SElement *se, int iCol) const {return Dof(0,0);}
+ virtual void elementMatrix(SElement *se, fullMatrix<dataMat> &m) const {m.scale(0.);}
+ virtual void elementVector(SElement *se, fullVector<dataVec> &m) const {m.scale(0.);}
+ public:
+ void addToMatrix(dofManager<dataVec> &dm,
+ fullMatrix<dataMat> &localMatrix,
+ std::vector<Dof> R,std::vector<Dof> C) const
+ {
+ dm.assemble(R, C, localMatrix);
+ }
+
+
+ void addToMatrix(dofManager<dataVec> &dm,
+ fullMatrix<dataMat> &localMatrix,
+ std::vector<Dof> R) const // symmetric version.
+ {
+ dm.assemble(R, localMatrix);
+ }
+
+};
+
+
+
+
+
#endif
diff --git a/Solver/functionSpace.cpp b/Solver/functionSpace.cpp
new file mode 100644
index 0000000000..853c23188f
--- /dev/null
+++ b/Solver/functionSpace.cpp
@@ -0,0 +1,14 @@
+//
+// C++ Implementation: functionSpace
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include "functionSpace.h"
+
+const SVector3 VectorLagrangeFunctionSpace::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index cc7ea643c9..3a540f69ec 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -2,12 +2,14 @@
#define _FUNCTION_SPACE_H_
#include "SVector3.h"
+#include "STensor3.h"
#include <vector>
#include <iterator>
#include "Numeric.h"
+#include "MElement.h"
+#include "dofManager.h"
-
-class STensor3{};
+//class STensor3{};
class SVoid{};
class basisFunction{
@@ -71,6 +73,7 @@ class FunctionSpace {
public:
virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0;
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)=0;
+// virtual groupOfElements* getSupport()=0;// probablement inutile
// virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads, STensor3 &invjac)=0;// on passe le jacobien que l'on veut ...
/* virtual int hessf(MElement *ele, double u, double v, double w,std::vector<HessType> &hesss);
virtual int divf(MElement *ele, double u, double v, double w,std::vector<DivType> &divs);
@@ -112,7 +115,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
int ndofs= ele->getNumVertices();
- grads.reserve(grads.size()+ndofs);
+// grads.reserve(grads.size()+ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
double jac[3][3];
@@ -136,12 +139,13 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
virtual int getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
{
int ndofs= ele->getNumVertices();
- keys.reserve(keys.size()+ndofs);
+// keys.reserve(keys.size()+ndofs);
for (int i=0;i<ndofs;++i)
keys.push_back(getLocalDof(ele,i));
}
};
+/*
template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T> // scalarFS * const vector (avec vecteur non const, peut etre utilise pour xfem directement
{
public :
@@ -163,41 +167,160 @@ public :
int nbdofs=valsd.size();
for (int i=0;i<nbdofs;++i) vals.push_back(multiplier*valsd[i]);
}
- virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads){};
+
+ virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+ {
+ std::vector<SVector3> gradsd;
+ ScalarFS->gradf(ele,u,v,w,gradsd);
+ int nbdofs=gradsd.size();
+ GradType val;
+ for (int i=0;i<nbdofs;++i)
+ {
+ tensprod(multiplier,gradsd[i],val);
+ grads.push_back(val);
+ }
+ };
virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele);}
virtual int getKeys(MElement *ele, Dof *keys){ ScalarFS->getKeys(ele,keys);}
virtual int getKeys(MElement *ele, std::vector<Dof> &keys) {ScalarFS->getKeys(ele,keys);}
};
+*/
+template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T>
+{
+public :
+ typedef typename TensorialTraits<T>::ValType ValType;
+ typedef typename TensorialTraits<T>::GradType GradType;
+ typedef typename TensorialTraits<T>::HessType HessType;
+ typedef typename TensorialTraits<T>::DivType DivType;
+ typedef typename TensorialTraits<T>::CurlType CurlType;
+protected :
+ std::vector<T> multipliers;
+ std::vector<int> comp;
+ FunctionSpace<double> *ScalarFS;
+public :
+ template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult, int comp_): ScalarFS(new T2(SFS))
+ {
+ multipliers.push_back(mult);comp.push_back(comp_);
+ }
+
+ template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult1, int comp1_,
+ const T& mult2, int comp2_): ScalarFS(new T2(SFS))
+ {
+ multipliers.push_back(mult1);multipliers.push_back(mult2);comp.push_back(comp1_);comp.push_back(comp2_);
+ }
+
+ template <class T2> ScalarToAnyFunctionSpace(const T2 &SFS, const T& mult1, int comp1_,
+ const T& mult2, int comp2_,const T& mult3, int comp3_): ScalarFS(new T2(SFS))
+ {
+ multipliers.push_back(mult1);multipliers.push_back(mult2);multipliers.push_back(mult3);
+ comp.push_back(comp1_);comp.push_back(comp2_);comp.push_back(comp3_);
+ }
-class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3> // it is a scalar lagrange times a constant vector.
+ virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;}
+ virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+ {
+ std::vector<double> valsd;
+ ScalarFS->f(ele,u,v,w,valsd);
+ int nbdofs=valsd.size();
+ int nbcomp=comp.size();
+ for (int j=0;j<nbcomp;++j)
+ {
+ for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]);
+ }
+ }
+
+ virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+ {
+ std::vector<SVector3> gradsd;
+ ScalarFS->gradf(ele,u,v,w,gradsd);
+ int nbdofs=gradsd.size();
+ int nbcomp=comp.size();
+ GradType val;
+ for (int j=0;j<nbcomp;++j)
+ {
+ for (int i=0;i<nbdofs;++i)
+ {
+ tensprod(multipliers[j],gradsd[i],val);
+ grads.push_back(val);
+ }
+ }
+ };
+ virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele)*comp.size();}
+ virtual int getKeys(MElement *ele, Dof *keys)
+ {
+ int nbcomp=comp.size();
+ int nk=ScalarFS->getNumKeys(ele);
+ Dof *kptr=keys;
+ ScalarFS->getKeys(ele,kptr);
+ kptr+=nk;
+ for (int j=1;j<nbcomp;++j)
+ {
+ for (int i=0;i<nk;++i)
+ kptr[i]=kptr[i-nk];
+ kptr+=nk;
+ }
+ kptr=keys;
+ for (int j=0;j<nbcomp;++j)
+ {
+ for (int i=0;i<nk;++i)
+ {
+ int i1,i2;
+ Dof::getTwoIntsFromType(kptr[i].getType(), i1,i2);
+ kptr[i]=Dof(kptr[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i2));
+ }
+ kptr+=nk;
+ }
+ }
+ virtual int getKeys(MElement *ele, std::vector<Dof> &keys)
+ {
+ int nk=ScalarFS->getNumKeys(ele);
+ std::vector<Dof> bufk;
+ bufk.reserve(nk);
+ ScalarFS->getKeys(ele,bufk);
+ int nbcomp=comp.size();
+ for (int j=0;j<nbcomp;++j)
+ {
+ for (int i=0;i<nk;++i)
+ {
+ int i1,i2;
+ Dof::getTwoIntsFromType(bufk[i].getType(), i1,i2);
+ keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i2)));
+ }
+ }
+ }
+};
+
+class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
{
public:
enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 };
+ static const SVector3 BasisVectors[3];
+
typedef TensorialTraits<SVector3>::ValType ValType;
typedef TensorialTraits<SVector3>::GradType GradType;
typedef TensorialTraits<SVector3>::HessType HessType;
typedef TensorialTraits<SVector3>::DivType DivType;
typedef TensorialTraits<SVector3>::CurlType CurlType;
- protected:
-// Along direction;
- public:
- VectorLagrangeFunctionSpace(int id,Along t) : ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),SVector3(0.,0.,0.) )//,direction(t)
- {
- multiplier[t]=1.;
- }
-// virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals);
-// virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
-// virtual int getNumKeys(MElement *ele);
-// virtual int getKeys(MElement *ele, Dof *keys);
-// virtual int getKeys(MElement *ele, std::vector<Dof> &keys);
+ VectorLagrangeFunctionSpace(int id) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
+ SVector3(1.,0.,0.),VECTOR_X, SVector3(0.,1.,0.),VECTOR_Y,SVector3(0.,0.,1.),VECTOR_Z)
+ {}
+ VectorLagrangeFunctionSpace(int id,Along comp1) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
+ BasisVectors[comp1],comp1)
+ {}
+ VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
+ BasisVectors[comp1],comp1, BasisVectors[comp2],comp2)
+ {}
+ VectorLagrangeFunctionSpace(int id,Along comp1,Along comp2, Along comp3) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpace(id),
+ BasisVectors[comp1],comp3, BasisVectors[comp2],comp2, BasisVectors[comp3],comp3)
+ {}
};
-
-
-
-
template<class T>
class CompositeFunctionSpace : public FunctionSpace<T>
{
@@ -235,8 +358,18 @@ class CompositeFunctionSpace : public FunctionSpace<T>
delete (*it);
}
- virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals) {}
- virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) {}
+ virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
+ {
+ for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+ (*it)->f(ele,u,v,w,vals);
+ }
+
+ virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+ {
+ for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+ (*it)->gradf(ele,u,v,w,grads);
+ }
+
virtual int getNumKeys(MElement *ele)
{
int ndofs=0;
@@ -244,6 +377,7 @@ class CompositeFunctionSpace : public FunctionSpace<T>
ndofs+=(*it)->getNumKeys(ele);
return ndofs;
}
+
virtual int getKeys(MElement *ele, Dof *keys)
{
Dof *kptr=keys;
@@ -280,4 +414,28 @@ class xFemFunctionSpace : public FunctionSpace<T>
};
+template<class T,class F>
+class FilteredFunctionSpace : public FunctionSpace<T>
+{
+ typedef typename TensorialTraits<T>::ValType ValType;
+ typedef typename TensorialTraits<T>::GradType GradType;
+ typedef typename TensorialTraits<T>::HessType HessType;
+ typedef typename TensorialTraits<T>::DivType DivType;
+ typedef typename TensorialTraits<T>::CurlType CurlType;
+ private:
+ FunctionSpace<T>* _spacebase;
+ F &_filter;
+
+ public:
+ ValType f(MElement *ele, double u, double v, double w);
+ GradType gradf(MElement *ele, double u, double v, double w);
+ int getNumKeys(MElement *ele);
+ void getKeys(MElement *ele, Dof *keys);
+ void getKeys(MElement *ele, std::vector<Dof> &keys);
+};
+
+
+
+
+
#endif
--
GitLab