From 587f74b57d227690d964ae444aed874d0aac21b6 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 8 Feb 2009 13:53:37 +0000
Subject: [PATCH] templatize gmshAssembler
---
Geo/GFaceCompound.cpp | 4 +-
Mesh/gmshSmoothHighOrder.cpp | 6 +-
Numeric/Makefile | 2 -
Numeric/gmshAssembler.cpp | 71 -------------------
Numeric/gmshAssembler.h | 90 +++++++++++++++++++-----
Numeric/gmshTermOfFormulation.cpp | 110 ------------------------------
Numeric/gmshTermOfFormulation.h | 80 +++++++++++++++-------
7 files changed, 132 insertions(+), 231 deletions(-)
delete mode 100644 Numeric/gmshAssembler.cpp
delete mode 100644 Numeric/gmshTermOfFormulation.cpp
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 08e1c2df2b..4269251d7e 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -43,7 +43,7 @@ class gmshGradientBasedDiffusivity : public gmshFunction
}
};
-static void fixEdgeToValue (GEdge *ed, double value, gmshAssembler &myAssembler)
+static void fixEdgeToValue(GEdge *ed, double value, gmshAssembler<double> &myAssembler)
{
myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value);
myAssembler.fixVertex(ed->getEndVertex()->mesh_vertices[0], 0, 1, value);
@@ -205,7 +205,7 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
#else
gmshLinearSystemFull<double> lsys;
#endif
- gmshAssembler myAssembler(&lsys);
+ gmshAssembler<double> myAssembler(&lsys);
gmshGradientBasedDiffusivity diffusivity(coordinates);
if (ITER > 0) diffusivity.setComponent(_isU);
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 09a4e5557e..59f611c8a5 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -95,7 +95,7 @@ void gmshHighOrderSmoother::smooth(GRegion *gr)
void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
{
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
- gmshAssembler myAssembler(lsys);
+ gmshAssembler<double> myAssembler(lsys);
gmshElasticityTerm El(0,1.0,.333,getTag());
std::vector<MElement*> v, layer;
@@ -434,7 +434,7 @@ void localHarmonicMapping(GModel *gm,
std::vector<MVertex*> &e) {
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
- gmshAssembler myAssembler(lsys);
+ gmshAssembler<double> myAssembler(lsys);
gmshFunction f(1.0);
gmshLaplaceTerm Laplace (gm,&f,0);
@@ -449,7 +449,7 @@ void localHarmonicMapping(GModel *gm,
lsys->systemSolve();
gmshLinearSystemGmm<double> *lsys1 = new gmshLinearSystemGmm<double>;
- gmshAssembler myAssembler1(lsys1);
+ gmshAssembler<double> myAssembler1(lsys1);
gmshLaplaceTerm Laplace1 (gm,&f,1);
myAssembler1.fixVertex ( n2 , 0 , 1 , -1.0);
diff --git a/Numeric/Makefile b/Numeric/Makefile
index caef387d32..0ead3c621f 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -13,8 +13,6 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
SRC = Numeric.cpp\
GmshMatrix.cpp\
- gmshAssembler.cpp\
- gmshTermOfFormulation.cpp\
gmshLaplace.cpp\
gmshElasticity.cpp\
EigSolve.cpp\
diff --git a/Numeric/gmshAssembler.cpp b/Numeric/gmshAssembler.cpp
deleted file mode 100644
index 19160913a5..0000000000
--- a/Numeric/gmshAssembler.cpp
+++ /dev/null
@@ -1,71 +0,0 @@
-// 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 "MVertex.h"
-#include "gmshAssembler.h"
-
-void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
- MVertex *vC, int iCompC, int iFieldC,
- double val)
-{
- if (!lsys->isAllocated()){
- lsys->allocate(numbering.size());
- }
-
- std::map<gmshDofKey, int>::iterator
- itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
- if (itR != numbering.end()){
- std::map<gmshDofKey, int>::iterator
- itC = numbering.find(gmshDofKey(vC, iCompC, iFieldC));
- if (itC != numbering.end()){
- lsys->addToMatrix(itR->second, itC->second, val);
- }
- else {
- std::map<gmshDofKey, double>::iterator
- itF = fixed.find(gmshDofKey(vC, iCompC, iFieldC));
- if (itF != fixed.end()){
- lsys->addToRightHandSide(itR->second, -val*itF->second);
- }
- else{
- std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator
- itConstrC = constraints.find(gmshDofKey(vC, iCompC, iFieldC));
- if (itConstrC != constraints.end()){
- for (unsigned int i = 0; i < itConstrC->second.size(); i++){
- gmshDofKey &dofKeyConstrC = itConstrC->second[i].first;
- double valConstrC = itConstrC->second[i].second;
- assemble(vR, iCompR, iFieldR,
- dofKeyConstrC.v, dofKeyConstrC.comp, dofKeyConstrC.field,
- val * valConstrC);
- }
- }
- }
- }
- }
- else{
- std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator
- itConstrR = constraints.find(gmshDofKey(vR, iCompR, iFieldR));
- if (itConstrR != constraints.end()){
- for (unsigned int i = 0; i < itConstrR->second.size(); i++){
- gmshDofKey &dofKeyConstrR = itConstrR->second[i].first;
- double valConstrR = itConstrR->second[i].second;
- assemble(dofKeyConstrR.v,dofKeyConstrR.comp, dofKeyConstrR.field,
- vC, iCompC, iFieldC,
- val * valConstrR);
- }
- }
- }
-}
-
-void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
- double val)
-{
- if (!lsys->isAllocated())lsys->allocate(numbering.size());
- std::map<gmshDofKey, int>::iterator
- itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
- if (itR != numbering.end()){
- lsys->addToRightHandSide(itR->second, val);
- }
-}
-
diff --git a/Numeric/gmshAssembler.h b/Numeric/gmshAssembler.h
index 158e68a7e3..195edbf328 100644
--- a/Numeric/gmshAssembler.h
+++ b/Numeric/gmshAssembler.h
@@ -30,18 +30,20 @@ struct gmshDofKey{
}
};
+template<class scalar>
class gmshAssembler {
- std::map<gmshDofKey, int> numbering;
- std::map<gmshDofKey, double> fixed;
- std::map<gmshDofKey, std::vector<std::pair<gmshDofKey,double> > > constraints;
- gmshLinearSystem<double> *lsys;
-public:
- gmshAssembler(gmshLinearSystem<double> *l) : lsys(l) {}
+ private:
+ std::map<gmshDofKey, int> numbering;
+ std::map<gmshDofKey, scalar> fixed;
+ std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > > constraints;
+ gmshLinearSystem<scalar> *lsys;
+ public:
+ gmshAssembler(gmshLinearSystem<scalar> *l) : lsys(l) {}
inline void constraintVertex(MVertex*v, int iComp, int iField,
std::vector<MVertex*> &verts,
- std::vector<double> &coeffs)
+ std::vector<scalar> &coeffs)
{
- std::vector<std::pair<gmshDofKey, double> > constraint;
+ std::vector<std::pair<gmshDofKey, scalar> > constraint;
gmshDofKey key(v, iComp, iField);
for (unsigned int i = 0; i < verts.size(); i++){
gmshDofKey key2(verts[i], iComp, iField);
@@ -73,15 +75,15 @@ public:
numbering[key] = size;
}
}
- inline void fixVertex(MVertex*v, int iComp, int iField, double val)
+ inline void fixVertex(MVertex*v, int iComp, int iField, scalar val)
{
fixed[gmshDofKey(v, iComp, iField)] = val;
}
- inline double getDofValue(MVertex*v, int iComp, int iField) const
+ inline scalar getDofValue(MVertex *v, int iComp, int iField) const
{
gmshDofKey key(v, iComp, iField);
{
- std::map<gmshDofKey, double>::const_iterator it = fixed.find(key);
+ typename std::map<gmshDofKey, scalar>::const_iterator it = fixed.find(key);
if (it != fixed.end()) return it->second;
}
{
@@ -90,13 +92,13 @@ public:
return lsys->getFromSolution(it->second);
}
{
- std::map<gmshDofKey, std::vector<std::pair<gmshDofKey,double> > >::
+ typename std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > >::
const_iterator itConstr = constraints.find(key);
if (itConstr != constraints.end()){
- double val = 0;
+ scalar val = 0.;
for (unsigned int i = 0; i < itConstr->second.size(); i++){
const gmshDofKey &dofKeyConstr = itConstr->second[i].first;
- double valConstr = itConstr->second[i].second;
+ scalar valConstr = itConstr->second[i].second;
val += getDofValue(dofKeyConstr.v, dofKeyConstr.comp, dofKeyConstr.field)
* valConstr;
}
@@ -106,10 +108,62 @@ public:
return 0.0;
}
void assemble(MVertex *vR , int iCompR, int iFieldR,
- MVertex *vC , int iCompC, int iFieldC,
- double val);
- void assemble(MVertex *vR , int iCompR, int iFieldR,
- double val);
+ MVertex *vC , int iCompC, int iFieldC, scalar val)
+ {
+ if (!lsys->isAllocated()) lsys->allocate(numbering.size());
+
+ std::map<gmshDofKey, int>::iterator
+ itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
+ if (itR != numbering.end()){
+ std::map<gmshDofKey, int>::iterator
+ itC = numbering.find(gmshDofKey(vC, iCompC, iFieldC));
+ if (itC != numbering.end()){
+ lsys->addToMatrix(itR->second, itC->second, val);
+ }
+ else {
+ typename std::map<gmshDofKey, scalar>::iterator
+ itF = fixed.find(gmshDofKey(vC, iCompC, iFieldC));
+ if (itF != fixed.end()){
+ lsys->addToRightHandSide(itR->second, -val*itF->second);
+ }
+ else{
+ typename std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > >::
+ iterator itConstrC = constraints.find(gmshDofKey(vC, iCompC, iFieldC));
+ if (itConstrC != constraints.end()){
+ for (unsigned int i = 0; i < itConstrC->second.size(); i++){
+ gmshDofKey &dofKeyConstrC = itConstrC->second[i].first;
+ scalar valConstrC = itConstrC->second[i].second;
+ assemble(vR, iCompR, iFieldR,
+ dofKeyConstrC.v, dofKeyConstrC.comp, dofKeyConstrC.field,
+ val * valConstrC);
+ }
+ }
+ }
+ }
+ }
+ else{
+ typename std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > >::
+ iterator itConstrR = constraints.find(gmshDofKey(vR, iCompR, iFieldR));
+ if (itConstrR != constraints.end()){
+ for (unsigned int i = 0; i < itConstrR->second.size(); i++){
+ gmshDofKey &dofKeyConstrR = itConstrR->second[i].first;
+ scalar valConstrR = itConstrR->second[i].second;
+ assemble(dofKeyConstrR.v,dofKeyConstrR.comp, dofKeyConstrR.field,
+ vC, iCompC, iFieldC,
+ val * valConstrR);
+ }
+ }
+ }
+ }
+ void assemble(MVertex *vR , int iCompR, int iFieldR, scalar val)
+ {
+ if (!lsys->isAllocated())lsys->allocate(numbering.size());
+ std::map<gmshDofKey, int>::iterator
+ itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
+ if (itR != numbering.end()){
+ lsys->addToRightHandSide(itR->second, val);
+ }
+ }
int sizeOfR() const { return numbering.size(); }
int sizeOfF() const { return fixed.size(); }
};
diff --git a/Numeric/gmshTermOfFormulation.cpp b/Numeric/gmshTermOfFormulation.cpp
deleted file mode 100644
index 4afb8d88ae..0000000000
--- a/Numeric/gmshTermOfFormulation.cpp
+++ /dev/null
@@ -1,110 +0,0 @@
-// 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 "Gmsh.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "GmshMatrix.h"
-#include "gmshTermOfFormulation.h"
-#include "gmshFunction.h"
-#include "gmshAssembler.h"
-
-gmshNodalFemTerm::~gmshNodalFemTerm ()
-{
-}
-
-void gmshNodalFemTerm::addDirichlet(int physical,
- int dim,
- int comp,
- int field,
- const gmshFunction & e,
- gmshAssembler &lsys)
-{
-}
-
-void gmshNodalFemTerm::addNeumann(int physical,
- int dim,
- int comp,
- int field,
- const gmshFunction & fct,
- gmshAssembler &lsys)
-{
-}
-
-void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys) const
-{
- if (_gm->getNumRegions()){
- for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
- addToMatrix(lsys, *it);
- }
- }
- else if(_gm->getNumFaces()){
- for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
- addToMatrix(lsys, *it);
- }
- }
-}
-
-void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys,const std::vector<MElement*> &v) const
-{
- for (unsigned int i = 0; i < v.size(); i++)
- addToMatrix(lsys, v[i]);
-}
-
-void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys,
- gmshMatrix<double> &localMatrix,
- MElement *e) const
-{
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
- for (int j = 0; j < nbR; j++){
- MVertex *vR;int iCompR,iFieldR;
- getLocalDofR (e, j, &vR, &iCompR, &iFieldR);
- for (int k = 0; k < nbC; k++){
- MVertex *vC;
- int iCompC, iFieldC;
- getLocalDofC(e, k, &vC, &iCompC, &iFieldC);
- lsys.assemble(vR, iCompR, iFieldR,
- vC, iCompC, iFieldC,
- localMatrix(j, k));
- }
- }
-}
-
-void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys, MElement *e) const
-{
- const int nbR = sizeOfR(e);
- const int nbC = sizeOfC(e);
- gmshMatrix<double> localMatrix (nbR, nbC);
- elementMatrix(e, localMatrix);
- addToMatrix(lsys, localMatrix, e);
-}
-
-void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys, GEntity *ge) const
-{
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement(i);
- addToMatrix(lsys, e);
- }
-}
-
-void gmshNodalFemTerm::addToRightHandSide(gmshAssembler &lsys) const
-{
- if (_gm->getNumRegions()){
- for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
- addToRightHandSide(lsys, *it);
- }
- }
- else if(_gm->getNumFaces()){
- for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
- addToRightHandSide(lsys, *it);
- }
- }
-}
-
-void gmshNodalFemTerm::addToRightHandSide(gmshAssembler &lsys, GEntity *ge) const
-{
- throw;
-}
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index 966a038965..d7b9fc8dbb 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -10,13 +10,9 @@
#include <map>
#include <vector>
#include "GmshMatrix.h"
-
-class GModel;
-class GEntity;
-class MElement;
-class MVertex;
-class gmshFunction;
-class gmshAssembler;
+#include "gmshAssembler.h"
+#include "GModel.h"
+#include "MElement.h"
class gmshTermOfFormulation {
protected:
@@ -24,8 +20,7 @@ class gmshTermOfFormulation {
public:
gmshTermOfFormulation(GModel *gm) : _gm(gm) {}
virtual ~gmshTermOfFormulation(){}
- virtual void addToMatrix(gmshAssembler&) const = 0;
- virtual void addToRightHandSide(gmshAssembler&) const = 0;
+ virtual void addToMatrix(gmshAssembler<double> &) const = 0;
};
// a nodal finite element term : variables are always defined at nodes
@@ -47,23 +42,58 @@ class gmshNodalFemTerm : public gmshTermOfFormulation {
}
public:
gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation(gm) {}
- virtual ~gmshNodalFemTerm ();
- // compute the element matrix
+ virtual ~gmshNodalFemTerm (){}
virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const = 0;
-
- void addToMatrix(gmshAssembler &J, MElement *e) const;
- void addToMatrix(gmshAssembler &J, GEntity *ge) const;
- void addToMatrix(gmshAssembler &J) const;
- void addToMatrix(gmshAssembler &J,const std::vector<MElement*> &) const;
- void addToMatrix(gmshAssembler &Jac, gmshMatrix<double> &localMatrix, MElement *e) const;
-
- void addDirichlet(int physical, int dim, int comp, int field, const gmshFunction &e,
- gmshAssembler &);
- void addNeumann(int physical, int dim, int icomp, int field, const gmshFunction &e,
- gmshAssembler &);
- void addToRightHandSide(gmshAssembler &J, GEntity *ge) const;
- void addToRightHandSide(gmshAssembler &r) const;
-
+ void addToMatrix(gmshAssembler<double> &lsys) const
+ {
+ if (_gm->getNumRegions()){
+ for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
+ addToMatrix(lsys, *it);
+ }
+ }
+ else if(_gm->getNumFaces()){
+ for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
+ addToMatrix(lsys, *it);
+ }
+ }
+ }
+ void addToMatrix(gmshAssembler<double> &lsys, GEntity *ge) const
+ {
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement(i);
+ addToMatrix(lsys, e);
+ }
+ }
+ void addToMatrix(gmshAssembler<double> &lsys, MElement *e) const
+ {
+ const int nbR = sizeOfR(e);
+ const int nbC = sizeOfC(e);
+ gmshMatrix<double> localMatrix (nbR, nbC);
+ elementMatrix(e, localMatrix);
+ addToMatrix(lsys, localMatrix, e);
+ }
+ void addToMatrix(gmshAssembler<double> &lsys, const std::vector<MElement*> &v) const
+ {
+ for (unsigned int i = 0; i < v.size(); i++)
+ addToMatrix(lsys, v[i]);
+ }
+ void addToMatrix(gmshAssembler<double> &lsys, gmshMatrix<double> &localMatrix,
+ MElement *e) const
+ {
+ const int nbR = sizeOfR(e);
+ const int nbC = sizeOfC(e);
+ for (int j = 0; j < nbR; j++){
+ MVertex *vR;
+ int iCompR, iFieldR;
+ getLocalDofR(e, j, &vR, &iCompR, &iFieldR);
+ for (int k = 0; k < nbC; k++){
+ MVertex *vC;
+ int iCompC, iFieldC;
+ getLocalDofC(e, k, &vC, &iCompC, &iFieldC);
+ lsys.assemble(vR, iCompR, iFieldR, vC, iCompC, iFieldC, localMatrix(j, k));
+ }
+ }
+ }
};
#endif
--
GitLab