new projection finite element formulation (jf pirenne)

Graphics/Makefile

@@ -77,7 +77,7 @@ drawContext${OBJEXT}: drawContext.cpp ../Common/GmshMessage.h drawContext.h \
   ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
   ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Post/PView.h ../Post/PViewOptions.h \
-  ../Post/ColorTable.h
+  ../Post/ColorTable.h gl2ps.h
 drawMesh${OBJEXT}: drawMesh.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Common/GmshMessage.h ../Common/GmshDefines.h \
   ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \

Mesh/Makefile

@@ -259,7 +259,8 @@ meshGFaceOptimize${OBJEXT}: meshGFaceOptimize.cpp meshGFaceOptimize.h \
   ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
-  ../Geo/Pair.h BackgroundMesh.h Generator.h
+  ../Geo/Pair.h BackgroundMesh.h ../Numeric/Numeric.h \
+  ../Numeric/GmshMatrix.h Generator.h
 meshGFaceQuadrilateralize${OBJEXT}: meshGFaceQuadrilateralize.cpp \
   meshGFaceQuadrilateralize.h ../Common/GmshMessage.h \
   ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \

Mesh/meshGFaceOptimize.cpp

@@ -244,8 +244,6 @@ void laplaceSmoothing(GFace *gf)
   }
 }
 
-extern void fourthPoint(double *p1, double *p2, double *p3, double *p4);
-
 double surfaceTriangleUV(MVertex *v1, MVertex *v2, MVertex *v3,           
                          const std::vector<double> &Us,
                          const std::vector<double> &Vs)

Numeric/Makefile

@@ -16,6 +16,7 @@ SRC = Numeric.cpp\
       gmshLaplace.cpp\
         gmshHelmholtz.cpp\
         gmshElasticity.cpp\
+        gmshProjection.cpp\
       EigSolve.cpp\
       FunctionSpace.cpp\
       GaussQuadratureLin.cpp\
@@ -103,6 +104,22 @@ gmshElasticity${OBJEXT}: gmshElasticity.cpp gmshElasticity.h \
   ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
   ../Common/Gmsh.h ../Common/GmshMessage.h Numeric.h
+gmshProjection${OBJEXT}: gmshProjection.cpp gmshProjection.h \
+  gmshTermOfFormulation.h GmshMatrix.h ../Common/GmshConfig.h \
+  ../Common/GmshMessage.h gmshFunction.h gmshAssembler.h \
+  gmshLinearSystem.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h \
+  ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
+  ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
+  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h \
+  ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h \
+  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/MElement.h \
+  ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+  ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
+  ../Common/Gmsh.h ../Common/GmshMessage.h Numeric.h
 EigSolve${OBJEXT}: EigSolve.cpp
 FunctionSpace${OBJEXT}: FunctionSpace.cpp FunctionSpace.h GmshMatrix.h \
   ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/GmshDefines.h

Numeric/gmshHelmholtz.cpp

@@ -10,7 +10,7 @@ void gmshHelmholtzTerm::elementMatrix(MElement *e,
                                       gmshMatrix<std::complex<double> > &m) const
 {
   int nbNodes = e->getNumVertices();
-  int integrationOrder = 2 * e->getPolynomialOrder() - 0;
+  int integrationOrder = 2 * e->getPolynomialOrder();
   int npts;
   IntPt *GP;
   double jac[3][3];
@@ -44,7 +44,7 @@ void gmshHelmholtzTerm::elementMatrix(MElement *e,
 	m(j, k) += ((Grads[j][0] * Grads[k][0] +
                      Grads[j][1] * Grads[k][1] +
                      Grads[j][2] * Grads[k][2]) 
-                    - kk * kk * sf[j] * sf[k]) * weight * detJ * 0.5;
+                    - kk * kk * sf[j] * sf[k]) * weight * detJ;
       }
     }
   }

Numeric/gmshProjection.cpp

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "gmshProjection.h"
+#include "Numeric.h"
+
+void gmshProjectionTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+{
+  int nbNodes = e->getNumVertices();
+  int integrationOrder = 2 * e->getPolynomialOrder();
+  int npts;
+  IntPt *GP;
+  double jac[3][3];
+  double sf[256];
+  e->getIntegrationPoints(integrationOrder, &npts, &GP);
+  m.set_all(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);   
+    SPoint3 p; e->pnt(u, v, w, p);
+    e->getShapeFunctions(u, v, w, sf);
+    for (int j = 0; j < nbNodes; j++){
+      for (int k = 0; k < nbNodes; k++){
+	m(j, k) += sf[j] * sf[k] * weight * detJ;
+      }
+    }
+  }
+} 

Numeric/gmshProjection.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _GMSH_PROJECTION_H_
+#define _GMSH_PROJECTION_H_
+
+#include "gmshTermOfFormulation.h"
+#include "gmshFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "GmshMatrix.h"
+
+class gmshProjectionTerm : public gmshNodalFemTerm<double> {
+ private:
+  const int _iField ;
+ protected:
+  virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
+  virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
+  void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
+  {
+    *vR = e->getVertex(iRow);
+    *iCompR = 0;
+    *iFieldR = _iField;
+  }
+ public:
+  gmshProjectionTerm(GModel *gm, int iField = 0) : 
+    gmshNodalFemTerm<double>(gm), _iField(iField){}
+  virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+};
+
+#endif

Numeric/gmshTermOfFormulation.h

@@ -124,6 +124,7 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
     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];
@@ -140,17 +141,9 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
 	  const double w = GP[ip].pt[2];
 	  const double weight = GP[ip].weight;
 	  const double detJ = e->getJacobian(u, v, w, jac);
-	  double x = 0;
-	  double y = 0;
-	  double z = 0;
-          double sf[nbNodes];
+          SPoint3 p; e->pnt(u, v, w, p);
           e->getShapeFunctions(u, v, w, sf);
-	  for (int k = 0; k < nbNodes; k++){  
-	    x += e->getVertex(k)->x() * sf[k];
-	    y += e->getVertex(k)->y() * sf[k];
-	    z += e->getVertex(k)->z() * sf[k];
-	  }
-	  const scalar FCT = fct (x,y,z);
+	  const scalar FCT = fct(p.x(), p.y(), p.z());
 	  for (int k = 0; k < nbNodes; k++){
 	    lsys.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
 	  }

Plugin/FiniteElement.cpp

@@ -3,27 +3,33 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <stdlib.h>
 #include "GmshConfig.h"
 #include "GModel.h"
 #include "FiniteElement.h"
 #include "gmshAssembler.h"
+#include "gmshProjection.h"
 #include "gmshLaplace.h"
 #include "gmshHelmholtz.h"
 #include "gmshLinearSystemGmm.h"
 
+#if defined(HAVE_MATH_EVAL)
+#include "matheval.h"
+#endif
+
 StringXNumber FiniteElementOptions_Number[] = {
-  {GMSH_FULLRC, "Parameter", NULL, 1.},
-  {GMSH_FULLRC, "Volume", NULL, 1.},
-  {GMSH_FULLRC, "Surface1", NULL, 2.},
-  {GMSH_FULLRC, "Value1", NULL, 0.},
-  {GMSH_FULLRC, "Surface2", NULL, 3.},
-  {GMSH_FULLRC, "Value2", NULL, 1.}
+  {GMSH_FULLRC, "Omega", NULL, 1.},
+  {GMSH_FULLRC, "Gamma1", NULL, 0.},
+  {GMSH_FULLRC, "Gamma1Value", NULL, 0.},
+  {GMSH_FULLRC, "Gamma2", NULL, 0.},
+  {GMSH_FULLRC, "Gamma2Value", NULL, 0.}
 };
 
 StringXString FiniteElementOptions_String[] = {
-  {GMSH_FULLRC, "Equation", NULL, "Laplace"},
-  {GMSH_FULLRC, "BC1", NULL, "Dirichlet"},
-  {GMSH_FULLRC, "BC2", NULL, "Dirichlet"},
+  {GMSH_FULLRC, "Equation", NULL, "Projection"},
+  {GMSH_FULLRC, "Parameter", NULL, "Sin(x*y)"},
+  {GMSH_FULLRC, "Gamma1BC", NULL, ""},
+  {GMSH_FULLRC, "Gamma2BC", NULL, ""},
 };
 
 extern "C"
@@ -78,6 +84,40 @@ void GMSH_FiniteElementPlugin::catchErrorMessage(char *errorMessage) const
   strcpy(errorMessage, "FiniteElement failed...");
 }
 
+template<class scalar>
+class gmshMathEvalFunction : public gmshFunction<scalar> {
+ private:
+  std::string _str;
+  void *_f;
+ public:
+  gmshMathEvalFunction(std::string str) : _str(str)
+  {
+#if defined(HAVE_MATH_EVAL)
+    _f = evaluator_create((char*)_str.c_str());
+    if(!_f) Msg::Error("Invalid expression '%s'", _str.c_str());
+#endif
+  }
+  virtual ~gmshMathEvalFunction()
+  {
+#if defined(HAVE_MATH_EVAL)
+    if(_f) evaluator_destroy(_f);
+#endif
+  }
+  virtual scalar operator () (double x, double y, double z) const 
+  { 
+#if defined(HAVE_MATH_EVAL)
+    if(_f){
+      char *names[] = {"x", "y", "z"};
+      double values[] = {x, y, z};
+      return evaluator_evaluate(_f, sizeof(names)/sizeof(names[0]), names, values);
+    }
+    return atof(_str.c_str());
+#else
+    return atof(_str.c_str());
+#endif
+  }
+};
+
 #if defined(HAVE_GMM)
 template<class scalar>
 class solver{
@@ -86,13 +126,13 @@ class solver{
   gmshAssembler<scalar> *myAssembler;
   solver()
   {
-    int volume = (int)FiniteElementOptions_Number[1].def;
-    int surface1 = (int)FiniteElementOptions_Number[2].def;
-    double value1 = FiniteElementOptions_Number[3].def;
-    int surface2 = (int)FiniteElementOptions_Number[4].def;
-    double value2 = FiniteElementOptions_Number[5].def;
-    std::string bc1 = FiniteElementOptions_String[1].def;
-    std::string bc2 = FiniteElementOptions_String[2].def;
+    int volume = (int)FiniteElementOptions_Number[0].def;
+    int surface1 = (int)FiniteElementOptions_Number[1].def;
+    double value1 = FiniteElementOptions_Number[2].def;
+    int surface2 = (int)FiniteElementOptions_Number[3].def;
+    double value2 = FiniteElementOptions_Number[4].def;
+    std::string bc1 = FiniteElementOptions_String[2].def;
+    std::string bc2 = FiniteElementOptions_String[3].def;
 
     lsys = new gmshLinearSystemGmm<scalar>;
     lsys->setPrec(1.e-10);
@@ -127,8 +167,8 @@ class solver{
 PView *GMSH_FiniteElementPlugin::execute(PView *v)
 {
   std::string equation = FiniteElementOptions_String[0].def;
-  double parameter = FiniteElementOptions_Number[0].def;
-  int volume = (int)FiniteElementOptions_Number[1].def;
+  std::string parameter = FiniteElementOptions_String[1].def;
+  int volume = (int)FiniteElementOptions_Number[0].def;
 
 #if defined(HAVE_GMM)
   GModel *m = GModel::current();
@@ -140,11 +180,25 @@ PView *GMSH_FiniteElementPlugin::execute(PView *v)
   m->getMeshVertices(volume, dim, vertices);
   std::map<int, std::vector<double> > data;
 
-  if(equation == "Laplace"){
+  if(equation == "Projection"){
+    solver<double> s;
+    if(!s.myAssembler->sizeOfR()) return 0;
+    gmshMathEvalFunction<double> par(parameter);
+    gmshProjectionTerm projection(m, 1);
+    for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
+      projection.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
+    projection.addNeumann(volume, dim, 0, 1, par, *s.myAssembler);
+    s.lsys->systemSolve();
+    for (unsigned int i = 0; i < vertices.size(); i++)
+      data[vertices[i]->getNum()].push_back
+        (s.myAssembler->getDofValue(vertices[i], 0, 1));
+    PView *pv = new PView("projection", "NodeData", m, data, 0.);
+  }
+  else if(equation == "Laplace"){
     solver<double> s;
     if(!s.myAssembler->sizeOfR()) return 0;
-    gmshFunction<double> diffusivity(parameter);
-    gmshLaplaceTerm laplace(m, &diffusivity, 1);
+    gmshMathEvalFunction<double> par(parameter);
+    gmshLaplaceTerm laplace(m, &par, 1);
     for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
       laplace.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
     s.lsys->systemSolve();
@@ -156,8 +210,8 @@ PView *GMSH_FiniteElementPlugin::execute(PView *v)
   else if(equation == "Helmholtz"){
     solver<std::complex<double> > s;
     if(!s.myAssembler->sizeOfR()) return 0;
-    gmshFunction<std::complex<double> > waveNumber(parameter);
-    gmshHelmholtzTerm helmholtz(m, &waveNumber, 1);
+    gmshMathEvalFunction<std::complex<double> > par(parameter);
+    gmshHelmholtzTerm helmholtz(m, &par, 1);
     for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
       helmholtz.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
     s.lsys->setGmres(1);

Plugin/Makefile

@@ -320,7 +320,7 @@ FiniteElement${OBJEXT}: FiniteElement.cpp ../Common/GmshConfig.h ../Geo/GModel.h
   ../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
   ../Post/PViewDataList.h ../Post/PViewData.h ../Numeric/GmshMatrix.h \
   ../Common/ListUtils.h ../Numeric/gmshAssembler.h \
-  ../Numeric/gmshLinearSystem.h ../Numeric/gmshLaplace.h \
+  ../Numeric/gmshLinearSystem.h ../Numeric/gmshProjection.h \
   ../Numeric/gmshTermOfFormulation.h ../Numeric/GmshMatrix.h \
   ../Numeric/gmshFunction.h ../Numeric/gmshAssembler.h ../Geo/MElement.h \
   ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
@@ -328,6 +328,8 @@ FiniteElement${OBJEXT}: FiniteElement.cpp ../Common/GmshConfig.h ../Geo/GModel.h
   ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
   ../Numeric/gmshFunction.h ../Common/Gmsh.h ../Common/GmshMessage.h \
+  ../Numeric/GmshMatrix.h ../Numeric/gmshLaplace.h \
+  ../Numeric/gmshTermOfFormulation.h ../Numeric/gmshFunction.h \
   ../Numeric/GmshMatrix.h ../Numeric/gmshHelmholtz.h \
   ../Numeric/gmshTermOfFormulation.h ../Numeric/gmshFunction.h \
   ../Numeric/GmshMatrix.h ../Numeric/gmshLinearSystemGmm.h \

doc/TODO.txt

-$Id: TODO.txt,v 1.20 2009-03-10 12:06:14 geuzaine Exp $
+$Id: TODO.txt,v 1.21 2009-03-17 07:36:33 geuzaine Exp $
+
+********************************************************************
+
+extrude along a given curve (use curve->closestPoint for each xyz in
+ExtrudeParams and extrude along the tangent? not that simple: this
+could lead to self-intersection if curvature >>)
 
 ********************************************************************
 
@@ -129,10 +135,6 @@ implement better algo to determine which axes to draw
 
 ********************************************************************
 
-extrude along a given function or along a parametric curve
-
-********************************************************************
-
 Custom range on "filled iso" 3D views produces ugly plots, where one
 sees inside the cut elements. (In the meantime, one can use
 Plugin(CutMap) with "ExtractVolume" set to 1 (or -1): this will

doc/VERSIONS.txt

-$Id: VERSIONS.txt,v 1.40 2009-03-16 07:23:48 geuzaine Exp $
+$Id: VERSIONS.txt,v 1.41 2009-03-17 07:36:33 geuzaine Exp $
 
-(?): removed GSL dependency (Gmsh now simply requires Blas and
+2.3.1 (?): removed GSL dependency (Gmsh now simply uses Blas and
 Lapack); new per-window visibility; added support for composite window
-printing and background images; fixes for string options in parser;
-fixed surface mesh orientation for Open Cascade models; fixed triangle
-orientations in Delaunay and Frontal algorithms.
+printing and background images; fixed string option affectation in
+parser; fixed surface mesh orientation for Open Cascade models; fixed
+random triangle orientations in Delaunay and Frontal algorithms.
 
 2.3.0 (Jan 23, 2009): major graphics and GUI code refactoring; new
 full-quad/hexa subdivision algorithm (removed Mesh.RecombineAlgo);