Skip to content
Snippets Groups Projects
Commit a19c7ee3 authored by Éric Béchet's avatar Éric Béchet
Browse files

No commit message

No commit message
parent 42fad48e
No related branches found
No related tags found
No related merge requests found
...@@ -15,10 +15,12 @@ ...@@ -15,10 +15,12 @@
#include <vector> #include <vector>
#include <complex> #include <complex>
#include <iostream> #include <iostream>
//#include <tr1/memory>
#include "dofManager.h" #include "dofManager.h"
#include "SVector3.h" #include "SVector3.h"
#include "MElement.h" #include "MElement.h"
#include "MVertex.h"
...@@ -31,7 +33,7 @@ template<class T> struct TensorialTraits ...@@ -31,7 +33,7 @@ template<class T> struct TensorialTraits
typedef T ValType; typedef T ValType;
typedef T GradType[3]; typedef T GradType[3];
typedef T HessType[3][3]; typedef T HessType[3][3];
static const int nb_basis_vectors=1; // par défaut considéré comme un scalaire static const int nb_basis_vectors=1; // par défaut, considéré comme un scalaire
}; };
template<> struct TensorialTraits<double> template<> struct TensorialTraits<double>
...@@ -56,13 +58,31 @@ template<class T> class Function // Fonction au sens EF du terme. ...@@ -56,13 +58,31 @@ template<class T> class Function // Fonction au sens EF du terme.
// renvoie valeur , gradient, hessien, etc...pour un element donné et un point donné // renvoie valeur , gradient, hessien, etc...pour un element donné et un point donné
{ {
public: public:
// typedef std::tr1::shared_ptr<Function<T> > FunctionPtr;
virtual ~Function(){} virtual ~Function(){}
virtual void GetVal (double uvw[],MElement *e, T& Val)const{};// =0; virtual void GetVal (double uvw[],MElement *e, T& Val)const=0;
virtual void GetGrad(double uvw[],MElement *e,typename TensorialTraits<T>::GradType &Grad) const {};//=0; virtual void GetGrad(double uvw[],MElement *e,typename TensorialTraits<T>::GradType &Grad) const =0;
virtual void GetHess(double uvw[],MElement *e,typename TensorialTraits<T>::HessType &Hess)const {};//=0; virtual void GetHess(double uvw[],MElement *e,typename TensorialTraits<T>::HessType &Hess)const =0;
}; };
class SpaceBase // renvoie des clefs de dofs
class LagrangeShapeFunction: public Function<double>
{
public:
virtual void GetVal(double uvw[],MElement *e, double& Val)
{
// double s[100];
// _e->getShapeFunctions(uvw[0], uvw[1], uvw[2], s);
}
};
class SpaceBase // renvoie des clefs des dofs
{ {
public: public:
virtual void getDofs(MElement *e,std::vector<Dof> &vecD)=0; virtual void getDofs(MElement *e,std::vector<Dof> &vecD)=0;
...@@ -73,8 +93,8 @@ template<class T> class Space : public SpaceBase // renvoie des clefs de dofs et ...@@ -73,8 +93,8 @@ template<class T> class Space : public SpaceBase // renvoie des clefs de dofs et
public: public:
Space(){}; Space(){};
virtual ~Space(){}; virtual ~Space(){};
virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF) {} virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF)=0;
virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF){};//=0; virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF)=0;
}; };
...@@ -94,7 +114,10 @@ public: ...@@ -94,7 +114,10 @@ public:
SpaceLagrange(int iField):_iField(iField){}; SpaceLagrange(int iField):_iField(iField){};
virtual ~SpaceLagrange(){}; virtual ~SpaceLagrange(){};
virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF){} virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF){}
virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF){};//=0; virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF)
{
};
virtual void getDofs(MElement *e,std::vector<Dof> &vecD) virtual void getDofs(MElement *e,std::vector<Dof> &vecD)
{ {
int ndofs= e->getNumVertices()*TensorialTraits<T>::nb_basis_vectors; int ndofs= e->getNumVertices()*TensorialTraits<T>::nb_basis_vectors;
...@@ -104,49 +127,36 @@ public: ...@@ -104,49 +127,36 @@ public:
}; };
template<class T> class Field : public Function<T> // renvoie des valeurs de champ (ff*valeurs dofs), gradient , etc... template<class T> class SpaceXfem : public Space<T> // approximation xfem (spacebase + enrich)
{ {
public: int _iEnrich;
Field(){}; Space<T> &_SpaceBase;
virtual ~Field(){}; Function<double> &enrich;
}; std::unary_function<MVertex*,bool> &test;
Dof getLocalDof(MElement *e, int i) const
template<class T1,class T2> class TermBilinear // terme associe a un "element" / pt de gauss (contribution élémentaire)
// typiquement celui ci stoque ce qui doit etre stocke. C'est la base configurable du code
// on doit associer cela a un allocateur qui renvoie un pointeur sur ce truc
// Soit tous les elements/pts de gauss ont le meme terme , allocateur unique (pas de stockage aux pts de
// gauss par exemple
// soit ils ont qqc a stockuer et sont tous distincts
// on peut utiliser des membres statiques pour ce qui est constant pour tous les instances
{ {
int iComp = i / e->getNumVertices();
int ithLocalVertex = i % e->getNumVertices();
MVertex * v= e->getVertex(ithLocalVertex);
if (test(v))
return Dof(e->getVertex(ithLocalVertex)->getNum(),
Dof::createTypeWithTwoInts(iComp, _iEnrich));
}
public: public:
typedef void Stype; SpaceXfem(int iEnrich,Space<T>& SpaceBase):_iEnrich(iEnrich), _SpaceBase(SpaceBase){};
// TermBilinear(){}; virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF){}
virtual double getTerm(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF,Stype* info) {}; virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF){};
virtual void Update(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {}; virtual void getDofs(MElement *e,std::vector<Dof> &vecD)
// toute fonctios utiles . prevoir un algorithme ad hoc pour l'appliquer sur ts les pts de gauss ...
};
class TermBilinearMeca : public TermBilinear<double,double>
{ {
public: _SpaceBase.getDofs(e,vecD);
// typedef Tensor Stype; }
// TermBilinear(){};
// virtual double getTerm(double uvw[],MElement &e,Function<double> &SF,Function<double> &TF,Stype* info) {};
// virtual void Update(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
// toute fonctios utiles . prevoir un algorithme
}; };
template<class T> class Field : public Function<T> // renvoie des valeurs de champ (ff*valeurs dofs), gradient , etc...
class TermBilinearMecaNL : public TermBilinearMeca
{ {
public: public:
typedef Tensor2 Stype; Field(){};
// TermBilinear(){}; virtual ~Field(){};
// virtual double getTerm(double uvw[],MElement &e,Function<double> &SF,Function<double> &TF,Stype* info) {};
// virtual void Update(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
// toute fonctios utiles . prevoir un algorithme
}; };
...@@ -156,12 +166,15 @@ class FormBilinearBase ...@@ -156,12 +166,15 @@ class FormBilinearBase
public: public:
template <class T> void allocate(T *p) {p=new T[10];std::cout << "10" << std::endl;} template <class T> void allocate(T *p) {p=new T[10];std::cout << "10" << std::endl;}
template <class T> void get(T* tab, T *p, int ndx) {p=&tab[ndx];} template <class T> void get(T* tab, T *p, int ndx) {p=&tab[ndx];}
virtual void getDofs(MElement *e)=0;
virtual void getFuncs(MElement *e,const double uvw[3])=0;
virtual void getGradFuncs(MElement *e,const double uvw[3])=0;
}; };
template <> void FormBilinearBase::allocate(void *p) {p=NULL; std::cout << "null" << std::endl;} template <> void FormBilinearBase::allocate(void *p) {p=NULL; std::cout << "null" << std::endl;}
template <> void FormBilinearBase::get(void *tab,void *p, int ndx) {p=NULL;} template <> void FormBilinearBase::get(void *tab,void *p, int ndx) {p=NULL;}
template <class Term> class FormBilinear : public FormBilinearBase template <class Term,class SpaceL,class SpaceR> class FormBilinear : public FormBilinearBase
// Renvoie des matrices élémentaires (ff) // Renvoie des matrices élémentaires (ff)
// Accessoirement stocke des pointeurs vers les termes pour chaque element // Accessoirement stocke des pointeurs vers les termes pour chaque element
// Doit etre initialisée AVANT toute opération (pour l'allocation) // Doit etre initialisée AVANT toute opération (pour l'allocation)
...@@ -171,17 +184,49 @@ template <class Term> class FormBilinear : public FormBilinearBase ...@@ -171,17 +184,49 @@ template <class Term> class FormBilinear : public FormBilinearBase
{ {
typedef typename Term::Stype Stype; // le truc a stocker typedef typename Term::Stype Stype; // le truc a stocker
typedef typename Term::Datamat Datamat; // type elementaire
Term instance; // une instance du terme. E.g. utile pour toute initialisation (calcul tenseur materiel, etc...) Term instance; // une instance du terme. E.g. utile pour toute initialisation (calcul tenseur materiel, etc...)
Stype *p; Stype *p;
MElement *e; MElement *e;
Function<double> *a; Function<double> *a;
Function<double> *b; Function<double> *b;
double uvw[3]; double uvw[3];
fullMatrix<Datamat> mat;
SpaceL *_spacel;
SpaceR *_spacer;
std::vector<Dof> Dofsl;
std::vector<Dof> Dofsr;
std::vector<Datamat> Funcl;
std::vector<Datamat> Funcr;
public: public:
FormBilinear(SpaceL &sl,SpaceR &sr) {}
void getMatrix(fullMatrix<Datamat> *v) {v=&mat;}
void func(void) { allocate(p);} void func(void) { allocate(p);}
void func2(void) { instance.getTerm(uvw,*e,*a,*b,p); } void func2(void) { }
void getDofs(MElement *e)
{
_spacel->getDofs(e,Dofsl);
_spacer->getDofs(e,Dofsr);
};
void getFuncs(MElement *e,const double uvw[3]) {};
void getGradFuncs(MElement *e,const double uvw[3]) {};
void Init(MElement *e) { getDofs(e);} // called once for each element
void Accumulate(MElement *e,const double uvw[3]) // called for every GP
{
// for (int i)
// for (int j)
// mat(i,j)+=instance.getTerm(uvw,*e,*a,*b,NULL,this);
instance.Init(e,this); // called once for each GP
if (instance.NeedFuncs()) getFuncs(e,uvw);
if (instance.NeedGradFuncs()) getGradFuncs(e,uvw);
Datamat result;
instance.getTerm(uvw,e,NULL,*a,*b,result); // called for every SF combination (times every GP, times every element)
}
}; };
//template <class Term,class SpaceL,class datamat=double> class FormBilinear<Term,SpaceL,SpaceR> : public FormBilinearBase {}; cas symétrique
class FormLinear{}; // renvoie des vecteurs élémentaires class FormLinear{}; // renvoie des vecteurs élémentaires
// on devrait pouvoir construire une forme lin à partir d'une forme bilin pour les pb "matrix free" // on devrait pouvoir construire une forme lin à partir d'une forme bilin pour les pb "matrix free"
...@@ -198,5 +243,54 @@ void Construct(FormZero &Z,Region &R,Integrator &I); ...@@ -198,5 +243,54 @@ void Construct(FormZero &Z,Region &R,Integrator &I);
void Assemble(FormZero &Z,Region &R,Assembler &A,Integrator &I); void Assemble(FormZero &Z,Region &R,Assembler &A,Integrator &I);
*/ */
template<class datamat, class T1,class T2> class TermBilinear // terme associe a un "element" / pt de gauss (contribution élémentaire)
// typiquement celui ci stoque ce qui doit etre stocke. C'est la base configurable du code
// on doit associer cela a un allocateur qui renvoie un pointeur sur ce truc
// Soit tous les elements/pts de gauss ont le meme terme , allocateur unique (pas de stockage aux pts de
// gauss par exemple
// soit ils ont qqc a stockuer et sont tous distincts
// on peut utiliser des membres statiques pour ce qui est constant pour tous les instances
{
public:
typedef void Stype;
typedef datamat Datamat;
// TermBilinear(){};
virtual void Init(MElement *e, FormBilinearBase * caller) //called once for every element
{
}
virtual bool NeedGradFuncs(void) { return false ;}// Query if getTerm will use the gradient of the SF
virtual bool NeedFuncs(void) { return false ;}// Query if it will need the value of the SF
virtual void getTerm(const double uvw[],MElement *e,Stype* info, Function<T1> &SF,Function<T2> &TF, datamat &res ) // called for every gauss point
{
};
virtual void Update(const double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
// toute fonctios utiles . prevoir un algorithme ad hoc pour l'appliquer sur ts les pts de gauss ...
};
class TermBilinearMeca : public TermBilinear<double, double,double>
{
public:
virtual bool NeedGradFuncs(void) { return true ;}// Query if getTerm will use the gradient of the SF
// typedef Tensor Stype;
// TermBilinear(){};
// virtual double getTerm(double uvw[],MElement &e,Function<double> &SF,Function<double> &TF,Stype* info) {};
// virtual void Update(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
// toute fonctios utiles . prevoir un algorithme
};
class TermBilinearMecaNL : public TermBilinearMeca
{
public:
typedef Tensor2 Stype;
// TermBilinear(){};
// virtual double getTerm(double uvw[],MElement &e,Function<double> &SF,Function<double> &TF,Stype* info) {};
// virtual void Update(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
// toute fonctios utiles . prevoir un algorithme
};
#endif // highlevel_H #endif // highlevel_H
...@@ -11,17 +11,16 @@ int main (int argc, char* argv[]) ...@@ -11,17 +11,16 @@ int main (int argc, char* argv[])
return -1; return -1;
} }
MElement *e;
FormBilinear<TermBilinearMeca> f; const double uvw[3]={0.,0.,0.};
SpaceLagrange<double> L(123);
FormBilinear<TermBilinearMeca,SpaceLagrange<double>,SpaceLagrange<double> > f(L,L);
f.func(); f.func();
f.func2(); f.Accumulate(e,uvw);
FormBilinear<TermBilinearMecaNL> fnl; FormBilinear<TermBilinearMecaNL,SpaceLagrange<double>,SpaceLagrange<double> > fnl(L,L);
fnl.func(); fnl.func();
fnl.func2();
SpaceLagrange<double> L(123);
// globals are still present in Gmsh // globals are still present in Gmsh
GmshInitialize(argc, argv); GmshInitialize(argc, argv);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment