style

- Make const and add `mutable` to appropriate member vars

examples/ex10.cpp

@@ -136,12 +136,12 @@ class ElasticEnergyCoefficient : public Coefficient
 private:
    HyperelasticModel  &model;
    const GridFunction &x;
-   DenseMatrix         J;
+   mutable DenseMatrix         J;
 
 public:
    ElasticEnergyCoefficient(HyperelasticModel &m, const GridFunction &x_)
       : model(m), x(x_) { }
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    virtual ~ElasticEnergyCoefficient() { }
 };
 
@@ -582,7 +582,7 @@ HyperelasticOperator::~HyperelasticOperator()
 
 
 double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    model.SetTransformation(T);
    x.GetVectorGradient(T, J);

examples/ex10p.cpp

@@ -141,12 +141,12 @@ class ElasticEnergyCoefficient : public Coefficient
 private:
    HyperelasticModel     &model;
    const ParGridFunction &x;
-   DenseMatrix            J;
+   mutable DenseMatrix            J;
 
 public:
    ElasticEnergyCoefficient(HyperelasticModel &m, const ParGridFunction &x_)
       : model(m), x(x_) { }
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    virtual ~ElasticEnergyCoefficient() { }
 };
 
@@ -661,7 +661,7 @@ HyperelasticOperator::~HyperelasticOperator()
 
 
 double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    model.SetTransformation(T);
    x.GetVectorGradient(T, J);

examples/ex17.cpp

@@ -60,7 +60,7 @@ protected:
    GridFunction *u; // displacement
    int si, sj; // component of the stress to evaluate, 0 <= si,sj < dim
 
-   DenseMatrix grad; // auxiliary matrix, used in Eval
+   mutable DenseMatrix grad; // auxiliary matrix, used in Eval
 
 public:
    StressCoefficient(Coefficient &lambda_, Coefficient &mu_)
@@ -69,7 +69,7 @@ public:
    void SetDisplacement(GridFunction &u_) { u = &u_; }
    void SetComponent(int i, int j) { si = i; sj = j; }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
 };
 
 // Simple GLVis visualization manager.
@@ -338,7 +338,7 @@ void InitDisplacement(const Vector &x, Vector &u)
 
 
 double StressCoefficient::Eval(ElementTransformation &T,
-                               const IntegrationPoint &ip)
+                               const IntegrationPoint &ip) const
 {
    MFEM_ASSERT(u != NULL, "displacement field is not set");
 

examples/ex17p.cpp

@@ -60,7 +60,7 @@ protected:
    GridFunction *u; // displacement
    int si, sj; // component of the stress to evaluate, 0 <= si,sj < dim
 
-   DenseMatrix grad; // auxiliary matrix, used in Eval
+   mutable DenseMatrix grad; // auxiliary matrix, used in Eval
 
 public:
    StressCoefficient(Coefficient &lambda_, Coefficient &mu_)
@@ -69,7 +69,7 @@ public:
    void SetDisplacement(GridFunction &u_) { u = &u_; }
    void SetComponent(int i, int j) { si = i; sj = j; }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
 };
 
 // Simple GLVis visualization manager.
@@ -377,7 +377,7 @@ void InitDisplacement(const Vector &x, Vector &u)
 
 
 double StressCoefficient::Eval(ElementTransformation &T,
-                               const IntegrationPoint &ip)
+                               const IntegrationPoint &ip) const
 {
    MFEM_ASSERT(u != NULL, "displacement field is not set");
 

examples/ex25.cpp

@@ -104,7 +104,7 @@ public:
    using VectorCoefficient::Eval;
 
    virtual void Eval(Vector &K, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       double x[3];
       Vector transip(x, 3);

examples/ex25p.cpp

@@ -103,7 +103,7 @@ public:
    using VectorCoefficient::Eval;
 
    virtual void Eval(Vector &K, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       double x[3];
       Vector transip(x, 3);

examples/ex36.cpp

@@ -53,7 +53,7 @@ public:
                                     double min_val_=-36)
       : u(&u_), obstacle(&obst_), min_val(min_val_) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
 };
 
 class ExponentialGridFunctionCoefficient : public Coefficient
@@ -69,7 +69,7 @@ public:
                                       double min_val_=0.0, double max_val_=1e6)
       : u(&u_), obstacle(&obst_), min_val(min_val_), max_val(max_val_) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
 };
 
 int main(int argc, char *argv[])
@@ -381,7 +381,7 @@ int main(int argc, char *argv[])
 }
 
 double LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
-                                              const IntegrationPoint &ip)
+                                              const IntegrationPoint &ip) const
 {
    MFEM_ASSERT(u != NULL, "grid function is not set");
 
@@ -390,7 +390,7 @@ double LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
 }
 
 double ExponentialGridFunctionCoefficient::Eval(ElementTransformation &T,
-                                                const IntegrationPoint &ip)
+                                                const IntegrationPoint &ip) const
 {
    MFEM_ASSERT(u != NULL, "grid function is not set");
 

examples/ex36p.cpp

@@ -53,7 +53,7 @@ public:
                                     double min_val_=-36)
       : u(&u_), obstacle(&obst_), min_val(min_val_) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
 };
 
 class ExponentialGridFunctionCoefficient : public Coefficient
@@ -69,7 +69,7 @@ public:
                                       double min_val_=0.0, double max_val_=1e6)
       : u(&u_), obstacle(&obst_), min_val(min_val_), max_val(max_val_) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
 };
 
 int main(int argc, char *argv[])
@@ -445,7 +445,7 @@ int main(int argc, char *argv[])
 }
 
 double LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
-                                              const IntegrationPoint &ip)
+                                              const IntegrationPoint &ip) const
 {
    MFEM_ASSERT(u != NULL, "grid function is not set");
 
@@ -454,7 +454,7 @@ double LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
 }
 
 double ExponentialGridFunctionCoefficient::Eval(ElementTransformation &T,
-                                                const IntegrationPoint &ip)
+                                                const IntegrationPoint &ip) const
 {
    MFEM_ASSERT(u != NULL, "grid function is not set");
 

examples/ex37.hpp

@@ -53,7 +53,7 @@ public:
 
 
    virtual double Eval(ElementTransformation &T,
-                       const IntegrationPoint &ip)
+                       const IntegrationPoint &ip) const
    {
       return fun(GridFunctionCoefficient::Eval(T, ip));
    }
@@ -84,7 +84,7 @@ public:
        fun(fun_) {}
 
    virtual double Eval(ElementTransformation &T,
-                       const IntegrationPoint &ip)
+                       const IntegrationPoint &ip) const
    {
       const double value1 = fun(GridFunctionCoefficient::Eval(T, ip));
       const double value2 = fun(OtherGridF_cf.Eval(T, ip));
@@ -108,7 +108,7 @@ public:
       : rho_filter(rho_filter_), min_val(min_val_), max_val(max_val_),
         exponent(exponent_) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       double val = rho_filter->GetValue(T, ip);
       double coeff = min_val + pow(val,exponent)*(max_val-min_val);
@@ -125,7 +125,7 @@ protected:
    Coefficient * mu=nullptr;
    GridFunction *u = nullptr; // displacement
    GridFunction *rho_filter = nullptr; // filter density
-   DenseMatrix grad; // auxiliary matrix, used in Eval
+   mutable DenseMatrix grad; // auxiliary matrix, used in Eval
    double exponent;
    double rho_min;
 
@@ -142,7 +142,7 @@ public:
       MFEM_ASSERT(rho_filter, "density field is not set");
    }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       double L = lambda->Eval(T, ip);
       double M = mu->Eval(T, ip);
@@ -177,7 +177,7 @@ public:
    using VectorCoefficient::Eval;
 
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       Vector xx; xx.SetSize(T.GetDimension());
       T.Transform(ip,xx);

examples/petsc/ex10p.cpp

@@ -170,7 +170,7 @@ private:
 public:
    ElasticEnergyCoefficient(HyperelasticModel &m, const ParGridFunction &x_)
       : model(m), x(x_) { }
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    virtual ~ElasticEnergyCoefficient() { }
 };
 
@@ -743,7 +743,7 @@ Solver* PreconditionerFactory::NewPreconditioner(const mfem::OperatorHandle& oh)
 }
 
 double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    model.SetTransformation(T);
    x.GetVectorGradient(T, J);

examples/sundials/ex10.cpp

@@ -195,7 +195,7 @@ private:
 public:
    ElasticEnergyCoefficient(HyperelasticModel &m, const GridFunction &x_)
       : model(m), x(x_) { }
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    virtual ~ElasticEnergyCoefficient() { }
 };
 
@@ -798,7 +798,7 @@ HyperelasticOperator::~HyperelasticOperator()
 
 
 double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    model.SetTransformation(T);
    x.GetVectorGradient(T, J);

examples/sundials/ex10p.cpp

@@ -200,7 +200,7 @@ private:
 public:
    ElasticEnergyCoefficient(HyperelasticModel &m, const ParGridFunction &x_)
       : model(m), x(x_) { }
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    virtual ~ElasticEnergyCoefficient() { }
 };
 
@@ -883,7 +883,7 @@ HyperelasticOperator::~HyperelasticOperator()
 
 
 double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    model.SetTransformation(T);
    x.GetVectorGradient(T, J);

fem/bilininteg.cpp

@@ -4310,7 +4310,7 @@ struct ShapeCoefficient : public VectorCoefficient
 
    using VectorCoefficient::Eval;
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       V.SetSize(vdim);
       fe.CalcPhysShape(T, V);
@@ -4351,9 +4351,9 @@ ScalarVectorProductInterpolator::AssembleElementMatrix2(
 
       VShapeCoefficient(Coefficient &q, const FiniteElement &fe_, int sdim)
          : MatrixCoefficient(fe_.GetDof(), sdim), Q(q), fe(fe_) { }
-
+      using MatrixCoefficient::Eval;
       virtual void Eval(DenseMatrix &M, ElementTransformation &T,
-                        const IntegrationPoint &ip)
+                        const IntegrationPoint &ip) const
       {
          M.SetSize(height, width);
          fe.CalcPhysVShape(T, M);
@@ -4383,14 +4383,16 @@ VectorScalarProductInterpolator::AssembleElementMatrix2(
    {
       VectorCoefficient &VQ;
       const FiniteElement &fe;
-      Vector vc, shape;
+
+      mutable Vector vc, shape;
 
       VecShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
          : MatrixCoefficient(fe_.GetDof(), vq.GetVDim()), VQ(vq), fe(fe_),
            vc(width), shape(height) { }
 
+      using MatrixCoefficient::Eval;
       virtual void Eval(DenseMatrix &M, ElementTransformation &T,
-                        const IntegrationPoint &ip)
+                        const IntegrationPoint &ip) const
       {
          M.SetSize(height, width);
          VQ.Eval(vc, T, ip);
@@ -4421,8 +4423,9 @@ ScalarCrossProductInterpolator::AssembleElementMatrix2(
    {
       VectorCoefficient &VQ;
       const FiniteElement &fe;
-      DenseMatrix vshape;
-      Vector vc;
+
+      mutable DenseMatrix vshape;
+      mutable Vector vc;
 
       VCrossVShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
          : VectorCoefficient(fe_.GetDof()), VQ(vq), fe(fe_),
@@ -4430,7 +4433,7 @@ ScalarCrossProductInterpolator::AssembleElementMatrix2(
 
       using VectorCoefficient::Eval;
       virtual void Eval(Vector &V, ElementTransformation &T,
-                        const IntegrationPoint &ip)
+                        const IntegrationPoint &ip) const
       {
          V.SetSize(vdim);
          VQ.Eval(vc, T, ip);
@@ -4463,8 +4466,9 @@ VectorCrossProductInterpolator::AssembleElementMatrix2(
    {
       VectorCoefficient &VQ;
       const FiniteElement &fe;
-      DenseMatrix vshape;
-      Vector vc;
+
+      mutable DenseMatrix vshape;
+      mutable Vector vc;
 
       VCrossVShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
          : MatrixCoefficient(fe_.GetDof(), vq.GetVDim()), VQ(vq), fe(fe_),
@@ -4473,8 +4477,9 @@ VectorCrossProductInterpolator::AssembleElementMatrix2(
          MFEM_ASSERT(width == 3, "");
       }
 
+      using MatrixCoefficient::Eval;
       virtual void Eval(DenseMatrix &M, ElementTransformation &T,
-                        const IntegrationPoint &ip)
+                        const IntegrationPoint &ip) const
       {
          M.SetSize(height, width);
          VQ.Eval(vc, T, ip);
@@ -4514,8 +4519,9 @@ struct VDotVShapeCoefficient : public VectorCoefficient
 {
    VectorCoefficient &VQ;
    const FiniteElement &fe;
-   DenseMatrix vshape;
-   Vector vc;
+
+   mutable DenseMatrix vshape;
+   mutable Vector vc;
 
    VDotVShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
       : VectorCoefficient(fe_.GetDof()), VQ(vq), fe(fe_),
@@ -4523,7 +4529,7 @@ struct VDotVShapeCoefficient : public VectorCoefficient
 
    using VectorCoefficient::Eval;
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       V.SetSize(vdim);
       VQ.Eval(vc, T, ip);

fem/coefficient.cpp

@@ -48,7 +48,7 @@ ElementTransformation *RefinedToCoarse(
    return coarse_T;
 }
 
-void Coefficient::Project(QuadratureFunction &qf)
+void Coefficient::Project(QuadratureFunction &qf) const
 {
    QuadratureSpaceBase &qspace = *qf.GetSpace();
    const int ne = qspace.GetNE();
@@ -68,13 +68,13 @@ void Coefficient::Project(QuadratureFunction &qf)
    }
 }
 
-void ConstantCoefficient::Project(QuadratureFunction &qf)
+void ConstantCoefficient::Project(QuadratureFunction &qf) const
 {
    qf = constant;
 }
 
 double PWConstCoefficient::Eval(ElementTransformation & T,
-                                const IntegrationPoint & ip)
+                                const IntegrationPoint & ip) const
 {
    int att = T.Attribute;
    return (constants(att-1));
@@ -112,7 +112,7 @@ void PWCoefficient::SetTime(double t)
 }
 
 double PWCoefficient::Eval(ElementTransformation &T,
-                           const IntegrationPoint &ip)
+                           const IntegrationPoint &ip) const
 {
    const int att = T.Attribute;
    std::map<int, Coefficient*>::const_iterator p = pieces.find(att);
@@ -127,7 +127,7 @@ double PWCoefficient::Eval(ElementTransformation &T,
 }
 
 double FunctionCoefficient::Eval(ElementTransformation & T,
-                                 const IntegrationPoint & ip)
+                                 const IntegrationPoint & ip) const
 {
    double x[3];
    Vector transip(x, 3);
@@ -145,42 +145,42 @@ double FunctionCoefficient::Eval(ElementTransformation & T,
 }
 
 double CartesianCoefficient::Eval(ElementTransformation & T,
-                                  const IntegrationPoint & ip)
+                                  const IntegrationPoint & ip) const
 {
    T.Transform(ip, transip);
    return transip[comp];
 }
 
 double CylindricalRadialCoefficient::Eval(ElementTransformation & T,
-                                          const IntegrationPoint & ip)
+                                          const IntegrationPoint & ip) const
 {
    T.Transform(ip, transip);
    return sqrt(transip[0] * transip[0] + transip[1] * transip[1]);
 }
 
 double CylindricalAzimuthalCoefficient::Eval(ElementTransformation & T,
-                                             const IntegrationPoint & ip)
+                                             const IntegrationPoint & ip) const
 {
    T.Transform(ip, transip);
    return atan2(transip[1], transip[0]);
 }
 
 double SphericalRadialCoefficient::Eval(ElementTransformation & T,
-                                        const IntegrationPoint & ip)
+                                        const IntegrationPoint & ip) const
 {
    T.Transform(ip, transip);
    return sqrt(transip * transip);
 }
 
 double SphericalAzimuthalCoefficient::Eval(ElementTransformation & T,
-                                           const IntegrationPoint & ip)
+                                           const IntegrationPoint & ip) const
 {
    T.Transform(ip, transip);
    return atan2(transip[1], transip[0]);
 }
 
 double SphericalPolarCoefficient::Eval(ElementTransformation & T,
-                                       const IntegrationPoint & ip)
+                                       const IntegrationPoint & ip) const
 {
    T.Transform(ip, transip);
    return atan2(sqrt(transip[0] * transip[0] + transip[1] * transip[1]),
@@ -188,7 +188,7 @@ double SphericalPolarCoefficient::Eval(ElementTransformation & T,
 }
 
 double GridFunctionCoefficient::Eval (ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    Mesh *gf_mesh = GridF->FESpace()->GetMesh();
    if (T.mesh->GetNE() == gf_mesh->GetNE())
@@ -203,7 +203,7 @@ double GridFunctionCoefficient::Eval (ElementTransformation &T,
    }
 }
 
-void GridFunctionCoefficient::Project(QuadratureFunction &qf)
+void GridFunctionCoefficient::Project(QuadratureFunction &qf) const
 {
    qf.ProjectGridFunction(*GridF);
 }
@@ -216,16 +216,16 @@ void TransformedCoefficient::SetTime(double t)
 }
 
 double TransformedCoefficient::Eval(ElementTransformation &T,
-                                    const IntegrationPoint &ip)
+                                    const IntegrationPoint &ip) const
 {
    if (Q2)
    {
-      return Transform2(Q1->Eval(T, ip, GetTime()),
-                        Q2->Eval(T, ip, GetTime()));
+      return Transform2(Q1->Eval(T, ip),
+                        Q2->Eval(T, ip));
    }
    else
    {
-      return Transform1(Q1->Eval(T, ip, GetTime()));
+      return Transform1(Q1->Eval(T, ip));
    }
 }
 
@@ -243,17 +243,17 @@ void DeltaCoefficient::SetDeltaCenter(const Vector& vcenter)
    sdim = vcenter.Size();
 }
 
-void DeltaCoefficient::GetDeltaCenter(Vector& vcenter)
+void DeltaCoefficient::GetDeltaCenter(Vector& vcenter) const
 {
    vcenter.SetSize(sdim);
    vcenter = center;
 }
 
 double DeltaCoefficient::EvalDelta(ElementTransformation &T,
-                                   const IntegrationPoint &ip)
+                                   const IntegrationPoint &ip) const
 {
    double w = Scale();
-   return weight ? weight->Eval(T, ip, GetTime())*w : w;
+   return weight ? weight->Eval(T, ip)*w : w;
 }
 
 void RestrictedCoefficient::SetTime(double t)
@@ -263,7 +263,7 @@ void RestrictedCoefficient::SetTime(double t)
 }
 
 void VectorCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                             const IntegrationRule &ir)
+                             const IntegrationRule &ir) const
 {
    Vector Mi;
    M.SetSize(vdim, ir.GetNPoints());
@@ -276,7 +276,7 @@ void VectorCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
    }
 }
 
-void VectorCoefficient::Project(QuadratureFunction &qf)
+void VectorCoefficient::Project(QuadratureFunction &qf) const
 {
    MFEM_VERIFY(vdim == qf.GetVDim(), "Wrong sizes.");
    QuadratureSpaceBase &qspace = *qf.GetSpace();
@@ -339,7 +339,7 @@ void PWVectorCoefficient::SetTime(double t)
 }
 
 void PWVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
-                               const IntegrationPoint &ip)
+                               const IntegrationPoint &ip) const
 {
    const int att = T.Attribute;
    std::map<int, VectorCoefficient*>::const_iterator p = pieces.find(att);
@@ -357,14 +357,20 @@ void PWVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
 }
 
 void PositionVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                     const IntegrationPoint &ip)
+                                     const IntegrationPoint &ip) const
 {
    V.SetSize(vdim);
    T.Transform(ip, V);
 }
 
+void VectorFunctionCoefficient::SetTime(double t)
+{
+   if (Q) { Q->SetTime(t); }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void VectorFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                     const IntegrationPoint &ip)
+                                     const IntegrationPoint &ip) const
 {
    double x[3];
    Vector transip(x, 3);
@@ -382,7 +388,7 @@ void VectorFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
    }
    if (Q)
    {
-      V *= Q->Eval(T, ip, GetTime());
+      V *= Q->Eval(T, ip);
    }
 }
 
@@ -421,7 +427,7 @@ VectorArrayCoefficient::~VectorArrayCoefficient()
 }
 
 void VectorArrayCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                  const IntegrationPoint &ip)
+                                  const IntegrationPoint &ip) const
 {
    V.SetSize(vdim);
    for (int i = 0; i < vdim; i++)
@@ -443,7 +449,7 @@ void VectorGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
 }
 
 void VectorGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                         const IntegrationPoint &ip)
+                                         const IntegrationPoint &ip) const
 {
    Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
    if (T.mesh->GetNE() == gf_mesh->GetNE())
@@ -459,7 +465,7 @@ void VectorGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
 }
 
 void VectorGridFunctionCoefficient::Eval(
-   DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
+   DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir) const
 {
    if (T.mesh == GridFunc->FESpace()->GetMesh())
    {
@@ -471,7 +477,7 @@ void VectorGridFunctionCoefficient::Eval(
    }
 }
 
-void VectorGridFunctionCoefficient::Project(QuadratureFunction &qf)
+void VectorGridFunctionCoefficient::Project(QuadratureFunction &qf) const
 {
    qf.ProjectGridFunction(*GridFunc);
 }
@@ -491,7 +497,7 @@ void GradientGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
 }
 
 void GradientGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                           const IntegrationPoint &ip)
+                                           const IntegrationPoint &ip) const
 {
    Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
    if (T.mesh->GetNE() == gf_mesh->GetNE())
@@ -507,7 +513,7 @@ void GradientGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
 }
 
 void GradientGridFunctionCoefficient::Eval(
-   DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
+   DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir) const
 {
    if (T.mesh == GridFunc->FESpace()->GetMesh())
    {
@@ -532,7 +538,7 @@ void CurlGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
 }
 
 void CurlGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                       const IntegrationPoint &ip)
+                                       const IntegrationPoint &ip) const
 {
    Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
    if (T.mesh->GetNE() == gf_mesh->GetNE())
@@ -554,7 +560,7 @@ DivergenceGridFunctionCoefficient::DivergenceGridFunctionCoefficient (
 }
 
 double DivergenceGridFunctionCoefficient::Eval(ElementTransformation &T,
-                                               const IntegrationPoint &ip)
+                                               const IntegrationPoint &ip) const
 {
    Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
    if (T.mesh->GetNE() == gf_mesh->GetNE())
@@ -582,10 +588,9 @@ void VectorDeltaCoefficient::SetDirection(const Vector &d_)
 }
 
 void VectorDeltaCoefficient::EvalDelta(
-   Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   Vector &V, ElementTransformation &T, const IntegrationPoint &ip) const
 {
    V = dir;
-   d.SetTime(GetTime());
    V *= d.EvalDelta(T, ip);
 }
 
@@ -596,12 +601,11 @@ void VectorRestrictedCoefficient::SetTime(double t)
 }
 
 void VectorRestrictedCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                       const IntegrationPoint &ip)
+                                       const IntegrationPoint &ip) const
 {
    V.SetSize(vdim);
    if (active_attr[T.Attribute-1])
    {
-      c->SetTime(GetTime());
       c->Eval(V, T, ip);
    }
    else
@@ -611,11 +615,10 @@ void VectorRestrictedCoefficient::Eval(Vector &V, ElementTransformation &T,
 }
 
 void VectorRestrictedCoefficient::Eval(
-   DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
+   DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir) const
 {
    if (active_attr[T.Attribute-1])
    {
-      c->SetTime(GetTime());
       c->Eval(M, T, ir);
    }
    else
@@ -625,7 +628,7 @@ void VectorRestrictedCoefficient::Eval(
    }
 }
 
-void MatrixCoefficient::Project(QuadratureFunction &qf, bool transpose)
+void MatrixCoefficient::Project(QuadratureFunction &qf, bool transpose) const
 {
    MFEM_VERIFY(qf.GetVDim() == height*width, "Wrong sizes.");
    QuadratureSpaceBase &qspace = *qf.GetSpace();
@@ -697,7 +700,7 @@ void PWMatrixCoefficient::SetTime(double t)
 }
 
 void PWMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
-                               const IntegrationPoint &ip)
+                               const IntegrationPoint &ip) const
 {
    const int att = T.Attribute;
    std::map<int, MatrixCoefficient*>::const_iterator p = pieces.find(att);
@@ -721,7 +724,7 @@ void MatrixFunctionCoefficient::SetTime(double t)
 }
 
 void MatrixFunctionCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
-                                     const IntegrationPoint &ip)
+                                     const IntegrationPoint &ip) const
 {
    double x[3];
    Vector transip(x, 3);
@@ -771,7 +774,7 @@ void MatrixFunctionCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
 
    if (Q)
    {
-      K *= Q->Eval(T, ip, GetTime());
+      K *= Q->Eval(T, ip);
    }
 }
 
@@ -826,7 +829,7 @@ void SymmetricMatrixCoefficient::ProjectSymmetric(QuadratureFunction &qf)
 
 
 void SymmetricMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    mat.SetSize(height);
    Eval(mat, T, ip);
@@ -847,7 +850,7 @@ void SymmetricMatrixFunctionCoefficient::SetTime(double t)
 
 void SymmetricMatrixFunctionCoefficient::Eval(DenseSymmetricMatrix &K,
                                               ElementTransformation &T,
-                                              const IntegrationPoint &ip)
+                                              const IntegrationPoint &ip) const
 {
    double x[3];
    Vector transip(x, 3);
@@ -871,7 +874,7 @@ void SymmetricMatrixFunctionCoefficient::Eval(DenseSymmetricMatrix &K,
 
    if (Q)
    {
-      K *= Q->Eval(T, ip, GetTime());
+      K *= Q->Eval(T, ip);
    }
 }
 
@@ -912,7 +915,7 @@ MatrixArrayCoefficient::~MatrixArrayCoefficient ()
 }
 
 void MatrixArrayCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
-                                  const IntegrationPoint &ip)
+                                  const IntegrationPoint &ip) const
 {
    K.SetSize(height, width);
    for (int i = 0; i < height; i++)
@@ -931,11 +934,10 @@ void MatrixRestrictedCoefficient::SetTime(double t)
 }
 
 void MatrixRestrictedCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
-                                       const IntegrationPoint &ip)
+                                       const IntegrationPoint &ip) const
 {
    if (active_attr[T.Attribute-1])
    {
-      c->SetTime(GetTime());
       c->Eval(K, T, ip);
    }
    else
@@ -989,7 +991,7 @@ void InnerProductCoefficient::SetTime(double t)
 }
 
 double InnerProductCoefficient::Eval(ElementTransformation &T,
-                                     const IntegrationPoint &ip)
+                                     const IntegrationPoint &ip) const
 {
    a->Eval(va, T, ip);
    b->Eval(vb, T, ip);
@@ -1013,7 +1015,7 @@ void VectorRotProductCoefficient::SetTime(double t)
 }
 
 double VectorRotProductCoefficient::Eval(ElementTransformation &T,
-                                         const IntegrationPoint &ip)
+                                         const IntegrationPoint &ip) const
 {
    a->Eval(va, T, ip);
    b->Eval(vb, T, ip);
@@ -1035,7 +1037,7 @@ void DeterminantCoefficient::SetTime(double t)
 }
 
 double DeterminantCoefficient::Eval(ElementTransformation &T,
-                                    const IntegrationPoint &ip)
+                                    const IntegrationPoint &ip) const
 {
    a->Eval(ma, T, ip);
    return ma.Det();
@@ -1092,7 +1094,7 @@ void VectorSumCoefficient::SetTime(double t)
 }
 
 void VectorSumCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                const IntegrationPoint &ip)
+                                const IntegrationPoint &ip) const
 {
    V.SetSize(A.Size());
    if (    ACoef) { ACoef->Eval(A, T, ip); }
@@ -1122,7 +1124,7 @@ void ScalarVectorProductCoefficient::SetTime(double t)
 }
 
 void ScalarVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                          const IntegrationPoint &ip)
+                                          const IntegrationPoint &ip) const
 {
    double sa = (a == NULL) ? aConst : a->Eval(T, ip);
    b->Eval(V, T, ip);
@@ -1141,7 +1143,7 @@ void NormalizedVectorCoefficient::SetTime(double t)
 }
 
 void NormalizedVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                       const IntegrationPoint &ip)
+                                       const IntegrationPoint &ip) const
 {
    a->Eval(V, T, ip);
    double nv = V.Norml2();
@@ -1166,7 +1168,7 @@ void VectorCrossProductCoefficient::SetTime(double t)
 }
 
 void VectorCrossProductCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                         const IntegrationPoint &ip)
+                                         const IntegrationPoint &ip) const
 {
    a->Eval(va, T, ip);
    b->Eval(vb, T, ip);
@@ -1194,7 +1196,7 @@ void MatrixVectorProductCoefficient::SetTime(double t)
 }
 
 void MatrixVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                          const IntegrationPoint &ip)
+                                          const IntegrationPoint &ip) const
 {
    a->Eval(ma, T, ip);
    b->Eval(vb, T, ip);
@@ -1203,7 +1205,7 @@ void MatrixVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
 }
 
 void IdentityMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                                     const IntegrationPoint &ip)
+                                     const IntegrationPoint &ip) const
 {
    M.SetSize(dim);
    M = 0.0;
@@ -1230,7 +1232,7 @@ void MatrixSumCoefficient::SetTime(double t)
 }
 
 void MatrixSumCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                                const IntegrationPoint &ip)
+                                const IntegrationPoint &ip) const
 {
    b->Eval(M, T, ip);
    if ( beta != 1.0 ) { M *= beta; }
@@ -1251,7 +1253,7 @@ MatrixProductCoefficient::MatrixProductCoefficient(MatrixCoefficient &A,
 }
 
 void MatrixProductCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                                    const IntegrationPoint &ip)
+                                    const IntegrationPoint &ip) const
 {
    a->Eval(ma, T, ip);
    b->Eval(mb, T, ip);
@@ -1279,7 +1281,7 @@ void ScalarMatrixProductCoefficient::SetTime(double t)
 
 void ScalarMatrixProductCoefficient::Eval(DenseMatrix &M,
                                           ElementTransformation &T,
-                                          const IntegrationPoint &ip)
+                                          const IntegrationPoint &ip) const
 {
    double sa = (a == NULL) ? aConst : a->Eval(T, ip);
    b->Eval(M, T, ip);
@@ -1298,7 +1300,7 @@ void TransposeMatrixCoefficient::SetTime(double t)
 
 void TransposeMatrixCoefficient::Eval(DenseMatrix &M,
                                       ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    a->Eval(M, T, ip);
    M.Transpose();
@@ -1320,7 +1322,7 @@ void InverseMatrixCoefficient::SetTime(double t)
 
 void InverseMatrixCoefficient::Eval(DenseMatrix &M,
                                     ElementTransformation &T,
-                                    const IntegrationPoint &ip)
+                                    const IntegrationPoint &ip) const
 {
    a->Eval(M, T, ip);
    M.Invert();
@@ -1340,7 +1342,7 @@ void OuterProductCoefficient::SetTime(double t)
 }
 
 void OuterProductCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                                   const IntegrationPoint &ip)
+                                   const IntegrationPoint &ip) const
 {
    a->Eval(va, T, ip);
    b->Eval(vb, T, ip);
@@ -1373,7 +1375,7 @@ void CrossCrossCoefficient::SetTime(double t)
 }
 
 void CrossCrossCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                                 const IntegrationPoint &ip)
+                                 const IntegrationPoint &ip) const
 {
    k->Eval(vk, T, ip);
    M.SetSize(vk.Size(), vk.Size());
@@ -1587,7 +1589,7 @@ void VectorQuadratureFunctionCoefficient::SetComponent(int index_, int length_)
 
 void VectorQuadratureFunctionCoefficient::Eval(Vector &V,
                                                ElementTransformation &T,
-                                               const IntegrationPoint &ip)
+                                               const IntegrationPoint &ip) const
 {
    QuadF.HostRead();
 
@@ -1616,7 +1618,7 @@ void VectorQuadratureFunctionCoefficient::Eval(Vector &V,
    return;
 }
 
-void VectorQuadratureFunctionCoefficient::Project(QuadratureFunction &qf)
+void VectorQuadratureFunctionCoefficient::Project(QuadratureFunction &qf) const
 {
    qf = QuadF;
 }
@@ -1628,7 +1630,7 @@ QuadratureFunctionCoefficient::QuadratureFunctionCoefficient(
 }
 
 double QuadratureFunctionCoefficient::Eval(ElementTransformation &T,
-                                           const IntegrationPoint &ip)
+                                           const IntegrationPoint &ip) const
 {
    QuadF.HostRead();
    Vector temp(1);
@@ -1641,7 +1643,7 @@ double QuadratureFunctionCoefficient::Eval(ElementTransformation &T,
    return temp[0];
 }
 
-void QuadratureFunctionCoefficient::Project(QuadratureFunction &qf)
+void QuadratureFunctionCoefficient::Project(QuadratureFunction &qf) const
 {
    qf = QuadF;
 }

fem/gridfunc.cpp

@@ -4449,7 +4449,7 @@ double ComputeElementLpDistance(double p, int i,
 
 
 double ExtrudeCoefficient::Eval(ElementTransformation &T,
-                                const IntegrationPoint &ip)
+                                const IntegrationPoint &ip) const
 {
    ElementTransformation *T_in =
       mesh_in->GetElementTransformation(T.ElementNo / n);

fem/gridfunc.hpp

@@ -874,7 +874,8 @@ private:
 public:
    ExtrudeCoefficient(Mesh *m, Coefficient &s, int n_)
       : n(n_), mesh_in(m), sol_in(s) { }
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   using Coefficient::Eval;
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    virtual ~ExtrudeCoefficient() { }
 };
 

mesh/mesh.cpp

@@ -13398,7 +13398,7 @@ NodeExtrudeCoefficient::NodeExtrudeCoefficient(const int dim, const int n_,
 }
 
 void NodeExtrudeCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                  const IntegrationPoint &ip)
+                                  const IntegrationPoint &ip) const
 {
    V.SetSize(vdim);
    T.Transform(ip, tip);

mesh/mesh.hpp

@@ -2550,13 +2550,14 @@ class NodeExtrudeCoefficient : public VectorCoefficient
 private:
    int n, layer;
    double p[2], s;
-   Vector tip;
+
+   mutable Vector tip;
 public:
    NodeExtrudeCoefficient(const int dim, const int n_, const double s_);
    void SetLayer(const int l) { layer = l; }
    using VectorCoefficient::Eval;
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip);
+                     const IntegrationPoint &ip) const;
    virtual ~NodeExtrudeCoefficient() { }
 };
 

miniapps/common/dist_solver.cpp

@@ -285,7 +285,7 @@ void HeatDistanceSolver::ComputeScalarDistance(Coefficient &zero_level_set,
 }
 
 double NormalizationDistanceSolver::NormalizationCoeff::
-Eval(ElementTransformation &T,const IntegrationPoint &ip)
+Eval(ElementTransformation &T,const IntegrationPoint &ip) const
 {
    T.SetIntPoint(&ip);
    Vector u_grad;

miniapps/common/dist_solver.hpp

@@ -90,7 +90,7 @@ private:
 
    public:
       NormalizationCoeff(ParGridFunction &u_gf) : u(u_gf) { }
-      virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+      virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    };
 
 public:
@@ -132,7 +132,7 @@ public:
 
    using VectorCoefficient::Eval;
 
-   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip) const
    {
       T.SetIntPoint(&ip);
 
@@ -153,7 +153,7 @@ public:
    PProductCoefficient(Coefficient& basec_, Coefficient& corrc_)
       : basef(basec_), corrf(corrc_) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       T.SetIntPoint(&ip);
       double u = basef.Eval(T,ip);

miniapps/common/mesh_extras.cpp

@@ -235,7 +235,7 @@ void AttrToMarker(int max_attr, const Array<int> &attrs, Array<int> &marker)
 }
 
 void KershawTransformation::Eval(Vector &V, ElementTransformation &T,
-                                 const IntegrationPoint &ip)
+                                 const IntegrationPoint &ip) const
 {
    V = 0.0;
    Vector pos(dim);

miniapps/common/mesh_extras.hpp

@@ -76,20 +76,20 @@ public:
    }
 
    // 1D transformation at the right boundary.
-   double right(const double eps, const double x)
+   double right(const double eps, const double x) const
    {
       return (x <= 0.5) ? (2-eps) * x : 1 + eps*(x-1);
    }
 
    // 1D transformation at the left boundary
-   double left(const double eps, const double x)
+   double left(const double eps, const double x) const
    {
       return 1-right(eps,1-x);
    }
 
    // Transition from a value of "a" for x=0, to a value of "b" for x=1.
    // Controlled through "smooth" parameter.
-   double step(const double a, const double b, double x)
+   double step(const double a, const double b, double x) const
    {
       if (x <= 0) { return a; }
       if (x >= 1) { return b; }
@@ -99,7 +99,7 @@ public:
    }
 
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip);
+                     const IntegrationPoint &ip) const;
 
    using VectorCoefficient::Eval;
 };

miniapps/dpg/util/pml.hpp

@@ -81,7 +81,7 @@ public:
    PmlCoefficient(double (*F)(const Vector &, CartesianPML *), CartesianPML * pml_)
       : pml(pml_), Function(F)
    {}
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       double x[3];
       Vector transip(x, 3);
@@ -102,7 +102,7 @@ public:
       : MatrixCoefficient(dim), pml(pml_), Function(F)
    {}
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       double x[3];
       Vector transip(x, 3);

miniapps/electromagnetics/joule_solver.cpp

@@ -902,7 +902,7 @@ void MagneticDiffusionEOperator::Debug(const char *base, double)
 }
 
 double JouleHeatingCoefficient::Eval(ElementTransformation &T,
-                                     const IntegrationPoint &ip)
+                                     const IntegrationPoint &ip) const
 {
    Vector E;
    double thisSigma;
@@ -930,7 +930,7 @@ MeshDependentCoefficient::MeshDependentCoefficient(
 }
 
 double MeshDependentCoefficient::Eval(ElementTransformation &T,
-                                      const IntegrationPoint &ip)
+                                      const IntegrationPoint &ip) const
 {
    // given the attribute, extract the coefficient value from the map
    std::map<int, double>::iterator it;
@@ -957,7 +957,7 @@ ScaledGFCoefficient::ScaledGFCoefficient(GridFunction *gf,
    : GridFunctionCoefficient(gf), mdc(input_mdc) {}
 
 double ScaledGFCoefficient::Eval(ElementTransformation &T,
-                                 const IntegrationPoint &ip)
+                                 const IntegrationPoint &ip) const
 {
    return mdc.Eval(T,ip) * GridFunctionCoefficient::Eval(T,ip);
 }

miniapps/electromagnetics/joule_solver.hpp

@@ -53,7 +53,7 @@ public:
    MeshDependentCoefficient(const std::map<int, double> &inputMap,
                             double scale = 1.0);
    MeshDependentCoefficient(const MeshDependentCoefficient &cloneMe);
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    void SetScaleFactor(const double &scale) { scaleFactor = scale; }
    virtual ~MeshDependentCoefficient()
    {
@@ -70,7 +70,7 @@ private:
    MeshDependentCoefficient mdc;
 public:
    ScaledGFCoefficient(GridFunction *gf, MeshDependentCoefficient &input_mdc);
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    void SetMDC(const MeshDependentCoefficient &input_mdc) { mdc = input_mdc; }
    virtual ~ScaledGFCoefficient() {}
 };
@@ -235,7 +235,7 @@ public:
    JouleHeatingCoefficient(const MeshDependentCoefficient &sigma_,
                            ParGridFunction &E_gf_)
       : E_gf(E_gf_), sigma(sigma_) {}
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const;
    virtual ~JouleHeatingCoefficient() {}
 };
 

miniapps/meshing/mesh-optimizer.hpp

@@ -139,7 +139,7 @@ public:
       : TMOPMatrixCoefficient(dim), metric(metric_id) { }
 
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       Vector pos(3);
       T.Transform(ip, pos);
@@ -250,7 +250,7 @@ public:
         hr_target_type(hr_target_type_) { }
 
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       Vector pos(3);
       T.Transform(ip, pos);

miniapps/meshing/reflector.cpp

@@ -68,13 +68,13 @@ public:
    { }
 
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip);
+                     const IntegrationPoint &ip) const;
 
    using VectorCoefficient::Eval;
 };
 
 void ReflectedCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                const IntegrationPoint &ip)
+                                const IntegrationPoint &ip) const
 {
    double x[3];
    Vector transip(x, 3);

miniapps/shifted/distance.cpp

@@ -138,7 +138,7 @@ public:
    ExactDistSphereLoc(ParGridFunction &d)
       : dist(d), dx(dist.ParFESpace()->GetParMesh()->GetElementSize(0)) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       Vector pos(T.GetDimension());
       T.Transform(ip, pos);

miniapps/shifted/extrapolator.hpp

@@ -93,7 +93,7 @@ public:
    using VectorCoefficient::Eval;
 
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       Vector grad_ls(vdim), n(vdim);
       ls_gf.GetGradient(T, grad_ls);
@@ -116,7 +116,7 @@ private:
 public:
    GradComponentCoeff(const ParGridFunction &u, int c) : u_gf(u), comp(c) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       Vector grad_u(T.GetDimension());
       u_gf.GetGradient(T, grad_u);
@@ -134,7 +134,7 @@ public:
    NormalGradCoeff(const ParGridFunction &u, VectorCoefficient &n)
       : u_gf(u), n_coeff(n) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       const int dim = T.GetDimension();
       Vector n(dim), grad_u(dim);
@@ -155,7 +155,7 @@ public:
                             const ParGridFunction &dy, VectorCoefficient &n)
       : du_dx(dx), du_dy(dy), n_coeff(n) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       const int dim = T.GetDimension();
       Vector n(dim), grad_u(dim);

miniapps/shifted/sbm_aux.hpp

@@ -143,7 +143,7 @@ public:
    Dist_Level_Set_Coefficient(int type_)
       : Coefficient(), type(type_) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       Vector x(3);
       T.Transform(ip, x);
@@ -166,7 +166,7 @@ public:
 
    int GetNLevelSets() { return dls.Size(); }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       MFEM_VERIFY(dls.Size() > 0,
                   "Add at least 1 Dist_level_Set_Coefficient to the Combo.");
@@ -192,7 +192,7 @@ public:
    using VectorCoefficient::Eval;
 
    virtual void Eval(Vector &p, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       Vector x;
       T.Transform(ip, x);

miniapps/shifted/sbm_solver.cpp

@@ -16,7 +16,7 @@ namespace mfem
 
 double ShiftedFunctionCoefficient::Eval(ElementTransformation & T,
                                         const IntegrationPoint & ip,
-                                        const Vector &D)
+                                        const Vector &D) const
 {
    if (constantcoefficient) { return constant; }
 
@@ -33,7 +33,7 @@ double ShiftedFunctionCoefficient::Eval(ElementTransformation & T,
 void ShiftedVectorFunctionCoefficient::Eval(Vector &V,
                                             ElementTransformation & T,
                                             const IntegrationPoint & ip,
-                                            const Vector &D)
+                                            const Vector &D) const
 {
    Vector transip;
    T.Transform(ip, transip);

miniapps/shifted/sbm_solver.hpp

@@ -33,7 +33,7 @@ public:
    ShiftedFunctionCoefficient(double constant_)
       : constant(constant_), constantcoefficient(true) { }
 
-   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) const
    {
       if (constantcoefficient) { return constant; }
 
@@ -45,7 +45,7 @@ public:
    /// Evaluate the coefficient at @a ip + @a D.
    double Eval(ElementTransformation &T,
                const IntegrationPoint &ip,
-               const Vector &D);
+               const Vector &D) const;
 };
 
 class ShiftedVectorFunctionCoefficient : public VectorCoefficient
@@ -60,7 +60,7 @@ public:
 
    using VectorCoefficient::Eval;
    virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+                     const IntegrationPoint &ip) const
    {
       Vector D(vdim);
       D = 0.;
@@ -71,7 +71,7 @@ public:
    void Eval(Vector &V,
              ElementTransformation &T,
              const IntegrationPoint &ip,
-             const Vector &D);
+             const Vector &D) const;
 };
 
 /// BilinearFormIntegrator for the high-order extension of shifted boundary

miniapps/tools/display-basis.cpp

@@ -79,12 +79,13 @@ public:
    Deformation(int dim, DefType dType, const DeformationData & data)
       : VectorCoefficient(dim), dim_(dim), dType_(dType), data_(data) {}
 
-   void Eval(Vector &v, ElementTransformation &T, const IntegrationPoint &ip);
+   void Eval(Vector &v, ElementTransformation &T,
+             const IntegrationPoint &ip) const;
    using VectorCoefficient::Eval;
 private:
-   void Def1D(const Vector & u, Vector & v);
-   void Def2D(const Vector & u, Vector & v);
-   void Def3D(const Vector & u, Vector & v);
+   void Def1D(const Vector & u, Vector & v) const;
+   void Def2D(const Vector & u, Vector & v) const;
+   void Def3D(const Vector & u, Vector & v) const;
 
    int     dim_;
    DefType dType_;
@@ -623,7 +624,7 @@ mapTypeStr(int mType)
 
 void
 Deformation::Eval(Vector &v, ElementTransformation &T,
-                  const IntegrationPoint &ip)
+                  const IntegrationPoint &ip) const
 {
    Vector u(dim_);
    T.Transform(ip, u);
@@ -643,7 +644,7 @@ Deformation::Eval(Vector &v, ElementTransformation &T,
 }
 
 void
-Deformation::Def1D(const Vector & u, Vector & v)
+Deformation::Def1D(const Vector & u, Vector & v) const
 {
    v = u;
    if ( dType_ == UNIFORM )
@@ -653,7 +654,7 @@ Deformation::Def1D(const Vector & u, Vector & v)
 }
 
 void
-Deformation::Def2D(const Vector & u, Vector & v)
+Deformation::Def2D(const Vector & u, Vector & v) const
 {
    switch (dType_)
    {
@@ -676,7 +677,7 @@ Deformation::Def2D(const Vector & u, Vector & v)
 }
 
 void
-Deformation::Def3D(const Vector & u, Vector & v)
+Deformation::Def3D(const Vector & u, Vector & v) const
 {
    switch (dType_)
    {

miniapps/tools/nodal-transfer.cpp

@@ -47,7 +47,7 @@ public:
 
    virtual
    double Eval(ElementTransformation &T,
-               const IntegrationPoint &ip)
+               const IntegrationPoint &ip) const
    {
       if (T.GetSpaceDim()==3)
       {