update removeinternalboundariesnotadjacentto

Quadrature function size fix

config/config.hpp

@@ -34,6 +34,12 @@
 // Macro needed to get defines like M_PI from <cmath>. (Visual Studio C++ only?)
 #define _USE_MATH_DEFINES
 #endif
+// On Cygwin the option -std=c++11 prevents the definition of M_PI. Defining
+// the following macro allows us to get M_PI and some needed functions, e.g.
+// posix_memalign(), strdup(), strerror_r().
+#ifdef __CYGWIN__
+#define _XOPEN_SOURCE 600
+#endif
 
 // Check dependencies:
 

fem/bilininteg.hpp

@@ -2050,8 +2050,6 @@ public:
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe,
                                          ElementTransformation &Trans);
-
-   void SetupPA(const FiniteElementSpace &fes);
 };
 
 /** Mass integrator (u, v) restricted to the boundary of a domain */

fem/bilininteg_vecmass.cpp

@@ -106,7 +106,7 @@ void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
 static void PAVectorMassApply2D(const int NE,
-                                const Array<double> &_B,
+                                const Array<double> &B_,
                                 const Array<double> &_Bt,
                                 const Vector &_op,
                                 const Vector &_x,
@@ -119,7 +119,7 @@ static void PAVectorMassApply2D(const int NE,
    constexpr int VDIM = 2;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto Bt = Reshape(_Bt.Read(), D1D, Q1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, NE);
    auto x = Reshape(_x.Read(), D1D, D1D, VDIM, NE);
@@ -203,7 +203,7 @@ static void PAVectorMassApply2D(const int NE,
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
 static void PAVectorMassApply3D(const int NE,
-                                const Array<double> &_B,
+                                const Array<double> &B_,
                                 const Array<double> &_Bt,
                                 const Vector &_op,
                                 const Vector &_x,
@@ -216,7 +216,7 @@ static void PAVectorMassApply3D(const int NE,
    constexpr int VDIM = 3;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto Bt = Reshape(_Bt.Read(), D1D, Q1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
    auto x = Reshape(_x.Read(), D1D, D1D, D1D, VDIM, NE);
@@ -381,7 +381,7 @@ void VectorMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
 static void PAVectorMassAssembleDiagonal2D(const int NE,
-                                           const Array<double> &_B,
+                                           const Array<double> &B_,
                                            const Array<double> &_Bt,
                                            const Vector &_op,
                                            Vector &_diag,
@@ -393,7 +393,7 @@ static void PAVectorMassAssembleDiagonal2D(const int NE,
    constexpr int VDIM = 2;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, NE);
    auto y = Reshape(_diag.ReadWrite(), D1D, D1D, VDIM, NE);
    MFEM_FORALL(e, NE,
@@ -433,7 +433,7 @@ static void PAVectorMassAssembleDiagonal2D(const int NE,
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
 static void PAVectorMassAssembleDiagonal3D(const int NE,
-                                           const Array<double> &_B,
+                                           const Array<double> &B_,
                                            const Array<double> &_Bt,
                                            const Vector &_op,
                                            Vector &_diag,
@@ -445,7 +445,7 @@ static void PAVectorMassAssembleDiagonal3D(const int NE,
    constexpr int VDIM = 3;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
    auto y = Reshape(_diag.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
    MFEM_FORALL(e, NE,

fem/coefficient.cpp

@@ -485,32 +485,32 @@ VectorSumCoefficient::VectorSumCoefficient(int dim)
 }
 
 VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &_A,
-                                           VectorCoefficient &_B,
+                                           VectorCoefficient &B_,
                                            double _alpha, double _beta)
    : VectorCoefficient(_A.GetVDim()),
-     ACoef(&_A), BCoef(&_B),
+     ACoef(&_A), BCoef(&B_),
      A(_A.GetVDim()), B(_A.GetVDim()),
      alphaCoef(NULL), betaCoef(NULL),
      alpha(_alpha), beta(_beta)
 {
-   MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
+   MFEM_ASSERT(_A.GetVDim() == B_.GetVDim(),
                "VectorSumCoefficient:  "
                "Arguments must have the same dimension.");
 }
 
 VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &_A,
-                                           VectorCoefficient &_B,
+                                           VectorCoefficient &B_,
                                            Coefficient &_alpha,
                                            Coefficient &_beta)
    : VectorCoefficient(_A.GetVDim()),
-     ACoef(&_A), BCoef(&_B),
+     ACoef(&_A), BCoef(&B_),
      A(_A.GetVDim()),
      B(_A.GetVDim()),
      alphaCoef(&_alpha),
      betaCoef(&_beta),
      alpha(0.0), beta(0.0)
 {
-   MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
+   MFEM_ASSERT(_A.GetVDim() == B_.GetVDim(),
                "VectorSumCoefficient:  "
                "Arguments must have the same dimension.");
 }
@@ -697,6 +697,11 @@ void OuterProductCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
    }
 }
 
+CrossCrossCoefficient::CrossCrossCoefficient(double A, VectorCoefficient &K)
+   : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(A), a(NULL), k(&K),
+     vk(K.GetVDim())
+{}
+
 CrossCrossCoefficient::CrossCrossCoefficient(Coefficient &A,
                                              VectorCoefficient &K)
    : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(0.0), a(&A), k(&K),

fem/coefficient.hpp

@@ -1171,8 +1171,8 @@ public:
                         double _alpha = 1.0, double _beta = 1.0);
 
    /** Constructor with scalar coefficients.
-       Result is _alpha * _A + _beta * _B */
-   VectorSumCoefficient(VectorCoefficient &_A, VectorCoefficient &_B,
+       Result is _alpha * _A + _beta * B_ */
+   VectorSumCoefficient(VectorCoefficient &_A, VectorCoefficient &B_,
                         Coefficient &_alpha, Coefficient &_beta);
 
    /// Reset the first vector coefficient
@@ -1201,7 +1201,7 @@ public:
    const Vector & GetA() const { return A; }
 
    /// Reset the second vector as a constant
-   void SetB(const Vector &_B) { B = _B; BCoef = NULL; }
+   void SetB(const Vector &B_) { B = B_; BCoef = NULL; }
    /// Return the second vector constant
    const Vector & GetB() const { return B; }
 

fem/conduitdatacollection.cpp

@@ -981,6 +981,11 @@ ConduitDataCollection::SaveRootFile(int num_domains,
    n_root["file_pattern"]     = MeshFilePattern(relay_protocol);
    n_root["tree_pattern"]     = "";
 
+   // Add the time, time step, and cycle
+   n_root["blueprint_index/mesh/state/time"] = time;
+   n_root["blueprint_index/mesh/state/time_step"] = time_step;
+   n_root["blueprint_index/mesh/state/cycle"] = cycle;
+
    relay::io::save(n_root, RootFileName(), root_proto);
 }
 

fem/convergence.cpp

@@ -110,7 +110,8 @@ void ConvergenceStudy::AddL2Error(GridFunction *gf,
 
 void ConvergenceStudy::AddGf(GridFunction *gf, Coefficient *scalar_u,
                              VectorCoefficient *grad,
-                             Coefficient *ell_coeff, double Nu)
+                             Coefficient *ell_coeff,
+                             JumpScaling jump_scaling)
 {
    cont_type = gf->FESpace()->FEColl()->GetContType();
 
@@ -140,7 +141,7 @@ void ConvergenceStudy::AddGf(GridFunction *gf, Coefficient *scalar_u,
 
    if (cont_type == mfem::FiniteElementCollection::DISCONTINUOUS && ell_coeff)
    {
-      double DGErr = gf->ComputeDGFaceJumpError(scalar_u,ell_coeff,Nu);
+      double DGErr = gf->ComputeDGFaceJumpError(scalar_u,ell_coeff,jump_scaling);
       DGFaceErrors.Append(DGErr);
       // Compute the rate of convergence by:
       // rate = log (||u - u_h|| / ||u - u_{h/2}||)/log(2)
@@ -270,26 +271,26 @@ void ConvergenceStudy::Print(bool relative, std::ostream &out)
             }
             out << "\n";
          }
-         if (cont_type == 3 && fcounter)
+      }
+      if (cont_type == 3 && fcounter)
+      {
+         out << " -------------------------------------------" << "\n";
+         out << "            DG Face Jump Error          " << "\n";
+         out << " -------------------------------------------"
+             << "\n";
+         out << std::right<< std::setw(11)<< "DOFs "<< std::setw(13);
+         out << "Error ";
+         out << std::setw(15) << "Rate " << "\n";
+         out << " -------------------------------------------"
+             << "\n";
+         out << std::setprecision(4);
+         for (int i =0; i<fcounter; i++)
          {
-            out << " -------------------------------------------" << "\n";
-            out << "            DG Face Jump Error          " << "\n";
-            out << " -------------------------------------------"
-                << "\n";
-            out << std::right<< std::setw(11)<< "DOFs "<< std::setw(13);
-            out << "Error ";
-            out << std::setw(15) << "Rate " << "\n";
-            out << " -------------------------------------------"
-                << "\n";
-            out << std::setprecision(4);
-            for (int i =0; i<fcounter; i++)
-            {
-               out << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
-                   << std::scientific << DGFaceErrors[i] << std::setw(13)
-                   << std::fixed << DGFaceRates[i] << "\n";
-            }
-            out << "\n";
+            out << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
+                << std::scientific << DGFaceErrors[i] << std::setw(13)
+                << std::fixed << DGFaceRates[i] << "\n";
          }
+         out << "\n";
       }
    }
 }

fem/convergence.hpp

@@ -59,7 +59,8 @@ private:
                    VectorCoefficient *vector_u);
    void AddGf(GridFunction *gf, Coefficient *scalar_u,
               VectorCoefficient *grad=nullptr,
-              Coefficient *ell_coeff=nullptr, double Nu=1.0);
+              Coefficient *ell_coeff=nullptr,
+              JumpScaling jump_scaling = {1.0, JumpScaling::ONE_OVER_H});
    void AddGf(GridFunction *gf, VectorCoefficient *vector_u,
               VectorCoefficient *curl, Coefficient *div);
    // returns the L2-norm of scalar_u or vector_u
@@ -75,9 +76,10 @@ public:
    /// DG face jumps parameters
    void AddL2GridFunction(GridFunction *gf, Coefficient *scalar_u,
                           VectorCoefficient *grad=nullptr,
-                          Coefficient *ell_coeff=nullptr, double Nu=1.0)
+                          Coefficient *ell_coeff=nullptr,
+                          JumpScaling jump_scaling = {1.0, JumpScaling::ONE_OVER_H})
    {
-      AddGf(gf, scalar_u, grad, ell_coeff, Nu);
+      AddGf(gf, scalar_u, grad, ell_coeff, jump_scaling);
    }
 
    /// Add H1 GridFunction, the exact solution and possibly its gradient

fem/gridfunc.cpp

@@ -2777,10 +2777,11 @@ double GridFunction::ComputeDivError(
 }
 
 double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
-                                            Coefficient *ell_coeff, double Nu,
+                                            Coefficient *ell_coeff,
+                                            class JumpScaling jump_scaling,
                                             const IntegrationRule *irs[])  const
 {
-   int fdof, dim, intorder, k;
+   int fdof, intorder, k;
    Mesh *mesh;
    const FiniteElement *fe;
    ElementTransformation *transf;
@@ -2791,20 +2792,24 @@ double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
    double error = 0.0;
 
    mesh = fes->GetMesh();
-   dim = mesh->Dimension();
 
    for (int i = 0; i < mesh->GetNumFaces(); i++)
    {
-      face_elem_transf = mesh->GetFaceElementTransformations(i, 5);
-      int i1 = face_elem_transf->Elem1No;
-      int i2 = face_elem_transf->Elem2No;
+      int i1, i2;
+      mesh->GetFaceElements(i, &i1, &i2);
+      double h = mesh->GetElementSize(i1);
       intorder = fes->GetFE(i1)->GetOrder();
       if (i2 >= 0)
+      {
          if ( (k = fes->GetFE(i2)->GetOrder()) > intorder )
          {
             intorder = k;
          }
+         h = std::min(h, mesh->GetElementSize(i2));
+      }
+      int p = intorder;
       intorder = 2 * intorder;  // <-------------
+      face_elem_transf = mesh->GetFaceElementTransformations(i, 5);
       const IntegrationRule *ir;
       if (irs)
       {
@@ -2875,8 +2880,9 @@ double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
       {
          const IntegrationPoint &ip = ir->IntPoint(j);
          transf->SetIntPoint(&ip);
-         error += (ip.weight * Nu * ell_coeff_val(j) *
-                   pow(transf->Weight(), 1.0-1.0/(dim-1)) *
+         double nu = jump_scaling.Eval(h, p);
+         error += (ip.weight * nu * ell_coeff_val(j) *
+                   transf->Weight() *
                    err_val(j) * err_val(j));
       }
    }
@@ -2884,6 +2890,15 @@ double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
    return (error < 0.0) ? -sqrt(-error) : sqrt(error);
 }
 
+double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
+                                            Coefficient *ell_coeff,
+                                            double Nu,
+                                            const IntegrationRule *irs[])  const
+{
+   return ComputeDGFaceJumpError(
+             exsol, ell_coeff, {Nu, JumpScaling::ONE_OVER_H}, irs);
+}
+
 double GridFunction::ComputeH1Error(Coefficient *exsol,
                                     VectorCoefficient *exgrad,
                                     Coefficient *ell_coef, double Nu,
@@ -2892,7 +2907,11 @@ double GridFunction::ComputeH1Error(Coefficient *exsol,
    double error1 = 0.0;
    double error2 = 0.0;
    if (norm_type & 1) { error1 = GridFunction::ComputeGradError(exgrad); }
-   if (norm_type & 2) { error2 = GridFunction::ComputeDGFaceJumpError(exsol,ell_coef,Nu); }
+   if (norm_type & 2)
+   {
+      error2 = GridFunction::ComputeDGFaceJumpError(
+                  exsol, ell_coef, {Nu, JumpScaling::ONE_OVER_H});
+   }
 
    return sqrt(error1 * error1 + error2 * error2);
 }
@@ -3670,7 +3689,7 @@ QuadratureFunction & QuadratureFunction::operator=(double value)
 
 QuadratureFunction & QuadratureFunction::operator=(const Vector &v)
 {
-   MFEM_ASSERT(qspace && v.Size() == qspace->GetSize(), "");
+   MFEM_ASSERT(qspace && v.Size() == this->Size(), "");
    Vector::operator=(v);
    return *this;
 }

fem/gridfunc.hpp

@@ -451,7 +451,17 @@ public:
    virtual double ComputeDivError(Coefficient *exdiv,
                                   const IntegrationRule *irs[] = NULL) const;
 
-   /// Returns the Face Jumps error for L2 elements
+   /// Returns the Face Jumps error for L2 elements. The error can be weighted
+   /// by a constant nu, by nu/h, or nu*p^2/h, depending on the value of
+   /// @a jump_scaling.
+   virtual double ComputeDGFaceJumpError(Coefficient *exsol,
+                                         Coefficient *ell_coeff,
+                                         class JumpScaling jump_scaling,
+                                         const IntegrationRule *irs[] = NULL)
+   const;
+
+   /// Returns the Face Jumps error for L2 elements, with 1/h scaling.
+   MFEM_DEPRECATED
    virtual double ComputeDGFaceJumpError(Coefficient *exsol,
                                          Coefficient *ell_coeff,
                                          double Nu,
@@ -664,6 +674,32 @@ public:
     derived class ParGridFunction */
 std::ostream &operator<<(std::ostream &out, const GridFunction &sol);
 
+/// Class used to specify how the jump terms in
+/// GridFunction::ComputeDGFaceJumpError are scaled.
+class JumpScaling
+{
+public:
+   enum JumpScalingType
+   {
+      CONSTANT,
+      ONE_OVER_H,
+      P_SQUARED_OVER_H
+   };
+private:
+   double nu;
+   JumpScalingType type;
+public:
+   JumpScaling(double nu_=1.0, JumpScalingType type_=CONSTANT)
+      : nu(nu_), type(type_) { }
+   double Eval(double h, int p) const
+   {
+      double val = nu;
+      if (type != CONSTANT) { val /= h; }
+      if (type == P_SQUARED_OVER_H) { val *= p*p; }
+      return val;
+   }
+};
+
 
 /** @brief Class representing a function through its values (scalar or vector)
     at quadrature points. */
@@ -752,7 +788,8 @@ public:
 
    /// Copy the data from @a v.
    /** The size of @a v must be equal to the size of the associated
-       QuadratureSpace #qspace. */
+       QuadratureSpace #qspace times the QuadratureFunction dimension
+       i.e. QuadratureFunction::Size(). */
    QuadratureFunction &operator=(const Vector &v);
 
    /// Copy assignment. Only the data of the base class Vector is copied.

fem/lininteg.cpp

@@ -10,8 +10,8 @@
 // CONTRIBUTING.md for details.
 
 
-#include <cmath>
 #include "fem.hpp"
+#include <cmath>
 
 namespace mfem
 {

fem/pgridfunc.cpp

@@ -657,12 +657,12 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
 
 double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
                                                Coefficient *ell_coeff,
-                                               double Nu,
+                                               JumpScaling jump_scaling,
                                                const IntegrationRule *irs[]) const
 {
    const_cast<ParGridFunction *>(this)->ExchangeFaceNbrData();
 
-   int fdof, dim, intorder, k;
+   int fdof, intorder, k;
    ElementTransformation *transf;
    Vector shape, el_dofs, err_val, ell_coeff_val;
    Array<int> vdofs;
@@ -670,7 +670,6 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
    double error = 0.0;
 
    ParMesh *mesh = pfes->GetParMesh();
-   dim = mesh->Dimension();
 
    std::map<int,int> local_to_shared;
    for (int i = 0; i < mesh->GetNSharedFaces(); ++i)
@@ -687,6 +686,7 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
       mesh->GetFaceElements(i, &iel1, &iel2);
       mesh->GetFaceInfos(i, &info1, &info2);
 
+      double h = mesh->GetElementSize(iel1);
       intorder = fes->GetFE(iel1)->GetOrder();
 
       FaceElementTransformations *face_elem_transf;
@@ -703,11 +703,10 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
          }
          shared_face = true;
          shared_face_factor = 0.5;
+         h = std::min(h, mesh->GetFaceNbrElementSize(iel2));
       }
       else
       {
-         face_elem_transf = mesh->GetFaceElementTransformations(i);
-
          if (iel2 >= 0)
          {
             fe2 = pfes->GetFE(iel2);
@@ -715,12 +714,15 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
             {
                intorder = k;
             }
+            h = std::min(h, mesh->GetElementSize(iel2));
          }
          else
          {
             fe2 = NULL;
          }
+         face_elem_transf = mesh->GetFaceElementTransformations(i);
       }
+      int p = intorder;
 
       intorder = 2 * intorder;  // <-------------
       const IntegrationRule *ir;
@@ -806,8 +808,9 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
       {
          const IntegrationPoint &ip = ir->IntPoint(j);
          transf->SetIntPoint(&ip);
-         error += shared_face_factor*(ip.weight * Nu * ell_coeff_val(j) *
-                                      pow(transf->Weight(), 1.0-1.0/(dim-1)) *
+         double nu = jump_scaling.Eval(h, p);
+         error += shared_face_factor*(ip.weight * nu * ell_coeff_val(j) *
+                                      transf->Weight() *
                                       err_val(j) * err_val(j));
       }
    }

fem/pgridfunc.hpp

@@ -310,7 +310,7 @@ public:
    /// Returns the Face Jumps error for L2 elements
    virtual double ComputeDGFaceJumpError(Coefficient *exsol,
                                          Coefficient *ell_coeff,
-                                         double Nu,
+                                         JumpScaling jump_scaling,
                                          const IntegrationRule *irs[]=NULL)
    const;
 

general/mem_manager.cpp

@@ -482,7 +482,8 @@ public:
       HostMemorySpace(),
       name(mm.GetUmpireAllocatorHostName()),
       rm(umpire::ResourceManager::getInstance()),
-      h_allocator(rm.isAllocator(name)? rm.getAllocator(name):
+      h_allocator((!std::strcmp(name, "HOST") || rm.isAllocator(name)) ?
+                  rm.getAllocator(name) :
                   rm.makeAllocator<umpire::strategy::DynamicPool>
                   (name, rm.getAllocator("HOST"))),
       strat(h_allocator.getAllocationStrategy()) { }
@@ -506,7 +507,8 @@ public:
       DeviceMemorySpace(),
       name(mm.GetUmpireAllocatorDeviceName()),
       rm(umpire::ResourceManager::getInstance()),
-      d_allocator(rm.isAllocator(name)? rm.getAllocator(name):
+      d_allocator((!std::strcmp(name, "DEVICE") || rm.isAllocator(name)) ?
+                  rm.getAllocator(name) :
                   rm.makeAllocator<umpire::strategy::DynamicPool>
                   (name, rm.getAllocator("DEVICE"))) { }
    void Alloc(Memory &base) { base.d_ptr = d_allocator.allocate(base.bytes); }

general/text.hpp

@@ -12,6 +12,7 @@
 #ifndef MFEM_TEXT
 #define MFEM_TEXT
 
+#include "../config/config.hpp"
 #include <istream>
 #include <iomanip>
 #include <sstream>
@@ -24,6 +25,8 @@ namespace mfem
 
 // Utilities for text parsing
 
+using std::to_string;
+
 /// Check if the stream starts with @a comment_char. If so skip it.
 inline void skip_comment_lines(std::istream &is, const char comment_char)
 {
@@ -47,18 +50,6 @@ inline void filter_dos(std::string &line)
    }
 }
 
-/// Convert an integer to an std::string.
-inline std::string to_string(int i)
-{
-   std::stringstream ss;
-   ss << i;
-
-   // trim leading spaces
-   std::string out_str = ss.str();
-   out_str = out_str.substr(out_str.find_first_not_of(" \t"));
-   return out_str;
-}
-
 /// Convert an integer to a 0-padded string with the given number of @a digits
 inline std::string to_padded_string(int i, int digits)
 {

linalg/kernels.hpp

@@ -12,11 +12,6 @@
 #ifndef MFEM_KERNELS_HPP
 #define MFEM_KERNELS_HPP
 
-#ifdef _WIN32
-#define _USE_MATH_DEFINES
-#include <cmath>
-#endif
-
 #include "../config/config.hpp"
 #include "../general/backends.hpp"
 #include "../general/globals.hpp"

linalg/vector.cpp

@@ -22,6 +22,10 @@
 #endif
 #endif
 
+#ifdef MFEM_USE_OPENMP
+#include <omp.h>
+#endif
+
 #include <iostream>
 #include <iomanip>
 #include <cmath>
@@ -1076,6 +1080,30 @@ double Vector::operator*(const Vector &v) const
 #ifdef MFEM_USE_OPENMP
    if (Device::Allows(Backend::OMP_MASK))
    {
+#define MFEM_USE_OPENMP_DETERMINISTIC_DOT
+#ifdef MFEM_USE_OPENMP_DETERMINISTIC_DOT
+      // By default, use a deterministic way of computing the dot product
+      static Vector th_dot;
+      #pragma omp parallel
+      {
+         const int nt = omp_get_num_threads();
+         #pragma omp master
+         th_dot.SetSize(nt);
+         const int tid    = omp_get_thread_num();
+         const int stride = (size + nt - 1)/nt;
+         const int start  = tid*stride;
+         const int stop   = std::min(start + stride, size);
+         double my_dot = 0.0;
+         for (int i = start; i < stop; i++)
+         {
+            my_dot += m_data[i] * v_data[i];
+         }
+         #pragma omp barrier
+         th_dot(tid) = my_dot;
+      }
+      return th_dot.Sum();
+#else
+      // The standard way of computing the dot product is non-deterministic
       double prod = 0.0;
       #pragma omp parallel for reduction(+:prod)
       for (int i = 0; i < size; i++)
@@ -1083,8 +1111,9 @@ double Vector::operator*(const Vector &v) const
          prod += m_data[i] * v_data[i];
       }
       return prod;
+#endif // MFEM_USE_OPENMP_DETERMINISTIC_DOT
    }
-#endif
+#endif // MFEM_USE_OPENMP
    if (Device::Allows(Backend::DEBUG_DEVICE))
    {
       const int N = size;

linalg/vector.hpp

@@ -42,10 +42,19 @@ namespace mfem
 inline int CheckFinite(const double *v, const int n);
 
 /// Define a shortcut for std::numeric_limits<double>::infinity()
+#ifndef __CYGWIN__
 inline double infinity()
 {
    return std::numeric_limits<double>::infinity();
 }
+#else
+// On Cygwin math.h defines a function 'infinity()' which will conflict with the
+// above definition if we have 'using namespace mfem;' and try to use something
+// like 'double a = infinity();'. This 'infinity()' function is non-standard and
+// is defined by the Newlib C standard library implementation used by Cygwin,
+// see https://en.wikipedia.org/wiki/Newlib, http://www.sourceware.org/newlib.
+using ::infinity;
+#endif
 
 /// Vector data type.
 class Vector

mesh/mesh.cpp

@@ -29,6 +29,7 @@
 #include <cstring>
 #include <ctime>
 #include <functional>
+#include <unordered_set>
 
 // Include the METIS header, if using version 5. If using METIS 4, the needed
 // declarations are inlined below, i.e. no header is needed.
@@ -71,10 +72,14 @@ void Mesh::GetElementCenter(int i, Vector &center)
    eltransf->Transform(Geometries.GetCenter(geom), center);
 }
 
-double Mesh::GetElementSize(int i, int type)
+double Mesh::GetElementSize(ElementTransformation *T, int type)
 {
    DenseMatrix J(Dim);
-   GetElementJacobian(i, J);
+
+   Geometry::Type geom = T->GetGeometryType();
+   T->SetIntPoint(&Geometries.GetCenter(geom));
+   Geometries.JacToPerfJac(geom, T->Jacobian(), J);
+
    if (type == 0)
    {
       return pow(fabs(J.Det()), 1./Dim);
@@ -89,6 +94,11 @@ double Mesh::GetElementSize(int i, int type)
    }
 }
 
+double Mesh::GetElementSize(int i, int type)
+{
+   return GetElementSize(GetElementTransformation(i), type);
+}
+
 double Mesh::GetElementSize(int i, const Vector &dir)
 {
    DenseMatrix J(Dim);
@@ -1463,6 +1473,13 @@ void Mesh::AddBdrQuadAsTriangles(const int *vi, int attr)
    }
 }
 
+int Mesh::AddBdrPoint(int v, int attr)
+{
+   CheckEnlarge(boundary, NumOfBdrElements);
+   boundary[NumOfBdrElements] = new Point(&v, attr);
+   return NumOfBdrElements++;
+}
+
 void Mesh::GenerateBoundaryElements()
 {
    int i, j;
@@ -10533,6 +10550,227 @@ void Mesh::RemoveInternalBoundaries()
    attribs.Copy(bdr_attributes);
 }
 
+void Mesh::RemoveInternalBoundaries(const Array<int> &_keep)
+{
+   if (NURBSext || ncmesh) { return; }
+
+   std::unordered_set<int> keep(_keep.GetData(),
+                                _keep.GetData() + _keep.Size());
+
+   int num_bdr_elem = 0;
+   int new_bel_to_edge_nnz = 0;
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      if (FaceIsInterior(GetBdrElementEdgeIndex(i)) && 
+         (keep.count(boundary[i]->GetAttribute()) == 0))
+      {
+         FreeElement(boundary[i]);
+      }
+      else
+      {
+         num_bdr_elem++;
+         if (Dim == 3)
+         {
+            new_bel_to_edge_nnz += bel_to_edge->RowSize(i);
+         }
+      }
+   }
+
+   if (num_bdr_elem == GetNBE()) { return; }
+
+   Array<Element *> new_boundary(num_bdr_elem);
+   Array<int> new_be_to_edge, new_be_to_face;
+   Table *new_bel_to_edge = NULL;
+   new_boundary.SetSize(0);
+   if (Dim == 2)
+   {
+      new_be_to_edge.Reserve(num_bdr_elem);
+   }
+   else if (Dim == 3)
+   {
+      new_be_to_face.Reserve(num_bdr_elem);
+      new_bel_to_edge = new Table;
+      new_bel_to_edge->SetDims(num_bdr_elem, new_bel_to_edge_nnz);
+   }
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      /// if it's not interior or is in the keep boundary list
+      if (!FaceIsInterior(GetBdrElementEdgeIndex(i)) || 
+          keep.count(boundary[i]->GetAttribute()))
+      {
+         new_boundary.Append(boundary[i]);
+         if (Dim == 2)
+         {
+            new_be_to_edge.Append(be_to_edge[i]);
+         }
+         else if (Dim == 3)
+         {
+            int row = new_be_to_face.Size();
+            new_be_to_face.Append(be_to_face[i]);
+            int *e = bel_to_edge->GetRow(i);
+            int ne = bel_to_edge->RowSize(i);
+            int *new_e = new_bel_to_edge->GetRow(row);
+            for (int j = 0; j < ne; j++)
+            {
+               new_e[j] = e[j];
+            }
+            new_bel_to_edge->GetI()[row+1] = new_bel_to_edge->GetI()[row] + ne;
+         }
+      }
+   }
+
+   NumOfBdrElements = new_boundary.Size();
+   mfem::Swap(boundary, new_boundary);
+
+   if (Dim == 2)
+   {
+      mfem::Swap(be_to_edge, new_be_to_edge);
+   }
+   else if (Dim == 3)
+   {
+      mfem::Swap(be_to_face, new_be_to_face);
+      delete bel_to_edge;
+      bel_to_edge = new_bel_to_edge;
+   }
+
+   Array<int> attribs(num_bdr_elem);
+   for (int i = 0; i < attribs.Size(); i++)
+   {
+      attribs[i] = GetBdrAttribute(i);
+   }
+   attribs.Sort();
+   attribs.Unique();
+   bdr_attributes.DeleteAll();
+   attribs.Copy(bdr_attributes);
+}
+
+void Mesh::RemoveInternalBoundariesNotAdjacentTo(const Array<int> &_regions)
+{
+   if (NURBSext || ncmesh) { return; }
+
+   std::unordered_set<int> regions(_regions.GetData(),
+                                   _regions.GetData() + _regions.Size());
+
+   int num_bdr_elem = 0;
+   int new_bel_to_edge_nnz = 0;
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      auto faceNo = GetBdrElementEdgeIndex(i);
+      bool interior = FaceIsInterior(faceNo);
+
+      auto face_info = faces_info[faceNo];
+      auto adj = regions.count(elements[face_info.Elem1No]->GetAttribute());
+      if (face_info.Elem2No > 0)
+         adj += regions.count(elements[face_info.Elem2No]->GetAttribute());
+
+      bool remove = interior & !adj;
+
+      // if (remove)
+      // {
+      //    std::cout << "remove face on: " << boundary[i]->GetAttribute()
+      //              << " adj to ("
+      //              << elements[face_info.Elem1No]->GetAttribute();
+      //    if (face_info.Elem2No > 0)
+      //    {
+      //       std::cout << ", "
+      //                 << elements[face_info.Elem2No]->GetAttribute() << ")\n";
+      //    }
+      //    else
+      //    {
+      //       std::cout << ")\n";
+      //    }
+      // }
+
+      if (remove)
+      {
+         FreeElement(boundary[i]);
+      }
+      else
+      {
+         num_bdr_elem++;
+         if (Dim == 3)
+         {
+            new_bel_to_edge_nnz += bel_to_edge->RowSize(i);
+         }
+      }
+   }
+
+   if (num_bdr_elem == GetNBE()) { return; }
+
+   Array<Element *> new_boundary(num_bdr_elem);
+   Array<int> new_be_to_edge, new_be_to_face;
+   Table *new_bel_to_edge = NULL;
+   new_boundary.SetSize(0);
+   if (Dim == 2)
+   {
+      new_be_to_edge.Reserve(num_bdr_elem);
+   }
+   else if (Dim == 3)
+   {
+      new_be_to_face.Reserve(num_bdr_elem);
+      new_bel_to_edge = new Table;
+      new_bel_to_edge->SetDims(num_bdr_elem, new_bel_to_edge_nnz);
+   }
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      auto faceNo = GetBdrElementEdgeIndex(i);
+      bool interior = FaceIsInterior(faceNo);
+
+      auto face_info = faces_info[faceNo];
+      auto adj = regions.count(elements[face_info.Elem1No]->GetAttribute());
+      if (face_info.Elem2No > 0)
+         adj += regions.count(elements[face_info.Elem2No]->GetAttribute());
+
+      bool keep = !interior || adj;
+
+      if (keep)
+      {
+         new_boundary.Append(boundary[i]);
+         if (Dim == 2)
+         {
+            new_be_to_edge.Append(be_to_edge[i]);
+         }
+         else if (Dim == 3)
+         {
+            int row = new_be_to_face.Size();
+            new_be_to_face.Append(be_to_face[i]);
+            int *e = bel_to_edge->GetRow(i);
+            int ne = bel_to_edge->RowSize(i);
+            int *new_e = new_bel_to_edge->GetRow(row);
+            for (int j = 0; j < ne; j++)
+            {
+               new_e[j] = e[j];
+            }
+            new_bel_to_edge->GetI()[row+1] = new_bel_to_edge->GetI()[row] + ne;
+         }
+      }
+   }
+
+   NumOfBdrElements = new_boundary.Size();
+   mfem::Swap(boundary, new_boundary);
+
+   if (Dim == 2)
+   {
+      mfem::Swap(be_to_edge, new_be_to_edge);
+   }
+   else if (Dim == 3)
+   {
+      mfem::Swap(be_to_face, new_be_to_face);
+      delete bel_to_edge;
+      bel_to_edge = new_bel_to_edge;
+   }
+
+   Array<int> attribs(num_bdr_elem);
+   for (int i = 0; i < attribs.Size(); i++)
+   {
+      attribs[i] = GetBdrAttribute(i);
+   }
+   attribs.Sort();
+   attribs.Unique();
+   bdr_attributes.DeleteAll();
+   attribs.Copy(bdr_attributes);
+}
+
 void Mesh::FreeElement(Element *E)
 {
 #ifdef MFEM_USE_MEMALLOC

mesh/mesh.hpp

@@ -471,6 +471,8 @@ protected:
    void GetElementData(const Array<Element*> &elem_array, int geom,
                        Array<int> &elem_vtx, Array<int> &attr) const;
 
+   double GetElementSize(ElementTransformation *T, int type = 0);
+
 public:
 
    Mesh() { SetEmpty(); }
@@ -555,6 +557,8 @@ public:
    int AddBdrQuad(const int *vi, int attr = 1);
    void AddBdrQuadAsTriangles(const int *vi, int attr = 1);
 
+   int AddBdrPoint(int v, int attr = 1);
+
    void GenerateBoundaryElements();
    /// Finalize the construction of a triangular Mesh.
    void FinalizeTriMesh(int generate_edges = 0, int refine = 0,
@@ -1302,6 +1306,12 @@ public:
        have two adjacent faces in 3D, or edges in 2D. */
    void RemoveInternalBoundaries();
 
+   /** Remove boundary elements that lie in the interior of the mesh, except 
+       for those with boundary attributes in `keep` */
+   void RemoveInternalBoundaries(const Array<int> &keep);
+
+   void RemoveInternalBoundariesNotAdjacentTo(const Array<int> &regions);
+
    /** @brief Get the size of the i-th element relative to the perfect
        reference element. */
    double GetElementSize(int i, int type = 0);

mesh/pmesh.cpp

@@ -1735,6 +1735,11 @@ void ParMesh::GetFaceNbrElementTransformation(
    }
 }
 
+double ParMesh::GetFaceNbrElementSize(int i, int type)
+{
+   return GetElementSize(GetFaceNbrElementTransformation(i), type);
+}
+
 void ParMesh::DeleteFaceNbrData()
 {
    if (!have_face_nbr_data)

mesh/pmesh.hpp

@@ -305,6 +305,10 @@ public:
       return &FaceNbrTransformation;
    }
 
+   /// Get the size of the i-th face neighbor element relative to the reference
+   /// element.
+   double GetFaceNbrElementSize(int i, int type=0);
+
    /// Return the number of shared faces (3D), edges (2D), vertices (1D)
    int GetNSharedFaces() const;
 

tests/convergence/prates.cpp

@@ -26,6 +26,7 @@
 //               mpirun -np 4 prates -m ../../data/inline-hex.mesh -sr 0 -pr 1 -prob 1 -o 3
 //               mpirun -np 4 prates -m ../../data/square-disc.mesh -sr 1 -pr 2 -prob 1 -o 2
 //               mpirun -np 4 prates -m ../../data/star.mesh -sr 1 -pr 2 -prob 3 -o 2
+//               mpirun -np 4 prates -m ../../data/star.mesh -sr 1 -pr 2 -prob 3 -o 2 -j 0
 //               mpirun -np 4 prates -m ../../data/inline-hex.mesh -sr 1 -pr 1 -prob 3 -o 2
 //
 // Description:  This example code demonstrates the use of MFEM to define and
@@ -80,6 +81,7 @@ int main(int argc, char *argv[])
    bool visualization = 1;
    int sr = 1;
    int pr = 1;
+   int jump_scaling_type = 1;
    double sigma = -1.0;
    double kappa = -1.0;
    OptionsParser args(argc, argv);
@@ -95,6 +97,9 @@ int main(int argc, char *argv[])
    args.AddOption(&kappa, "-k", "--kappa",
                   "One of the two DG penalty parameters, should be positive."
                   " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&jump_scaling_type, "-j", "--jump-scaling",
+                  "Scaling of the jump error for DG methods: "
+                  "0: no scaling, 1: 1/h, 2: p^2/h");
    args.AddOption(&sr, "-sr", "--serial_ref",
                   "Number of serial refinements.");
    args.AddOption(&pr, "-pr", "--parallel_ref",
@@ -285,12 +290,15 @@ int main(int argc, char *argv[])
       delete solver;
 
       x = *X;
+      JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
+                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                        : JumpScaling::CONSTANT);
       switch (prob)
       {
          case 0: rates.AddH1GridFunction(&x,scalar_u,gradu); break;
          case 1: rates.AddHcurlGridFunction(&x,vector_u,curlu); break;
          case 2: rates.AddHdivGridFunction(&x,vector_u,divu);  break;
-         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one); break;
+         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one,js); break;
       }
 
       delete X;

tests/convergence/rates.cpp

@@ -26,6 +26,7 @@
 //               rates -m ../../data/inline-hex.mesh -sr 1 -prob 1 -o 2
 //               rates -m ../../data/square-disc.mesh -sr 2 -prob 1 -o 1
 //               rates -m ../../data/star.mesh -sr 2 -prob 3 -o 2
+//               rates -m ../../data/star.mesh -sr 2 -prob 3 -o 2 -j 0
 //               rates -m ../../data/inline-hex.mesh -sr 1 -prob 3 -o 1
 //
 // Description:  This example code demonstrates the use of MFEM to define and
@@ -73,6 +74,7 @@ int main(int argc, char *argv[])
    int order = 1;
    bool visualization = 1;
    int sr = 1;
+   int jump_scaling_type = 1;
    double sigma = -1.0;
    double kappa = -1.0;
    OptionsParser args(argc, argv);
@@ -88,6 +90,9 @@ int main(int argc, char *argv[])
    args.AddOption(&kappa, "-k", "--kappa",
                   "One of the two DG penalty parameters, should be positive."
                   " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&jump_scaling_type, "-j", "--jump-scaling",
+                  "Scaling of the jump error for DG methods: "
+                  "0: no scaling, 1: 1/h, 2: p^2/h");
    args.AddOption(&sr, "-sr", "--serial_ref",
                   "Number of serial refinements.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -218,12 +223,16 @@ int main(int argc, char *argv[])
          PCG(A, M, b, x, 0, 500, 1e-12, 0.0);
       }
 
+      JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
+                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                        : JumpScaling::CONSTANT);
+
       switch (prob)
       {
          case 0: rates.AddH1GridFunction(&x,scalar_u,gradu); break;
          case 1: rates.AddHcurlGridFunction(&x,vector_u,curlu); break;
          case 2: rates.AddHdivGridFunction(&x,vector_u,divu);  break;
-         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one); break;
+         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one,js); break;
       }
 
       if (l==sr) break;

tests/unit/linalg/test_matrix_block.cpp

@@ -10,9 +10,9 @@
 // CONTRIBUTING.md for details.
 
 #include <fstream>
-using namespace std;
-
 #include "mfem.hpp"
+
+using namespace std;
 using namespace mfem;
 
 #include "unit_tests.hpp"

tests/unit/miniapps/test_sedov.cpp

@@ -98,14 +98,14 @@ struct Tensors1D
 
 template<int DIM, int D1D, int Q1D, int L1D, int H1D, int NBZ =1> static
 void kSmemForceMult2D(const int NE,
-                      const Array<double> &_B,
+                      const Array<double> &B_,
                       const Array<double> &_Bt,
                       const Array<double> &_Gt,
                       const DenseTensor &_sJit,
                       const Vector &_e,
                       Vector &_v)
 {
-   auto b = Reshape(_B.Read(), Q1D, L1D);
+   auto b = Reshape(B_.Read(), Q1D, L1D);
    auto bt = Reshape(_Bt.Read(), H1D, Q1D);
    auto gt = Reshape(_Gt.Read(), H1D, Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*NE*2*2),
@@ -241,14 +241,14 @@ void kSmemForceMult2D(const int NE,
 
 template<int DIM, int D1D, int Q1D, int L1D, int H1D> static
 void kSmemForceMult3D(const int NE,
-                      const Array<double> &_B,
+                      const Array<double> &B_,
                       const Array<double> &_Bt,
                       const Array<double> &_Gt,
                       const DenseTensor &_sJit,
                       const Vector &_e,
                       Vector &_v)
 {
-   auto b = Reshape(_B.Read(), Q1D, L1D);
+   auto b = Reshape(B_.Read(), Q1D, L1D);
    auto bt = Reshape(_Bt.Read(), H1D, Q1D);
    auto gt = Reshape(_Gt.Read(), H1D, Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*Q1D*NE*3*3),
@@ -493,14 +493,14 @@ static void kForceMult(const int DIM,
 template<int DIM, int D1D, int Q1D, int L1D, int H1D, int NBZ =1> static
 void kSmemForceMultTranspose2D(const int NE,
                                const Array<double> &_Bt,
-                               const Array<double> &_B,
+                               const Array<double> &B_,
                                const Array<double> &_G,
                                const DenseTensor &_sJit,
                                const Vector &_v,
                                Vector &_e)
 {
    MFEM_VERIFY(D1D==H1D,"");
-   auto b = Reshape(_B.Read(), Q1D,H1D);
+   auto b = Reshape(B_.Read(), Q1D,H1D);
    auto g = Reshape(_G.Read(), Q1D,H1D);
    auto bt = Reshape(_Bt.Read(), L1D,Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*NE*2*2),
@@ -633,14 +633,14 @@ void kSmemForceMultTranspose2D(const int NE,
 template<int DIM, int D1D, int Q1D, int L1D, int H1D> static
 void kSmemForceMultTranspose3D(const int NE,
                                const Array<double> &_Bt,
-                               const Array<double> &_B,
+                               const Array<double> &B_,
                                const Array<double> &_G,
                                const DenseTensor &_sJit,
                                const Vector &_v,
                                Vector &_e)
 {
    MFEM_VERIFY(D1D==H1D,"");
-   auto b = Reshape(_B.Read(), Q1D,H1D);
+   auto b = Reshape(B_.Read(), Q1D,H1D);
    auto g = Reshape(_G.Read(), Q1D,H1D);
    auto bt = Reshape(_Bt.Read(), L1D,Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*Q1D*NE*3*3),