From fffbd52f72f0b5dac45d722399eb33152fd0df1e Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 8 Feb 2009 19:22:47 +0000
Subject: [PATCH] added simple Helhmoltz solver + demo finite element plugin
---
Geo/GFaceCompound.cpp | 4 +-
Geo/Makefile | 17 ++--
Mesh/Makefile | 14 +--
Mesh/gmshSmoothHighOrder.cpp | 2 +-
Numeric/Makefile | 40 ++++----
Numeric/gmshElasticity.h | 2 +-
Numeric/gmshFunction.h | 7 +-
Numeric/gmshHelmholtz.cpp | 52 +++++++++++
Numeric/gmshHelmholtz.h | 37 ++++++++
Numeric/gmshLaplace.h | 4 +-
Plugin/FiniteElement.cpp | 174 +++++++++++++++++++++++++++++++++++
Plugin/FiniteElement.h | 30 ++++++
Plugin/Makefile | 31 ++++++-
Plugin/PluginManager.cpp | 3 +
14 files changed, 370 insertions(+), 47 deletions(-)
create mode 100644 Numeric/gmshHelmholtz.cpp
create mode 100644 Numeric/gmshHelmholtz.h
create mode 100644 Plugin/FiniteElement.cpp
create mode 100644 Plugin/FiniteElement.h
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 4269251d7e..5869bd998d 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -12,7 +12,7 @@
#include "gmshLinearSystemGmm.h"
#include "gmshLinearSystemFull.h"
-class gmshGradientBasedDiffusivity : public gmshFunction
+class gmshGradientBasedDiffusivity : public gmshFunction<double>
{
private:
MElement *_current;
@@ -274,7 +274,7 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
Msg::Debug("Creating term %d dofs numbered %d fixed",
myAssembler.sizeOfR(), myAssembler.sizeOfF());
- gmshLaplaceTerm laplace (model(), &diffusivity, 1);
+ gmshLaplaceTerm laplace(model(), &diffusivity, 1);
it = _compound.begin();
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
diff --git a/Geo/Makefile b/Geo/Makefile
index ffdc7096fe..241c249759 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -99,15 +99,16 @@ GFaceCompound${OBJEXT}: GFaceCompound.cpp ../Common/GmshConfig.h GFaceCompound.h
../Numeric/gmshAssembler.h ../Numeric/gmshLinearSystem.h \
../Numeric/gmshLaplace.h ../Numeric/gmshTermOfFormulation.h \
../Numeric/GmshMatrix.h ../Common/GmshMessage.h \
+ ../Numeric/gmshAssembler.h ../Geo/GModel.h ../Geo/GVertex.h \
+ ../Geo/GEdge.h ../Geo/GFace.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 \
../Numeric/gmshFunction.h ../Common/Gmsh.h ../Common/GmshMessage.h \
- ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.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 ../Numeric/GmshMatrix.h ../Numeric/Numeric.h \
- ../Numeric/GmshMatrix.h ../Common/Octree.h ../Common/OctreeInternals.h \
+ ../Numeric/GmshMatrix.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+ ../Common/Octree.h ../Common/OctreeInternals.h \
../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
../Numeric/gmshLinearSystemFull.h ../Numeric/gmshLinearSystem.h \
../Numeric/GmshMatrix.h
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 369729d83f..5da0e80b3f 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -129,13 +129,13 @@ gmshSmoothHighOrder${OBJEXT}: gmshSmoothHighOrder.cpp HighOrder.h \
meshGFaceDelaunayInsertion.h gmshSmoothHighOrder.h \
../Numeric/gmshAssembler.h ../Numeric/gmshLinearSystem.h \
../Numeric/gmshLaplace.h ../Numeric/gmshTermOfFormulation.h \
- ../Numeric/GmshMatrix.h ../Numeric/gmshFunction.h ../Common/Gmsh.h \
- ../Common/GmshMessage.h ../Numeric/GmshMatrix.h \
- ../Numeric/gmshElasticity.h ../Numeric/gmshTermOfFormulation.h \
- ../Numeric/GmshMatrix.h ../Numeric/gmshLinearSystemGmm.h \
- ../Numeric/gmshLinearSystem.h ../Common/Context.h ../Geo/CGNSOptions.h \
- ../Mesh/meshPartitionOptions.h ../Numeric/Numeric.h \
- ../Numeric/GmshMatrix.h
+ ../Numeric/GmshMatrix.h ../Numeric/gmshAssembler.h \
+ ../Numeric/gmshFunction.h ../Common/Gmsh.h ../Common/GmshMessage.h \
+ ../Numeric/GmshMatrix.h ../Numeric/gmshElasticity.h \
+ ../Numeric/gmshTermOfFormulation.h ../Numeric/GmshMatrix.h \
+ ../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
+ ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
+ ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
meshGEdge${OBJEXT}: meshGEdge.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index a32b9d2795..ef5e22824f 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -435,7 +435,7 @@ void localHarmonicMapping(GModel *gm,
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
gmshAssembler<double> myAssembler(lsys);
- gmshFunction f(1.0);
+ gmshFunction<double> f(1.0);
gmshLaplaceTerm Laplace(gm, &f, 0);
myAssembler.fixVertex ( n1 , 0 , 0 , -1.0);
diff --git a/Numeric/Makefile b/Numeric/Makefile
index 0ead3c621f..9ab2bb5e4d 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -14,6 +14,7 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
SRC = Numeric.cpp\
GmshMatrix.cpp\
gmshLaplace.cpp\
+ gmshHelmholtz.cpp\
gmshElasticity.cpp\
EigSolve.cpp\
FunctionSpace.cpp\
@@ -55,10 +56,9 @@ Numeric${OBJEXT}: Numeric.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
Numeric.h GmshMatrix.h
GmshMatrix${OBJEXT}: GmshMatrix.cpp ../Common/GmshConfig.h GmshMatrix.h \
../Common/GmshMessage.h
-gmshAssembler${OBJEXT}: gmshAssembler.cpp ../Geo/MVertex.h ../Geo/SPoint2.h \
- ../Geo/SPoint3.h gmshAssembler.h gmshLinearSystem.h
-gmshTermOfFormulation${OBJEXT}: gmshTermOfFormulation.cpp ../Common/Gmsh.h \
- ../Common/GmshMessage.h ../Geo/GModel.h ../Geo/GVertex.h \
+gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
+ GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.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 \
@@ -70,27 +70,27 @@ gmshTermOfFormulation${OBJEXT}: gmshTermOfFormulation.cpp ../Common/Gmsh.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 \
- ../Common/GmshConfig.h ../Numeric/Gauss.h GmshMatrix.h \
- gmshTermOfFormulation.h gmshFunction.h gmshAssembler.h \
- gmshLinearSystem.h
-gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
- GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
- gmshFunction.h ../Common/Gmsh.h ../Common/GmshMessage.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 \
+ ../Numeric/Gauss.h gmshFunction.h ../Common/Gmsh.h \
+ ../Common/GmshMessage.h Numeric.h
+gmshHelmholtz${OBJEXT}: gmshHelmholtz.cpp gmshHelmholtz.h \
+ gmshTermOfFormulation.h GmshMatrix.h ../Common/GmshConfig.h \
+ ../Common/GmshMessage.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 Numeric.h
+ ../Numeric/Gauss.h gmshFunction.h ../Common/Gmsh.h \
+ ../Common/GmshMessage.h Numeric.h
gmshElasticity${OBJEXT}: gmshElasticity.cpp gmshElasticity.h \
gmshTermOfFormulation.h GmshMatrix.h ../Common/GmshConfig.h \
- ../Common/GmshMessage.h ../Common/Gmsh.h ../Common/GmshMessage.h \
+ ../Common/GmshMessage.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 \
@@ -102,7 +102,7 @@ gmshElasticity${OBJEXT}: gmshElasticity.cpp gmshElasticity.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 Numeric.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
diff --git a/Numeric/gmshElasticity.h b/Numeric/gmshElasticity.h
index b0089f2f23..348545483a 100644
--- a/Numeric/gmshElasticity.h
+++ b/Numeric/gmshElasticity.h
@@ -13,7 +13,7 @@
#include "GmshMatrix.h"
class gmshElasticityTerm : public gmshNodalFemTerm<double> {
- double _E,_nu;
+ double _E, _nu;
int _iField;
protected:
virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
diff --git a/Numeric/gmshFunction.h b/Numeric/gmshFunction.h
index ed11c8876a..04478989ed 100644
--- a/Numeric/gmshFunction.h
+++ b/Numeric/gmshFunction.h
@@ -6,12 +6,13 @@
#ifndef _GMSH_FUNCTION_H_
#define _GMSH_FUNCTION_H_
+template <class scalar>
class gmshFunction {
- double _val;
+ scalar _val;
public :
- gmshFunction(double val = 0) : _val(val) {}
+ gmshFunction(scalar val=0) : _val(val) {}
virtual ~gmshFunction(){}
- virtual double operator () (double x, double y, double z) const { return _val; }
+ virtual scalar operator () (double x, double y, double z) const { return _val; }
};
#endif
diff --git a/Numeric/gmshHelmholtz.cpp b/Numeric/gmshHelmholtz.cpp
new file mode 100644
index 0000000000..132fa0648c
--- /dev/null
+++ b/Numeric/gmshHelmholtz.cpp
@@ -0,0 +1,52 @@
+// 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 "gmshHelmholtz.h"
+#include "Numeric.h"
+
+void gmshHelmholtzTerm::elementMatrix(MElement *e,
+ gmshMatrix<std::complex<double> > &m) const
+{
+ int nbNodes = e->getNumVertices();
+ int integrationOrder = 2 * e->getPolynomialOrder() - 0;
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ double invjac[3][3];
+ double sf[256], Grads[256][3], grads[256][3];
+ 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);
+ const std::complex<double> kk = (*_waveNumber)(p.x(), p.y(), p.z());
+ inv3x3(jac, invjac) ;
+ e->getShapeFunctions(u, v, w, sf);
+ e->getGradShapeFunctions(u, v, w, grads);
+ for (int j = 0; j < nbNodes; j++){
+ Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ invjac[0][2] * grads[j][2];
+ Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][2] * grads[j][2];
+ Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ invjac[2][2] * grads[j][2];
+ }
+ for (int j = 0; j < nbNodes; j++){
+ for (int k = 0; k < nbNodes; k++){
+ 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;
+ }
+ }
+ }
+
+}
diff --git a/Numeric/gmshHelmholtz.h b/Numeric/gmshHelmholtz.h
new file mode 100644
index 0000000000..f767f24a31
--- /dev/null
+++ b/Numeric/gmshHelmholtz.h
@@ -0,0 +1,37 @@
+// 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_HELMHOLTZ_H_
+#define _GMSH_HELMHOLTZ_H_
+
+#include <complex>
+#include "gmshTermOfFormulation.h"
+#include "gmshFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "GmshMatrix.h"
+
+class gmshHelmholtzTerm : public gmshNodalFemTerm<std::complex<double> > {
+ private:
+ const gmshFunction<std::complex<double> > *_waveNumber;
+ 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:
+ gmshHelmholtzTerm(GModel *gm, gmshFunction<std::complex<double> > *waveNumber,
+ int iField = 0)
+ : gmshNodalFemTerm<std::complex<double> >(gm), _waveNumber(waveNumber), _iField(iField){}
+ virtual void elementMatrix(MElement *e, gmshMatrix<std::complex<double> > &m) const;
+};
+
+#endif
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index 41188e8dc3..3a462294f9 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -15,7 +15,7 @@
class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
private:
- const gmshFunction *_diffusivity;
+ const gmshFunction<double> *_diffusivity;
const int _iField ;
protected:
virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
@@ -27,7 +27,7 @@ class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
*iFieldR = _iField;
}
public:
- gmshLaplaceTerm(GModel *gm, gmshFunction *diffusivity, int iField = 0) :
+ gmshLaplaceTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0) :
gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField){}
virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
};
diff --git a/Plugin/FiniteElement.cpp b/Plugin/FiniteElement.cpp
new file mode 100644
index 0000000000..ce281332e7
--- /dev/null
+++ b/Plugin/FiniteElement.cpp
@@ -0,0 +1,174 @@
+// 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 "GmshConfig.h"
+#include "GModel.h"
+#include "FiniteElement.h"
+#include "gmshAssembler.h"
+#include "gmshLaplace.h"
+#include "gmshHelmholtz.h"
+#include "gmshLinearSystemGmm.h"
+
+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.}
+};
+
+StringXString FiniteElementOptions_String[] = {
+ {GMSH_FULLRC, "Equation", NULL, "Laplace"},
+ {GMSH_FULLRC, "BC1", NULL, "Dirichlet"},
+ {GMSH_FULLRC, "BC2", NULL, "Dirichlet"},
+};
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterFiniteElementPlugin()
+ {
+ return new GMSH_FiniteElementPlugin();
+ }
+}
+
+void GMSH_FiniteElementPlugin::getName(char *name) const
+{
+ strcpy(name, "Finite Element");
+}
+
+void GMSH_FiniteElementPlugin::getInfos(char *author, char *copyright,
+ char *help_text) const
+{
+ strcpy(author, "C. Geuzaine, J.-F. Remacle");
+ strcpy(copyright, "C. Geuzaine, J.-F. Remacle");
+ strcpy(help_text,
+ "Plugin(FiniteElement) solves simple PDEs\n"
+ "using the finite element method. This is only\n"
+ "intended as a demonstration!\n"
+ "Plugin(FiniteElement) is creates a new view.\n");
+}
+
+int GMSH_FiniteElementPlugin::getNbOptions() const
+{
+ return sizeof(FiniteElementOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_FiniteElementPlugin::getOption(int iopt)
+{
+ return &FiniteElementOptions_Number[iopt];
+}
+
+int GMSH_FiniteElementPlugin::getNbOptionsStr() const
+{
+ return sizeof(FiniteElementOptions_String) / sizeof(StringXString);
+}
+
+StringXString *GMSH_FiniteElementPlugin::getOptionStr(int iopt)
+{
+ return &FiniteElementOptions_String[iopt];
+}
+
+void GMSH_FiniteElementPlugin::catchErrorMessage(char *errorMessage) const
+{
+ strcpy(errorMessage, "FiniteElement failed...");
+}
+
+#if defined(HAVE_GMM)
+template<class scalar>
+class solver{
+ public:
+ gmshLinearSystemGmm<scalar> *lsys;
+ 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;
+
+ lsys = new gmshLinearSystemGmm<scalar>;
+ lsys->setPrec(1.e-10);
+ lsys->setNoisy(2);
+ myAssembler = new gmshAssembler<scalar>(lsys);
+
+ GModel *m = GModel::current();
+ std::vector<MVertex*> vertices;
+
+ if(bc1 == "Dirichlet"){
+ m->getMeshVertices(surface1, 2, vertices);
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ myAssembler->fixVertex(vertices[i], 0, 1, value1);
+ }
+ if(bc2 == "Dirichlet"){
+ m->getMeshVertices(surface2, 2, vertices);
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ myAssembler->fixVertex(vertices[i], 0, 1, value2);
+ }
+ m->getMeshVertices(volume, 3, vertices);
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ myAssembler->numberVertex(vertices[i], 0, 1);
+
+ Msg::Info("Assembler: %d dofs, %d fixed",
+ myAssembler->sizeOfR(), myAssembler->sizeOfF());
+ }
+};
+#endif
+
+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;
+
+#if defined(HAVE_GMM)
+ GModel *m = GModel::current();
+ std::map<int, std::vector<GEntity*> > groups[4];
+ m->getPhysicalGroups(groups);
+
+ std::vector<MVertex*> vertices;
+ m->getMeshVertices(volume, 3, vertices);
+ std::map<int, std::vector<double> > data;
+
+ if(equation == "Laplace"){
+ solver<double> s;
+ gmshFunction<double> diffusivity(parameter);
+ gmshLaplaceTerm laplace(m, &diffusivity, 1);
+ for(unsigned int i = 0; i < groups[3][volume].size(); i++)
+ laplace.addToMatrix(*s.myAssembler, groups[3][volume][i]);
+ 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("laplace", "NodeData", m, data, 0.);
+ }
+ else if(equation == "Helmholtz"){
+ solver<std::complex<double> > s;
+ std::complex<double> k(parameter, 0.1);
+ gmshFunction<std::complex<double> > waveNumber(k);
+ gmshHelmholtzTerm helmholtz(m, &waveNumber, 1);
+ for(unsigned int i = 0; i < groups[3][volume].size(); i++)
+ helmholtz.addToMatrix(*s.myAssembler, groups[3][volume][i]);
+ s.lsys->setGmres(1);
+ 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).real());
+ PView *pv = new PView("helmholtz", "NodeData", m, data, 0.);
+ data.clear();
+ for (unsigned int i = 0; i < vertices.size(); i++)
+ data[vertices[i]->getNum()].push_back
+ (s.myAssembler->getDofValue(vertices[i], 0, 1).imag());
+ pv->addStep(m, data, 1.);
+ }
+ else{
+ Msg::Error("Unknown equation to solve");
+ }
+#endif
+ return 0;
+}
diff --git a/Plugin/FiniteElement.h b/Plugin/FiniteElement.h
new file mode 100644
index 0000000000..7eaa1b7948
--- /dev/null
+++ b/Plugin/FiniteElement.h
@@ -0,0 +1,30 @@
+// 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 _FINITE_ELEMENT_H_
+#define _FINITE_ELEMENT_H_
+
+#include "Plugin.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterFiniteElementPlugin();
+}
+
+class GMSH_FiniteElementPlugin : public GMSH_PostPlugin
+{
+public:
+ GMSH_FiniteElementPlugin(){}
+ void getName(char *name) const;
+ void getInfos(char *author, char *copyright, char *helpText) const;
+ void catchErrorMessage(char *errorMessage) const;
+ int getNbOptions() const;
+ StringXNumber* getOption(int iopt);
+ int getNbOptionsStr() const;
+ StringXString* getOptionStr(int iopt);
+ PView *execute(PView *);
+};
+
+#endif
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 0cd3036a2d..cc9f78fe04 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -9,7 +9,8 @@ LIB = ../lib/libGmshPlugin${LIBEXT}
INC = ${DASH}I../Common ${DASH}I../Graphics ${DASH}I../Geo\
${DASH}I../Mesh ${DASH}I../Post ${DASH}I../Fltk ${DASH}I../Numeric\
- ${DASH}I../contrib/ANN/include ${DASH}I../contrib/MathEval
+ ${DASH}I../contrib/ANN/include ${DASH}I../contrib/MathEval\
+ ${DASH}I../contrib/gmm
CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
@@ -31,7 +32,8 @@ SRC = Plugin.cpp PluginManager.cpp\
Integrate.cpp Gradient.cpp Curl.cpp Divergence.cpp\
Annotate.cpp Remove.cpp\
Probe.cpp\
- HarmonicToTime.cpp ModulusPhase.cpp
+ HarmonicToTime.cpp ModulusPhase.cpp\
+ FiniteElement.cpp
OBJ = ${SRC:.cpp=${OBJEXT}}
@@ -77,7 +79,7 @@ PluginManager${OBJEXT}: PluginManager.cpp ../Common/GmshConfig.h Plugin.h \
LongitudeLatitude.h Triangulate.h Warp.h SphericalRaise.h \
Eigenvectors.h Eigenvalues.h Lambda2.h Evaluate.h ../Post/OctreePost.h \
../Common/Octree.h ../Common/OctreeInternals.h Probe.h FieldView.h \
- GSHHS.h ../Common/Context.h ../Geo/CGNSOptions.h \
+ GSHHS.h FiniteElement.h ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/meshPartitionOptions.h
Levelset${OBJEXT}: Levelset.cpp Levelset.h Plugin.h ../Common/Options.h \
../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
@@ -306,3 +308,26 @@ ModulusPhase${OBJEXT}: ModulusPhase.cpp ModulusPhase.h Plugin.h \
../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
../Post/PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h
+FiniteElement${OBJEXT}: FiniteElement.cpp ../Common/GmshConfig.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 FiniteElement.h Plugin.h ../Common/Options.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/gmshTermOfFormulation.h ../Numeric/GmshMatrix.h \
+ ../Numeric/gmshAssembler.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 ../Numeric/gmshFunction.h ../Common/Gmsh.h \
+ ../Common/GmshMessage.h ../Numeric/GmshMatrix.h \
+ ../Numeric/gmshHelmholtz.h ../Numeric/gmshTermOfFormulation.h \
+ ../Numeric/gmshFunction.h ../Numeric/GmshMatrix.h \
+ ../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index ccdecc79e9..b599f8e69e 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -39,6 +39,7 @@
#include "Probe.h"
#include "FieldView.h"
#include "GSHHS.h"
+#include "FiniteElement.h"
#include "Context.h"
#if !defined(HAVE_NO_DLL)
@@ -223,6 +224,8 @@ void PluginManager::registerDefaultPlugins()
allPlugins.insert(std::pair<const char*, GMSH_Plugin*>
("CutParametric", GMSH_RegisterCutParametricPlugin()));
#endif
+ allPlugins.insert(std::pair<const char*, GMSH_Plugin*>
+ ("FiniteElement", GMSH_RegisterFiniteElementPlugin()));
}
#if defined(HAVE_FLTK)
--
GitLab