Merge branch 'master' into mprotect_trace_mmcheck

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

CMakeLists.txt

@@ -633,7 +633,7 @@ set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
     NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
     ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CALIPER CODIPACK
     BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
-    ALGOIM ENZYME CUDA::cudart)
+    ALGOIM ENZYME CUDA::cudart LIBBACKTRACE)
 
 # Add all created targets and *_FOUND libraries in the variables TPL_TARGETS and
 # TPL_LIBRARIES, respectively.
@@ -657,7 +657,7 @@ list(REMOVE_DUPLICATES TPL_LIBRARIES)
 list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_INCLUDE_DIRS)
 # message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")
-
+ 
 message(STATUS "MFEM shared library: BUILD_SHARED_LIBS = ${BUILD_SHARED_LIBS}")
 message(STATUS "MFEM build type: CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 message(STATUS "MFEM version: v${MFEM_VERSION_STRING}")
@@ -724,7 +724,7 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})
 # Declaring the library
 mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 target_compile_features(mfem PUBLIC cxx_std_${CMAKE_CXX_STANDARD})
-# message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
+message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
 target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
 if (TPL_TARGETS)
   add_dependencies(mfem ${TPL_TARGETS})

config/cmake/modules/FindLIBBACKTRACE.cmake

@@ -0,0 +1,85 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+#[=======================================================================[.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.
+#]=======================================================================]
+
+include(MfemCmakeUtilities)
+include(FindPackageHandleStandardArgs)
+
+# https://github.com/ianlancetaylor/libbacktrace
+
+if(${CMAKE_HOST_SYSTEM_NAME} MATCHES "Darwin") ################################
+  # Find libbacktrace.
+  # Defines the following variables:
+  #   - LIBBACKTRACE_FOUND
+  #   - LIBBACKTRACE_LIBRARIES    (if needed)
+  #   - LIBBACKTRACE_INCLUDE_DIRS (if needed)
+  message(STATUS "[🟠 LIBBACKTRACE 🟠] Looking for LIBBACKTRACE ...")
+
+  set(BACKTRACE_DIR "$ENV{HOME}/usr/local/backtrace" CACHE PATH "")
+
+  # set(LIBBACKTRACE_SKIP_LOOKING_MSG TRUE)
+  mfem_find_package(LIBBACKTRACE LIBBACKTRACE BACKTRACE_DIR 
+                    "include" backtrace.h
+                    "lib" backtrace
+                    "Paths to headers required by LIBBACKTRACE." 
+                    "Libraries required by LIBBACKTRACE.")
+
+elseif(${CMAKE_HOST_SYSTEM_NAME} MATCHES "Linux") #############################
+  find_package(LIBBACKTRACE QUIET NO_MODULE)
+  # 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()
+  find_package_handle_standard_args(LIBBACKTRACE DEFAULT_MSG
+      LIBBACKTRACE_FOUND
+      LIBBACKTRACE_LIBRARIES
+      LIBBACKTRACE_INCLUDE_DIRS)
+else() 
+  message(STATUS "[🟠 LIBBACKTRACE 🟠] Unsupported platform!")
+endif()
+
+message(STATUS "[🟠 LIBBACKTRACE 🟠] LIBBACKTRACE_FOUND: ${LIBBACKTRACE_FOUND}")
+message(STATUS "[🟠 LIBBACKTRACE 🟠] LIBBACKTRACE_INCLUDE_DIRS: ${LIBBACKTRACE_INCLUDE_DIRS}")
+message(STATUS "[🟠 LIBBACKTRACE 🟠] LIBBACKTRACE_LIBRARIES: ${LIBBACKTRACE_LIBRARIES}")
+message(STATUS "[🟠 LIBBACKTRACE 🟠] LIBBACKTRACE_FIND_VERSION: ${LIBBACKTRACE_FIND_VERSION}")
+message(STATUS "[🟠 LIBBACKTRACE 🟠] LIBBACKTRACE_FIND_COMPONENTS: ${LIBBACKTRACE_FIND_COMPONENTS}")
+
+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")

examples/ex1.cpp

@@ -68,12 +68,25 @@
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>
+#include <unistd.h>
 
 using namespace std;
 using namespace mfem;
 
+#include <miniapps/tools/mem_check.hpp>
+
 int main(int argc, char *argv[])
 {
+   // Initialize the memory manager checker
+   mfem::MemoryManagerCheck::Init(argv[0]);
+
+   pid_t pid = getpid();
+   uid_t uid = getuid();
+   mfem::out << "\x1b[33m"
+             << "pid: " << pid
+             << ", uid: " << uid << "\x1b[m" <<
+             std::endl;
+
    // 1. Parse command-line options.
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
@@ -128,14 +141,14 @@ int main(int argc, char *argv[])
    //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
    //    largest number that gives a final mesh with no more than 50,000
    //    elements.
-   {
+   /*{
       int ref_levels =
          (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh.UniformRefinement();
       }
-   }
+   }*/
 
    // 5. Define a finite element space on the mesh. Here we use continuous
    //    Lagrange finite elements of the specified order. If order < 1, we

examples/mprotect.d

@@ -0,0 +1,9 @@
+#!/usr/sbin/dtrace -s
+
+#pragma D option quiet
+
+pid$target:libsystem_kernel.dylib:mprotect:entry
+{
+    printf("mprotect called: pid=%d, addr=%p, len=%d, prot=%d\n",
+           pid, arg0, arg1, arg2);
+}

examples/test_mprotect.cpp

@@ -0,0 +1,32 @@
+#include <sys/mman.h>
+#include <iostream>
+#include <unistd.h>
+
+int main()
+{
+   // Allocate memory
+   void *addr = mmap(nullptr, 4096, PROT_READ | PROT_WRITE,
+                     MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
+   if (addr == MAP_FAILED)
+   {
+      std::cerr << "mmap failed\n";
+      return 1;
+   }
+
+   std::cout << "mprotect called, PID: " << getpid() << std::endl;
+   sleep(10);
+
+   // Call mprotect
+   if (mprotect(addr, 4096, PROT_READ) == -1)
+   {
+      std::cerr << "mprotect failed\n";
+      return 1;
+   }
+
+   // Keep the process running to allow DTrace to attach
+   sleep(30);
+
+   // Clean up
+   munmap(addr, 4096);
+   return 0;
+}

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++;

general/nvtx.hpp

@@ -0,0 +1 @@
+../../stash/debug/nvtx.hpp

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,31 @@ 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)
+  message(STATUS "[🟠 LIBBACKTRACE 🟠] Building mprotect_trace miniapp")
+  message(STATUS "[🟠 LIBBACKTRACE 🟠] PROJECT_BINARY_DIR: ${PROJECT_BINARY_DIR}")
+  
+   find_package(LIBBACKTRACE REQUIRED)
+   
+   # MPROTECT_TRACE library
+   add_library(mprotect_trace SHARED mprotect_trace.cpp)
+   # remove the preceding "lib"
+   set_target_properties(mprotect_trace PROPERTIES PREFIX "")
+   target_include_directories(mprotect_trace PRIVATE ${PROJECT_BINARY_DIR}/config)
+   target_link_libraries(mprotect_trace PRIVATE LIBBACKTRACE::LIBBACKTRACE mfem)
+   install(TARGETS mprotect_trace DESTINATION miniapps)
+
+   # MEM CHECK library
+   add_library(mem_check SHARED mem_check.cpp)
+   set_target_properties(mem_check PROPERTIES PREFIX "")
+  #  target_include_directories(mem_check PRIVATE /opt/homebrew/opt/fmt/include)
+   target_include_directories(mem_check PRIVATE ${PROJECT_BINARY_DIR}/config)
+   target_link_libraries(mem_check PRIVATE mfem LIBBACKTRACE::LIBBACKTRACE)
+  #  target_link_libraries(mem_check PRIVATE -L/opt/homebrew/opt/fmt/lib -lfmt)
+   install(TARGETS mem_check DESTINATION miniapps)
+else()
+  message(STATUS "[🟠 LIBBACKTRACE 🟠] Skipping mprotect_trace miniapp")
+  message(STATUS "[🟠 LIBBACKTRACE 🟠] Set MFEM_ENABLE_LIBBACKTRACE=ON to build mprotect_trace")
+  message(STATUS "[🟠 LIBBACKTRACE 🟠] Set CMAKE_CXX_STANDARD=17 or higher to build mprotect_trace")
+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/mem_check.cpp

@@ -0,0 +1,727 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+///////////////////////////////////////////////////////////////////////////////
+// Environment variables:
+//   - MM_TRACE: trace functions: malloc, free, calloc, realloc, memalign
+//   - MM_DEBUG: print debug information
+//   - MM_CHECK: check memory manager
+//   - MM_ALL: print full backtrace
+//   - MM_ARGS: print function arguments
+//   - MM_ORG: switch to org mode, print '*' instead of ' '
+
+///////////////////////////////////////////////////////////////////////////////
+// On macOS, we need to set the following flags to allow dynamic lookup,
+// 'dsymutil ex1' is also required to generate the debug symbols
+// - -g -O1 -fno-inline-functions -fno-omit-frame-pointer
+// - -Wl,-undefined,dynamic_lookup
+
+///////////////////////////////////////////////////////////////////////////////
+// For example 1:
+// - add '#include <miniapps/tools/mem_check.hpp>'
+// - add 'mfem::MemoryManagerCheck::Init(argv[0]);'
+// - make ex1 && dsymutil ex1
+// - tput reset && MM_TRACE=1 MM_DEBUG=1 DYLD_INSERT_LIBRARIES=../miniapps/tools/mem_check.dylib ./ex1
+
+#include <cassert>
+#include <cstdlib>
+#include <cxxabi.h>
+#include <dlfcn.h>
+#include <iostream>
+#include <cstring>
+
+#include <backtrace.h>
+#include <backtrace-supported.h>
+
+#if BACKTRACE_SUPPORTED != 1
+#error "Backtrace not supported! See output file backtrace-supported.h for details."
+#endif
+
+#include "mem_check.hpp"
+
+#include "config/config.hpp" // IWYU pragma: keep
+#include "general/mem_manager.hpp"
+
+// #undef NVTX_COLOR
+// #define NVTX_COLOR nvtx::kLawnGreen
+// #include "general/nvtx.hpp"
+#define dbg(...)
+
+///////////////////////////////////////////////////////////////////////////////
+// static helper functions
+static const char *strrnchr(const char *s, const unsigned char c, int n)
+{
+   assert(s);
+   size_t len = strlen(s);
+   char *p = (char *)s + len - 1;
+   for (; n; n--, p--, len--)
+   {
+      for (; len; p--, len--) { if (*p == c) { break; } }
+      if (!len) { return nullptr; }
+      if (n == 1) { return p; }
+   }
+   return nullptr;
+}
+
+static const char *getFilename(const char *filename, const unsigned char delim,
+                               const int n)
+{
+   assert(filename);
+   const char *f = strrnchr(filename, delim, n);
+   return f ? f + 1 : filename;
+}
+
+static void err_callback(void *data, const char *msg, int errnum)
+{
+   dbg("error: {} errnum: #{}", msg, errnum);
+   assert(false);
+}
+
+static const char *cxx_demangle(const char *mangled_name)
+{
+   int status;
+   assert(mangled_name);
+   const char *demangled_name =
+      abi::__cxa_demangle(mangled_name, nullptr, nullptr, &status);
+   const bool succeeded = status == 0;
+   const bool memory_allocation_failure_occurred = status == -1;
+   const bool one_argument_is_invalid = status == -3;
+   assert(not one_argument_is_invalid);
+   if (memory_allocation_failure_occurred)
+   {
+      printf("[demangle] memory_allocation_failure_occurred!");
+      fflush(nullptr);
+      assert(false);
+   }
+   return (succeeded) ? demangled_name : mangled_name;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+class backtrace_data
+{
+public:
+   backtrace_state *state = nullptr;
+   char *function = nullptr;
+   const char *filename = nullptr;
+   int lineno = -1;
+   uintptr_t address = 0x0;
+   bool dump = false;
+   bool main = false;
+   bool mm = false;
+   bool mfem = false;
+   bool got = false;
+   int depth = 0;
+   static constexpr int STACK_MAX = 32768;
+   char stack[STACK_MAX];
+   bool skip = false;
+
+   backtrace_data(backtrace_state *state): state(state) {}
+
+   ~backtrace_data() { if (function) { free((void*) function); } }
+
+   void ini(const bool dmp = false)
+   {
+      if (function) { free((void*) function); }
+      function = nullptr;
+      lineno = -1;
+      address = 0x0;
+      dump = dmp;
+      main = false;
+      mm = false;
+      mfem = false;
+      got = false;
+      depth = 0;
+      stack[0] = 0;
+      skip = false;
+   }
+
+   /// Returns 0 to continue tracing
+   inline int continue_tracing() { return main ? 1 : 0; }
+
+   void update(const char *demangled, uintptr_t PC,
+               const char *path_name = nullptr, const int line = -1)
+   {
+      assert(demangled);
+      // assert(path_name);
+      // update the context depth
+      depth += 1;
+      dbg("{}", demangled);
+      if (path_name) { dbg("{}:{}", path_name, line); }
+      // check if we have reached the 'main' function
+      main = !strncmp(demangled, "main", 4);
+      dbg("main:\x1b[33m {}", main);
+      // MFEM namespace test
+      const bool mfem_namespace = strncmp(demangled, "mfem::", 6) == 0;
+      dbg("mfem_namespace:\x1b[33m {}", mfem_namespace);
+      mfem |= mfem_namespace;
+      // Test if we are doing an allocation from inside the mfem::mem_manager
+      // If it is the case, set 'skip' to true
+      if (path_name && !skip)
+      {
+         const char *file = getFilename(path_name, '/', 2);
+         dbg("file:\x1b[33m {}", file);
+         const bool mm_hpp_file = strncmp(file, "general/mem_manager.hpp", 13) == 0;
+         const bool mm_cpp_file = strncmp(file, "general/mem_manager.cpp", 13) == 0;
+         // We want to filter out the allocations done in the mem_manager
+         mm |= mm_cpp_file or mm_hpp_file;
+         dbg("mem_manager:\x1b[33m {}", mm);
+         // Skip where the 'New' and 'Delete' are done
+         // Not sure these are the ones to skip!! ❌❌
+         const size_t len = std::strlen(demangled);
+         const char *new_suffix = ">::New(int, mfem::MemoryType)";
+         const char *del_suffix = ">::Delete(int, mfem::MemoryType)";
+         const size_t new_suffix_len = strlen(new_suffix);
+         const size_t del_suffix_len = strlen(del_suffix);
+         const bool New = (new_suffix_len <= len) &&
+                          strncmp(demangled + len - new_suffix_len, new_suffix,
+                                  new_suffix_len) == 0;
+         const bool Delete = (del_suffix_len <= len) &&
+                             strncmp(demangled + len - del_suffix_len, del_suffix,
+                                     del_suffix_len) == 0;
+         skip |= New or Delete;
+         dbg("skip:\x1b[33m {}", skip);
+      }
+      // Record the stack
+      static char path_file[PATH_MAX];
+      const int n_char_printed =
+         snprintf(path_file, PATH_MAX, "\n\t%s %s:%d", demangled,
+                  path_name ? path_name : "???", line);
+      assert(n_char_printed < STACK_MAX);
+      dbg("path_file:\x1b[33m {}", path_file);
+      strcat(stack, path_file);
+      // if we already have caught the function name, return
+      if (function != nullptr)
+      {
+         assert(filename); // we should have recorded the filename
+         // assert(lineno >= 0); // and set the line number
+         dbg("function:\x1b[33m {}, returning", function);
+         return;
+      }
+      function = strdup(demangled);
+      // assert(path_name);
+      filename = path_name ? strdup(path_name) : strdup("???");
+      address = PC, lineno = line;
+      dbg("\x1b[31m[{}] {}:{}", filename, function, line);
+   }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+class backrace_check
+{
+   backtrace_data *data {};
+   backtrace_state *state {};
+   // number of frames to skip:
+   // '1' will go down to 'mm_alloc'
+   // '3' should be sufficient
+   static constexpr int SKIP = 3;
+
+public:
+   backrace_check()  = default;
+
+   ~backrace_check() { delete data; }
+
+   void ini(const char *argv0)
+   {
+      assert(argv0);
+      state = backtrace_create_state(argv0,
+                                     BACKTRACE_SUPPORTS_THREADS,
+                                     err_callback, nullptr);
+      data = new backtrace_data(state);
+   }
+
+   int backtrace(const bool dump = false)
+   {
+      if (!state or !data) { return EXIT_FAILURE; } // not ready
+      data->ini(dump); // flush the data
+      // int code =
+      backtrace_simple(state, SKIP, simple_callback, err_callback, data);
+      // printf("\x1b[33m[%d]\x1b[m",code);
+      return EXIT_SUCCESS;
+   }
+
+private:
+   static void sym_callback(void *data, uintptr_t pc, const char *symname,
+                            uintptr_t symval, uintptr_t symsize)
+   {
+      if (!symname) { return; }
+      auto *ctx = static_cast<backtrace_data *>(data);
+      const char *demangled = cxx_demangle(symname);
+      ctx->update(demangled, pc);
+   }
+
+   static int filter(const bool debug, const char *demangled)
+   {
+      if (!debug) { return 0; }
+      assert(demangled);
+      printf("\nFiltering OUT '%s'!\n", demangled);
+      return 0;
+   }
+
+   static int full_callback(void *data, uintptr_t pc, const char *filename,
+                            int lineno, const char *function)
+   {
+      dbg("filename: {}, lineno: {}, function: {}",
+          filename ? filename : "?",
+          lineno,
+          function ? function : "?");
+      static const bool all = getenv("MM_ALL");
+      auto *ctx = static_cast<backtrace_data *>(data);
+      const bool debug = ctx->dump or all;
+
+      if (!function) // symbol hit
+      {
+         if (debug)
+         {
+            printf("\n[full_callback:symbol] filename:%s, lineno=%d, pc=0x%lx",
+                   filename ? filename : "???", lineno, pc);
+         }
+         return backtrace_syminfo(ctx->state, pc, sym_callback, err_callback,
+                                  data);
+      }
+
+      const char *demangled = cxx_demangle(function);
+      assert(demangled);
+
+      if (debug)
+      {
+         printf("\n\t\x1b[32m[full_callback:function] %s\x1b[m", demangled);
+      }
+
+      // Filtering
+      if (strncmp("std::", demangled, 5) == 0) { return filter(all, demangled); }
+      if (strncmp("__gnu_cxx::", demangled, 11) == 0) { return filter(all, demangled); }
+
+      // Debug if ALL
+      if (debug) { printf("\n\t"); }
+
+      // Update context
+      assert(filename);
+      ctx->update(demangled, pc, filename, lineno);
+
+      if (debug) { printf("%s:%d %s", filename, lineno, demangled); }
+      return EXIT_SUCCESS;
+   }
+
+   static int simple_callback(void *data, uintptr_t pc)
+   {
+      auto *ctx = static_cast<backtrace_data *>(data);
+      backtrace_pcinfo(ctx->state, pc, full_callback, err_callback, data);
+      return ctx->continue_tracing(); // returns 0 to continue tracing
+   }
+
+public:
+   bool mm() { return data->mm; }
+   bool mfem() { return data->mfem; }
+   bool skip() { return data->skip; }
+   int depth() { return data->depth; }
+   uintptr_t address() { return data->address; }
+   const char* function() { return data->function; }
+   const char* filename() { return data->filename; }
+   int lineno() { return data->lineno; }
+   char *stack() { return data->stack; }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+using mm_stack_t = std::unordered_map<const void *, const char *>;
+static mm_stack_t *mm_stack_map = nullptr;
+
+static backrace_check *bt_check = nullptr;
+static bool trace = false, hooked = false, dlsymd = false;
+
+namespace mfem
+{
+
+///////////////////////////////////////////////////////////////////////////////
+// MemoryManagerCheck
+MemoryManagerCheck::MemoryManagerCheck(const char *argv0)
+{
+   hooked = false;
+   dbg();
+   bt_check = new backrace_check();
+   mm_stack_map = new mm_stack_t();
+   bt_check->ini(argv0);
+   hooked = true;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+MemoryManagerCheck::~MemoryManagerCheck()
+{
+   hooked = false;
+   dbg();
+   delete bt_check, bt_check = nullptr;
+   delete mm_stack_map, mm_stack_map = nullptr;
+   hooked = true;
+}
+
+} // namespace mfem
+
+///////////////////////////////////////////////////////////////////////////////
+static void mmCheck(const void *ptr, const bool is_new, const bool dump)
+{
+   static const bool all = getenv("MM_ALL") != nullptr;
+   static const bool debug = getenv("MM_DEBUG") != nullptr;
+   static const bool mfem_debug = getenv("MFEM_DEBUG") != nullptr;
+   static const bool args = getenv("MM_ARGS") != nullptr;
+   static const bool org_mode = getenv("MM_ORG") != nullptr;
+   static const bool mm_check = getenv("MM_CHECK") != nullptr;
+
+   const auto is_del = not is_new;
+   const auto tab = org_mode ? "*" : "  ";
+
+   // now backtracing if ready
+   if (!bt_check) { return; }
+
+   if (bt_check->backtrace(dump) != EXIT_SUCCESS) { return; }
+
+   // we have the stack, checking the backtrace
+   dbg();
+
+   const bool from_memory_manager = bt_check->mm();     // memory manager
+   dbg("from_memory_manager:\x1b[35m {}", from_memory_manager);
+   const bool mfem_namespace = bt_check->mfem();    // mfem namespace
+   dbg("mfem_namespace:\x1b[35m {}", mfem_namespace);
+   const bool skip = bt_check->skip();
+   dbg("skip:\x1b[35m {}", skip);
+   const int depth = bt_check->depth();
+   const int frames = depth - (org_mode ? -1 : 1);
+   const auto address = bt_check->address();
+   const char *function = bt_check->function() ? bt_check->function() : "???";
+   const char *filename = bt_check->filename() ? bt_check->filename() : "???";
+   const int lineno = bt_check->lineno();
+   // dbg("\x1b[35m[{}] {}:{}", filename, function, lineno);
+
+   const std::string demangled_function(function);
+   const size_t first_parenthesis = demangled_function.find_first_of('(');
+   const std::string no_args_demangled_function =
+      demangled_function.substr(0, first_parenthesis);
+   const std::string display_function =
+      args ? demangled_function : no_args_demangled_function;
+   dbg("display_function:\x1b[35m {}", display_function);
+
+   const size_t first_3A = display_function.find_first_of(':');
+   const size_t first_3C = display_function.find_first_of('<');
+   const size_t first_5B = display_function.find_first_of('[');
+   assert(first_3A <= (first_5B < 0) ? first_3A : first_5B);
+   const size_t first_3AC = ((first_3A ^ first_3C) < 0)
+                            ? std::max(first_3A, first_3C)
+                            : std::min(first_3A, first_3C);
+   const std::string root = (first_3A != first_3C)
+                            ? display_function.substr(0, first_3AC)
+                            : display_function;
+   dbg("root:\x1b[35m {}", root);
+   const int color = address % (256 - 46) + 46;
+
+   if (debug && !mfem_debug)
+   {
+      if (all) { std::cout << std::endl; }
+      // Generating tabs
+      for (int k = 0; k < frames; ++k) { std::cout << tab; }
+      // Bold outputing
+      if (!org_mode) { std::cout << "\x1b[38;5;" << color << ";1m"; }
+      else { std::cout << " "; }
+      mfem::out << "[" << (filename ? filename : "")
+                << ":" << lineno
+                << ":" << display_function
+                << "]\x1b[m";
+   }
+
+   if (skip and debug) { printf(" skip!"); }
+   if (skip or not mm_check) { return; }
+
+   assert(ptr);
+
+   // should return false when the memory manager does not exist,
+   // or is not configured and also when then maps are not valid.
+   const bool known = mfem::mm.IsKnown((void *)ptr);
+
+   dbg("known:\x1b[35m {}", known);
+
+   if (debug) { printf(" %sMFEM", mfem_namespace ? "" : "Not "); }
+
+   static auto error = [](const void *ptr)
+   {
+      assert(bt_check), assert(mm_stack_map);
+      printf("\nStack:%s\n", bt_check->stack());
+      printf("\nFirst:%s\n", mm_stack_map->at(ptr));
+      fflush(nullptr);
+      assert(false);
+   };
+
+   // all this logic needs to be adapted to the new memory manager
+   if (mfem_namespace)
+   {
+      dbg("in mfem namespace");
+      if (from_memory_manager)
+      {
+         dbg("in mfem namespace, from memory manager");
+         if (known and is_new)
+         {
+            printf("\nTrying to 'insert' (%p), known by the MM!", ptr);
+            return error(ptr);
+         }
+         else if (not known and is_del) // when device is not configured ?!
+         {
+            printf("\nTrying to 'erase' (%p), not known by the MM!", ptr);
+            return error(ptr);
+         }
+         else
+         {
+            if (debug) { printf(", known: ok"); }
+         }
+      }
+      else // not from memory manager
+      {
+         dbg("in mfem namespace, not from memory manager");
+         if (debug) { printf(", !MM"); }
+         if (known and is_new)
+         {
+            printf("\nTrying to 'new' (%p), known by the MM!", ptr);
+            return error(ptr);
+         }
+         else if (known and is_del)
+         {
+            printf("\nTrying to 'delete' (%p), known by the MM!", ptr);
+            return error(ptr);
+         }
+         else if (not known)
+         {
+            // Should we allow user to do allocations ?!
+            printf("\nTrying to new/del (%p), not known by the MM!", ptr);
+            return error(ptr);
+         }
+         else
+         {
+            if (debug) { printf(" unknown: ok/error ?"); }
+            return error(ptr); // 🔥🔥🔥 return error for now
+         }
+      }
+   }
+   else // not mfem namespace
+   {
+      // from memory manager but not mfem namespace: possible ? 🔥
+      if (from_memory_manager)
+      {
+         if (known and is_new)
+         {
+            printf("\nTrying to 'insert' (%p), known by the MM!", ptr);
+            return error(ptr);
+         }
+         else if (not known and is_del)
+         {
+            printf("Trying to 'erase' (%p), not known by the MM!", ptr);
+            return error(ptr);
+         }
+         else
+         {
+            if (debug) { printf(" known: ok"); }
+         }
+      }
+      else // not from memory manager, not mfem namespace
+      {
+         if (known and is_new)
+         {
+            printf("\nTrying to 'new' (%p), known by the MM!", ptr);
+            return error(ptr);
+         }
+         if (known and is_del)
+         {
+            printf("Trying to 'delete' (%p), known by the MM!", ptr);
+            return error(ptr);
+         }
+         if (debug) { printf(", unknown: leave it"); } // 🔥 to the user
+      }
+   }
+
+   static auto mmAdd = [](const void *ptr,
+                          const bool is_new,
+                          const char *stack)
+   {
+      assert(ptr);
+      assert(stack);
+      const bool known = mm_stack_map->find(ptr) != mm_stack_map->end();
+      if (is_new and not known)
+      {
+         dbg("Add to stack map");
+         mm_stack_map->emplace(ptr, strdup(stack));
+      }
+      if (known and !is_new) // 🔥 to check ?!
+      {
+         dbg("Remove from stack map");
+         mm_stack_map->erase(ptr);
+      }
+   };
+
+   mmAdd(ptr, is_new, bt_check->stack());
+   fflush(nullptr);
+   // assert(false);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Linux: LD_PRELOAD + dlsym(RTLD_NEXT, ...)
+// MacOS: DYLD_INSERT_LIBRARIES + DYLD_INTERPOSE
+// Apple macOS specifics to interpose functions, instead of dlsym NEXT
+#ifdef __APPLE__
+#define DYLD_INTERPOSE(_to,_from) \
+   __attribute__((used)) static struct{ const void* to; const void* from; } \
+   _interpose_##_from __attribute__ ((section ("__DATA,__interpose"))) = \
+   { (const void*)(unsigned long)&_to, (const void*)(unsigned long)&_from };
+#else
+#define DYLD_INTERPOSE(...)
+#endif
+
+///////////////////////////////////////////////////////////////////////////////
+// LD_PRELOADed functions: malloc, free, calloc, realloc and memalign
+using free_t = void (void *);
+using malloc_t = void *(size_t);
+using calloc_t = void *(size_t, size_t);
+using realloc_t = void *(void *, size_t);
+#ifdef __APPLE__
+using memalign_t = int (void **memptr, size_t alignment, size_t size);
+#else
+using memalign_t = void *(size_t, size_t);
+#endif
+using mm_t = std::unordered_map<void *, size_t>;
+
+static free_t *_free = nullptr;
+static malloc_t *_malloc = nullptr;
+static calloc_t *_calloc = nullptr;
+static realloc_t *_realloc = nullptr;
+static memalign_t *_memalign = nullptr;
+
+///////////////////////////////////////////////////////////////////////////////
+// With APPLE, don't dlsym the functions, as it will use the interposed one
+static void _init()
+{
+   if (getenv("MM_TRACE")) { trace = true; }
+#ifdef __APPLE__ // use the system ones
+   _free = (free_t *)free, assert(_free);
+   _malloc = (malloc_t *)malloc, assert(_malloc);
+   _calloc = (calloc_t *)calloc, assert(_calloc);
+   _realloc = (realloc_t *)realloc, assert(_realloc);
+   _memalign = (memalign_t *)posix_memalign, assert(_memalign);
+#else
+   _free = (free_t *)dlsym(RTLD_NEXT, "free"), assert(_free);
+   _malloc = (malloc_t *)dlsym(RTLD_NEXT, "malloc"), assert(_malloc);
+   _calloc = (calloc_t *)dlsym(RTLD_NEXT, "calloc"), assert(_calloc);
+   _realloc = (realloc_t *)dlsym(RTLD_NEXT, "realloc"), assert(_realloc);
+   _memalign = (memalign_t *)dlsym(RTLD_NEXT, "memalign"), assert(_memalign);
+#endif
+   assert(_free and _malloc and _calloc and _realloc and _memalign);
+   hooked = true, dlsymd = true;
+}
+
+// ALLOC //////////////////////////////////////////////////////////////////////
+#ifdef __APPLE__
+void *mm_malloc(size_t size) // Red
+#else
+void *malloc(size_t size) // Red
+#endif
+{
+   if (!_malloc) { _init(); }
+   if (!hooked) { return _malloc(size); }
+   hooked = false;
+   void *ptr = _malloc(size);
+   assert(ptr);
+   if (trace) { printf("\n\x1b[31m[malloc] %p (%ld)\x1b[m", ptr, size); }
+   mmCheck(ptr, true, false); // new, dont show full stack
+   hooked = true;
+   return ptr;
+}
+DYLD_INTERPOSE(mm_malloc, malloc)
+
+// FREE ///////////////////////////////////////////////////////////////////////
+void mm_free(void *ptr) // Green
+{
+   if (!_free) { _init(); }
+   if (!hooked) { return _free(ptr); }
+   if (!ptr) { return; }
+   hooked = false;
+   if (trace) { printf("\n\x1b[32m[free] %p\x1b[m", ptr); }
+   mmCheck(ptr, false, false); // delete, dont show full stack
+   _free(ptr);
+   hooked = true;
+}
+DYLD_INTERPOSE(mm_free, free)
+
+// CALLOC /////////////////////////////////////////////////////////////////////
+void *mm_calloc(size_t nmemb, size_t size) // Yellow
+{
+   if (not dlsymd) // if we are not yet dlsym'ed, just do it ourselves
+   {
+      static const size_t MEM_MAX = 8192;
+      static char mem[MEM_MAX];
+      static size_t m = 0;
+      const size_t bytes = nmemb * size;
+      void *ptr = &mem[m];
+      m += bytes;
+      assert(m < MEM_MAX);
+      for (size_t k = 0; k < bytes; k += 1) { *(((char *)ptr) + k) = 0; }
+      return ptr;
+   }
+   if (!hooked) { return _calloc(nmemb, size); }
+   hooked = false;
+   void *ptr = _calloc(nmemb, size);
+   if (trace) { printf("\n\x1b[33m[calloc] %p (%ld)\x1b[m", ptr, size); }
+   mmCheck(ptr, true, false); // new, dont show full stack
+   hooked = true;
+   return ptr;
+}
+DYLD_INTERPOSE(mm_calloc, calloc)
+
+// REALLOC ////////////////////////////////////////////////////////////////////
+void *mm_realloc(void *ptr, size_t size) // Blue
+{
+   if (!_realloc) { _init(); }
+   if (!hooked) { return _realloc(ptr, size); }
+   hooked = false;
+   void *nptr = _realloc(ptr, size);
+   assert(nptr);
+   if (trace) { printf("\n\x1b[34;7m[realloc] %p(%ld)\x1b[m", nptr, size); }
+   mmCheck(nptr, true, false); // new, dont show full stack
+   hooked = true;
+   return nptr;
+}
+DYLD_INTERPOSE(mm_realloc, realloc)
+
+// MEMALIGN ///////////////////////////////////////////////////////////////////
+#ifdef __APPLE__
+int mm_memalign(void **memptr, size_t alignment, size_t size) // Magenta
+#else
+void *memalign(size_t alignment, size_t size) // Magenta
+#endif
+{
+   if (!_memalign) { _init(); }
+   if (!hooked)
+   {
+#ifdef __APPLE__
+      return _memalign(memptr, alignment, size);
+#else
+      return _memalign(alignment, size);
+#endif
+   }
+   hooked = false;
+#ifdef __APPLE__
+   const int rtn = _memalign(memptr, alignment, size);
+   void *ptr = *memptr;
+#else
+   void *ptr = _memalign(alignment, size);
+#endif
+   assert(ptr);
+   if (trace) { printf("\n\x1b[35;7m[memalign] %p(%ld)\x1b[m", ptr, size); }
+   mmCheck(ptr, true, false); // new, dont show full stack
+   hooked = true;
+#ifdef __APPLE__
+   return rtn;
+#else
+   return ptr;
+#endif
+}
+DYLD_INTERPOSE(mm_memalign, posix_memalign)

miniapps/tools/mem_check.hpp

@@ -0,0 +1,35 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+#pragma once
+
+namespace mfem
+{
+
+// The Memory Manager checker singleton class
+class MemoryManagerCheck
+{
+   MemoryManagerCheck(const char *argv0);
+   ~MemoryManagerCheck();
+
+   static MemoryManagerCheck& Singleton(const char *arg0)
+   {
+      static MemoryManagerCheck instance(arg0);
+      return instance;
+   }
+
+public:
+   static void Init(const char *arg0) { Singleton(arg0); }
+
+   MemoryManagerCheck(const MemoryManagerCheck&) = delete;
+   MemoryManagerCheck& operator=(const MemoryManagerCheck&) = delete;
+};
+
+} // namespace mfem

miniapps/tools/mprotect_trace.cpp

@@ -0,0 +1,338 @@
+/**
+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 <cassert>
+#include <execinfo.h>
+#include <csignal>
+#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>
+#include <thread>
+#include <atomic>
+
+#include "general/mem_manager.hpp"
+
+#if BACKTRACE_SUPPORTED != 1
+#error "Backtrace not supported! See output file backtrace-supported.h for details."
+#endif
+
+#ifdef __APPLE__
+#define DYLD_SYMBOL(name) name
+#define DYLD_INTERPOSE(_to,_from) \
+   __attribute__((used)) static struct{ const void* to; const void* from; } \
+   _interpose_##_from __attribute__ ((section ("__DATA,__interpose"))) = \
+   { (const void*)(unsigned long)&_to, (const void*)(unsigned long)&_from };
+#else
+#define DYLD_INTERPOSE(...)
+#define DYLD_SYMBOL(name) dlsym(RTLD_NEXT, #name);
+#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.
+// static std::map<void*, ProtectionInfo> g_protected;
+static std::map<void*, void*> 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[PATH_MAX];
+   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);
+}
+
+using mprotect_t = int (void*, size_t, int);
+
+int computeValue()
+{
+   // Simulate some work
+   std::this_thread::sleep_for(std::chrono::milliseconds(100));
+   return 42;
+}
+
+std::atomic<int> result{0}; // Atomic variable to store the result
+
+/// Our custom mprotect
+/// 🔥🔥🔥 WIP, as on macos and apple hardware, mprotect cannot be used anymore.
+/// 🔥🔥🔥 Trying the dtrace approach...
+#ifdef __APPLE__
+extern "C"
+int mm_mprotect(void* addr, size_t len, int prot)
+#else
+extern "C" int mprotect(void* addr, size_t len, int prot)
+#endif
+{
+   static mprotect_t *_mprotect = (mprotect_t *) mprotect;
+   assert(true); // ✅
+   // static std::map<void*, int> map; // ✅
+
+   assert(getuid() == 501);
+   /*{
+      return _mprotect(addr, len, prot);
+   }*/
+
+
+   assert(false);
+   const bool known = mfem::mm.IsKnown(addr);
+   assert(false);
+   // assert(known);
+
+   {
+      // auto worker = []() { /*computeValue();*/ };  // ❌
+      // Launch the thread
+      // std::thread t(worker);//, std::ref(result)); // ❌
+      // Join the thread to ensure the result is set
+      // t.join();
+   }
+
+   {
+      // std::lock_guard<std::mutex> lock(g_mutex); // ✅
+      // printf("\n\x1b.\x1b[m"); // ❌
+      // map[addr] = prot; // ❌
+   }
+   // const int rtn = _mprotect(addr, len, prot);
+   // assert(rtn == 0);
+   // printf("\n\x1b[31m[mprotect] %p (%ld, %x)\x1b[m", addr, len, prot);
+
+   // std::vector<std::string> stack;
+   // collect_trace(stack); // ❌
+   // auto timepoint = std::chrono::steady_clock::now(); // ✅
+   // (void) timepoint;
+   {
+      // std::lock_guard<std::mutex> lock(g_mutex);
+      // g_protected[addr] = addr;
+      /*{
+         addr,
+         // len, prot,
+         // stack,
+         // timepoint
+      };*/
+   }
+   return _mprotect(addr, len, prot);
+   // return rtn;
+   // return EXIT_SUCCESS;
+}
+DYLD_INTERPOSE(mm_mprotect, mprotect)
+
+/// 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, 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);
    }
 }