Small bug fix

fem/bilininteg.cpp

@@ -1061,7 +1061,7 @@ void DiffusionIntegrator::ComputeElementFlux
 
 double DiffusionIntegrator::ComputeFluxEnergy
 ( const FiniteElement &fluxelem, ElementTransformation &Trans,
-  Vector &flux, Vector* d_energy)
+  Vector &flux, bool with_coef, Vector* d_energy)
 {
    int nd = fluxelem.GetDof();
    int dim = fluxelem.GetDim();
@@ -1107,18 +1107,38 @@ double DiffusionIntegrator::ComputeFluxEnergy
       if (MQ)
       {
          MQ->Eval(M, Trans, ip);
+         if (with_coef) { M.Invert(); }
          energy += w * M.InnerProduct(pointflux, pointflux);
       }
       else if (VQ)
       {
          VQ->Eval(D, Trans, ip);
+         if (with_coef)
+         {
+            Vector Dinv(D.Size());
+            Dinv = 1.0;
+            Dinv /= D;
+            D = Dinv;
+         }
+
          D *= pointflux;
+
          energy += w * (D * pointflux);
       }
       else
       {
          double e = (pointflux * pointflux);
-         if (Q) { e *= Q->Eval(Trans, ip); }
+         if (Q)
+         {
+            if (with_coef)
+            {
+               e /= Q->Eval(Trans, ip);
+            }
+            else
+            {
+               e *= Q->Eval(Trans, ip);
+            }
+         }
          energy += w * e;
       }
 
@@ -1958,7 +1978,8 @@ void CurlCurlIntegrator
 
 double CurlCurlIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
                                              ElementTransformation &Trans,
-                                             Vector &flux, Vector *d_energy)
+                                             Vector &flux, bool with_coef,
+                                             Vector *d_energy)
 {
    int nd = fluxelem.GetDof();
    int dim = fluxelem.GetDim();
@@ -2851,7 +2872,8 @@ void ElasticityIntegrator::ComputeElementFlux(
 
 double ElasticityIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
                                                ElementTransformation &Trans,
-                                               Vector &flux, Vector *d_energy)
+                                               Vector &flux, bool with_coef,
+                                               Vector *d_energy)
 {
    const int dof = fluxelem.GetDof();
    const int dim = fluxelem.GetDim();

fem/bilininteg.hpp

@@ -234,6 +234,8 @@ public:
                          position of the mesh element.
        @param[in] flux   "Flux" coefficients representing the expansion of the
                          "flux" function in the basis of @a fluxelem.
+       @param[in] wcoef  If true, @a flux includes the coefficient of this
+                         integrator.
        @param[out] d_energy  If not NULL, the given Vector should be set to
                              represent directional energy split that can be used
                              for anisotropic error estimation.
@@ -241,7 +243,8 @@ public:
     */
    virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                     ElementTransformation &Trans,
-                                    Vector &flux, Vector *d_energy = NULL)
+                                    Vector &flux, bool wcoef = true,
+                                    Vector *d_energy = NULL)
    { return 0.0; }
 
    virtual ~BilinearFormIntegrator() { }
@@ -2033,7 +2036,8 @@ public:
 
    virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                     ElementTransformation &Trans,
-                                    Vector &flux, Vector *d_energy = NULL);
+                                    Vector &flux, bool wcoef = true,
+                                    Vector *d_energy = NULL);
 
    using BilinearFormIntegrator::AssemblePA;
 
@@ -2452,7 +2456,8 @@ public:
 
    virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                     ElementTransformation &Trans,
-                                    Vector &flux, Vector *d_energy = NULL);
+                                    Vector &flux, bool wcoef = true,
+                                    Vector *d_energy = NULL);
 
    using BilinearFormIntegrator::AssemblePA;
    virtual void AssemblePA(const FiniteElementSpace &fes);
@@ -2785,7 +2790,8 @@ public:
        s_yz in 3D. */
    virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                     ElementTransformation &Trans,
-                                    Vector &flux, Vector *d_energy = NULL);
+                                    Vector &flux, bool wcoef = true,
+                                    Vector *d_energy = NULL);
 };
 
 /** Integrator for the DG form:

fem/gridfunc.cpp

@@ -4003,7 +4003,7 @@ double ZZErrorEstimator(BilinearFormIntegrator &blfi,
          fl -= fla;
 
          double err = blfi.ComputeFluxEnergy(*ffes->GetFE(i), *Transf, fl,
-                                             (aniso_flags ? &d_xyz : NULL));
+                                             with_coeff, (aniso_flags ? &d_xyz : NULL));
 
          error_estimates(i) = std::sqrt(err);
          total_error += err;