From 457909676ebb95687e7e78d7eb0ab7eb29b937f8 Mon Sep 17 00:00:00 2001 From: Eric Bechet <eric.bechet@ulg.ac.be> Date: Wed, 13 Jan 2010 15:47:44 +0000 Subject: [PATCH] no more mVertex* - related members in function spaces. Use Mpoint instead. --- Solver/elasticitySolver.cpp | 19 ++-- Solver/function.h | 7 ++ Solver/functionSpace.h | 160 ++++----------------------------- Solver/solverAlgorithms.h | 29 ------ Solver/solverField.h | 25 +----- Solver/terms.h | 15 ---- contrib/arc/conge.dat | 5 +- contrib/arc/mainElasticity.cpp | 104 +++++++++++++++++++-- 8 files changed, 133 insertions(+), 231 deletions(-) diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp index 35051215ff..e461f4522e 100644 --- a/Solver/elasticitySolver.cpp +++ b/Solver/elasticitySolver.cpp @@ -16,7 +16,7 @@ #include "solverAlgorithms.h" #include "quadratureRules.h" #include "solverField.h" - +#include "MPoint.h" #if defined(HAVE_POST) #include "PView.h" #include "PViewData.h" @@ -164,14 +164,11 @@ void elasticitySolver::solve() // give priority to fixations : when a dof is fixed, it cannot be // numbered afterwards - + std::cout << "Dirichlet BC"<< std::endl; for (unsigned int i = 0; i < allDirichlet.size(); i++) { FilterDofComponent filter(allDirichlet[i]._comp); - if (allDirichlet[i].onWhat==BoundaryCondition::ON_VERTEX) - FixNodalDofs(*LagSpace,allDirichlet[i].g->vbegin(),allDirichlet[i].g->vend(),*pAssembler,allDirichlet[i]._f,filter); - else - FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter); + FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter); } // we number the dofs : when a dof is numbered, it cannot be numbered @@ -185,14 +182,11 @@ void elasticitySolver::solve() // First build the force vector GaussQuadrature Integ_Boundary(GaussQuadrature::Val); - + std::cout << "Neumann BC"<< std::endl; for (unsigned int i = 0; i < allNeumann.size(); i++) { LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f); - if (allNeumann[i].onWhat==BoundaryCondition::ON_VERTEX) - Assemble(Lterm,*LagSpace,allNeumann[i].g->vbegin(),allNeumann[i].g->vend(),*pAssembler); - else - Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler); + Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler); } // bulk material law @@ -279,7 +273,8 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName) for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it) { SVector3 val; - Field.f(*it,val); + MPoint p(*it); + Field.f(&p,0,0,0,val); std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2); data[(*it)->getNum()]=vec; } diff --git a/Solver/function.h b/Solver/function.h index adeabe4873..51916a46d1 100644 --- a/Solver/function.h +++ b/Solver/function.h @@ -51,6 +51,13 @@ public : inline bool somethingDependOnMe() { return !_dependOnMe.empty(); } + inline int howManyDependOnMe() { + return _dependOnMe.size(); + } + inline int howManyDoIDependOn() { + return _iDependOn.size(); + } + }; // dataCache when the value is a double diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h index bfe5df1273..1b7880f3f4 100644 --- a/Solver/functionSpace.h +++ b/Solver/functionSpace.h @@ -10,66 +10,40 @@ #include "MElement.h" #include "dofManager.h" -class SVoid{}; - -class basisFunction{ - public: - virtual void f(MElement *ele, double u, double v, double w, double *s) = 0; - virtual void df(MElement *ele, double u, double v, double w, double *s) = 0; -}; - -class scalarPolynomialBasisFunction : public basisFunction{ - private: - polynomialBasis *_basis; - public: - virtual void f(MElement *ele, double u, double v, double w, double *s); - virtual void df(MElement *ele, double u, double v, double w, double *s); -}; - -class vectorPolynomialBasisFunction : public basisFunction{ - private: - polynomialBasis *_basis[3]; - public: - virtual void f(MElement *ele, double u, double v, double w, double *s); - virtual void df(MElement *ele, double u, double v, double w, double *s); -}; +//class SVoid{}; template<class T> struct TensorialTraits { typedef T ValType; typedef T GradType[3]; - typedef T HessType[3][3]; +/* typedef T HessType[3][3]; typedef SVoid DivType; - typedef SVoid CurlType; + typedef SVoid CurlType;*/ }; template<> struct TensorialTraits<double> { typedef double ValType; typedef SVector3 GradType; - typedef STensor3 HessType; +/* typedef STensor3 HessType; typedef SVoid DivType; - typedef SVoid CurlType; + typedef SVoid CurlType;*/ }; template<> struct TensorialTraits<SVector3> { typedef SVector3 ValType; typedef STensor3 GradType; - typedef SVoid HessType; +/* typedef SVoid HessType; typedef double DivType; - typedef SVector3 CurlType; + typedef SVector3 CurlType;*/ }; - - class FunctionSpaceBase { public: virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0; - virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs - virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)=0; }; template<class T> @@ -78,23 +52,10 @@ class FunctionSpace : public FunctionSpaceBase 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;*/ - - public: - virtual void f(MVertex *ver, std::vector<ValType> &vals) {} // used for neumann BC virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0; virtual void 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); - virtual int curlf(MElement *ele, double u, double v, double w,std::vector<CurlType> &curls);*/ virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0; - virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs - virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)=0; }; class ScalarLagrangeFunctionSpace : public FunctionSpace<double> @@ -102,14 +63,16 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> public: typedef TensorialTraits<double>::ValType ValType; typedef TensorialTraits<double>::GradType GradType; - typedef TensorialTraits<double>::HessType HessType; - typedef TensorialTraits<double>::DivType DivType; - typedef TensorialTraits<double>::CurlType CurlType; - protected: - std::vector<basisFunction*> basisFunctions; + protected: int _iField; // field number (used to build dof keys) + private: + virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) + { + keys.push_back(Dof(ver->getNum(), _iField)); + } + public: ScalarLagrangeFunctionSpace(int i=0):_iField(i) {} virtual int getId(void) const {return _iField;}; @@ -119,8 +82,8 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> int curpos=vals.size(); vals.resize(curpos+ndofs); ele->getShapeFunctions(u, v, w, &(vals[curpos])); - }; - virtual void f(MVertex *ver, std::vector<ValType> &vals) {vals.push_back(1.0);} // used for neumann BC + } + virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) { int ndofs= ele->getNumVertices(); @@ -137,7 +100,8 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2], invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2] )); - }; + } + virtual int getNumKeys(MElement *ele) {return ele->getNumVertices();} virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ... @@ -147,11 +111,6 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> for (int i=0;i<ndofs;++i) getKeys(ele->getVertex(i),keys); } - virtual int getNumKeys(MVertex *ver) { return 1;} - virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) - { - keys.push_back(Dof(ver->getNum(), _iField)); - }; }; @@ -160,9 +119,6 @@ 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; @@ -187,19 +143,6 @@ public : } virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;} - virtual void f(MVertex *ver, std::vector<ValType> &vals) - { - std::vector<double> valsd; - ScalarFS->f(ver,valsd); - int nbdofs=valsd.size(); - int nbcomp=comp.size(); - int curpos=vals.size(); - vals.reserve(curpos+nbcomp*nbdofs); - for (int j=0;j<nbcomp;++j) - { - for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]); - } - } // used for neumann BC virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) { @@ -257,28 +200,6 @@ public : } } } - - virtual int getNumKeys(MVertex *ver) { return ScalarFS->getNumKeys(ver)*comp.size();} - virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) - { - int nk=ScalarFS->getNumKeys(ver); - std::vector<Dof> bufk; - bufk.reserve(nk); - ScalarFS->getKeys(ver,bufk); - int nbdofs=bufk.size(); - int nbcomp=comp.size(); - int curpos=keys.size(); - keys.reserve(curpos+nbcomp*nbdofs); - 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],i1))); - } - } - }; }; class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3> @@ -289,9 +210,6 @@ class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3> enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 }; 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; 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) @@ -317,9 +235,6 @@ class CompositeFunctionSpace : 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; typedef typename std::vector<FunctionSpace<T>* >::iterator iterFS; protected: @@ -347,13 +262,6 @@ class CompositeFunctionSpace : public FunctionSpace<T> for (iterFS it=_spaces.begin(); it!=_spaces.end();++it) delete (*it); } - - virtual void f(MVertex *ver, std::vector<ValType> &vals) - { - for (iterFS it=_spaces.begin(); it!=_spaces.end();++it) - (*it)->f(ver,vals); - } - virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals) { for (iterFS it=_spaces.begin(); it!=_spaces.end();++it) @@ -379,20 +287,6 @@ class CompositeFunctionSpace : public FunctionSpace<T> for (iterFS it=_spaces.begin(); it!=_spaces.end();++it) (*it)->getKeys(ele,keys); } - - virtual int getNumKeys(MVertex *ver) - { - int ndofs=0; - for (iterFS it=_spaces.begin(); it!=_spaces.end();++it) - ndofs+=(*it)->getNumKeys(ver); - return ndofs; - } - - virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) - { - for (iterFS it=_spaces.begin(); it!=_spaces.end();++it) - (*it)->getKeys(ver,keys); - }; }; template<class T> @@ -401,21 +295,15 @@ class xFemFunctionSpace : 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; private: FunctionSpace<T>* _spacebase; // Function<double>* enrichment; public: - virtual void f(MVertex *ver, std::vector<ValType> &vals); + virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals); virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads); virtual int getNumKeys(MElement *ele); virtual void getKeys(MElement *ele, std::vector<Dof> &keys); - virtual int getNumKeys(MVertex *ver); - virtual void getKeys(MVertex *ver, std::vector<Dof> &keys); - }; @@ -424,25 +312,15 @@ 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: - virtual void f(MVertex *ver, std::vector<ValType> &vals); virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals); virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads); virtual int getNumKeys(MElement *ele); virtual void getKeys(MElement *ele, std::vector<Dof> &keys); - virtual int getNumKeys(MVertex *ver); - virtual void getKeys(MVertex *ver, std::vector<Dof> &keys); }; - - - - #endif diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h index 1e0cdc5590..d605063fad 100644 --- a/Solver/solverAlgorithms.h +++ b/Solver/solverAlgorithms.h @@ -65,20 +65,6 @@ template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,Func } } -template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler) -{ - fullVector<typename Assembler::dataMat> localVector; - std::vector<Dof> R; - for (Iterator it = itbegin;it!=itend; ++it) - { - MVertex *v = *it; - R.clear(); - term.get(v,localVector); - space.getKeys(v,R); - assembler.assemble(R, localVector); - } -} - template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler) { fullVector<typename Assembler::dataMat> localVector; @@ -177,20 +163,6 @@ template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e } } -template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MVertex *v,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter) -{ - std::vector<Dof> R; - space.getKeys(v,R); - for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd) - { - Dof key=*itd; - if (filter(key)) - { - assembler.fixDof(key, fct(v->x(),v->y(),v->z())); - } - } -} - template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter) { for (Iterator it=itbegin;it!=itend;++it) @@ -211,5 +183,4 @@ template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &spac } } - #endif// _SOLVERALGORITHMS_H_ diff --git a/Solver/solverField.h b/Solver/solverField.h index d2d96ad570..1b8e2e9484 100644 --- a/Solver/solverField.h +++ b/Solver/solverField.h @@ -21,7 +21,7 @@ #include "functionSpace.h" template<class T> -class SolverField : public FunctionSpace<T> // being able to use it instead of a real function space is interesting (nbkeys=1, dofs undefined (or could be defined by element-wise ) +class SolverField : public FunctionSpace<T> // being able to use it instead of a real function space is interesting (nbkeys=1, explicit keys/dofs undefined (or could be defined element-wise ) { public: typedef typename TensorialTraits<T>::ValType ValType; @@ -38,26 +38,6 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) {Msg::Error("getKeys for SolverField should'nt be called");} public: - virtual void f(MVertex *ver, std::vector<ValType> &vals) - { - ValType val; - f(ver,val); - vals.push_back(val); - } - - virtual void f(MVertex *ver, ValType &val) - { - std::vector<Dof> D; - std::vector<ValType> SFVals; - std::vector<double> DMVals; - fs->getKeys(ver,D); - dm->getDofValue(D,DMVals); - fs->f(ver,SFVals); - val=ValType(); - for (int i=0;i<D.size();++i) - val+=SFVals[i]*DMVals[i]; - } - virtual void f(MElement *ele, double u, double v, double w, ValType &val) { std::vector<Dof> D; @@ -99,9 +79,6 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a } }; - - - /* class Formulation diff --git a/Solver/terms.h b/Solver/terms.h index de33311651..cce9fd5a3c 100644 --- a/Solver/terms.h +++ b/Solver/terms.h @@ -44,7 +44,6 @@ class LinearTermBase public: virtual ~LinearTermBase() {} virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0; - virtual void get(MVertex *ver,fullVector<double> &m) =0; }; template<class T1> class LinearTerm : public LinearTermBase @@ -309,20 +308,6 @@ template<class T1> class LoadTerm : public LinearTerm<T1> } } } - - virtual void get(MVertex *ver,fullVector<double> &m) - { - double nbFF=LinearTerm<T1>::space1.getNumKeys(ver); - double jac[3][3]; - m.resize(nbFF); - std::vector<typename TensorialTraits<T1>::ValType> Vals; - LinearTerm<T1>::space1.f(ver, Vals); - typename TensorialTraits<T1>::ValType load=Load(ver->x(),ver->y(),ver->z()); - for (int j = 0; j < nbFF ; ++j) - { - m(j)=dot(Vals[j],load); - } - } }; diff --git a/contrib/arc/conge.dat b/contrib/arc/conge.dat index 1b2aab2281..dc3ffaf1fa 100644 --- a/contrib/arc/conge.dat +++ b/contrib/arc/conge.dat @@ -5,10 +5,9 @@ EdgeDisplacement 8 1 0 EdgeDisplacement 8 2 0 #NodeDisplacement 20 0 0 #NodeDisplacement 20 1 0 -#NodeDisplacement 20 2 0 #NodeDisplacement 21 1 0 -#NodeDisplacement 21 2 0 +#NodeForce 21 100e6 0 0 EdgeForce 9 100e6 0 0 #FaceForce 7 0 -1.0e10 0 -#NodeForce 20 -1e7 0 0 +#NodeForce 20 -1e8 0 0 diff --git a/contrib/arc/mainElasticity.cpp b/contrib/arc/mainElasticity.cpp index d7a911468b..d654c0ca62 100644 --- a/contrib/arc/mainElasticity.cpp +++ b/contrib/arc/mainElasticity.cpp @@ -5,16 +5,106 @@ #include "highlevel.h" #include "groupOfElements.h" #include <iterator> +#include "function.h" +#include "fullMatrix.h" +/* +class functionAdd : public function +{ + private: + std::vector<std::string> strvec; + class data : public dataCacheDouble + { + private: + std::vector<dataCacheDouble *> dcvec; + public: + data(const functionAdd * fm,dataCacheMap *m) : + dataCacheDouble(m->getNbEvaluationPoints(),1) + { + for (int i=0;i<fm->strvec.size();++i) + dcvec.push_back(&(m->get(fm->strvec[i],this))); + } + void _eval() + { + _value(0, 0)=0; + for (int i=0;i<dcvec.size();++i) { _value(0, 0) += (*(dcvec[i]))(0,0);} + } + ~data() + { + } + }; + public: + dataCacheDouble *newDataCache(dataCacheMap *m) + { + return new data(this,m); + } + functionAdd(std::string _a,std::string _b) {strvec.push_back(_a);strvec.push_back(_b);} + void addNewTerm(std::string _a) { strvec.push_back(_a);} +}; + +*/ + int main (int argc, char* argv[]) { - + if (argc != 2){ printf("usage : elasticity input_file_name\n"); return -1; } +/* + fullMatrix<double> a(1,1); + a(0,0)=1.0; + fullMatrix<double> b(1,1); + b(0,0)=2.0; + fullMatrix<double> c(1,1); + c(0,0)=4.0; + fullMatrix<double> d(1,1); + d(0,0)=8.0; + fullMatrix<double> e(1,1); + e(0,0)=16.0; + + std::cout << argc << std::endl; + functionConstant fc_a(&a); + functionConstant fc_b(&b); + functionConstant fc_c(&c); + functionConstant fc_d(&d); + functionConstant fc_e(&e); + + functionAdd fa_ab(fc_a.getName(),fc_b.getName()); + functionAdd fa_cd(fc_c.getName(),fc_d.getName()); + functionAdd fa_abcd("a+b","c+d"); + + function::add("a+b", &fa_ab ); + function::add("c+d", &fa_cd); + function::add("a+b+c+d", &fa_abcd); + + dataCacheMap m(1); + std::cout << "start" << std::endl; + dataCacheDouble &dc_a=m.get(fc_a.getName()); + dataCacheDouble &dc_abcd=m.get("a+b+c+d"); + fa_abcd.addNewTerm(fc_e.getName()); + std::cout << "nb depsabcd=" << dc_abcd.howManyDoIDependOn() << std::endl; + std::cout << "nb depsa=" << dc_a.howManyDependOnMe() << std::endl; + std::cout << "a" << std::endl; + std::cout << dc_a(0,0) << std::endl; + std::cout << "a+b+c+d" << std::endl; + std::cout << dc_abcd(0,0) << std::endl; + std::cout << "dca.set(b)" << std::endl; + dc_a.set(b); + std::cout << "a+b+c+d" << std::endl; + std::cout << dc_abcd(0,0) << std::endl; +*/ + +/* + functionMult fm("axbxcxd","axbxcxd"); + dataCacheDouble *res; + res=fm.newDataCache(&m); + std::cout << "*res" << std::endl; + std::cout << (*res)(0,0) << std::endl; +*/ + +// return(0); - GmshInitialize(argc, argv); // globals are still present in Gmsh @@ -38,7 +128,7 @@ int main (int argc, char* argv[]) // stop gmsh GmshFinalize(); - + } @@ -60,13 +150,13 @@ int main (int argc, char* argv[]) ScalarLagrangeFunctionSpace L(100); std::cout << L.getNumKeys(e) << "fonctions de formes L" << std::endl; L.getKeys(e,dofs); - for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl; + for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl; dofs.clear(); L.f(e,0.1,0.1,0,vals); L.gradf(e,0.1,0.1,0,grads); std::copy(vals.begin(),vals.end(),output); std::cout << std::endl; for (std::vector<SVector3>::iterator it=grads.begin();it!=grads.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; } - + VectorLagrangeFunctionSpace L1(100,VectorLagrangeFunctionSpace::VECTOR_X); VectorLagrangeFunctionSpace L2(100,VectorLagrangeFunctionSpace::VECTOR_Y); std::cout << L2.getNumKeys(e) << "fonctions de formes L2" << std::endl; @@ -81,7 +171,7 @@ int main (int argc, char* argv[]) std::cout << P123.getNumKeys(e) << "fonctions de formes P123" << std::endl; P123.getKeys(e,dofs); std::cout << dofs.size() << std::endl; - for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl; + for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl; vals2.clear(); grads2.clear(); @@ -99,7 +189,7 @@ int main (int argc, char* argv[]) FormBilinear<TermBilinearMecaNL,ScalarLagrangeFunctionSpace,ScalarLagrangeFunctionSpace > fnl(L,L); fnl.func(); - + */ -- GitLab