Skip to content
Snippets Groups Projects
Commit 587f74b5 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

templatize gmshAssembler

parent 6b28c99d
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GFaceCompound.cpp
+ 2
− 2
View file @ 587f74b5
...@@ -43,7 +43,7 @@ class gmshGradientBasedDiffusivity : public gmshFunction ...@@ -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->getBeginVertex()->mesh_vertices[0], 0, 1, value);
myAssembler.fixVertex(ed->getEndVertex()->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 ...@@ -205,7 +205,7 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
#else #else
gmshLinearSystemFull<double> lsys; gmshLinearSystemFull<double> lsys;
#endif #endif
gmshAssembler myAssembler(&lsys); gmshAssembler<double> myAssembler(&lsys);
gmshGradientBasedDiffusivity diffusivity(coordinates); gmshGradientBasedDiffusivity diffusivity(coordinates);
if (ITER > 0) diffusivity.setComponent(_isU); if (ITER > 0) diffusivity.setComponent(_isU);
... ...
......
Mesh/gmshSmoothHighOrder.cpp
+ 3
− 3
View file @ 587f74b5
...@@ -95,7 +95,7 @@ void gmshHighOrderSmoother::smooth(GRegion *gr) ...@@ -95,7 +95,7 @@ void gmshHighOrderSmoother::smooth(GRegion *gr)
void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all) void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
{ {
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>; gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
gmshAssembler myAssembler(lsys); gmshAssembler<double> myAssembler(lsys);
gmshElasticityTerm El(0,1.0,.333,getTag()); gmshElasticityTerm El(0,1.0,.333,getTag());
std::vector<MElement*> v, layer; std::vector<MElement*> v, layer;
...@@ -434,7 +434,7 @@ void localHarmonicMapping(GModel *gm, ...@@ -434,7 +434,7 @@ void localHarmonicMapping(GModel *gm,
std::vector<MVertex*> &e) { std::vector<MVertex*> &e) {
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>; gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
gmshAssembler myAssembler(lsys); gmshAssembler<double> myAssembler(lsys);
gmshFunction f(1.0); gmshFunction f(1.0);
gmshLaplaceTerm Laplace (gm,&f,0); gmshLaplaceTerm Laplace (gm,&f,0);
...@@ -449,7 +449,7 @@ void localHarmonicMapping(GModel *gm, ...@@ -449,7 +449,7 @@ void localHarmonicMapping(GModel *gm,
lsys->systemSolve(); lsys->systemSolve();
gmshLinearSystemGmm<double> *lsys1 = new gmshLinearSystemGmm<double>; gmshLinearSystemGmm<double> *lsys1 = new gmshLinearSystemGmm<double>;
gmshAssembler myAssembler1(lsys1); gmshAssembler<double> myAssembler1(lsys1);
gmshLaplaceTerm Laplace1 (gm,&f,1); gmshLaplaceTerm Laplace1 (gm,&f,1);
myAssembler1.fixVertex ( n2 , 0 , 1 , -1.0); myAssembler1.fixVertex ( n2 , 0 , 1 , -1.0);
... ...
......
Numeric/Makefile
+ 0
− 2
View file @ 587f74b5
...@@ -13,8 +13,6 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE} ...@@ -13,8 +13,6 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
SRC = Numeric.cpp\ SRC = Numeric.cpp\
GmshMatrix.cpp\ GmshMatrix.cpp\
gmshAssembler.cpp\
gmshTermOfFormulation.cpp\
gmshLaplace.cpp\ gmshLaplace.cpp\
gmshElasticity.cpp\ gmshElasticity.cpp\
EigSolve.cpp\ EigSolve.cpp\
... ...
......
Numeric/gmshAssembler.cpp deleted 100644 → 0
+ 0
− 71
View file @ 6b28c99d
// 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);
}
}
Numeric/gmshAssembler.h
+ 72
− 18
View file @ 587f74b5
...@@ -30,18 +30,20 @@ struct gmshDofKey{ ...@@ -30,18 +30,20 @@ struct gmshDofKey{
} }
}; };
template<class scalar>
class gmshAssembler { class gmshAssembler {
private:
std::map<gmshDofKey, int> numbering; std::map<gmshDofKey, int> numbering;
std::map<gmshDofKey, double> fixed; std::map<gmshDofKey, scalar> fixed;
std::map<gmshDofKey, std::vector<std::pair<gmshDofKey,double> > > constraints; std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > > constraints;
gmshLinearSystem<double> *lsys; gmshLinearSystem<scalar> *lsys;
public: public:
gmshAssembler(gmshLinearSystem<double> *l) : lsys(l) {} gmshAssembler(gmshLinearSystem<scalar> *l) : lsys(l) {}
inline void constraintVertex(MVertex*v, int iComp, int iField, inline void constraintVertex(MVertex*v, int iComp, int iField,
std::vector<MVertex*> &verts, 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); gmshDofKey key(v, iComp, iField);
for (unsigned int i = 0; i < verts.size(); i++){ for (unsigned int i = 0; i < verts.size(); i++){
gmshDofKey key2(verts[i], iComp, iField); gmshDofKey key2(verts[i], iComp, iField);
...@@ -73,15 +75,15 @@ public: ...@@ -73,15 +75,15 @@ public:
numbering[key] = size; 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; 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); 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; if (it != fixed.end()) return it->second;
} }
{ {
...@@ -90,13 +92,13 @@ public: ...@@ -90,13 +92,13 @@ public:
return lsys->getFromSolution(it->second); 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); const_iterator itConstr = constraints.find(key);
if (itConstr != constraints.end()){ if (itConstr != constraints.end()){
double val = 0; scalar val = 0.;
for (unsigned int i = 0; i < itConstr->second.size(); i++){ for (unsigned int i = 0; i < itConstr->second.size(); i++){
const gmshDofKey &dofKeyConstr = itConstr->second[i].first; 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) val += getDofValue(dofKeyConstr.v, dofKeyConstr.comp, dofKeyConstr.field)
* valConstr; * valConstr;
} }
...@@ -106,10 +108,62 @@ public: ...@@ -106,10 +108,62 @@ public:
return 0.0; return 0.0;
} }
void assemble(MVertex *vR , int iCompR, int iFieldR, void assemble(MVertex *vR , int iCompR, int iFieldR,
MVertex *vC , int iCompC, int iFieldC, MVertex *vC , int iCompC, int iFieldC, scalar val)
double val); {
void assemble(MVertex *vR , int iCompR, int iFieldR, if (!lsys->isAllocated()) lsys->allocate(numbering.size());
double val);
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 sizeOfR() const { return numbering.size(); }
int sizeOfF() const { return fixed.size(); } int sizeOfF() const { return fixed.size(); }
}; };
... ...
......
Numeric/gmshTermOfFormulation.cpp deleted 100644 → 0
+ 0
− 110
View file @ 6b28c99d
// 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;
}
Numeric/gmshTermOfFormulation.h
+ 55
− 25
View file @ 587f74b5
...@@ -10,13 +10,9 @@ ...@@ -10,13 +10,9 @@
#include <map> #include <map>
#include <vector> #include <vector>
#include "GmshMatrix.h" #include "GmshMatrix.h"
#include "gmshAssembler.h"
class GModel; #include "GModel.h"
class GEntity; #include "MElement.h"
class MElement;
class MVertex;
class gmshFunction;
class gmshAssembler;
class gmshTermOfFormulation { class gmshTermOfFormulation {
protected: protected:
...@@ -24,8 +20,7 @@ class gmshTermOfFormulation { ...@@ -24,8 +20,7 @@ class gmshTermOfFormulation {
public: public:
gmshTermOfFormulation(GModel *gm) : _gm(gm) {} gmshTermOfFormulation(GModel *gm) : _gm(gm) {}
virtual ~gmshTermOfFormulation(){} virtual ~gmshTermOfFormulation(){}
virtual void addToMatrix(gmshAssembler&) const = 0; virtual void addToMatrix(gmshAssembler<double> &) const = 0;
virtual void addToRightHandSide(gmshAssembler&) const = 0;
}; };
// a nodal finite element term : variables are always defined at nodes // a nodal finite element term : variables are always defined at nodes
...@@ -47,23 +42,58 @@ class gmshNodalFemTerm : public gmshTermOfFormulation { ...@@ -47,23 +42,58 @@ class gmshNodalFemTerm : public gmshTermOfFormulation {
} }
public: public:
gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation(gm) {} gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation(gm) {}
virtual ~gmshNodalFemTerm (); virtual ~gmshNodalFemTerm (){}
// compute the element matrix
virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const = 0; virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const = 0;
void addToMatrix(gmshAssembler<double> &lsys) const
void addToMatrix(gmshAssembler &J, MElement *e) const; {
void addToMatrix(gmshAssembler &J, GEntity *ge) const; if (_gm->getNumRegions()){
void addToMatrix(gmshAssembler &J) const; for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
void addToMatrix(gmshAssembler &J,const std::vector<MElement*> &) const; addToMatrix(lsys, *it);
void addToMatrix(gmshAssembler &Jac, gmshMatrix<double> &localMatrix, MElement *e) const; }
}
void addDirichlet(int physical, int dim, int comp, int field, const gmshFunction &e, else if(_gm->getNumFaces()){
gmshAssembler &); for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
void addNeumann(int physical, int dim, int icomp, int field, const gmshFunction &e, addToMatrix(lsys, *it);
gmshAssembler &); }
void addToRightHandSide(gmshAssembler &J, GEntity *ge) const; }
void addToRightHandSide(gmshAssembler &r) const; }
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 #endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment