From 2ec0553a729292aa0aafc2245cf1c63ea26464a1 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Thu, 17 Sep 2009 12:05:21 +0000
Subject: [PATCH] more solver stuff
---
Solver/elasticitySolver.h | 54 ++++++++++++++++++++++
Solver/helmholtzTerm.h | 94 +++++++++++++++++++++++++++++++++++++++
Solver/laplaceTerm.h | 35 +++++++++++++++
3 files changed, 183 insertions(+)
create mode 100644 Solver/elasticitySolver.h
create mode 100644 Solver/helmholtzTerm.h
create mode 100644 Solver/laplaceTerm.h
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
new file mode 100644
index 0000000000..276d3db572
--- /dev/null
+++ b/Solver/elasticitySolver.h
@@ -0,0 +1,54 @@
+#ifndef _ELASTICITY_SOLVER_H_
+#define _ELASTICITY_SOLVER_H_
+
+#include <map>
+#include <string>
+#include "SVector3.h"
+#include "dofManager.h"
+
+class GModel;
+
+// an elastic solver ...
+namespace gsolver {
+ class elasticitySolver
+ {
+ protected:
+ GModel *pModel;
+ int _dim,_tag;
+ dofManager<double,double> *pAssembler;
+ // young modulus and poisson coefficient per physical
+ std::map<int, std::pair<double,double> > elasticConstants;
+ // imposed nodal forces
+ std::map<int, SVector3 > nodalForces;
+ // imposed line forces
+ std::map<int, SVector3 > lineForces;
+ // imposed face forces
+ std::map<int, SVector3 > faceForces;
+ // imposed volume forces
+ std::map<int, SVector3 > volumeForces;
+ // imposed nodal displacements
+ std::map<std::pair<int,int>, double> nodalDisplacements;
+ // imposed edge displacements
+ std::map<std::pair<int,int>, double> edgeDisplacements;
+ // imposed face displacements
+ std::map<std::pair<int,int>, double> faceDisplacements;
+ public:
+ elasticitySolver (int tag) : _tag(tag){}
+ void addNodalForces (int iNode, const SVector3 & f )
+ {nodalForces[iNode] = f;}
+ void addNodalDisplacement (int iNode, int dir, double val )
+ {nodalDisplacements[std::make_pair(iNode,dir)] = val;}
+ void addElasticConstants (double e, double nu, int physical)
+ {elasticConstants[physical] = std::make_pair(e,nu);}
+ void setMesh (const std::string &meshFileName);
+ virtual void solve ();
+ void readInputFile (const std::string &meshFileName);
+ // PView * buildDisplacementView (const std::string &postFileName);
+ // PView * buildVonMisesView(const std::string &postFileName);
+ // std::pair<PView *, PView*> buildErrorEstimateView (const std::string &errorFileName, double, int);
+ // std::pair<PView *, PView*> buildErrorEstimateView (const std::string &errorFileName, const gmshElasticityData &ref, double, int);
+ }; // end of class definition
+}// end of namespace
+#endif
+
+
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
new file mode 100644
index 0000000000..c2a6c3179e
--- /dev/null
+++ b/Solver/helmholtzTerm.h
@@ -0,0 +1,94 @@
+// 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 _HELMHOLTZ_H_
+#define _HELMHOLTZ_H_
+
+#include <assert.h>
+#include "termOfFormulation.h"
+#include "gmshFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "GmshMatrix.h"
+#include "Numeric.h"
+
+namespace gsolver {
+ // \nabla \cdot k \nabla U + a U
+ template<class scalar>
+ class helmoltzTerm : public femTerm<scalar,scalar> {
+ protected:
+ const gmshFunction<scalar> *_k, *_a;
+ const int _iFieldR;
+ int _iFieldC ;
+ public:
+ // one dof per vertex (nodal fem)
+ virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
+ virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
+ Dof getLocalDofR(MElement *e, int iRow) const
+ {return Dof (e->getVertex(iRow)->getNum(),_iFieldR);}
+ Dof getLocalDofC(MElement *e, int iRow) const
+ {return Dof (e->getVertex(iRow)->getNum(),_iFieldC);}
+ public:
+ helmoltzTerm(GModel *gm, int iFieldR, int iFieldC,
+ gmshFunction<scalar> *k,
+ gmshFunction<scalar> *a) :
+ femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR),_iFieldC(iFieldC)
+ {}
+ virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const{
+ // compute integration rule
+ const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() : 2 * (e->getPolynomialOrder() - 1);
+ int npts;IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ // get the number of nodes
+ const int nbNodes = e->getNumVertices();
+ // assume a maximum of 100 nodes
+ assert (nbNodes < 100);
+ double jac[3][3];
+ double invjac[3][3];
+ double Grads[100][3], grads[100][3];
+ double sf[100];
+ // set the local matrix to 0
+ m.set_all(0.);
+ // loop over integration points
+ for (int i = 0; i < npts; i++){
+ // compute stuff at this point
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
+ const double weightDetJ = GP[i].weight * e->getJacobian(u, v, w, jac);
+ SPoint3 p; e->pnt(u, v, w, p);
+ const scalar K = (*_k)(p.x(), p.y(), p.z());
+ const scalar A = _a ? (*_a)(p.x(), p.y(), p.z()) : 0.0;
+ inv3x3(jac, invjac) ;
+ e->getGradShapeFunctions(u, v, w, grads);
+ if (_a) e->getShapeFunctions(u, v, w, sf);
+ 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];
+ if (!_a) sf[j] = 0;
+ }
+ for (int j = 0; j < nbNodes; j++){
+ for (int k = 0; k <= j; k++){
+ m(j, k) += (K* (Grads[j][0] * Grads[k][0] +
+ Grads[j][1] * Grads[k][1] +
+ Grads[j][2] * Grads[k][2]) + A * sf[j] * sf[k]) *
+ weightDetJ;
+ }
+ }
+ }
+ for (int j = 0; j < nbNodes; j++)
+ for (int k = 0; k < j; k++)
+ m(k, j) = m(j, k);
+
+ }// end of element matrix
+ };// end of the class definition
+}// end of namespace
+
+#endif
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
new file mode 100644
index 0000000000..a1c0394133
--- /dev/null
+++ b/Solver/laplaceTerm.h
@@ -0,0 +1,35 @@
+// 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 _LAPLACE_H_
+#define _LAPLACE_H_
+
+#include "helmholtzTerm.h"
+
+namespace gsolver {
+ // \nabla \cdot k \nabla U
+ template<class scalar>
+ class laplaceTerm : public helmholtzTerm<scalar> {
+ gmshFunction<scalar> *ONE;
+ public:
+ laplaceTerm(GModel *gm, int iFieldR) :
+ helmholtzTerm<scalar>(gm,iFieldR,iFieldR,0,0)
+ {
+ ONE = new gmshFunction<scalar>(1.0);
+ _k = ONE;
+ }
+ laplaceTerm(GModel *gm, int iFieldR, gmshFunction<scalar> *k) :
+ helmholtzTerm<scalar>(gm,iFieldR,iFieldR,k,0),ONE(0)
+ {}
+ laplaceTerm(GModel *gm, int iFieldR, int iFieldC, gmshFunction<scalar> *k) :
+ helmholtzTerm<scalar>(gm,iFieldR,iFieldC,k,0),ONE(0)
+ {}
+ virtual ~laplaceTerm(){
+ if(ONE) delete ONE;
+ }
+ };// end of the class definition
+}// end of namespace
+
+#endif
--
GitLab