update removeinternalboundariesnotadjacentto

Quadrature function size fix

CMakeLists.txt

@@ -14,7 +14,7 @@ set(USER_CONFIG "${CMAKE_CURRENT_SOURCE_DIR}/config/user.cmake" CACHE PATH
   "Path to optional user configuration file.")
 
 # Require C++11 and disable compiler-specific extensions
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 11)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 

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;
 

miniapps/CMakeLists.txt

@@ -26,4 +26,3 @@ add_subdirectory(tools)
 add_subdirectory(toys)
 add_subdirectory(nurbs)
 add_subdirectory(gslib)
-add_subdirectory(variationalform)

miniapps/variationalform/CMakeLists.txt

@@ -1,15 +0,0 @@
-add_library(parvariationalform parvariationalform.cpp)
-target_include_directories(parvariationalform PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
-target_link_libraries(parvariationalform mfem)
-
-add_executable(parvariationalform_example parvariationalform_example.cpp)
-target_include_directories(parvariationalform_example PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
-target_link_libraries(parvariationalform_example parvariationalform mfem)
-
-add_executable(parvariationalform_ex1 parvariationalform_ex1.cpp)
-target_include_directories(parvariationalform_ex1 PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
-target_link_libraries(parvariationalform_ex1 parvariationalform mfem)
-
-add_executable(parvariationalform_unit_tests test_tensor_ad.cpp)
-target_include_directories(parvariationalform_unit_tests PRIVATE ${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_SOURCE_DIR}/tests/unit)
-target_link_libraries(parvariationalform_unit_tests parvariationalform mfem)

miniapps/variationalform/genericintegrator.hpp

@@ -1,23 +0,0 @@
-#include "mfem.hpp"
-
-#pragma once
-
-namespace mfem
-{
-class GenericIntegrator
-{
-public:
-   GenericIntegrator(const IntegrationRule *ir = nullptr) : IntRule(ir) {}
-
-   virtual void Setup(const FiniteElementSpace &) = 0;
-
-   virtual void Apply(const Vector &x, Vector &y) const = 0;
-
-   virtual void ApplyGradient(const Vector &x,
-                              const Vector &v,
-                              Vector &y) const = 0;
-
-protected:
-   const IntegrationRule *IntRule;
-};
-} // namespace mfem

miniapps/variationalform/parvariationalform.cpp

@@ -1,94 +0,0 @@
-#include "parvariationalform.hpp"
-#include "qfuncintegrator.hpp"
-
-namespace mfem
-{
-ParVariationalForm::ParVariationalForm(ParFiniteElementSpace *f)
-   : Operator(f->GetTrueVSize()), fes(f), P(f->GetProlongationMatrix()),
-     grad(*this)
-{
-   G = fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
-   MFEM_ASSERT(G, "Some GetElementRestriction error");
-   x_local.SetSize(G->Height(), Device::GetMemoryType());
-   v_local.SetSize(G->Height(), Device::GetMemoryType());
-   y_local.SetSize(G->Height(), Device::GetMemoryType());
-}
-
-void ParVariationalForm::Mult(const Vector &x, Vector &y) const
-{
-   px.SetSize(P->Height());
-   py.SetSize(P->Height());
-
-   P->Mult(x, px);
-   G->Mult(px, x_local);
-
-   y_local = 0.0;
-   for (int i = 0; i < domain_integrators.Size(); ++i)
-   {
-      // y += F(x)
-      domain_integrators[i]->Apply(x_local, y_local);
-   }
-
-   G->MultTranspose(y_local, py);
-   P->MultTranspose(py, y);
-
-   y.HostReadWrite();
-   for (int i = 0; i < ess_tdof_list.Size(); i++)
-   {
-      y(ess_tdof_list[i]) = 0.0;
-   }
-}
-
-void ParVariationalForm::GradientMult(const Vector &v, Vector &y) const
-{
-   px.SetSize(P->Height());
-   py.SetSize(P->Height());
-   pv.SetSize(P->Height());
-
-   P->Mult(v, pv);
-   G->Mult(pv, v_local);
-
-   P->Mult(x_lin, px);
-   G->Mult(px, x_local);
-
-   y_local = 0.0;
-   for (int i = 0; i < domain_integrators.Size(); ++i)
-   {
-      // y += dF(x)/dx * v
-      if (is_linear)
-      {
-         // take care of RHS
-         // domain_integrators[i]->Apply(v_local, y_local);
-      }
-      domain_integrators[i]->ApplyGradient(x_local, v_local, y_local);
-   }
-
-   G->MultTranspose(y_local, py);
-   P->MultTranspose(py, y);
-
-   y.HostReadWrite();
-   for (int i = 0; i < ess_tdof_list.Size(); i++)
-   {
-      y(ess_tdof_list[i]) = v(ess_tdof_list[i]);
-   }
-}
-
-Operator &ParVariationalForm::GetGradient(const Vector &x) const
-{
-   x_lin = x;
-   return grad;
-}
-
-HypreParMatrix *ParVariationalForm::GetGradientMatrix(const Vector &x)
-{
-   delete gradient_matrix;
-   gradient_matrix = new HypreParMatrix;
-   return gradient_matrix;
-}
-
-void ParVariationalForm::SetEssentialBC(const Array<int> &ess_attr)
-{
-   fes->GetEssentialTrueDofs(ess_attr, ess_tdof_list);
-}
-
-} // namespace mfem

miniapps/variationalform/parvariationalform.hpp

@@ -1,109 +0,0 @@
-#include "mfem.hpp"
-#include "genericintegrator.hpp"
-#include "qfuncintegrator.hpp"
-
-#pragma once
-
-namespace mfem
-{
-class ParVariationalForm : public Operator
-{
-   class Gradient : public Operator
-   {
-   public:
-      Gradient(ParVariationalForm &f) : Operator(f.Height()), form(f){};
-
-      void Mult(const Vector &x, Vector &y) const override
-      {
-         form.GradientMult(x, y);
-      }
-
-   private:
-      ParVariationalForm &form;
-   };
-
-public:
-   ParVariationalForm(ParFiniteElementSpace *f);
-
-   void Mult(const Vector &x, Vector &y) const override;
-
-   void AddDomainIntegrator(GenericIntegrator *i)
-   {
-      domain_integrators.Append(i);
-      i->Setup(*fes);
-   }
-
-   template<typename integrator_type,
-            typename qfunc_type,
-            typename... qfunc_args_type>
-   void AddDomainIntegrator(qfunc_type f, qfunc_args_type const &... fargs)
-   {
-      if constexpr (std::is_same_v<integrator_type, DomainLFIntegrator>)
-      {
-         auto i = new QFunctionIntegrator(
-            [&](auto... args) {
-               auto du = std::get<1>(std::tuple{args...});
-               return std::tuple{f(args...), decltype(du){}};
-            },
-            0,
-            *fes->GetParMesh(),
-            fargs...);
-
-         domain_integrators.Append(i);
-         i->Setup(*fes);
-      }
-
-      if constexpr (std::is_same_v<integrator_type, DiffusionIntegrator>)
-      {
-         auto i = new QFunctionIntegrator(
-            [&](auto... args) {
-               auto du = std::get<1>(std::tuple{args...});
-               return std::tuple{0.0, f(args...) * du};
-            },
-            0,
-            *fes->GetParMesh(),
-            fargs...);
-
-         domain_integrators.Append(i);
-         i->Setup(*fes);
-      }
-   }
-
-   void AssumeLinear() { is_linear = true; }
-
-   // Return an Operator that provides a Mult(x, y) which is the MatVec of the
-   // gradient of the ParVariationalForm wrt x. Acts as a "passthrough" to
-   // ::GradientMult in order to satisfy mfem interfaces.
-   Operator &GetGradient(const Vector &x) const override;
-
-   // Return an assmbled parallel matrix which represents the gradient of the
-   // ParVariationalForm wrt to x.
-   HypreParMatrix *GetGradientMatrix(const Vector &x);
-
-   void SetEssentialBC(const Array<int> &ess_attr);
-
-protected:
-   // y = F'(x_lin) * v
-   void GradientMult(const Vector &v, Vector &y) const;
-
-   bool is_linear = false;
-   ParFiniteElementSpace *fes;
-   Array<GenericIntegrator *> domain_integrators;
-   Array<int> ess_tdof_list;
-
-   // T -> L
-   const Operator *P;
-
-   // L -> E
-   const Operator *G;
-
-   mutable Vector x_local, y_local, v_local, px, py, pv;
-
-   // State to build the Gradient on, single source of the true state.
-   mutable Vector x_lin;
-
-   mutable Gradient grad;
-
-   HypreParMatrix *gradient_matrix = nullptr;
-};
-} // namespace mfem

miniapps/variationalform/parvariationalform_ex1.cpp

@@ -1,164 +0,0 @@
-#include "mfem.hpp"
-#include "parvariationalform.hpp"
-#include "tensor.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   const char *mesh_file = "../data/inline-quad.mesh";
-   int order = 1;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = true;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
-   args.AddOption(&order,
-                  "-o",
-                  "--order",
-                  "Finite element order (polynomial degree) or -1 for"
-                  " isoparametric space.");
-   args.AddOption(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-
-   {
-      int ref_levels = (int) floor(log(10. / mesh.GetNE()) / log(2.) / dim);
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-   {
-      int par_ref_levels = 2;
-      for (int l = 0; l < par_ref_levels; l++)
-      {
-         pmesh.UniformRefinement();
-      }
-   }
-
-   FiniteElementCollection *fec;
-   bool delete_fec;
-   if (order > 0)
-   {
-      fec = new H1_FECollection(order, dim);
-      delete_fec = true;
-   }
-   else if (pmesh.GetNodes())
-   {
-      fec = pmesh.GetNodes()->OwnFEC();
-      delete_fec = false;
-      if (myid == 0)
-      {
-         cout << "Using isoparametric FEs: " << fec->Name() << endl;
-      }
-   }
-   else
-   {
-      fec = new H1_FECollection(order = 1, dim);
-      delete_fec = true;
-   }
-   ParFiniteElementSpace fespace(&pmesh, fec);
-   HYPRE_Int size = fespace.GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of finite element unknowns: " << size << endl;
-   }
-
-   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
-   ess_bdr = 1;
-
-   FunctionCoefficient u_excoeff([&](const Vector &coords) {
-      double x = coords(0);
-      double y = coords(1);
-      return x * x + y * y;
-   });
-
-   ParGridFunction x(&fespace);
-
-   x.ProjectBdrCoefficient(u_excoeff, ess_bdr);
-
-   ParVariationalForm form(&fespace);
-
-   auto b_coeff = [&](auto u, auto du, auto x) {
-      return 4.0 * (1.0 + 2.0 * x[0] * x[0] + 2.0 * x[1] * x[1]);
-   };
-
-   form.AddDomainIntegrator<DomainLFIntegrator>(b_coeff);
-
-   auto a_coeff = [&](auto u, auto du, auto x) { return 1.0 + u; };
-
-   form.AddDomainIntegrator<DiffusionIntegrator>(a_coeff);
-
-   form.SetEssentialBC(ess_bdr);
-
-   CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-6);
-   cg.SetMaxIter(2000);
-   cg.SetPrintLevel(1);
-
-   NewtonSolver newton(MPI_COMM_WORLD);
-   newton.SetOperator(form);
-   newton.SetSolver(cg);
-   newton.SetPrintLevel(1);
-   newton.SetRelTol(1e-8);
-   newton.SetMaxIter(100);
-
-   Vector zero;
-   Vector X;
-   x.GetTrueDofs(X);
-   newton.Mult(zero, X);
-
-   x.Distribute(X);
-
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int visport = 19916;
-      socketstream sol_sock(vishost, visport);
-      sol_sock << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock.precision(8);
-      sol_sock << "solution\n" << pmesh << x << flush;
-   }
-
-   if (delete_fec)
-   {
-      delete fec;
-   }
-   MPI_Finalize();
-
-   return 0;
-}

miniapps/variationalform/parvariationalform_example.cpp

@@ -1,118 +0,0 @@
-#include "mfem.hpp"
-#include "parvariationalform.hpp"
-#include "qfuncintegrator.hpp"
-#include "tensor.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   const char *mesh_file = "../data/inline-quad.mesh";
-   int order = 1;
-   int refinements = 0;
-   double p = 5.0;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&refinements, "-r", "--ref", "");
-   args.AddOption(&order, "-o", "--order", "");
-
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   Mesh mesh(mesh_file, 1, 1);
-
-   int dim = mesh.Dimension();
-   {
-      for (int l = 0; l < refinements; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-
-   auto fec = H1_FECollection(order, dim);
-   ParFiniteElementSpace fespace(&pmesh, &fec);
-
-   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
-   ess_bdr = 1;
-
-   FunctionCoefficient u_excoeff([&](const Vector &coords) {
-      double x = coords(0);
-      double y = coords(1);
-
-      // return 1.0 - pow(sqrt(x * x + y * y), p / (p - 1.0));
-      return x * x + y * y;
-   });
-
-   ParGridFunction x(&fespace);
-   x.Randomize();
-
-   x.ProjectBdrCoefficient(u_excoeff, ess_bdr);
-
-   ParVariationalForm form(&fespace);
-
-   auto plaplacian = new QFunctionIntegrator([&](auto u, auto du, auto x) {
-      // auto f0 = -1.0;
-      // auto f1 = pow(norm(du), p - 2.0) * du;
-      auto f0 = 4.0 * (1.0 + 2.0 * x[0] * x[0] + 2.0 * x[1] * x[1]);
-      auto f1 = (1.0 + u) * du;
-      return std::tuple{f0, f1};
-   }, 0, pmesh);
-
-   form.AddDomainIntegrator(plaplacian);
-
-   form.SetEssentialBC(ess_bdr);
-
-   CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-6);
-   cg.SetMaxIter(2000);
-   cg.SetPrintLevel(1);
-
-   NewtonSolver newton(MPI_COMM_WORLD);
-   newton.SetOperator(form);
-   newton.SetSolver(cg);
-   newton.SetPrintLevel(1);
-   newton.SetRelTol(1e-8);
-   newton.SetMaxIter(100);
-
-   Vector zero;
-   Vector X;
-   x.GetTrueDofs(X);
-   newton.Mult(zero, X);
-
-   x.Distribute(X);
-
-   // x.ProjectCoefficient(u_excoeff);
-   
-   char vishost[] = "localhost";
-   int visport = 19916;
-   socketstream sol_sock(vishost, visport);
-   sol_sock << "parallel " << num_procs << " " << myid << "\n";
-   sol_sock.precision(8);
-   sol_sock << "solution\n" << pmesh << x << flush;
-
-   MPI_Finalize();
-
-   return 0;
-}

miniapps/variationalform/qfuncintegrator.hpp

@@ -1,404 +0,0 @@
-#include "mfem.hpp"
-#include "../../general/forall.hpp"
-#include "genericintegrator.hpp"
-#include "tensor.hpp"
-
-#pragma once
-
-namespace mfem
-{
-template<typename T>
-struct supported_type
-{
-   static constexpr bool value = false;
-};
-
-template<>
-struct supported_type<ParMesh>
-{
-   static constexpr bool value = true;
-};
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-class QFunctionIntegrator : public GenericIntegrator
-{
-public:
-   QFunctionIntegrator(qfunc_type f,
-                       qfunc_grad_type f_grad,
-                       qfunc_args_type const &... fargs);
-
-   QFunctionIntegrator(qfunc_type f, qfunc_args_type const &... fargs);
-
-   void Setup(const FiniteElementSpace &fes) override;
-
-   void Apply(const Vector &, Vector &) const override;
-
-   // y += F'(x) * v
-   void ApplyGradient(const Vector &x,
-                      const Vector &v,
-                      Vector &y) const override;
-
-protected:
-   template<int D1D, int Q1D>
-   void Apply2D(const Vector &u_in_, Vector &y_) const;
-
-   template<int D1D, int Q1D>
-   void ApplyGradient2D(const Vector &u_in_,
-                        const Vector &v_in_,
-                        Vector &y_) const;
-
-   auto EvaluateFargValue(const Mesh &m, const double qx, const double qy) const;
-
-   const FiniteElementSpace *fespace;
-   const DofToQuad *maps;        ///< Not owned
-   const GeometricFactors *geom; ///< Not owned
-   int dim, ne, nq, dofs1D, quad1D;
-
-   // Geometric factors
-   Vector J_;
-   Vector W_;
-
-   qfunc_type qf;
-   qfunc_grad_type qf_grad;
-   std::tuple<qfunc_args_type...> qf_farg_values;
-};
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
-   QFunctionIntegrator(qfunc_type f,
-                       qfunc_grad_type df,
-                       qfunc_args_type const &... fargs)
-   : GenericIntegrator(nullptr), maps(nullptr), geom(nullptr), qf(f),
-     qf_grad(df), qf_farg_values(std::tuple{fargs...})
-{
-   static_assert((supported_type<qfunc_args_type>::value && ...),
-                 "Type not supported for parameter expansion. See "
-                 "documentation for supported types.");
-}
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
-   QFunctionIntegrator(qfunc_type f, qfunc_args_type const &... fargs)
-   : GenericIntegrator(nullptr), maps(nullptr), geom(nullptr), qf(f),
-     qf_grad(qfunc_grad_type{}), qf_farg_values(std::tuple{fargs...})
-{
-   static_assert((supported_type<qfunc_args_type>::value && ...),
-                 "Type not supported for parameter expansion. See "
-                 "documentation for supported types.");
-}
-
-template<typename qfunc_type, typename... qfunc_args_type>
-QFunctionIntegrator(qfunc_type, qfunc_args_type const &...)
-   -> QFunctionIntegrator<qfunc_type, int, qfunc_args_type const &...>;
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::Setup(
-   const FiniteElementSpace &fes)
-{
-   // Assuming the same element type
-   fespace = &fes;
-   Mesh *mesh = fes.GetMesh();
-   if (mesh->GetNE() == 0)
-   {
-      return;
-   }
-   const FiniteElement &el = *fes.GetFE(0);
-   ElementTransformation *T = mesh->GetElementTransformation(0);
-   const IntegrationRule *ir = nullptr;
-   if (!IntRule)
-   {
-      IntRule = &IntRules.Get(el.GetGeomType(), el.GetOrder() * 2);
-   }
-   ir = IntRule;
-
-   dim = mesh->Dimension();
-   ne = fes.GetMesh()->GetNE();
-   nq = ir->GetNPoints();
-   geom = mesh->GetGeometricFactors(*ir,
-                                    GeometricFactors::COORDINATES
-                                       | GeometricFactors::JACOBIANS);
-   maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
-   dofs1D = maps->ndof;
-   quad1D = maps->nqpt;
-   //    pa_data.SetSize(ne * nq, Device::GetDeviceMemoryType());
-
-   W_.SetSize(nq, Device::GetDeviceMemoryType());
-   W_.GetMemory().CopyFrom(ir->GetWeights().GetMemory(), nq);
-
-   // J.SetSize(ne * nq, Device::GetDeviceMemoryType());
-   J_ = geom->J;
-}
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-auto QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
-   EvaluateFargValue(const Mesh &m, const double qx, const double qy) const
-{
-   Vector trip(3);
-   trip = 0.0;
-   ElementTransformation *tr = const_cast<Mesh &>(m).GetElementTransformation(0);
-   tr->Transform(IntRule->IntPoint(qx + quad1D * qy), trip);
-   return tensor<double, 3>{
-      {trip(0), trip(1), m.SpaceDimension() == 2 ? 0.0 : trip(2)}};
-}
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::Apply(
-   const Vector &x, Vector &y) const
-{
-   if (dim == 2)
-   {
-      switch ((dofs1D << 4) | quad1D)
-      {
-      case 0x22:
-         return Apply2D<2, 2>(x, y);
-      default:
-         MFEM_ASSERT(false, "NOPE");
-      }
-   }
-}
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-template<int D1D, int Q1D>
-void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::Apply2D(
-   const Vector &u_in_, Vector &y_) const
-{
-   int NE = ne;
-
-   auto v1d = Reshape(maps->B.Read(), Q1D, D1D);
-   auto dv1d_dX = Reshape(maps->G.Read(), Q1D, D1D);
-   // (NQ x SDIM x DIM x NE)
-   auto J = Reshape(J_.Read(), Q1D, Q1D, 2, 2, NE);
-   auto W = Reshape(W_.Read(), Q1D, Q1D);
-   auto u = Reshape(u_in_.Read(), D1D, D1D, NE);
-   auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
-
-   // MFEM_FORALL(e, NE, {
-   for (int e = 0; e < NE; e++)
-   {
-      // loop over quadrature points
-      for (int qy = 0; qy < Q1D; ++qy)
-      {
-         for (int qx = 0; qx < Q1D; ++qx)
-         {
-            double u_q = 0.0;
-            double du_dX_q[2] = {0.0};
-            for (int ix = 0; ix < D1D; ix++)
-            {
-               for (int iy = 0; iy < D1D; iy++)
-               {
-                  u_q += u(ix, iy, e) * v1d(qx, ix) * v1d(qy, iy);
-                  du_dX_q[0] += u(ix, iy, e) * dv1d_dX(qx, ix) * v1d(qy, iy);
-                  du_dX_q[1] += u(ix, iy, e) * v1d(qx, ix) * dv1d_dX(qy, iy);
-               }
-            }
-
-            // du_dx_q = invJ^T * du_dX_q
-            //         = (adjJ^T * du_dX_q) / detJ
-            double J_q[2][2] = {{J(qx, qy, 0, 0, e),
-                                 J(qx, qy, 0, 1, e)}, // J_q[0][0], J_q[0][1]
-                                {J(qx, qy, 1, 0, e),
-                                 J(qx, qy, 1, 1, e)}}; // J_q[1][0], J_q[1][1]
-
-            double detJ_q = (J_q[0][0] * J_q[1][1]) - (J_q[0][1] * J_q[1][0]);
-
-            double adjJ[2][2] = {{J_q[1][1], -J_q[0][1]},
-                                 {-J_q[1][0], J_q[0][0]}};
-
-            tensor<double, 2> du_dx_q
-               = {(adjJ[0][0] * du_dX_q[0] + adjJ[1][0] * du_dX_q[1]) / detJ_q,
-                  (adjJ[0][1] * du_dX_q[0] + adjJ[1][1] * du_dX_q[1]) / detJ_q};
-
-            auto processed_qf_farg_values = std::apply(
-               [=](auto &... a) {
-                  return std::make_tuple(u_q,
-                                         du_dx_q,
-                                         EvaluateFargValue(a, qx, qy)...);
-               },
-               qf_farg_values);
-
-            auto [f0, f1] = std::apply(qf, processed_qf_farg_values);
-
-            double f0_X = f0 * detJ_q;
-
-            // f1_X = invJ * f1 * detJ
-            //      = adjJ * f1
-            double f1_X[2] = {
-               adjJ[0][0] * f1[0] + adjJ[0][1] * f1[1],
-               adjJ[1][0] * f1[0] + adjJ[1][1] * f1[1],
-            };
-
-            for (int ix = 0; ix < D1D; ix++)
-            {
-               for (int iy = 0; iy < D1D; iy++)
-               {
-                  // accumulate v * f0 + dot(dv_dx, f1)
-                  y(ix, iy, e) += (f0_X * v1d(qx, ix) * v1d(qy, iy)
-                                   + f1_X[0] * dv1d_dX(qx, ix) * v1d(qy, iy)
-                                   + f1_X[1] * dv1d_dX(qy, iy) * v1d(qx, ix))
-                                  * W(qx, qy);
-               }
-            }
-         }
-      }
-   }
-   // });
-}
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
-   ApplyGradient(const Vector &x, const Vector &v, Vector &y) const
-{
-   ApplyGradient2D<2, 2>(x, v, y);
-}
-
-template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
-template<int D1D, int Q1D>
-void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
-   ApplyGradient2D(const Vector &u_in_, const Vector &v_in_, Vector &y_) const
-{
-   int NE = ne;
-
-   auto v1d = Reshape(maps->B.Read(), Q1D, D1D);
-   auto dv1d_dX = Reshape(maps->G.Read(), Q1D, D1D);
-   // (NQ x SDIM x DIM x NE)
-   auto J = Reshape(J_.Read(), Q1D, Q1D, 2, 2, NE);
-   auto W = Reshape(W_.Read(), Q1D, Q1D);
-   auto u = Reshape(u_in_.Read(), D1D, D1D, NE);
-   auto v = Reshape(v_in_.Read(), D1D, D1D, NE);
-   auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
-
-   for (int e = 0; e < NE; e++)
-   {
-      // loop over quadrature points
-      for (int qy = 0; qy < Q1D; ++qy)
-      {
-         for (int qx = 0; qx < Q1D; ++qx)
-         {
-            double u_q = 0.0;
-            double du_dX_q[2] = {0.0};
-            double v_q = 0.0;
-            double dv_dX_q[2] = {0.0};
-            for (int ix = 0; ix < D1D; ix++)
-            {
-               for (int iy = 0; iy < D1D; iy++)
-               {
-                  u_q += u(ix, iy, e) * v1d(qx, ix) * v1d(qy, iy);
-
-                  du_dX_q[0] += u(ix, iy, e) * dv1d_dX(qx, ix) * v1d(qy, iy);
-                  du_dX_q[1] += u(ix, iy, e) * v1d(qx, ix) * dv1d_dX(qy, iy);
-
-                  v_q += v(ix, iy, e) * v1d(qx, ix) * v1d(qy, iy);
-
-                  dv_dX_q[0] += v(ix, iy, e) * dv1d_dX(qx, ix) * v1d(qy, iy);
-                  dv_dX_q[1] += v(ix, iy, e) * v1d(qx, ix) * dv1d_dX(qy, iy);
-               }
-            }
-
-            // du_dx_q = invJ^T * du_dX_q
-            //         = (adjJ^T * du_dX_q) / detJ
-            double J_q[2][2] = {{J(qx, qy, 0, 0, e),
-                                 J(qx, qy, 0, 1, e)}, // J_q[0][0], J_q[0][1]
-                                {J(qx, qy, 1, 0, e),
-                                 J(qx, qy, 1, 1, e)}}; // J_q[1][0], J_q[1][1]
-
-            double detJ_q = (J_q[0][0] * J_q[1][1]) - (J_q[0][1] * J_q[1][0]);
-
-            double adjJ[2][2] = {{J_q[1][1], -J_q[0][1]},
-                                 {-J_q[1][0], J_q[0][0]}};
-
-            tensor<double, 2> du_dx_q
-               = {(adjJ[0][0] * du_dX_q[0] + adjJ[1][0] * du_dX_q[1]) / detJ_q,
-                  (adjJ[0][1] * du_dX_q[0] + adjJ[1][1] * du_dX_q[1]) / detJ_q};
-
-            double dv_dx_q[2]
-               = {(adjJ[0][0] * dv_dX_q[0] + adjJ[1][0] * dv_dX_q[1]) / detJ_q,
-                  (adjJ[0][1] * dv_dX_q[0] + adjJ[1][1] * dv_dX_q[1]) / detJ_q};
-
-            // compute dF(u, du)/du
-            auto processed_qf_farg_values_u = std::apply(
-               [=](auto &... a) {
-                  return std::make_tuple(derivative_wrt(u_q),
-                                         du_dx_q,
-                                         EvaluateFargValue(a, qx, qy)...);
-               },
-               qf_farg_values);
-
-            auto [f0u, f1u] = std::apply(qf, processed_qf_farg_values_u);
-
-            double f00 = 0.0;
-            if constexpr (std::is_same_v<decltype(f0u), double>)
-            {
-               f00 = 0.0;
-            }
-            else
-            {
-               f00 = f0u.gradient;
-            }
-
-            tensor<double, 2> f10;
-            if constexpr (std::is_same_v<decltype(f1u), tensor<double, 2>>)
-            {
-               f10 = {0.0, 0.0};
-            }
-            else
-            {
-               f10 = {f1u[0].gradient, f1u[1].gradient};
-            }
-
-            // compute dF(u, du)/ddu
-            auto processed_qf_farg_values_du = std::apply(
-               [=](auto &... a) {
-                  return std::make_tuple(u_q,
-                                         derivative_wrt(du_dx_q),
-                                         EvaluateFargValue(a, qx, qy)...);
-               },
-               qf_farg_values);
-
-            auto [f0du, f1du] = std::apply(qf, processed_qf_farg_values_du);
-
-            tensor<double, 2> f01;
-            if constexpr (std::is_same_v<decltype(f0du), double>)
-            {
-               f01 = {0.0, 0.0};
-            }
-            else
-            {
-               f01 = f0du.gradient;
-            }
-            tensor<double, 2, 2> f11 = {
-               {{f1du[0].gradient[0], f1du[0].gradient[1]},
-                {f1du[1].gradient[0], f1du[1].gradient[1]}}};
-
-            double W0 = f00 * v_q + f01[0] * dv_dx_q[0] + f01[1] * dv_dx_q[1];
-
-            double W1[2] = {f10[0] * v_q + f11[0][0] * dv_dx_q[0]
-                               + f11[0][1] * dv_dx_q[1],
-                            f10[1] * v_q + f11[1][0] * dv_dx_q[0]
-                               + +f11[1][1] * dv_dx_q[1]};
-
-            double W0_X = W0 * detJ_q;
-
-            // W1_X = invJ * W1 * detJ
-            //      = adjJ * W1
-            double W1_X[2] = {
-               adjJ[0][0] * W1[0] + adjJ[0][1] * W1[1],
-               adjJ[1][0] * W1[0] + adjJ[1][1] * W1[1],
-            };
-
-            for (int ix = 0; ix < D1D; ix++)
-            {
-               for (int iy = 0; iy < D1D; iy++)
-               {
-                  // @TODO: proper comment
-                  y(ix, iy, e) += (W0_X * v1d(qx, ix) * v1d(qy, iy)
-                                   + W1_X[0] * dv1d_dX(qx, ix) * v1d(qy, iy)
-                                   + W1_X[1] * dv1d_dX(qy, iy) * v1d(qx, ix))
-                                  * W(qx, qy);
-               }
-            }
-         }
-      }
-   }
-}
-
-} // namespace mfem

miniapps/variationalform/tensor.hpp

@@ -1,347 +0,0 @@
-// tensor
-#include <iostream>
-
-#pragma once
-
-template<typename T, int... n>
-struct tensor;
-
-template<typename T>
-struct tensor<T, 1>
-{
-   static constexpr int shape[1] = {1};
-   operator T() { return value; }
-   T value;
-};
-
-template<typename T, int n>
-struct tensor<T, n>
-{
-   static constexpr int shape[1] = {n};
-   constexpr auto &operator[](int i) { return value[i]; };
-   constexpr auto operator[](int i) const { return value[i]; };
-   T value[n];
-};
-
-template<typename T, int first, int... rest>
-struct tensor<T, first, rest...>
-{
-   static constexpr int shape[1 + sizeof...(rest)] = {first, rest...};
-   constexpr auto &operator[](int i) { return value[i]; };
-   constexpr auto operator[](int i) const { return value[i]; };
-   tensor<T, rest...> value[first];
-};
-
-template<int n>
-constexpr int product(int (&values)[n])
-{
-   int p = 1;
-   for (int i = 0; i < n; i++)
-   {
-      p *= values[i];
-   }
-   return p;
-}
-
-template<typename S, typename T, int... n>
-auto operator+(tensor<S, n...> A, tensor<T, n...> B)
-{
-   tensor<decltype(S{} + T{}), n...> C{};
-   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
-   {
-      C[i] = A[i] + B[i];
-   }
-   return C;
-}
-
-template<typename S, typename T, int... n>
-auto operator-(tensor<S, n...> A, tensor<T, n...> B)
-{
-   tensor<decltype(S{} - T{}), n...> C{};
-   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
-   {
-      C[i] = A[i] - B[i];
-   }
-   return C;
-}
-
-template<typename S, typename T, int... n>
-auto operator*(S scale, tensor<T, n...> A)
-{
-   tensor<decltype(S{} * T{}), n...> C{};
-   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
-   {
-      C[i] = scale * A[i];
-   }
-   return C;
-}
-
-template<typename S, typename T, int... n>
-auto operator*(tensor<T, n...> A, S scale)
-{
-   tensor<decltype(T{} * S{}), n...> C{};
-   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
-   {
-      C[i] = A[i] * scale;
-   }
-   return C;
-}
-
-template<typename S, typename T, int... n>
-auto operator/(S scale, tensor<T, n...> A)
-{
-   tensor<decltype(S{} / T{}), n...> C{};
-   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
-   {
-      C[i] = scale / A[i];
-   }
-   return C;
-}
-
-template<typename S, typename T, int... n>
-auto operator/(tensor<T, n...> A, S scale)
-{
-   tensor<decltype(T{} / S{}), n...> C{};
-   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
-   {
-      C[i] = A[i] / scale;
-   }
-   return C;
-}
-
-template<typename S, typename T, int m, int n, int p>
-auto dot(tensor<S, m, n> A, tensor<T, n, p> B)
-{
-   tensor<decltype(S{} * T{}), m, p> AB{};
-   for (int i = 0; i < m; i++)
-   {
-      for (int j = 0; j < p; j++)
-      {
-         for (int k = 0; k < n; k++)
-         {
-            AB[i][j] = AB[i][j] + A[i][k] * B[k][j];
-         }
-      }
-   }
-   return AB;
-}
-
-template<typename T, int m, int n>
-auto inner(tensor<T, m, n> A, tensor<T, m, n> B)
-{
-   double value = 0.0;
-   for (int i = 0; i < m; i++)
-   {
-      for (int j = 0; j < n; j++)
-      {
-         value += A[i][j] * B[i][j];
-      }
-   }
-   return value;
-}
-
-auto inner(double a, double b)
-{
-   return a * b;
-}
-
-template<typename T, int n>
-auto tr(tensor<T, n, n> A)
-{
-   T trA{};
-   for (int i = 0; i < n; i++)
-   {
-      trA = trA + A[i][i];
-   }
-   return trA;
-}
-
-template<int dim>
-constexpr tensor<double, dim, dim> Identity()
-{
-   tensor<double, dim, dim> I{};
-   for (int i = 0; i < dim; i++)
-   {
-      for (int j = 0; j < dim; j++)
-      {
-         I[i][j] = (i == j);
-      }
-   }
-   return I;
-}
-
-template<typename T, int m, int n>
-auto transpose(const tensor<T, m, n> &A)
-{
-   tensor<T, n, m> AT{};
-   for (int i = 0; i < n; i++)
-   {
-      for (int j = 0; j < m; j++)
-      {
-         AT[i][j] = A[j][i];
-      }
-   }
-   return AT;
-}
-
-template<typename T, int n>
-auto norm(const tensor<T, n> &A)
-{
-   T r = {};
-   for (int i = 0; i < n; i++)
-   {
-      r = r + A[i] * A[i];
-   }
-   return pow(r, 0.5);
-}
-
-template<typename gradient_type>
-struct dual
-{
-   double value;
-   gradient_type gradient;
-};
-
-template<typename gradient_type>
-auto operator+(dual<gradient_type> a, double b)
-{
-   return dual<gradient_type>{a.value + b, a.gradient};
-}
-
-template<typename gradient_type>
-auto operator+(double a, dual<gradient_type> b)
-{
-   return dual<gradient_type>{a + b.value, b.gradient};
-}
-
-template<typename gradient_type>
-auto operator+(dual<gradient_type> a, dual<gradient_type> b)
-{
-   return dual<gradient_type>{a.value + b.value, a.gradient + b.gradient};
-}
-
-template<typename gradient_type>
-auto operator*(dual<gradient_type> a, double b)
-{
-   return dual<gradient_type>{a.value * b, a.gradient * b};
-}
-
-template<typename gradient_type>
-auto operator*(double a, dual<gradient_type> b)
-{
-   return dual<gradient_type>{a * b.value, a * b.gradient};
-}
-
-template<typename gradient_type>
-auto operator*(dual<gradient_type> a, dual<gradient_type> b)
-{
-   return dual<gradient_type>{a.value * b.value,
-                              b.value * a.gradient + a.value * b.gradient};
-}
-
-template<typename gradient_type>
-auto cos(dual<gradient_type> a)
-{
-   return dual<gradient_type>{cos(a.value), -a.gradient * sin(a.value)};
-}
-
-template<typename gradient_type>
-auto exp(dual<gradient_type> a)
-{
-   return dual<gradient_type>{exp(a.value), exp(a.value)};
-}
-
-template<typename gradient_type>
-auto log(dual<gradient_type> a)
-{
-   return dual<gradient_type>{log(a.value), a.gradient / a.value};
-}
-
-template<typename gradient_type>
-auto pow(dual<gradient_type> a, dual<gradient_type> b)
-{
-   double value = pow(a.value, b.value);
-   return dual<gradient_type>{value,
-                              value
-                                 * (a.gradient * (b.value / a.value)
-                                    + b.gradient * log(a.value))};
-}
-
-template<typename gradient_type>
-auto pow(double a, dual<gradient_type> b)
-{
-   double value = pow(a, b.value);
-   return dual<gradient_type>{value, value * b.gradient * log(a)};
-}
-
-template<typename gradient_type>
-auto pow(dual<gradient_type> a, double b)
-{
-   double value = pow(a.value, b);
-   return dual<gradient_type>{value, value * a.gradient * b / a.value};
-}
-
-template<typename T, int... n>
-auto &operator<<(std::ostream &out, dual<T> A)
-{
-   out << '(' << A.value << ' ' << A.gradient << ')';
-   return out;
-}
-
-template<typename T, int... n>
-auto &operator<<(std::ostream &out, tensor<T, n...> A)
-{
-   out << '{' << A[0];
-   for (int i = 1; i < tensor<T, n...>::shape[0]; i++)
-   {
-      out << ", " << A[i];
-   }
-   out << '}';
-   return out;
-}
-
-auto derivative_wrt(double a)
-{
-   return dual<double>{a, 1};
-}
-
-template<typename T, int m>
-auto derivative_wrt(tensor<T, m> A)
-{
-   tensor<dual<tensor<double, m>>, m> A_dual{};
-   for (int i = 0; i < m; i++)
-   {
-      A_dual[i].value = A[i];
-      A_dual[i].gradient[i] = 1.0;
-   }
-   return A_dual;
-}
-
-template<typename T, int m, int n>
-auto derivative_wrt(tensor<T, m, n> A)
-{
-   tensor<dual<tensor<double, m, n>>, m, n> A_dual{};
-   for (int i = 0; i < m; i++)
-   {
-      for (int j = 0; j < n; j++)
-      {
-         A_dual[i][j].value = A[i][j];
-         A_dual[i][j].gradient[i][j] = 1.0;
-      }
-   }
-   return A_dual;
-}
-
-template<typename grad_type, int nrows, int ncols>
-auto directional_derivative(tensor<dual<grad_type>, nrows, ncols> A, grad_type n)
-{
-   tensor<double, nrows, ncols> dA_dn{};
-   for (int i = 0; i < nrows; i++)
-   {
-      for (int j = 0; j < ncols; j++)
-      {
-         dA_dn[i][j] = inner(A[i][j].gradient, n);
-      }
-   }
-   return dA_dn;
-}

miniapps/variationalform/test_tensor_ad.cpp

@@ -1,69 +0,0 @@
-#include "mfem.hpp"
-#include "tensor.hpp"
-
-#define CATCH_CONFIG_MAIN
-#include "catch.hpp"
-
-using namespace mfem;
-
-TEST_CASE("Tensor tests", "[tensor]")
-{
-   SECTION("norm")
-   {
-      tensor<double, 5> a = {{1.0, 2.0, 3.0, 4.0, 5.0}};
-      REQUIRE(norm(a) - sqrt(55) == Approx(0.0));
-   }
-}
-
-TEST_CASE("Dual number tensor tests", "[DualNumber]")
-{
-   auto eps = std::numeric_limits<double>::epsilon();
-   double x = 0.5;
-
-   SECTION("cos")
-   {
-      auto xd = cos(derivative_wrt(x));
-      REQUIRE(abs(-sin(x) - xd.gradient) == Approx(0.0));
-   }
-
-   SECTION("exp")
-   {
-      auto xd = exp(derivative_wrt(x));
-      REQUIRE(abs(exp(x) - xd.gradient) == Approx(0.0));
-   }
-
-   SECTION("log")
-   {
-      auto xd = log(derivative_wrt(x));
-      REQUIRE(abs(1.0 / x - xd.gradient) == Approx(0.0));
-   }
-
-   SECTION("pow")
-   {
-      // f(x) = x^3/2
-      auto xd = pow(derivative_wrt(x), 1.5);
-      REQUIRE(abs(1.5 * pow(x, 0.5) - xd.gradient) == Approx(0.0));
-   }
-
-   SECTION("mixed operations")
-   {
-      auto xd = derivative_wrt(x);
-      auto r = cos(xd) * cos(xd);
-      REQUIRE(abs(-2.0 * sin(x) * cos(x) - r.gradient) == Approx(0.0));
-
-      r = exp(xd) * cos(xd);
-      REQUIRE(abs(exp(x) * (cos(x) - sin(x)) - r.gradient) < eps);
-
-      r = log(xd) * cos(xd);
-      REQUIRE(abs((cos(x) / x - log(x) * sin(x)) - r.gradient) < eps);
-
-      r = exp(xd) * pow(xd, 1.5);
-      REQUIRE(abs((exp(x) * (pow(x, 1.5) + 1.5 * pow(x, 0.5))) - r.gradient)
-              < eps);
-
-      tensor<double, 2> vx = {{0.5, 0.25}};
-      tensor<double, 2> vre = {{0.894427190999916, 0.4472135954999579}};
-      auto vr = norm(derivative_wrt(vx));
-      REQUIRE(norm(vr.gradient - vre) < eps);
-   }
-}

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),