diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index e7e5575caee7d26bf71b5550f87648f53429e877..b69c5ef74732b007a3a98186d5a5795f690cf6f1 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -15,6 +15,7 @@
#include "terms.h"
#include "solverAlgorithms.h"
#include "quadratureRules.h"
+#include "solverField.h"
#if defined(HAVE_POST)
#include "PView.h"
@@ -285,9 +286,32 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
return pv;
}
-PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
+PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
{
- std::map<int, std::vector<double> > data;
+ std::map<int, std::vector<double> > data;
+ GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+ double energ=0;
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ IsotropicElasticTerm<SolverField<SVector3> ,SolverField<SVector3> > Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ ScalarTermOne One;
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
+ {
+ MElement *e=*it;
+ fullMatrix<double> localMatrix;
+ IntPt *GP;
+ int npts=Integ_Bulk.getIntPoints(e,&GP);
+ Eterm.get(e,npts,GP,localMatrix);
+ double vol=0;
+ One.get(e,npts,GP,vol);
+ std::vector<double> vec;
+ vec.push_back(localMatrix(0,0)/vol);
+ energ+=localMatrix(0,0);
+ data[(*it)->getNum()]=vec;
+ }
+ }
+ std::cout<< "elastic energy=" << energ << std::endl;
PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
return pv;
}
@@ -300,7 +324,7 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName
return 0;
}
-PView* elasticitySolver::buildVonMisesView(const std::string &postFileName)
+PView* elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
{
Msg::Error("Post-pro module not available");
return 0;
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index d9b82a8fe77b09495bf6c1c4a46c5949601e297c..ef5109181f3ade19094980b8dad58d5e0152cd29 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -74,7 +74,8 @@ class elasticitySolver
void setMesh(const std::string &meshFileName);
virtual void solve();
virtual PView *buildDisplacementView(const std::string &postFileName);
- virtual PView *buildVonMisesView(const std::string &postFileName);
+ virtual PView *buildElasticEnergyView(const std::string &postFileName);
+ // virtual PView *buildVonMisesView(const std::string &postFileName);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int);
// std::pair<PView *, PView*> buildErrorEstimateView
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 659037c2f4a3414d72a0648c25d704b8e3658de8..bfe5df12734158ab16ad27cc652455844188906d 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -67,9 +67,9 @@ class FunctionSpaceBase
{
public:
virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0;
virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
- virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
};
template<class T>
@@ -83,18 +83,18 @@ class FunctionSpace : public FunctionSpaceBase
typedef typename TensorialTraits<T>::CurlType CurlType;*/
public:
- virtual int f(MVertex *ver, std::vector<ValType> &vals) {} // used for neumann BC
- 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 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 int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0;
virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
- virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
};
class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
@@ -113,15 +113,15 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
public:
ScalarLagrangeFunctionSpace(int i=0):_iField(i) {}
virtual int getId(void) const {return _iField;};
- virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+ virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
int ndofs= ele->getNumVertices();
int curpos=vals.size();
vals.resize(curpos+ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
};
- virtual int f(MVertex *ver, std::vector<ValType> &vals) {vals.push_back(1.0);} // used for neumann BC
- virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+ 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();
grads.reserve(grads.size()+ndofs);
@@ -140,7 +140,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
};
virtual int getNumKeys(MElement *ele) {return ele->getNumVertices();}
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
{
int ndofs= ele->getNumVertices();
keys.reserve(keys.size()+ndofs);
@@ -148,7 +148,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
getKeys(ele->getVertex(i),keys);
}
virtual int getNumKeys(MVertex *ver) { return 1;}
- virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
{
keys.push_back(Dof(ver->getNum(), _iField));
};
@@ -187,7 +187,7 @@ public :
}
virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;}
- virtual int f(MVertex *ver, std::vector<ValType> &vals)
+ virtual void f(MVertex *ver, std::vector<ValType> &vals)
{
std::vector<double> valsd;
ScalarFS->f(ver,valsd);
@@ -201,7 +201,7 @@ public :
}
} // used for neumann BC
- virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+ virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
std::vector<double> valsd;
ScalarFS->f(ele,u,v,w,valsd);
@@ -215,7 +215,7 @@ public :
}
}
- virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+ virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
std::vector<SVector3> gradsd;
ScalarFS->gradf(ele,u,v,w,gradsd);
@@ -236,7 +236,7 @@ public :
virtual int getNumKeys(MElement *ele) {return ScalarFS->getNumKeys(ele)*comp.size();}
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys)
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys)
{
int nk=ScalarFS->getNumKeys(ele);
@@ -259,7 +259,7 @@ public :
}
virtual int getNumKeys(MVertex *ver) { return ScalarFS->getNumKeys(ver)*comp.size();}
- virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
{
int nk=ScalarFS->getNumKeys(ver);
std::vector<Dof> bufk;
@@ -348,19 +348,19 @@ class CompositeFunctionSpace : public FunctionSpace<T>
delete (*it);
}
- virtual int f(MVertex *ver, std::vector<ValType> &vals)
+ virtual void f(MVertex *ver, std::vector<ValType> &vals)
{
for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
(*it)->f(ver,vals);
}
- virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &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)
(*it)->f(ele,u,v,w,vals);
}
- virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+ virtual void 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);
@@ -374,7 +374,7 @@ class CompositeFunctionSpace : public FunctionSpace<T>
return ndofs;
}
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys)
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys)
{
for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
(*it)->getKeys(ele,keys);
@@ -388,7 +388,7 @@ class CompositeFunctionSpace : public FunctionSpace<T>
return ndofs;
}
- virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
{
for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
(*it)->getKeys(ver,keys);
@@ -408,13 +408,13 @@ class xFemFunctionSpace : public FunctionSpace<T>
FunctionSpace<T>* _spacebase;
// Function<double>* enrichment;
public:
- virtual int f(MVertex *ver, std::vector<ValType> &vals);
- 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 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 int getKeys(MVertex *ver, std::vector<Dof> &keys);
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
};
@@ -432,13 +432,13 @@ class FilteredFunctionSpace : public FunctionSpace<T>
F &_filter;
public:
- virtual int f(MVertex *ver, std::vector<ValType> &vals);
- 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 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 int getKeys(MVertex *ver, std::vector<Dof> &keys);
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
};
diff --git a/Solver/terms.h b/Solver/terms.h
index b5734be749fba741b64dbafdb670a3e57bf96e0c..dff4c69fd0bcab317f0796aa2b0b653fea6579e9 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -21,46 +21,6 @@
#include "functionSpace.h"
#include "groupOfElements.h"
-
-
-// evaluation of a field ???
-template<class T>
-class Field {
- protected:
- typedef typename TensorialTraits<T>::ValType ValType;
- typedef typename TensorialTraits<T>::GradType GradType;
- dofManager<double> *dm; //
-/* typedef typename TensorialTraits<T>::HessType HessType;
- typedef typename TensorialTraits<T>::DivType DivType;
- typedef typename TensorialTraits<T>::CurlType CurlType;*/
-
- public:
- virtual int f(MElement *ele, double u, double v, double w, ValType &val)=0;
- virtual int gradf(MElement *ele, double u, double v, double w,GradType &grad)=0;
-// 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 int getKeys(MElement *ele, Dof *keys)=0; // may be faster once the number of dofs is known
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
-};
-
-
-
-
-
-
-class Formulation
-{
- std::vector<FunctionSpace<double>* > scalarfs;
- std::vector<FunctionSpace<SVector3>* > vectorfs;
- std::vector<groupOfElements* > groups;
- std::vector<std::pair<MElement*,std::vector<groupOfElements&> > > links;
- dofManager<double> *dm; //
-
-};
-
class BilinearTermBase
{
public :
@@ -74,14 +34,13 @@ template<class S1,class S2> class BilinearTerm : public BilinearTermBase
S2& space2;
public :
BilinearTerm(S1& space1_,S2& space2_) : space1(space1_),space2(space2_) {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
};
class LinearTermBase
{
public:
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
- virtual void get(MVertex *v,fullVector<double> &m) =0;
+ virtual void get(MVertex *ver,fullVector<double> &m) =0;
};
template<class S1> class LinearTerm : public LinearTermBase
@@ -90,7 +49,6 @@ template<class S1> class LinearTerm : public LinearTermBase
S1& space1;
public :
LinearTerm(S1& space1_) : space1(space1_) {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) =0;
};
class ScalarTermBase
@@ -102,7 +60,22 @@ class ScalarTermBase
class ScalarTerm : public ScalarTermBase
{
public :
- virtual void get(MElement *ele,int npts,IntPt *GP,double &v) =0;
+};
+
+class ScalarTermOne : public ScalarTerm
+{
+ public :
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
+ {
+ double jac[3][3];
+ val=0;
+ for (int i = 0; i < npts; i++)
+ {
+ const double u = GP[i].pt[0];const double v = GP[i].pt[1];const double w = GP[i].pt[2];
+ const double weight = GP[i].weight;const double detJ = ele->getJacobian(u, v, w, jac);
+ val+=weight*detJ;
+ }
+ }
};
@@ -117,7 +90,7 @@ template<class S1,class S2> class LaplaceTerm : public BilinearTerm<S1,S2>
{
Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
}
- virtual void get(MVertex *v,fullMatrix<double> &m)
+ virtual void get(MVertex *ver,fullMatrix<double> &m)
{
Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
}