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

No commit message 

No commit message
parent 31bd17d0
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Solver/elasticitySolver.cpp
+ 27
− 3
View file @ f8ee7c1c
......@@ -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;
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;
......
Solver/elasticitySolver.h
+ 2
− 1
View file @ f8ee7c1c
......@@ -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
......
Solver/functionSpace.h
+ 30
− 30
View file @ f8ee7c1c
......@@ -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);
};
......
Solver/terms.h
+ 18
− 45
View file @ f8ee7c1c
......@@ -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");
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment