From cb6804905a716dda81fdbbdbe7c5100f27da7bdc Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 13 Nov 2009 13:17:03 +0000
Subject: [PATCH] - simplified dofManager using traits (allows to remove second
template argument) - renamed functionSpace into polynomialBasis
(functionSpace will be much higher level)
---
Geo/GFaceCompound.cpp | 10 ++--
Geo/MElement.cpp | 4 +-
Geo/MElement.h | 4 +-
Geo/MElementCut.h | 8 +--
Geo/MLine.cpp | 12 ++--
Geo/MLine.h | 2 +-
Geo/MQuadrangle.cpp | 22 ++++----
Geo/MQuadrangle.h | 2 +-
Geo/MTetrahedron.cpp | 22 ++++----
Geo/MTetrahedron.h | 2 +-
Geo/MTriangle.cpp | 22 ++++----
Geo/MTriangle.h | 2 +-
Mesh/HighOrder.cpp | 12 ++--
Mesh/highOrderSmoother.cpp | 6 +-
Mesh/highOrderSmoother.h | 2 +-
Mesh/qualityMeasures.cpp | 2 +-
Numeric/CMakeLists.txt | 2 +-
...{functionSpace.cpp => polynomialBasis.cpp} | 56 +++++++++----------
.../{functionSpace.h => polynomialBasis.h} | 28 +++++-----
Post/PViewDataGModel.cpp | 2 +-
Post/PViewDataList.cpp | 4 +-
Solver/convexCombinationTerm.h | 4 +-
Solver/crossConfTerm.h | 4 +-
Solver/dgGroupOfElements.cpp | 6 +-
Solver/dgGroupOfElements.h | 8 +--
Solver/dofManager.h | 32 +++++++++--
Solver/eigenSolver.cpp | 2 +-
Solver/eigenSolver.h | 4 +-
Solver/elasticitySolver.cpp | 4 +-
Solver/elasticitySolver.h | 2 +-
Solver/elasticityTerm.h | 6 +-
Solver/femTerm.h | 17 +++---
Solver/helmholtzTerm.h | 4 +-
Solver/multiscaleLaplace.cpp | 2 +-
34 files changed, 174 insertions(+), 147 deletions(-)
rename Numeric/{functionSpace.cpp => polynomialBasis.cpp} (94%)
rename Numeric/{functionSpace.h => polynomialBasis.h} (83%)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 973133fa95..6afed65789 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -16,7 +16,7 @@
#include "Octree.h"
#include "SBoundingBox3d.h"
#include "SPoint3.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
#include "robustPredicates.h"
#include "MElementCut.h"
#include "GEntity.h"
@@ -34,7 +34,7 @@
#include "discreteFace.h"
#include "eigenSolver.h"
-static void fixEdgeToValue(GEdge *ed, double value, dofManager<double, double> &myAssembler)
+static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
{
myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value);
myAssembler.fixVertex(ed->getEndVertex()->mesh_vertices[0], 0, 1, value);
@@ -869,7 +869,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
{
- dofManager<double, double> myAssembler(_lsys);
+ dofManager<double> myAssembler(_lsys);
simpleFunction<double> ONE(1.0);
if(_type == UNITCIRCLE){
@@ -976,7 +976,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
void GFaceCompound::parametrize_conformal() const
{
- dofManager<double, double> myAssembler(_lsys);
+ dofManager<double> myAssembler(_lsys);
std::vector<MVertex*> ordered;
std::vector<double> coords;
@@ -1059,7 +1059,7 @@ void GFaceCompound::compute_distance() const
double L = norm(SVector3(bbox.max(), bbox.min()));
double mu = L/28;
simpleFunction<double> DIFF(mu * mu), MONE(1.0);
- dofManager<double, double> myAssembler(_lsys);
+ dofManager<double> myAssembler(_lsys);
distanceTerm distance(model(), 1, &DIFF, &MONE);
std::vector<MVertex*> ordered;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 01b7398403..ca3ba8a384 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -110,7 +110,7 @@ double MElement::rhoShapeMeasure()
void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
{
- const functionSpace* fs = getFunctionSpace(o);
+ const polynomialBasis* fs = getFunctionSpace(o);
if(fs) fs->f(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
@@ -118,7 +118,7 @@ void MElement::getShapeFunctions(double u, double v, double w, double s[], int o
void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],
int o)
{
- const functionSpace* fs = getFunctionSpace(o);
+ const polynomialBasis* fs = getFunctionSpace(o);
if(fs) fs->df(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index a35759b733..44f6c2b73c 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -14,7 +14,7 @@
#include "MVertex.h"
#include "MEdge.h"
#include "MFace.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
#include "Gauss.h"
class GFace;
@@ -173,7 +173,7 @@ class MElement
virtual std::string getInfoString();
// get the function space for the element
- virtual const functionSpace* getFunctionSpace(int order=-1) const { return 0; }
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
// return the interpolating nodal shape functions evaluated at point
// (u,v,w) in parametric coordinates (if order == -1, use the
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index fa37382c5a..a9f355c15d 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -105,7 +105,7 @@ class MPolyhedron : public MElement {
return vol;
}
virtual int getVolumeSign() { return (getVolume() >= 0) ? 1 : -1; }
- virtual const functionSpace* getFunctionSpace(int order=-1) const
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const
{
return _orig->getFunctionSpace(order);
}
@@ -207,7 +207,7 @@ class MPolygon : public MElement {
}
virtual int getType() const { return TYPE_POLYG; }
virtual int getTypeForMSH() const { return MSH_POLYG_; }
- virtual const functionSpace* getFunctionSpace(int order=-1) const
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const
{
return _orig->getFunctionSpace(order);
}
@@ -251,7 +251,7 @@ class MTriangleBorder : public MTriangle {
if(_domains[1]) return _domains[1]->getParent();
return NULL;
}
- virtual const functionSpace* getFunctionSpace(int order=-1) const
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const
{
return getParent()->getFunctionSpace(order);
}
@@ -287,7 +287,7 @@ class MLineBorder : public MLine {
if(_domains[1]) return _domains[1]->getParent();
return NULL;
}
- virtual const functionSpace* getFunctionSpace(int order=-1) const
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const
{
return getParent()->getFunctionSpace(order);
}
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index cdd91eb117..a4f906ffa2 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -8,16 +8,16 @@
#include "Context.h"
#include "qualityMeasures.h"
-const functionSpace* MLine::getFunctionSpace(int o) const
+const polynomialBasis* MLine::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 1: return &functionSpaces::find(MSH_LIN_2);
- case 2: return &functionSpaces::find(MSH_LIN_3);
- case 3: return &functionSpaces::find(MSH_LIN_4);
- case 4: return &functionSpaces::find(MSH_LIN_5);
- case 5: return &functionSpaces::find(MSH_LIN_6);
+ case 1: return &polynomialBases::find(MSH_LIN_2);
+ case 2: return &polynomialBases::find(MSH_LIN_3);
+ case 3: return &polynomialBases::find(MSH_LIN_4);
+ case 4: return &polynomialBases::find(MSH_LIN_5);
+ case 5: return &polynomialBases::find(MSH_LIN_6);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index c78fe34652..74542824aa 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -63,7 +63,7 @@ class MLine : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
- virtual const functionSpace* getFunctionSpace(int o=-1) const;
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 3daf1324c1..aec951b0c4 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -8,7 +8,7 @@
#include "Context.h"
#include "qualityMeasures.h"
#define SQU(a) ((a)*(a))
-const functionSpace* MQuadrangle::getFunctionSpace(int o) const
+const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
@@ -16,19 +16,19 @@ const functionSpace* MQuadrangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 1: return &functionSpaces::find(MSH_QUA_4);
- case 2: return &functionSpaces::find(MSH_QUA_8);
- case 3: return &functionSpaces::find(MSH_QUA_12);
- case 4: return &functionSpaces::find(MSH_QUA_16I);
- case 5: return &functionSpaces::find(MSH_QUA_20);
+ case 1: return &polynomialBases::find(MSH_QUA_4);
+ case 2: return &polynomialBases::find(MSH_QUA_8);
+ case 3: return &polynomialBases::find(MSH_QUA_12);
+ case 4: return &polynomialBases::find(MSH_QUA_16I);
+ case 5: return &polynomialBases::find(MSH_QUA_20);
}
}
switch (order) {
- case 1: return &functionSpaces::find(MSH_QUA_4);
- case 2: return &functionSpaces::find(MSH_QUA_9);
- case 3: return &functionSpaces::find(MSH_QUA_16);
- case 4: return &functionSpaces::find(MSH_QUA_25);
- case 5: return &functionSpaces::find(MSH_QUA_36);
+ case 1: return &polynomialBases::find(MSH_QUA_4);
+ case 2: return &polynomialBases::find(MSH_QUA_9);
+ case 3: return &polynomialBases::find(MSH_QUA_16);
+ case 4: return &polynomialBases::find(MSH_QUA_25);
+ case 5: return &polynomialBases::find(MSH_QUA_36);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
return 0;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index d6457574b0..e4e5cbb373 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -115,7 +115,7 @@ class MQuadrangle : public MElement {
virtual const char *getStringForPOS() const { return "SQ"; }
virtual const char *getStringForBDF() const { return "CQUAD4"; }
virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
- virtual const functionSpace* getFunctionSpace(int o=-1) const;
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 9875212bff..e79da9bb50 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -83,7 +83,7 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3])
sys3x3(mat, b, uvw, &det);
}
-const functionSpace* MTetrahedron::getFunctionSpace(int o) const
+const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
@@ -91,21 +91,21 @@ const functionSpace* MTetrahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return &functionSpaces::find(MSH_TET_4);
- case 2: return &functionSpaces::find(MSH_TET_10);
- case 3: return &functionSpaces::find(MSH_TET_20);
- case 4: return &functionSpaces::find(MSH_TET_34);
- case 5: return &functionSpaces::find(MSH_TET_52);
+ case 1: return &polynomialBases::find(MSH_TET_4);
+ case 2: return &polynomialBases::find(MSH_TET_10);
+ case 3: return &polynomialBases::find(MSH_TET_20);
+ case 4: return &polynomialBases::find(MSH_TET_34);
+ case 5: return &polynomialBases::find(MSH_TET_52);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return &functionSpaces::find(MSH_TET_4);
- case 2: return &functionSpaces::find(MSH_TET_10);
- case 3: return &functionSpaces::find(MSH_TET_20);
- case 4: return &functionSpaces::find(MSH_TET_35);
- case 5: return &functionSpaces::find(MSH_TET_56);
+ case 1: return &polynomialBases::find(MSH_TET_4);
+ case 2: return &polynomialBases::find(MSH_TET_10);
+ case 3: return &polynomialBases::find(MSH_TET_20);
+ case 4: return &polynomialBases::find(MSH_TET_35);
+ case 5: return &polynomialBases::find(MSH_TET_56);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index dec9e6e20a..bc98392d9a 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -130,7 +130,7 @@ class MTetrahedron : public MElement {
virtual double distoShapeMeasure();
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
- virtual const functionSpace* getFunctionSpace(int o=-1) const;
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index cd97bce754..45b1b501c7 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -41,7 +41,7 @@ double MTriangle::gammaShapeMeasure()
#endif
}
-const functionSpace* MTriangle::getFunctionSpace(int o) const
+const polynomialBasis* MTriangle::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
@@ -49,21 +49,21 @@ const functionSpace* MTriangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 1: return &functionSpaces::find(MSH_TRI_3);
- case 2: return &functionSpaces::find(MSH_TRI_6);
- case 3: return &functionSpaces::find(MSH_TRI_9);
- case 4: return &functionSpaces::find(MSH_TRI_12);
- case 5: return &functionSpaces::find(MSH_TRI_15I);
+ case 1: return &polynomialBases::find(MSH_TRI_3);
+ case 2: return &polynomialBases::find(MSH_TRI_6);
+ case 3: return &polynomialBases::find(MSH_TRI_9);
+ case 4: return &polynomialBases::find(MSH_TRI_12);
+ case 5: return &polynomialBases::find(MSH_TRI_15I);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return &functionSpaces::find(MSH_TRI_3);
- case 2: return &functionSpaces::find(MSH_TRI_6);
- case 3: return &functionSpaces::find(MSH_TRI_10);
- case 4: return &functionSpaces::find(MSH_TRI_15);
- case 5: return &functionSpaces::find(MSH_TRI_21);
+ case 1: return &polynomialBases::find(MSH_TRI_3);
+ case 2: return &polynomialBases::find(MSH_TRI_6);
+ case 3: return &polynomialBases::find(MSH_TRI_10);
+ case 4: return &polynomialBases::find(MSH_TRI_15);
+ case 5: return &polynomialBases::find(MSH_TRI_21);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index e59c1b159a..67172ff358 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -121,7 +121,7 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
- virtual const functionSpace* getFunctionSpace(int o=-1) const;
+ virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 6c361c397f..eecf866a62 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -21,7 +21,7 @@
#include "Numeric.h"
#include "Context.h"
#include "fullMatrix.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
#define SQU(a) ((a)*(a))
@@ -521,15 +521,15 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
switch (nPts){
case 2 :
- points = functionSpaces::find(MSH_TRI_10).points;
+ points = polynomialBases::find(MSH_TRI_10).points;
start = 9;
break;
case 3 :
- points = functionSpaces::find(MSH_TRI_15).points;
+ points = polynomialBases::find(MSH_TRI_15).points;
start = 12;
break;
case 4 :
- points = functionSpaces::find(MSH_TRI_21).points;
+ points = polynomialBases::find(MSH_TRI_21).points;
start = 15;
break;
default :
@@ -614,11 +614,11 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
switch (nPts){
case 3 :
- points = functionSpaces::find(MSH_TET_35).points;
+ points = polynomialBases::find(MSH_TET_35).points;
start = 34;
break;
case 4 :
- points = functionSpaces::find(MSH_TET_56).points;
+ points = polynomialBases::find(MSH_TET_56).points;
start = 52;
break;
default :
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index dc56882ddd..e6e5168e53 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -434,7 +434,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
lsys->setGmres(1);
lsys->setPrec(5.e-8);
#endif
- dofManager<double,double> myAssembler(lsys);
+ dofManager<double> myAssembler(lsys);
elasticityTerm El(0, 1.0, .333, getTag());
std::vector<MElement*> layer, v;
@@ -540,7 +540,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
GFace *gf,
- dofManager<double,double> &myAssembler,
+ dofManager<double> &myAssembler,
std::set<MVertex*> &verticesToMove,
elasticityTerm &El)
{
@@ -609,7 +609,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
lsys->setGmres(1);
lsys->setPrec(5.e-8);
#endif
- dofManager<double,double> myAssembler(lsys);
+ dofManager<double> myAssembler(lsys);
elasticityTerm El(0, 1.0, .333, getTag());
std::vector<MElement*> layer, v;
diff --git a/Mesh/highOrderSmoother.h b/Mesh/highOrderSmoother.h
index 8b16f060d2..3b4379fbc0 100644
--- a/Mesh/highOrderSmoother.h
+++ b/Mesh/highOrderSmoother.h
@@ -44,7 +44,7 @@ public:
}
void smooth(std::vector<MElement*> & );
double smooth_metric_(std::vector<MElement*> &, GFace *gf,
- dofManager<double,double> &myAssembler,
+ dofManager<double> &myAssembler,
std::set<MVertex*> &verticesToMove,
elasticityTerm &El);
void smooth_metric(std::vector<MElement*> &, GFace *gf );
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index c35eb72848..a4bc0b5b97 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -9,7 +9,7 @@
#include "MTriangle.h"
#include "MTetrahedron.h"
#include "Numeric.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
#include "GmshMessage.h"
double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index bd8066edd2..e2e2451b4d 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -6,7 +6,7 @@
set(SRC
Numeric.cpp
fullMatrix.cpp
- functionSpace.cpp
+ polynomialBasis.cpp
GaussQuadratureLin.cpp
GaussQuadratureTri.cpp
GaussQuadratureQuad.cpp
diff --git a/Numeric/functionSpace.cpp b/Numeric/polynomialBasis.cpp
similarity index 94%
rename from Numeric/functionSpace.cpp
rename to Numeric/polynomialBasis.cpp
index fa2f604eb4..6a82c3809b 100644
--- a/Numeric/functionSpace.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -9,7 +9,7 @@
#include "GmshDefines.h"
#include "GmshMessage.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
static fullMatrix<double> generate1DMonomials(int order)
{
@@ -733,7 +733,7 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
}
}
-static void generate3dFaceClosure(functionSpace::clCont &closure, int order)
+static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
{
for (int iRotate = 0; iRotate < 3; iRotate++){
for (int iSign = 1; iSign >= -1; iSign -= 2){
@@ -746,7 +746,7 @@ static void generate3dFaceClosure(functionSpace::clCont &closure, int order)
}
}
-static void generate2dEdgeClosure(functionSpace::clCont &closure, int order)
+static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order)
{
closure.clear();
// first give edge nodes of the three edges in direct order
@@ -764,14 +764,14 @@ static void generate2dEdgeClosure(functionSpace::clCont &closure, int order)
closure.push_back(c3); closure.push_back(c4); closure.push_back(c5);
}
-std::map<int, functionSpace> functionSpaces::fs;
+std::map<int, polynomialBasis> polynomialBases::fs;
-const functionSpace &functionSpaces::find(int tag)
+const polynomialBasis &polynomialBases::find(int tag)
{
- std::map<int, functionSpace>::const_iterator it = fs.find(tag);
+ std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
if (it != fs.end()) return it->second;
- functionSpace F;
+ polynomialBasis F;
switch (tag){
case MSH_PNT:
@@ -801,77 +801,77 @@ const functionSpace &functionSpaces::find(int tag)
case MSH_TRI_3 :
F.monomials = generatePascalTriangle(1);
F.points = gmshGeneratePointsTriangle(1, false);
- generate2dEdgeClosure (F.edgeClosure,1);
+ generate2dEdgeClosure(F.edgeClosure, 1);
break;
case MSH_TRI_6 :
F.monomials = generatePascalTriangle(2);
F.points = gmshGeneratePointsTriangle(2, false);
- generate2dEdgeClosure (F.edgeClosure,2);
+ generate2dEdgeClosure(F.edgeClosure, 2);
break;
case MSH_TRI_9 :
F.monomials = generatePascalSerendipityTriangle(3);
F.points = gmshGeneratePointsTriangle(3, true);
- generate2dEdgeClosure (F.edgeClosure,3);
+ generate2dEdgeClosure(F.edgeClosure, 3);
break;
case MSH_TRI_10 :
F.monomials = generatePascalTriangle(3);
F.points = gmshGeneratePointsTriangle(3, false);
- generate2dEdgeClosure (F.edgeClosure,3);
+ generate2dEdgeClosure(F.edgeClosure, 3);
break;
case MSH_TRI_12 :
F.monomials = generatePascalSerendipityTriangle(4);
F.points = gmshGeneratePointsTriangle(4, true);
- generate2dEdgeClosure (F.edgeClosure,4);
+ generate2dEdgeClosure(F.edgeClosure, 4);
break;
case MSH_TRI_15 :
F.monomials = generatePascalTriangle(4);
F.points = gmshGeneratePointsTriangle(4, false);
- generate2dEdgeClosure (F.edgeClosure,4);
+ generate2dEdgeClosure(F.edgeClosure, 4);
break;
case MSH_TRI_15I :
F.monomials = generatePascalSerendipityTriangle(5);
F.points = gmshGeneratePointsTriangle(5, true);
- generate2dEdgeClosure (F.edgeClosure,5);
+ generate2dEdgeClosure(F.edgeClosure, 5);
break;
case MSH_TRI_21 :
F.monomials = generatePascalTriangle(5);
F.points = gmshGeneratePointsTriangle(5, false);
- generate2dEdgeClosure (F.edgeClosure,5);
+ generate2dEdgeClosure(F.edgeClosure, 5);
break;
case MSH_TET_4 :
F.monomials = generatePascalTetrahedron(1);
F.points = gmshGeneratePointsTetrahedron(1, false);
- generate3dFaceClosure (F.faceClosure,1);
+ generate3dFaceClosure(F.faceClosure, 1);
break;
case MSH_TET_10 :
F.monomials = generatePascalTetrahedron(2);
F.points = gmshGeneratePointsTetrahedron(2, false);
- generate3dFaceClosure (F.faceClosure,2);
+ generate3dFaceClosure(F.faceClosure, 2);
break;
case MSH_TET_20 :
F.monomials = generatePascalTetrahedron(3);
F.points = gmshGeneratePointsTetrahedron(3, false);
- generate3dFaceClosure (F.faceClosure,3);
+ generate3dFaceClosure(F.faceClosure, 3);
break;
case MSH_TET_35 :
F.monomials = generatePascalTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, false);
- generate3dFaceClosure (F.faceClosure,4);
+ generate3dFaceClosure(F.faceClosure, 4);
break;
case MSH_TET_34 :
F.monomials = generatePascalSerendipityTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, true);
- generate3dFaceClosure (F.faceClosure,4);
+ generate3dFaceClosure(F.faceClosure, 4);
break;
case MSH_TET_52 :
F.monomials = generatePascalSerendipityTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, true);
- generate3dFaceClosure (F.faceClosure,5);
+ generate3dFaceClosure(F.faceClosure, 5);
break;
case MSH_TET_56 :
F.monomials = generatePascalTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, false);
- generate3dFaceClosure (F.faceClosure,5);
+ generate3dFaceClosure(F.faceClosure, 5);
break;
case MSH_QUA_4 :
F.monomials = generatePascalQuad(1);
@@ -913,7 +913,7 @@ const functionSpace &functionSpaces::find(int tag)
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.monomials = generatePascalTetrahedron(1);
F.points = gmshGeneratePointsTetrahedron(1, false);
- generate3dFaceClosure (F.faceClosure,1);
+ generate3dFaceClosure(F.faceClosure, 1);
break;
}
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
@@ -921,18 +921,18 @@ const functionSpace &functionSpaces::find(int tag)
return fs[tag];
}
-std::map<std::pair<int, int>, fullMatrix<double> > functionSpaces::injector;
+std::map<std::pair<int, int>, fullMatrix<double> > polynomialBases::injector;
-const fullMatrix<double> &functionSpaces::findInjector(int tag1, int tag2)
+const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
{
std::pair<int,int> key(tag1,tag2);
std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it = injector.find(key);
if (it != injector.end()) return it->second;
- const functionSpace& fs1 = find(tag1);
- const functionSpace& fs2 = find(tag2);
+ const polynomialBasis& fs1 = find(tag1);
+ const polynomialBasis& fs2 = find(tag2);
- fullMatrix<double> inj(fs1.points.size1(),fs2.points.size1());
+ fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
double sf[256];
diff --git a/Numeric/functionSpace.h b/Numeric/polynomialBasis.h
similarity index 83%
rename from Numeric/functionSpace.h
rename to Numeric/polynomialBasis.h
index 6f8ed6be74..01b8f846c3 100644
--- a/Numeric/functionSpace.h
+++ b/Numeric/polynomialBasis.h
@@ -3,19 +3,19 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef _FUNCTION_SPACE_H_
-#define _FUNCTION_SPACE_H_
+#ifndef _POLYNOMIAL_BASIS_H_
+#define _POLYNOMIAL_BASIS_H_
#include <math.h>
#include <map>
#include <vector>
#include "fullMatrix.h"
-// presently thos function spaces are only for simplices
-// should be extended to other elements like quads and hexes
-struct functionSpace
+// presently those function spaces are only for simplices and quads;
+// should be extended to other elements like hexes
+struct polynomialBasis
{
- typedef std::vector<std::vector<int> > clCont;
+ typedef std::vector<std::vector<int> > clCont;
clCont faceClosure;
clCont edgeClosure;
fullMatrix<double> points;
@@ -23,11 +23,13 @@ struct functionSpace
fullMatrix<double> coefficients;
// for a given face/edge, with both a sign and a rotation,
// give an ordered list of nodes on this face/edge
- inline const std::vector<int> & getFaceClosure (int iFace, int iSign, int iRot)const{
- return faceClosure[iFace+4*(iSign==1?0:1)+8*iRot];
+ inline const std::vector<int> &getFaceClosure(int iFace, int iSign, int iRot) const
+ {
+ return faceClosure[iFace + 4 * (iSign == 1 ? 0 : 1) + 8 * iRot];
}
- inline const std::vector<int> & getEdgeClosure (int iEdge, int iSign) const{
- return edgeClosure[iSign == 1 ? iEdge : 3+iEdge];
+ inline const std::vector<int> &getEdgeClosure(int iEdge, int iSign) const
+ {
+ return edgeClosure[iSign == 1 ? iEdge : 3 + iEdge];
}
inline void evaluateMonomials(double u, double v, double w, double p[]) const
{
@@ -108,13 +110,13 @@ struct functionSpace
}
};
-class functionSpaces
+class polynomialBases
{
private:
- static std::map<int, functionSpace> fs;
+ static std::map<int, polynomialBasis> fs;
static std::map<std::pair<int, int>, fullMatrix<double> > injector;
public :
- static const functionSpace &find(int);
+ static const polynomialBasis &find(int);
static const fullMatrix<double> &findInjector(int, int);
};
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index c82b6eb396..4f1bf69a30 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -93,7 +93,7 @@ bool PViewDataGModel::finalize()
void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
{
int type = e->getType();
- const functionSpace *fs = e->getFunctionSpace();
+ const polynomialBasis *fs = e->getFunctionSpace();
if(fs){
if(e->getPolynomialOrder() > 1)
setInterpolationMatrices(type, fs->coefficients, fs->monomials,
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 1a9f11a5f6..31d3d44b31 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -6,7 +6,7 @@
#include "PViewDataList.h"
#include "GmshMessage.h"
#include "GmshDefines.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
#include "Numeric.h"
#include "SmoothData.h"
#include "Context.h"
@@ -784,7 +784,7 @@ void PViewDataList::setOrder2(int type)
case TYPE_PRI: typeMSH = MSH_PRI_18; break;
case TYPE_PYR: typeMSH = MSH_PYR_14; break;
}
- const functionSpace *fs = &functionSpaces::find(typeMSH);
+ const polynomialBasis *fs = &polynomialBases::find(typeMSH);
if(!fs){
Msg::Error("Could not find function space for element type %d", typeMSH);
return;
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index 6a2088756c..2e3e068f5b 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -15,13 +15,13 @@
#include "fullMatrix.h"
#include "Numeric.h"
-class convexCombinationTerm : public femTerm<double,double> {
+class convexCombinationTerm : public femTerm<double> {
protected:
const simpleFunction<double> *_k;
const int _iField;
public:
convexCombinationTerm(GModel *gm, int iField, simpleFunction<double> *k)
- : femTerm<double,double>(gm), _iField(iField), _k(k) {}
+ : femTerm<double>(gm), _iField(iField), _k(k) {}
virtual int sizeOfR(SElement *se) const
{
return se->getMeshElement()->getNumVertices();
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index 52521de117..8f9ce458f9 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -14,7 +14,7 @@
#include "fullMatrix.h"
#include "Numeric.h"
-class crossConfTerm : public femTerm<double, double> {
+class crossConfTerm : public femTerm<double> {
protected:
const simpleFunction<double> *_diffusivity;
const int _iFieldR;
@@ -22,7 +22,7 @@ class crossConfTerm : public femTerm<double, double> {
public:
crossConfTerm(GModel *gm, int iFieldR, int iFieldC,
simpleFunction<double> *diffusivity)
- : femTerm<double, double>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
+ : femTerm<double>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
_diffusivity(diffusivity) {}
virtual int sizeOfR(SElement *se) const
{
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 30c269ab42..de8bfe3ddb 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -1,6 +1,6 @@
#include "dgGroupOfElements.h"
#include "MElement.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
#include "Numeric.h"
#include "MTriangle.h"
#include "MLine.h"
@@ -20,9 +20,9 @@ static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
}
static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
- const functionSpace *fsFace,
+ const polynomialBasis *fsFace,
const fullMatrix<double> &intgFace,
- const functionSpace *fsElement,
+ const polynomialBasis *fsElement,
const std::vector <int> *closure) {
int npts=intgFace.size1();
fullMatrix<double> *m = new fullMatrix<double> (npts, 4);
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 6277359e34..f00b80a7b4 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -16,7 +16,7 @@
class MElement;
class MFace;
class MEdge;
-class functionSpace;
+class polynomialBasis;
class dgElement {
MElement *_element;
@@ -42,7 +42,7 @@ class dgGroupOfElements {
// the ONLY function space that is used to
// inerpolated the fields (may be different to the
// one of the elements)
- const functionSpace &_fs;
+ const polynomialBasis &_fs;
// Ni integration points, matrix of size Ni x 4 (u,v,w,weight)
fullMatrix<double> *_integration;
// collocation matrix that maps vertices to integration points.
@@ -108,9 +108,9 @@ class dgGroupOfFaces {
void addFace(const MFace &topoFace, int iElLeft, int iElRight);
void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
void computeFaceNormals();
- // Two functionSpaces for left and right elements
+ // Two polynomialBases for left and right elements
// the group has always the same types for left and right
- const functionSpace *_fsLeft,*_fsRight, *_fsFace;
+ const polynomialBasis *_fsLeft,*_fsRight, *_fsFace;
// N elements in the group
std::vector<int>_left, _right;
std::vector<MElement *>_faces;
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 9ff945e0ff..fd2d2fccb4 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -8,6 +8,7 @@
#include <vector>
#include <string>
+#include <complex>
#include <map>
#include "MVertex.h"
#include "linearSystem.h"
@@ -43,21 +44,42 @@ class Dof{
}
};
-template<class dataVec, class dataMat>
+template<class T> struct dofTraits
+{
+ typedef T VecType;
+ typedef T MatType;
+};
+
+template<> struct dofTraits<fullVector<double> >
+{
+ typedef fullVector<double> VecType;
+ typedef fullMatrix<double> MatType;
+};
+
+template<> struct dofTraits<fullVector<std::complex<double> > >
+{
+ typedef fullVector<std::complex<double> > VecType;
+ typedef fullMatrix<std::complex<double> > MatType;
+};
+
+template<class T>
class DofAffineConstraint{
public:
- std::vector<std::pair<Dof, dataMat> > linear;
- dataVec shift;
+ std::vector<std::pair<Dof, typename dofTraits<T>::MatType> > linear;
+ typename dofTraits<T>::VecType shift;
};
// data=float, double, complex<double>, smallVec<double>, smallVec<complex<double> >...
-template <class dataVec, class dataMat>
+template <class T>
class dofManager{
private:
+ typedef typename dofTraits<T>::VecType dataVec;
+ typedef typename dofTraits<T>::MatType dataMat;
+
// general affine constraint on sub-blocks, treated by adding
// equations:
// dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
- std::map<std::pair<Dof, dataMat>, DofAffineConstraint<dataVec, dataMat> > constraints;
+ std::map<std::pair<Dof, dataMat>, DofAffineConstraint<dataVec> > constraints;
// fixations on full blocks, treated by eliminating equations:
// DofVec = dataVec
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index ddfade0e24..8f771bdbbc 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -9,7 +9,7 @@
#include <slepceps.h>
-eigenSolver::eigenSolver(dofManager<double, double> *manager, std::string A,
+eigenSolver::eigenSolver(dofManager<double> *manager, std::string A,
std::string B) : _A(0), _B(0)
{
if(A.size()){
diff --git a/Solver/eigenSolver.h b/Solver/eigenSolver.h
index 0ca589d14e..ee3e8d9ff8 100644
--- a/Solver/eigenSolver.h
+++ b/Solver/eigenSolver.h
@@ -24,7 +24,7 @@ class eigenSolver{
std::vector<std::vector<std::complex<double> > > _eigenVectors;
void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
public:
- eigenSolver(dofManager<double, double> *manager, std::string A,
+ eigenSolver(dofManager<double> *manager, std::string A,
std::string B="");
void solve(int numEigenValues=0, std::string which="");
int getNumEigenValues(){ return _eigenValues.size(); }
@@ -38,7 +38,7 @@ class eigenSolver{
private:
std::vector<std::complex<double> > _dummy;
public:
- eigenSolver(dofManager<double, double> *manager, std::string A,
+ eigenSolver(dofManager<double> *manager, std::string A,
std::string B=""){}
void solve(int numEigenValues=0, std::string which="")
{
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 73f99fe40c..4358067e65 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -59,7 +59,7 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName
#endif
-static double vonMises(dofManager<double,double> *a, MElement *e,
+static double vonMises(dofManager<double> *a, MElement *e,
double u, double v, double w,
double E, double nu, int _tag)
{
@@ -195,7 +195,7 @@ void elasticitySolver::solve()
linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
lsys->setNoisy(2);
#endif
- pAssembler = new dofManager<double, double>(lsys);
+ pAssembler = new dofManager<double>(lsys);
std::map<int, std::vector<GEntity*> > groups[4];
pModel->getPhysicalGroups(groups);
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 10fa60f342..c0fedc0bb2 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -29,7 +29,7 @@ class elasticitySolver{
protected:
GModel *pModel;
int _dim, _tag;
- dofManager<double, double> *pAssembler;
+ dofManager<double> *pAssembler;
// young modulus and poisson coefficient per physical
std::vector<elasticField> elasticFields;
// imposed nodal forces
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index 1c0292cf88..09e3242901 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -12,7 +12,7 @@
#include "SElement.h"
#include "fullMatrix.h"
-class elasticityTerm : public femTerm<double, double> {
+class elasticityTerm : public femTerm<double> {
protected:
double _E, _nu;
int _iFieldR, _iFieldC;
@@ -49,9 +49,9 @@ class elasticityTerm : public femTerm<double, double> {
}
public:
elasticityTerm(GModel *gm, double E, double nu, int fieldr, int fieldc)
- : femTerm<double, double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldc) {}
+ : femTerm<double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldc) {}
elasticityTerm(GModel *gm, double E, double nu, int fieldr)
- : femTerm<double, double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldr) {}
+ : femTerm<double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldr) {}
void setVector(const SVector3 &f) { _volumeForce = f; }
void elementMatrix(SElement *se, fullMatrix<double> &m) const;
void elementVector(SElement *se, fullVector<double> &m) const;
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 5c35bbe4c0..3536058089 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -18,8 +18,11 @@
// a nodal finite element term : variables are always defined at nodes
// of the mesh
-template<class dataVec, class dataMat>
+template<class T>
class femTerm {
+ private:
+ typedef typename dofTraits<T>::VecType dataVec;
+ typedef typename dofTraits<T>::MatType dataMat;
protected:
GModel *_gm;
public:
@@ -46,7 +49,7 @@ class femTerm {
// add the contribution from all the elements in the intersection
// of two element groups L and C
- void addToMatrix(dofManager<dataVec, dataMat> &dm,
+ void addToMatrix(dofManager<dataVec> &dm,
groupOfElements &L,
groupOfElements &C) const
{
@@ -61,7 +64,7 @@ class femTerm {
}
// add the contribution from a single element to the dof manager
- void addToMatrix(dofManager<dataVec, dataMat> &dm, SElement *se) const
+ void addToMatrix(dofManager<dataVec> &dm, SElement *se) const
{
const int nbR = sizeOfR(se);
const int nbC = sizeOfC(se);
@@ -69,7 +72,7 @@ class femTerm {
elementMatrix(se, localMatrix);
addToMatrix(dm, localMatrix, se);
}
- void addToMatrix(dofManager<dataVec, dataMat> &dm,
+ void addToMatrix(dofManager<dataVec> &dm,
fullMatrix<dataMat> &localMatrix,
SElement *se) const
{
@@ -116,7 +119,7 @@ class femTerm {
}
void dirichletNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &e,
- dofManager<dataVec, dataMat> &dm)
+ dofManager<dataVec> &dm)
{
std::vector<MVertex *> v;
GModel *m = _gm;
@@ -126,7 +129,7 @@ class femTerm {
}
void neumannNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &fct,
- dofManager<dataVec, dataMat> &dm)
+ dofManager<dataVec> &dm)
{
std::map<int, std::vector<GEntity*> > groups[4];
GModel *m = _gm;
@@ -161,7 +164,7 @@ class femTerm {
}
}
- void addToRightHandSide(dofManager<dataVec, dataMat> &dm, groupOfElements &C) const
+ void addToRightHandSide(dofManager<dataVec> &dm, groupOfElements &C) const
{
groupOfElements::elementContainer::const_iterator it = C.begin();
for ( ; it != C.end(); ++it){
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 09af7cf50f..9d4401e203 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -17,7 +17,7 @@
// \nabla \cdot k \nabla U - a U
template<class scalar>
-class helmholtzTerm : public femTerm<scalar, scalar> {
+class helmholtzTerm : public femTerm<scalar> {
protected:
const simpleFunction<scalar> *_k, *_a;
const int _iFieldR;
@@ -25,7 +25,7 @@ class helmholtzTerm : public femTerm<scalar, scalar> {
public:
helmholtzTerm(GModel *gm, int iFieldR, int iFieldC, simpleFunction<scalar> *k,
simpleFunction<scalar> *a)
- : femTerm<scalar, scalar>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
+ : femTerm<scalar>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
_k(k), _a(a) {}
// one dof per vertex (nodal fem)
virtual int sizeOfR(SElement *se) const
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index e1d4b7b93c..0f5c09f579 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -120,7 +120,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
std::map<MVertex*,SPoint2> solution;
// PARAMETRIZE ---------------------------------
for (int step =0 ; step<2 ; step++){
- dofManager<double, double> myAssembler(_lsys);
+ dofManager<double> myAssembler(_lsys);
for(std::map<MVertex*,SPoint2>::iterator it = level.coordinates.begin();
it != level.coordinates.end(); ++it){
MVertex *v = it->first;
--
GitLab