diff --git a/contrib/arc/highlevel.h b/contrib/arc/highlevel.h
index 8ed6d0ab493de289190a13368d91a1b63af9c6a0..70a497433d43967699f0320e6dfec8fa9ade87a2 100644
--- a/contrib/arc/highlevel.h
+++ b/contrib/arc/highlevel.h
@@ -15,10 +15,12 @@
#include <vector>
#include <complex>
#include <iostream>
+//#include <tr1/memory>
#include "dofManager.h"
#include "SVector3.h"
#include "MElement.h"
+#include "MVertex.h"
@@ -31,7 +33,7 @@ template<class T> struct TensorialTraits
typedef T ValType;
typedef T GradType[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>
@@ -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é
{
public:
- virtual ~Function(){}
- 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 GetHess(double uvw[],MElement *e,typename TensorialTraits<T>::HessType &Hess)const {};//=0;
+// typedef std::tr1::shared_ptr<Function<T> > FunctionPtr;
+ virtual ~Function(){}
+ 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 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:
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
public:
Space(){};
virtual ~Space(){};
- 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 getDofsandSFs(MElement *e,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:
SpaceLagrange(int iField):_iField(iField){};
virtual ~SpaceLagrange(){};
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)
{
int ndofs= e->getNumVertices()*TensorialTraits<T>::nb_basis_vectors;
@@ -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)
{
+ int _iEnrich;
+ Space<T> &_SpaceBase;
+ Function<double> &enrich;
+ std::unary_function<MVertex*,bool> &test;
+ Dof getLocalDof(MElement *e, int i) const
+ {
+ 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:
- Field(){};
- virtual ~Field(){};
-};
-
-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
-{
-public:
- typedef void Stype;
-// TermBilinear(){};
- virtual double getTerm(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF,Stype* info) {};
- virtual void Update(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>
-{
-public:
-// 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
+ SpaceXfem(int iEnrich,Space<T>& SpaceBase):_iEnrich(iEnrich), _SpaceBase(SpaceBase){};
+ virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF){}
+ virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF){};
+ virtual void getDofs(MElement *e,std::vector<Dof> &vecD)
+ {
+ _SpaceBase.getDofs(e,vecD);
+ }
};
-
-class TermBilinearMecaNL : public TermBilinearMeca
+template<class T> class Field : public Function<T> // renvoie des valeurs de champ (ff*valeurs dofs), gradient , etc...
{
-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
+public:
+ Field(){};
+ virtual ~Field(){};
};
@@ -156,12 +166,15 @@ class FormBilinearBase
public:
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];}
+ 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::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)
// Accessoirement stocke des pointeurs vers les termes pour chaque element
// Doit etre initialisée AVANT toute opération (pour l'allocation)
@@ -171,17 +184,49 @@ template <class Term> class FormBilinear : public FormBilinearBase
{
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...)
Stype *p;
MElement *e;
Function<double> *a;
Function<double> *b;
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:
+ FormBilinear(SpaceL &sl,SpaceR &sr) {}
+ void getMatrix(fullMatrix<Datamat> *v) {v=&mat;}
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
// 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);
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
diff --git a/contrib/arc/mainElasticity.cpp b/contrib/arc/mainElasticity.cpp
index e39b9fe4118b9e1d338f93d22a1881b111f6991d..99b08ca69a7c04d80e117d95e7f74148c2b2892d 100644
--- a/contrib/arc/mainElasticity.cpp
+++ b/contrib/arc/mainElasticity.cpp
@@ -11,18 +11,17 @@ int main (int argc, char* argv[])
return -1;
}
-
- FormBilinear<TermBilinearMeca> f;
+ MElement *e;
+ const double uvw[3]={0.,0.,0.};
+ SpaceLagrange<double> L(123);
+ FormBilinear<TermBilinearMeca,SpaceLagrange<double>,SpaceLagrange<double> > f(L,L);
f.func();
- f.func2();
+ f.Accumulate(e,uvw);
- FormBilinear<TermBilinearMecaNL> fnl;
+ FormBilinear<TermBilinearMecaNL,SpaceLagrange<double>,SpaceLagrange<double> > fnl(L,L);
fnl.func();
- fnl.func2();
-
- SpaceLagrange<double> L(123);
-
+
// globals are still present in Gmsh
GmshInitialize(argc, argv);