diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index eb86ec0d4b697e34d1f9bf5123d79940085a7089..262a8329b0154283ddb242fd18b6eb1d23eb352e 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -17,20 +17,51 @@
#include "quadratureRules.h"
#include "solverField.h"
#include "MPoint.h"
+#include "DILevelset.h"
#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
#endif
+static void printLinearSystem(linearSystemCSRTaucs<double> *lsys)
+{
+ int *startIndex;
+ int *columns;
+ double *values;
+ lsys->getMatrix(startIndex, columns, values);
+ for(int i = 0; i < lsys->getNbUnk(); i++){
+ std::cout<<"a ";
+ for(int j = 0; j < lsys->getNbUnk(); j++){
+ double val = 0.;
+ for(int k = startIndex[i]; k < startIndex[i+1]; k++){
+ if(columns[k] == j) {val = values[k]; break;}
+ }
+ std::cout<< val <<" ";
+ }
+ std::cout<<std::endl;
+ }std::cout<<std::endl;
+ for(int i = 0; i < lsys->getNbUnk(); i++) {
+ double val; lsys->getFromRightHandSide(i,val);
+ std::cout<<"b "<<val<<std::endl;
+ }std::cout<<std::endl;
+ for(int i = 0; i < lsys->getNbUnk(); i++) {
+ double val; lsys->getFromSolution(i,val);
+ std::cout<<"x "<<val<<std::endl;
+ }
+}
+
void elasticitySolver::setMesh(const std::string &meshFileName)
{
pModel = new GModel();
pModel->readMSH(meshFileName.c_str());
_dim = pModel->getNumRegions() ? 3 : 2;
+
if (LagSpace) delete LagSpace;
- if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
- if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+ if (_dim==3) LagSpace=new VectorLagrangeFunctionSpaceOfParent(_tag);
+ if (_dim==2) LagSpace=new VectorLagrangeFunctionSpaceOfParent(_tag,VectorLagrangeFunctionSpaceOfParent::VECTOR_X,VectorLagrangeFunctionSpaceOfParent::VECTOR_Y);
+ if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
+ LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1);
}
void elasticitySolver::readInputFile(const std::string &fn)
@@ -47,12 +78,26 @@ void elasticitySolver::readInputFile(const std::string &fn)
field.g = new groupOfElements (_dim, physical);
elasticFields.push_back(field);
}
+ if (!strcmp(what, "LagrangeMultipliers")){
+ LagrangeMultiplierField field;
+ int physical;
+ double d1, d2, d3, val;
+ if(fscanf(f, "%d %lf %lf %lf %lf %lf %d", &physical, &field._tau,
+ &d1, &d2, &d3, &val, &field._tag) != 7) return;
+ field._tag = _tag;
+ field._d = SVector3(d1, d2, d3);
+ field._f = simpleFunction<double>(val);
+ field.g = new groupOfElements (_dim - 1, physical);
+ LagrangeMultiplierFields.push_back(field);
+ }
else if (!strcmp(what, "Void")){
- // elasticField field;
- // create enrichment ...
- // create the group ...
- // assign a tag
- // elasticFields.push_back(field);
+ elasticField field;
+ int physical;
+ if(fscanf(f, "%d", &physical) != 1) return;
+ field._E = field._nu = 0;
+ field.g = new groupOfElements (_dim, physical);
+ field._tag = 0;
+ elasticFields.push_back(field);
}
else if (!strcmp(what, "NodeDisplacement")){
double val;
@@ -139,6 +184,20 @@ void elasticitySolver::readInputFile(const std::string &fn)
if(fscanf(f, "%s", name) != 1) return;
setMesh(name);
}
+ else if (!strcmp(what, "CutMeshPlane")){
+ double a, b, c, d;
+ if(fscanf(f, "%lf %lf %lf %lf", &a, &b, &c, &d) != 4) return;
+ int tag=1;
+ gLevelsetPlane ls(a,b,c,d,tag);
+ pModel = pModel->buildCutGModel(&ls);
+ }
+ else if (!strcmp(what, "CutMeshSphere")){
+ double x, y, z, r;
+ if(fscanf(f, "%lf %lf %lf %lf", &x, &y, &z, &r) != 4) return;
+ int tag=1;
+ gLevelsetSphere ls(x,y,z,r,tag);
+ pModel = pModel->buildCutGModel(&ls);
+ }
else {
Msg::Error("Invalid input : %s", what);
// return;
@@ -164,23 +223,30 @@ void elasticitySolver::solve()
// give priority to fixations : when a dof is fixed, it cannot be
// numbered afterwards
- std::cout << "Dirichlet BC"<< std::endl;
+ // Dirichlet conditions
for (unsigned int i = 0; i < allDirichlet.size(); i++)
{
FilterDofComponent filter(allDirichlet[i]._comp);
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
- // again with another number.
+ // LagrangeMultipliers
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i)
+ {
+ NumberDofs(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(), LagrangeMultiplierFields[i].g->end(), *pAssembler);
+ }
+ // Elastic Fields
for (unsigned int i = 0; i < elasticFields.size(); ++i)
{
- NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
+ if(elasticFields[i]._E != 0.)
+ NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
}
-
- // Now we start the assembly process
- // First build the force vector
-
+ // Voids
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ if(elasticFields[i]._E == 0.)
+ FixVoidNodalDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(), *pAssembler);
+ }
+ // Neumann conditions
GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
std::cout << "Neumann BC"<< std::endl;
for (unsigned int i = 0; i < allNeumann.size(); i++)
@@ -188,20 +254,37 @@ void elasticitySolver::solve()
LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
}
+ // Assemble cross term, laplace term and rhs term for LM
+ GaussQuadrature Integ_LagrangeMult(GaussQuadrature::ValVal);
+ GaussQuadrature Integ_Laplace(GaussQuadrature::GradGrad);
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++)
+ {
+ LagrangeMultiplierTerm LagTerm(*LagSpace, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i]._d);
+ Assemble(LagTerm, *LagSpace, *LagrangeMultiplierSpace,
+ LagrangeMultiplierFields[i].g->begin(), LagrangeMultiplierFields[i].g->end(), Integ_LagrangeMult, *pAssembler);
-// bulk material law
+ LaplaceTerm<double,double> LapTerm(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i]._tau);
+ Assemble(LapTerm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(), LagrangeMultiplierFields[i].g->end(), Integ_Laplace, *pAssembler);
+ LoadTerm<double> Lterm(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i]._f);
+ Assemble(Lterm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(), LagrangeMultiplierFields[i].g->end(), Integ_Boundary, *pAssembler);
+ }
+ // Assemble elastic term for
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
-// LaplaceTerm<SVector3,SVector3> Eterm(*LagSpace);
Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
}
+
printf("nDofs=%d\n",pAssembler->sizeOfR());
+ printf("nFixed=%d\n",pAssembler->sizeOfF());
printf("-- done assembling!\n");
lsys->systemSolve();
printf("-- done solving!\n");
+
+ printLinearSystem(lsys);
+
double energ=0;
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
@@ -214,7 +297,6 @@ void elasticitySolver::solve()
}
-
#if defined(HAVE_POST)
static double vonMises(dofManager<double> *a, MElement *e,
@@ -259,42 +341,91 @@ static double vonMises(dofManager<double> *a, MElement *e,
PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
{
+ std::cout << "build Displacement View"<< std::endl;
std::set<MVertex*> v;
+ std::map<MVertex*,MElement*> vCut;
for (unsigned int i = 0; i < elasticFields.size(); ++i)
{
- for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
- {
+ if(elasticFields[i]._E == 0.) continue;
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it){
MElement *e=*it;
- for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
+ if(e->getParent()) {
+ for (int j = 0; j < e->getNumVertices(); ++j) {
+ if(vCut.find(e->getVertex(j)) == vCut.end())
+ vCut[e->getVertex(j)] = e->getParent();
+ }
+ }
+ else {
+ for (int j = 0; j < e->getNumVertices(); ++j)
+ v.insert(e->getVertex(j));
+ }
}
}
std::map<int, std::vector<double> > data;
SolverField<SVector3> Field(pAssembler, LagSpace);
- for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
- {
+ for (std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it){
SVector3 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;
}
+ for (std::map<MVertex*,MElement*>::iterator it = vCut.begin(); it != vCut.end(); ++it){
+ SVector3 val;
+ double uvw[3];
+ double xyz[3] = {it->first->x(), it->first->y(), it->first->z()};
+ it->second->xyz2uvw(xyz, uvw);
+ Field.f(it->second,uvw[0],uvw[1],uvw[2],val);
+ std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
+ data[it->first->getNum()]=vec;
+ }
+ PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
+ return pv;
+}
+
+PView* elasticitySolver::buildLagrangeMultiplierView (const std::string &postFileName)
+{
+ std::cout << "build Lagrange Multiplier View"<< std::endl;
+ if(!LagrangeMultiplierSpace) return new PView();
+ std::set<MVertex*> v;
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i)
+ {
+ for(groupOfElements::elementContainer::const_iterator it = LagrangeMultiplierFields[i].g->begin(); it != LagrangeMultiplierFields[i].g->end(); ++it)
+ {
+ MElement *e = *it;
+ for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
+ }
+ }
+ std::map<int, std::vector<double> > data;
+ SolverField<double> Field(pAssembler, LagrangeMultiplierSpace);
+ for(std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
+ {
+ double val;
+ MPoint p(*it);
+ Field.f(&p, 0, 0, 0, val);
+ std::vector<double> vec;
+ vec.push_back(val);
+ data[(*it)->getNum()] = vec;
+ }
PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
return pv;
}
PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
{
+ std::cout << "build Elastic Energy View"<< std::endl;
std::map<int, std::vector<double> > data;
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); ++i)
{
+ if(elasticFields[i]._E == 0.) continue;
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
ScalarTermConstant One(1.0);
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
{
- MElement *e=*it;
+ MElement *e = *it;
double energ;
double vol;
IntPt *GP;
@@ -303,7 +434,7 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
One.get(e,npts,GP,vol);
std::vector<double> vec;
vec.push_back(energ/vol);
- data[(*it)->getNum()]=vec;
+ data[e->getNum()]=vec;
}
}
PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
@@ -311,7 +442,7 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
}
PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
-{
+{std::cout << "build elastic view"<< std::endl;
std::map<int, std::vector<double> > data;
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); ++i)
@@ -337,7 +468,13 @@ PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
#else
-PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
+PView* elasticitySolver::buildDisplacementView(const std::string &postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
+
+PView* elasticitySolver::buildLagrangeMultiplierView (const std::string &postFileName)
{
Msg::Error("Post-pro module not available");
return 0;
@@ -349,5 +486,10 @@ PView* elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
return 0;
}
+PView* elasticitySolver::buildVonMisesView(const std::string &postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
#endif
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index fe1d228d5b61a1c7070a748ce3cda3ee77d410f8..b70f3c25383b7b45574ffc45b16cf3dbad0476b7 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -17,6 +17,15 @@ class GModel;
class PView;
class groupOfElements;
+struct LagrangeMultiplierField {
+ int _tag;
+ groupOfElements *g;
+ double _tau;
+ SVector3 _d;
+ simpleFunction<double> _f;
+ LagrangeMultiplierField() : g(0),_tag(0){}
+};
+
struct elasticField {
int _tag; // tag for the dofManager
groupOfElements *g; // support for this field
@@ -54,25 +63,29 @@ class elasticitySolver
int _dim, _tag;
dofManager<double> *pAssembler;
FunctionSpace<SVector3> *LagSpace;
+ FunctionSpace<double> *LagrangeMultiplierSpace;
// young modulus and poisson coefficient per physical
std::vector<elasticField> elasticFields;
+ std::vector<LagrangeMultiplierField> LagrangeMultiplierFields;
// neumann BC
std::vector<neumannBC> allNeumann;
// dirichlet BC
std::vector<dirichletBC> allDirichlet;
public:
- elasticitySolver(int tag) : _tag(tag),LagSpace(0),pAssembler(0) {}
+ elasticitySolver(int tag) : _tag(tag),LagSpace(0),pAssembler(0),LagrangeMultiplierSpace(0) {}
virtual ~elasticitySolver()
{
if (LagSpace) delete LagSpace;
+ if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
if (pAssembler) delete pAssembler;
}
void readInputFile(const std::string &meshFileName);
- void setMesh(const std::string &meshFileName);
+ virtual void setMesh(const std::string &meshFileName);
virtual void solve();
virtual PView *buildDisplacementView(const std::string &postFileName);
+ virtual PView *buildLagrangeMultiplierView(const std::string &posFileName);
virtual PView *buildElasticEnergyView(const std::string &postFileName);
virtual PView *buildVonMisesView(const std::string &postFileName);
// std::pair<PView *, PView*> buildErrorEstimateView
diff --git a/Solver/functionSpace.cpp b/Solver/functionSpace.cpp
index fed1efa9045ff667d19c0f1529b96ac7497f6b52..47067cd4ad8199db126da2018e7779bdb33ecb77 100644
--- a/Solver/functionSpace.cpp
+++ b/Solver/functionSpace.cpp
@@ -12,4 +12,4 @@
#include "functionSpace.h"
const SVector3 VectorLagrangeFunctionSpace::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};
-
+const SVector3 VectorLagrangeFunctionSpaceOfParent::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 53ea4afd43807e5aa2aa78beb1a753c7fc0cac3e..e5bb028499d6d57f8da6aaa1eaa57948f8703a23 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -54,12 +54,9 @@ class FunctionSpace : public FunctionSpaceBase
typedef typename TensorialTraits<T>::GradType GradType;
typedef typename TensorialTraits<T>::HessType HessType;
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 void gradfuvw(MElement *ele, double u, double v, double w,
- std::vector<GradType> &grads) {} // should return to pure virtual once all is done.
- virtual void hessfuvw(MElement *ele, double u, double v, double w,
- std::vector<HessType> &hess) = 0;
+ virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) = 0;
+ virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) {} // should return to pure virtual once all is done.
+ virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess) = 0;
virtual int getNumKeys(MElement *ele) = 0; // if one needs the number of dofs
virtual void getKeys(MElement *ele, std::vector<Dof> &keys) = 0;
};
@@ -85,17 +82,24 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
virtual int getId(void) const {return _iField;}
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
- if (ele->getParent()) ele = ele->getParent();
+ if(ele->getParent()) {
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ ele->movePointFromParentSpaceToElementSpace(u, v, w);
+ }
+ }
int ndofs = ele->getNumVertices();
int curpos = vals.size();
vals.resize(curpos + ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
}
-
// Fonction renvoyant un vecteur contenant le grandient de chaque FF
virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
- if (ele->getParent()) ele = ele->getParent();
+ if(ele->getParent()) {
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ ele->movePointFromParentSpaceToElementSpace(u, v, w);
+ }
+ }
int ndofs = ele->getNumVertices();
grads.reserve(grads.size() + ndofs);
double gradsuvw[256][3];
@@ -110,9 +114,14 @@ 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]));
}
- // Fonction renvoyant un vecteur contenant le hessien [][] de chaque FF dans l'espace ISOPARAMETRIQUE
+ // Fonction renvoyant un vecteur contenant le hessien [][] de chaque FF dans l'espace ISOPARAMETRIQUE
virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess)
{
+ if(ele->getParent()) {
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ ele->movePointFromParentSpaceToElementSpace(u, v, w);
+ }
+ }
int ndofs = ele->getNumVertices(); // ATTENTION RETOURNE LE NBBRE DE NOEUDS ET PAS LE NBRE DE DDL
hess.reserve(hess.size() + ndofs); // permet de mettre les composantes suivantes à la suite du vecteur
double hessuvw[256][3][3];
@@ -128,27 +137,26 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
- if (ele->getParent()) ele = ele->getParent();
+ if(ele->getParent()) {
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ ele->movePointFromParentSpaceToElementSpace(u, v, w);
+ }
+ }
int ndofs = ele->getNumVertices();
grads.reserve(grads.size() + ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
for (int i = 0; i < ndofs; ++i)
- grads.push_back(GradType(
- gradsuvw[i][0],
- gradsuvw[i][1],
- gradsuvw[i][2]));
+ grads.push_back(GradType(gradsuvw[i][0], gradsuvw[i][1], gradsuvw[i][2]));
}
virtual int getNumKeys(MElement *ele)
{
- if (ele->getParent()) ele = ele->getParent();
return ele->getNumVertices();
}
virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
{
- if (ele->getParent()) ele = ele->getParent();
int ndofs = ele->getNumVertices();
keys.reserve(keys.size() + ndofs);
for (int i = 0; i < ndofs; ++i)
@@ -156,6 +164,92 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
}
};
+class ScalarLagrangeFunctionSpaceOfParent : public FunctionSpace<double>
+{
+ public:
+ typedef TensorialTraits<double>::ValType ValType;
+ typedef TensorialTraits<double>::GradType GradType;
+ typedef TensorialTraits<double>::HessType HessType;
+
+ 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:
+ ScalarLagrangeFunctionSpaceOfParent(int i = 0) : _iField(i) {}
+ virtual int getId(void) const {return _iField;}
+ virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+ {
+ if(ele->getParent()) ele = ele->getParent();
+ int ndofs = ele->getNumVertices();
+ int curpos = vals.size();
+ vals.resize(curpos + ndofs);
+ ele->getShapeFunctions(u, v, w, &(vals[curpos]));
+ }
+
+ // Fonction renvoyant un vecteur contenant le grandient de chaque FF
+ virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
+ {
+ if(ele->getParent()) ele = ele->getParent();
+ int ndofs = ele->getNumVertices();
+ grads.reserve(grads.size() + ndofs);
+ double gradsuvw[256][3];
+ ele->getGradShapeFunctions(u, v, w, gradsuvw);
+ double jac[3][3];
+ double invjac[3][3];
+ const double detJ = ele->getJacobian(u, v, w, jac); // redondant : on fait cet appel a l'exterieur
+ inv3x3(jac, invjac);
+ for (int i = 0; i < ndofs; ++i)
+ grads.push_back(GradType(
+ invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2],
+ 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]));
+ }
+ // Fonction renvoyant un vecteur contenant le hessien [][] de chaque FF dans l'espace ISOPARAMETRIQUE
+ virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess)
+ {
+ if(ele->getParent()) ele = ele->getParent();
+ int ndofs = ele->getNumVertices(); // ATTENTION RETOURNE LE NBBRE DE NOEUDS ET PAS LE NBRE DE DDL
+ hess.reserve(hess.size() + ndofs); // permet de mettre les composantes suivantes à la suite du vecteur
+ double hessuvw[256][3][3];
+ ele->getHessShapeFunctions(u, v, w, hessuvw);
+ HessType hesst;
+ for (int i = 0; i < ndofs; ++i){
+ hesst(0,0) = hessuvw[i][0][0]; hesst(0,1) = hessuvw[i][0][1]; hesst(0,2) = hessuvw[i][0][2];
+ hesst(1,0) = hessuvw[i][1][0]; hesst(1,1) = hessuvw[i][1][1]; hesst(1,2) = hessuvw[i][1][2];
+ hesst(2,0) = hessuvw[i][2][0]; hesst(2,1) = hessuvw[i][2][1]; hesst(2,2) = hessuvw[i][2][2];
+ hess.push_back(hesst);
+ }
+ }
+ virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
+ {
+ if(ele->getParent()) ele = ele->getParent();
+ int ndofs = ele->getNumVertices();
+ grads.reserve(grads.size() + ndofs);
+ double gradsuvw[256][3];
+ ele->getGradShapeFunctions(u, v, w, gradsuvw);
+ for (int i = 0; i < ndofs; ++i)
+ grads.push_back(GradType(gradsuvw[i][0], gradsuvw[i][1], gradsuvw[i][2]));
+ }
+ virtual int getNumKeys(MElement *ele)
+ {
+ if(ele->getParent()) ele = ele->getParent();
+ return ele->getNumVertices();
+ }
+
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys) // appends ...
+ {
+ if(ele->getParent()) ele = ele->getParent();
+ int ndofs = ele->getNumVertices();
+ keys.reserve(keys.size() + ndofs);
+ for (int i = 0; i < ndofs; ++i)
+ getKeys(ele->getVertex(i), keys);
+ }
+};
template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T>
{
@@ -295,6 +389,31 @@ class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
{}
};
+class VectorLagrangeFunctionSpaceOfParent : public ScalarToAnyFunctionSpace<SVector3>
+{
+ protected:
+ static const SVector3 BasisVectors[3];
+ public:
+ enum Along { VECTOR_X = 0, VECTOR_Y = 1, VECTOR_Z = 2 };
+ typedef TensorialTraits<SVector3>::ValType ValType;
+ typedef TensorialTraits<SVector3>::GradType GradType;
+ VectorLagrangeFunctionSpaceOfParent(int id) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfParent(id),
+ SVector3(1.,0.,0.), VECTOR_X, SVector3(0.,1.,0.), VECTOR_Y, SVector3(0.,0.,1.), VECTOR_Z)
+ {}
+ VectorLagrangeFunctionSpaceOfParent(int id,Along comp1) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfParent(id),
+ BasisVectors[comp1], comp1)
+ {}
+ VectorLagrangeFunctionSpaceOfParent(int id,Along comp1,Along comp2) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfParent(id),
+ BasisVectors[comp1], comp1, BasisVectors[comp2], comp2)
+ {}
+ VectorLagrangeFunctionSpaceOfParent(int id,Along comp1,Along comp2, Along comp3) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeFunctionSpaceOfParent(id),
+ BasisVectors[comp1], comp1, BasisVectors[comp2], comp2, BasisVectors[comp3], comp3)
+ {}
+};
template<class T>
class CompositeFunctionSpace : public FunctionSpace<T>
@@ -489,7 +608,7 @@ template <class T> int xFemFunctionSpace<T>::getNumKeys(MElement *ele)
MElement * elep;
if (ele->getParent()) elep = ele->getParent();
else elep = ele;
- int nbdofs = xFemFunctionSpace<T>::_spacebase->getNumKeys(ele);
+ int nbdofs = xFemFunctionSpace<T>::_spacebase->getNumKeys(elep);
return nbdofs;
}
@@ -617,7 +736,4 @@ template <class T,class F> void FilteredFunctionSpace<T,F>::getKeys(MElement *el
}
-
-
-
#endif
diff --git a/Solver/quadratureRules.h b/Solver/quadratureRules.h
index 4ae9438c0edeba1d25f41b4acd54a3d9252107da..12c12b3dc006f418276875438f1d1a67605ed527 100644
--- a/Solver/quadratureRules.h
+++ b/Solver/quadratureRules.h
@@ -56,7 +56,7 @@ class GaussQuadrature : public QuadratureBase
integrationOrder=2*geoorder;
break;
case GradGrad :
- integrationOrder=3*(geoorder-1);
+ integrationOrder=3*(geoorder-1)+1;
break;
default : integrationOrder=1;
}
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index ac501538e771c2aa840db4a65e5c1cc60569f2ef..33f5a2b88665b3f20f998fdbc9b85378207089d9 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -38,42 +38,44 @@ template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,Fu
}
}
+template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler) // symmetric
+{
+ fullMatrix<typename Assembler::dataMat> localMatrix;
+ std::vector<Dof> R;
+ IntPt *GP;
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localMatrix);
+ space.getKeys(e,R);
+ assembler.assemble(R, localMatrix);
+}
+
template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,
- FunctionSpaceBase &shapeFcts,
- FunctionSpaceBase &testFcts,
- Iterator itbegin,
- Iterator itend,
- QuadratureBase &integrator,
- Assembler &assembler) // non symmetric
+ FunctionSpaceBase &shapeFcts,
+ FunctionSpaceBase &testFcts,
+ Iterator itbegin,
+ Iterator itend,
+ QuadratureBase &integrator,
+ Assembler &assembler) // non symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
std::vector<Dof> R;
std::vector<Dof> C;
- for (Iterator it = itbegin;it!=itend; ++it)
+ for (Iterator it = itbegin; it != itend; ++it)
{
MElement *e = *it;
R.clear();
C.clear();
IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localMatrix);
- shapeFcts.getKeys(e,R);
- testFcts.getKeys(e,C);
+ int npts = integrator.getIntPoints(e, &GP);
+ term.get(e, npts, GP, localMatrix);
+ shapeFcts.getKeys(e, R);
+ testFcts.getKeys(e, C);
assembler.assemble(R, C, localMatrix);
+ assembler.assemble(C, R, localMatrix.transpose());
}
}
-template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler) // symmetric
-{
- fullMatrix<typename Assembler::dataMat> localMatrix;
- std::vector<Dof> R;
- IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localMatrix);
- space.getKeys(e,R);
- assembler.assemble(R, localMatrix);
-}
template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler)
{
@@ -125,7 +127,6 @@ template<class Iterator,class dataMat> void Assemble(ScalarTermBase &term,MEleme
}
-
template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &dofs,std::vector<typename Assembler::dataVec> &vals)
{
int nbff=dofs.size();
@@ -141,7 +142,7 @@ class FilterDof
virtual bool operator()(Dof key)=0;
};
-class FilterDofTrivial
+class FilterDofTrivial : public FilterDof
{
public:
virtual bool operator()(Dof key) {return true;}
@@ -197,6 +198,13 @@ template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &sp
}
}
+template<class Iterator,class Assembler> void FixVoidNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
+{
+ FilterDofTrivial filter;
+ simpleFunction<typename Assembler::dataVec> fct(0.);
+ FixNodalDofs(space, itbegin, itend, assembler, fct, filter);
+}
+
template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
{
for (Iterator it=itbegin;it!=itend;++it)
diff --git a/Solver/terms.h b/Solver/terms.h
index f629943724b894fefb515d4c9bae263c0e1eda86..7a4ba1bdfd0fa8f912472c5d2d3fb86098e9dd30 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -35,7 +35,7 @@ template<class T1,class T2> class BilinearTerm : public BilinearTermBase
FunctionSpace<T1>& space1;
FunctionSpace<T2>& space2;
public :
- BilinearTerm(FunctionSpace<T1>& space1_,FunctionSpace<T1>& space2_) : space1(space1_),space2(space2_) {}
+ BilinearTerm(FunctionSpace<T1>& space1_,FunctionSpace<T2>& space2_) : space1(space1_),space2(space2_) {}
virtual ~BilinearTerm() {}
};
@@ -129,8 +129,9 @@ template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symmetric
{
+ double diffusivity;
public :
- LaplaceTerm(FunctionSpace<T1>& space1_) : BilinearTerm<T1,T1>(space1_,space1_)
+ LaplaceTerm(FunctionSpace<T1>& space1_, double diff=1) : BilinearTerm<T1,T1>(space1_,space1_), diffusivity(diff)
{}
virtual ~LaplaceTerm() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
@@ -149,7 +150,7 @@ template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symm
{
for (int k = j; k < nbFF; k++)
{
- double contrib=weight * detJ * dot(Grads[j],Grads[k]);
+ double contrib=weight * detJ * dot(Grads[j],Grads[k]) * diffusivity;
m(j,k)+=contrib;
if (j!=k) m(k,j)+=contrib;
}
@@ -161,8 +162,6 @@ template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symm
}; // class
-
-
class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{
protected :
@@ -202,8 +201,8 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
virtual ~IsotropicElasticTerm() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
{
- if (ele->getParent()) ele=ele->getParent();
- if (sym)
+ if (ele->getParent()) ele=ele->getParent();
+ if (sym)
{
int nbFF = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
double jac[3][3];
@@ -217,7 +216,6 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
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);
std::vector<TensorialTraits<SVector3>::GradType> Grads;
BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads); // a optimiser ??
@@ -285,7 +283,6 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
inline double dot(const double &a, const double &b)
{ return a*b; }
-
template<class T1> class LoadTerm : public LinearTerm<T1>
{
simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
@@ -295,29 +292,56 @@ template<class T1> class LoadTerm : public LinearTerm<T1>
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
{
- if (ele->getParent()) ele=ele->getParent();
- int nbFF=LinearTerm<T1>::space1.getNumKeys(ele);
+ if (ele->getParent()) ele=ele->getParent();
+ int nbFF=LinearTerm<T1>::space1.getNumKeys(ele);
double jac[3][3];
m.resize(nbFF);
m.scale(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);
+ 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);
std::vector<typename TensorialTraits<T1>::ValType> Vals;
- LinearTerm<T1>::space1.f(ele,u, v, w, Vals);
+ LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
SPoint3 p;
ele->pnt(u, v, w, p);
typename TensorialTraits<T1>::ValType load=Load(p.x(),p.y(),p.z());
for (int j = 0; j < nbFF ; ++j)
{
- m(j)+=dot(Vals[j],load)*weight*detJ;
+ m(j) += dot(Vals[j], load) * weight * detJ;
}
}
}
};
-
-
+class LagrangeMultiplierTerm : public BilinearTerm<SVector3,double>
+{
+ SVector3 _d;
+ public :
+ LagrangeMultiplierTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<double>& space2_, const SVector3 &d) :
+ BilinearTerm<SVector3,double>(space1_, space2_) {for(int i=0; i < 3; i++) _d(i) = d(i);}
+ virtual ~LagrangeMultiplierTerm() {}
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
+ {
+ int nbFF1 = BilinearTerm<SVector3,double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+ int nbFF2 = BilinearTerm<SVector3,double>::space2.getNumKeys(ele); //nbVertices of boundary
+ double jac[3][3];
+ m.resize(nbFF1, nbFF2);
+ m.setAll(0.);
+ for (int i = 0; i < npts; i++) {
+ double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
+ const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+ std::vector<TensorialTraits<SVector3>::ValType> Vals;
+ std::vector<TensorialTraits<double>::ValType> ValsT;
+ BilinearTerm<SVector3,double>::space1.f(ele, u, v, w, Vals);
+ BilinearTerm<SVector3,double>::space2.f(ele, u, v, w, ValsT);
+ for (int j = 0; j < nbFF1; j++) {
+ for (int k = 0; k < nbFF2; k++) {
+ m(j, k) += dot(Vals[j], _d) * ValsT[k] * weight * detJ;
+ }
+ }
+ }
+ }
+};
#endif// _TERMS_H_