diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index cd1f1665e2dfcf313f5e42cd38158fc5bd8f2d60..f20b610a095b8b809b82de1ae22aa3f9aaa6ec1f 100644 --- a/Geo/GFaceCompound.cpp +++ b/Geo/GFaceCompound.cpp @@ -12,7 +12,7 @@ static void fixEdgeToValue (GEdge *ed, 0, 1, value); myAssembler.fixVertex(ed->getEndVertex()->mesh_vertices[0], 0, 1, value); - for (int i = 0 ; i < ed->mesh_vertices.size(); i++){ + for (unsigned int i = 0 ; i < ed->mesh_vertices.size(); i++){ myAssembler.fixVertex(ed->mesh_vertices[i],0, 1, value); } } @@ -20,7 +20,7 @@ static void fixEdgeToValue (GEdge *ed, static void fixEdgeToValueX (GEdge *ed, gmshAssembler &myAssembler) { - for (int i = 0 ; i < ed->lines.size(); i++){ + for (unsigned int i = 0 ; i < ed->lines.size(); i++){ myAssembler.fixVertex(ed->lines[i]->getVertex(0), 0, 1, ed->lines[i]->getVertex(0)->x()); myAssembler.fixVertex(ed->lines[i]->getVertex(1), 0, 1, ed->lines[i]->getVertex(1)->x()); } @@ -29,7 +29,7 @@ static void fixEdgeToValueX (GEdge *ed, static void fixEdgeToValueY (GEdge *ed, gmshAssembler &myAssembler) { - for (int i = 0 ; i < ed->lines.size(); i++){ + for (unsigned int i = 0 ; i < ed->lines.size(); i++){ myAssembler.fixVertex(ed->lines[i]->getVertex(0), 0, 1, ed->lines[i]->getVertex(0)->y()); myAssembler.fixVertex(ed->lines[i]->getVertex(1), 0, 1, ed->lines[i]->getVertex(1)->y()); } @@ -74,7 +74,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound, std::list<GEdge*> &V1) : GFace(m,tag),_compound(compound),_U0(U0),_U1(U1),_V0(V0),_V1(V1), oct(0) { - printf("%d %d %d %d \n",_U0.size(),_U1.size(),_V0.size(),_V1.size()); + //printf("%d %d %d %d \n",_U0.size(),_U1.size(),_V0.size(),_V1.size()); getBoundingEdges(); } @@ -130,7 +130,7 @@ void GFaceCompound::parametrize (bool _isU) const std::list<GFace*> :: const_iterator it = _compound.begin(); for ( ; it != _compound.end() ; ++it){ - for ( int i=0;i<(*it)->triangles.size() ; ++i){ + for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; myAssembler.numberVertex(t->getVertex(0),0, 1); myAssembler.numberVertex(t->getVertex(1),0, 1); @@ -148,10 +148,10 @@ void GFaceCompound::parametrize (bool _isU) const it = _compound.begin(); for ( ; it != _compound.end() ; ++it){ - for ( int i=0;i<(*it)->triangles.size() ; ++i){ + for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; - double uu[3],vv[3]; - for (int j=0;j<3;j++){ + double uu[3], vv[3]; + for (int j = 0; j < 3; j++){ MVertex *v = t->getVertex(j); double value = myAssembler.getDofValue(v,0,1); std::map<MVertex*,SPoint2>::iterator itf = coordinates.find(v); @@ -160,7 +160,7 @@ void GFaceCompound::parametrize (bool _isU) const coordinates[v] = p; } else{ - if(_isU)itf->second[0]= value; + if(_isU) itf->second[0]= value; else itf->second[1]= value; uu[j] = itf->second[0]; vv[j] = itf->second[1]; @@ -174,14 +174,14 @@ void GFaceCompound::parametrize (bool _isU) const GPoint GFaceCompound::point(double par1, double par2) const { parametrize(); - double U,V; + double U, V; MTriangle *t; getTriangle (par1, par2, &t, U,V); - SPoint3 p(0,0,0); - if (!t)return GPoint(p.x(),p.y(),p.z(),this); - t->pnt(U,V,0,p); - double par[2] = {par1,par2}; - return GPoint(p.x(),p.y(),p.z(),this,par); + SPoint3 p(0, 0, 0); + if (!t) return GPoint(p.x(), p.y(), p.z(), this); + t->pnt(U, V, 0, p); + double par[2] = {par1, par2}; + return GPoint(p.x(), p.y(), p.z(), this, par); } /* @@ -191,27 +191,27 @@ GPoint GFaceCompound::point(double par1, double par2) const Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const { parametrize(); - double U,V; + double U, V; MTriangle *t; - getTriangle (param.x(), param.y(), &t, U,V); + getTriangle(param.x(), param.y(), &t, U,V); double jac[3][3]; - t->getJacobian(U,V,0,jac); - return Pair<SVector3, SVector3>(SVector3(jac[0][0],jac[0][1],jac[0][2]), - SVector3(jac[1][0],jac[1][1],jac[1][2])); + t->getJacobian(U, V, 0, jac); + return Pair<SVector3, SVector3>(SVector3(jac[0][0], jac[0][1], jac[0][2]), + SVector3(jac[1][0], jac[1][1], jac[1][2])); } -static void GFaceCompoundBB(void *a, double*mmin, double*mmax) +static void GFaceCompoundBB(void *a, double*mmin, double*mmax) { GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a; - mmin[0] = std::min(std::min(t->p1.x(),t->p2.x()),t->p3.x()); - mmin[1] = std::min(std::min(t->p1.y(),t->p2.y()),t->p3.y()); - mmax[0] = std::max(std::max(t->p1.x(),t->p2.x()),t->p3.x()); - mmax[1] = std::max(std::max(t->p1.y(),t->p2.y()),t->p3.y()); + mmin[0] = std::min(std::min(t->p1.x(), t->p2.x()), t->p3.x()); + mmin[1] = std::min(std::min(t->p1.y(), t->p2.y()), t->p3.y()); + mmax[0] = std::max(std::max(t->p1.x(), t->p2.x()), t->p3.x()); + mmax[1] = std::max(std::max(t->p1.y(), t->p2.y()), t->p3.y()); mmin[2] = -1; mmax[2] = 1; } -static void GFaceCompoundCentroid(void *a, double*c) +static void GFaceCompoundCentroid(void *a, double*c) { GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a; c[0] = (t->p1.x() + t->p2.x() + t->p3.x())/3.0; @@ -219,7 +219,7 @@ static void GFaceCompoundCentroid(void *a, double*c) c[2] = 0.0; } -static int GFaceCompoundInEle(void *a, double*c) +static int GFaceCompoundInEle(void *a, double*c) { GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a; double M[2][2],R[2],X[2]; @@ -234,17 +234,15 @@ static int GFaceCompoundInEle(void *a, double*c) R[0] = (c[0] - p0.x()); R[1] = (c[1] - p0.y()); sys2x2(M,R,X); - if (X[0] > -eps && - X[1] > -eps && - 1.-X[0]-X[1] > -eps){ + if (X[0] > -eps && X[1] > -eps && 1. - X[0] - X[1] > -eps){ return 1; } return 0; } -void GFaceCompound::getTriangle (double u, double v, - MTriangle **t, - double &_u, double &_v) const +void GFaceCompound::getTriangle(double u, double v, + MTriangle **t, + double &_u, double &_v) const { parametrize(); @@ -254,7 +252,7 @@ void GFaceCompound::getTriangle (double u, double v, *t = tt->t; double M[2][2],X[2],R[2]; - const double eps = 1.e-6; +// const double eps = 1.e-6; const SPoint2 p0 = tt->p1; const SPoint2 p1 = tt->p2; const SPoint2 p2 = tt->p3; @@ -310,7 +308,7 @@ void GFaceCompound::buildOct() const int count = 0; std::list<GFace*> :: const_iterator it = _compound.begin(); for ( ; it != _compound.end() ; ++it){ - for ( int i=0;i<(*it)->triangles.size() ; ++i){ + for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; for (int j=0;j<3;j++){ std::map<MVertex*,SPoint2>::const_iterator itj = @@ -344,7 +342,7 @@ void GFaceCompound::buildOct() const fprintf(xyzv,"View \"\"{\n"); for ( ; it != _compound.end() ; ++it){ - for ( int i=0;i<(*it)->triangles.size() ; ++i){ + for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *t = (*it)->triangles[i]; std::map<MVertex*,SPoint2>::const_iterator it0 = coordinates.find(t->getVertex(0)); diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp index 80491f5b8b3be5885c02f1a7c9ca770ddce15d67..6608363e637864cf4dbbd2d3d0ccf6f33124abd4 100644 --- a/Mesh/HighOrder.cpp +++ b/Mesh/HighOrder.cpp @@ -937,8 +937,6 @@ void getFaceVertices(GRegion *gr, MElement *ele, face.getVertex(2), hoEdgeNodes,nPts+1); - double X(0),Y(0),Z(0); - for (int k=start;k<points.size1();k++) { double t1 = points(k,0); @@ -1612,8 +1610,8 @@ void localHarmonicMapping(GModel *gm, myAssembler.fixVertex ( n2 , 0 , 0 , -1.0); myAssembler.fixVertex ( n3 , 0 , 0 , 1.0); myAssembler.fixVertex ( n4 , 0 , 0 , 1.0); - for (int i=0;i<e1.size() ; i++) myAssembler.fixVertex ( e1[i] , 0 , 0, -1.0); - for (int i=0;i<e3.size() ; i++) myAssembler.fixVertex ( e3[i] , 0 , 0, 1.0); + for (unsigned int i = 0; i < e1.size(); i++) myAssembler.fixVertex(e1[i], 0, 0, -1.0); + for (unsigned int i = 0; i < e3.size(); i++) myAssembler.fixVertex(e3[i], 0, 0, 1.0); Laplace.addToMatrix(myAssembler,t1); Laplace.addToMatrix(myAssembler,t2); lsys->systemSolve(); @@ -1626,8 +1624,8 @@ void localHarmonicMapping(GModel *gm, myAssembler1.fixVertex ( n3 , 0 , 1 , -1.0); myAssembler1.fixVertex ( n4 , 0 , 1 , 1.0); myAssembler1.fixVertex ( n1 , 0 , 1 , 1.0); - for (int i=0;i<e2.size() ; i++) myAssembler1.fixVertex ( e2[i] , 0 , 1, -1.0); - for (int i=0;i<e4.size() ; i++) myAssembler1.fixVertex ( e4[i] , 0 , 1, 1.0); + for (unsigned int i = 0; i < e2.size(); i++) myAssembler1.fixVertex(e2[i], 0, 1, -1.0); + for (unsigned int i = 0; i < e4.size(); i++) myAssembler1.fixVertex(e4[i], 0, 1, 1.0); Laplace1.addToMatrix(myAssembler1,t1); Laplace1.addToMatrix(myAssembler1,t2); lsys1->systemSolve(); @@ -1638,7 +1636,7 @@ void localHarmonicMapping(GModel *gm, // this can be done by evaluating the - for (int i=0;i<e.size() ; i++){ + for (unsigned int i = 0; i < e.size(); i++){ MVertex *v = e[i]; const double U = myAssembler.getDofValue (v, 0 ,0); const double V = myAssembler1.getDofValue (v, 0 ,1); @@ -1660,7 +1658,7 @@ void getParametricCoordnates ( GFace *gf, std::vector<MVertex*> &e, std::vector<SPoint2> ¶m){ param.clear(); - for (int i=0;i<e.size();i++){ + for (unsigned int i = 0; i < e.size(); i++){ double U,V; parametricCoordinates(e[i] , gf, U, V); param.push_back(SPoint2(U,V)); diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp index acbd09cc23e3c119957525061f615771aa39df82..c6d68fb56dde7538fa2df5289d29c3d96715e552 100644 --- a/Mesh/gmshSmoothHighOrder.cpp +++ b/Mesh/gmshSmoothHighOrder.cpp @@ -14,11 +14,11 @@ void getDistordedElements ( const std::vector<MElement*> & v, double &minD){ d.clear(); minD = 1; - for (int i=0;i<v.size() ; i++){ + for (unsigned int i = 0; i < v.size(); i++){ const double disto = v[i]->distoShapeMeasure(); - if ( disto < threshold) + if (disto < threshold) d.push_back(v[i]); - minD = std::min(minD,disto); + minD = std::min(minD, disto); } } @@ -26,10 +26,10 @@ void addOneLayer ( const std::vector<MElement*> & v, std::vector<MElement*> & d , std::vector<MElement*> & layer ){ std::set<MVertex*> all; - for (int i=0;i<d.size() ; i++){ + for (unsigned int i = 0; i < d.size(); i++){ MElement *e = d[i]; int n = e->getNumPrimaryVertices(); - for (int j=0;j<n;j++){ + for (int j = 0; j < n; j++){ all.insert(e->getVertex(j)); } } @@ -38,13 +38,13 @@ void addOneLayer ( const std::vector<MElement*> & v, std::sort(d.begin(), d.end()); - for (int i=0;i<v.size() ; i++){ + for (unsigned int i = 0; i < v.size(); i++){ MElement *e = v[i]; bool found = std::binary_search(d.begin(), d.end(), e); // element is not yet there if (!found){ int n = e->getNumPrimaryVertices(); - for (int j=0;j<n;j++){ + for (int j = 0; j < n; j++){ MVertex *vert = e->getVertex(j); if (all.find(vert) != all.end()){ layer.push_back(e); @@ -90,7 +90,7 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all) { // printf("%d elements / %d distorted min Disto = %g\n",all.size(),v.size(), minD); - if (!v.size())return; + if (!v.size()) return; const int nbLayers = 2; for (int i=0;i<nbLayers;i++){ @@ -99,19 +99,19 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all) { } // 3 -> .4 - printf("%d elements after adding %d layers\n",v.size(),nbLayers); + printf("%d elements after adding %d layers\n", v.size(), nbLayers); - addOneLayer ( all, v, layer); + addOneLayer(all, v, layer); // printf("%d elements in the next layer\n",layer.size()); - for (int i=0;i<layer.size() ; i++){ - for (int j=0;j<layer[i]->getNumVertices(); j++){ + for (unsigned int i = 0; i < layer.size(); i++){ + for (int j = 0; j < layer[i]->getNumVertices(); j++){ MVertex *vert = layer[i]->getVertex(j); - myAssembler.fixVertex ( vert , 0 , getTag() , 0); - myAssembler.fixVertex ( vert , 1 , getTag() , 0); - myAssembler.fixVertex ( vert , 2 , getTag() , 0); + myAssembler.fixVertex(vert, 0, getTag(), 0); + myAssembler.fixVertex(vert, 1, getTag(), 0); + myAssembler.fixVertex(vert, 2, getTag(), 0); } } @@ -121,8 +121,8 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all) { // printf("%d vertices \n", _displ.size()); - for (int i=0;i<v.size() ; i++){ - for (int j=0;j<v[i]->getNumVertices(); j++){ + for (unsigned int i = 0; i < v.size(); i++){ + for (int j = 0; j < v[i]->getNumVertices(); j++){ MVertex *vert = v[i]->getVertex(j); // printf("%d %d %d v\n",i,j,v[i]->getNumVertices()); it = _displ.find(vert); @@ -151,12 +151,12 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all) { } // number the other DOFs - for (int i=0;i<v.size() ; i++){ - for (int j=0;j<v[i]->getNumVertices(); j++){ + for (unsigned int i = 0; i < v.size(); i++){ + for (int j = 0; j < v[i]->getNumVertices(); j++){ MVertex *vert = v[i]->getVertex(j); - myAssembler.numberVertex ( vert , 0 , getTag() ); - myAssembler.numberVertex ( vert , 1 , getTag() ); - myAssembler.numberVertex ( vert , 2 , getTag() ); + myAssembler.numberVertex(vert, 0, getTag()); + myAssembler.numberVertex(vert, 1, getTag()); + myAssembler.numberVertex(vert, 2, getTag()); // gather all vertices that are supposed to move verticesToMove[vert] = SVector3(0.0,0.0,0.0); } diff --git a/Numeric/gmshAssembler.h b/Numeric/gmshAssembler.h index 529506e92ec98e7ec3651312b246e4410f3c8510..18d769affb4b4a8d268520835c58bb374ef846ec 100644 --- a/Numeric/gmshAssembler.h +++ b/Numeric/gmshAssembler.h @@ -35,23 +35,23 @@ class gmshAssembler { public: gmshAssembler (gmshLinearSystem *l) : lsys(l){} inline void constraintVertex(MVertex*v, int iComp, int iField, - std::vector<MVertex *>&verts, + std::vector<MVertex*> &verts, std::vector<double> &coeffs) { - std::vector<std::pair<gmshDofKey,double> > constraint; - gmshDofKey key (v,iComp,iField); - for (int i=0;i<verts.size();i++){ + std::vector<std::pair<gmshDofKey, double> > constraint; + gmshDofKey key(v, iComp, iField); + for (unsigned int i = 0; i < verts.size(); i++){ gmshDofKey key2 (verts[i],iComp,iField); - constraint.push_back(std::make_pair(key2,coeffs[i])); + constraint.push_back(std::make_pair(key2, coeffs[i])); } constraints[key] = constraint; } inline void numberVertex(MVertex*v, int iComp, int iField) { - gmshDofKey key (v,iComp,iField); - if (fixed.find(key) != fixed.end())return; - if (constraints.find(key) != constraints.end())return; - std::map<gmshDofKey, int> :: iterator it = numbering.find(key); + gmshDofKey key(v, iComp, iField); + if (fixed.find(key) != fixed.end()) return; + if (constraints.find(key) != constraints.end()) return; + std::map<gmshDofKey, int>::iterator it = numbering.find(key); if (it == numbering.end()){ int size = numbering.size(); numbering[key] = size; @@ -59,41 +59,41 @@ public: } inline void numberVertex(MVertex*v, int iComp, int iField, int iField2) { - gmshDofKey key (v,iComp,iField); - gmshDofKey key2 (v,iComp,iField2); - if (fixed.find(key) != fixed.end())return; - if (fixed.find(key2) != fixed.end())return; - if (constraints.find(key) != constraints.end())return; + gmshDofKey key(v, iComp, iField); + gmshDofKey key2(v, iComp, iField2); + if (fixed.find(key) != fixed.end()) return; + if (fixed.find(key2) != fixed.end()) return; + if (constraints.find(key) != constraints.end()) return; std::map<gmshDofKey, int> :: iterator it = numbering.find(key); if (it == numbering.end()){ int size = numbering.size(); numbering[key] = size; } } - inline void fixVertex (MVertex*v, int iComp, int iField, double val) + inline void fixVertex(MVertex*v, int iComp, int iField, double val) { - fixed[gmshDofKey(v,iComp,iField)] = val; + fixed[gmshDofKey(v, iComp, iField)] = val; } - inline double getDofValue (MVertex*v, int iComp, int iField) const + inline double 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); - if (it != fixed.end())return it->second; + std::map<gmshDofKey, double>::const_iterator it = fixed.find(key); + if (it != fixed.end()) return it->second; } { - std::map<gmshDofKey, int> :: const_iterator it = numbering.find(key); + std::map<gmshDofKey, int>::const_iterator it = numbering.find(key); if (it != numbering.end()) - return lsys->getFromSolution (it->second); + return lsys->getFromSolution(it->second); } { std::map<gmshDofKey, std::vector<std::pair<gmshDofKey,double> > >:: const_iterator itConstr = constraints.find(key); if (itConstr != constraints.end()){ double val = 0; - for (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; - double valConstr = itConstr->second[i].second; + double valConstr = itConstr->second[i].second; val += getDofValue(dofKeyConstr.v, dofKeyConstr.comp, dofKeyConstr.field) * valConstr; } diff --git a/Numeric/gmshFunction.h b/Numeric/gmshFunction.h index b11288f3adc2a07074e00c57138fca9472bfa340..a3187794863d04c4eeaf0c7b376850a3b70a6d0a 100644 --- a/Numeric/gmshFunction.h +++ b/Numeric/gmshFunction.h @@ -5,6 +5,7 @@ class gmshFunction { double _val; public : gmshFunction(double val = 0) : _val(val) {} + virtual ~gmshFunction(){} virtual double operator () (double x, double y, double z) const { return _val; } }; diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h index 13eff4cd1f11d6a4d0c7b29adca23c7bdc2aac55..bd830cc0e558bbcde25adf4eaa7eea285629930e 100644 --- a/Numeric/gmshLinearSystemGmm.h +++ b/Numeric/gmshLinearSystemGmm.h @@ -11,11 +11,10 @@ #include <gmm.h> class gmshLinearSystemGmm : public gmshLinearSystem { - gmm::row_matrix< gmm::wsvector<double> > *_a; - std::vector<double> *_x; - std::vector<double> *_b; + gmm::row_matrix<gmm::wsvector<double> > *_a; + std::vector<double> *_b, *_x; public : - gmshLinearSystemGmm () : _a(0),_b(0),_x(0) {} + gmshLinearSystemGmm () : _a(0), _b(0), _x(0) {} virtual bool isAllocated () const {return _a != 0;} virtual void allocate (int _nbRows) { @@ -58,7 +57,7 @@ public : } virtual void zeroRightHandSide () { - for (int i=0;i<_b->size();i++)(*_b)[i] = 0; + for (unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0; } virtual int systemSolve () { @@ -68,6 +67,7 @@ public : //iter.set_noisy(2); //gmm::gmres(*_a, *_x, *_b, P, 100, iter); // execute the GMRES algorithm gmm::cg(*_a, *_x, *_b, P, iter); // execute the CG algorithm + return 1; } }; diff --git a/Numeric/gmshTermOfFormulation.cpp b/Numeric/gmshTermOfFormulation.cpp index 5a8f4d14f9aae75802af6b9b7aca69cbfbc3475f..12df17436f88362a011c95151376c28b4bf28780 100644 --- a/Numeric/gmshTermOfFormulation.cpp +++ b/Numeric/gmshTermOfFormulation.cpp @@ -47,8 +47,8 @@ void gmshNodalFemTerm::addToMatrix (gmshAssembler &lsys) const void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys,const std::vector<MElement*> &v) const { - for (int i=0;i<v.size();i++) - addToMatrix (lsys,v[i]); + for (unsigned int i = 0; i < v.size(); i++) + addToMatrix(lsys, v[i]); } void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys, diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h index 2ddfa7e6e41d629b783500833273c2a7ac04b92b..689d6ffae206918d9a243d1aaced2174e1498db9 100644 --- a/Numeric/gmshTermOfFormulation.h +++ b/Numeric/gmshTermOfFormulation.h @@ -20,6 +20,7 @@ protected: GModel *_gm; public: gmshTermOfFormulation(GModel *gm) : _gm(gm) {} + virtual ~gmshTermOfFormulation(){} virtual void addToMatrix(gmshAssembler&) const = 0; virtual void addToRightHandSide(gmshAssembler&) const = 0; };