Merge branch 'master' into mprotect_trace

Fixed NULL pointer segfault in GetSurfaceFittingErrors

CHANGELOG

@@ -17,6 +17,9 @@ Discretization improvements
   Vector and VectorFE, also NURBS versions. Optionally different types of
   projections can be selected, default behaviour has not changed.
 
+- Added methods to estimate function extremum using piecewise linear bounds +
+  recursive subdivision.
+
 Meshing improvements
 --------------------
 - Improved support for 1D NURBS meshes with variable order, including using

config/cmake/modules/FindLIBBACKTRACE.cmake

@@ -0,0 +1,55 @@
+#[=======================================================================[.rst:
+FindLIBBACKTRACE
+-------
+
+Finds the LIBBACKTRACE library.
+
+Result Variables
+^^^^^^^^^^^^^^^^
+
+This will define the following variables:
+
+``LIBBACKTRACE_FOUND``
+  True if the system has the LIBBACKTRACE library.
+``LIBBACKTRACE_INCLUDE_DIRS``
+  Include directories needed to use LIBBACKTRACE.
+``LIBBACKTRACE_LIBRARIES``
+  Libraries needed to link to LIBBACKTRACE.
+#]=======================================================================]
+
+find_package(LIBBACKTRACE QUIET NO_MODULE)
+if(LIBBACKTRACE_FOUND)
+    message("LIBBACKTRACE ${LIBBACKTRACE_FIND_VERSION} found.")
+    if(LIBBACKTRACE_FIND_COMPONENTS)
+        message("LIBBACKTRACE components found:")
+        message("${LIBBACKTRACE_FIND_COMPONENTS}")
+    endif()
+    return()
+endif()
+
+set(LIBBACKTRACE_FOUND TRUE)
+
+## Find headers and libraries
+# should be in CMAKE_PREFIX_PATH if libunwind is
+# loaded with spack.
+find_library(LIBBACKTRACE_LIBRARIES libbacktrace.so)
+find_path(LIBBACKTRACE_INCLUDE_DIRS backtrace.h)
+if(NOT LIBBACKTRACE_LIBRARIES OR NOT LIBBACKTRACE_INCLUDE_DIRS)
+    set(LIBBACKTRACE_FOUND FALSE)
+endif()
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(LIBBACKTRACE DEFAULT_MSG
+                                  LIBBACKTRACE_FOUND
+                                  LIBBACKTRACE_LIBRARIES
+                                  LIBBACKTRACE_INCLUDE_DIRS)
+
+if (LIBBACKTRACE_FOUND)
+   add_library(LIBBACKTRACE::LIBBACKTRACE UNKNOWN IMPORTED)
+   set_target_properties(LIBBACKTRACE::LIBBACKTRACE PROPERTIES
+       INTERFACE_INCLUDE_DIRECTORIES "${LIBBACKTRACE_INCLUDE_DIRS}"
+       IMPORTED_LOCATION ${LIBBACKTRACE_LIBRARIES}
+       )
+endif()
+
+mark_as_advanced(LIBBACKTRACE_LIBRARIES LIBBACKTRACE_INCLUDE_DIRS LIBBACKTRACE_FOUND)

config/defaults.cmake

@@ -71,6 +71,7 @@ option(MFEM_USE_BENCHMARK "Enable Google Benchmark" OFF)
 option(MFEM_USE_PARELAG "Enable ParELAG" OFF)
 option(MFEM_USE_TRIBOL "Enable Tribol" OFF)
 option(MFEM_USE_ENZYME "Enable Enzyme" OFF)
+option(MFEM_USE_LIBBACKTRACE "Enable Libbacktrace" OFF)
 
 # Optional overrides for autodetected MPIEXEC and MPIEXEC_NUMPROC_FLAG
 # set(MFEM_MPIEXEC "mpirun" CACHE STRING "Command for running MPI tests")

fem/bounds.cpp

@@ -39,8 +39,8 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
    b_type = b_type_i;
    cp_type = cp_type_i;
    tol = tol_i;
-   lbound.SetSize(nb, ncp);
-   ubound.SetSize(nb, ncp);
+   lbound.SetSize(ncp, nb);
+   ubound.SetSize(ncp, nb);
    nodes.SetSize(nb);
    weights.SetSize(nb);
    control_points.SetSize(ncp);
@@ -125,21 +125,25 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
       {
          if (j == 0)
          {
-            lbound(i, j) = bv(i);
-            ubound(i, j) = bv(i);
+            lbound(j,i) = bv(i);
+            ubound(j,i) = bv(i);
          }
          else if (j == ncp-1)
          {
-            lbound(i, j) = bv(i);
-            ubound(i, j) = bv(i);
+            lbound(j,i) = bv(i);
+            ubound(j,i) = bv(i);
          }
          else
          {
             vals(0)  = bv(i);
             vals(1) =  bmv(i) +  dm*bdmv(i);
             vals(2) =  bpv(i) +  dp*bdpv(i);
-            lbound(i, j) = vals.Min()-tol; // tolerance for good measure
-            ubound(i, j) = vals.Max()+tol; // tolerance for good measure
+            lbound(j,i) = vals.Min()-tol; // tolerance for good measure
+            ubound(j,i) = vals.Max()+tol; // tolerance for good measure
+            if (b_type == 2)
+            {
+               lbound(j,i) = std::max(lbound(j,i),0_r);
+            }
          }
       }
    }
@@ -273,8 +277,7 @@ void PLBound::Get1DBounds(const Vector &coeff, Vector &intmin,
    intmax.SetSize(ncp);
    intmin = 0.0;
    intmax = 0.0;
-   Vector coeffm(nb);
-   coeffm = 0.0;
+   Vector coeffm;
 
    real_t a0 = 0.0;
    real_t a1 = 0.0;
@@ -302,6 +305,8 @@ void PLBound::Get1DBounds(const Vector &coeff, Vector &intmin,
    // compute L2 projection for linear bases: a0 + a1*x
    if (proj)
    {
+      coeffm.SetSize(nb);
+      coeffm = 0.0;
       for (int i = 0; i < nb; i++)
       {
          x = 2.0*nodes_int(i)-1;
@@ -342,8 +347,8 @@ void PLBound::Get1DBounds(const Vector &coeff, Vector &intmin,
       real_t c = coeffm(i);
       for (int j = 0; j < ncp; j++)
       {
-         intmin(j) += min(lbound(i,j)*c, ubound(i,j)*c);
-         intmax(j) += max(lbound(i,j)*c, ubound(i,j)*c);
+         intmin(j) += min(lbound(j,i)*c, ubound(j,i)*c);
+         intmax(j) += max(lbound(j,i)*c, ubound(j,i)*c);
       }
    }
 }
@@ -474,10 +479,10 @@ void PLBound::Get2DBounds(const Vector &coeff, Vector &intmin,
          real_t w1 = intmaxT(id2++);
          for (int k = 0; k < ncp; k++) // kth row
          {
-            vals(0) = w0*lbound(j,k);
-            vals(1) = w0*ubound(j,k);
-            vals(2) = w1*lbound(j,k);
-            vals(3) = w1*ubound(j,k);
+            vals(0) = w0*lbound(k,j);
+            vals(1) = w0*ubound(k,j);
+            vals(2) = w1*lbound(k,j);
+            vals(3) = w1*ubound(k,j);
             intmin(k*ncp+i) += vals.Min();
             intmax(k*ncp+i) += vals.Max();
          }
@@ -553,17 +558,17 @@ void PLBound::Get3DBounds(const Vector &coeff, Vector &intmin,
             for (int i = 0; i < nb; i++)
             {
                x = 2.0*nodes(i)-1; // x-coordinate
-               minBounds(i) -= a0V(j) + a1V(j)*x;
-               maxBounds(i) -= a0V(j) + a1V(j)*x;
+               minNodalVals(i) -= a0V(j) + a1V(j)*x;
+               maxNodalVals(i) -= a0V(j) + a1V(j)*x;
             }
             // Compute Bernstein coefficients
             LUFactors lu(basisMatLU.GetData(), lu_ip.GetData());
-            lu.Solve(nb, 1, minBounds.GetData());
-            lu.Solve(nb, 1, maxBounds.GetData());
+            lu.Solve(nb, 1, minNodalVals.GetData());
+            lu.Solve(nb, 1, maxNodalVals.GetData());
             for (int i = 0; i < nb; i++)
             {
-               intminT(i*ncp2+j) = minBounds(i);
-               intmaxT(i*ncp2+j) = maxBounds(i);
+               intminT(i*ncp2+j) = minNodalVals(i);
+               intmaxT(i*ncp2+j) = maxNodalVals(i);
             }
          }
       }
@@ -617,10 +622,10 @@ void PLBound::Get3DBounds(const Vector &coeff, Vector &intmin,
          real_t w1 = intmaxT(id2++);
          for (int k = 0; k < ncp; k++) // kth slice
          {
-            vals(0) = w0*lbound(j,k);
-            vals(1) = w0*ubound(j,k);
-            vals(2) = w1*lbound(j,k);
-            vals(3) = w1*ubound(j,k);
+            vals(0) = w0*lbound(k,j);
+            vals(1) = w0*ubound(k,j);
+            vals(2) = w1*lbound(k,j);
+            vals(3) = w1*ubound(k,j);
             intmin(k*ncp2+i) += vals.Min();
             intmax(k*ncp2+i) += vals.Max();
          }
@@ -653,7 +658,8 @@ void PLBound::SetupBernsteinBasisMat(DenseMatrix &basisMat,
                                      Vector &nodesBern) const
 {
    const int nbern = nodesBern.Size();
-   L2_SegmentElement el(nbern-1, 2); // we use L2 to leverage lexicographic order
+   L2_SegmentElement el(nbern-1, 2);
+   // we use L2 to leverage lexicographic order
    Array<int> ordering = el.GetLexicographicOrdering();
    basisMat.SetSize(nbern, nbern);
    Vector shape(nbern);
@@ -666,6 +672,39 @@ void PLBound::SetupBernsteinBasisMat(DenseMatrix &basisMat,
    }
 }
 
+DenseMatrix PLBound::GetBoundingMatrix(int dim, bool is_lower) const
+{
+   if (dim > 1)
+   {
+      const int ncpd = static_cast<int>(std::pow(ncp, dim));
+      const int nbd = static_cast<int>(std::pow(nb, dim));
+      DenseMatrix boundND(ncpd, nbd);
+      Vector phimin, phimax, col;
+      Vector coeffs(nbd);
+      coeffs = 0.0;
+      for (int j = 0; j < nbd; j++)
+      {
+         coeffs(j) = 1.0;
+         boundND.GetColumnReference(j, col);
+         GetNDBounds(dim, coeffs, phimin, phimax);
+         col = is_lower ? phimin : phimax;
+         coeffs(j) = 0.0;
+      }
+      return boundND;
+   }
+   return is_lower ? lbound : ubound;
+}
+
+DenseMatrix PLBound::GetLowerBoundMatrix(int dim) const
+{
+   return GetBoundingMatrix(dim, true);
+}
+
+DenseMatrix PLBound::GetUpperBoundMatrix(int dim) const
+{
+   return GetBoundingMatrix(dim, false);
+}
+
 constexpr int PLBound::min_ncp_gl_x[2][11];
 constexpr int PLBound::min_ncp_gll_x[2][11];
 constexpr int PLBound::min_ncp_pos_x[2][11];
@@ -716,4 +755,4 @@ void PLBound::Print(std::ostream &outp) const
    ubound.Print(outp);
 }
 
-}
+}

fem/bounds.hpp

@@ -19,14 +19,18 @@ namespace mfem
 {
 
 /** @name Piecewise linear bounds of bases
-    \brief Piecewise linear bounds of bases can be used to compute bounds on the grid function in each element. The bounds for the bases are constructed based on the following parameters:
+    \brief Piecewise linear bounds of bases can be used to compute bounds on
+    the grid function in each element. The bounds for the bases are constructed
+    based on the following parameters:
 
     (i) @b nb: number of bases/nodes in 1D (i.e. polynomial order+1),
 
-    (ii) @b b_type: bases type, 0 - Lagrange interpolants on Gauss-Legendre nodes, 1 - Lagrange interpolants on Gauss-Lobatto-Legendre nodes, and
+    (ii) @b b_type: bases type, 0 - Lagrange interpolants on Gauss-Legendre
+    nodes, 1 - Lagrange interpolants on Gauss-Lobatto-Legendre nodes, and
     2 - Positive/Bernstein bases on uniformly distributed nodes,
 
-    (iii) @b ncp: number of control points used to construct the piecewise linear bounds
+    (iii) @b ncp: number of control points used to construct the piecewise
+    linear bounds
 
     (iv) @b cp_type: control point distribution. 0 - GL + end-points,
                     1 - Chebyshev.
@@ -35,7 +39,9 @@ namespace mfem
 
     If the user does not specify @b ncp and @b cp_type, the minimum value of
     @b ncp is used that would bound the bases for the @b cp_type. We default
-    to @b cp_type = 0 as it requires fewer number of points to bound the bases. Typically, @b ncp = 2 @b nb is sufficient to get fairly compact bounds, and increasing @b ncp results in tighter bounds.
+    to @b cp_type = 0 as it requires fewer number of points to bound the bases.
+    Typically, @b ncp = 2 @b nb is sufficient to get fairly compact bounds, and
+    increasing @b ncp results in tighter bounds.
 
     Finally, only tensor-product elements are currently supported.
 
@@ -54,7 +60,7 @@ private:
    bool proj = true; // Use linear projection to compute bounds.
    real_t tol = 0.0; // offset bounds to avoid round-off errors
    Vector nodes, weights, control_points;
-   DenseMatrix lbound, ubound; // nb x ncp matrices with bounds of all bases
+   DenseMatrix lbound, ubound; // ncp x nb matrices with bounds of all bases
    // Some auxillary storage for computing the bounds with Bernstein
    DenseMatrix basisMatNodes; // Bernstein bases at equispaced nodes
    DenseMatrix basisMatInt;   // Bernstein bases at GLL nodes
@@ -80,6 +86,9 @@ private:
       {3,5,8,9,11,12,13,13,14,15,16}
    };
 
+   /// Helper function to extract lower or upper bounding matrix
+   DenseMatrix GetBoundingMatrix(int dim, bool is_lower) const;
+
 public:
    // Constructor
    PLBound(const int nb_i, const int ncp_i, const int b_type_i,
@@ -92,40 +101,82 @@ public:
    PLBound(const FiniteElementSpace *fes,
            const int ncp_i = -1, const int cp_type_i = 0);
 
-   // Get minimum number of control points needed to bound the given bases
+   /// Get minimum number of control points needed to bound the given bases
    int GetMinimumPointsForGivenBases(int nb_i, int b_type_i,
                                      int cp_type_i) const;
 
-   // Print information about the bounds
+   /// Print information about the bounds
    void Print(std::ostream &outp = mfem::out) const;
 
-   // Enable (default) or disable linear projection before bounding.
-   // This projection increases the computational cost but results in tighter
-   // bounds.
+   /** @brief Enable (default) or disable linear projection before bounding.
+    *
+    *  @details This projection increases the computational cost but results in
+    *  tighter bounds.
+    */
    void SetProjectionFlagForBounding(bool proj_) { proj = proj_; }
 
-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 1D/2D/3D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 1D/2D/3D.
+    *
+    *  @param[in] rdim    The spatial dimension of the element (1, 2, or 3).
+    *  @param[in] coeff   The vector of lexicographically-ordered coefficients.
+    *                     Should be of size nb^rdim, where nb is the number of
+    *                     bases/nodes in 1D. These coefficients must correspond
+    *                     to the bases type and number of bases, used in the
+    *                     constructor of PLBound.
+    *
+    *  @param[out] intmin The vector of minimum bound for all control points.
+    *  @param[out] intmax The vector of maximum bound for all control points.
+    *                     Both intmin and intmax are of size ncp^rdim, where
+    *                     ncp is the number of control points in 1D, and are
+    *                     ordered lexicographically.
+    */
    void GetNDBounds(const int rdim, const Vector &coeff,
                     Vector &intmin, Vector &intmax) const;
 
    /// Get number of control points used to compute the bounds.
    int GetNControlPoints() const { return ncp; }
+
+   /// Get 1D control point locations (lexicographic order) in [0,1].
+   const Vector &GetControlPoints() const { return control_points; }
+
+   /** @brief Get lower and upper bounding matrix (ncp^dim x nb^dim)
+    *
+    *  @details The matrices can be used to compute the bounds at control points
+    *           by a simple matrix-vector product with the
+    *           lexicographically-ordered nodal coefficients.
+    *           The resulting output is also lexicographically-ordered.
+    *
+    *  @note These matrices do not account for the linear projection step that
+    *        is optionally done in GetNDBounds before bounding the function.
+    */
+   ///@{
+   DenseMatrix GetLowerBoundMatrix(int dim = 1) const;
+   DenseMatrix GetUpperBoundMatrix(int dim = 1) const;
+   ///@}
+
 private:
-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 1D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 1D.
+    *  See GetNDBounds for details of the input and output parameters.
+    */
    void Get1DBounds(const Vector &coeff, Vector &intmin, Vector &intmax) const;
 
-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 2D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 2D.
+    *  See GetNDBounds for details of the input and output parameters.
+    */
    void Get2DBounds(const Vector &coeff, Vector &intmin, Vector &intmax) const;
 
-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 3D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 3D.
+    *  See GetNDBounds for details of the input and output parameters.
+    */
    void Get3DBounds(const Vector &coeff, Vector &intmin, Vector &intmax) const;
 
-   /// Setup matrix used to compute values at given 1D locations in [0,1]
-   /// for Bernstein bases.
+   /** @brief Setup matrix used to compute values at given 1D locations in [0,1]
+    *  for Bernstein bases.
+    */
    void SetupBernsteinBasisMat(DenseMatrix &basisMat, Vector &nodesBern) const;
 
    void Setup(const int nb_i, const int ncp_i, const int b_type_i,

fem/fe/fe_h1.cpp

@@ -589,7 +589,7 @@ void H1_TriangleElement::CalcHessian(const IntegrationPoint &ip,
    Vector   shape_x(p + 1),   shape_y(p + 1),   shape_l(p + 1);
    Vector  dshape_x(p + 1),  dshape_y(p + 1),  dshape_l(p + 1);
    Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
-   DenseMatrix ddu(dof, dim);
+   DenseMatrix ddu(dof, (dim*(dim+1))/2);
 #endif
 
    poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);

fem/gridfunc.cpp

@@ -30,6 +30,7 @@
 #include <cmath>
 #include <iostream>
 #include <algorithm>
+#include <queue>
 
 namespace mfem
 {
@@ -5117,6 +5118,103 @@ void GridFunction::GetElementBoundsAtControlPoints(const int elem,
    }
 }
 
+void GridFunction::GetElementBoundsAtControlPoints(const int elem,
+                                                   const PLBound &plb,
+                                                   const Vector &ref_range,
+                                                   const int vdim,
+                                                   Vector &lower, Vector &upper,
+                                                   Vector &control_pos) const
+{
+   const FiniteElement *fe = fes->GetFE(elem);
+   const IntegrationRule ir_in = fe->GetNodes();
+   IntegrationRule ir_new(ir_in.GetNPoints());
+   const int dim = fes->GetMesh()->Dimension();
+   const L2_FECollection *l2fec = dynamic_cast<const L2_FECollection *>
+                                  (fes->FEColl());
+
+   const TensorBasisElement *tbe =
+      dynamic_cast<const TensorBasisElement *>(fe);
+   MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
+
+   const Array<int> &dof_map = tbe->GetDofMap();
+   bool lexico = (dof_map.Size() == 0);
+   bool bern = (tbe->GetBasisType() == BasisType::Positive);
+   bool h1   = (l2fec == nullptr);
+
+   Vector loc_data; // gridfunction values
+   // Construct an integration rule to evaluate the gridfunction in
+   // subinterval.
+   for (int i = 0; i < ir_in.GetNPoints(); i++)
+   {
+      IntegrationPoint &ip_new = ir_new.IntPoint(i);
+      const IntegrationPoint &ip_old =
+         ir_in.IntPoint((lexico || bern) ? i : dof_map[i]);
+      Vector ip_coord(dim);
+      ip_old.Get(ip_coord.GetData(), dim);
+      for (int d = 0; d < dim; d++)
+      {
+         ip_coord(d) = ref_range(d) +
+                       (ref_range(dim+d) - ref_range(d)) * ip_coord(d);
+      }
+      ip_new.Set(ip_coord.GetData(), dim);
+   }
+   GetValues(elem, ir_new, loc_data, vdim);
+   // At this point, the loc_data contains function values ordered
+   // lexicographically, unless we are using Bernstein bases.
+   // For Bernstein, we need to project and get coefficients first.
+
+   // For bernstein, we get coefficients corresponding to these function values
+   if (bern)
+   {
+      int bt = 4; // BasisType::ClosedUniform
+      int o = fe->GetOrder();
+      DenseMatrix projmat;
+      NodalTensorFiniteElement *ntfe = nullptr;
+      if (dim == 1)
+      {
+         if (h1) { ntfe = new H1_SegmentElement(o, bt); }
+         else    { ntfe = new L2_SegmentElement(o, bt); }
+      }
+      else if (dim == 2)
+      {
+         if (h1) { ntfe = new H1_QuadrilateralElement(o, bt); }
+         else    { ntfe = new L2_QuadrilateralElement(o, bt); }
+      }
+      else if (dim == 3)
+      {
+         if (h1) { ntfe = new H1_HexahedronElement(o, bt); }
+         else    { ntfe = new L2_HexahedronElement(o, bt); }
+      }
+      // projection matrix from H1 to Positive
+      ElementTransformation *eltran = fes->GetElementTransformation(elem);
+      fe->Project(*ntfe, *eltran, projmat);
+      Vector loc_data_temp(loc_data.Size());
+      projmat.Mult(loc_data, loc_data_temp);
+      for (int i = 0; i < dof_map.Size(); i++)
+      {
+         loc_data(i) = loc_data_temp(dof_map[i]);
+      }
+      if (dof_map.Size() == 0) { loc_data = loc_data_temp; }
+      delete ntfe;
+   }
+
+   // Get bounds at control points
+   plb.GetNDBounds(dim, loc_data, lower, upper);
+
+   // Save control point positions
+   int ncp = plb.GetNControlPoints();
+   control_pos.SetSize(dim * ncp);
+   const Vector control_pos_1D = plb.GetControlPoints();
+   for (int i = 0; i < ncp; i++)
+   {
+      for (int d = 0; d < dim; d++)
+      {
+         control_pos(i + d*ncp) =
+            ref_range(d) + (ref_range(dim+d)-ref_range(d))*control_pos_1D(i);
+      }
+   }
+}
+
 void GridFunction::GetElementBounds(const int elem, const PLBound &plb,
                                     Vector &lower, Vector &upper,
                                     const int vdim) const
@@ -5197,6 +5295,467 @@ PLBound GridFunction::GetBounds(Vector &lower, Vector &upper,
    return plb;
 }
 
+struct IntervalNode
+{
+   real_t val_min;
+   real_t val_max;
+   Array<IntervalNode *> child;
+   IntervalNode(real_t vmin, real_t vmax)
+      : val_min(vmin), val_max(vmax)
+   {
+      child.SetSize(0);
+   }
+   void AddChild(IntervalNode *ch) { child.Append(ch); }
+   real_t GetChildMinLower()
+   {
+      if (child.Size() == 0)
+      {
+         return val_min;
+      }
+      real_t valmin = numeric_limits<real_t>::max();
+      for (int i = 0; i < child.Size(); i++)
+      {
+         real_t candidate = child[i]->GetChildMinLower();
+         valmin = std::min(valmin, candidate);
+      }
+      return valmin;
+   }
+   real_t GetChildMinUpper()
+   {
+      if (child.Size() == 0)
+      {
+         return val_max;
+      }
+      real_t valmax = numeric_limits<real_t>::max();
+      for (int i = 0; i < child.Size(); i++)
+      {
+         real_t candidate = child[i]->GetChildMinUpper();
+         valmax = std::min(valmax, candidate);
+      }
+      return valmax;
+   }
+   real_t GetChildMaxLower()
+   {
+      if (child.Size() == 0)
+      {
+         return val_min;
+      }
+      real_t valmin = numeric_limits<real_t>::lowest();
+      for (int i = 0; i < child.Size(); i++)
+      {
+         real_t candidate = child[i]->GetChildMaxLower();
+         valmin = std::max(valmin, candidate);
+      }
+      return valmin;
+   }
+   real_t GetChildMaxUpper()
+   {
+      if (child.Size() == 0)
+      {
+         return val_max;
+      }
+      real_t valmax = numeric_limits<real_t>::lowest();
+      for (int i = 0; i < child.Size(); i++)
+      {
+         real_t candidate = child[i]->GetChildMaxUpper();
+         valmax = std::max(valmax, candidate);
+      }
+      return valmax;
+   }
+   void DeleteChildren()
+   {
+      for (int i = 0; i < child.Size(); i++)
+      {
+         child[i]->DeleteChildren();
+         delete child[i];
+      }
+      child.SetSize(0);
+   }
+};
+
+struct SearchInterval
+{
+   Vector ref_range;
+   int depth;
+   IntervalNode *node;
+   SearchInterval(const Vector &ref_range_in, int d, IntervalNode *n)
+      : ref_range(ref_range_in), depth(d), node(n)
+   { }
+};
+
+struct IntervalCompareMin
+{
+   bool operator()(const SearchInterval *a, const SearchInterval *b) const
+   {
+      return a->node->val_min > b->node->val_min;
+   }
+};
+
+struct IntervalCompareMax
+{
+   bool operator()(const SearchInterval *a, const SearchInterval *b) const
+   {
+      return a->node->val_max < b->node->val_max;
+   }
+};
+
+std::pair<real_t, real_t> GridFunction::EstimateFunctionMinimum(
+   const int elem, const PLBound &plb, const int vdim,
+   const int max_depth, const real_t tol) const
+{
+   real_t min_threshold = std::numeric_limits<real_t>::max();
+   return EstimateFunctionMinimum(elem, plb, vdim, max_depth, tol,
+                                  min_threshold);
 }
 
+std::pair<real_t, real_t> GridFunction::EstimateFunctionMinimum(
+   const int elem, const PLBound &plb, const int vdim,
+   const int max_depth, const real_t tol, real_t &min_threshold) const
+{
+   const int dim = this->FESpace()->GetMesh()->Dimension();
+   const int ncp = plb.GetNControlPoints();
+   Vector pos_range(2*dim); pos_range = 0.0;
+   for (int d = 0; d < dim; d++) { pos_range(d+dim) = 1.0; }
+   Vector lower, upper, cp_ref_loc;
 
+   GetElementBoundsAtControlPoints(elem, plb, lower, upper, vdim);
+   real_t val_min = lower.Min();
+   real_t val_max = upper.Min();
+
+   min_threshold = std::min(min_threshold, val_max);
+
+   // Pruning: if the element's lower bound is greater than the current global
+   // upper bound, this element cannot contain the global minimum.
+   if (val_min >= min_threshold)
+   {
+      return std::make_pair(val_min, val_max);
+   }
+
+   if (val_min == val_max || max_depth == 0)
+   {
+      min_threshold = std::min(min_threshold, val_min);
+      return std::make_pair(val_min, val_max);
+   }
+   real_t abs_tol = tol*(val_max-val_min);
+
+   IntervalNode *initial_node = new IntervalNode(val_min, val_max);
+   SearchInterval *initial_interval = new SearchInterval(pos_range, 0,
+                                                         initial_node);
+
+   std::priority_queue<SearchInterval*,
+       std::vector<SearchInterval*>, IntervalCompareMin> pq;
+   pq.push(initial_interval);
+
+   real_t min_upper_bound = upper.Min();
+   real_t min_lower_bound = lower.Min();
+
+   while (!pq.empty())
+   {
+      SearchInterval *current = pq.top();
+      pq.pop();
+      int curr_depth = current->depth;
+
+      // Reached max depth or this interval cannot contain the global minimum
+      if (current->node->val_min >= min_threshold || curr_depth >= max_depth)
+      {
+         delete current;
+         continue;
+      }
+
+      min_lower_bound = initial_node->GetChildMinLower();
+      if (min_upper_bound - min_lower_bound < abs_tol)
+      {
+         delete current;
+         break;
+      }
+
+      // Subdivide the interval and get bounds on it
+      GetElementBoundsAtControlPoints(elem, plb, current->ref_range,
+                                      vdim, lower, upper, cp_ref_loc);
+
+      // process the bounds and create sub-intervals
+      for (int k = 0; k < (dim == 3 ? ncp-1 : 1); k++)
+      {
+         for (int j = 0; j < (dim >= 2 ? ncp-1 : 1); j++)
+         {
+            for (int i = 0; i < ncp-1; i++)
+            {
+               real_t lv = 0.0, uv = 0.0;
+               if (dim == 1)
+               {
+                  lv = std::min(lower(i), lower(i+1));
+                  uv = std::min(upper(i), upper(i+1));
+               }
+               else if (dim == 2)
+               {
+                  lv = std::min({lower(i + j*ncp), lower((i+1) + j*ncp),
+                                 lower(i + (j+1)*ncp),
+                                 lower((i+1) + (j+1)*ncp)});
+                  uv = std::min({upper(i + j*ncp), upper((i+1) + j*ncp),
+                                 upper(i + (j+1)*ncp),
+                                 upper((i+1) + (j+1)*ncp)});
+               }
+               else if (dim == 3)
+               {
+                  lv = std::min({lower(i + j*ncp + k*ncp*ncp),
+                                 lower((i+1) + j*ncp + k*ncp*ncp),
+                                 lower(i + (j+1)*ncp + k*ncp*ncp),
+                                 lower((i+1) + (j+1)*ncp + k*ncp*ncp),
+                                 lower(i + j*ncp + (k+1)*ncp*ncp),
+                                 lower((i+1) + j*ncp + (k+1)*ncp*ncp),
+                                 lower(i + (j+1)*ncp + (k+1)*ncp*ncp),
+                                 lower((i+1) + (j+1)*ncp + (k+1)*ncp*ncp)});
+                  uv = std::min({upper(i + j*ncp + k*ncp*ncp),
+                                 upper((i+1) + j*ncp + k*ncp*ncp),
+                                 upper(i + (j+1)*ncp + k*ncp*ncp),
+                                 upper((i+1) + (j+1)*ncp + k*ncp*ncp),
+                                 upper(i + j*ncp + (k+1)*ncp*ncp),
+                                 upper((i+1) + j*ncp + (k+1)*ncp*ncp),
+                                 upper(i + (j+1)*ncp + (k+1)*ncp*ncp),
+                                 upper((i+1) + (j+1)*ncp + (k+1)*ncp*ncp)});
+               }
+               IntervalNode *child_node = new IntervalNode(lv, uv);
+               current->node->AddChild(child_node);
+
+               if (lv < min_threshold)
+               {
+                  min_upper_bound = std::min(min_upper_bound, uv);
+                  min_threshold = std::min(min_threshold, uv);
+                  if (curr_depth < max_depth)
+                  {
+                     pos_range(0) = cp_ref_loc(i);
+                     pos_range(0+dim) = cp_ref_loc(i+1);
+                     if (dim >= 2)
+                     {
+                        pos_range(1) = cp_ref_loc(ncp + j);
+                        pos_range(1+dim) = cp_ref_loc(ncp + j+1);
+                     }
+                     if (dim == 3)
+                     {
+                        pos_range(2) = cp_ref_loc(2*ncp + k);
+                        pos_range(2+dim) = cp_ref_loc(2*ncp + k+1);
+                     }
+                     SearchInterval *child_interval =
+                        new  SearchInterval(pos_range, curr_depth + 1,
+                                            child_node);
+                     pq.push(child_interval);
+                  }
+               }
+            }
+         }
+      }
+      delete current;
+   }
+
+   // clean up remaining intervals in queue
+   while (!pq.empty())
+   {
+      delete pq.top();
+      pq.pop();
+   }
+
+   min_lower_bound = initial_node->GetChildMinLower();
+   initial_node->DeleteChildren();
+   delete initial_node;
+
+   min_threshold = std::min(min_threshold, min_lower_bound);
+   return std::make_pair(min_lower_bound, min_upper_bound);
+}
+
+std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
+   const int elem, const PLBound &plb, const int vdim,
+   const int max_depth, const real_t tol) const
+{
+   real_t max_threshold = std::numeric_limits<real_t>::lowest();
+   return EstimateFunctionMaximum(elem, plb, vdim, max_depth, tol,
+                                  max_threshold);
+}
+
+std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
+   const int elem, const PLBound &plb, const int vdim,
+   const int max_depth, const real_t tol, real_t &max_threshold) const
+{
+   const int dim = this->FESpace()->GetMesh()->Dimension();
+   const int ncp = plb.GetNControlPoints();
+   Vector pos_range(2*dim); pos_range = 0.0;
+   for (int d = 0; d < dim; d++) { pos_range(d+dim) = 1.0; }
+   Vector lower, upper, cp_ref_loc;
+
+   GetElementBoundsAtControlPoints(elem, plb, lower, upper, vdim);
+   real_t val_min = lower.Max();
+   real_t val_max = upper.Max();
+
+   max_threshold = std::max(max_threshold, val_min);
+
+   // Pruning: if the element's upper bound is less than the current global
+   // lower bound, this element cannot contain the global maximum.
+   if (val_max <= max_threshold)
+   {
+      return std::make_pair(val_min, val_max);
+   }
+
+   if (val_min == val_max || max_depth == 0)
+   {
+      max_threshold = std::max(max_threshold, val_max);
+      return std::make_pair(val_min, val_max);
+   }
+   real_t abs_tol = tol*(val_max-val_min);
+
+   IntervalNode *initial_node = new IntervalNode(val_min, val_max);
+   SearchInterval *initial_interval = new SearchInterval(pos_range, 0,
+                                                         initial_node);
+
+   std::priority_queue<SearchInterval*,
+       std::vector<SearchInterval*>, IntervalCompareMax> pq;
+   pq.push(initial_interval);
+
+   real_t max_lower_bound = val_min;
+   real_t max_upper_bound = val_max;
+
+   while (!pq.empty())
+   {
+      SearchInterval *current = pq.top();
+      pq.pop();
+      int curr_depth = current->depth;
+
+      // Reached max depth or this interval cannot contain the global maximum.
+      if (current->node->val_max <= max_threshold || curr_depth >= max_depth)
+      {
+         delete current;
+         continue;
+      }
+
+      max_upper_bound = initial_node->GetChildMaxUpper();
+      if (max_upper_bound - max_lower_bound < abs_tol)
+      {
+         delete current;
+         break;
+      }
+
+      // Subdivide the interval and get bounds on it
+      GetElementBoundsAtControlPoints(elem, plb, current->ref_range,
+                                      vdim, lower, upper, cp_ref_loc);
+
+      // process the bounds and create sub-intervals
+      for (int k = 0; k < (dim == 3 ? ncp-1 : 1); k++)
+      {
+         for (int j = 0; j < (dim >= 2 ? ncp-1 : 1); j++)
+         {
+            for (int i = 0; i < ncp-1; i++)
+            {
+               real_t lv = 0.0, uv = 0.0;
+               if (dim == 1)
+               {
+                  lv = std::max(lower(i), lower(i+1));
+                  uv = std::max(upper(i), upper(i+1));
+               }
+               else if (dim == 2)
+               {
+                  lv = std::max({lower(i + j*ncp), lower((i+1) + j*ncp),
+                                 lower(i + (j+1)*ncp),
+                                 lower((i+1) + (j+1)*ncp)});
+                  uv = std::max({upper(i + j*ncp), upper((i+1) + j*ncp),
+                                 upper(i + (j+1)*ncp),
+                                 upper((i+1) + (j+1)*ncp)});
+               }
+               else if (dim == 3)
+               {
+                  lv = std::max({lower(i + j*ncp + k*ncp*ncp),
+                                 lower((i+1) + j*ncp + k*ncp*ncp),
+                                 lower(i + (j+1)*ncp + k*ncp*ncp),
+                                 lower((i+1) + (j+1)*ncp + k*ncp*ncp),
+                                 lower(i + j*ncp + (k+1)*ncp*ncp),
+                                 lower((i+1) + j*ncp + (k+1)*ncp*ncp),
+                                 lower(i + (j+1)*ncp + (k+1)*ncp*ncp),
+                                 lower((i+1) + (j+1)*ncp + (k+1)*ncp*ncp)});
+                  uv = std::max({upper(i + j*ncp + k*ncp*ncp),
+                                 upper((i+1) + j*ncp + k*ncp*ncp),
+                                 upper(i + (j+1)*ncp + k*ncp*ncp),
+                                 upper((i+1) + (j+1)*ncp + k*ncp*ncp),
+                                 upper(i + j*ncp + (k+1)*ncp*ncp),
+                                 upper((i+1) + j*ncp + (k+1)*ncp*ncp),
+                                 upper(i + (j+1)*ncp + (k+1)*ncp*ncp),
+                                 upper((i+1) + (j+1)*ncp + (k+1)*ncp*ncp)});
+               }
+               IntervalNode *child_node = new IntervalNode(lv, uv);
+               current->node->AddChild(child_node);
+
+               if (uv > max_threshold)
+               {
+                  max_lower_bound = std::max(max_lower_bound, lv);
+                  max_threshold = std::max(max_threshold, lv);
+                  if (curr_depth < max_depth)
+                  {
+                     pos_range(0) = cp_ref_loc(i);
+                     pos_range(0+dim) = cp_ref_loc(i+1);
+                     if (dim >= 2)
+                     {
+                        pos_range(1) = cp_ref_loc(ncp + j);
+                        pos_range(1+dim) = cp_ref_loc(ncp + j+1);
+                     }
+                     if (dim == 3)
+                     {
+                        pos_range(2) = cp_ref_loc(2*ncp + k);
+                        pos_range(2+dim) = cp_ref_loc(2*ncp + k+1);
+                     }
+                     SearchInterval *child_interval =
+                        new SearchInterval(pos_range, curr_depth + 1,
+                                           child_node);
+                     pq.push(child_interval);
+                  }
+               }
+            }
+         }
+      }
+      delete current;
+   }
+   // clean up remaining intervals in queue
+   while (!pq.empty())
+   {
+      delete pq.top();
+      pq.pop();
+   }
+
+   max_upper_bound = initial_node->GetChildMaxUpper();
+   initial_node->DeleteChildren();
+   delete initial_node;
+   max_threshold = std::max(max_threshold, max_upper_bound);
+
+   return std::make_pair(max_lower_bound, max_upper_bound);
+}
+
+std::pair<real_t, real_t> GridFunction::EstimateFunctionMinimum(
+   const int vdim, const PLBound &plb, const int max_depth,
+   const real_t tol) const
+{
+   real_t global_min_lower = std::numeric_limits<real_t>::max();
+   real_t global_min_upper = std::numeric_limits<real_t>::max();
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      std::pair<real_t, real_t> min_pair =
+         EstimateFunctionMinimum(i, plb, vdim, max_depth, tol,
+                                 global_min_lower);
+      global_min_upper = std::min(global_min_upper, min_pair.second);
+   }
+   return std::make_pair(global_min_lower, global_min_upper);
+}
+
+std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
+   const int vdim, const PLBound &plb, const int max_depth,
+   const real_t tol) const
+{
+   real_t global_max_lower = std::numeric_limits<real_t>::lowest();
+   real_t global_max_upper = std::numeric_limits<real_t>::lowest();
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      std::pair<real_t, real_t> max_pair =
+         EstimateFunctionMaximum(i, plb, vdim, max_depth, tol,
+                                 global_max_upper);
+      global_max_lower = std::max(global_max_lower, max_pair.first);
+   }
+   return std::make_pair(global_max_lower, global_max_upper);
+}
+
+}

fem/gridfunc.hpp

@@ -564,6 +564,70 @@ protected:
    /// P-refinement version of Update().
    void UpdatePRef();
 
+   /** @brief Estimate the minimum value of the GridFunction in element @a elem
+    *  if it is below a certain @a min_threshold.
+    *
+    *  @details For a given element \p elem and grid function component \p vdim
+    *  an estimate of the function minimum is the minimum of the piecewise
+    *  linear lower bound obtained using the given PLBound object. The actual
+    *  minimum is between [minimum lower bound, minimum upper bound]. We
+    *  improve the estimate of the function minimum by recursively
+    *  subdividing the interval with the lowest lower bound, and computing
+    *  bounds on the sub-intervals.
+    *  This process continues until (i) the maximum recursion depth is reached
+    *  or (ii) the difference between the minimum upper bound and minimum lower
+    *  bound is less than a certain tolerance (\p tol * [initial maximum
+    *  upper bound - initial minimum lower bound]).
+    *  The function also terminates if the lowest minima estimate is found
+    *  to be above the given threshold \p min_threshold. This is useful when
+    *  we are interested in computing the global minimum of the function
+    *  over all elements. In this case we can reject elements where the lowest
+    *  bound is above the current global minimum. In case the function
+    *  minimum on the element is below the global minimum, we update
+    *  \p min_threshold.
+    *
+    *  We return a pair of values that bracket the actual minimum, i.e.
+    *  [min_lower_bound, min_upper_bound].
+    */
+   std::pair<real_t,real_t> EstimateFunctionMinimum(const int elem,
+                                                    const PLBound &plb,
+                                                    const int vdim,
+                                                    const int max_depth,
+                                                    const real_t tol,
+                                                    real_t &min_threshold)const;
+
+   /** @brief Estimate the maximum value of the GridFunction in element @a elem
+    *  if it is below a certain @a max_threshold.
+    *
+    *  @details For a given element \p elem and grid function component \p vdim
+    *  an estimate of the function maximum is the maximum of the piecewise
+    *  linear upper bound obtained using the given PLBound object. The actual
+    *  maximum is between [maximum lower bound, maximum upper bound]. We
+    *  improve the estimate of the function maximum by recursively
+    *  subdividing the interval with the highest upper bound, and computing
+    *  bounds on the sub-intervals.
+    *  This process continues until (i) the maximum recursion depth is reached
+    *  or (ii) the difference between the maximum upper bound and maximum lower
+    *  bound is less than a certain tolerance (\p tol * [initial maximum
+    *  upper bound - initial maximum lower bound]).
+    *  The function also terminates if the highest maxima estimate is found
+    *  to be below the given threshold \p max_threshold. This is useful when
+    *  we are interested in computing the global maximum of the function
+    *  over all elements. In this case we can reject elements where the upper
+    *  bound is below the current global maximum. In case the function
+    *  maximum on the element is above the global maximum, we update
+    *  \p max_threshold.
+    *
+    *  We return a pair of values that bracket the actual maximum, i.e.
+    *  [max_lower_bound, max_upper_bound].
+    */
+   std::pair<real_t,real_t> EstimateFunctionMaximum(const int elem,
+                                                    const PLBound &plb,
+                                                    const int vdim,
+                                                    const int max_depth,
+                                                    const real_t tol,
+                                                    real_t &max_threshold)const;
+
 public:
    /** @brief For each vdof, counts how many elements contain the vdof,
        as containment is determined by FiniteElementSpace::GetElementVDofs(). */
@@ -1662,21 +1726,21 @@ public:
     */
    ///@{
    /// Computes the \ref PLBound for the gridfunction with number of control
-   /// points based on @a ref_factor, and returns the overall bounds for each
-   /// vdim (across all elements) in @b lower and @b upper. We also return the
+   /// points based on \p ref_factor, and returns the overall bounds for each
+   /// vdim (across all elements) in \p lower and \p upper. We also return the
    /// PLBound object used to compute the bounds.
-   /// We compute the bounds for each vdim if @a vdim < 1.
+   /// We compute the bounds for each vdim if \p vdim < 1.
    /// Note: For most cases, this method/interface will be sufficient.
    virtual PLBound GetBounds(Vector &lower, Vector &upper,
                              const int ref_factor=1, const int vdim=-1) const;
 
    /// Computes the \ref PLBound for the gridfunction with number of control
-   /// points based on @a ref_factor, and returns the bounds for each element
-   /// ordered byVDim:
+   /// points based on \p ref_factor, and returns the bounds for each element
+   /// ordered byNodes:
    /// lower_{0,0}, lower_{1,0}, ..., lower_{ne-1,0},
    /// lower_{0,1}, ..., lower_{ne-1,vdim-1}. We also return the
    /// PLBound object used to compute the bounds.
-   /// We compute the bounds for each vdim if @a vdim < 1.
+   /// We compute the bounds for each vdim if \p vdim < 1.
    PLBound GetElementBounds(Vector &lower, Vector &upper,
                             const int ref_factor=1, const int vdim=-1) const;
 
@@ -1687,6 +1751,18 @@ public:
                                         Vector &lower, Vector &upper,
                                         const int vdim = -1) const;
 
+   /** @brief Gets the bounds on given reference range inside an element.
+    *
+    *  @details @a ref_range is a vector of size 2*dim that specifies the
+    *  lower and upper limits in each dimension of the reference element.
+    *  For example, in 2D, ref_range = [rmin, smin, rmax, smax].
+    */
+   void GetElementBoundsAtControlPoints(const int elem, const PLBound &plb,
+                                        const Vector &ref_range,
+                                        const int vdim,
+                                        Vector &lower, Vector &upper,
+                                        Vector &control_pos) const;
+
    /// Compute bounds on the grid function for the given element.
    /// The bounds are stored in @b lower and @b upper.
    void GetElementBounds(const int elem, const PLBound &plb,
@@ -1694,11 +1770,45 @@ public:
                          const int vdim = -1) const;
 
    /// Compute bounds on the grid function for all the elements. The bounds
-   /// are returned in @b lower and @b upper, ordered byVDim:
+   /// are returned in @b lower and @b upper, ordered byNodes:
    /// lower_{0,0}, lower_{1,0}, ..., lower_{ne-1,0},
    /// lower_{0,1}, ..., lower_{ne-1,vdim-1}
    void GetElementBounds(const PLBound &plb, Vector &lower, Vector &upper,
                          const int vdim=-1) const;
+
+   /** @brief Estimate the minimum value of the GridFunction in element @a elem.
+    *
+    *  @details See the protected version of EstimateFunctionMinimum for
+    *  details.
+    */
+   std::pair<real_t, real_t> EstimateFunctionMinimum(const int elem,
+                                                     const PLBound &plb,
+                                                     const int vdim,
+                                                     const int max_depth,
+                                                     const real_t tol) const;
+
+   /** @brief Estimate the minimum value of the GridFunction in element @a elem.
+    *
+    *  @details See the protected version of EstimateFunctionMaximum for
+    *  details.
+    */
+   std::pair<real_t, real_t> EstimateFunctionMaximum(const int elem,
+                                                     const PLBound &plb,
+                                                     const int vdim,
+                                                     const int max_depth,
+                                                     const real_t tol) const;
+
+   /** @brief Estimate the GridFunction minimum across all elements. */
+   virtual std::pair<real_t,real_t> EstimateFunctionMinimum(const int vdim,
+                                                            const PLBound &plb,
+                                                            const int max_depth,
+                                                            const real_t tol) const;
+
+   /** @brief Estimate the GridFunction maximum across all elements. */
+   virtual std::pair<real_t,real_t> EstimateFunctionMaximum(const int vdim,
+                                                            const PLBound &plb,
+                                                            const int max_depth,
+                                                            const real_t tol) const;
    ///@}
 
    /// Destroys grid function.

fem/integ/bilininteg_vecdiffusion_pa.hpp

@@ -171,15 +171,15 @@ template<int DIM, int T_SDIM, int T_D1D, int T_Q1D>
 VectorDiffusionIntegrator::ApplyKernelType
 VectorDiffusionIntegrator::ApplyPAKernels::Kernel()
 {
-   if (DIM == 2)
+   if constexpr (DIM == 2)
    {
       return internal::SmemPAVectorDiffusionApply2D<T_SDIM, T_D1D, T_Q1D>;
    }
-   else if (DIM == 3)
+   else if constexpr (DIM == 3)
    {
       return internal::SmemPAVectorDiffusionApply3D<T_SDIM, T_D1D, T_Q1D>;
    }
-   else { MFEM_ABORT("Unsupported kernel"); }
+   MFEM_ABORT("Unsupported kernel");
 }
 
 inline VectorDiffusionIntegrator::ApplyKernelType

fem/integ/bilininteg_vecmass_pa.hpp

@@ -182,15 +182,15 @@ template<int DIM, int T_D1D, int T_Q1D>
 VectorMassIntegrator::VectorMassAddMultPAType
 VectorMassIntegrator::VectorMassAddMultPA::Kernel()
 {
-   if (DIM == 2)
+   if constexpr (DIM == 2)
    {
       return internal::SmemPAVectorMassApply2D<T_D1D,T_Q1D>;
    }
-   else if (DIM == 3)
+   else if constexpr (DIM == 3)
    {
       return internal::SmemPAVectorMassApply3D<T_D1D, T_Q1D>;
    }
-   else { MFEM_ABORT("Unsupported kernel"); }
+   MFEM_ABORT("Unsupported kernel");
 }
 
 inline VectorMassIntegrator::VectorMassAddMultPAType

fem/integ/lininteg_domain_kernels.hpp

@@ -301,18 +301,14 @@ template <int DIM, int T_D1D, int T_Q1D>
 DomainLFIntegrator::AssembleKernelType
 DomainLFIntegrator::AssembleKernels::Kernel()
 {
-   switch (DIM)
-   {
-      case 1:
-         return DLFEvalAssemble1D<T_D1D, T_Q1D>;
-      case 2:
-         return DLFEvalAssemble2D<T_D1D, T_Q1D>;
-      case 3:
-         return DLFEvalAssemble3D<T_D1D, T_Q1D>;
-   }
+   if constexpr (DIM == 1) { return DLFEvalAssemble1D<T_D1D, T_Q1D>; }
+   if constexpr (DIM == 2) { return DLFEvalAssemble2D<T_D1D, T_Q1D>; }
+   if constexpr (DIM == 3) { return DLFEvalAssemble3D<T_D1D, T_Q1D>; }
    MFEM_ABORT("");
 }
+
 /// \endcond DO_NOT_DOCUMENT
 
 } // namespace mfem
-#endif
+
+#endif // MFEM_LININTEG_DOMAIN_KERNELS_HPP

fem/pgridfunc.cpp

@@ -1568,6 +1568,39 @@ PLBound ParGridFunction::GetBounds(Vector &lower, Vector &upper,
    return plb;
 }
 
+std::pair<real_t, real_t> ParGridFunction::EstimateFunctionMinimum(
+   const int vdim, const PLBound &plb, const int max_depth,
+   const real_t tol) const
+{
+   std::pair<real_t, real_t> minmax =
+      GridFunction::EstimateFunctionMinimum(vdim, plb, max_depth, tol);
+
+   real_t glob_min_lower = minmax.first;
+   real_t glob_min_upper = minmax.second;
+   MPI_Allreduce(MPI_IN_PLACE, &glob_min_lower, 1,
+                 MFEM_MPI_REAL_T, MPI_MIN, pfes->GetComm());
+   MPI_Allreduce(MPI_IN_PLACE, &glob_min_upper, 1,
+                 MFEM_MPI_REAL_T, MPI_MIN, pfes->GetComm());
+
+   return std::make_pair(glob_min_lower, glob_min_upper);
+}
+
+std::pair<real_t, real_t> ParGridFunction::EstimateFunctionMaximum(
+   const int vdim, const PLBound &plb, const int max_depth,
+   const real_t tol) const
+{
+   std::pair<real_t, real_t> minmax =
+      GridFunction::EstimateFunctionMaximum(vdim, plb, max_depth, tol);
+
+   real_t glob_max_lower = minmax.first;
+   real_t glob_max_upper = minmax.second;
+   MPI_Allreduce(MPI_IN_PLACE, &glob_max_lower, 1,
+                 MFEM_MPI_REAL_T, MPI_MAX, pfes->GetComm());
+   MPI_Allreduce(MPI_IN_PLACE, &glob_max_upper, 1,
+                 MFEM_MPI_REAL_T, MPI_MAX, pfes->GetComm());
+   return std::make_pair(glob_max_lower, glob_max_upper);
+}
+
 } // namespace mfem
 
 #endif // MFEM_USE_MPI

fem/pgridfunc.hpp

@@ -609,6 +609,18 @@ public:
    PLBound GetBounds(Vector &lower, Vector &upper,
                      const int ref_factor=1, const int vdim=-1) const override;
 
+   /** @brief Estimate the GridFunction minimum across all elements. */
+   std::pair<real_t, real_t> EstimateFunctionMinimum(const int vdim,
+                                                     const PLBound &plb,
+                                                     const int max_depth,
+                                                     const real_t tol) const override;
+
+   /** @brief Estimate the GridFunction maximum across all elements. */
+   std::pair<real_t, real_t> EstimateFunctionMaximum(const int vdim,
+                                                     const PLBound &plb,
+                                                     const int max_depth,
+                                                     const real_t tol) const override;
+
    /** Save the local portion of the ParGridFunction. This differs from the
        serial GridFunction::Save in that it takes into account the signs of
        the local dofs. */

fem/qinterp/det.hpp

@@ -334,17 +334,16 @@ template<int DIM, int SDIM, int D1D, int Q1D>
 QuadratureInterpolator::DetKernelType
 QuadratureInterpolator::DetKernels::Kernel()
 {
-   if (DIM == 1)
+   if constexpr (DIM == 1)
    {
-      if (SDIM == 1) { return internal::quadrature_interpolator::Det1D; }
-      else if (SDIM == 2) { return internal::quadrature_interpolator::Det1DSurface<D1D, Q1D, 2>; }
-      else if (SDIM == 3) { return internal::quadrature_interpolator::Det1DSurface<D1D, Q1D, 3>; }
-      else { MFEM_ABORT(""); }
+      if constexpr (SDIM == 1) { return internal::quadrature_interpolator::Det1D; }
+      else if constexpr (SDIM == 2) { return internal::quadrature_interpolator::Det1DSurface<D1D, Q1D, 2>; }
+      else if constexpr (SDIM == 3) { return internal::quadrature_interpolator::Det1DSurface<D1D, Q1D, 3>; }
    }
-   else if (DIM == 2 && SDIM == 2) { return internal::quadrature_interpolator::Det2D<D1D, Q1D>; }
-   else if (DIM == 2 && SDIM == 3) { return internal::quadrature_interpolator::Det2DSurface<D1D, Q1D>; }
-   else if (DIM == 3) { return internal::quadrature_interpolator::Det3D<D1D, Q1D>; }
-   else { MFEM_ABORT(""); }
+   else if constexpr (DIM == 2 && SDIM == 2) { return internal::quadrature_interpolator::Det2D<D1D, Q1D>; }
+   else if constexpr (DIM == 2 && SDIM == 3) { return internal::quadrature_interpolator::Det2DSurface<D1D, Q1D>; }
+   else if constexpr (DIM == 3) { return internal::quadrature_interpolator::Det3D<D1D, Q1D>; }
+   MFEM_ABORT("");
 }
 
 /// @endcond

fem/qinterp/eval.hpp

@@ -203,10 +203,10 @@ template<int DIM, QVectorLayout Q_LAYOUT,
 QuadratureInterpolator::TensorEvalKernelType
 QuadratureInterpolator::TensorEvalKernels::Kernel()
 {
-   if (DIM == 1) { return internal::quadrature_interpolator::Values1D<Q_LAYOUT>; }
-   else if (DIM == 2) { return internal::quadrature_interpolator::Values2D<Q_LAYOUT, VDIM, D1D, Q1D, NBZ>; }
-   else if (DIM == 3) { return internal::quadrature_interpolator::Values3D<Q_LAYOUT, VDIM, D1D, Q1D>; }
-   else { MFEM_ABORT(""); }
+   if constexpr (DIM == 1) { return internal::quadrature_interpolator::Values1D<Q_LAYOUT>; }
+   else if constexpr (DIM == 2) { return internal::quadrature_interpolator::Values2D<Q_LAYOUT, VDIM, D1D, Q1D, NBZ>; }
+   else if constexpr (DIM == 3) { return internal::quadrature_interpolator::Values3D<Q_LAYOUT, VDIM, D1D, Q1D>; }
+   MFEM_ABORT("");
 }
 
 /// @endcond

fem/qinterp/eval_hdiv.hpp

@@ -453,8 +453,15 @@ QuadratureInterpolator::TensorEvalHDivKernels::Kernel()
 {
    using namespace internal::quadrature_interpolator;
    static_assert(DIM == 2 || DIM == 3, "only DIM=2 and DIM=3 are implemented!");
-   if (DIM == 2) { return EvalHDiv2D<Q_LAYOUT, FLAGS, D1D, Q1D>; }
-   return EvalHDiv3D<Q_LAYOUT, FLAGS, D1D, Q1D>;
+   if constexpr (DIM == 2)
+   {
+      return EvalHDiv2D<Q_LAYOUT, FLAGS, D1D, Q1D>;
+   }
+   else if constexpr (DIM == 3)
+   {
+      return EvalHDiv3D<Q_LAYOUT, FLAGS, D1D, Q1D>;
+   }
+   MFEM_ABORT("only DIM=2 and DIM=3 are implemented!");
 }
 
 /// @endcond

fem/qinterp/grad.hpp

@@ -592,10 +592,10 @@ template<int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int VDIM, int D1D,
 QuadratureInterpolator::GradKernelType
 QuadratureInterpolator::GradKernels::Kernel()
 {
-   if (DIM == 1) { return internal::quadrature_interpolator::Derivatives1D<Q_LAYOUT, GRAD_PHYS>; }
-   else if (DIM == 2) { return internal::quadrature_interpolator::Derivatives2D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D, NBZ>; }
-   else if (DIM == 3) { return internal::quadrature_interpolator::Derivatives3D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D>; }
-   else { MFEM_ABORT(""); }
+   if constexpr (DIM == 1) { return internal::quadrature_interpolator::Derivatives1D<Q_LAYOUT, GRAD_PHYS>; }
+   else if constexpr (DIM == 2) { return internal::quadrature_interpolator::Derivatives2D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D, NBZ>; }
+   else if constexpr (DIM == 3) { return internal::quadrature_interpolator::Derivatives3D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D>; }
+   MFEM_ABORT("");
 }
 
 template<int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int VDIM, int D1D,
@@ -603,10 +603,10 @@ template<int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int VDIM, int D1D,
 QuadratureInterpolator::CollocatedGradKernelType
 QuadratureInterpolator::CollocatedGradKernels::Kernel()
 {
-   if (DIM == 1) { return internal::quadrature_interpolator::CollocatedDerivatives1D<Q_LAYOUT, GRAD_PHYS>; }
-   else if (DIM == 2) { return internal::quadrature_interpolator::CollocatedDerivatives2D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, NBZ>; }
-   else if (DIM == 3) { return internal::quadrature_interpolator::CollocatedDerivatives3D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D>; }
-   else { MFEM_ABORT(""); }
+   if constexpr (DIM == 1) { return internal::quadrature_interpolator::CollocatedDerivatives1D<Q_LAYOUT, GRAD_PHYS>; }
+   else if constexpr (DIM == 2) { return internal::quadrature_interpolator::CollocatedDerivatives2D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, NBZ>; }
+   else if constexpr (DIM == 3) { return internal::quadrature_interpolator::CollocatedDerivatives3D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D>; }
+   MFEM_ABORT("");
 }
 
 /// @endcond

fem/quadinterpolator.cpp

@@ -752,10 +752,10 @@ template <int DIM, int VDIM, int ND, int NQ>
 EvalKernel QuadratureInterpolator::EvalKernels::Kernel()
 {
    using namespace internal::quadrature_interpolator;
-   if (DIM == 1) { return Eval1D; }
-   else if (DIM == 2) { return Eval2D<VDIM,ND,NQ>; }
-   else if (DIM == 3) { return Eval3D<VDIM,ND,NQ>; }
-   else { MFEM_ABORT(""); }
+   if constexpr (DIM == 1) { return Eval1D; }
+   else if constexpr (DIM == 2) { return Eval2D<VDIM,ND,NQ>; }
+   else if constexpr (DIM == 3) { return Eval3D<VDIM,ND,NQ>; }
+   MFEM_ABORT("");
 }
 
 template <int DIM>

fem/tmop.cpp

@@ -4102,8 +4102,11 @@ void TMOP_Integrator::GetSurfaceFittingErrors(const Vector &d_loc,
 #ifdef MFEM_USE_MPI
       // Don't count the overlapping DOFs in parallel.
       // The pfes might be ordered byVDIM, while the loop goes consecutively.
-      const int dof_i = pfes->DofToVDof(i, 0);
-      if (parallel && pfes->GetLocalTDofNumber(dof_i) < 0) { continue; }
+      if (parallel)
+      {
+         const int dof_i = pfes->DofToVDof(i, 0);
+         if (pfes->GetLocalTDofNumber(dof_i) < 0) { continue; }
+      }
 #endif
 
       dof_cnt++;

mesh/mesh.cpp

@@ -15751,9 +15751,18 @@ Mesh PartitionMPI(int dim, int mpi_cnt, int elem_per_mpi, bool print,
 {
    MFEM_VERIFY(dim > 1, "Not implemented for 1D meshes.");
 
-   auto factor = [&](int N)
+   // Closest int divisor to the cubit root, going down.
+   auto factor3 = [](int N)
    {
-      for (int i = static_cast<int>(sqrt(N)); i > 0; i--)
+      for (int i = static_cast<int>(round(cbrt(N))); i > 0; i--)
+      {  if (N % i == 0) { return i; } }
+      return 1;
+   };
+
+   // Closest int divisor to the square root, going down.
+   auto factor2 = [](int N)
+   {
+      for (int i = static_cast<int>(round(sqrt(N))); i > 0; i--)
       { if (N % i == 0) { return i; } }
       return 1;
    };
@@ -15777,22 +15786,22 @@ Mesh PartitionMPI(int dim, int mpi_cnt, int elem_per_mpi, bool print,
    int el0_x, el0_y, el0_z;
    if (dim == 2)
    {
-      mpi_x = factor(mpi_cnt);
+      mpi_x = factor2(mpi_cnt);
       mpi_y = mpi_cnt / mpi_x;
 
       // Switch order for better balance.
-      el0_y = factor(el0);
+      el0_y = factor2(el0);
       el0_x = el0 / el0_y;
    }
    else
    {
-      mpi_x = factor(mpi_cnt);
-      mpi_y = factor(mpi_cnt / mpi_x);
+      mpi_x = factor3(mpi_cnt);
+      mpi_y = factor2(mpi_cnt / mpi_x);
       mpi_z = mpi_cnt / mpi_x / mpi_y;
 
       // Switch order for better balance.
-      el0_z = factor(el0);
-      el0_y = factor(el0 / el0_z);
+      el0_z = factor3(el0);
+      el0_y = factor2(el0 / el0_z);
       el0_x = el0 / el0_y / el0_z;
    }
 

miniapps/nurbs/nurbs_ex10p.cpp

@@ -264,7 +264,7 @@ int main(int argc, char *argv[])
    //    parallel degrees of freedom, with offsets given by array true_offset.
    FiniteElementCollection *fec = nullptr;
    NURBSExtension *NURBSext = nullptr;
-   if (mesh->NURBSext)
+   if (pmesh->NURBSext)
    {
       NURBSext = new NURBSExtension(pmesh->NURBSext, order);
       fec = new NURBSFECollection(order);

miniapps/tools/CMakeLists.txt

@@ -67,3 +67,16 @@ if (MFEM_USE_MPI)
    add_mfem_miniapp(nodal-transfer
      MAIN nodal-transfer.cpp LIBRARIES mfem)
 endif()
+
+if(MFEM_ENABLE_LIBBACKTRACE AND CMAKE_CXX_STANDARD GREATER_EQUAL 17)
+   find_package(LIBBACKTRACE REQUIRED)
+   
+   add_library(mprotect_trace SHARED mprotect_trace.cpp)
+   # remove the preceding "lib"
+   set_target_properties(mprotect_trace PROPERTIES PREFIX "")
+   target_link_libraries(mprotect_trace PRIVATE
+      LIBBACKTRACE::LIBBACKTRACE
+   )
+   install(TARGETS mprotect_trace
+   DESTINATION miniapps)
+endif()

miniapps/tools/gridfunction-bounds.cpp

@@ -23,6 +23,8 @@
 // (2) Dzanic et al., "A method for bounding high-order finite element
 //     functions: Applications to mesh validity and bounds-preserving limiters".
 //
+// We also use a recursive subdivision strategy to compute tighter estimate of
+// the function extremum.
 //
 // Compile with: make gridfunction-bounds
 //
@@ -31,9 +33,6 @@
 //  mpirun -np 4 gridfunction-bounds -nb 100 -ref 5 -bt 2 -l2
 
 #include "mfem.hpp"
-#include <memory>
-#include <iostream>
-#include <fstream>
 
 using namespace mfem;
 using namespace std;
@@ -56,6 +55,8 @@ int main (int argc, char *argv[])
    int  b_type           = -1;
    bool continuous       = true;
    int  nbrute           = 0;
+   int  rec_depth        = 4;
+   real_t rel_tol        = 1e-4;
 
    // Parse command-line options.
    OptionsParser args(argc, argv);
@@ -83,6 +84,12 @@ int main (int argc, char *argv[])
    args.AddOption(&nbrute, "-nb", "--nbrute",
                   "Brute force search for minimum in an array of nxnxn points "
                   "in each element.");
+   args.AddOption(&rec_depth, "-rd", "--rec-depth",
+                  "Maximum depth for recursive subdivision to compute function "
+                  "extremum.");
+   args.AddOption(&rel_tol, "-rt", "--rel-tol",
+                  "Relative tolerance for termination of recursive "
+                  "subdivision.");
    args.ParseCheck();
 
    Mesh mesh(mesh_file, 1, 1, false);
@@ -151,7 +158,19 @@ int main (int argc, char *argv[])
    ParGridFunction lowerb(&fes_pc), upperb(&fes_pc);
 
    // Compute bounds
-   pfunc_proj->GetElementBounds(lowerb, upperb, ref);
+   PLBound plb = pfunc_proj->GetElementBounds(lowerb, upperb, ref);
+
+   // Compute minimum and maximum bounds via recursion
+   Vector bound_rec_min(vdim), bound_rec_max(vdim);
+   for (int d = 0; d < vdim; d++)
+   {
+      auto min_interval = pfunc_proj->EstimateFunctionMinimum(d, plb, rec_depth,
+                                                              rel_tol);
+      auto max_interval = pfunc_proj->EstimateFunctionMaximum(d, plb, rec_depth,
+                                                              rel_tol);
+      bound_rec_min(d) = min_interval.first;
+      bound_rec_max(d) = max_interval.second;
+   }
 
    Vector bound_min(vdim), bound_max(vdim);
    for (int d = 0; d < vdim; d++)
@@ -236,17 +255,31 @@ int main (int argc, char *argv[])
       {
          for (int d = 0; d < vdim; d++)
          {
-            cout << "Brute force and bounding comparison for component " <<
+            cout << "Compare function extremum for component " <<
                  d << endl;
-            cout << "Brute force minimum and minimum bound: " << global_min(d)
-                 << " " <<  bound_min(d) << endl;
-
-            cout << "Brute force maximum and maximum bound: " << global_max(d)
-                 << " " <<  bound_max(d) << endl;
-
-            cout << "The difference in bounds is: " <<
-                 global_min(d)-bound_min(d) << " " <<
-                 bound_max(d)-global_max(d) << endl;
+            constexpr int w = 20;
+            cout << left << setw(w) << " "
+                 << setw(w)  << "Brute force"
+                 << setw(w) << "PL Bound"
+                 << setw(w) << "PL Bound + recursion" << endl
+                 << left << setw(w) << "Minimum: "
+                 << setw(w) << global_min(d)
+                 << setw(w) << bound_min(d)
+                 << setw(w) << bound_rec_min(d) << endl
+                 << left << setw(w) << "Difference: "
+                 << setw(w) << "-"
+                 << setw(w) << global_min(d)-bound_min(d)
+                 << setw(w) << global_min(d)-bound_rec_min(d) << endl;
+            cout << endl
+                 << left << setw(w) << "Maximum: "
+                 << setw(w) << global_max(d)
+                 << setw(w) << bound_max(d)
+                 << setw(w) << bound_rec_max(d) << endl
+                 << left << setw(w) << "Difference: "
+                 << setw(w) << "-"
+                 << setw(w) << bound_max(d)-global_max(d)
+                 << setw(w) << bound_rec_max(d)-global_max(d) << endl;
+            cout << endl;
          }
       }
    }
@@ -255,10 +288,19 @@ int main (int argc, char *argv[])
    {
       for (int d = 0; d < vdim; d++)
       {
-         cout << "Minimum bound for component " << d << " is " <<
-              bound_min(d) << endl;
-         cout << "Maximum bound for component " << d << " is " <<
-              bound_max(d) << endl;
+         cout << "Compare function extremum for component " <<
+              d << endl;
+         constexpr int w = 20;
+         cout << left << setw(w) << " "
+              << setw(w) << "PL Bound"
+              << setw(w) << "PL Bound + recursion" << endl
+              << left << setw(w) << "Minimum: "
+              << setw(w) << bound_min(d)
+              << setw(w) << bound_rec_min(d) << endl;
+         cout << endl
+              << left << setw(w) << "Maximum: "
+              << setw(w) << bound_max(d)
+              << setw(w) << bound_rec_max(d) << endl;
       }
    }
 

miniapps/tools/mprotect_trace.cpp

@@ -0,0 +1,266 @@
+/**
+This library helps in finding the sources of memory validity errors especially
+with the "debug" device backend. It works by preloading a modified mprotect
+which tracks every address that mprotect is called, and stores the stack
+trace of that call. When a SIGSEGV with code
+
+This should be built as a shared library with position independent code.
+
+To use, build this as a shared library "mprotect_trace.so" and set LD_PRELOAD=/path/to/mprotect_trace.so
+If it is in the same directory as your present directory, you may need to use LD_PRELOAD=./mprotect_trace.so
+
+Requires libbacktrace which can be obtained from
+https://github.com/ianlancetaylor/libbacktrace
+or
+https://github.com/gcc-mirror/gcc/tree/master/libbacktrace
+or from a package manager like spack.
+ */
+
+#include <dlfcn.h>
+#include <execinfo.h>
+#include <signal.h>
+#include <unistd.h>
+#include <sys/mman.h>
+#include <map>
+#include <mutex>
+#include <iostream>
+#include <backtrace.h>
+#include <backtrace-supported.h>
+#include <cstring>
+#include <vector>
+#include <cxxabi.h>
+#include <sstream>
+#include <chrono>
+
+#if BACKTRACE_SUPPORTED != 1
+#error "Backtrace not supported! See output file backtrace-supported.h for details."
+#endif
+
+/// Try to undo the name mangling. Return mangled name if unsuccesful.
+std::string demangle(const char* mangled)
+{
+   int status = 0;
+   char* demangled = abi::__cxa_demangle(mangled, nullptr, nullptr, &status);
+   std::string result = (status == 0 && demangled) ? demangled : mangled;
+   free(demangled);
+   return result;
+}
+
+/// Fallback filename to library if source wasn't available
+std::string fallback_file_name(uintptr_t pc)
+{
+   Dl_info info;
+   if (dladdr((void*)pc, &info))
+   {
+      // use ostringstream because it's easier to convert to from hex.
+      return std::string (info.dli_fname ? info.dli_fname : "??");
+   }
+   else
+   {
+      return "??";
+   }
+}
+
+/// Fallback function name from dladdr when nicer name from backtrace fails.
+std::string fallback_func_info(uintptr_t pc)
+{
+   Dl_info info;
+   std::ostringstream ss;
+   if (dladdr((void*)pc, &info))
+   {
+      if (info.dli_sname)
+      {
+         ss << " in " << abi::__cxa_demangle(info.dli_sname, nullptr, nullptr, nullptr);
+      }
+      else
+      {
+         ss << "pc " << (void*)pc;
+      }
+   }
+   else
+   {
+      ss << "pc " << (void*)pc;
+   }
+   return ss.str();
+}
+
+/// Tracks information about mprotect.
+struct ProtectionInfo
+{
+   /// Address that mprotect was called on.
+   void* addr;
+   /// The length of the protected buffer.
+   size_t len;
+   /// The protection code used in mprotect.
+   int prot;
+   /// The stack trace.
+   std::vector<std::string> trace;
+   /// The timespace of the mprotect call, used to differentiate calls to overlapping memory spans.
+   std::chrono::steady_clock::time_point timestamp;
+};
+
+static std::mutex g_mutex;
+
+/// the keys of this map are the addresses that mprotect was called on.s
+static std::map<void*, ProtectionInfo> g_protected;
+
+void error_callback(void*, const char* msg, int errnum)
+{
+   std::cerr << "libbacktrace error: " << msg << " (" << errnum << ")\n";
+}
+
+/// Get the string representation of the trace.
+int mprotect_backtrace_full_callback(void* data, uintptr_t pc,
+                                     const char* filename, int lineno, const char* function)
+{
+   char buf[512];
+   std::string demangled_name;
+   std::string fname_str;
+   if (!function)
+   {
+      demangled_name = fallback_func_info(pc);
+   }
+   else
+   {
+      demangled_name = demangle( function);
+   }
+   if (!filename)
+   {
+      fname_str = fallback_file_name(pc);
+   }
+   else
+   {
+      fname_str = filename;
+   }
+
+   snprintf(buf, sizeof(buf), "%s:%d in %s",fname_str.c_str(), lineno,
+            demangled_name.c_str());
+   ((std::vector<std::string>*)data)->emplace_back(buf);
+   return 0;
+}
+
+/// This gets called if libbacktrace ran into a problem.
+void mprotect_backtrace_error_callback(void* data, const char* msg, int)
+{
+   ((std::vector<std::string>*)data)->emplace_back(std::string("???: ") + msg);
+}
+
+/// Call this to try to generate a trace.
+void collect_trace(std::vector<std::string>& out)
+{
+   static backtrace_state* state = backtrace_create_state(nullptr, 1,
+                                                          error_callback, nullptr);
+   backtrace_full(state, 0,
+                  mprotect_backtrace_full_callback,
+                  mprotect_backtrace_error_callback,
+                  &out);
+}
+
+/// Our custom mprotect
+extern "C" int mprotect(void* addr, size_t len, int prot)
+{
+   static auto real_mprotect = (int(*)(void*, size_t, int))dlsym(RTLD_NEXT,
+                                                                 "mprotect");
+
+   std::vector<std::string> stack;
+   collect_trace(stack);
+   auto timepoint = std::chrono::steady_clock::now();
+
+   {
+      std::lock_guard<std::mutex> lock(g_mutex);
+      g_protected[addr] = { addr, len, prot, stack, timepoint};
+   }
+
+   return real_mprotect(addr, len, prot);
+}
+
+/// Return human readable protection code
+/// There are some codes I didn't write a conversion for.
+std::string prot_to_string(int prot)
+{
+   std::string prot_string;
+   if (prot == PROT_NONE) { prot_string += " PROT_NONE "; }
+   if (prot & PROT_READ) { prot_string += " PROT_READ "; }
+   if (prot & PROT_WRITE) { prot_string += " PROT_WRITE "; }
+   if (prot & PROT_EXEC) { prot_string += " PROT_EXEC "; }
+   return prot_string;
+}
+
+void segv_handler(int /*sig*/, siginfo_t* info, void*)
+{
+   void* fault_addr = info->si_addr;
+   auto code = info->si_code;
+   std::cerr << "Caught SIGSEGV with code " << code <<  " at address " <<
+             fault_addr << ".\n";
+   if (code == SEGV_ACCERR)
+   {
+      std::cerr <<
+                "This is a permission violation which was likely caused by accessing memory protected by mprotect.\n";
+   }
+   else
+   {
+      std::cerr <<
+                "This is not a permission violation which means the mprotect trace is likely not useful.\n";
+   }
+
+   std::vector<std::string> access_trace;
+   collect_trace(access_trace);
+
+   std::cerr << "\nACCESS STACK TRACE (this caused the violation):\n";
+   for (const auto& s : access_trace) { std::cerr << "  " << s << "\n"; }
+   {
+      std::lock_guard<std::mutex> lock(g_mutex);
+      // loop through addresses that m_protect was called on and see
+      // if the offending address fell in any of their ranges.
+      // It's possible that this has multiple hits if data was freed
+      // and another mprotect was called in a shifted location. This is dealt
+      // with here using timestamps from a high resolution clock. The most
+      // recent call is used.
+      static std::map<std::chrono::steady_clock::time_point, ProtectionInfo>
+      faulty_candidates;
+      static std::map<std::chrono::steady_clock::time_point, bool>
+      multiple_candidates;
+      for (const auto& [base, pi] : g_protected)
+      {
+         if (fault_addr >= base && fault_addr < (char*)base + pi.len)
+         {
+            multiple_candidates[pi.timestamp] = faulty_candidates.count(pi.timestamp) > 0;
+            faulty_candidates[pi.timestamp] = pi;
+         }
+      }
+      if (!faulty_candidates.empty())
+      {
+         // reverse the map from newest to oldest time. Extract most recent time.
+         const auto most_recent_time = faulty_candidates.rbegin()->first;
+         const auto &pi = faulty_candidates[most_recent_time];
+         const auto multiple_matches_found = multiple_candidates[most_recent_time];
+         std::cerr <<
+                   "\nPROTECTION STACK TRACE (this altered the access permissions most recently):\n";
+         std::cerr << "mprotect set the following permissions: " << prot_to_string(
+                      pi.prot) << "\n";
+         std::cerr << "Full protection number: " << pi.prot << "\n";
+         for (const auto& s : pi.trace) { std::cerr << "  " << s << "\n"; }
+         if (multiple_matches_found)
+         {
+            std::cerr <<
+                      "WARNING : Multiple protection stack traces found. Only showing one, but it may not be right. Maybe this tool needs a finer resolution clock.\n";
+         }
+      }
+      else
+      {
+         std::cerr << "No offending mprotect stack trace could be found\n";
+      }
+   }
+
+   exit(1);
+}
+
+// This code gets run before main.
+__attribute__((constructor))
+void install_handler()
+{
+   struct sigaction sa {};
+   sa.sa_flags = SA_SIGINFO;
+   sa.sa_sigaction = segv_handler;
+   sigaction(SIGSEGV, &sa, nullptr);
+}

tests/unit/fem/test_col_lag_der.cpp

@@ -151,6 +151,8 @@ TEST_CASE("Collocated Derivative Kernels", "[QuadratureInterpolator]")
       auto L = GENERATE(QVectorLayout::byNODES, QVectorLayout::byVDIM);
       auto P = GENERATE(true, false);
 
+      CAPTURE(L, P);
+
       const int nd = maps.ndof;
       const int nq = maps.nqpt;
 
@@ -163,7 +165,15 @@ TEST_CASE("Collocated Derivative Kernels", "[QuadratureInterpolator]")
       CGK::Run(dim, L, P, vdim, nd, nelem, maps.G.Read(), geom->J.Read(),
                evec_values.Read(), col_der.Write(), sdim, vdim, nd);
 
+      const real_t max_norm = qp_der.Normlinf();
+
       qp_der -= col_der;
-      REQUIRE(qp_der.Normlinf() == MFEM_Approx(0.0, 1e-10, 1e-10));
+
+      const real_t abs_err = qp_der.Normlinf();
+      const real_t rel_err = max_norm > 0_r ?
+                             abs_err/max_norm :
+                             abs_err > 0_r ? mfem::infinity() : 0_r;
+      CAPTURE(rel_err, max_norm);
+      CHECK(rel_err <= 1e-13);
    }
 }