*** empty log message ***

Makefile

@@ -20,7 +20,7 @@ GMSH_DATE = `date "+%Y%m%d"`
 GMSH_API = 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/GFaceCompound.h\
-           Geo/MVertex.h Geo/MEdge.h Geo/MFace.h Geo/MElement.h\
+           Geo/MVertex.h Geo/MEdge.h Geo/MFace.h Geo/MElement.h Geo/MTriangle.h Geo/MQuadrangle.h Geo/MTetrahedron.h Geo/MHexahedron.h Geo/MPyramid.h Geo/MPrism.h Geo/MLine.h \
            Geo/discreteVertex.h Geo/discreteEdge.h Geo/discreteFace.h Geo/discreteRegion.h\
            Geo/SPoint2.h Geo/SPoint3.h Geo/SVector3.h Geo/SBoundingBox3d.h\
            Geo/Pair.h Geo/Range.h\

Numeric/gmshElasticity.cpp

@@ -74,3 +74,29 @@ void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
     m.gemm(BTH, B, weight * detJ, 1.);
   } 
 }
+
+void gmshElasticityTerm::elementVector(MElement *e, gmshVector<double> &m) const{
+  int nbNodes = e->getNumVertices();
+  int integrationOrder = 2 * e->getPolynomialOrder();
+  int npts;
+  IntPt *GP;
+  double jac[3][3];
+  double ff[256];
+  e->getIntegrationPoints(integrationOrder, &npts, &GP);
+  
+  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 = e->getJacobian(u, v, w, jac);
+    e->getShapeFunctions(u, v, w, ff);
+    for (int j = 0; j < nbNodes; j++){
+      m(j)           += _f.x() *ff[j] * weight * detJ *.5;
+      m(j+nbNodes)   += _f.y() *ff[j] * weight * detJ *.5;
+      m(j+2*nbNodes) += _f.z() *ff[j] * weight * detJ *.5;
+    }
+  } 
+}

Numeric/gmshElasticity.h

@@ -15,6 +15,7 @@
 class gmshElasticityTerm : public gmshNodalFemTerm<double> {
   double _E, _nu;
   int _iField;
+  SVector3 _f;
  protected:
   virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
   virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
@@ -28,7 +29,9 @@ class gmshElasticityTerm : public gmshNodalFemTerm<double> {
  public:
  gmshElasticityTerm(GModel *gm, double E, double nu, int iField = 1) : 
   gmshNodalFemTerm<double>(gm), _E(E), _nu(nu), _iField(iField){}
+  void setVector(const SVector3 &f) {_f = f;}
   void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+  void elementVector(MElement *e, gmshVector<double> &m) const;
 };
 
 #endif

Numeric/gmshTermOfFormulation.h

@@ -47,6 +47,9 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
   gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation<scalar>(gm) {}
   virtual ~gmshNodalFemTerm (){}
   virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const = 0;
+  virtual void elementVector(MElement *e, gmshVector<scalar> &m) const {
+    m.scale(0.0);
+  }
   void addToMatrix(gmshAssembler<scalar> &lsys) const
   {
     GModel *m = gmshTermOfFormulation<scalar>::_gm;
@@ -151,6 +154,34 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
       }
     }
   }
+  void addToRightHandSide (gmshAssembler<scalar> &lsys) const {
+    GModel *m = gmshTermOfFormulation<scalar>::_gm;
+    if (m->getNumRegions()){
+      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
+	addToRightHandSide(lsys,*it);
+      }
+    }
+    else if(m->getNumFaces()){
+      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+	addToRightHandSide(lsys,*it);
+      }
+    }  
+  }    
+  void addToRightHandSide (gmshAssembler<scalar> &lsys, GEntity *ge) const {
+    for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+      MElement *e = ge->getMeshElement (i);
+      int nbR = sizeOfR(e);
+      gmshVector<scalar> V (nbR);
+      elementVector (e, V);
+      // assembly
+      for (int j=0;j<nbR;j++){
+	MVertex *vR;int iCompR,iFieldR;
+	getLocalDofR (e,j,&vR,&iCompR,&iFieldR);
+	lsys.assemble(vR,iCompR,iFieldR,V(j));
+      }
+    }
+  }
+
 };
 
 #endif

Post/PViewDataGModelIO.cpp

@@ -120,7 +120,6 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
   GModel *model = _steps[0]->getModel();
 
   binary = true;
-
   if(!model->writeMSH(fileName, 2.0, binary, true)) return false;
 
   // append data