Try to fix CI failures on Windows by including windows.h before psapi.h

* class Memory; some methods were still using int for sizes
* class Vector, for size and indexing
* all 1D "forall" macros and function templates
  - HIP cannot launch kernels with >= 2^32 total threads
  - CUDA seems to support kernel launches with >= 2^32 total threads
* class DeviceTensor; individual dimensions still use int, however, the
  total 1D index computation uses bigint
* element and face geometric factors use bigint sizes for memory allocations

Left some FIXME comments to be addressed later.

CHANGELOG

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

config/config.hpp

@@ -62,6 +62,11 @@ constexpr real_t operator""_r(unsigned long long v)
    return static_cast<real_t>(v);
 }
 
+// MFEM bigint type
+
+/// MFEM's "big" integer type.
+typedef long long int bigint;
+
 } // namespace mfem
 
 // Return value for main function in examples that should be skipped by testing

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_base.cpp

@@ -663,58 +663,59 @@ const
    #pragma omp critical (DofToQuad)
 #endif
    {
-      // If the new Dof2Quad is already present, e.g. added in a previous call
-      // or added by another omp thread, return.
+      // Do not run if the new Dof2Quad is already present, e.g. added in a
+      // previous call or added by another omp thread.
       if (DofToQuad::SearchArray(dof2quad_array, ir,
-                                 DofToQuad::LEXICOGRAPHIC_FULL))
-      { return; }
-
-      // Undo the native ordering which is what FiniteElement::GetDofToQuad
-      // returns.
-      auto *d2q_new = new DofToQuad(d2q);
-      d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
-      const int nqpt = ir.GetNPoints();
-
-      const int b_dim = (range_type == VECTOR) ? dim : 1;
-
-      for (int i = 0; i < nqpt; i++)
+                                 DofToQuad::LEXICOGRAPHIC_FULL) == nullptr)
       {
-         for (int d = 0; d < b_dim; d++)
+         // Undo the native ordering which is what FiniteElement::GetDofToQuad
+         // returns.
+         auto *d2q_new = new DofToQuad(d2q);
+         d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
+         const int nqpt = ir.GetNPoints();
+
+         const int b_dim = (range_type == VECTOR) ? dim : 1;
+
+         for (int i = 0; i < nqpt; i++)
          {
-            for (int j = 0; j < dof; j++)
+            for (int d = 0; d < b_dim; d++)
             {
-               const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
-               d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
-               d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
+                  d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
+                  d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               }
             }
          }
-      }
 
-      const int g_dim = [this]()
-      {
-         switch (deriv_type)
+         const int g_dim = [this]()
          {
-            case GRAD: return dim;
-            case DIV: return 1;
-            case CURL: return cdim;
-            default: return 0;
-         }
-      }();
-
-      for (int i = 0; i < nqpt; i++)
-      {
-         for (int d = 0; d < g_dim; d++)
-         {
-            for (int j = 0; j < dof; j++)
+            switch (deriv_type)
             {
-               const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
-               d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
-               d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               case GRAD: return dim;
+               case DIV: return 1;
+               case CURL: return cdim;
+               default: return 0;
+            }
+         }();
+
+         for (int i = 0; i < nqpt; i++)
+         {
+            for (int d = 0; d < g_dim; d++)
+            {
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
+                  d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
+                  d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               }
             }
          }
+
+         dof2quad_array.Append(d2q_new);
       }
 
-      dof2quad_array.Append(d2q_new);
    }
 }
 

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/fespace.cpp

@@ -282,14 +282,7 @@ int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs_) const
 void FiniteElementSpace::AdjustVDofs(Array<int> &vdofs)
 {
    int n = vdofs.Size(), *vdof = vdofs;
-   for (int i = 0; i < n; i++)
-   {
-      int j;
-      if ((j = vdof[i]) < 0)
-      {
-         vdof[i] = -1-j;
-      }
-   }
+   for (int i = 0; i < n; i++) { vdof[i] = UnsignIndex(vdof[i]); }
 }
 
 void FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs,
@@ -483,13 +476,14 @@ void FiniteElementSpace::ReorderElementToDofTable()
    for (int k = 0, dof_counter = 0; k < nnz; k++)
    {
       const int sdof = J[k]; // signed dof
-      const int dof = (sdof < 0) ? -1-sdof : sdof;
+      const int dof = UnsignIndex(sdof);
       int new_dof = dof_marker[dof];
       if (new_dof < 0)
       {
          dof_marker[dof] = new_dof = dof_counter++;
       }
-      J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
+      // Preserve the sign of sdof
+      J[k] = (sdof < 0) ? FlipIndexSign(new_dof) : new_dof;
    }
 }
 
@@ -547,7 +541,7 @@ void MarkDofs(const Array<int> &dofs, Array<int> &mark_array)
 {
    for (auto d : dofs)
    {
-      mark_array[d >= 0 ? d : -1 - d] = -1;
+      mark_array[UnsignIndex(d)] = -1;
    }
 }
 
@@ -931,7 +925,7 @@ void FiniteElementSpace::AddDependencies(
             if (std::abs(coef) > 1e-12)
             {
                const int mdof = master_dofs[j];
-               if (mdof != sdof && mdof != (-1-sdof))
+               if (mdof != sdof && mdof != FlipIndexSign(sdof))
                {
                   deps.Add(sdof, mdof, coef);
                }
@@ -1024,7 +1018,7 @@ int FiniteElementSpace::GetDegenerateFaceDofs(int index, Array<int> &dofs,
    // FiniteElementSpace::AddDependencies.
 
    Array<int> edof;
-   int order = GetEdgeDofs(-1 - index, edof, variant);
+   int order = GetEdgeDofs(FlipIndexSign(index), edof, variant);
 
    int nv = fec->DofForGeometry(Geometry::POINT);
    int ne = fec->DofForGeometry(Geometry::SEGMENT);
@@ -1710,8 +1704,8 @@ SparseMatrix *FiniteElementSpace::RefinementMatrix_main(
 
          for (int i = 0; i < fine_ldof; i++)
          {
-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);
 
             if (!mark[m])
             {
@@ -1772,7 +1766,7 @@ SparseMatrix *FiniteElementSpace::VariableOrderRefinementMatrix(
          for (int i = 0; i < fine_ldof; i++)
          {
             const int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int m = UnsignIndex(r);
 
             if (!mark[m])
             {
@@ -2482,8 +2476,8 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
          {
             if (!std::isfinite(lR(i, 0))) { continue; }
 
-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);
 
             if (is_dg || !mark[m])
             {
@@ -3201,7 +3195,7 @@ void FiniteElementSpace::CalcEdgeFaceVarOrders(
             else
             {
                // degenerate face (i.e., edge-face constraint)
-               slave_orders |= edge_orders[-1 - slave.index];
+               slave_orders |= edge_orders[FlipIndexSign(slave.index)];
             }
          }
 

fem/fespace.hpp

@@ -1150,7 +1150,7 @@ public:
 
    /// Helper to return the DOF associated with a sign encoded DOF
    static inline int DecodeDof(int dof)
-   { return (dof >= 0) ? dof : (-1 - dof); }
+   { return UnsignIndex(dof); }
 
    /// Helper to determine the DOF and sign of a sign encoded DOF
    static inline int DecodeDof(int dof, real_t& sign)

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/lor/lor.cpp

@@ -158,15 +158,16 @@ void LORBase::ConstructLocalDofPermutation(Array<int> &perm_) const
                int i;
                i = dofmap_lor[off_lor + i1 + i2*2];
                int s1 = i < 0 ? -1 : 1;
-               int idof_lor = vdof_lor[absdof(i)];
+               int idof_lor = vdof_lor[UnsignIndex(i)];
                i = dofmap_ho[off_ho + i1*n1 + i2*n2];
                int s2 = i < 0 ? -1 : 1;
-               int idof_ho = vdof_ho[absdof(i)];
+               int idof_ho = vdof_ho[UnsignIndex(i)];
                int s3 = idof_lor < 0 ? -1 : 1;
                int s4 = idof_ho < 0 ? -1 : 1;
                int s = s1*s2*s3*s4;
-               i = absdof(idof_ho);
-               perm_[absdof(idof_lor)] = s < 0 ? -1-absdof(i) : absdof(i);
+               i = UnsignIndex(idof_ho);
+               perm_[UnsignIndex(idof_lor)] = s < 0 ? -1-UnsignIndex(i) :
+                                              UnsignIndex(i);
             }
          }
       };
@@ -232,7 +233,7 @@ void LORBase::ConstructDofPermutation() const
          int j = l_perm[i];
          int s = j < 0 ? -1 : 1;
          int t_i = pfes_lor->GetLocalTDofNumber(i);
-         int t_j = pfes_ho->GetLocalTDofNumber(absdof(j));
+         int t_j = pfes_ho->GetLocalTDofNumber(UnsignIndex(j));
          // Either t_i and t_j both -1, or both non-negative
          if ((t_i < 0 && t_j >=0) || (t_j < 0 && t_i >= 0))
          {

fem/lor/lor.hpp

@@ -57,8 +57,6 @@ private:
    /// values (after temporarily changing them for LOR assembly).
    void ResetIntegrationRules(GetIntegratorsFn get_integrators);
 
-   static inline int absdof(int i) { return i < 0 ? -1-i : i; }
-
 protected:
    enum FESpaceType { H1, ND, RT, L2, INVALID };
 

fem/pfespace.cpp

@@ -424,7 +424,7 @@ void ParFiniteElementSpace::GetGroupComm(
             {
                if (ind[l] < 0)
                {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                   if (g_ldof_sign)
                   {
                      (*g_ldof_sign)[dofs[l]] = -1;
@@ -462,7 +462,7 @@ void ParFiniteElementSpace::GetGroupComm(
             {
                if (ind[l] < 0)
                {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                   if (g_ldof_sign)
                   {
                      (*g_ldof_sign)[dofs[l]] = -1;
@@ -500,7 +500,7 @@ void ParFiniteElementSpace::GetGroupComm(
             {
                if (ind[l] < 0)
                {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                   if (g_ldof_sign)
                   {
                      (*g_ldof_sign)[dofs[l]] = -1;
@@ -538,16 +538,16 @@ void ParFiniteElementSpace::ApplyLDofSigns(Array<int> &dofs) const
    {
       if (dofs[i] < 0)
       {
-         if (ldof_sign[-1-dofs[i]] < 0)
+         if (ldof_sign[FlipIndexSign(dofs[i])] < 0)
          {
-            dofs[i] = -1-dofs[i];
+            dofs[i] = FlipIndexSign(dofs[i]);
          }
       }
       else
       {
          if (ldof_sign[dofs[i]] < 0)
          {
-            dofs[i] = -1-dofs[i];
+            dofs[i] = FlipIndexSign(dofs[i]);
          }
       }
    }
@@ -699,7 +699,8 @@ void ParFiniteElementSpace::GetSharedEdgeDofs(
       for (int i = 0; i < dofs.Size(); i++)
       {
          const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = di >= 0 ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
       }
    }
 }
@@ -723,7 +724,8 @@ void ParFiniteElementSpace::GetSharedTriangleDofs(
       for (int i = 0; i < dofs.Size(); i++)
       {
          const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = di >= 0 ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
       }
    }
 }
@@ -747,7 +749,8 @@ void ParFiniteElementSpace::GetSharedQuadrilateralDofs(
       for (int i = 0; i < dofs.Size(); i++)
       {
          const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = (di >= 0) ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
       }
    }
 }
@@ -1487,7 +1490,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
          GetElementVDofs(my_elems[i], ldofs);
          for (int j = 0; j < ldofs.Size(); j++)
          {
-            int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
+            int ldof = UnsignIndex(ldofs[j]);
 
             if (ldof_marker[ldof] != fn)
             {
@@ -1548,7 +1551,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
          GetElementVDofs(my_elems[i], ldofs);
          for (int j = 0; j < ldofs.Size(); j++)
          {
-            int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
+            int ldof = UnsignIndex(ldofs[j]);
 
             if (ldof_marker[ldof] != fn)
             {
@@ -1573,14 +1576,15 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
 
       for (int i = 0; i < num_ldofs; i++)
       {
-         int ldof = (ldofs_fn[i] >= 0 ? ldofs_fn[i] : -1-ldofs_fn[i]);
+         int ldof = UnsignIndex(ldofs_fn[i]);
          ldof_marker[ldof] = i;
       }
 
       for ( ; j < j_end; j++)
       {
-         int ldof = (send_J[j] >= 0 ? send_J[j] : -1-send_J[j]);
-         send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] : -1-ldof_marker[ldof]);
+         const int ldof = UnsignIndex(send_J[j]);
+         send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] :
+                      FlipIndexSign(ldof_marker[ldof]));
       }
    }
 
@@ -1672,12 +1676,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
    {
       for (int j_end = face_nbr_ldof.GetI()[fn+1]; j < j_end; j++)
       {
-         int ldof = face_nbr_ldof.GetJ()[j];
-         if (ldof < 0)
-         {
-            ldof = -1-ldof;
-         }
-
+         const int ldof = UnsignIndex(face_nbr_ldof.GetJ()[j]);
          face_nbr_glob_dof_map[j] = dof_face_nbr_offsets[fn] + ldof;
       }
    }
@@ -1721,7 +1720,7 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
    MFEM_ASSERT(Nonconforming() && i >= pmesh->GetNumFaces(), "");
    int el1, el2, inf1, inf2;
    pmesh->GetFaceElements(i, &el1, &el2);
-   el2 = -1 - el2;
+   el2 = FlipIndexSign(el2);
    pmesh->GetFaceInfos(i, &inf1, &inf2);
    MFEM_ASSERT(0 <= el2 && el2 < face_nbr_element_dof.Size(), "");
    const int nd = face_nbr_element_dof.RowSize(el2);
@@ -1737,7 +1736,8 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
    for (int j = 0; j < vdofs.Size(); j++)
    {
       const int ldof = vdofs[j];
-      vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] : -1-vol_vdofs[-1-ldof];
+      vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] :
+                 FlipIndexSign(vol_vdofs[FlipIndexSign(ldof)]);
    }
 }
 
@@ -2061,8 +2061,8 @@ void ParFiniteElementSpace::GetGhostFaceDofs(const MeshId &face_id,
 
          for (int j = 0; j < ne; j++)
          {
-            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
-                             /*         */ : (-1 - (first + (-1 - ind[j])));
+            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j]) :
+                             FlipIndexSign(first + FlipIndexSign(ind[j]));
          }
       }
       else
@@ -2072,8 +2072,8 @@ void ParFiniteElementSpace::GetGhostFaceDofs(const MeshId &face_id,
          const int *ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[i]);
          for (int j = 0; j < ne; j++)
          {
-            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
-                             /*         */ : (-1 - (first + (-1 - ind[j])));
+            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j]) :
+                             FlipIndexSign(first + FlipIndexSign(ind[j]));
          }
       }
    }
@@ -2866,7 +2866,7 @@ void NeighborRowMessage::Encode(int rank)
 
             if (ind && (edof = ind[edof]) < 0)
             {
-               edof = -1 - edof;
+               edof = FlipIndexSign(edof);
                s = -1;
             }
 
@@ -3067,10 +3067,10 @@ void NeighborRowMessage::Decode(int rank)
 
          // If edof arrived with a negative index, flip it, and the scaling.
          real_t s = (edof < 0) ? -1.0 : 1.0;
-         edof = (edof < 0) ? -1 - edof : edof;
+         edof = UnsignIndex(edof);
          if (ind && (edof = ind[edof]) < 0)
          {
-            edof = -1 - edof;
+            edof = FlipIndexSign(edof);
             s *= -1.0;
          }
 
@@ -3121,10 +3121,10 @@ void NeighborRowMessage::Decode(int rank)
 
             // If edof arrived with a negative index, flip it, and the scaling.
             s = (edof < 0) ? -1.0 : 1.0;
-            edof = (edof < 0) ? -1 - edof : edof;
+            edof = UnsignIndex(edof);
             if (ind && (edof = ind[edof]) < 0)
             {
-               edof = -1 - edof;
+               edof = FlipIndexSign(edof);
                s *= -1.0;
             }
 
@@ -4405,12 +4405,9 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
          {
             for (int j = 0; j < dofs.Size(); j++)
             {
-               int row = DofToVDof(dofs[j], vd);
-               if (row < 0) { row = -1 - row; }
-
-               int col = DofToVDof(old_dofs[j], vd, old_ndofs);
-               if (col < 0) { col = -1 - col; }
-
+               const int row = UnsignIndex(DofToVDof(dofs[j], vd));
+               const int col = UnsignIndex(DofToVDof(old_dofs[j], vd,
+                                                     old_ndofs));
                i_diag[row] = col;
             }
          }
@@ -4435,9 +4432,7 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
       {
          for (int j = 0; j < dofs.Size(); j++)
          {
-            int row = DofToVDof(dofs[j], vd);
-            if (row < 0) { row = -1 - row; }
-
+            const int row = UnsignIndex(DofToVDof(dofs[j], vd));
             if (i_diag[row] == i_diag[row+1]) // diag row empty?
             {
                i_offd[row] = old_dofs[j + vd * dofs.Size()];
@@ -4546,9 +4541,9 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
    {
       const Embedding &emb = dtrans.embeddings[k];
 
-      int fine_rank = old_ranks[k];
-      int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
-                        : old_pncmesh->ElementRank(emb.parent);
+      const int fine_rank = old_ranks[k];
+      const int coarse_rank = (emb.parent < 0) ? FlipIndexSign(emb.parent)
+                              : old_pncmesh->ElementRank(emb.parent);
 
       if (coarse_rank != MyRank && fine_rank == MyRank)
       {
@@ -4636,8 +4631,8 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
             {
                if (!std::isfinite(lR(i, 0))) { continue; }
 
-               int r = DofToVDof(dofs[i], vd);
-               int m = (r >= 0) ? r : (-1 - r);
+               const int r = DofToVDof(dofs[i], vd);
+               const int m = UnsignIndex(r);
 
                if (is_dg || !mark[m])
                {
@@ -4686,8 +4681,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
             {
                if (!std::isfinite(lR(i, 0))) { continue; }
 
-               int r = DofToVDof(dofs[i], vd);
-               int m = (r >= 0) ? r : (-1 - r);
+               const int m = UnsignIndex(DofToVDof(dofs[i], vd));
 
                if (is_dg || !mark[m])
                {

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/qfunction.hpp

@@ -271,10 +271,7 @@ inline void QuadratureFunction::GetValues(
    const int s_offset = qspace->Offset(idx);
    const int sl_size = qspace->Offset(idx + 1) - s_offset;
    // Make the values matrix memory an alias of the quadrature function memory
-   Memory<real_t> &values_mem = values.GetMemory();
-   values_mem.Delete();
-   values_mem.MakeAlias(GetMemory(), vdim*s_offset, vdim*sl_size);
-   values.SetSize(vdim, sl_size);
+   values.MakeRef(GetMemory(), vdim*s_offset, vdim, sl_size);
 }
 
 inline void QuadratureFunction::GetValues(

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/restriction.cpp

@@ -844,8 +844,6 @@ void ConformingFaceRestriction::ComputeGatherIndices(
    gather_offsets[0] = 0;
 }
 
-static inline int absdof(int i) { return i < 0 ? -1-i : i; }
-
 void ConformingFaceRestriction::SetFaceDofsScatterIndices(
    const Mesh::FaceInformation &face,
    const int face_index,
@@ -868,9 +866,9 @@ void ConformingFaceRestriction::SetFaceDofsScatterIndices(
    {
       const int lex_volume_dof = face_map[face_dof];
       const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof]; // signed
-      const int volume_dof = absdof(s_volume_dof);
+      const int volume_dof = UnsignIndex(s_volume_dof);
       const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
-      const int global_dof = absdof(s_global_dof);
+      const int global_dof = UnsignIndex(s_global_dof);
       const int restriction_dof = face_dofs*face_index + face_dof;
       scatter_indices[restriction_dof] = s_global_dof;
       ++gather_offsets[global_dof + 1];
@@ -897,10 +895,10 @@ void ConformingFaceRestriction::SetFaceDofsGatherIndices(
    {
       const int lex_volume_dof = face_map[face_dof];
       const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof];
-      const int volume_dof = absdof(s_volume_dof);
+      const int volume_dof = UnsignIndex(s_volume_dof);
       const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
       const int sgn = (s_global_dof >= 0) ? 1 : -1;
-      const int global_dof = absdof(s_global_dof);
+      const int global_dof = UnsignIndex(s_global_dof);
       const int restriction_dof = face_dofs*face_index + face_dof;
       const int s_restriction_dof = (sgn >= 0) ? restriction_dof : -1 -
                                     restriction_dof;

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/array.cpp

@@ -23,7 +23,7 @@ namespace mfem
 template <class T>
 void Array<T>::Print(std::ostream &os, int width) const
 {
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       os << data[i];
       if ( !((i+1) % width) || i+1 == size )
@@ -44,7 +44,7 @@ void Array<T>::Save(std::ostream &os, int fmt) const
    {
       os << size << '\n';
    }
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       os << operator[](i) << '\n';
    }
@@ -55,11 +55,11 @@ void Array<T>::Load(std::istream &in, int fmt)
 {
    if (fmt == 0)
    {
-      int new_size;
+      bigint new_size;
       in >> new_size;
       SetSize(new_size);
    }
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       in >> operator[](i);
    }
@@ -71,7 +71,7 @@ T Array<T>::Max() const
    MFEM_ASSERT(size > 0, "Array is empty with size " << size);
 
    T max = operator[](0);
-   for (int i = 1; i < size; i++)
+   for (bigint i = 1; i < size; i++)
    {
       if (max < operator[](i))
       {
@@ -88,7 +88,7 @@ T Array<T>::Min() const
    MFEM_ASSERT(size > 0, "Array is empty with size " << size);
 
    T min = operator[](0);
-   for (int i = 1; i < size; i++)
+   for (bigint i = 1; i < size; i++)
    {
       if (operator[](i) < min)
       {
@@ -104,7 +104,7 @@ template <class T>
 void Array<T>::PartialSum()
 {
    T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       sum+=operator[](i);
       operator[](i) = sum;
@@ -116,9 +116,8 @@ void Array<T>::Abs()
 {
    static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
    const bool useDevice = UseDevice();
-   const int N = size;
    auto y = ReadWrite(useDevice);
-   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(useDevice, size, [=] MFEM_HOST_DEVICE (bigint i)
    {
       y[i] = std::abs(y[i]);
    });
@@ -129,7 +128,7 @@ template <class T>
 T Array<T>::Sum() const
 {
    T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       sum+=operator[](i);
    }
@@ -141,7 +140,7 @@ template <class T>
 int Array<T>::IsSorted() const
 {
    T val_prev = operator[](0), val;
-   for (int i = 1; i < size; i++)
+   for (bigint i = 1; i < size; i++)
    {
       val=operator[](i);
       if (val < val_prev)
@@ -159,7 +158,7 @@ bool Array<T>::IsConstant() const
 {
    if (size < 2) { return true; }
    const T v0 = data[0];
-   for (int i = 1; i < size; i++)
+   for (bigint i = 1; i < size; i++)
    {
       if (data[i] != v0)
       {

general/array.hpp

@@ -50,9 +50,9 @@ protected:
    /// Pointer to data
    Memory<T> data;
    /// Size of the array
-   int size;
+   bigint size;
 
-   inline void GrowSize(int minsize);
+   inline void GrowSize(bigint minsize);
 
    static_assert(std::is_trivial<T>::value, "type T must be trivial");
 
@@ -68,11 +68,11 @@ public:
    inline Array(MemoryType mt) : data(mt), size(0) { }
 
    /// Creates array of @a asize elements
-   explicit inline Array(int asize)
+   explicit inline Array(bigint asize)
       : size(asize) { if (asize > 0) { data.New(asize); } }
 
    /// Creates array of @a asize elements with a given MemoryType
-   inline Array(int asize, MemoryType mt)
+   inline Array(bigint asize, MemoryType mt)
       : data(mt), size(asize) { if (asize > 0) { data.New(asize, mt); } }
 
    /** @brief Creates array using an externally allocated host pointer @a data_
@@ -81,7 +81,7 @@ public:
 
        When @a own_data is true, the pointer @a data_ must be allocated with
        MemoryType given by MemoryManager::GetHostMemoryType(). */
-   inline Array(T *data_, int asize, bool own_data = false)
+   inline Array(T *data_, bigint asize, bool own_data = false)
    { data.Wrap(data_, asize, own_data); size = asize; }
 
    /// Copy constructor: deep copy from @a src
@@ -93,7 +93,7 @@ public:
    inline Array(const Array<CT> &src);
 
    /// Construct an Array from a C-style array of static length
-   template <typename CT, int N>
+   template <typename CT, bigint N>
    explicit inline Array(const CT (&values)[N]);
 
    /// Construct an Array from a braced initializer list of convertible type
@@ -114,10 +114,22 @@ public:
    Array<T> &operator=(const Array<T> &src) { src.Copy(*this); return *this; }
 
    /// Move assignment operator
+   /** If *this is a non-owning view (e.g., from MakeRef()), the data is copied
+       so that the base is also modified. */
    Array<T> &operator=(Array<T> &&src)
    {
       if (this == &src) { return *this; }
-      Swap(src);  // Swap does not use move assignment!
+      // If *this is a non-owning view (alias), and its capacity is sufficient
+      // to contain src, then copy into *this so that the alias's base memory is
+      // modified.
+      if (!OwnsData() && Capacity() >= src.Size())
+      {
+         *this = src; // Copy assignment.
+      }
+      else
+      {
+         Swap(src); // Swap the pointers only.
+      }
       src.DeleteAll();
       return *this;
    }
@@ -163,47 +175,47 @@ public:
    void MakeDataOwner() const { data.SetHostPtrOwner(true); }
 
    /// Return the logical size of the array.
-   inline int Size() const { return size; }
+   inline bigint Size() const { return size; }
 
    /// Change the logical size of the array, keep existing entries.
-   inline void SetSize(int nsize);
+   inline void SetSize(bigint nsize);
 
-   /// Same as SetSize(int) plus initialize new entries with 'initval'.
-   inline void SetSize(int nsize, const T &initval);
+   /// Same as SetSize(bigint) plus initialize new entries with 'initval'.
+   inline void SetSize(bigint nsize, const T &initval);
 
    /** @brief Resize the array to size @a nsize using MemoryType @a mt. Note
        that unlike the other versions of SetSize(), the current content of the
        array is not preserved. */
-   inline void SetSize(int nsize, MemoryType mt);
+   inline void SetSize(bigint nsize, MemoryType mt);
 
    /** Maximum number of entries the array can store without allocating more
        memory. */
-   inline int Capacity() const { return data.Capacity(); }
+   inline bigint Capacity() const { return data.Capacity(); }
 
    /// Ensures that the allocated size is at least the given size.
-   inline void Reserve(int capacity)
+   inline void Reserve(bigint capacity)
    { if (capacity > Capacity()) { GrowSize(capacity); } }
 
    /// Reference access to the ith element.
-   inline T & operator[](int i);
+   inline T & operator[](bigint i);
 
    /// Const reference access to the ith element.
-   inline const T &operator[](int i) const;
+   inline const T &operator[](bigint i) const;
 
    /// Append element 'el' to array, resize if necessary.
-   inline int Append(const T & el);
+   inline bigint Append(const T & el);
 
    /// STL-like push_back. Append element 'el' to array, resize if necessary.
    void push_back(const T &el) { Append(el); }
 
    /// Append another array to this array, resize if necessary.
-   inline int Append(const T *els, int nels);
+   inline bigint Append(const T *els, bigint nels);
 
    /// Append another array to this array, resize if necessary.
-   inline int Append(const Array<T> &els) { return Append(els, els.Size()); }
+   inline bigint Append(const Array<T> &els) { return Append(els, els.Size()); }
 
    /// Prepend an 'el' to the array, resize if necessary.
-   inline int Prepend(const T &el);
+   inline bigint Prepend(const T &el);
 
    /// Return the last element in the array.
    inline T &Last();
@@ -212,13 +224,13 @@ public:
    inline const T &Last() const;
 
    /// Append element when it is not yet in the array, return index.
-   inline int Union(const T & el);
+   inline bigint Union(const T & el);
 
    /// Return the first index where 'el' is found; return -1 if not found.
-   inline int Find(const T &el) const;
+   inline bigint Find(const T &el) const;
 
    /// Do bisection search for 'el' in a sorted array; return -1 if not found.
-   inline int FindSorted(const T &el) const;
+   inline bigint FindSorted(const T &el) const;
 
    /// Delete the last entry of the array.
    inline void DeleteLast() { if (size > 0) { size--; } }
@@ -241,16 +253,19 @@ public:
    /// Make this Array a reference to a pointer.
    /** When @a own_data is true, the pointer @a data_ must be allocated with
        MemoryType given by MemoryManager::GetHostMemoryType(). */
-   inline void MakeRef(T *data_, int size_, bool own_data = false);
+   inline void MakeRef(T *data_, bigint size_, bool own_data = false);
 
    /// Make this Array a reference to a pointer.
    /** When @a own_data is true, the pointer @a data_ must be allocated with
        MemoryType given by @a mt. */
-   inline void MakeRef(T *data_, int size, MemoryType mt, bool own_data);
+   inline void MakeRef(T *data_, bigint size, MemoryType mt, bool own_data);
 
    /// Make this Array a reference to 'master'.
    inline void MakeRef(const Array &master);
 
+   /// Make this Array a reference to the given sub-Memory of @a base.
+   inline void MakeRef(Memory<T> &base, bigint offset, bigint size_);
+
    /// Reset the Array to use the given external Memory @a mem and size @a s.
    /** If @a own_mem is false, the Array will not own any of the pointers of
        @a mem.
@@ -258,7 +273,7 @@ public:
        Note that when @a own_mem is true, the @a mem object can be destroyed
        immediately by the caller but `mem.Delete()` should NOT be called since
        the Array object takes ownership of all pointers owned by @a mem. */
-   inline void NewMemoryAndSize(const Memory<T> &mem, int s, bool own_mem);
+   inline void NewMemoryAndSize(const Memory<T> &mem, bigint s, bool own_mem);
 
    /**
     * @brief Permute the array using the provided indices. Sorts the indices
@@ -274,7 +289,7 @@ public:
    inline void Permute(const I &indices) { Permute(I(indices)); }
 
    /// Copy sub array starting from @a offset out to the provided @a sa.
-   inline void GetSubArray(int offset, int sa_size, Array<T> &sa) const;
+   inline void GetSubArray(bigint offset, bigint sa_size, Array<T> &sa) const;
 
    /// Prints array to stream with width elements per row.
    void Print(std::ostream &out = mfem::out, int width = 4) const;
@@ -297,7 +312,7 @@ public:
 
    /** @brief Set the Array size to @a new_size and read that many entries from
        the stream @a in. */
-   void Load(int new_size, std::istream &in)
+   void Load(bigint new_size, std::istream &in)
    { SetSize(new_size); Load(in, 1); }
 
    /** @brief Find the maximal element in the array, using the comparison
@@ -320,7 +335,7 @@ public:
    void Unique()
    {
       T* end = std::unique((T*)data, data + size);
-      SetSize((int)(end - data));
+      SetSize((bigint)(end - data));
    }
 
    /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
@@ -406,7 +421,7 @@ template <class T>
 inline bool operator==(const Array<T> &LHS, const Array<T> &RHS)
 {
    if ( LHS.Size() != RHS.Size() ) { return false; }
-   for (int i=0; i<LHS.Size(); i++)
+   for (bigint i=0; i<LHS.Size(); i++)
    {
       if ( LHS[i] != RHS[i] ) { return false; }
    }
@@ -436,13 +451,13 @@ public:
    Array2D() { M = N = 0; }
 
    /// Construct an m x n 2D array.
-   Array2D(int m, int n) : array1d(m*n) { M = m; N = n; }
+   Array2D(int m, int n) : array1d(bigint(m)*n) { M = m; N = n; }
 
    Array2D(const Array2D &) = default;
    Array2D(Array2D &&) = default;
 
    /// Set the 2D array size to m x n.
-   void SetSize(int m, int n) { array1d.SetSize(m*n); M = m; N = n; }
+   void SetSize(int m, int n) { array1d.SetSize(bigint(m)*n); M = m; N = n; }
 
    int NumRows() const { return M; }
    int NumCols() const { return N; }
@@ -547,16 +562,16 @@ public:
 
    /// Construct a 3D array of size n1 x n2 x n3.
    Array3D(int n1, int n2, int n3)
-      : array1d(n1*n2*n3) { N2 = n2; N3 = n3; }
+      : array1d(bigint(n1)*n2*n3) { N2 = n2; N3 = n3; }
 
    /// Set the 3D array size to n1 x n2 x n3.
    void SetSize(int n1, int n2, int n3)
-   { array1d.SetSize(n1*n2*n3); N2 = n2; N3 = n3; }
+   { array1d.SetSize(bigint(n1)*n2*n3); N2 = n2; N3 = n3; }
 
    /// Get the 3D array size in the first dimension.
    int GetSize1() const
    {
-      const int size = array1d.Size();
+      const bigint size = array1d.Size();
       return size == 0 ? 0 : size / (N2 * N3);
    }
 
@@ -764,7 +779,7 @@ inline Array<T>::Array(const Array<CT> &src)
    : size(src.Size())
 {
    size > 0 ? data.New(size) : data.Reset();
-   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
+   for (bigint i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
 }
 
 template <typename T>
@@ -775,7 +790,7 @@ inline Array<T>::Array(std::initializer_list<CT> values) : Array(values.size())
    std::copy(values.begin(), values.end(), begin());
 }
 
-template <typename T> template <typename CT, int N>
+template <typename T> template <typename CT, bigint N>
 inline Array<T>::Array(const CT (&values)[N]) : Array(N)
 {
    std::copy(values, values + N, begin());
@@ -789,9 +804,9 @@ inline void Array<T>::Swap(Array &other)
 }
 
 template <class T>
-inline void Array<T>::GrowSize(int minsize)
+inline void Array<T>::GrowSize(bigint minsize)
 {
-   const int nsize = std::max(minsize, 2 * data.Capacity());
+   const bigint nsize = std::max(minsize, 2 * data.Capacity());
    Memory<T> p(nsize, data.GetMemoryType());
    p.CopyFrom(data, size);
    p.UseDevice(data.UseDevice());
@@ -814,7 +829,7 @@ template <typename T>
 template <typename I>
 inline void Array<T>::Permute(I &&indices)
 {
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       auto current = i;
       while (i != indices[current])
@@ -832,12 +847,12 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
    SetSize(src.Size());
-   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
+   for (bigint i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
    return *this;
 }
 
 template <class T>
-inline void Array<T>::SetSize(int nsize)
+inline void Array<T>::SetSize(bigint nsize)
 {
    MFEM_ASSERT( nsize>=0, "Size must be non-negative.  It is " << nsize );
    if (nsize > Capacity())
@@ -848,7 +863,7 @@ inline void Array<T>::SetSize(int nsize)
 }
 
 template <class T>
-inline void Array<T>::SetSize(int nsize, const T &initval)
+inline void Array<T>::SetSize(bigint nsize, const T &initval)
 {
    MFEM_ASSERT( nsize>=0, "Size must be non-negative.  It is " << nsize );
    if (nsize > size)
@@ -857,7 +872,7 @@ inline void Array<T>::SetSize(int nsize, const T &initval)
       {
          GrowSize(nsize);
       }
-      for (int i = size; i < nsize; i++)
+      for (bigint i = size; i < nsize; i++)
       {
          data[i] = initval;
       }
@@ -866,7 +881,7 @@ inline void Array<T>::SetSize(int nsize, const T &initval)
 }
 
 template <class T>
-inline void Array<T>::SetSize(int nsize, MemoryType mt)
+inline void Array<T>::SetSize(bigint nsize, MemoryType mt)
 {
    MFEM_ASSERT(nsize >= 0, "invalid new size: " << nsize);
    if (mt == data.GetMemoryType())
@@ -893,7 +908,7 @@ inline void Array<T>::SetSize(int nsize, MemoryType mt)
 }
 
 template <class T>
-inline T &Array<T>::operator[](int i)
+inline T &Array<T>::operator[](bigint i)
 {
    MFEM_ASSERT( i>=0 && i<size,
                 "Access element " << i << " of array, size = " << size );
@@ -901,7 +916,7 @@ inline T &Array<T>::operator[](int i)
 }
 
 template <class T>
-inline const T &Array<T>::operator[](int i) const
+inline const T &Array<T>::operator[](bigint i) const
 {
    MFEM_ASSERT( i>=0 && i<size,
                 "Access element " << i << " of array, size = " << size );
@@ -909,7 +924,7 @@ inline const T &Array<T>::operator[](int i) const
 }
 
 template <class T>
-inline int Array<T>::Append(const T &el)
+inline bigint Array<T>::Append(const T &el)
 {
    SetSize(size+1);
    data[size-1] = el;
@@ -917,12 +932,12 @@ inline int Array<T>::Append(const T &el)
 }
 
 template <class T>
-inline int Array<T>::Append(const T *els, int nels)
+inline bigint Array<T>::Append(const T *els, bigint nels)
 {
-   const int old_size = size;
+   const bigint old_size = size;
 
    SetSize(size + nels);
-   for (int i = 0; i < nels; i++)
+   for (bigint i = 0; i < nels; i++)
    {
       data[old_size+i] = els[i];
    }
@@ -930,10 +945,10 @@ inline int Array<T>::Append(const T *els, int nels)
 }
 
 template <class T>
-inline int Array<T>::Prepend(const T &el)
+inline bigint Array<T>::Prepend(const T &el)
 {
    SetSize(size+1);
-   for (int i = size-1; i > 0; i--)
+   for (bigint i = size-1; i > 0; i--)
    {
       data[i] = data[i-1];
    }
@@ -956,9 +971,9 @@ inline const T &Array<T>::Last() const
 }
 
 template <class T>
-inline int Array<T>::Union(const T &el)
+inline bigint Array<T>::Union(const T &el)
 {
-   int i = 0;
+   bigint i = 0;
    while ((i < size) && (data[i] != el)) { i++; }
    if (i == size)
    {
@@ -968,9 +983,9 @@ inline int Array<T>::Union(const T &el)
 }
 
 template <class T>
-inline int Array<T>::Find(const T &el) const
+inline bigint Array<T>::Find(const T &el) const
 {
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       if (data[i] == el) { return i; }
    }
@@ -978,18 +993,18 @@ inline int Array<T>::Find(const T &el) const
 }
 
 template <class T>
-inline int Array<T>::FindSorted(const T &el) const
+inline bigint Array<T>::FindSorted(const T &el) const
 {
    const T *begin = data, *end = begin + size;
    const T* first = std::lower_bound(begin, end, el);
    if (first == end || !(*first == el)) { return  -1; }
-   return (int)(first - begin);
+   return (bigint)(first - begin);
 }
 
 template <class T>
 inline void Array<T>::DeleteFirst(const T &el)
 {
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       if (data[i] == el)
       {
@@ -1012,8 +1027,8 @@ inline void Array<T>::DeleteAt(const Array<int> &indices)
    Array<int> sorted_indices(indices);
    sorted_indices.Sort();
 
-   int rm_count = 0;
-   for (int i = 0; i < size; i++)
+   bigint rm_count = 0;
+   for (bigint i = 0; i < size; i++)
    {
       if (rm_count < sorted_indices.Size() && i == sorted_indices[rm_count])
       {
@@ -1050,7 +1065,7 @@ inline void Array<T>::Copy(Array &copy) const
 }
 
 template <class T>
-inline void Array<T>::MakeRef(T *data_, int size_, bool own_data)
+inline void Array<T>::MakeRef(T *data_, bigint size_, bool own_data)
 {
    data.Delete();
    data.Wrap(data_, size_, own_data);
@@ -1058,7 +1073,8 @@ inline void Array<T>::MakeRef(T *data_, int size_, bool own_data)
 }
 
 template <class T>
-inline void Array<T>::MakeRef(T *data_, int size_, MemoryType mt, bool own_data)
+inline void Array<T>::MakeRef(T *data_, bigint size_, MemoryType mt,
+                              bool own_data)
 {
    data.Delete();
    data.Wrap(data_, size_, mt, own_data);
@@ -1073,9 +1089,17 @@ inline void Array<T>::MakeRef(const Array &master)
    data.MakeAlias(master.GetMemory(), 0, size);
 }
 
+template <class T>
+inline void Array<T>::MakeRef(Memory<T> &base, bigint offset, bigint size_)
+{
+   data.Delete();
+   size = size_;
+   data.MakeAlias(base, offset, size_);
+}
+
 template <class T>
 inline void Array<T>::NewMemoryAndSize(
-   const Memory<T> &mem, int s, bool own_mem)
+   const Memory<T> &mem, bigint s, bool own_mem)
 {
    data.Delete();
    size = s;
@@ -1090,10 +1114,11 @@ inline void Array<T>::NewMemoryAndSize(
 }
 
 template <class T>
-inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
+inline void Array<T>::GetSubArray(bigint offset, bigint sa_size,
+                                  Array<T> &sa) const
 {
    sa.SetSize(sa_size);
-   for (int i = 0; i < sa_size; i++)
+   for (bigint i = 0; i < sa_size; i++)
    {
       sa[i] = (*this)[offset+i];
    }
@@ -1102,7 +1127,7 @@ inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 template <class T>
 inline void Array<T>::operator=(const T &a)
 {
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       data[i] = a;
    }

general/device.cpp

@@ -24,6 +24,18 @@
 #include <map>
 #include <sstream>
 #include <iomanip>
+#if defined(__linux__)
+#include <sys/resource.h>  // getrusage
+#include <unistd.h>        // sysconf
+#include <cstdio>          // fopen, fscanf, fclose
+#elif defined(__APPLE__)
+#include <mach/mach_init.h>  // mach_task_self
+#include <mach/task.h>       // task_info
+#elif defined(_WIN32)
+#include <windows.h>
+#include <psapi.h>  // GetProcessMemoryInfo
+#pragma comment(lib, "psapi.lib")
+#endif
 
 namespace mfem
 {
@@ -718,6 +730,52 @@ void Device::DeviceMem(size_t *free, size_t *total)
 #endif
 }
 
+// static method
+void Device::HostMem(size_t *rss_p, size_t *maxrss_p)
+{
+   size_t rss = 0, maxrss = 0;
+
+#if defined(__linux__)
+   struct rusage usage;
+   if (getrusage(RUSAGE_SELF, &usage)) { usage.ru_maxrss = 0; }
+   maxrss = 1024*usage.ru_maxrss;
+   static const long PAGE_SIZE = sysconf(_SC_PAGESIZE);
+   FILE *statm = fopen("/proc/self/statm", "r");
+   if (statm)
+   {
+      // Values are measured in pages, see Table 1-3 at
+      // https://www.kernel.org/doc/Documentation/filesystems/proc.txt
+      long rss_pages;
+      if (fscanf(statm, "%*d %ld", &rss_pages) == EOF) { rss_pages = 0; }
+      fclose(statm);
+      rss = rss_pages * PAGE_SIZE;
+   }
+#elif defined(__APPLE__)
+   struct mach_task_basic_info info;
+   mach_msg_type_number_t count = MACH_TASK_BASIC_INFO_COUNT;
+   if (task_info(mach_task_self(), MACH_TASK_BASIC_INFO,
+                 (task_info_t)&info, &count) == KERN_SUCCESS)
+   {
+      rss = info.resident_size;
+      maxrss = info.resident_size_max;
+   }
+#elif defined(_WIN32)
+   PROCESS_MEMORY_COUNTERS mem_counters;
+   if (GetProcessMemoryInfo(GetCurrentProcess(),
+                            &mem_counters,
+                            sizeof(mem_counters)))
+   {
+      // Reference:
+      // https://learn.microsoft.com/en-us/windows/win32/api/psapi/ns-psapi-process_memory_counters
+      rss = mem_counters.WorkingSetSize;
+      maxrss = mem_counters.PeakWorkingSetSize;
+   }
+#endif
+
+   *rss_p = rss;
+   *maxrss_p = maxrss;
+}
+
 std::string Device::GetUUID(const int device_id)
 {
    std::stringstream res;

general/device.hpp

@@ -321,6 +321,10 @@ public:
 
    /** @brief Gets the @a free and @a total memory on the device. */
    static void DeviceMem(size_t *free, size_t *total);
+
+   /** @brief Gets the @a rss (resident set size) and @a maxrss (maximum
+       resident set size) memory of the host process, in bytes. */
+   static void HostMem(size_t *rss, size_t *maxrss);
 };
 
 
@@ -349,14 +353,14 @@ inline MemoryClass GetMemoryClass(const Memory<T> &mem, bool on_dev)
     HostMemoryClass, otherwise. */
 /** Also, if @a on_dev = true, the device flag of @a mem will be set. */
 template <typename T>
-inline const T *Read(const Memory<T> &mem, int size, bool on_dev = true)
+inline const T *Read(const Memory<T> &mem, bigint size, bool on_dev = true)
 {
    return mem.Read(GetMemoryClass(mem, on_dev), size);
 }
 
-/** @brief Shortcut to Read(const Memory<T> &mem, int size, false) */
+/** @brief Shortcut to Read(const Memory<T> &mem, bigint size, false) */
 template <typename T>
-inline const T *HostRead(const Memory<T> &mem, int size)
+inline const T *HostRead(const Memory<T> &mem, bigint size)
 {
    return mfem::Read(mem, size, false);
 }
@@ -366,14 +370,14 @@ inline const T *HostRead(const Memory<T> &mem, int size)
     HostMemoryClass, otherwise. */
 /** Also, if @a on_dev = true, the device flag of @a mem will be set. */
 template <typename T>
-inline T *Write(Memory<T> &mem, int size, bool on_dev = true)
+inline T *Write(Memory<T> &mem, bigint size, bool on_dev = true)
 {
    return mem.Write(GetMemoryClass(mem, on_dev), size);
 }
 
-/** @brief Shortcut to Write(const Memory<T> &mem, int size, false) */
+/** @brief Shortcut to Write(const Memory<T> &mem, bigint size, false) */
 template <typename T>
-inline T *HostWrite(Memory<T> &mem, int size)
+inline T *HostWrite(Memory<T> &mem, bigint size)
 {
    return mfem::Write(mem, size, false);
 }
@@ -383,14 +387,14 @@ inline T *HostWrite(Memory<T> &mem, int size)
     HostMemoryClass, otherwise. */
 /** Also, if @a on_dev = true, the device flag of @a mem will be set. */
 template <typename T>
-inline T *ReadWrite(Memory<T> &mem, int size, bool on_dev = true)
+inline T *ReadWrite(Memory<T> &mem, bigint size, bool on_dev = true)
 {
    return mem.ReadWrite(GetMemoryClass(mem, on_dev), size);
 }
 
-/** @brief Shortcut to ReadWrite(Memory<T> &mem, int size, false) */
+/** @brief Shortcut to ReadWrite(Memory<T> &mem, bigint size, false) */
 template <typename T>
-inline T *HostReadWrite(Memory<T> &mem, int size)
+inline T *HostReadWrite(Memory<T> &mem, bigint size)
 {
    return mfem::ReadWrite(mem, size, false);
 }

general/error.hpp

@@ -68,7 +68,8 @@ void mfem_error(const char *msg = NULL);
 __attribute__((enzyme_inactive))
 #endif
 void mfem_warning(const char *msg = NULL);
-}
+
+} // namespace mfem
 
 #ifndef _MFEM_FUNC_NAME
 #ifndef _MSC_VER

general/forall.hpp

@@ -176,10 +176,10 @@ private:
 // with CUDA/HIP language. Otherwise, this macro is a no-op.
 #if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #define MFEM_GPU_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
-                                       (int i) {__VA_ARGS__})
+                                       (bigint i) {__VA_ARGS__})
 #elif defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_GPU_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
-                                        (int i) {__VA_ARGS__})
+                                        (bigint i) {__VA_ARGS__})
 #else
 #define MFEM_GPU_FORALL(i, N,...) do { } while (false)
 #endif
@@ -189,7 +189,7 @@ private:
 
 // The MFEM_FORALL wrapper
 #define MFEM_FORALL(i,N,...) \
-   ForallWrap<1>(true,N,[=] MFEM_HOST_DEVICE (int i) {__VA_ARGS__})
+   ForallWrap<1>(true,N,[=] MFEM_HOST_DEVICE (bigint i) {__VA_ARGS__})
 
 // MFEM_FORALL with a 2D CUDA block
 #define MFEM_FORALL_2D(i,N,X,Y,BZ,...) \
@@ -208,16 +208,16 @@ private:
 // example the functions in vector.cpp, where we don't want to use the mfem
 // device for operations on small vectors.
 #define MFEM_FORALL_SWITCH(use_dev,i,N,...) \
-   ForallWrap<1>(use_dev,N,[=] MFEM_HOST_DEVICE (int i) {__VA_ARGS__})
+   ForallWrap<1>(use_dev,N,[=] MFEM_HOST_DEVICE (bigint i) {__VA_ARGS__})
 
 
 /// OpenMP backend
 template <typename HBODY>
-void OmpWrap(const int N, HBODY &&h_body)
+void OmpWrap(const bigint N, HBODY &&h_body)
 {
 #ifdef MFEM_USE_OPENMP
    #pragma omp parallel for
-   for (int k = 0; k < N; k++)
+   for (bigint k = 0; k < N; k++)
    {
       h_body(k);
    }
@@ -296,7 +296,7 @@ using hip_threads_z =
 
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_CUDA) && defined(__CUDACC__)
 template <typename DBODY>
-void RajaCuWrap1D(const int N, DBODY &&d_body)
+void RajaCuWrap1D(const bigint N, DBODY &&d_body)
 {
    //true denotes asynchronous kernel
    RAJA::forall<RAJA::cuda_exec<MFEM_CUDA_BLOCKS,true>>(RAJA::RangeSegment(0,N),
@@ -364,7 +364,7 @@ template <>
 struct RajaCuWrap<1>
 {
    template <typename DBODY>
-   static void run(const int N, DBODY &&d_body,
+   static void run(const bigint N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
       RajaCuWrap1D(N, d_body);
@@ -397,7 +397,7 @@ struct RajaCuWrap<3>
 
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_HIP) && defined(__HIP__)
 template <typename DBODY>
-void RajaHipWrap1D(const int N, DBODY &&d_body)
+void RajaHipWrap1D(const bigint N, DBODY &&d_body)
 {
    //true denotes asynchronous kernel
    RAJA::forall<RAJA::hip_exec<MFEM_HIP_BLOCKS,true>>(RAJA::RangeSegment(0,N),
@@ -465,7 +465,7 @@ template <>
 struct RajaHipWrap<1>
 {
    template <typename DBODY>
-   static void run(const int N, DBODY &&d_body,
+   static void run(const bigint N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
       RajaHipWrap1D(N, d_body);
@@ -500,7 +500,7 @@ struct RajaHipWrap<3>
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_OPENMP)
 
 template <typename HBODY>
-void RajaOmpWrap(const int N, HBODY &&h_body)
+void RajaOmpWrap(const bigint N, HBODY &&h_body)
 {
    RAJA::forall<RAJA::omp_parallel_for_exec>(RAJA::RangeSegment(0,N), h_body);
 }
@@ -546,7 +546,7 @@ void RajaOmpWrap3D(const int Nx, const int Ny, const int Nz, HBODY &&h_body)
 
 /// RAJA sequential loop backend
 template <typename HBODY>
-void RajaSeqWrap(const int N, HBODY &&h_body)
+void RajaSeqWrap(const bigint N, HBODY &&h_body)
 {
 #ifdef MFEM_USE_RAJA
 
@@ -571,9 +571,9 @@ void RajaSeqWrap(const int N, HBODY &&h_body)
 #if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 
 template <typename BODY> __global__ static
-void CuKernel1D(const int N, BODY body)
+void CuKernel1D(const bigint N, BODY body)
 {
-   const int k = blockDim.x*blockIdx.x + threadIdx.x;
+   const bigint k = bigint(blockDim.x)*blockIdx.x + threadIdx.x;
    if (k >= N) { return; }
    body(k);
 }
@@ -612,10 +612,10 @@ static void CuKernel3DLaunchBounds(const int N, BODY body)
 }
 
 template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
-void CuWrap1D(const int N, DBODY &&d_body)
+void CuWrap1D(const bigint N, DBODY &&d_body)
 {
    if (N==0) { return; }
-   const int GRID = (N+BLCK-1)/BLCK;
+   const unsigned int GRID = (N+BLCK-1)/BLCK;
    CuKernel1D<<<GRID,BLCK>>>(N, d_body);
    MFEM_GPU_CHECK(cudaGetLastError());
 }
@@ -676,7 +676,7 @@ template <int MAX_THREADS_PER_BLOCK>
 struct CuWrap<1, MAX_THREADS_PER_BLOCK>
 {
    template <typename DBODY>
-   static void run(const int N, DBODY &&d_body,
+   static void run(const bigint N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
       CuWrap1D<MFEM_CUDA_BLOCKS>(N, d_body);
@@ -735,9 +735,9 @@ struct CuWrap<3, MAX_THREADS_PER_BLOCK>
 #if defined(MFEM_USE_HIP) && defined(__HIP__)
 
 template <typename BODY> __global__ static
-void HipKernel1D(const int N, BODY body)
+void HipKernel1D(const bigint N, BODY body)
 {
-   const int k = hipBlockDim_x*hipBlockIdx_x + hipThreadIdx_x;
+   const bigint k = bigint(hipBlockDim_x)*hipBlockIdx_x + hipThreadIdx_x;
    if (k >= N) { return; }
    body(k);
 }
@@ -775,10 +775,10 @@ static void HipKernel3DLaunchBounds(const int N, BODY body)
 }
 
 template <int BLCK = MFEM_HIP_BLOCKS, typename DBODY>
-void HipWrap1D(const int N, DBODY &&d_body)
+void HipWrap1D(const bigint N, DBODY &&d_body)
 {
    if (N==0) { return; }
-   const int GRID = (N+BLCK-1)/BLCK;
+   const unsigned int GRID = (N+BLCK-1)/BLCK;
    hipLaunchKernelGGL(HipKernel1D,GRID,BLCK,0,nullptr,N,d_body);
    MFEM_GPU_CHECK(hipGetLastError());
 }
@@ -839,7 +839,7 @@ template <int MAX_THREADS_PER_BLOCK>
 struct HipWrap<1, MAX_THREADS_PER_BLOCK>
 {
    template <typename DBODY>
-   static void run(const int N, DBODY &&d_body,
+   static void run(const bigint N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
       HipWrap1D<MFEM_HIP_BLOCKS>(N, d_body);
@@ -897,7 +897,7 @@ struct HipWrap<3, MAX_THREADS_PER_BLOCK>
 /// Forall host & device kernel dispatch
 template <int DIM, int MAX_THREADS_PER_BLOCK = 0,
           typename d_lambda, typename h_lambda>
-inline void ForallWrap(const bool use_dev, const int N,
+inline void ForallWrap(const bool use_dev, const bigint N,
                        d_lambda &&d_body, h_lambda &&h_body,
                        const int X=0, const int Y=0, const int Z=0,
                        const int G=0)
@@ -963,13 +963,13 @@ backend_cpu:
    // Handle Backend::CPU. This is also a fallback for any allowed backends not
    // handled above, e.g. OCCA_CPU with configuration 'occa-cpu,cpu', or
    // OCCA_OMP with configuration 'occa-omp,cpu'.
-   for (int k = 0; k < N; k++) { h_body(k); }
+   for (bigint k = 0; k < N; k++) { h_body(k); }
 }
 
 ///////////////////////////////////////////////////////////////////////////////
 /// Forall host & device kernel wrappers
 template <int DIM, typename lambda>
-inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
+inline void ForallWrap(const bool use_dev, const bigint N, lambda &&body,
                        const int X=0, const int Y=0, const int Z=0,
                        const int G=0)
 {
@@ -977,7 +977,7 @@ inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
 }
 
 template <int DIM, int MAX_THREADS_PER_BLOCK, typename lambda>
-inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
+inline void ForallWrap(const bool use_dev, const bigint N, lambda &&body,
                        const int X=0, const int Y=0, const int Z=0,
                        const int G=0)
 {
@@ -987,14 +987,14 @@ inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
 ///////////////////////////////////////////////////////////////////////////////
 // forall interfaces
 template<typename lambda>
-inline void forall(int N, lambda &&body) { ForallWrap<1>(true, N, body); }
+inline void forall(bigint N, lambda &&body) { ForallWrap<1>(true, N, body); }
 
 template<typename lambda>
 inline void forall(int Nx, int Ny, lambda &&body)
 {
    if (Device::Allows(Backend::DEVICE_MASK))
    {
-      mfem::forall(Nx * Ny, [=] MFEM_HOST_DEVICE(int idx)
+      mfem::forall(bigint(Nx) * Ny, [=] MFEM_HOST_DEVICE(bigint idx)
       {
          int j = idx / Nx;
          int i = idx % Nx;
@@ -1030,12 +1030,12 @@ inline void forall(int Nx, int Ny, int Nz, lambda &&body)
 {
    if (Device::Allows(Backend::DEVICE_MASK))
    {
-      mfem::forall(Nx * Ny * Nz, [=] MFEM_HOST_DEVICE(int idx)
+      mfem::forall(bigint(Nx) * Ny * Nz, [=] MFEM_HOST_DEVICE(bigint idx)
       {
          int i = idx % Nx;
-         int j = idx / Nx;
-         int k = j / Ny;
-         j = j % Ny;
+         bigint jk = idx / Nx;
+         int k = jk / Ny;
+         int j = jk % Ny;
          body(i, j, k);
       });
    }
@@ -1067,7 +1067,7 @@ inline void forall(int Nx, int Ny, int Nz, lambda &&body)
 }
 
 template<typename lambda>
-inline void forall_switch(bool use_dev, int N, lambda &&body)
+inline void forall_switch(bool use_dev, bigint N, lambda &&body)
 {
    ForallWrap<1>(use_dev, N, body);
 }

general/globals.hpp

@@ -113,6 +113,10 @@ void SetGlobalMPI_Comm(MPI_Comm comm);
 /// to suppress the warning.
 const char* GetEnv(const char* name);
 
+/// Signed indices i -> -1 - i are used as a convention to encode orientation.
+inline MFEM_HOST_DEVICE int FlipIndexSign(int i) { return -1 - i; }
+inline MFEM_HOST_DEVICE int UnsignIndex(int i) { return i < 0 ? -1 - i : i; }
+
 } // namespace mfem
 
 #endif

general/mem_manager.cpp

@@ -150,8 +150,8 @@ template void Memory<int>::PrintFlags() const;
 template void Memory<real_t>::PrintFlags() const;
 
 // Instantiate Memory<T>::CompareHostAndDevice for T = int and T = real_t.
-template int Memory<int>::CompareHostAndDevice(int size) const;
-template int Memory<real_t>::CompareHostAndDevice(int size) const;
+template int Memory<int>::CompareHostAndDevice(bigint size) const;
+template int Memory<real_t>::CompareHostAndDevice(bigint size) const;
 
 
 namespace internal
@@ -759,7 +759,9 @@ private:
    {
       switch (mt)
       {
-         case MT::HOST_DEBUG: return new MmuHostMemorySpace();
+         case MT::HOST_DEBUG:
+            if (GetEnv("MFEM_MMU_STD")) { return new StdHostMemorySpace(); }
+            return new MmuHostMemorySpace();
 #ifdef MFEM_USE_UMPIRE
          case MT::HOST_UMPIRE:
             return new UmpireHostMemorySpace(
@@ -788,7 +790,9 @@ private:
          case MT::DEVICE_UMPIRE: return new NoDeviceMemorySpace();
          case MT::DEVICE_UMPIRE_2: return new NoDeviceMemorySpace();
 #endif
-         case MT::DEVICE_DEBUG: return new MmuDeviceMemorySpace();
+         case MT::DEVICE_DEBUG:
+            if (GetEnv("MFEM_MMU_STD")) { return new StdDeviceMemorySpace(); }
+            return new MmuDeviceMemorySpace();
          case MT::DEVICE:
          {
 #if defined(MFEM_USE_CUDA)

general/mem_manager.hpp

@@ -197,7 +197,7 @@ protected:
    /** The type of the pointer is given by the field #h_mt; it can be any type
        from MemoryClass::HOST. */
    T *h_ptr;
-   int capacity; ///< Size of the allocated memory
+   bigint capacity; ///< Size of the allocated memory
    MemoryType h_mt; ///< Host memory type
    mutable unsigned flags; ///< Bit flags defined from the #FlagMask enum
    // 'flags' is mutable so that it can be modified in Set{Host,Device}PtrOwner,
@@ -236,7 +236,7 @@ public:
    /// Allocate host memory for @a size entries.
    /** The allocation uses the current host memory type returned by
        MemoryManager::GetHostMemoryType(). */
-   explicit Memory(int size) { New(size); }
+   explicit Memory(bigint size) { New(size); }
 
    /// Creates a new empty Memory object with host MemoryType @a mt.
    explicit Memory(MemoryType mt) { Reset(mt); }
@@ -245,19 +245,20 @@ public:
        @a mt. */
    /** The newly allocated memory is not initialized, however the given
        MemoryType is still set as valid. */
-   Memory(int size, MemoryType mt) { New(size, mt); }
+   Memory(bigint size, MemoryType mt) { New(size, mt); }
 
    /** @brief Allocate memory for @a size entries with the given host MemoryType
        @a h_mt and device MemoryType @a d_mt. */
    /** The newly allocated memory is not initialized. The host pointer is set as
        valid. */
-   Memory(int size, MemoryType h_mt, MemoryType d_mt) { New(size, h_mt, d_mt); }
+   Memory(bigint size, MemoryType h_mt, MemoryType d_mt)
+   { New(size, h_mt, d_mt); }
 
    /** @brief Wrap an externally allocated host pointer, @a ptr with the current
        host memory type returned by MemoryManager::GetHostMemoryType(). */
    /** The parameter @a own determines whether @a ptr will be deleted when the
        method Delete() is called. */
-   explicit Memory(T *ptr, int size, bool own) { Wrap(ptr, size, own); }
+   explicit Memory(T *ptr, bigint size, bool own) { Wrap(ptr, size, own); }
 
    /// Wrap an externally allocated pointer, @a ptr, of the given MemoryType.
    /** The new memory object will have the given MemoryType set as valid.
@@ -267,13 +268,13 @@ public:
 
        The parameter @a own determines whether @a ptr will be deleted when the
        method Delete() is called. */
-   Memory(T *ptr, int size, MemoryType mt, bool own)
+   Memory(T *ptr, bigint size, MemoryType mt, bool own)
    { Wrap(ptr, size, mt, own); }
 
    /** @brief Alias constructor. Create a Memory object that points inside the
        Memory object @a base. */
    /** The new Memory object uses the same MemoryType(s) as @a base. */
-   Memory(const Memory &base, int offset, int size)
+   Memory(const Memory &base, bigint offset, bigint size)
    { MakeAlias(base, offset, size); }
 
    /// Destructor: default.
@@ -319,7 +320,7 @@ public:
    { flags = use_dev ? (flags | USE_DEVICE) : (flags & ~USE_DEVICE); }
 
    /// Return the size of the allocated memory.
-   int Capacity() const { return capacity; }
+   bigint Capacity() const { return capacity; }
 
    /// Reset the memory to be empty, ensuring that Delete() will be a no-op.
    /** This is the Memory class equivalent to setting a pointer to NULL, see
@@ -339,7 +340,7 @@ public:
    /** @brief Allocate host memory for @a size entries with the current host
        memory type returned by MemoryManager::GetHostMemoryType(). */
    /** @note The current memory is NOT deleted by this method. */
-   inline void New(int size);
+   inline void New(bigint size);
 
    /// Allocate memory for @a size entries with the given MemoryType.
    /** The newly allocated memory is not initialized, however the given
@@ -353,7 +354,7 @@ public:
        to be the dual of @a mt, see MemoryManager::GetDualMemoryType().
 
        @note The current memory is NOT deleted by this method. */
-   inline void New(int size, MemoryType mt);
+   inline void New(bigint size, MemoryType mt);
 
    /** @brief Allocate memory for @a size entries with the given host MemoryType
        @a h_mt and device MemoryType @a d_mt. */
@@ -361,7 +362,7 @@ public:
        valid.
 
        @note The current memory is NOT deleted by this method. */
-   inline void New(int size, MemoryType h_mt, MemoryType d_mt);
+   inline void New(bigint size, MemoryType h_mt, MemoryType d_mt);
 
    /** @brief Wrap an externally allocated host pointer, @a ptr with the current
        host memory type returned by MemoryManager::GetHostMemoryType(). */
@@ -369,7 +370,7 @@ public:
        method Delete() is called.
 
        @note The current memory is NOT deleted by this method. */
-   inline void Wrap(T *ptr, int size, bool own);
+   inline void Wrap(T *ptr, bigint size, bool own);
 
    /// Wrap an externally allocated pointer, @a ptr, of the given MemoryType.
    /** The new memory object will have the given MemoryType set as valid.
@@ -381,7 +382,7 @@ public:
        method Delete() is called.
 
        @note The current memory is NOT deleted by this method. */
-   inline void Wrap(T *ptr, int size, MemoryType mt, bool own);
+   inline void Wrap(T *ptr, bigint size, MemoryType mt, bool own);
 
    /** Wrap an externally pair of allocated pointers, @a h_ptr and @a d_ptr,
        of the given host MemoryType @a h_mt. */
@@ -405,14 +406,14 @@ public:
          - SetDevicePtrOwner.
 
        @note The current memory is NOT deleted by this method. */
-   inline void Wrap(T *h_ptr, T *d_ptr, int size, MemoryType h_mt, bool own,
+   inline void Wrap(T *h_ptr, T *d_ptr, bigint size, MemoryType h_mt, bool own,
                     bool valid_host = false, bool valid_device = true);
 
    /// Create a memory object that points inside the memory object @a base.
    /** The new Memory object uses the same MemoryType(s) as @a base.
 
        @note The current memory is NOT deleted by this method. */
-   inline void MakeAlias(const Memory &base, int offset, int size);
+   inline void MakeAlias(const Memory &base, bigint offset, bigint size);
 
    /// Set the device MemoryType to be used by the Memory object.
    /** If the specified @a d_mt is not a device MemoryType, i.e. not one of the
@@ -437,10 +438,10 @@ public:
    inline void DeleteDevice(bool copy_to_host = true);
 
    /// Array subscript operator for host memory.
-   inline T &operator[](int idx);
+   inline T &operator[](bigint idx);
 
    /// Array subscript operator for host memory, const version.
-   inline const T &operator[](int idx) const;
+   inline const T &operator[](bigint idx) const;
 
    /// Direct access to the host memory as T* (implicit conversion).
    /** When the type T is const-qualified, this method can be used only if the
@@ -492,11 +493,11 @@ public:
        Read() or Write() should be used instead of this method.
 
        The parameter @a size must not exceed the Capacity(). */
-   inline T *ReadWrite(MemoryClass mc, int size);
+   inline T *ReadWrite(MemoryClass mc, bigint size);
 
    /// Get read-only access to the memory with the given MemoryClass.
    /** The parameter @a size must not exceed the Capacity(). */
-   inline const T *Read(MemoryClass mc, int size) const;
+   inline const T *Read(MemoryClass mc, bigint size) const;
 
    /// Get write-only access to the memory with the given MemoryClass.
    /** The parameter @a size must not exceed the Capacity().
@@ -504,7 +505,7 @@ public:
        The contents of the returned pointer is undefined, unless it was
        validated by a previous call to Read() or ReadWrite() with
        the same MemoryClass. */
-   inline T *Write(MemoryClass mc, int size);
+   inline T *Write(MemoryClass mc, bigint size);
 
    /// Copy the host/device pointer validity flags from @a other to @a *this.
    /** This method synchronizes the pointer validity flags of two Memory objects
@@ -521,7 +522,7 @@ public:
        of the base incorrect. Calling this method will ensure that @a base is
        up-to-date. Note that this is achieved by moving/copying @a *this (if
        necessary), and not @a base. */
-   inline void SyncAlias(const Memory &base, int alias_size) const;
+   inline void SyncAlias(const Memory &base, bigint alias_size) const;
 
    /** @brief Return a MemoryType that is currently valid. If both the host and
        the device pointers are currently valid, then the device memory type is
@@ -544,20 +545,20 @@ public:
    /// Copy @a size entries from @a src to @a *this.
    /** The given @a size should not exceed the Capacity() of the source @a src
        and the destination, @a *this. */
-   inline void CopyFrom(const Memory &src, int size);
+   inline void CopyFrom(const Memory &src, bigint size);
 
    /// Copy @a size entries from the host pointer @a src to @a *this.
    /** The given @a size should not exceed the Capacity() of @a *this. */
-   inline void CopyFromHost(const T *src, int size);
+   inline void CopyFromHost(const T *src, bigint size);
 
    /// Copy @a size entries from @a *this to @a dest.
    /** The given @a size should not exceed the Capacity() of @a *this and the
        destination, @a dest. */
-   inline void CopyTo(Memory &dest, int size) const;
+   inline void CopyTo(Memory &dest, bigint size) const;
 
    /// Copy @a size entries from @a *this to the host pointer @a dest.
    /** The given @a size should not exceed the Capacity() of @a *this. */
-   inline void CopyToHost(T *dest, int size) const;
+   inline void CopyToHost(T *dest, bigint size) const;
 
    /// Print the internal flags.
    /** This method can be useful for debugging. It is explicitly instantiated
@@ -567,7 +568,7 @@ public:
    /// If both the host and the device data are valid, compare their contents.
    /** This method can be useful for debugging. It is explicitly instantiated
        for Memory<T> with T = int and T = real_t. */
-   inline int CompareHostAndDevice(int size) const;
+   inline int CompareHostAndDevice(bigint size) const;
 
 private:
    // GCC 4.8 workaround: max_align_t is not in std.
@@ -956,7 +957,7 @@ inline void Memory<T>::Reset(MemoryType host_mt)
 }
 
 template <typename T>
-inline void Memory<T>::New(int size)
+inline void Memory<T>::New(bigint size)
 {
    capacity = size;
    flags = OWNS_HOST | VALID_HOST;
@@ -966,7 +967,7 @@ inline void Memory<T>::New(int size)
 }
 
 template <typename T>
-inline void Memory<T>::New(int size, MemoryType mt)
+inline void Memory<T>::New(bigint size, MemoryType mt)
 {
    capacity = size;
    const size_t bytes = size*sizeof(T);
@@ -978,7 +979,8 @@ inline void Memory<T>::New(int size, MemoryType mt)
 }
 
 template <typename T>
-inline void Memory<T>::New(int size, MemoryType host_mt, MemoryType device_mt)
+inline void Memory<T>::New(bigint size, MemoryType host_mt,
+                           MemoryType device_mt)
 {
    capacity = size;
    const size_t bytes = size*sizeof(T);
@@ -989,7 +991,7 @@ inline void Memory<T>::New(int size, MemoryType host_mt, MemoryType device_mt)
 }
 
 template <typename T>
-inline void Memory<T>::Wrap(T *ptr, int size, bool own)
+inline void Memory<T>::Wrap(T *ptr, bigint size, bool own)
 {
    h_ptr = ptr;
    capacity = size;
@@ -1011,7 +1013,7 @@ inline void Memory<T>::Wrap(T *ptr, int size, bool own)
 }
 
 template <typename T>
-inline void Memory<T>::Wrap(T *ptr, int size, MemoryType mt, bool own)
+inline void Memory<T>::Wrap(T *ptr, bigint size, MemoryType mt, bool own)
 {
    capacity = size;
    if (IsHostMemory(mt))
@@ -1036,7 +1038,7 @@ inline void Memory<T>::Wrap(T *ptr, int size, MemoryType mt, bool own)
 }
 
 template <typename T>
-inline void Memory<T>::Wrap(T *h_ptr_, T *d_ptr, int size, MemoryType h_mt_,
+inline void Memory<T>::Wrap(T *h_ptr_, T *d_ptr, bigint size, MemoryType h_mt_,
                             bool own, bool valid_host, bool valid_device)
 {
    h_mt = h_mt_;
@@ -1053,7 +1055,7 @@ inline void Memory<T>::Wrap(T *h_ptr_, T *d_ptr, int size, MemoryType h_mt_,
 }
 
 template <typename T>
-inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
+inline void Memory<T>::MakeAlias(const Memory &base, bigint offset, bigint size)
 {
    MFEM_ASSERT(0 <= offset, "invalid offset = " << offset);
    MFEM_ASSERT(0 <= size, "invalid size = " << size);
@@ -1140,7 +1142,7 @@ inline void Memory<T>::DeleteDevice(bool copy_to_host)
 }
 
 template <typename T>
-inline T &Memory<T>::operator[](int idx)
+inline T &Memory<T>::operator[](bigint idx)
 {
    MFEM_ASSERT((flags & VALID_HOST) && !(flags & VALID_DEVICE),
                "invalid host pointer access");
@@ -1148,7 +1150,7 @@ inline T &Memory<T>::operator[](int idx)
 }
 
 template <typename T>
-inline const T &Memory<T>::operator[](int idx) const
+inline const T &Memory<T>::operator[](bigint idx) const
 {
    MFEM_ASSERT((flags & VALID_HOST), "invalid host pointer access");
    return h_ptr[idx];
@@ -1189,7 +1191,7 @@ inline Memory<T>::operator const U*() const
 }
 
 template <typename T>
-inline T *Memory<T>::ReadWrite(MemoryClass mc, int size)
+inline T *Memory<T>::ReadWrite(MemoryClass mc, bigint size)
 {
    const size_t bytes = size * sizeof(T);
    if (!(flags & Registered))
@@ -1202,7 +1204,7 @@ inline T *Memory<T>::ReadWrite(MemoryClass mc, int size)
 }
 
 template <typename T>
-inline const T *Memory<T>::Read(MemoryClass mc, int size) const
+inline const T *Memory<T>::Read(MemoryClass mc, bigint size) const
 {
    const size_t bytes = size * sizeof(T);
    if (!(flags & Registered))
@@ -1215,7 +1217,7 @@ inline const T *Memory<T>::Read(MemoryClass mc, int size) const
 }
 
 template <typename T>
-inline T *Memory<T>::Write(MemoryClass mc, int size)
+inline T *Memory<T>::Write(MemoryClass mc, bigint size)
 {
    const size_t bytes = size * sizeof(T);
    if (!(flags & Registered))
@@ -1242,7 +1244,7 @@ inline void Memory<T>::Sync(const Memory &other) const
 }
 
 template <typename T>
-inline void Memory<T>::SyncAlias(const Memory &base, int alias_size) const
+inline void Memory<T>::SyncAlias(const Memory &base, bigint alias_size) const
 {
    // Assuming that if *this is registered then base is also registered.
    MFEM_ASSERT(!(flags & Registered) || (base.flags & Registered),
@@ -1279,7 +1281,7 @@ inline bool Memory<T>::DeviceIsValid() const
 }
 
 template <typename T>
-inline void Memory<T>::CopyFrom(const Memory &src, int size)
+inline void Memory<T>::CopyFrom(const Memory &src, bigint size)
 {
    MFEM_VERIFY(src.capacity>=size && capacity>=size, "Incorrect size");
    if (size <= 0) { return; }
@@ -1300,7 +1302,7 @@ inline void Memory<T>::CopyFrom(const Memory &src, int size)
 }
 
 template <typename T>
-inline void Memory<T>::CopyFromHost(const T *src, int size)
+inline void Memory<T>::CopyFromHost(const T *src, bigint size)
 {
    MFEM_VERIFY(capacity>=size, "Incorrect size");
    if (size <= 0) { return; }
@@ -1321,13 +1323,13 @@ inline void Memory<T>::CopyFromHost(const T *src, int size)
 }
 
 template <typename T>
-inline void Memory<T>::CopyTo(Memory &dest, int size) const
+inline void Memory<T>::CopyTo(Memory &dest, bigint size) const
 {
    dest.CopyFrom(*this, size);
 }
 
 template <typename T>
-inline void Memory<T>::CopyToHost(T *dest, int size) const
+inline void Memory<T>::CopyToHost(T *dest, bigint size) const
 {
    MFEM_VERIFY(capacity>=size, "Incorrect size");
    if (size <= 0) { return; }
@@ -1359,7 +1361,7 @@ inline void Memory<T>::PrintFlags() const
 }
 
 template <typename T>
-inline int Memory<T>::CompareHostAndDevice(int size) const
+inline int Memory<T>::CompareHostAndDevice(bigint size) const
 {
    if (!(flags & VALID_HOST) || !(flags & VALID_DEVICE)) { return 0; }
    return MemoryManager::CompareHostAndDevice_(h_ptr, size*sizeof(T), flags);

general/stable3d.cpp

@@ -12,6 +12,7 @@
 
 #include "error.hpp"
 #include "stable3d.hpp"
+#include <limits>
 
 using namespace std;
 
@@ -90,6 +91,8 @@ int STable3D::Push (int r, int c, int f)
    node->Prev   = Rows[r];
    Rows[r] = node;
 
+   MFEM_VERIFY(NElem != std::numeric_limits<int>::max(),
+               "integer overflow error");
    NElem++;
    return (NElem-1);
 }

general/table.cpp

@@ -18,6 +18,7 @@
 #include "../general/mem_manager.hpp"
 #include <iostream>
 #include <iomanip>
+#include <limits>
 
 namespace mfem
 {
@@ -27,6 +28,8 @@ using namespace std;
 Table::Table(const Table &table1,
              const Table &table2, int offset)
 {
+   MFEM_VERIFY(!table1.UsingBigI() && !table2.UsingBigI(), "");
+
    MFEM_ASSERT(table1.size == table2.size,
                "Tables have different sizes can not merge.");
    size = table1.size;
@@ -60,6 +63,9 @@ Table::Table(const Table &table1,
              const Table &table2, int offset2,
              const Table &table3, int offset3)
 {
+   MFEM_VERIFY(!table1.UsingBigI() && !table2.UsingBigI() &&
+               !table3.UsingBigI(), "");
+
    MFEM_ASSERT(table1.size == table2.size,
                "Tables have different sizes can not merge.");
    MFEM_ASSERT(table1.size == table3.size,
@@ -98,17 +104,30 @@ Table::Table(const Table &table1,
 
 Table::Table (int dim, int connections_per_row)
 {
-   int i, j, sum = dim * connections_per_row;
+   bigint sum = bigint(dim) * connections_per_row;
 
    size = dim;
-   I.SetSize(size+1);
-   J.SetSize(sum);
-
-   I[0] = 0;
-   for (i = 1; i <= size; i++)
+   if (int(sum) == sum)
    {
-      I[i] = I[i-1] + connections_per_row;
-      for (j = I[i-1]; j < I[i]; j++) { J[j] = -1; }
+      I.SetSize(size+1);
+      J.SetSize(sum);
+      I[0] = 0;
+      for (int i = 1; i <= size; i++)
+      {
+         I[i] = I[i-1] + connections_per_row;
+         for (int j = I[i-1]; j < I[i]; j++) { J[j] = -1; }
+      }
+   }
+   else
+   {
+      bigI.SetSize(size+1);
+      J.SetSize(sum);
+      bigI[0] = 0;
+      for (int i = 1; i <= size; i++)
+      {
+         bigI[i] = bigI[i-1] + connections_per_row;
+         for (bigint j = bigI[i-1]; j < bigI[i]; j++) { J[j] = -1; }
+      }
    }
 }
 
@@ -139,77 +158,153 @@ void Table::MakeI(int nrows)
 
 void Table::MakeJ()
 {
-   int i, j, k;
+   bigint nnz;
 
-   for (k = i = 0; i < size; i++)
+   if (!UsingBigI())
    {
-      j = I[i], I[i] = k, k += j;
+      nnz = 0;
+      int nnz_int = 0;
+      for (int i = 0; i < size; i++)
+      {
+         const int row_size = I[i];
+         I[i] = nnz_int;
+         nnz_int += row_size;
+         nnz += row_size;
+         if (nnz_int != nnz) // check for overflow
+         {
+            bigI.SetSize(size+1);
+            for (int j = 0; j <= i; j++) { bigI[j] = I[j]; }
+            for (i++ ; i < size; i++)
+            {
+               bigI[i] = nnz;
+               nnz += I[i];
+            }
+            bigI[size] = nnz;
+            I.DeleteAll();
+            goto I_is_updated;
+         }
+      }
+      I[size] = nnz_int;
+      nnz = nnz_int;
+   I_is_updated: ;
+   }
+   else
+   {
+      nnz = 0;
+      for (int i = 0; i < size; i++)
+      {
+         const bigint row_size = bigI[i];
+         bigI[i] = nnz;
+         nnz += row_size;
+      }
+      bigI[size] = nnz;
    }
 
-   J.SetSize(I[size]=k);
+   J.SetSize(nnz);
 }
 
 void Table::AddConnections(int r, const int *c, int nc)
 {
-   int *jp = J+I[r];
+   int *jp = GetRow(r);
 
    for (int i = 0; i < nc; i++)
    {
       jp[i] = c[i];
    }
-   I[r] += nc;
+   UsingBigI() ? bigI[r] += nc : I[r] += nc;
 }
 
 void Table::ShiftUpI()
 {
-   for (int i = size; i > 0; i--)
+   if (!UsingBigI())
    {
-      I[i] = I[i-1];
+      for (int i = size; i > 0; i--)
+      {
+         I[i] = I[i-1];
+      }
+      I[0] = 0;
+   }
+   else
+   {
+      for (int i = size; i > 0; i--)
+      {
+         bigI[i] = bigI[i-1];
+      }
+      bigI[0] = 0;
    }
-   I[0] = 0;
 }
 
 void Table::SetSize(int dim, int connections_per_row)
 {
-   SetDims (dim, dim * connections_per_row);
+   SetDims(dim, bigint(dim) * connections_per_row);
 
    if (size > 0)
    {
-      I[0] = 0;
-      for (int i = 0, j = 0; i < size; i++)
+      if (!UsingBigI())
       {
-         int end = I[i] + connections_per_row;
-         I[i+1] = end;
-         for ( ; j < end; j++) { J[j] = -1; }
+         I[0] = 0;
+         for (int i = 0, j = 0; i < size; i++)
+         {
+            int end = I[i] + connections_per_row;
+            I[i+1] = end;
+            for ( ; j < end; j++) { J[j] = -1; }
+         }
+      }
+      else
+      {
+         bigint j = 0;
+         bigI[0] = 0;
+         for (int i = 0; i < size; i++)
+         {
+            bigint end = bigI[i] + connections_per_row;
+            bigI[i+1] = end;
+            for ( ; j < end; j++) { J[j] = -1; }
+         }
       }
    }
 }
 
-void Table::SetDims(int rows, int nnz)
+void Table::SetDims(int rows, bigint nnz)
 {
-   int j;
-
-   j = (I) ? (I[size]) : (0);
-   if (size != rows)
+   const bool new_use_big_i = (bigint(int(nnz)) != nnz);
+   if (size != rows || new_use_big_i != UsingBigI())
    {
       size = rows;
-      (rows >= 0) ? I.SetSize(rows+1) : I.DeleteAll();
+      if (new_use_big_i != UsingBigI())
+      {
+         UsingBigI() ? bigI.DeleteAll() : I.DeleteAll();
+      }
+      if (size >= 0)
+      {
+         new_use_big_i ? bigI.SetSize(size+1) : I.SetSize(size+1);
+      }
+      else
+      {
+         new_use_big_i ? bigI.DeleteAll() : I.DeleteAll();
+      }
    }
 
-   if (j != nnz)
-   {
-      (nnz > 0) ? J.SetSize(nnz) : J.DeleteAll();
-   }
+   (nnz > 0) ? J.SetSize(nnz) : J.DeleteAll();
 
    if (size >= 0)
    {
-      I[0] = 0;
-      I[size] = nnz;
+      if (!UsingBigI())
+      {
+         I[0] = 0;
+         I[size] = int(nnz);
+      }
+      else
+      {
+         bigI[0] = 0;
+         bigI[size] = nnz;
+      }
    }
 }
 
 int Table::operator()(int i, int j) const
 {
+   MFEM_VERIFY(!UsingBigI(), "");
+
    if ( i>=size || i<0 )
    {
       return -1;
@@ -236,7 +331,8 @@ void Table::GetRow(int i, Array<int> &row) const
                << size << ')');
 
    HostReadJ();
-   HostReadI();
+   if (UsingBigI()) { HostReadBigI(); }
+   else { HostReadI(); }
 
    row.SetSize(RowSize(i));
    row.Assign(GetRow(i));
@@ -244,14 +340,25 @@ void Table::GetRow(int i, Array<int> &row) const
 
 void Table::SortRows()
 {
-   for (int r = 0; r < size; r++)
+   if (!UsingBigI())
    {
-      std::sort(J + I[r], J + I[r+1]);
+      for (int r = 0; r < size; r++)
+      {
+         std::sort(J.GetData()+I[r], J.GetData()+I[r+1]);
+      }
+   }
+   else
+   {
+      for (int r = 0; r < size; r++)
+      {
+         std::sort(J.GetData()+bigI[r], J.GetData()+bigI[r+1]);
+      }
    }
 }
 
 void Table::SetIJ(int *newI, int *newJ, int newsize)
 {
+   if (UsingBigI()) { bigI.DeleteAll(); }
    if (newsize >= 0)
    {
       size = newsize;
@@ -262,6 +369,8 @@ void Table::SetIJ(int *newI, int *newJ, int newsize)
 
 int Table::Push(int i, int j)
 {
+   MFEM_VERIFY(!UsingBigI(), "");
+
    MFEM_ASSERT(i >=0 &&
                i<size, "Index out of bounds.  i = " << i << " size " << size);
 
@@ -286,6 +395,8 @@ int Table::Push(int i, int j)
 
 void Table::Finalize()
 {
+   MFEM_VERIFY(!UsingBigI(), "");
+
    int i, j, end, sum = 0, n = 0, newI = 0;
 
    for (i=0; i<I[size]; i++)
@@ -324,14 +435,16 @@ void Table::MakeFromList(int nrows, const Array<Connection> &list)
    Clear();
 
    size = nrows;
-   int nnz = list.Size();
+   const bigint nnz = list.Size();
+   const bool use_big_i = (bigint(int(nnz)) != nnz);
 
-   I.SetSize(size+1);
+   use_big_i ? bigI.SetSize(size+1) : I.SetSize(size+1);
    J.SetSize(nnz);
 
-   for (int i = 0, k = 0; i <= size; i++)
+   bigint k = 0;
+   for (int i = 0; i <= size; i++)
    {
-      I[i] = k;
+      use_big_i ? bigI[i] = k : I[i] = int(k);
       while (k < nnz && list[k].from == i)
       {
          J[k] = list[k].to;
@@ -342,30 +455,25 @@ void Table::MakeFromList(int nrows, const Array<Connection> &list)
 
 int Table::Width() const
 {
-   int width = -1, nnz = (size >= 0) ? I[size] : 0;
-   for (int k = 0; k < nnz; k++)
-   {
-      if (J[k] > width) { width = J[k]; }
-   }
-   return width + 1;
+   return (J.Size() > 0) ? J.Max() + 1 : 0;
 }
 
 void Table::Print(std::ostream & os, int width) const
 {
-   int i, j;
-
-   for (i = 0; i < size; i++)
+   for (int i = 0; i < size; i++)
    {
       os << "[row " << i << "]\n";
-      for (j = I[i]; j < I[i+1]; j++)
+      const int row_size = RowSize(i);
+      const int *row = GetRow(i);
+      for (int j = 0; j < row_size; j++)
       {
-         os << setw(5) << J[j];
-         if ( !((j+1-I[i]) % width) )
+         os << setw(5) << row[j];
+         if ( !((j+1) % width) )
          {
             os << '\n';
          }
       }
-      if ((j-I[i]) % width)
+      if (row_size % width)
       {
          os << '\n';
       }
@@ -378,9 +486,11 @@ void Table::PrintMatlab(std::ostream & os) const
 
    for (i = 0; i < size; i++)
    {
-      for (j = I[i]; j < I[i+1]; j++)
+      const int row_size = RowSize(i);
+      const int *row = GetRow(i);
+      for (j = 0; j < row_size; j++)
       {
-         os << i << " " << J[j] << " 1. \n";
+         os << i << " " << row[j] << " 1. \n";
       }
    }
 
@@ -391,35 +501,51 @@ void Table::Save(std::ostream &os) const
 {
    os << size << '\n';
 
-   for (int i = 0; i <= size; i++)
+   if (!UsingBigI())
    {
-      os << I[i] << '\n';
+      I.Save(os, 1);
    }
-   for (int i = 0, nnz = I[size]; i < nnz; i++)
+   else
    {
-      os << J[i] << '\n';
+      bigI.Save(os, 1);
    }
+   J.Save(os, 1);
 }
 
 void Table::Load(std::istream &in)
 {
+   Clear();
+
    in >> size;
    I.SetSize(size+1);
    for (int i = 0; i <= size; i++)
    {
-      in >> I[i];
-   }
-   int nnz = I[size];
-   J.SetSize(nnz);
-   for (int j = 0; j < nnz; j++)
-   {
-      in >> J[j];
+      bigint big_offset;
+      in >> big_offset;
+      const int offset = int(big_offset);
+      if (bigint(offset) != big_offset)
+      {
+         // switch to using bigI instead of I
+         bigI.SetSize(size+1);
+         for (int j = 0; j < i; j++) { bigI[j] = I[j]; }
+         I.DeleteAll();
+         bigI[i] = big_offset;
+         for (i++; i <= size; i++)
+         {
+            in >> bigI[i];
+         }
+         break;
+      }
+      I[i] = offset;
    }
+   J.SetSize(UsingBigI() ? bigI[size] : I[size]);
+   J.Load(in, 1);
 }
 
 void Table::Clear()
 {
    I.DeleteAll();
+   bigI.DeleteAll();
    J.DeleteAll();
    size = -1;
 }
@@ -436,50 +562,72 @@ void Table::Swap(Table & other)
 
 std::size_t Table::MemoryUsage() const
 {
-   if (size < 0 || I == NULL) { return 0; }
-   return (size+1 + I[size]) * sizeof(int);
+   return I.MemoryUsage() + bigI.MemoryUsage() + J.MemoryUsage();
 }
 
-void Transpose(const Table &A, Table &At, int ncols_A_)
+template <typename TI, typename TJ>
+void TransposeImpl(const TI *i_A, const TJ *j_A, const TJ nrows_A,
+                   const TJ ncols_A, const TI nnz_A, TI *i_At, TJ *j_At)
 {
-   const int *i_A     = A.GetI();
-   const int *j_A     = A.GetJ();
-   const int  nrows_A = A.Size();
-   const int  ncols_A = (ncols_A_ < 0) ? A.Width() : ncols_A_;
-   const int  nnz_A   = i_A[nrows_A];
-
-   At.SetDims (ncols_A, nnz_A);
-
-   int *i_At = At.GetI();
-   int *j_At = At.GetJ();
-
-   for (int i = 0; i <= ncols_A; i++)
+   for (TJ i = 0; i <= ncols_A; i++)
    {
       i_At[i] = 0;
    }
-   for (int i = 0; i < nnz_A; i++)
+   for (TI i = 0; i < nnz_A; i++)
    {
       i_At[j_A[i]+1]++;
    }
-   for (int i = 1; i < ncols_A; i++)
+   for (TJ i = 1; i < ncols_A; i++)
    {
       i_At[i+1] += i_At[i];
    }
 
-   for (int i = 0; i < nrows_A; i++)
+   for (TJ i = 0; i < nrows_A; i++)
    {
-      for (int j = i_A[i]; j < i_A[i+1]; j++)
+      for (TI j = i_A[i]; j < i_A[i+1]; j++)
       {
          j_At[i_At[j_A[j]]++] = i;
       }
    }
-   for (int i = ncols_A; i > 0; i--)
+   for (TJ i = ncols_A; i > 0; i--)
    {
       i_At[i] = i_At[i-1];
    }
    i_At[0] = 0;
 }
 
+void Transpose(const Table &A, Table &At, int ncols_A_)
+{
+   const int *j_A     = A.HostReadJ();
+   const int  nrows_A = A.Size();
+   const int  ncols_A = (ncols_A_ < 0) ? A.Width() : ncols_A_;
+
+   if (!A.UsingBigI())
+   {
+      const int *i_A   = A.HostReadI();
+      const int  nnz_A = i_A[nrows_A];
+
+      At.SetDims(ncols_A, nnz_A);
+
+      int *i_At = At.HostWriteI();
+      int *j_At = At.HostWriteJ();
+
+      TransposeImpl(i_A, j_A, nrows_A, ncols_A, nnz_A, i_At, j_At);
+   }
+   else
+   {
+      const bigint *i_A   = A.HostReadBigI();
+      const bigint  nnz_A = i_A[nrows_A];
+
+      At.SetDims(ncols_A, nnz_A);
+
+      bigint *i_At = At.HostWriteBigI();
+      int    *j_At = At.HostWriteJ();
+
+      TransposeImpl(i_A, j_A, nrows_A, ncols_A, nnz_A, i_At, j_At);
+   }
+}
+
 Table * Transpose(const Table &A)
 {
    Table * At = new Table;
@@ -504,6 +652,8 @@ void Transpose(const Array<int> &A, Table &At, int ncols_A_)
 
 void Mult(const Table &A, const Table &B, Table &C)
 {
+   MFEM_VERIFY(!A.UsingBigI() && !B.UsingBigI(), "");
+
    int  i, j, k, l, m;
    const int *i_A     = A.GetI();
    const int *j_A     = A.GetJ();
@@ -641,6 +791,8 @@ int DSTable::Push_(int r, int c)
    n->Index  = NumEntries;
    n->Prev   = Rows[r];
    Rows[r]   = n;
+   MFEM_VERIFY(NumEntries != std::numeric_limits<int>::max(),
+               "integer overflow error");
    return (NumEntries++);
 }
 

general/table.hpp

@@ -42,6 +42,7 @@ struct Connection
 class Table
 {
 protected:
+   // FIXME: this member can mess up the default move ctor ?!!
    int size; ///< The number of TYPE I elements.
 
    /// @name Arrays for the connectivity information in the CSR storage.
@@ -50,6 +51,10 @@ protected:
    /// The length of the I array is 'size + 1',
    Array<int> I;
 
+   /** @brief Alternative to the I array. Used when the number of connections
+       overflows the int type. */
+   Array<bigint> bigI;
+
    /// @brief The length of the J array is equal to the number of connections
    /// between TYPE I and TYPE II elements.
    Array<int> J;
@@ -83,21 +88,25 @@ public:
    /// @name Used together with the default constructor
    /// @{
    void MakeI(int nrows);
-   void AddAColumnInRow(int r) { I[r]++; }
-   void AddColumnsInRow(int r, int ncol) { I[r] += ncol; }
+   void AddAColumnInRow(int r) { UsingBigI() ? bigI[r]++ : I[r]++; }
+   void AddColumnsInRow (int r, int ncol)
+   { UsingBigI() ? bigI[r] += ncol : I[r] += ncol; }
    void MakeJ();
-   void AddConnection(int r, int c) { J[I[r]++] = c; }
+   void AddConnection (int r, int c)
+   { UsingBigI() ? J[bigI[r]++] = c : J[I[r]++] = c; }
    void AddConnections(int r, const int *c, int nc);
    void ShiftUpI();
    /// @}
 
+   bool UsingBigI() const { return !bigI.IsEmpty(); }
+
    /// Set the size and the number of connections for the table.
    void SetSize(int dim, int connections_per_row);
 
    /// @brief Set the rows and the number of all connections for the table.
    ///
    /// Does NOT initialize the whole array I ! (I[0]=0 and I[rows]=nnz only)
-   void SetDims(int rows, int nnz);
+   void SetDims(int rows, bigint nnz);
 
    /// Returns the number of TYPE I elements.
    inline int Size() const { return size; }
@@ -107,7 +116,7 @@ public:
    /// If Finalize() is not called, it returns the number of possible
    /// connections established by the used constructor. Otherwise, it is exactly
    /// the number of established connections after calling Finalize(). */
-   inline int Size_of_connections() const { return J.Size(); }
+   inline bigint Size_of_connections() const { return J.Size(); }
 
    /// @brief Returns index of the connection between element i of TYPE I and
    /// element j of TYPE II.
@@ -119,27 +128,94 @@ public:
    /// Return row i in array row (the Table must be finalized)
    void GetRow(int i, Array<int> &row) const;
 
-   int RowSize(int i) const { return I[i+1] - I[i]; }
+   int RowSize(int i) const
+   { return UsingBigI() ? int(bigI[i+1]-bigI[i]): I[i+1]-I[i]; }
 
-   const int *GetRow(int i) const { return J.GetMemory() + I[i]; }
-   int *GetRow(int i) { return J.GetMemory() + I[i]; }
+   const int *GetRow(int i) const
+   { return UsingBigI() ? J.GetData()+bigI[i] : J.GetData()+I[i]; }
+
+   int *GetRow(int i)
+   { return UsingBigI() ? J.GetData()+bigI[i] : J.GetData()+I[i]; }
+
+   int *GetI()
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.GetData();
+   }
 
-   int *GetI() { return I.GetData(); }
    int *GetJ() { return J.GetData(); }
-   const int *GetI() const { return I.GetData(); }
+
+   const int *GetI() const
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.GetData();
+   }
+
    const int *GetJ() const { return J.GetData(); }
 
-   Memory<int> &GetIMemory() { return I.GetMemory(); }
+   Memory<int> &GetIMemory()
+   { MFEM_ASSERT(!UsingBigI(), ""); return I.GetMemory(); }
+
    Memory<int> &GetJMemory() { return J.GetMemory(); }
-   const Memory<int> &GetIMemory() const { return I.GetMemory(); }
+
+   const Memory<int> &GetIMemory() const
+   { MFEM_ASSERT(!UsingBigI(), ""); return I.GetMemory(); }
+
    const Memory<int> &GetJMemory() const { return J.GetMemory(); }
 
-   const int *ReadI(bool on_dev = true) const { return I.Read(on_dev); }
-   int *WriteI(bool on_dev = true) { return I.Write(on_dev); }
-   int *ReadWriteI(bool on_dev = true) { return I.ReadWrite(on_dev); }
-   const int *HostReadI() const { return I.HostRead(); }
-   int *HostWriteI() { return I.HostWrite(); }
-   int *HostReadWriteI() { return I.HostReadWrite(); }
+   const int *ReadI(bool on_dev = true) const
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.Read(on_dev);
+   }
+
+   int *WriteI(bool on_dev = true)
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.Write(on_dev);
+   }
+
+   int *ReadWriteI(bool on_dev = true)
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.ReadWrite(on_dev);
+   }
+
+   const int *HostReadI() const
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.HostRead();
+   }
+
+   int *HostWriteI()
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.HostWrite();
+   }
+
+   int *HostReadWriteI()
+   {
+      MFEM_ASSERT(!UsingBigI(), "");
+      return I.HostReadWrite();
+   }
+
+   const bigint *HostReadBigI() const
+   {
+      MFEM_ASSERT(UsingBigI(), "");
+      return bigI.HostRead();
+   }
+
+   bigint *HostWriteBigI()
+   {
+      MFEM_ASSERT(UsingBigI(), "");
+      return bigI.HostWrite();
+   }
+
+   bigint *HostReadWriteBigI()
+   {
+      MFEM_ASSERT(UsingBigI(), "");
+      return bigI.HostReadWrite();
+   }
 
    const int *ReadJ(bool on_dev = true) const { return J.Read(on_dev); }
    int *WriteJ(bool on_dev = true) { return J.Write(on_dev); }
@@ -181,7 +257,7 @@ public:
    int Width() const;
 
    /// Releases ownership of and null-ifies the data.
-   void LoseData() { size = -1; I.LoseData(); J.LoseData(); }
+   void LoseData() { size = -1; I.LoseData(); bigI.LoseData(); J.LoseData(); }
 
    /// Prints the table to the stream @a out.
    void Print(std::ostream & out = mfem::out, int width = 4) const;

linalg/CMakeLists.txt

@@ -55,7 +55,6 @@ list(APPEND HDRS
   dinvariants.hpp
   dtensor.hpp
   dual.hpp
-  eigensolver.hpp
   filteredsolver.hpp
   handle.hpp
   invariants.hpp

linalg/densemat.hpp

@@ -82,6 +82,15 @@ public:
       }
    }
 
+   /// Make the DenseMatrix to reference the given sub-Memory of @a base.
+   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
+       not delete the @a base Memory. */
+   void MakeRef(Memory<real_t> &base, int offset, int h, int w)
+   {
+      data.MakeRef(base, offset, h*w);
+      height = h; width = w;
+   }
+
    /// Change the data array and the size of the DenseMatrix.
    /** The DenseMatrix does not assume ownership of the data array, i.e. it will
        not delete the data array @a d. */

linalg/dtensor.hpp

@@ -24,13 +24,14 @@ class TensorInd
 {
 public:
    MFEM_HOST_DEVICE
-   static inline int result(const int* sizes, T first, Args... args)
+   static inline bigint result(const int* sizes, T first, Args... args)
    {
 #if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
       MFEM_ASSERT(first<sizes[N-1],"Trying to access out of boundary.");
 #endif
-      return static_cast<int>(first + sizes[N - 1] * TensorInd < N + 1, Dim, Args... >
-                              ::result(sizes, args...));
+      return static_cast<bigint>(
+                first + sizes[N - 1] * TensorInd < N + 1, Dim, Args... >
+                ::result(sizes, args...));
    }
 };
 
@@ -40,13 +41,13 @@ class TensorInd<Dim, Dim, T, Args...>
 {
 public:
    MFEM_HOST_DEVICE
-   static inline int result(const int* sizes, T first, Args... args)
+   static inline bigint result(const int* sizes, T first, Args... args)
    {
 #if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
       MFEM_ASSERT(first<static_cast<T>(sizes[Dim-1]),
                   "Trying to access out of boundary.");
 #endif
-      return static_cast<int>(first);
+      return static_cast<bigint>(first);
    }
 };
 
@@ -57,7 +58,7 @@ class Init
 {
 public:
    MFEM_HOST_DEVICE
-   static inline int result(int* sizes, T first, Args... args)
+   static inline bigint result(int* sizes, T first, Args... args)
    {
       sizes[N - 1] = first;
       return first * Init < N + 1, Dim, Args... >::result(sizes, args...);
@@ -70,10 +71,10 @@ class Init<Dim, Dim, T, Args...>
 {
 public:
    MFEM_HOST_DEVICE
-   static inline int result(int* sizes, T first, Args... args)
+   static inline bigint result(int* sizes, T first, Args... args)
    {
       sizes[Dim - 1] = first;
-      return first;
+      return static_cast<bigint>(first);
    }
 };
 
@@ -83,7 +84,7 @@ template<int Dim, typename Scalar = real_t>
 class DeviceTensor
 {
 protected:
-   int capacity;
+   bigint capacity;
    Scalar *data;
    int sizes[Dim];
 
@@ -99,8 +100,7 @@ public:
    {
       static_assert(sizeof...(args) == Dim, "Wrong number of arguments");
       // Initialize sizes, and compute the number of values
-      const long int nb = Init<1, Dim, Args...>::result(sizes, args...);
-      capacity = nb;
+      capacity = Init<1, Dim, Args...>::result(sizes, args...);
       data = (capacity > 0) ? data_ : nullptr;
    }
 
@@ -122,7 +122,7 @@ public:
    }
 
    /// Subscript operator where the tensor is viewed as a 1D array.
-   MFEM_HOST_DEVICE inline Scalar& operator[](int i) const
+   MFEM_HOST_DEVICE inline Scalar& operator[](bigint i) const
    {
       return data[i];
    }

linalg/eigensolver.hpp

@@ -1,203 +0,0 @@
-// 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.
-
-/**
- * @file eigensolver.hpp
- *
- * @brief This file contains a common interface for all eigensolver classes
- */
-
-#ifndef MFEM_EIGENSOLVER
-#define MFEM_EIGENSOLVER
-
-#ifdef MFEM_HYPRE
-#include "hypre.hpp"
-#endif
-
-#ifdef MFEM_SLEPC
-#include "slepc.hpp"
-#endif
-
-namespace mfem
-{
-
-enum class EigenSolverType
-{
-   HYPRE,
-   SLEPC,
-   INVALID_TYPE
-};
-
-/// Provides base class for MFEM Eigensolvers
-class EigenSolverBase
-{
-public:
-   EigenSolverBase() {}
-
-   /// Destructor
-   virtual ~EigenSolverBase() = default;
-
-   /// Solves the eigenvalue problem
-   virtual void Solve() = 0;
-
-   /// Set the required number of modes
-   virtual void SetNumModes(int num_Modes)
-   {
-      numModes=num_Modes;
-   }
-
-   /// @brief  Set the operator to the eigenvalue problem
-   /// @param A - operator
-   virtual void SetOperator(Operator& A) = 0;
-
-   /// @brief Sets operators for the generalized eigenvalue problem
-   /// @param A - operator
-   /// @param M - mass matrix
-   virtual void SetOperator(Operator& A, Operator& M)
-   {
-      MFEM_ABORT("Generalized eigensolver is not supported!");
-   }
-
-   /// Optional method - sets preconditioner for the
-   /// eigenvalue solver.
-   virtual void SetPreconditioner(Solver& precond)
-   {
-      MFEM_ABORT("Preconditioner is not supported!");
-   }
-
-   /// Returns the converged eigenvalues
-   virtual void GetEigenvalues(Array<real_t>& eigen_vals) = 0;
-
-   /// Returns the vec_index eigenvector.
-   virtual void GetEigenvector(int vec_index, Vector& vector) = 0;
-
-   /// Returns the eigensolver type.
-   EigenSolverType GetSolverType() { return eigSolverType; }
-
-protected:
-   int numModes = 0;
-   EigenSolverType eigSolverType = EigenSolverType::INVALID_TYPE;
-};
-
-#ifdef MFEM_HYPRE
-class EigenSolverHypreLOBPCG : public EigenSolverBase
-{
-public:
-   EigenSolverHypreLOBPCG(MPI_Comm comm)
-   {
-      eigenSolver = std::make_unique<HypreLOBPCG>(comm);
-      eigSolverType = EigenSolverType::HYPRE;
-   }
-
-   ~EigenSolverHypreLOBPCG() {}
-
-   void Solve() override { eigenSolver->Solve(); }
-   void SetNumModes(int num_Modes) override
-   {
-      eigenSolver->SetNumModes(num_Modes);
-      numModes = num_Modes;
-   }
-
-   void SetOperator(Operator& A) override { eigenSolver->SetOperator(A); }
-
-   void SetOperator(Operator& A, Operator& M) override
-   {
-      eigenSolver->SetOperator(A);
-      eigenSolver->SetMassMatrix(M);
-   }
-
-   void SetPreconditioner(Solver& precond) override { eigenSolver->SetPreconditioner(precond); }
-   void GetEigenvalues(Array<real_t>& eigen_vals) override { eigenSolver->GetEigenvalues(eigen_vals); }
-   void GetEigenvector(int vec_index, Vector& vector) override
-   {
-      const HypreParVector& eigenvec = eigenSolver->GetEigenvector(vec_index);
-      vector = eigenvec;
-   }
-
-   void SetTol(real_t tol) { eigenSolver->SetTol(tol); }
-   void SetRelTol(real_t rel_tol) { eigenSolver->SetRelTol(rel_tol); }
-   void SetMaxIter(int max_iter) { eigenSolver->SetMaxIter(max_iter); }
-   void SetPrintLevel(int logging) { eigenSolver->SetPrintLevel(logging); }
-   void SetRandomSeed(int seed) { eigenSolver->SetRandomSeed(seed); }
-   void SetPrecondUsageMode(int usage_mode) { eigenSolver->SetPrecondUsageMode(usage_mode); }
-
-private:
-   std::unique_ptr<HypreLOBPCG> eigenSolver = nullptr;
-};
-#endif
-
-#ifdef MFEM_SLEPC
-class EigenSolverSlepc : public EigenSolverBase
-{
-public:
-   EigenSolverSlepc(MPI_Comm comm)
-   {
-      eigSolverType = EigenSolverType::SLEPC;
-      eigenSolver = std::make_unique<SlepcEigenSolver>(comm);
-
-      eigenSolver->SetWhichEigenpairs(SlepcEigenSolver::TARGET_REAL);
-      eigenSolver->SetTarget(0.0);
-      eigenSolver->SetSpectralTransformation(SlepcEigenSolver::SHIFT_INVERT);
-   }
-
-   ~EigenSolverSlepc() {}
-
-   void Solve() override { eigenSolver->Solve(); }
-   void SetNumModes(int num_Modes) override
-   {
-      eigenSolver->SetNumModes(num_Modes);
-      numModes = num_Modes;
-   }
-   /// @brief  Set the operator to the slepc eigenvalue problem. This method deep copies data to create a PetscParMatrix
-   /// @param A - operator, must be of type HypreParMatrix.
-   void SetOperator(Operator& A) override
-   {
-      petscMatA = std::make_unique<PetscParMatrix>
-                  (dynamic_cast<HypreParMatrix*>(&A));
-      eigenSolver->SetOperator(*petscMatA);
-   }
-   /// @brief  Set the operators to the slepc eigenvalue problem. This method deep copies data to create a PetscParMatrix
-   /// @param A - operator, must be of type HypreParMatrix.
-   /// @param M - operator, must be of type HypreParMatrix.
-   void SetOperator(Operator& A, Operator& M) override
-   {
-      petscMatA = std::make_unique<PetscParMatrix>
-                  (dynamic_cast<const HypreParMatrix*>(&A));
-      petscMatM = std::make_unique<PetscParMatrix>
-                  (dynamic_cast<const HypreParMatrix*>(&M));
-
-      eigenSolver->SetOperators(*petscMatA, *petscMatM);
-   }
-   void SetPreconditioner([[maybe_unused]] Solver& precond) override {}
-   void GetEigenvalues(Array<real_t>& eigen_vals) override
-   {
-      eigen_vals.SetSize(numModes);
-      for (int ik = 0; ik < numModes; ik++)
-      {
-         eigenSolver->GetEigenvalue(static_cast<unsigned int>(ik), eigen_vals[ik]);
-      }
-   }
-   void GetEigenvector( int vec_index, Vector& vector) override
-   { eigenSolver->GetEigenvector(vec_index, vector); }
-
-   void SetTol(real_t tol) { eigenSolver->SetTol(tol); }
-   void SetMaxIter(int max_iter) { eigenSolver->SetMaxIter(max_iter); }
-
-private:
-   std::unique_ptr<SlepcEigenSolver> eigenSolver = nullptr;
-   std::unique_ptr<PetscParMatrix> petscMatA = nullptr;
-   std::unique_ptr<PetscParMatrix> petscMatM = nullptr;
-};
-#endif
-
-}  // namespace mfem
-
-#endif

linalg/particlevector.cpp

@@ -16,7 +16,8 @@ namespace mfem
 
 void ParticleVector::GrowSize(int min_num_vectors, bool keep_data)
 {
-   const int nsize = std::max(min_num_vectors*vdim, 2 * data.Capacity());
+   const bigint nsize = std::max(bigint(min_num_vectors)*vdim,
+                                 2 * data.Capacity());
    Memory<real_t> p(nsize, data.GetMemoryType());
    if (keep_data) { p.CopyFrom(data, size); }
    p.UseDevice(data.UseDevice());

linalg/vector.cpp

@@ -107,7 +107,7 @@ static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
 
 Vector::Vector(const Vector &v)
 {
-   const int s = v.Size();
+   const bigint s = v.Size();
    size = s;
    if (s > 0)
    {
@@ -126,10 +126,8 @@ Vector::Vector(Vector &&v)
 
 void Vector::Load(std::istream **in, int np, int *dim)
 {
-   int i, j, s;
-
-   s = 0;
-   for (i = 0; i < np; i++)
+   bigint s = 0;
+   for (int i = 0; i < np; i++)
    {
       s += dim[i];
    }
@@ -137,10 +135,10 @@ void Vector::Load(std::istream **in, int np, int *dim)
    SetSize(s);
    HostWrite();
 
-   int p = 0;
-   for (i = 0; i < np; i++)
+   bigint p = 0;
+   for (int i = 0; i < np; i++)
    {
-      for (j = 0; j < dim[i]; j++)
+      for (int j = 0; j < dim[i]; j++)
       {
          *in[i] >> data[p++];
          // Clang's libc++ sets the failbit when (correctly) parsing subnormals,
@@ -153,12 +151,12 @@ void Vector::Load(std::istream **in, int np, int *dim)
    }
 }
 
-void Vector::Load(std::istream &in, int Size)
+void Vector::Load(std::istream &in, bigint Size)
 {
    SetSize(Size);
    HostWrite();
 
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       in >> data[i];
       // Clang's libc++ sets the failbit when (correctly) parsing subnormals,
@@ -170,12 +168,12 @@ void Vector::Load(std::istream &in, int Size)
    }
 }
 
-real_t &Vector::Elem(int i)
+real_t &Vector::Elem(bigint i)
 {
    return operator()(i);
 }
 
-const real_t &Vector::Elem(int i) const
+const real_t &Vector::Elem(bigint i) const
 {
    return operator()(i);
 }
@@ -187,7 +185,7 @@ real_t Vector::operator*(const real_t *v) const
 #ifdef MFEM_USE_LEGACY_OPENMP
    #pragma omp parallel for reduction(+:dot)
 #endif
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       dot += data[i] * v[i];
    }
@@ -232,18 +230,22 @@ Vector &Vector::operator=(Vector &&v)
 Vector &Vector::operator=(real_t value)
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = value; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] = value;
+   });
    return *this;
 }
 
 Vector &Vector::operator*=(real_t c)
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] *= c; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] *= c;
+   });
    return *this;
 }
 
@@ -252,20 +254,24 @@ Vector &Vector::operator*=(const Vector &v)
    MFEM_ASSERT(size == v.size, "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
    const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] *= x[i]; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] *= x[i];
+   });
    return *this;
 }
 
 Vector &Vector::operator/=(real_t c)
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    const real_t m = 1.0/c;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] *= m; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] *= m;
+   });
    return *this;
 }
 
@@ -274,19 +280,23 @@ Vector &Vector::operator/=(const Vector &v)
    MFEM_ASSERT(size == v.size, "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
    const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] /= x[i]; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] /= x[i];
+   });
    return *this;
 }
 
 Vector &Vector::operator-=(real_t c)
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] -= c; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] -= c;
+   });
    return *this;
 }
 
@@ -295,19 +305,23 @@ Vector &Vector::operator-=(const Vector &v)
    MFEM_ASSERT(size == v.size, "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
    const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] -= x[i]; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] -= x[i];
+   });
    return *this;
 }
 
 Vector &Vector::operator+=(real_t c)
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += c; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] += c;
+   });
    return *this;
 }
 
@@ -316,10 +330,12 @@ Vector &Vector::operator+=(const Vector &v)
    MFEM_ASSERT(size == v.size, "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
    const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += x[i]; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] += x[i];
+   });
    return *this;
 }
 
@@ -329,11 +345,13 @@ Vector &Vector::Add(const real_t a, const Vector &Va)
 
    if (a != 0.0)
    {
-      const int N = size;
       const bool use_dev = UseDevice() || Va.UseDevice();
       const auto x = Va.Read(use_dev);
       auto y = ReadWrite(use_dev);
-      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += a * x[i]; });
+      mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+      {
+         y[i] += a * x[i];
+      });
    }
    return *this;
 }
@@ -343,58 +361,69 @@ Vector &Vector::Set(const real_t a, const Vector &Va)
    MFEM_ASSERT(size == Va.size, "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || Va.UseDevice();
-   const int N = size;
    const auto x = Va.Read(use_dev);
    auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = a * x[i]; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] = a * x[i];
+   });
    return *this;
 }
 
-void Vector::SetVector(const Vector &v, int offset)
+void Vector::SetVector(const Vector &v, bigint offset)
 {
    MFEM_ASSERT(v.Size() + offset <= size, "invalid sub-vector");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   const int vs = v.Size();
+   const bigint vs = v.Size();
    const auto vp = v.Read(use_dev);
    // Use read+write access for *this - we only modify some of its entries
    auto p = ReadWrite(use_dev) + offset;
-   mfem::forall_switch(use_dev, vs, [=] MFEM_HOST_DEVICE (int i) { p[i] = vp[i]; });
+   mfem::forall_switch(use_dev, vs, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      p[i] = vp[i];
+   });
 }
 
-void Vector::AddSubVector(const Vector &v, int offset)
+void Vector::AddSubVector(const Vector &v, bigint offset)
 {
    MFEM_ASSERT(v.Size() + offset <= size, "invalid sub-vector");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   const int vs = v.Size();
+   const bigint vs = v.Size();
    const auto vp = v.Read(use_dev);
    auto p = ReadWrite(use_dev) + offset;
-   mfem::forall_switch(use_dev, vs, [=] MFEM_HOST_DEVICE (int i) { p[i] += vp[i]; });
+   mfem::forall_switch(use_dev, vs, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      p[i] += vp[i];
+   });
 }
 
 void Vector::Neg()
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = -y[i]; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] = -y[i];
+   });
 }
 
 void Vector::Reciprocal()
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = 1.0/y[i]; });
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] = 1.0/y[i];
+   });
 }
 
 void Vector::Abs()
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
    {
       y[i] = std::abs(y[i]);
    });
@@ -403,9 +432,8 @@ void Vector::Abs()
 void Vector::Pow(const real_t p)
 {
    const bool use_dev = UseDevice();
-   const int N = size;
    auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (bigint i)
    {
       y[i] = std::pow(y[i], p);
    });
@@ -418,15 +446,18 @@ void add(const Vector &v1, const Vector &v2, Vector &v)
 
 #if !defined(MFEM_USE_LEGACY_OPENMP)
    const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
-   const int N = v.size;
+   const bigint N = v.size;
    // Note: get read access first, in case v is the same as v1/v2.
    const auto x1 = v1.Read(use_dev);
    const auto x2 = v2.Read(use_dev);
    auto y = v.Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = x1[i] + x2[i]; });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      y[i] = x1[i] + x2[i];
+   });
 #else
    #pragma omp parallel for
-   for (int i = 0; i < v.size; i++)
+   for (bigint i = 0; i < v.size; i++)
    {
       v.data[i] = v1.data[i] + v2.data[i];
    }
@@ -450,21 +481,21 @@ void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v)
    {
 #if !defined(MFEM_USE_LEGACY_OPENMP)
       const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
-      const int N = v.size;
+      const bigint N = v.size;
       // Note: get read access first, in case v is the same as v1/v2.
       const auto d_x = v1.Read(use_dev);
       const auto d_y = v2.Read(use_dev);
       auto d_z = v.Write(use_dev);
-      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
       {
          d_z[i] = d_x[i] + alpha * d_y[i];
       });
 #else
       const real_t *v1p = v1.data, *v2p = v2.data;
       real_t *vp = v.data;
-      const int s = v.size;
+      const bigint s = v.size;
       #pragma omp parallel for
-      for (int i = 0; i < s; i++)
+      for (bigint i = 0; i < s; i++)
       {
          vp[i] = v1p[i] + alpha*v2p[i];
       }
@@ -489,12 +520,12 @@ void add(const real_t a, const Vector &x, const Vector &y, Vector &z)
    {
 #if !defined(MFEM_USE_LEGACY_OPENMP)
       const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
-      const int N = x.size;
+      const bigint N = x.size;
       // Note: get read access first, in case z is the same as x/y.
       const auto xd = x.Read(use_dev);
       const auto yd = y.Read(use_dev);
       auto zd = z.Write(use_dev);
-      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
       {
          zd[i] = a * (xd[i] + yd[i]);
       });
@@ -502,9 +533,9 @@ void add(const real_t a, const Vector &x, const Vector &y, Vector &z)
       const real_t *xp = x.data;
       const real_t *yp = y.data;
       real_t       *zp = z.data;
-      const int      s = x.size;
+      const bigint   s = x.size;
       #pragma omp parallel for
-      for (int i = 0; i < s; i++)
+      for (bigint i = 0; i < s; i++)
       {
          zp[i] = a * (xp[i] + yp[i]);
       }
@@ -544,12 +575,12 @@ void add(const real_t a, const Vector &x,
    {
 #if !defined(MFEM_USE_LEGACY_OPENMP)
       const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
-      const int N = x.size;
+      const bigint N = x.size;
       // Note: get read access first, in case z is the same as x/y.
       const auto xd = x.Read(use_dev);
       const auto yd = y.Read(use_dev);
       auto zd = z.Write(use_dev);
-      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
       {
          zd[i] = a * xd[i] + b * yd[i];
       });
@@ -557,9 +588,9 @@ void add(const real_t a, const Vector &x,
       const real_t *xp = x.data;
       const real_t *yp = y.data;
       real_t       *zp = z.data;
-      const int      s = x.size;
+      const bigint   s = x.size;
       #pragma omp parallel for
-      for (int i = 0; i < s; i++)
+      for (bigint i = 0; i < s; i++)
       {
          zp[i] = a * xp[i] + b * yp[i];
       }
@@ -574,12 +605,12 @@ void subtract(const Vector &x, const Vector &y, Vector &z)
 
 #if !defined(MFEM_USE_LEGACY_OPENMP)
    const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
-   const int N = x.size;
+   const bigint N = x.size;
    // Note: get read access first, in case z is the same as x/y.
    const auto xd = x.Read(use_dev);
    const auto yd = y.Read(use_dev);
    auto zd = z.Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
    {
       zd[i] = xd[i] - yd[i];
    });
@@ -587,9 +618,9 @@ void subtract(const Vector &x, const Vector &y, Vector &z)
    const real_t *xp = x.data;
    const real_t *yp = y.data;
    real_t       *zp = z.data;
-   const int     s = x.size;
+   const bigint   s = x.size;
    #pragma omp parallel for
-   for (int i = 0; i < s; i++)
+   for (bigint i = 0; i < s; i++)
    {
       zp[i] = xp[i] - yp[i];
    }
@@ -613,12 +644,12 @@ void subtract(const real_t a, const Vector &x, const Vector &y, Vector &z)
    {
 #if !defined(MFEM_USE_LEGACY_OPENMP)
       const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
-      const int N = x.size;
+      const bigint N = x.size;
       // Note: get read access first, in case z is the same as x/y.
       const auto xd = x.Read(use_dev);
       const auto yd = y.Read(use_dev);
       auto zd = z.Write(use_dev);
-      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
       {
          zd[i] = a * (xd[i] - yd[i]);
       });
@@ -626,9 +657,9 @@ void subtract(const real_t a, const Vector &x, const Vector &y, Vector &z)
       const real_t *xp = x.data;
       const real_t *yp = y.data;
       real_t       *zp = z.data;
-      const int      s = x.size;
+      const bigint   s = x.size;
       #pragma omp parallel for
-      for (int i = 0; i < s; i++)
+      for (bigint i = 0; i < s; i++)
       {
          zp[i] = a * (xp[i] - yp[i]);
       }
@@ -655,12 +686,12 @@ void Vector::median(const Vector &lo, const Vector &hi)
                "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || lo.UseDevice() || hi.UseDevice();
-   const int N = size;
+   const bigint N = size;
    // Note: get read access first, in case *this is the same as lo/hi.
    const auto l = lo.Read(use_dev);
    const auto h = hi.Read(use_dev);
    auto m = Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
    {
       if (m[i] < l[i])
       {
@@ -675,13 +706,13 @@ void Vector::median(const Vector &lo, const Vector &hi)
 
 void Vector::GetSubVector(const Array<int> &dofs, Vector &elemvect) const
 {
-   const int n = dofs.Size();
+   const bigint n = dofs.Size();
    elemvect.SetSize(n);
    const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
    const auto d_X = Read(use_dev);
    const auto d_dofs = dofs.Read(use_dev);
    auto d_y = elemvect.Write(use_dev);
-   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (bigint i)
    {
       const int dof_i = d_dofs[i];
       d_y[i] = dof_i >= 0 ? d_X[dof_i] : -d_X[-dof_i-1];
@@ -691,8 +722,8 @@ void Vector::GetSubVector(const Array<int> &dofs, Vector &elemvect) const
 void Vector::GetSubVector(const Array<int> &dofs, real_t *elem_data) const
 {
    HostRead();
-   const int n = dofs.Size();
-   for (int i = 0; i < n; i++)
+   const bigint n = dofs.Size();
+   for (bigint i = 0; i < n; i++)
    {
       const int j = dofs[i];
       elem_data[i] = (j >= 0) ? data[j] : -data[-1-j];
@@ -702,11 +733,11 @@ void Vector::GetSubVector(const Array<int> &dofs, real_t *elem_data) const
 void Vector::SetSubVector(const Array<int> &dofs, const real_t value)
 {
    const bool use_dev = UseDevice() || dofs.UseDevice();
-   const int n = dofs.Size();
+   const bigint n = dofs.Size();
    // Use read+write access for *this - we only modify some of its entries
    auto d_X = ReadWrite(use_dev);
    const auto d_dofs = dofs.Read(use_dev);
-   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (bigint i)
    {
       const int j = d_dofs[i];
       if (j >= 0)
@@ -723,7 +754,7 @@ void Vector::SetSubVector(const Array<int> &dofs, const real_t value)
 void Vector::SetSubVectorHost(const Array<int> &dofs, const real_t value)
 {
    HostReadWrite();
-   for (int i = 0; i < dofs.Size(); ++i)
+   for (bigint i = 0; i < dofs.Size(); ++i)
    {
       const int j = dofs[i];
       if (j >= 0)
@@ -744,12 +775,12 @@ void Vector::SetSubVector(const Array<int> &dofs, const Vector &elemvect)
                << ", length of elemvect is " << elemvect.Size());
 
    const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
-   const int n = dofs.Size();
+   const bigint n = dofs.Size();
    // Use read+write access for X - we only modify some of its entries
    auto d_X = ReadWrite(use_dev);
    const auto d_y = elemvect.Read(use_dev);
    const auto d_dofs = dofs.Read(use_dev);
-   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (bigint i)
    {
       const int dof_i = d_dofs[i];
       if (dof_i >= 0)
@@ -767,10 +798,10 @@ void Vector::SetSubVector(const Array<int> &dofs, real_t *elem_data)
 {
    // Use read+write access because we overwrite only part of the data.
    HostReadWrite();
-   const int n = dofs.Size();
-   for (int i = 0; i < n; i++)
+   const bigint n = dofs.Size();
+   for (bigint i = 0; i < n; i++)
    {
-      const int j= dofs[i];
+      const int j = dofs[i];
       if (j >= 0)
       {
          operator()(j) = elem_data[i];
@@ -789,11 +820,11 @@ void Vector::AddElementVector(const Array<int> &dofs, const Vector &elemvect)
                ", length of elemvect is " << elemvect.Size());
 
    const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
-   const int n = dofs.Size();
+   const bigint n = dofs.Size();
    const auto d_y = elemvect.Read(use_dev);
    const auto d_dofs = dofs.Read(use_dev);
    auto d_X = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (bigint i)
    {
       const int j = d_dofs[i];
       if (j >= 0)
@@ -810,8 +841,8 @@ void Vector::AddElementVector(const Array<int> &dofs, const Vector &elemvect)
 void Vector::AddElementVector(const Array<int> &dofs, real_t *elem_data)
 {
    HostReadWrite();
-   const int n = dofs.Size();
-   for (int i = 0; i < n; i++)
+   const bigint n = dofs.Size();
+   for (bigint i = 0; i < n; i++)
    {
       const int j = dofs[i];
       if (j >= 0)
@@ -833,11 +864,11 @@ void Vector::AddElementVector(const Array<int> &dofs, const real_t a,
                ", length of elemvect is " << elemvect.Size());
 
    const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
-   const int n = dofs.Size();
+   const bigint n = dofs.Size();
    const auto d_x = elemvect.Read(use_dev);
    const auto d_dofs = dofs.Read(use_dev);
    auto d_y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (int i)
+   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (bigint i)
    {
       const int j = d_dofs[i];
       if (j >= 0)
@@ -854,17 +885,26 @@ void Vector::AddElementVector(const Array<int> &dofs, const real_t a,
 void Vector::SetSubVectorComplement(const Array<int> &dofs, const real_t val)
 {
    const bool use_dev = UseDevice() || dofs.UseDevice();
-   const int n = dofs.Size();
-   const int N = size;
+   const bigint n = dofs.Size();
+   const bigint N = size;
    Vector dofs_vals(n, use_dev ?
                     Device::GetDeviceMemoryType() :
                     Device::GetHostMemoryType());
    auto d_data = ReadWrite(use_dev);
    auto d_dofs_vals = dofs_vals.Write(use_dev);
    const auto d_dofs = dofs.Read(use_dev);
-   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (int i) { d_dofs_vals[i] = d_data[d_dofs[i]]; });
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { d_data[i] = val; });
-   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (int i) { d_data[d_dofs[i]] = d_dofs_vals[i]; });
+   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      d_dofs_vals[i] = d_data[d_dofs[i]];
+   });
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      d_data[i] = val;
+   });
+   mfem::forall_switch(use_dev, n, [=] MFEM_HOST_DEVICE (bigint i)
+   {
+      d_data[d_dofs[i]] = d_dofs_vals[i];
+   });
 }
 
 void Vector::Print(std::ostream &os, int width) const
@@ -903,7 +943,6 @@ void Vector::Print(adios2stream &os,
 
 void Vector::Print_HYPRE(std::ostream &os) const
 {
-   int i;
    std::ios::fmtflags old_fmt = os.flags();
    os.setf(std::ios::scientific);
    std::streamsize old_prec = os.precision(14);
@@ -911,7 +950,7 @@ void Vector::Print_HYPRE(std::ostream &os) const
    os << size << '\n';  // number of rows
 
    data.Read(MemoryClass::HOST, size);
-   for (i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       os << ZeroSubnormal(data[i]) << '\n';
    }
@@ -931,7 +970,7 @@ void Vector::PrintMathematica(std::ostream & os) const
    os << "{\n";
 
    data.Read(MemoryClass::HOST, size);
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       os << "Internal`StringToMReal[\"" << ZeroSubnormal(data[i]) << "\"]";
       if (i < size - 1) { os << ','; }
@@ -959,7 +998,7 @@ void Vector::Randomize(int seed)
    srand((unsigned)seed);
 
    HostWrite();
-   for (int i = 0; i < size; i++)
+   for (bigint i = 0; i < size; i++)
    {
       data[i] = rand_real();
    }
@@ -978,6 +1017,7 @@ real_t Vector::Norml2() const
    res.first = 0;
    res.second = 0;
    // first compute sum (|m_data|/scale)^2
+   // FIXME: bigint support
    reduce(size, res, [=] MFEM_HOST_DEVICE(int i, value_type &r)
    {
       real_t n = fabs(m_data[i]);
@@ -1007,6 +1047,7 @@ real_t Vector::Normlinf() const
 
    real_t res = 0;
    const auto m_data = Read(UseDevice());
+   // FIXME: bigint support
    reduce(size, res, [=] MFEM_HOST_DEVICE(int i, real_t &r)
    {
       r = fmax(r, fabs(m_data[i]));
@@ -1021,6 +1062,7 @@ real_t Vector::Norml1() const
 
    real_t res = 0;
    const auto m_data = Read(UseDevice());
+   // FIXME: bigint support
    reduce(size, res, [=] MFEM_HOST_DEVICE(int i, real_t &r)
    {
       r += fabs(m_data[i]);
@@ -1050,6 +1092,7 @@ real_t Vector::Normlp(real_t p) const
       res.second = 0;
       const auto m_data = Read(UseDevice());
       // first compute sum (|m_data|/scale)^p
+      // FIXME: bigint support
       reduce(size, res, [=] MFEM_HOST_DEVICE(int i, value_type &r)
       {
          real_t n = fabs(m_data[i]);
@@ -1097,6 +1140,7 @@ real_t Vector::operator*(const Vector &v) const
    const auto compute_dot = [&]()
    {
       real_t res = 0;
+      // FIXME: bigint support
       reduce(size, res, [=] MFEM_HOST_DEVICE (int i, real_t &r)
       {
          r += m_data[i] * v_data[i];
@@ -1122,11 +1166,11 @@ real_t Vector::operator*(const Vector &v) const
          #pragma omp master
          th_dot.SetSize(nt);
          const int tid = omp_get_thread_num();
-         const int stride = (size + nt - 1) / nt;
-         const int start = tid * stride;
-         const int stop = std::min(start + stride, size);
+         const bigint stride = (size + nt - 1) / nt;
+         const bigint start = tid * stride;
+         const bigint stop = std::min(start + stride, size);
          real_t my_dot = 0.0;
-         for (int i = start; i < stop; i++)
+         for (bigint i = start; i < stop; i++)
          {
             my_dot += m_data[i] * v_data[i];
          }
@@ -1138,7 +1182,7 @@ real_t Vector::operator*(const Vector &v) const
       // The standard way of computing the dot product is non-deterministic
       real_t prod = 0.0;
       #pragma omp parallel for reduction(+ : prod)
-      for (int i = 0; i < size; i++)
+      for (bigint i = 0; i < size; i++)
       {
          prod += m_data[i] * v_data[i];
       }
@@ -1168,6 +1212,7 @@ real_t Vector::Min() const
    const auto compute_min = [&]()
    {
       real_t res = infinity();
+      // FIXME: bigint support
       reduce(size, res, [=] MFEM_HOST_DEVICE(int i, real_t &r)
       {
          r = fmin(r, m_data[i]);
@@ -1185,7 +1230,7 @@ real_t Vector::Min() const
    {
       real_t minimum = m_data[0];
       #pragma omp parallel for reduction(min:minimum)
-      for (int i = 0; i < size; i++)
+      for (bigint i = 0; i < size; i++)
       {
          minimum = std::min(minimum, m_data[i]);
       }
@@ -1214,6 +1259,7 @@ real_t Vector::Max() const
    const auto compute_max = [&]()
    {
       real_t res = -infinity();
+      // FIXME: bigint support
       reduce(size, res, [=] MFEM_HOST_DEVICE(int i, real_t &r)
       {
          r = fmax(r, m_data[i]);
@@ -1231,7 +1277,7 @@ real_t Vector::Max() const
    {
       real_t maximum = m_data[0];
       #pragma omp parallel for reduction(max : maximum)
-      for (int i = 0; i < size; i++)
+      for (bigint i = 0; i < size; i++)
       {
          maximum = fmax(maximum, m_data[i]);
       }
@@ -1249,6 +1295,7 @@ real_t Vector::Sum() const
 
    real_t res = 0;
    const auto m_data = Read(UseDevice());
+   // FIXME: bigint support
    reduce(size, res, [=] MFEM_HOST_DEVICE(int i, real_t &r)
    {
       r += m_data[i];
@@ -1266,10 +1313,13 @@ void Vector::DeleteAt(const Array<int> &indices)
       // extra entry for number of selected out
       Array<int> workspace(size + 1);
       const auto d_flag = workspace.Write(use_dev);
-      mfem::forall_switch(use_dev, size,
-      [=] MFEM_HOST_DEVICE(int i) { d_flag[i] = true; });
+      mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE(bigint i)
+      {
+         d_flag[i] = true;
+      });
       const auto d_indices = indices.Read(use_dev);
-      mfem::forall_switch(use_dev, indices.Size(), [=] MFEM_HOST_DEVICE(int i)
+      mfem::forall_switch(
+         use_dev, indices.Size(), [=] MFEM_HOST_DEVICE(bigint i)
       {
          // fine as long as indices are unique; to support non-unique indices
          // assignment to d_flag must be atomic

linalg/vector.hpp

@@ -40,7 +40,7 @@ namespace mfem
 
 /** Count the number of entries in an array of doubles for which isfinite
     is false, i.e. the entry is a NaN or +/-Inf. */
-inline int CheckFinite(const real_t *v, const int n);
+inline bigint CheckFinite(const real_t *v, const bigint n);
 
 /// Define a shortcut for std::numeric_limits<double>::infinity()
 #ifndef __CYGWIN__
@@ -83,7 +83,7 @@ class Vector
 protected:
 
    Memory<real_t> data;
-   int size;
+   bigint size;
 
 public:
 
@@ -99,30 +99,30 @@ public:
 
    /// @brief Creates vector of size s.
    /// @warning Entries are not initialized to zero!
-   explicit Vector(int s);
+   explicit Vector(bigint s);
 
    /// Creates a vector referencing an array of doubles, owned by someone else.
    /** The pointer @a data_ can be NULL. The data array can be replaced later
        with SetData(). */
-   Vector(real_t *data_, int size_)
+   Vector(real_t *data_, bigint size_)
    { data.Wrap(data_, size_, false); size = size_; }
 
    /** @brief Create a Vector referencing a sub-vector of the Vector @a base
        starting at the given offset, @a base_offset, and size @a size_. */
-   Vector(Vector &base, int base_offset, int size_)
+   Vector(Vector &base, bigint base_offset, bigint size_)
       : data(base.data, base_offset, size_), size(size_) { }
 
    /// Create a Vector of size @a size_ using MemoryType @a mt.
-   Vector(int size_, MemoryType mt)
+   Vector(bigint size_, MemoryType mt)
       : data(size_, mt), size(size_) { }
 
    /** @brief Create a Vector of size @a size_ using host MemoryType @a h_mt and
        device MemoryType @a d_mt. */
-   Vector(int size_, MemoryType h_mt, MemoryType d_mt)
+   Vector(bigint size_, MemoryType h_mt, MemoryType d_mt)
       : data(size_, h_mt, d_mt), size(size_) { }
 
    /// Create a vector from a statically sized C-style array of convertible type
-   template <typename CT, int N>
+   template <typename CT, bigint N>
    explicit Vector(const CT (&values)[N]) : Vector(N)
    { std::copy(values, values + N, begin()); }
 
@@ -130,7 +130,7 @@ public:
    template <typename CT, typename std::enable_if<
                 std::is_convertible<CT,real_t>::value,bool>::type = true>
    explicit Vector(std::initializer_list<CT> values) :
-      Vector(static_cast<int> (values.size()))
+      Vector(static_cast<bigint> (values.size()))
    { std::copy(values.begin(), values.end(), begin()); }
 
    /// Enable execution of Vector operations using the mfem::Device.
@@ -151,10 +151,10 @@ public:
    void Load(std::istream ** in, int np, int * dim);
 
    /// Load a vector from an input stream.
-   void Load(std::istream &in, int Size);
+   void Load(std::istream &in, bigint Size);
 
    /// Load a vector from an input stream, reading the size from the stream.
-   void Load(std::istream &in) { int s; in >> s; Load(in, s); }
+   void Load(std::istream &in) { bigint s; in >> s; Load(in, s); }
 
    /// @brief Resize the vector to size @a s.
    /** If the new size is less than or equal to Capacity() then the internal
@@ -164,16 +164,18 @@ public:
        @warning In the second case above (new size greater than current one),
        the vector will allocate new data array, even if it did not own the
        original data! Also, new entries are not initialized! */
-   void SetSize(int s);
+   void SetSize(bigint s);
 
    /// Resize the vector to size @a s using MemoryType @a mt.
-   void SetSize(int s, MemoryType mt);
+   void SetSize(bigint s, MemoryType mt);
 
    /// Resize the vector to size @a s using the MemoryType of @a v.
-   void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
+   void SetSize(bigint s, const Vector &v)
+   { SetSize(s, v.GetMemory().GetMemoryType()); }
 
-   /// Update \ref Capacity() to @a res (if less than current), keeping existing entries.
-   void Reserve(int res);
+   /** @brief Update \ref Capacity() to @a res (if less than current), keeping
+       existing entries. */
+   void Reserve(bigint res);
 
    /// Delete entries at @a indices and resize vector accordingly.
    /// @warning Indices must be unique!
@@ -188,13 +190,14 @@ public:
        also used as the new Capacity().
        @warning This method should be called only when OwnsData() is false.
        @sa NewDataAndSize(). */
-   void SetDataAndSize(real_t *d, int s) { data.Wrap(d, s, false); size = s; }
+   void SetDataAndSize(real_t *d, bigint s)
+   { data.Wrap(d, s, false); size = s; }
 
    /// Set the Vector data and size, deleting the old data, if owned.
    /** The Vector does not assume ownership of the new data. The new size is
        also used as the new Capacity().
        @sa SetDataAndSize(). */
-   void NewDataAndSize(real_t *d, int s)
+   void NewDataAndSize(real_t *d, bigint s)
    {
       data.Delete();
       SetDataAndSize(d, s);
@@ -209,14 +212,15 @@ public:
        the Vector object takes ownership of all pointers owned by @a mem.
 
        @sa NewDataAndSize(). */
-   inline void NewMemoryAndSize(const Memory<real_t> &mem, int s, bool own_mem);
+   inline void NewMemoryAndSize(const Memory<real_t> &mem, bigint s,
+                                bool own_mem);
 
    /// Reset the Vector to be a reference to a sub-vector of @a base.
-   inline void MakeRef(Vector &base, int offset, int size);
+   inline void MakeRef(Vector &base, bigint offset, bigint size);
 
    /** @brief Reset the Vector to be a reference to a sub-vector of @a base
        without changing its current size. */
-   inline void MakeRef(Vector &base, int offset);
+   inline void MakeRef(Vector &base, bigint offset);
 
    /// Set the Vector data (host pointer) ownership flag.
    void MakeDataOwner() const { data.SetHostPtrOwner(true); }
@@ -231,11 +235,11 @@ public:
    { data.DeleteDevice(copy_to_host); }
 
    /// Returns the size of the vector.
-   inline int Size() const { return size; }
+   inline bigint Size() const { return size; }
 
    /// Return the size of the currently allocated data array.
    /** It is always true that Capacity() >= Size(). */
-   inline int Capacity() const { return data.Capacity(); }
+   inline bigint Capacity() const { return data.Capacity(); }
 
    /// Return a pointer to the beginning of the Vector data.
    /** @warning This method should be used with caution as it gives write access
@@ -286,26 +290,26 @@ public:
    inline real_t *StealData() { real_t *p; StealData(&p); return p; }
 
    /// Access Vector entries. Index i = 0 .. size-1.
-   real_t &Elem(int i);
+   real_t &Elem(bigint i);
 
    /// Read only access to Vector entries. Index i = 0 .. size-1.
-   const real_t &Elem(int i) const;
+   const real_t &Elem(bigint i) const;
 
    /// Access Vector entries using () for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline real_t &operator()(int i);
+   inline real_t &operator()(bigint i);
 
    /// Read only access to Vector entries using () for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline const real_t &operator()(int i) const;
+   inline const real_t &operator()(bigint i) const;
 
    /// Access Vector entries using [] for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline real_t &operator[](int i) { return (*this)(i); }
+   inline real_t &operator[](bigint i) { return (*this)(i); }
 
    /// Read only access to Vector entries using [] for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline const real_t &operator[](int i) const { return (*this)(i); }
+   inline const real_t &operator[](bigint i) const { return (*this)(i); }
 
    /// Dot product with a `double *` array.
    /// This function always executes on the CPU. A HostRead() will be called if
@@ -357,10 +361,10 @@ public:
    Vector &Set(const real_t a, const Vector &x);
 
    /// (*this)[i + offset] = v[i]
-   void SetVector(const Vector &v, int offset);
+   void SetVector(const Vector &v, bigint offset);
 
    /// (*this)[i + offset] += v[i]
-   void AddSubVector(const Vector &v, int offset);
+   void AddSubVector(const Vector &v, bigint offset);
 
    /// (*this) = -(*this)
    void Neg();
@@ -511,7 +515,7 @@ public:
 
    /** @brief Count the number of entries in the Vector for which isfinite
        is false, i.e. the entry is a NaN or +/-Inf. */
-   int CheckFinite() const { return mfem::CheckFinite(HostRead(), size); }
+   bigint CheckFinite() const { return mfem::CheckFinite(HostRead(), size); }
 
    /// Destroys vector.
    virtual ~Vector();
@@ -561,17 +565,17 @@ inline bool IsFinite(const real_t &val)
 #endif
 }
 
-inline int CheckFinite(const real_t *v, const int n)
+inline bigint CheckFinite(const real_t *v, const bigint n)
 {
-   int bad = 0;
-   for (int i = 0; i < n; i++)
+   bigint bad = 0;
+   for (bigint i = 0; i < n; i++)
    {
       if (!IsFinite(v[i])) { bad++; }
    }
    return bad;
 }
 
-inline Vector::Vector(int s)
+inline Vector::Vector(bigint s)
 {
    MFEM_ASSERT(s>=0,"Unexpected negative size.");
    size = s;
@@ -581,7 +585,7 @@ inline Vector::Vector(int s)
    }
 }
 
-inline void Vector::SetSize(int s)
+inline void Vector::SetSize(bigint s)
 {
    if (s == size)
    {
@@ -601,7 +605,7 @@ inline void Vector::SetSize(int s)
    data.UseDevice(use_dev);
 }
 
-inline void Vector::SetSize(int s, MemoryType mt)
+inline void Vector::SetSize(bigint s, MemoryType mt)
 {
    if (mt == data.GetMemoryType())
    {
@@ -630,7 +634,7 @@ inline void Vector::SetSize(int s, MemoryType mt)
    data.UseDevice(use_dev);
 }
 
-inline void Vector::Reserve(int res)
+inline void Vector::Reserve(bigint res)
 {
    if (res > Capacity())
    {
@@ -642,7 +646,7 @@ inline void Vector::Reserve(int res)
    }
 }
 
-inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
+inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, bigint s,
                                      bool own_mem)
 {
    data.Delete();
@@ -657,14 +661,14 @@ inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
    }
 }
 
-inline void Vector::MakeRef(Vector &base, int offset, int s)
+inline void Vector::MakeRef(Vector &base, bigint offset, bigint s)
 {
    data.Delete();
    size = s;
    data.MakeAlias(base.GetMemory(), offset, s);
 }
 
-inline void Vector::MakeRef(Vector &base, int offset)
+inline void Vector::MakeRef(Vector &base, bigint offset)
 {
    data.Delete();
    data.MakeAlias(base.GetMemory(), offset, size);
@@ -678,7 +682,7 @@ inline void Vector::Destroy()
    data.UseDevice(use_dev);
 }
 
-inline real_t &Vector::operator()(int i)
+inline real_t &Vector::operator()(bigint i)
 {
    MFEM_ASSERT(data && i >= 0 && i < size,
                "index [" << i << "] is out of range [0," << size << ")");
@@ -686,7 +690,7 @@ inline real_t &Vector::operator()(int i)
    return data[i];
 }
 
-inline const real_t &Vector::operator()(int i) const
+inline const real_t &Vector::operator()(bigint i) const
 {
    MFEM_ASSERT(data && i >= 0 && i < size,
                "index [" << i << "] is out of range [0," << size << ")");
@@ -712,11 +716,11 @@ inline Vector::~Vector()
    data.Delete();
 }
 
-inline real_t DistanceSquared(const real_t *x, const real_t *y, const int n)
+inline real_t DistanceSquared(const real_t *x, const real_t *y, const bigint n)
 {
    real_t d = 0.0;
 
-   for (int i = 0; i < n; i++)
+   for (bigint i = 0; i < n; i++)
    {
       d += (x[i]-y[i])*(x[i]-y[i]);
    }
@@ -724,7 +728,7 @@ inline real_t DistanceSquared(const real_t *x, const real_t *y, const int n)
    return d;
 }
 
-inline real_t Distance(const real_t *x, const real_t *y, const int n)
+inline real_t Distance(const real_t *x, const real_t *y, const bigint n)
 {
    return std::sqrt(DistanceSquared(x, y, n));
 }

mesh/mesh.cpp

@@ -493,8 +493,7 @@ void Mesh::GetBdrElementTransformation(int i,
          {
             for (int j = 0; j < n; j++)
             {
-               int idx = vdofs[n*k+j];
-               pm(k,j) = nodes((idx<0)? -1-idx:idx);
+               pm(k,j) = nodes(UnsignIndex(vdofs[n*k+j]));
             }
          }
          ElTr->SetFE(bdr_el);
@@ -1356,7 +1355,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
                face.element[0].conformity = ElementConformity::Coincident;
                face.element[1].conformity = ElementConformity::Coincident;
                face.element[1].location = ElementLocation::FaceNbr;
-               face.element[1].index = -1 - e2;
+               face.element[1].index = FlipIndexSign(e2);
                face.element[1].orientation = inf2%64;
             }
          }
@@ -1379,7 +1378,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
                face.element[1].location = ElementLocation::FaceNbr;
                face.element[0].conformity = ElementConformity::Coincident;
                face.element[1].conformity = ElementConformity::Superset;
-               face.element[1].index = -1 - e2;
+               face.element[1].index = FlipIndexSign(e2);
                face.element[1].orientation = inf2%64;
             }
             face.point_matrix = nc_faces_info[ncface].PointMatrix;
@@ -1405,7 +1404,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
          face.element[1].location = ElementLocation::FaceNbr;
          face.element[0].conformity = ElementConformity::Superset;
          face.element[1].conformity = ElementConformity::Coincident;
-         face.element[1].index = -1 - e2;
+         face.element[1].index = FlipIndexSign(e2);
          face.element[1].orientation = inf2%64;
          face.point_matrix = nc_faces_info[ncface].PointMatrix;
       }
@@ -1438,7 +1437,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
          break;
       case FaceInfoTag::SharedConforming:
          res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
          res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
          res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
          break;
@@ -1448,7 +1447,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
          break;
       case FaceInfoTag::SharedSlaveNonconforming:
          res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
          res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
          res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
          break;
@@ -1456,7 +1455,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
          break;
       case FaceInfoTag::GhostSlave:
          res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
          res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
          res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
          break;
@@ -6514,7 +6513,7 @@ void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_ukv)
          input >> edge_to_ukv[j] >> v[0] >> v[1];
          if (v[0] > v[1])
          {
-            edge_to_ukv[j] = -1 - edge_to_ukv[j];
+            edge_to_ukv[j] = FlipIndexSign(edge_to_ukv[j]);
          }
       }
    }
@@ -6551,9 +6550,6 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
    const int NP = NumOfElements; // number of patches
    const int NPKV = NP * dim;    // number of patch knotvectors
    constexpr int notset = -9999999;
-   // Sign convention
-   auto flipSign = [](int i) { return -1 - i; };
-   auto unSign = [](int i) { return (i < 0) ? -1 - i : i; };
    // Local edge index -> dimension convention
    auto edge_to_dim = [](int i) { return (i < 8) ? ((i & 1) ? 1 : 0) : 2; };
 
@@ -6569,7 +6565,7 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
       {
          GetElementVertices(i, v);
          // Sign is based on the edge's vertex indices
-         edge_to_ukv[i] = (v[1] > v[0]) ? i : flipSign(i);
+         edge_to_ukv[i] = (v[1] > v[0]) ? i : FlipIndexSign(i);
          ukv_to_rpkv[i] = i;
       }
       return;
@@ -6619,14 +6615,14 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
          // We've set this edge already - link this index to it
          if (edge_to_pkv[edge] != notset)
          {
-            const int pkv_other = unSign(edge_to_pkv[edge]);
+            const int pkv_other = UnsignIndex(edge_to_pkv[edge]);
             unite(pkv, pkv_other);
          }
          else
          {
             GetEdgeVertices(edge, v);
             // Sign is based on the edge's vertex indices
-            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : flipSign(pkv);
+            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : FlipIndexSign(pkv);
          }
       }
    }
@@ -6653,10 +6649,10 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
    edge_to_ukv.SetSize(NumOfEdges);
    for (int i = 0; i < NumOfEdges; i++)
    {
-      const int pkv = unSign(edge_to_pkv[i]);
+      const int pkv = UnsignIndex(edge_to_pkv[i]);
       const int rpkv = pkv_to_rpkv[pkv];
       const int ukv = rpkv_to_ukv[rpkv];
-      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? flipSign(ukv) : ukv;
+      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? FlipIndexSign(ukv) : ukv;
    }
 
    CorrectPatchTopoOrientations(edge_to_ukv);
@@ -6667,9 +6663,6 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
    const int dim = Dimension(); // Topological (not physical) dimension
    if (dim == 1) { return; }
 
-   // Sign convention
-   auto flipSign = [](int i) { return -1 - i; };
-
    const Table *face2elem = GetFaceToElementTable();
    Array<int> pfaces, orient;
    Array<int> fe, feo;
@@ -6688,7 +6681,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          for (auto e : fe)
          {
             const int skv = edge_to_ukv[e];
-            if (skv == kv || flipSign(skv) == kv) { hasKV = true; }
+            if (skv == kv || FlipIndexSign(skv) == kv) { hasKV = true; }
          }
          if (hasKV)
          {
@@ -6718,7 +6711,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       };
    }
 
-   Array<int> ukvs((dim==2) ? 4 : 12);
+   Array<int> ukvs((dim == 2) ? 4 : 12);
    Array<int> pe, oe;
    bool initKV = false;
 
@@ -6732,7 +6725,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       for (int i = 0; i < pe.Size(); i++)
       {
          ukvs[i] = edge_to_ukv[pe[i]];
-         ukvs[i] = (oe[i] < 0) ? flipSign(ukvs[i]) : ukvs[i];
+         ukvs[i] = (oe[i] < 0) ? FlipIndexSign(ukvs[i]) : ukvs[i];
       }
 
       // Find the direction with this kv.
@@ -6740,12 +6733,19 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       for (int d=0; d<dim; ++d) // Loop over directions.
       {
          const int skv = edge_to_ukv[pe[dir_edges[d][0]]];
-         if (skv == kv || flipSign(skv) == kv)
+         if (skv == kv || FlipIndexSign(skv) == kv)
          {
-            thisDir = d;
+            for (auto e : dir_edges[d])
+               if (!edgeSet[pe[e]])
+               {
+                  thisDir = d;
+               }
          }
       }
-      MFEM_VERIFY(thisDir >= 0, "");
+      if (thisDir == -1)
+      {
+         return false;
+      }
 
       // For this direction, find any edge already set. If no edge is set, we
       // arbitrarily take the first.
@@ -6777,12 +6777,12 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          }
 
          const int edge = pe[i];
-         if ((dim == 2 && ukvs[i] != flipSign(ukvs[ref_edge0])) ||
-             (dim == 3 && ukvs[i] == flipSign(ukvs[ref_edge0])))
+         if ((dim == 2 && ukvs[i] != FlipIndexSign(ukvs[ref_edge0])) ||
+             (dim == 3 && ukvs[i] == FlipIndexSign(ukvs[ref_edge0])))
          {
             // Flip the sign of this edge
-            MFEM_VERIFY(!edgeSet[edge], "");
-            edge_to_ukv[edge] = flipSign(edge_to_ukv[edge]);
+            MFEM_ASSERT(!edgeSet[edge], "");
+            edge_to_ukv[edge] = FlipIndexSign(edge_to_ukv[edge]);
          }
 
          edgeSet[edge] = true;
@@ -6827,10 +6827,11 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       int unsetDim = -1;
       for (int d=0; d<dim; ++d) // Loop over dimensions.
       {
-         if (!edgeSet[pe[dir_edges[d][0]]])
-         {
-            unsetDim = d;
-         }
+         for (auto e : dir_edges[d])
+            if (!edgeSet[pe[e]])
+            {
+               unsetDim = d;
+            }
       }
 
       if (unsetDim == -1)
@@ -6839,9 +6840,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          continue;
       }
 
-      const int kv_signed = edge_to_ukv[pe[dir_edges[unsetDim][0]]];
-      const int kv = kv_signed < 0 ? flipSign(kv_signed) : kv_signed;
-      MFEM_VERIFY(!edgeSet[pe[dir_edges[unsetDim][0]]], "");
+      const int kv = UnsignIndex(edge_to_ukv[pe[dir_edges[unsetDim][0]]]);
 
       initKV = false;
 
@@ -6891,6 +6890,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       }
    }
 
+#ifdef MFEM_DEBUG
    bool allSet = true;
    for (auto eset : edgeSet)
    {
@@ -6899,7 +6899,8 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          allSet = false;
       }
    }
-   MFEM_VERIFY(allSet && unset.size() == 0, "Some edge is not set");
+   MFEM_ASSERT(allSet && unset.size() == 0, "Some edge is not set");
+#endif
 
    delete face2elem;
 }
@@ -6941,7 +6942,7 @@ void Mesh::LoadNonconformingPatchTopo(std::istream &input,
 
       if (v[0] > v[1])
       {
-         ukv = -1 - ukv;
+         ukv = FlipIndexSign(ukv);
       }
       edge_to_ukv[j] = ukv;
    }
@@ -12438,11 +12439,7 @@ void Mesh::PrintTopoEdges(std::ostream &os, const Array<int> &e_to_k,
    for (int i = 0; i < NumOfEdges; i++)
    {
       edge_vertex->GetRow(i, vert);
-      int ki = e_to_k[i];
-      if (ki < 0)
-      {
-         ki = -1 - ki;
-      }
+      const int ki = UnsignIndex(e_to_k[i]);
 
       if (vmap)
       {
@@ -15237,17 +15234,17 @@ void GeometricFactors::Compute(const GridFunction &nodes,
                         Device::GetDeviceMemoryType();
    if (computed_factors & GeometricFactors::COORDINATES)
    {
-      X.SetSize(vdim*NQ*NE, my_d_mt); // NQ x SDIM x NE
+      X.SetSize(bigint(vdim)*NQ*NE, my_d_mt); // NQ x SDIM x NE
       eval_flags |= QuadratureInterpolator::VALUES;
    }
    if (computed_factors & GeometricFactors::JACOBIANS)
    {
-      J.SetSize(dim*vdim*NQ*NE, my_d_mt); // NQ x SDIM x DIM x NE
+      J.SetSize(bigint(dim)*vdim*NQ*NE, my_d_mt); // NQ x SDIM x DIM x NE
       eval_flags |= QuadratureInterpolator::DERIVATIVES;
    }
    if (computed_factors & GeometricFactors::DETERMINANTS)
    {
-      detJ.SetSize(NQ*NE, my_d_mt); // NQ x NE
+      detJ.SetSize(bigint(NQ)*NE, my_d_mt); // NQ x NE
       eval_flags |= QuadratureInterpolator::DETERMINANTS;
    }
 
@@ -15265,7 +15262,7 @@ void GeometricFactors::Compute(const GridFunction &nodes,
 
    if (elem_restr) // Always true as of 2021-04-27
    {
-      Vector Enodes(vdim*ND*NE, my_d_mt);
+      Vector Enodes(bigint(vdim)*ND*NE, my_d_mt);
       elem_restr->Mult(nodes, Enodes);
       qi->Mult(Enodes, eval_flags, X, J, detJ);
    }
@@ -15307,22 +15304,22 @@ FaceGeometricFactors::FaceGeometricFactors(const Mesh *mesh,
 
    if (flags & FaceGeometricFactors::COORDINATES)
    {
-      X.SetSize(vdim*NQ*NF, my_d_mt);
+      X.SetSize(bigint(vdim)*NQ*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::VALUES;
    }
    if (flags & FaceGeometricFactors::JACOBIANS)
    {
-      J.SetSize(vdim*(mesh->Dimension() - 1)*NQ*NF, my_d_mt);
+      J.SetSize(bigint(vdim)*(mesh->Dimension() - 1)*NQ*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::DERIVATIVES;
    }
    if (flags & FaceGeometricFactors::DETERMINANTS)
    {
-      detJ.SetSize(NQ*NF, my_d_mt);
+      detJ.SetSize(bigint(NQ)*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::DETERMINANTS;
    }
    if (flags & FaceGeometricFactors::NORMALS)
    {
-      normal.SetSize(vdim*NQ*NF, my_d_mt);
+      normal.SetSize(bigint(vdim)*NQ*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::NORMALS;
    }
 
@@ -15751,9 +15748,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 +15783,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;
    }
 

mesh/mesh_operators.cpp

@@ -143,7 +143,7 @@ int ThresholdRefiner::ApplyImpl(Mesh &mesh)
    if (num_marked_elements == 0LL) { return STOP; }
 
    mesh.GeneralRefinement(marked_elements, non_conforming, nc_limit);
-   return CONTINUE + REFINED;
+   return static_cast<int>(CONTINUE) + static_cast<int>(REFINED);
 }
 
 void ThresholdRefiner::Reset()
@@ -162,7 +162,7 @@ int ThresholdDerefiner::ApplyImpl(Mesh &mesh)
    const Vector &local_err = estimator.GetLocalErrors();
    bool derefs = mesh.DerefineByError(local_err, threshold, nc_limit, op);
 
-   return derefs ? CONTINUE + DEREFINED : NONE;
+   return derefs ? static_cast<int>(CONTINUE) + static_cast<int>(DEREFINED) : NONE;
 }
 
 
@@ -290,7 +290,7 @@ int CoefficientRefiner::PreprocessMesh(Mesh &mesh, int max_it)
    }
    delete l2fes;
    delete gf;
-   return CONTINUE + REFINED;
+   return static_cast<int>(CONTINUE) + static_cast<int>(REFINED);
 
 }
 
@@ -310,7 +310,7 @@ int Rebalancer::ApplyImpl(Mesh &mesh)
    if (pmesh && pmesh->Nonconforming())
    {
       pmesh->Rebalance();
-      return CONTINUE + REBALANCED;
+      return static_cast<int>(CONTINUE) + static_cast<int>(REBALANCED);
    }
 #endif
    return NONE;

mesh/ncmesh.cpp

@@ -3542,7 +3542,7 @@ void NCMesh::TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
 
             // create a slave face record with a degenerate point matrix
             face_list.slaves.Append(
-               Slave(-1 - enode.edge_index,
+               Slave(FlipIndexSign(enode.edge_index),
                      eid[0].element, eid[0].local, Geometry::SQUARE));
             Slave &sl = face_list.slaves.Last();
 
@@ -3589,7 +3589,7 @@ void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
          // non-slave edge is really a (face-)slave itself.
          const MeshId &eid = *eid_and_type.id;
          face_list.slaves.Append(
-            Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
+            Slave(FlipIndexSign(eid.index), eid.element, eid.local, Geometry::TRIANGLE));
 
          int v0index = nodes[vn0].vert_index;
          int v1index = nodes[vn1].vert_index;

mesh/ncnurbs.cpp

@@ -93,7 +93,7 @@ void NCNURBSExtension::GetMasterEdgeEntities(
       }
       else
       {
-         const int auxEdge = -1 - edge_i;
+         const int auxEdge = FlipIndexSign(edge_i);
          GetAuxEdgeVertices(auxEdge, sverts);
       }
 
@@ -159,7 +159,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
             {
                if (edge < 0)
                {
-                  sideAuxEdges[s].Append(-1 - edge);
+                  sideAuxEdges[s].Append(FlipIndexSign(edge));
                }
                else
                {
@@ -456,7 +456,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                                        == afverts[j], "");
                         }
 
-                        ori_f2 = -1 - ori_f2;
+                        ori_f2 = FlipIndexSign(ori_f2);
                      }
                      else
                      {
@@ -468,7 +468,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                      }
 
                      facePairs.emplace_back(FacePairInfo{fverts[vMinID], f,
-                     SlaveFaceInfo{-1 - afid, ori_f2,
+                     SlaveFaceInfo{FlipIndexSign(afid), ori_f2,
                         {fki(vMinID,0), fki(vMinID,1)},
                         {
                            fki((vMinID + 2) % 4,0) - fki(vMinID,0),
@@ -509,7 +509,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                      auxFaces.push_back(auxFace);
 
                      facePairs.emplace_back(FacePairInfo{fverts[vMinID], f,
-                     SlaveFaceInfo{-1 - auxFaceId, ori_f,
+                     SlaveFaceInfo{FlipIndexSign(auxFaceId), ori_f,
                         {fki(vMinID,0), fki(vMinID,1)},
                         {
                            fki((vMinID + 2) % 4,0) - fki(vMinID,0),
@@ -622,7 +622,7 @@ void NCNURBSExtension::GetAuxFaceEdges(int auxFace, Array<int> &edges) const
       }
       else // Auxiliary edge
       {
-         edges[i] = -1 - auxv2e.at(edge_v);
+         edges[i] = FlipIndexSign(auxv2e.at(edge_v));
       }
    }
 }
@@ -633,7 +633,7 @@ int OffsetHelper(int i, int j, const Array<int> &a, const Array<int> &b)
 {
    if (i < 0)
    {
-      return b[-1 - i + j];
+      return b[FlipIndexSign(i) + j];
    }
    else if (i + j < a.Size())
    {
@@ -679,7 +679,7 @@ void NCNURBSExtension::GetMasterEdgeDofs(bool dof, int me,
       }
       else // Auxiliary edge
       {
-         GetAuxEdgeVertices(-1 - slaveId, svert);
+         GetAuxEdgeVertices(FlipIndexSign(slaveId), svert);
       }
 
       bool reverse = false;
@@ -872,7 +872,7 @@ void ReorderArray2D(int i0, int j0, const Array2D<int> &a,
 // Set a quadrilateral vertex index permutation for a given orientation.
 void GetVertexOrdering(int ori, std::array<int, 4> &perm)
 {
-   const int oriAbs = ori < 0 ? -1 - ori : ori;
+   const int oriAbs = UnsignIndex(ori);
 
    for (int i=0; i<4; ++i)
    {
@@ -1094,7 +1094,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
       if (slaveId < 0)
       {
          // Auxiliary face
-         const int auxFace = -1 - slaveId;
+         const int auxFace = FlipIndexSign(slaveId);
 
          // Set slave face entity dimensions.
          if (dof)
@@ -1171,7 +1171,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
             }
             else
             {
-               const int auxEdge = -1 - edge;
+               const int auxEdge = FlipIndexSign(edge);
                GetAuxEdgeVertices(auxEdge, evert);
             }
             MFEM_ASSERT(evert[0] == vstart || evert[1] == vstart, "");
@@ -1184,7 +1184,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
          // dimensions of the master face, by using ori.
          int e1 = -1, e2 = -1;
          {
-            const int aori = ori < 0 ? -1 - ori : ori;
+            const int aori = UnsignIndex(ori);
             if (aori % 2 == 0)
             {
                e1 = 0;
@@ -1416,14 +1416,15 @@ void NCNURBSExtension::ProcessVertexToKnot2D(const VertexToKnotSpan &v2k,
          {
             // Create a new auxiliary edge
             auxv2e[childPair] = auxEdges.size();
-            auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ?
-                                                parentEdge : -1 - parentEdge,
+            auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ? parentEdge :
+                                                FlipIndexSign(parentEdge),
             {childPair.first, childPair.second},
             {newParentEdge ? 0 : prevKI, ks}});
          }
       }
 
-      const int childEdge = childPairTopo ? v2e[childPair] : -1 - auxv2e[childPair];
+      const int childEdge = childPairTopo ? v2e[childPair] :
+                            FlipIndexSign(auxv2e[childPair]);
 
       // Check whether this is the final vertex in this parent edge. Note that
       // the logic for comparing (pv[0],pv[1]) to the next parents assumes the
@@ -1460,14 +1461,15 @@ void NCNURBSExtension::ProcessVertexToKnot2D(const VertexToKnotSpan &v2k,
 
                // -1 denotes `ne` at endpoint
                auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ?
-                                                   -1 - parentEdge : parentEdge,
+                                                   FlipIndexSign(parentEdge) :
+                                                   parentEdge,
                {finalChildPair.first, finalChildPair.second},
                {ks, -1}});
             }
          }
 
          const int finalChildEdge = finalChildPairTopo ? v2e[finalChildPair] :
-                                    -1 - auxv2e[finalChildPair];
+                                    FlipIndexSign(auxv2e[finalChildPair]);
          edgePairs.emplace_back(-1, -1, finalChildEdge, parentEdge);
       }
 
@@ -1805,7 +1807,7 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
                   auxFaces.push_back(auxFace);
                   facePairs.emplace_back(
                      FacePairInfo{cv[0], parentFace,
-                                  SlaveFaceInfo{-1 - auxv2f[childPair],
+                                  SlaveFaceInfo{FlipIndexSign(auxv2f[childPair]),
                                                 0, {i0, j0}, {d0, d1}}});
                }
             }
@@ -2111,7 +2113,7 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
                      auxv2e[childPair] = auxEdges.size();
                      auxEdges.emplace_back(AuxiliaryEdge{pv0 < pv1 ?
                                                          parentEdge :
-                                                         -1 - parentEdge,
+                                                         FlipIndexSign(parentEdge),
                      {childPair.first, childPair.second},
                      {knotIndex0, knotIndex1}});
                   }
@@ -2131,7 +2133,8 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
 
                   const EdgePairInfo ep_e((e_idx == n_d - de) ? -1 : tv,
                                           (e_idx == n_d - de) ? -1 : tvki,
-                                          -1 - auxv2e[childPair], parentEdge);
+                                          FlipIndexSign(auxv2e[childPair]),
+                                          parentEdge);
 
                   const bool unset = !edgePairs[edgePairOS[parentEdge] + e_idx].isSet;
                   if (unset)
@@ -2226,7 +2229,7 @@ void NCNURBSExtension::GetAuxFaceToPatchTable(Array2D<int> &auxface2patch)
                if (s < 0)
                {
                   // Auxiliary face.
-                  const int aux = -1 - s;
+                  const int aux = FlipIndexSign(s);
                   if (auxface2patch(aux, 0) >= 0)
                   {
                      if (auxface2patch(aux, 1) != -1) { consistent = false; }
@@ -2316,7 +2319,7 @@ void NCNURBSExtension::UpdateAuxiliaryKnotSpans(const Array<int> &rf)
    for (auto auxEdge : auxEdges)
    {
       const int p = auxEdge.parent;
-      const int parent = p < 0 ? -1 - p : p;
+      const int parent = UnsignIndex(p);
       const int kv = KnotInd(parent);
       for (int i=0; i<2; ++i)
       {
@@ -2382,14 +2385,8 @@ int NCNURBSExtension::AuxiliaryEdgeNE(int aux_edge)
    const int signedParentEdge = auxEdges[aux_edge].parent;
    const int ki0 = auxEdges[aux_edge].ksi[0];
    const int ki1raw = auxEdges[aux_edge].ksi[1];
-   int ki1 = ki1raw;
-   if (ki1raw == -1)
-   {
-      const bool rev = signedParentEdge < 0;
-      const int parentEdge = rev ? -1 - signedParentEdge : signedParentEdge;
-      ki1 = KnotVec(parentEdge)->GetNE();
-   }
-
+   const int ki1 = ki1raw == -1 ? KnotVec(UnsignIndex(signedParentEdge))->GetNE()
+                   : ki1raw;
    return ki1 - ki0;
 }
 
@@ -2403,7 +2400,7 @@ void NCNURBSExtension::SlaveEdgeToParent(int se, int parent,
    Array<int> sev(2);
    if (se < 0) // Auxiliary edge
    {
-      for (int i=0; i<2; ++i) { sev[i] = auxEdges[-1 - se].v[i]; }
+      for (int i=0; i<2; ++i) { sev[i] = auxEdges[FlipIndexSign(se)].v[i]; }
    }
    else
    {
@@ -2459,7 +2456,7 @@ void NCNURBSExtension::GetMasterEdgePieceOffsets(int mid, Array<int> &os)
       }
       else
       {
-         nes = AuxiliaryEdgeNE(-1 - s);
+         nes = AuxiliaryEdgeNE(FlipIndexSign(s));
       }
 
       os[i+1] = os[i] + nes;
@@ -2565,7 +2562,7 @@ int NCNURBSExtension::SetPatchFactors(int p)
                }
                else // Aux edge
                {
-                  const int aux_edge = -1 - s;
+                  const int aux_edge = FlipIndexSign(s);
                   if (auxef[aux_edge].Size() == 0)
                   {
                      auxef[aux_edge].SetSize(AuxiliaryEdgeNE(aux_edge));
@@ -2611,7 +2608,7 @@ int NCNURBSExtension::SetPatchFactors(int p)
    }
 
    MFEM_VERIFY(consistent, "");
-   return partialChange ? -1 - dirSet : dirSet;
+   return partialChange ? FlipIndexSign(dirSet) : dirSet;
 }
 
 void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
@@ -2715,7 +2712,7 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
                if (s < 0)
                {
                   // Auxiliary face.
-                  const int aux = -1 - s;
+                  const int aux = FlipIndexSign(s);
                   for (int i=0; i<2; ++i)
                   {
                      const int patch = auxface2patch(aux, i);
@@ -2763,7 +2760,7 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
 
          const int dirSetSigned = SetPatchFactors(p);
          const bool partialChange = dirSetSigned < 0;
-         const int dirSet = partialChange ? -1 - dirSetSigned : dirSetSigned;
+         const int dirSet = UnsignIndex(dirSetSigned);
          const bool changed = (patchState[p] != dirSet) || partialChange;
          patchState[p] = dirSet;
 
@@ -2806,8 +2803,8 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
             {
                const int dirSetSigned_i = SetPatchFactors(i);
                const bool partialChange_i = dirSetSigned_i < 0;
-               const int dirSet_i = partialChange_i ? -1 - dirSetSigned_i :
-                                    dirSetSigned_i;
+               const int dirSet_i = partialChange_i ?
+                                    FlipIndexSign(dirSetSigned_i) : dirSetSigned_i;
                const bool changed_i = (patchState[i] != dirSet_i) ||
                                       partialChange_i;
                patchState[p] = dirSet_i;
@@ -3027,7 +3024,7 @@ int GetFaceOrientation(const Mesh *mesh, const int face,
 
    // Check whether ordering is reversed.
    const bool rev = verts[(s + 1) % 4] != fverts[1];
-   if (rev) { s = -1 - s; } // Reversed order is encoded by the sign.
+   if (rev) { s = FlipIndexSign(s); } // Reversed order is encoded by the sign.
    return s;
 }
 
@@ -3040,7 +3037,7 @@ int GetFaceOrientation(const Mesh *mesh, const int face,
 // see GetFaceOrientation.
 bool Reorder2D(int ori, std::array<int, 2> &s0)
 {
-   const int shift = ori < 0 ? -1 - ori : ori;
+   const int shift = UnsignIndex(ori);
 
    // Shift is an F1 index in the counter-clockwise ordering of 4 quad vertices.
    // Now find the (i,j) indices of this index, with i,j in {0,1}.
@@ -3064,7 +3061,7 @@ void GetInverseShiftedDimensions2D(int signedShift, int sm, int sn, int &m,
                                    int &n)
 {
    const bool rev = (signedShift < 0);
-   const int shift = rev ? -1 - signedShift : signedShift;
+   const int shift = UnsignIndex(signedShift);
    MFEM_ASSERT(0 <= shift && shift < 4, "");
 
    // We consider 8 cases for the possible values of rev and shift.
@@ -3136,7 +3133,7 @@ void GetShiftedGridPoints2D(int m, int n, int i, int j, int signedShift,
                             int& sm, int& sn, int& si, int& sj)
 {
    const bool rev = (signedShift < 0);
-   const int shift = rev ? -1 - signedShift : signedShift;
+   const int shift = UnsignIndex(signedShift);
    MFEM_ASSERT(0 <= shift && shift < 4, "");
 
    // (0,0) <= (i,j) < (m,n) are old indices, and old vertex [shift] maps
@@ -3798,8 +3795,7 @@ void NCNURBSExtension::GenerateOffsets()
       const int signedParentEdge = auxEdges[e].parent;
       const int ki0 = auxEdges[e].ksi[0];
       const int ki1raw = auxEdges[e].ksi[1];
-      const bool rev = signedParentEdge < 0;
-      const int parentEdge = rev ? -1 - signedParentEdge : signedParentEdge;
+      const int parentEdge = UnsignIndex(signedParentEdge);
       const int masterNE = KnotVec(parentEdge)->GetNE();
       const int ki1 = ki1raw == -1 ? masterNE : ki1raw;
       const int perEdgeCP = GetNCPperEdge(KnotVec(e));

mesh/nurbs.cpp

@@ -3723,10 +3723,7 @@ bool NURBSExtension::CheckPatches()
       for (int i = 0; i < edges.Size(); i++)
       {
          edges[i] = edge_to_ukv[edges[i]];
-         if (oedge[i] < 0)
-         {
-            edges[i] = -1 - edges[i];
-         }
+         if (oedge[i] < 0) { edges[i] = FlipIndexSign(edges[i]); }
       }
 
       // In 2d - opposite edges must be same knotvector with opposite sign.
@@ -3740,7 +3737,7 @@ bool NURBSExtension::CheckPatches()
       // {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7} for Geometry::CUBE in 3D
       // See fem/geom.cpp for these definitions.
       if ((dim == 2 &&
-           (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
+           (edges[0] != FlipIndexSign(edges[2]) || edges[1] != FlipIndexSign(edges[3]))) ||
 
           (dim == 3 &&
            (edges[0] != edges[2] || edges[0] != edges[4] ||
@@ -3769,7 +3766,7 @@ void NURBSExtension::CheckBdrPatches()
          edges[i] = edge_to_ukv[edges[i]];
          if (oedge[i] < 0)
          {
-            edges[i] = -1 - edges[i];
+            edges[i] = FlipIndexSign(edges[i]);
          }
       }
 
@@ -4766,14 +4763,13 @@ void NURBSExtension::GenerateBdrElementDofTable()
    SetPatchToBdrElements();
 
    int *dof = bel_dof->GetJ();
-   int ndof = bel_dof->Size_of_connections();
+   const int ndof = bel_dof->Size_of_connections();
    for (int i = 0; i < ndof; i++)
    {
-      int idx = dof[i];
+      const int idx = dof[i];
       if (idx < 0)
       {
-         dof[i] = -1 - (activeDof[-1-idx] - 1);
-         dof[i] = -activeDof[-1-idx];
+         dof[i] = -activeDof[FlipIndexSign(idx)];
       }
       else
       {
@@ -4858,12 +4854,12 @@ void NURBSExtension::Generate2DBdrElementDofTable()
                   for (int ii = 0; ii <= ord0; ii++)
                   {
                      conn.to = DofMap(p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]);
-                     if (s == -1) { conn.to = -1 -conn.to; }
+                     if (s == -1) { conn.to = FlipIndexSign(conn.to); }
                      bel_dof_list.Append(conn);
                   }
                }
                bel_to_patch[lbe] = b;
-               bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
+               bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : FlipIndexSign(i);
                lbe++;
             }
             gbe++;
@@ -4936,14 +4932,14 @@ void NURBSExtension::Generate3DBdrElementDofTable()
                            {
                               const int ii_ = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
                               conn.to = DofMap(p2g(ii_, jj_));
-                              if (s == -1) { conn.to = -1 -conn.to; }
+                              if (s == -1) { conn.to = FlipIndexSign(conn.to); }
                               bel_dof_list.Append(conn);
                            }
                         }
                      }
                      bel_to_patch[lbe] = b;
-                     bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
-                     bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
+                     bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : FlipIndexSign(i);
+                     bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : FlipIndexSign(j);
                      lbe++;
                   }
                   gbe++;

mesh/nurbs.hpp

@@ -220,7 +220,7 @@ public:
        @a u.
        The main purpose of this function is its use in FindInterpolant.
        Use GetBotella instead for each shape function separately, perhaps in
-       conjuction with GetSpan and GetRefPoint.*/
+       conjunction with GetSpan and GetRefPoint.*/
    MFEM_DEPRECATED void FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const;
 
    /** @brief Global curve interpolation through the points @a x (overwritten).
@@ -1398,8 +1398,7 @@ inline const real_t &NURBSPatch::operator()(int i, int j, int k, int l) const
 
 inline int NURBSExtension::KnotInd(int edge) const
 {
-   const int kv = edge_to_ukv[edge];
-   return kv >= 0 ? kv : -1 - kv;
+   return UnsignIndex(edge_to_ukv[edge]);
 }
 
 inline int NURBSExtension::KnotSign(int edge) const
@@ -1429,7 +1428,7 @@ const
    else
    {
       *okv = -oedge;
-      return knotVectors[-1-kv];
+      return knotVectors[FlipIndexSign(kv)];
    }
 }
 

mesh/pncmesh.cpp

@@ -400,7 +400,7 @@ void ParNCMesh::MakeSharedList(const NCList &list, NCList &shared)
          }
          else // special case: prism edge-face constraint
          {
-            if (entity_owner[1][-1-si] != MyRank)
+            if (entity_owner[1][FlipIndexSign(si)] != MyRank)
             {
                master_flag |= 0x2;
             }
@@ -571,9 +571,10 @@ void ParNCMesh::CalculatePMatrixGroups()
       ranks.SetSize(0);
       for (int j = master_face.slaves_begin; j < master_face.slaves_end; j++)
       {
-         int si = face_list.slaves[j].index;
-         int owner = (si >= 0) ? entity_owner[2][si] // standard face dependency
-                     /*     */ : entity_owner[1][-1 - si]; // prism edge-face dep
+         const int si = face_list.slaves[j].index;
+         const int owner =
+            (si >= 0) ? entity_owner[2][si] : // standard face dependency
+            entity_owner[1][FlipIndexSign(si)]; // prism edge-face dep
          ranks.Append(groups[owner][0]);
       }
       ranks.Sort();
@@ -1181,7 +1182,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
       if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }
 
       Mesh::FaceInfo &fi = pmesh.faces_info[cf.index];
-      fi.Elem2No = -1 - fnbr_index[e[0]->index - NElements];
+      fi.Elem2No = FlipIndexSign(fnbr_index[e[0]->index - NElements]);
 
       if (Dim == 3)
       {
@@ -1270,7 +1271,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
                // In other words, side 1 IS the side that generated the face.
             }
             MFEM_ASSERT(fi.Elem2No >= NElements, "");
-            fi.Elem2No = -1 - fnbr_index[fi.Elem2No - NElements];
+            fi.Elem2No = FlipIndexSign(fnbr_index[fi.Elem2No - NElements]);
 
             const DenseMatrix* pm = full_list.point_matrices[sf.geom][sf.matrix];
             if (!sloc && Dim == 3)
@@ -2286,7 +2287,7 @@ void ParNCMesh::Derefine(const Array<int> &derefs)
       if (element_type[index] == 0)
       {
          // this coarse element will get pruned, encode who owns it now
-         index = -1 - elements[coarse[i]].rank;
+         index = FlipIndexSign(elements[coarse[i]].rank);
       }
       transforms.embeddings[i].parent = index;
    }

mesh/submesh/pncsubmesh.cpp

@@ -78,7 +78,7 @@ ParNCSubMesh::ParNCSubMesh(ParSubMesh& submesh, const ParNCMesh &parent,
 #ifdef MFEM_DEBUG
    // Check all processors have the same number of roots
    {
-      int p[2] = {root_state.Size(), -root_state.Size()};
+      int p[2] = {(int)root_state.Size(), -(int)root_state.Size()};
       MPI_Allreduce(MPI_IN_PLACE, p, 2, MPI_INT, MPI_MIN, submesh.GetComm());
       MFEM_ASSERT(p[0] == -p[1], "Ranks must agree on number of root elements: min "
                   << p[0] << " max " << -p[1] << " local " << root_state.Size() << " MyRank " <<
@@ -154,4 +154,4 @@ ParNCSubMesh::ParNCSubMesh(ParSubMesh& submesh, const ParNCMesh &parent,
 
 } // namespace mfem
 
-#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MPI

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/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;
       }
    }
 

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);
    }
 }

tests/unit/fem/test_derefine.cpp

@@ -213,7 +213,7 @@ void test_derefine_L2_element(int order, Element::Type el_type, int basis_type)
    double eps = 1.e-3;
 
    // limit to max 20 dofs for efficiency in 3D
-   int test_ndofs = std::min(coarse_soln_v.Size(), 20);
+   int test_ndofs = std::min(coarse_soln_v.Size(), bigint(20));
    for (int i = 0; i < test_ndofs; i++)
    {
       for (int f = -1; f <= 1; f += 2)

tests/unit/general/test_array.cpp

@@ -125,6 +125,59 @@ TEST_CASE("Array stl-interactions", "[Array]")
    }
 }
 
+
+TEST_CASE("Array move assignment to view (MakeRef)", "[Array]")
+{
+   constexpr int n = 5;
+
+   // Helper function that returns an Array by value (rvalue)
+   auto make_sequence_array = [](int n)
+   {
+      Array<int> arr(n);
+      for (int i = 0; i < n; i++)
+      {
+         arr[i] = i + 1;  // 1, 2, 3, ...
+      }
+      return arr;
+   };
+
+   // Create backing storage
+   Array<int> backing1(n), backing2(n);
+   backing1 = 0.0; backing2 = 0.0;
+
+   // Create a view into the backing storage
+   Array<int> view1, view2;
+   view1.MakeRef(backing1);
+   view2.MakeRef(backing2);
+
+   auto seq_array = make_sequence_array(n);
+   view1 = seq_array; // copy assign
+   view2 = std::move(seq_array); // move assign
+
+   CHECK(seq_array.Size() == 0); // seq_array is invalidated by the move
+
+   // Both assignments should be semantically equivalent.
+   for (int i = 0; i < n; i++)
+   {
+      CHECK(backing1[i] == backing2[i]);
+   }
+
+   // Create backing storage as raw array
+   int backing3[n] = {0, 0, 0, 0, 0};
+
+   // Create a view into the backing storage
+   Array<int> view3;
+   view3.MakeRef(backing3, n);
+
+   // Assign from rvalue
+   view3 = std::move(view2);
+
+   for (int i = 0; i < n; i++)
+   {
+      CHECK(backing3[i] == backing1[i]);
+   }
+}
+
 TEST_CASE("Array delete at indices", "[Array],[GPU]")
 {
    for (int use_dev = 0; use_dev < 2; use_dev++)

tests/unit/linalg/test_matrix_dense.cpp

@@ -961,3 +961,38 @@ TEST_CASE("NNLS", "[DenseMatrix]")
 }
 
 #endif // if MFEM_USE_LAPACK
+
+TEST_CASE("DenseTensor slice copy and move assign equivalency",
+          "[DenseMatrix][DenseTensor]")
+{
+   auto fill_matrix = [](int n)
+   {
+      DenseMatrix M(n, n);
+      M = 0.0;
+      for (int i = 0; i < n; i++)
+         for (int j = 0; j < n; j++)
+         {
+            M(i, j) = i+j*n+1;
+         }
+      return M;
+   };
+
+   constexpr int n = 3;
+   constexpr int k = 2;
+
+   DenseTensor tensor1(n, n, k), tensor2(n, n, k);
+   tensor1 = 0.0; tensor2 = 0.0;
+
+   DenseMatrix temp = fill_matrix(n);
+   tensor1(0) = temp; // copy assign
+   tensor2(0) = std::move(temp); // move assign
+
+   // Check that the tensor was actually updated
+   for (int i = 0; i < n; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         CHECK(tensor1(i, j, 0) == tensor2(i, j, 0));
+      }
+   }
+}

tests/unit/mesh/test_nurbs.cpp

@@ -222,13 +222,13 @@ TEST_CASE("Greville, Botella and Demko points", "[NURBS]")
       mfem::out<<"Knotvector : "; kvp.Print(mfem::out);
 
       kvp.GetGreville(xi);
-      mfem::out<<"Greville points : "; xi.Print(std::cout,999);
+      mfem::out<<"Greville points : "; xi.Print(mfem::out,999);
 
       kvp.GetBotella(xi);
-      mfem::out<<"Botella  points : "; xi.Print(std::cout,999);
+      mfem::out<<"Botella  points : "; xi.Print(mfem::out,999);
 
       kvp.GetDemko(xi);
-      mfem::out<<"Demko    points : "; xi.Print(std::cout,999);
+      mfem::out<<"Demko    points : "; xi.Print(mfem::out,999);
    }
 
    KnotVector kv(3, Vector({0.0, 0.3, 0.3, 0.3, 0.6, 1.0}));