- add LagrangeMultipliers

- split FunctionSpace / FunctionSpaceOfParent
- add Voids

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
 

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

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)};

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

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;
     }

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)

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_