thermic solver

CMakeLists.txt

@@ -109,6 +109,7 @@ set(GMSH_API
   Numeric/Numeric.h Numeric/GaussIntegration.h Numeric/polynomialBasis.h
     Numeric/JacobianBasis.h Numeric/MetricBasis.h Numeric/bezierBasis.h Numeric/fullMatrix.h
     Numeric/simpleFunction.h Numeric/cartesian.h Numeric/ElementType.h
+    Numeric/BasisFactory.h
   Geo/GModel.h Geo/GEntity.h Geo/GPoint.h Geo/GVertex.h Geo/GEdge.h 
     Geo/GFace.h Geo/GRegion.h Geo/GEdgeLoop.h Geo/GEdgeCompound.h 
     Geo/GFaceCompound.h Geo/GRegionCompound.h Geo/GRbf.h Geo/MVertex.h
@@ -131,6 +132,7 @@ set(GMSH_API
     Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h 
     Solver/linearSystemFull.h Solver/elasticitySolver.h Solver/sparsityPattern.h 
     Solver/groupOfElements.h Solver/linearSystemPETSc.h Solver/linearSystemMUMPS.h
+    Solver/thermicSolver.h
   Post/PView.h Post/PViewData.h Plugin/PluginManager.h Post/OctreePost.h 
   Post/PViewDataList.h Post/PViewDataGModel.h Post/PViewOptions.h Post/ColorTable.h
    Numeric/nodalBasis.h

Solver/CMakeLists.txt

@@ -12,6 +12,7 @@ set(SRC
   groupOfElements.cpp
   elasticityTerm.cpp
 elasticitySolver.cpp
+  thermicSolver.cpp
   SElement.cpp
   eigenSolver.cpp
   multiscaleLaplace.cpp

Solver/elasticitySolver.h

@@ -22,7 +22,7 @@ struct LagrangeMultiplierField {
   groupOfElements *g;
   double _tau;
   SVector3 _d;
-  simpleFunction<double> _f;
+  simpleFunction<double> *_f;
   LagrangeMultiplierField() : _tag(0), g(0){}
 };
 
@@ -79,8 +79,11 @@ class elasticitySolver
   elasticitySolver(GModel *model, int tag);
 
   void addDirichletBC(int dim, int entityId, int component, double value);
+  void addDirichletBC(int dim, std::string phys, int component, double value);
   void addNeumannBC(int dim, int entityId, const std::vector<double> value);
+  void addNeumannBC(int dim, std::string phys, const std::vector<double> value);
   void addElasticDomain(int tag, double e, double nu);
+  void addElasticDomain(std::string phys, double e, double nu);
 
   virtual ~elasticitySolver()
   {
@@ -95,12 +98,15 @@ class elasticitySolver
   void solve();
   void postSolve();
   void exportKb();
-  void getSolutionOnElement(MElement *el, fullMatrix<double> &sol);
+  void computeEffectiveStiffness(std::vector<double> stiff);
+  void computeEffectiveStrain(std::vector<double> strain);
   virtual PView *buildDisplacementView(const std::string postFileName);
+  virtual PView *buildStrainView(const std::string postFileName);
   virtual PView *buildStressesView(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);
+  virtual PView *buildVolumeView(const std::string postFileName);
   // std::pair<PView *, PView*> buildErrorEstimateView
   //   (const std::string &errorFileName, double, int);
   // std::pair<PView *, PView*> buildErrorEstimateView

Solver/femTerm.h

@@ -118,23 +118,13 @@ class femTerm {
       dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
   }
 
-  void neumannNodalBC(int physical, int dim, int comp, int field,
+  void neumannNodalBC(MElement *e, int comp, int field,
                       const simpleFunction<dataVec> &fct,
                       dofManager<dataVec> &dm)
   {
-    std::map<int, std::vector<GEntity*> > groups[4];
-    GModel *m = _gm;
-    m->getPhysicalGroups(groups);
-    std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
-    if (it == groups[dim].end()) return;
     double jac[3][3];
     double sf[256];
-    for (unsigned int i = 0; i < it->second.size(); ++i){
-      GEntity *ge = it->second[i];
-      for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
-        MElement *e = ge->getMeshElement(j);
     int integrationOrder = 2 * e->getPolynomialOrder();
-        int nbNodes = e->getNumVertices();
     int npts;
     IntPt *GP;
     e->getIntegrationPoints(integrationOrder, &npts, &GP);
@@ -147,10 +137,46 @@ class femTerm {
       SPoint3 p; e->pnt(u, v, w, p);
       e->getShapeFunctions(u, v, w, sf);
       const dataVec FCT = fct(p.x(), p.y(), p.z());
-          for (int k = 0; k < nbNodes; k++){
-            dm.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
+      for (int k = 0; k < e->getNumShapeFunctions(); k++){
+        dm.assemble(e->getShapeFunctionNode(k), comp, field, detJ * weight * sf[k] * FCT);
+      }
     }
   }
+
+  void neumannNodalBC(int physical, int dim, int comp, int field,
+                      const simpleFunction<dataVec> &fct,
+                      dofManager<dataVec> &dm)
+  {
+    std::map<int, std::vector<GEntity*> > groups[4];
+    GModel *m = _gm;
+    m->getPhysicalGroups(groups);
+    std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
+    if (it == groups[dim].end()) return;
+    for (unsigned int i = 0; i < it->second.size(); ++i){
+      GEntity *ge = it->second[i];
+      for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
+        MElement *e = ge->getMeshElement(j);
+        neumannNodalBC(e, comp, field, fct, dm);
+      }
+    }
+  }
+  void neumannNormalNodalBC(int physical, int dim, int field,
+                            const simpleFunction<dataVec> &fct,
+                            dofManager<dataVec> &dm)
+  {
+    std::map<int, std::vector<GEntity*> > groups[4];
+    GModel *m = _gm;
+    m->getPhysicalGroups(groups);
+    std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
+    if (it == groups[dim].end()) return;
+    for (unsigned int i = 0; i < it->second.size(); ++i){
+      GEntity *ge = it->second[i];
+      for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
+        MElement *e = ge->getMeshElement(j);
+        
+        neumannNodalBC(e, 0, field, fct, dm);
+        neumannNodalBC(e, 1, field, fct, dm);
+        neumannNodalBC(e, 2, field, fct, dm);
       }
     }
   }

Solver/functionSpace.h

@@ -115,7 +115,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
   virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
   {
     if(ele->getParent()) {
-      if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+      if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+         ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
         ele->movePointFromParentSpaceToElementSpace(u, v, w);
       }
     }
@@ -128,7 +129,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
   virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
   {
     if(ele->getParent()) {
-      if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+      if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+         ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
         ele->movePointFromParentSpaceToElementSpace(u, v, w);
       }
     }
@@ -150,7 +152,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
   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...
+      if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+         ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
         ele->movePointFromParentSpaceToElementSpace(u, v, w);
       }
     }
@@ -169,7 +172,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
   virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
   {
     if(ele->getParent()) {
-      if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+      if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+         ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
         ele->movePointFromParentSpaceToElementSpace(u, v, w);
       }
     }

Solver/terms.cpp

@@ -142,31 +142,6 @@ void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<do
   }
 }
 
-void LagrangeMultiplierTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
-{
-  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;
-      }
-    }
-  }
-}
-
 void LagMultTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
 {
   int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent

Solver/terms.h

@@ -34,8 +34,14 @@ inline double dot(const double &a, const double &b)
 
 inline int delta(int i,int j) {if (i==j) return 1; else return 0;}
 
-
-
+inline void setDirection(double &a, const double &b)
+{
+  a = b;
+}
+inline void setDirection(SVector3 &a, const SVector3 &b)
+{
+  for(int i = 0; i < 3; i++) a(i) = b(i);
+}
 
 template<class T2=double> class  ScalarTermBase
 {
@@ -234,9 +240,9 @@ public :
 template<class T1> class LoadTerm : public LinearTerm<T1>
 {
 protected:
-  simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
+  simpleFunction<typename TensorialTraits<T1>::ValType> *Load;
 public :
-  LoadTerm(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :
+  LoadTerm(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> *Load_) :
       LinearTerm<T1>(space1_), Load(Load_) {}
   virtual ~LoadTerm() {}
   virtual LinearTermBase<double>* clone () const { return new LoadTerm<T1>(LinearTerm<T1>::space1,Load);}
@@ -244,19 +250,19 @@ public :
   virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<double> > &vv) const {};
 };
 
-class LagrangeMultiplierTerm : public BilinearTerm<SVector3, double>
+template<class T1> class LagrangeMultiplierTerm : public BilinearTerm<T1, double>
 {
-  SVector3 _d;
+  T1 _d;
 public :
-  LagrangeMultiplierTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<double>& space2_, const SVector3 &d) :
-      BilinearTerm<SVector3, double>(space1_, space2_)
+  LagrangeMultiplierTerm(FunctionSpace<T1>& space1_, FunctionSpace<double>& space2_, const T1 &d) :
+      BilinearTerm<T1, double>(space1_, space2_)
   {
-    for (int i = 0; i < 3; i++) _d(i) = d(i);
+    setDirection(_d, d);
   }
   virtual ~LagrangeMultiplierTerm() {}
   virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
   virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< fullMatrix<double> > &vm) const{};
-  virtual BilinearTermBase* clone () const {return new LagrangeMultiplierTerm(BilinearTerm<SVector3, double>::space1,BilinearTerm<SVector3, double>::space2,_d);}
+  virtual BilinearTermBase* clone () const {return new LagrangeMultiplierTerm(BilinearTerm<T1, double>::space1, BilinearTerm<T1, double>::space2, _d);}
 };
 
 class LagMultTerm : public BilinearTerm<SVector3, SVector3>
@@ -276,9 +282,9 @@ template<class T1> class LoadTermOnBorder : public LinearTerm<T1>
 {
 private :
   double _eqfac;
-  simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
+  simpleFunction<typename TensorialTraits<T1>::ValType> *Load;
 public :
-  LoadTermOnBorder(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_, double eqfac = 1.0) :
+  LoadTermOnBorder(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> *Load_, double eqfac = 1.0) :
       LinearTerm<T1>(space1_), _eqfac(eqfac), Load(Load_) {}
   virtual ~LoadTermOnBorder() {}
   virtual LinearTermBase<double>* clone () const { return new LoadTermOnBorder<T1>(LinearTerm<T1>::space1,Load,_eqfac);}

Solver/terms.hpp

@@ -7,7 +7,7 @@
 //   Eric Bechet
 //
 
-#include "terms.h"
+//#include "terms.h"
 
 template<class T2> void LinearTermBase<T2>::get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &vec) const
 {
@@ -101,7 +101,7 @@ template<class T1> void LoadTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fu
     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());
+    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;
@@ -191,25 +191,49 @@ template<class T1> void GradTerm<T1>::get(MElement *ele, int npts, IntPt *GP, st
   }
 }
 
-
+template<class T1> void LagrangeMultiplierTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+  int nbFF1 = BilinearTerm<T1, double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+  int nbFF2 = BilinearTerm<T1, 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<typename TensorialTraits<T1>::ValType> Vals;
+    std::vector<TensorialTraits<double>::ValType> ValsT;
+    BilinearTerm<T1,double>::space1.f(ele, u, v, w, Vals);
+    BilinearTerm<T1,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;
+      }
+    }
+  }
+}
 
 template<class T1> void LoadTermOnBorder<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const
 {
-  MElement *elep;
-  if (ele->getParent()) elep = 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];
+    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<typename TensorialTraits<T1>::ValType> Vals;
     LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
+    if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+       ele->getTypeForMSH() == MSH_POLYG_B)
+      ele->movePointFromParentSpaceToElementSpace(u, v, w);
     SPoint3 p;
     ele->pnt(u, v, w, p);
-    typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+    typename TensorialTraits<T1>::ValType load = (*Load)(p.x(), p.y(), p.z());
     for(int j = 0; j < nbFF ; ++j){
       m(j) += _eqfac * dot(Vals[j], load) * weight * detJ;
     }

Solver/thermicSolver.cpp

+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include <string.h>
+#include "GmshConfig.h"
+#include "thermicSolver.h"
+#include "linearSystemCSR.h"
+#include "linearSystemPETSc.h"
+#include "linearSystemGMM.h"
+#include "linearSystemFull.h"
+#include "Numeric.h"
+#include "GModel.h"
+#include "functionSpace.h"
+#include "terms.h"
+#include "solverAlgorithms.h"
+#include "quadratureRules.h"
+#include "solverField.h"
+#include "MPoint.h"
+#include "gmshLevelset.h"
+#if defined(HAVE_POST)
+#include "PView.h"
+#include "PViewData.h"
+#endif
+
+void thermicSolver::setMesh(const std::string &meshFileName)
+{
+  pModel = new GModel();
+  pModel->readMSH(meshFileName.c_str());
+  _dim = pModel->getNumRegions() ? 3 : 2;
+
+  if (LagSpace) delete LagSpace;
+  LagSpace = new ScalarLagrangeFunctionSpace(_tag);
+
+  if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
+  LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement(_tag + 1);
+}
+
+void thermicSolver::solve()
+{
+  linearSystemFull<double> *lsys = new linearSystemFull<double>;
+/*#if defined(HAVE_TAUCS)
+  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+  linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#else
+  linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+  lsys->setNoisy(2);
+#endif*/
+  assemble(lsys);
+  lsys->systemSolve();
+  printf("-- done solving!\n");
+}
+
+void thermicSolver::cutMesh(gLevelset *ls)
+{
+  pModel = pModel->buildCutGModel(ls);
+  pModel->writeMSH("cutMesh.msh");
+}
+
+void thermicSolver::setThermicDomain(int phys, double k)
+{
+  thermicField field;
+  field._k = k;
+  field._tag = _tag;
+  field.g = new groupOfElements (_dim, phys);
+  thermicFields.push_back(field);
+}
+
+void thermicSolver::setLagrangeMultipliers(int phys, double tau, int tag, simpleFunction<double> *f)
+{
+  LagrangeMultiplierFieldT field;
+  field._tau = tau;
+  field._tag = tag;
+  field._f = f;
+  field.g = new groupOfElements (_dim - 1, phys);
+  LagrangeMultiplierFields.push_back(field);
+}
+
+void thermicSolver::setEdgeTemp(int edge, simpleFunction<double> *f)
+{
+  dirichletBCT diri;
+  diri.g = new groupOfElements (1, edge);
+  diri._f = f;
+  diri._tag = edge;
+  diri.onWhat = BoundaryConditionT::ON_EDGE;
+  allDirichlet.push_back(diri);
+}
+
+void thermicSolver::assemble(linearSystem<double> *lsys)
+{
+  if (pAssembler) delete pAssembler;
+  pAssembler = new dofManager<double>(lsys);
+
+  // we first do all fixations. the behavior of the dofManager is to
+  // give priority to fixations : when a dof is fixed, it cannot be
+  // numbered afterwards
+
+  // Dirichlet conditions
+  for (unsigned int i = 0; i < allDirichlet.size(); i++){
+    FilterDofTrivial filter;
+    FixNodalDofs(*LagSpace, allDirichlet[i].g->begin(), allDirichlet[i].g->end(),
+                 *pAssembler, *allDirichlet[i]._f, filter);
+  }
+  // LagrangeMultipliers
+  for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
+    NumberDofs(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(),
+               LagrangeMultiplierFields[i].g->end(), *pAssembler);
+  }
+  // Thermic Fields
+  for (unsigned int i = 0; i < thermicFields.size(); ++i){
+    NumberDofs(*LagSpace, thermicFields[i].g->begin(), thermicFields[i].g->end(), *pAssembler);
+  }
+  // Neumann conditions
+  GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
+  for (unsigned int i = 0; i < allNeumann.size(); i++){
+    std::cout <<  "Neumann BC" << std::endl;
+    LoadTerm<double> 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++){
+    printf("Lagrange Mult Lag\n");
+    LagrangeMultiplierTerm<double> LagTerm(*LagSpace, *LagrangeMultiplierSpace, 1.);
+    Assemble(LagTerm, *LagSpace, *LagrangeMultiplierSpace,
+             LagrangeMultiplierFields[i].g->begin(), LagrangeMultiplierFields[i].g->end(),
+             Integ_LagrangeMult, *pAssembler);
+    printf("Lagrange Mult Lap\n");
+    LaplaceTerm<double, double> LapTerm(*LagrangeMultiplierSpace,
+                                        -LagrangeMultiplierFields[i]._tau);
+    Assemble(LapTerm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(),
+             LagrangeMultiplierFields[i].g->end(), Integ_Laplace, *pAssembler);
+    printf("Lagrange Mult Load\n");
+    LoadTermOnBorder<double> Lterm(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i]._f);
+    Assemble(Lterm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(),
+             LagrangeMultiplierFields[i].g->end(), Integ_Boundary, *pAssembler);
+  }
+  // Assemble thermic term
+  GaussQuadrature Integ_Bulk(GaussQuadrature::ValVal);
+  for (unsigned int i = 0; i < thermicFields.size(); i++){
+    printf("Thermic Term\n");
+    LaplaceTerm<double, double> Tterm(*LagSpace, thermicFields[i]._k);
+    Assemble(Tterm, *LagSpace, thermicFields[i].g->begin(), thermicFields[i].g->end(),
+             Integ_Bulk, *pAssembler);
+  }
+
+  /*for (int i = 0;i<pAssembler->sizeOfR();i++){
+    for (int j = 0;j<pAssembler->sizeOfR();j++){
+      double d; lsys->getFromMatrix(i, j, d);
+      printf("%g ", d);
+    }
+    double d; lsys->getFromRightHandSide(i, d);
+    printf(" |  %g\n", d);
+  }*/
+
+  printf("nDofs=%d\n", pAssembler->sizeOfR());
+  printf("nFixed=%d\n", pAssembler->sizeOfF());
+}
+
+double thermicSolver::computeL2Norm(simpleFunction<double> *sol) {
+  double val = 0.0;
+  SolverField<double> solField(pAssembler, LagSpace);
+  for (unsigned int i = 0; i < thermicFields.size(); ++i){
+    for (groupOfElements::elementContainer::const_iterator it = thermicFields[i].g->begin(); it != thermicFields[i].g->end(); ++it){
+      MElement *e = *it;
+//printf("element (%g,%g) (%g,%g) (%g,%g)\n",e->getVertex(0)->x(),e->getVertex(0)->y(),e->getVertex(1)->x(),e->getVertex(1)->y(),e->getVertex(2)->x(),e->getVertex(2)->y());
+      int npts;
+      IntPt *GP;
+      double jac[3][3];
+      int integrationOrder = 2 * (e->getPolynomialOrder()+5);
+      e->getIntegrationPoints(integrationOrder, &npts, &GP);
+      for (int j = 0; j < npts; j++){
+        double u = GP[j].pt[0];
+        double v = GP[j].pt[1];
+        double w = GP[j].pt[2];
+        double weight = GP[j].weight;
+        double detJ = fabs(e->getJacobian (u, v, w, jac));
+        SPoint3 p;
+        e->pnt(u, v, w, p);
+	double FEMVALUE;
+	solField.f(e, u, v, w, FEMVALUE);
+        double diff = (*sol)(p.x(), p.y(), p.z()) - FEMVALUE;
+        val += diff * diff * detJ * weight;
+//printf("(%g %g) : detJ=%g we=%g FV=%g sol=%g diff=%g\n",p.x(),p.y(),detJ,weight,FEMVALUE,(*sol)(p.x(), p.y(), p.z()),diff);
+      }
+    }
+  }
+printf("L2Norm = %g\n",sqrt(val));
+  return sqrt(val);
+}
+
+double thermicSolver::computeLagNorm(int tag, simpleFunction<double> *sol) {
+  double val = 0.0, val2 = 0.0;
+  SolverField<double> solField(pAssembler, LagrangeMultiplierSpace);
+  for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
+    if(tag != LagrangeMultiplierFields[i]._tag) continue;
+    for(groupOfElements::elementContainer::const_iterator it = LagrangeMultiplierFields[i].g->begin();
+        it != LagrangeMultiplierFields[i].g->end(); ++it){
+      MElement *e = *it;
+printf("element (%g,%g) (%g,%g)\n",e->getVertex(0)->x(),e->getVertex(0)->y(),e->getVertex(1)->x(),e->getVertex(1)->y());
+      int npts;
+      IntPt *GP;
+      double jac[3][3];
+      int integrationOrder = 2 * (e->getPolynomialOrder()+1);
+      e->getIntegrationPoints(integrationOrder, &npts, &GP);
+      for (int j = 0; j < npts; j++){
+        double u = GP[j].pt[0];
+        double v = GP[j].pt[1];
+        double w = GP[j].pt[2];
+        double weight = GP[j].weight;
+        double detJ = fabs(e->getJacobian (u, v, w, jac));
+        SPoint3 p;
+        e->getParent()->pnt(u, v, w, p);
+	double FEMVALUE;
+	solField.f(e, u, v, w, FEMVALUE);
+        double diff = (*sol)(p.x(), p.y(), p.z()) - FEMVALUE;
+        val += diff * diff * detJ * weight;
+        val2 += (*sol)(p.x(), p.y(), p.z()) * (*sol)(p.x(), p.y(), p.z()) * detJ * weight;
+printf("(%g %g) : u,v=(%g,%g) detJ=%g we=%g FV=%g sol=%g diff=%g\n",p.x(),p.y(),u,v,detJ,weight,FEMVALUE,(*sol)(p.x(), p.y(), p.z()),diff);
+      }
+    }
+  }
+printf("LagNorm = %g\n",sqrt(val/val2));
+  return sqrt(val/val2);
+}
+
+#if defined(HAVE_POST)
+
+PView* thermicSolver::buildTemperatureView (const std::string postFileName)
+{
+  std::cout <<  "build Temperature View"<< std::endl;
+  std::set<MVertex*> v;
+  std::map<MVertex*, MElement*> vCut;
+  for (unsigned int i = 0; i < thermicFields.size(); ++i){
+    for (groupOfElements::elementContainer::const_iterator it = thermicFields[i].g->begin(); it != thermicFields[i].g->end(); ++it){
+      MElement *e = *it;
+      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<double> Field(pAssembler, LagSpace);
+  for (std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it){
+    double val;
+    MPoint p(*it);
+    Field.f(&p, 0, 0, 0, val); //printf("valv=%g\n",val);
+    std::vector<double> vec; vec.push_back(val);
+    data[(*it)->getNum()] = vec;
+  }
+  for (std::map<MVertex*, MElement*>::iterator it = vCut.begin(); it != vCut.end(); ++it){
+    double 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); //printf("valvc=%g\n",val);
+    std::vector<double> vec; vec.push_back(val);
+    data[it->first->getNum()] = vec;
+  }
+  PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0, 1);
+  return pv;
+}
+
+PView* thermicSolver::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, 1);
+  return pv;
+}
+
+#else
+PView* thermicSolver::buildTemperatureView(const std::string postFileName)
+{
+  Msg::Error("Post-pro module not available");
+  return 0;
+}
+
+PView* thermicSolver::buildLagrangeMultiplierView (const std::string postFileName)
+{
+  Msg::Error("Post-pro module not available");
+  return 0;
+}
+#endif
+
+
+
+
+

Solver/thermicSolver.h

+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _THERMIC_SOLVER_H_
+#define _THERMIC_SOLVER_H_
+
+#include <map>
+#include <string>
+#include "SVector3.h"
+#include "dofManager.h"
+#include "simpleFunction.h"
+#include "functionSpace.h"
+
+class GModel;
+class PView;
+class groupOfElements;
+class gLevelset;
+
+struct LagrangeMultiplierFieldT {
+  int _tag;
+  groupOfElements *g;
+  double _tau;
+  simpleFunction<double> *_f;
+  LagrangeMultiplierFieldT() : _tag(0), g(0){}
+};
+
+struct thermicField {
+  int _tag; // tag for the dofManager
+  groupOfElements *g; // support for this field
+  double _k; // diffusivity
+  thermicField () : _tag(0), g(0){}
+};
+
+struct BoundaryConditionT
+{
+  int _tag; // tag for the dofManager
+  enum location{UNDEF, ON_VERTEX, ON_EDGE, ON_FACE, ON_VOLUME};
+  location onWhat; // on vertices or elements
+  groupOfElements *g; // support for this BC
+  BoundaryConditionT() : _tag(0), onWhat(UNDEF), g(0) {}
+};
+
+struct dirichletBCT : public BoundaryConditionT
+{
+  simpleFunction<double> *_f;
+  dirichletBCT ():BoundaryConditionT(), _f(0){}
+};
+
+struct neumannBCT  : public BoundaryConditionT
+{
+  simpleFunction<double> *_f;
+  neumannBCT () : BoundaryConditionT(), _f(0){}
+};
+// a thermic solver ...
+class thermicSolver
+{
+ protected:
+  GModel *pModel;
+  int _dim, _tag;
+  dofManager<double> *pAssembler;
+  FunctionSpace<double> *LagSpace;
+  FunctionSpace<double> *LagrangeMultiplierSpace;
+
+  // young modulus and poisson coefficient per physical
+  std::vector<thermicField> thermicFields;
+  std::vector<LagrangeMultiplierFieldT> LagrangeMultiplierFields;
+  // neumann BC
+  std::vector<neumannBCT> allNeumann;
+  // dirichlet BC
+  std::vector<dirichletBCT> allDirichlet;
+
+ public:
+  thermicSolver(int tag) : _tag(tag), pAssembler(0), LagSpace(0), LagrangeMultiplierSpace(0) {}
+
+  virtual ~thermicSolver()
+  {
+    if (LagSpace) delete LagSpace;
+    if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
+    if (pAssembler) delete pAssembler;
+  }
+  void assemble (linearSystem<double> *lsys);
+  virtual void setMesh(const std::string &meshFileName);
+  void cutMesh(gLevelset *ls);
+  void setThermicDomain(int phys, double k);
+  void setLagrangeMultipliers(int phys, double tau, int tag, simpleFunction<double> *f);
+  void setEdgeTemp(int edge, simpleFunction<double> *f);
+  void solve();
+  virtual PView *buildTemperatureView(const std::string postFileName);
+  virtual PView *buildLagrangeMultiplierView(const std::string posFileName);
+  double computeL2Norm(simpleFunction<double> *f);
+  double computeLagNorm(int tag, simpleFunction<double> *f);
+  // std::pair<PView *, PView*> buildErrorEstimateView
+  //   (const std::string &errorFileName, double, int);
+  // std::pair<PView *, PView*> buildErrorEstimateView
+  //   (const std::string &errorFileName, const elasticityData &ref, double, int);
+};
+
+#endif

wrappers/gmshpy/gmshPost.i

@@ -16,6 +16,7 @@
   #include "PViewFactory.h"
   #include "PViewData.h"
   #include "PViewAsSimpleFunction.h"
+  #include "PViewDataGModel.h"
 #endif
 %}
 
@@ -34,6 +35,7 @@ namespace std {
 %include "simpleFunction.h"
 %template(simpleFunctionDouble) simpleFunction<double>;
 %include "PViewAsSimpleFunction.h"
+%include "PViewDataGModel.h"
 %include "Plugin.h"
 %include "PluginManager.h"
 #endif

wrappers/gmshpy/gmshSolver.i

@@ -10,6 +10,7 @@
 #if defined(HAVE_SOLVER)
   #include "dofManager.h"
   #include "elasticitySolver.h"
+  #include "thermicSolver.h"
   #include "frameSolver.h"
   #include "linearSystem.h"
   #include "linearSystemCSR.h"
@@ -24,6 +25,7 @@
 %include "dofManager.h"
 %template(dofManagerDouble) dofManager<double>;
 %include "elasticitySolver.h"
+%include "thermicSolver.h"
 %include "frameSolver.h"
 %include "linearSystem.h"
 %template(linearSystemDouble) linearSystem<double>;