Merge branch 'master' into dfem-transpose

Add if constexpr to Kernel Specializations

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

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/dfem/doperator.hpp

@@ -84,29 +84,33 @@ public:
       const std::vector<derivative_action_t> &derivative_actions,
       const FieldDescriptor &direction,
       const int &daction_l_size,
-      const std::vector<derivative_action_t> &derivative_actions_transpose,
+      const int &derivative_action_tr_l_size,
+      const std::vector<derivative_action_t> &derivative_tr_actions,
       const FieldDescriptor &transpose_direction,
       const int &daction_transpose_l_size,
       const std::vector<Vector *> &solutions_l,
       const std::vector<Vector *> &parameters_l,
       const restriction_callback_t &restriction_callback,
       const std::function<void(Vector &, Vector &)> &prolongation_transpose,
+      const std::function<void(Vector &, Vector &)> &tr_prolongation_transpose,
       const std::vector<assemble_derivative_sparsematrix_callback_t>
       &assemble_derivative_sparsematrix_callbacks,
-      const std::vector<assemble_derivative_hypreparmatrix_callback_t>
-      &assemble_derivative_hypreparmatrix_callbacks) :
+      const assemble_derivative_hypreparmatrix_callback_t
+      &assemble_derivative_hypreparmatrix_callback) :
       Operator(height, width),
       derivative_actions(derivative_actions),
       direction(direction),
       daction_l(daction_l_size),
       daction_l_size(daction_l_size),
-      derivative_actions_transpose(derivative_actions_transpose),
+      derivative_action_tr_l_size(derivative_action_tr_l_size),
+      derivative_tr_actions(derivative_tr_actions),
       transpose_direction(transpose_direction),
       prolongation_transpose(prolongation_transpose),
+      tr_prolongation_transpose(tr_prolongation_transpose),
       assemble_derivative_sparsematrix_callbacks(
          assemble_derivative_sparsematrix_callbacks),
-      assemble_derivative_hypreparmatrix_callbacks(
-         assemble_derivative_hypreparmatrix_callbacks)
+      assemble_derivative_hypreparmatrix_callback(
+         assemble_derivative_hypreparmatrix_callback)
    {
       std::vector<Vector> s_l(solutions_l.size());
       for (size_t i = 0; i < s_l.size(); i++)
@@ -156,18 +160,18 @@ public:
    /// direction_t on T-dofs.
    void MultTranspose(const Vector &direction_t, Vector &result_t) const override
    {
-      MFEM_ASSERT(!derivative_actions_transpose.empty(),
+      MFEM_ASSERT(!derivative_tr_actions.empty(),
                   "derivative can't be used to be multiplied in transpose mode");
 
-      daction_l.SetSize(width);
+      daction_l.SetSize(derivative_action_tr_l_size);
       daction_l = 0.0;
 
       prolongation(transpose_direction, direction_t, direction_l);
-      for (const auto &f : derivative_actions_transpose)
+      for (const auto &f : derivative_tr_actions)
       {
          f(fields_e, direction_l, daction_l);
       }
-      prolongation_transpose(daction_l, result_t);
+      tr_prolongation_transpose(daction_l, result_t);
    };
 
    /// @brief Assemble the derivative operator into a SparseMatrix.
@@ -183,6 +187,10 @@ public:
       {
          f(fields_e, A);
       }
+
+      // SparseMatrix A is finalized after all callbacks have contributed to
+      // it.
+      A->Finalize();
    }
 
    /// @brief Assemble the derivative operator into a HypreParMatrix.
@@ -191,13 +199,7 @@ public:
    /// be an uninitialized object.
    void Assemble(HypreParMatrix *&A)
    {
-      MFEM_ASSERT(!assemble_derivative_hypreparmatrix_callbacks.empty(),
-                  "derivative can't be assembled into a HypreParMatrix");
-
-      for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
-      {
-         f(fields_e, A);
-      }
+      assemble_derivative_hypreparmatrix_callback(fields_e, A);
    }
 
 private:
@@ -212,10 +214,12 @@ private:
 
    const int daction_l_size;
 
+   const int derivative_action_tr_l_size;
+
    /// Transpose Derivative action callbacks. Depending on the requested
    /// derivatives in DifferentiableOperator the callbacks represent certain
    /// combinations of actions of derivatives of the forward operator.
-   std::vector<derivative_action_t> derivative_actions_transpose;
+   std::vector<derivative_action_t> derivative_tr_actions;
 
    FieldDescriptor transpose_direction;
 
@@ -225,13 +229,15 @@ private:
 
    std::function<void(Vector &, Vector &)> prolongation_transpose;
 
+   std::function<void(Vector &, Vector &)> tr_prolongation_transpose;
+
    /// Callbacks that assemble derivatives into a SparseMatrix.
    std::vector<assemble_derivative_sparsematrix_callback_t>
    assemble_derivative_sparsematrix_callbacks;
 
    /// Callbacks that assemble derivatives into a HypreParMatrix.
-   std::vector<assemble_derivative_hypreparmatrix_callback_t>
-   assemble_derivative_hypreparmatrix_callbacks;
+   assemble_derivative_hypreparmatrix_callback_t
+   assemble_derivative_hypreparmatrix_callback;
 };
 
 /// Class representing a differentiable operator which acts on solution and
@@ -457,7 +463,10 @@ public:
          dir_l = s_l[derivative_idx];
       }
 
-      derivative_setup_callbacks[derivative_id][0](fields_e, dir_l);
+      for (size_t i = 0; i < derivative_setup_callbacks[derivative_id].size(); i++)
+      {
+         derivative_setup_callbacks[derivative_id][i](fields_e, dir_l);
+      }
 
       return std::make_shared<DerivativeOperator>(
                 height,
@@ -465,15 +474,17 @@ public:
                 derivative_action_callbacks[derivative_id],
                 fields[derivative_idx],
                 residual_l.Size(),
-                daction_transpose_callbacks[derivative_id],
+                derivative_action_tr_l_size[derivative_id],
+                derivative_action_tr_callbacks[derivative_id],
                 fields[test_space_field_idx],
                 GetVSize(fields[test_space_field_idx]),
                 sol_l,
                 par_l,
                 restriction_callback,
                 prolongation_transpose,
+                derivative_tr_prolongation_transpose[derivative_id],
                 assemble_derivative_sparsematrix_callbacks[derivative_id],
-                assemble_derivative_hypreparmatrix_callbacks[derivative_id]);
+                assemble_derivative_hypreparmatrix_callback[derivative_id]);
    }
 
 private:
@@ -486,13 +497,15 @@ private:
    std::map<size_t,
        std::vector<derivative_action_t>> derivative_action_callbacks;
    std::map<size_t,
-       std::vector<derivative_action_t>> daction_transpose_callbacks;
+       std::vector<derivative_action_t>> derivative_action_tr_callbacks;
+   std::map<size_t,
+       std::function<void(Vector &, Vector &)>> derivative_tr_prolongation_transpose;
+   std::map<size_t, int> derivative_action_tr_l_size;
    std::map<size_t,
        std::vector<assemble_derivative_sparsematrix_callback_t>>
        assemble_derivative_sparsematrix_callbacks;
-   std::map<size_t,
-       std::vector<assemble_derivative_hypreparmatrix_callback_t>>
-       assemble_derivative_hypreparmatrix_callbacks;
+   std::map<size_t, assemble_derivative_hypreparmatrix_callback_t>
+   assemble_derivative_hypreparmatrix_callback;
 
    std::vector<FieldDescriptor> solutions;
    std::vector<FieldDescriptor> parameters;
@@ -510,7 +523,7 @@ private:
    std::function<void(Vector &, Vector &)> output_restriction_transpose;
    restriction_callback_t restriction_callback;
 
-   std::map<size_t, Vector> derivative_qp_caches;
+   std::map<size_t, std::vector<Vector>> derivative_qp_caches;
 
    std::map<size_t, size_t> assembled_vector_sizes;
 
@@ -769,9 +782,10 @@ void DifferentiableOperator::AddIntegrator(
    auto input_size_on_qp =
       get_input_size_on_qp(inputs, std::make_index_sequence<num_inputs> {});
 
+   // printf("calculate shmem action info\n");
    auto action_shmem_info =
       get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>
-      (input_dtq_maps, output_dtq_maps, fields, num_entities, inputs, num_qp,
+      (input_dtq_maps, output_dtq_maps, fields, num_entities, num_qp,
        input_size_on_qp, residual_size_on_qp, element_dof_ordering);
 
    Vector shmem_cache(action_shmem_info.total_size);
@@ -894,9 +908,10 @@ void DifferentiableOperator::AddIntegrator(
          const int da_size_on_qp =
             GetSizeOnQP<entity_t>(output_fop, fields[test_space_field_idx]);
 
+         // printf("calculate shmem derivative action info\n");
          auto shmem_info =
             get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>(
-               input_dtq_maps, output_dtq_maps, fields, num_entities, inputs,
+               input_dtq_maps, output_dtq_maps, fields, num_entities,
                num_qp, input_size_on_qp, residual_size_on_qp,
                element_dof_ordering, d_field_idx);
 
@@ -968,9 +983,12 @@ void DifferentiableOperator::AddIntegrator(
          // Quadrature point local derivative cache for each element, with data
          // layout:
          // [test_vdim, test_op_dim, trial_vdim, trial_op_dim, qp, num_entities].
-         derivative_qp_caches[derivative_id] = Vector(test_vdim * test_op_dim *
-                                                      trial_vdim *
-                                                      total_trial_op_dim * num_qp * num_entities);
+         derivative_qp_caches[derivative_id].push_back(
+            Vector(test_vdim * test_op_dim * trial_vdim * total_trial_op_dim * num_qp *
+                   num_entities));
+
+         const int cache_index = this->derivative_qp_caches[derivative_id].size() - 1;
+
          // Create local references for MSVC lambda capture compatibility
          auto& fields_ref = this->fields;
          auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_id];
@@ -1012,6 +1030,8 @@ void DifferentiableOperator::AddIntegrator(
                trial_vdim,
                inputs_trial_op_dim,
 
+               qpdc_idx = cache_index,
+
                // capture by ref:
                &qpdc_mem = derivative_qp_caches_ref
             ](std::vector<Vector> &f_e, const Vector &dir_l) mutable
@@ -1024,7 +1044,7 @@ void DifferentiableOperator::AddIntegrator(
                                                shmem_info.direction_size,
                                                num_entities);
 
-            auto qpdc = Reshape(qpdc_mem.ReadWrite(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].ReadWrite(), test_vdim, test_op_dim,
                                 trial_vdim, total_trial_op_dim, num_qp, num_entities);
 
             auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1095,6 +1115,7 @@ void DifferentiableOperator::AddIntegrator(
                inputs_trial_op_dim,
                total_trial_op_dim,
                trial_vdim,
+               qpdc_idx = cache_index,
                // capture by ref:
                &qpdc_mem = derivative_qp_caches_ref,
                &or_transpose
@@ -1112,7 +1133,7 @@ void DifferentiableOperator::AddIntegrator(
                                                shmem_info.direction_size,
                                                num_entities);
 
-            auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
                                 trial_vdim, total_trial_op_dim, num_qp, num_entities);
 
             auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1157,6 +1178,154 @@ void DifferentiableOperator::AddIntegrator(
             or_transpose(derivative_action_e, der_action_l);
          });
 
+         // This prevents Sum/Identity. These are invalid
+         // as input FieldOperators anyways.
+         constexpr auto dummy_fop = Value<0> {};
+
+         auto [input_rt,
+               input_e_sz] = get_restriction_transpose<entity_t>
+                             (fields[d_field_idx],
+                              element_dof_ordering, dummy_fop);
+
+         const auto input_restriction_transpose = input_rt;
+
+         derivative_tr_prolongation_transpose[derivative_id] =
+            get_prolongation_transpose(
+               fields[d_field_idx], dummy_fop, mesh.GetComm());
+
+         const auto d_tr_field_idx = test_space_field_idx;
+         const auto direction_tr = fields[d_tr_field_idx];
+
+         auto output_size_on_qp =
+            get_input_size_on_qp(outputs, std::make_index_sequence<num_outputs> {});
+
+         const int residual_tr_size_on_qp = trial_vdim * total_trial_op_dim;
+
+         auto shmem_tr_info =
+            get_shmem_info<entity_t, num_fields, num_outputs, num_inputs>(
+               output_dtq_maps, input_dtq_maps, fields, num_entities,
+               num_qp, output_size_on_qp, residual_tr_size_on_qp,
+               element_dof_ordering, test_space_field_idx);
+
+         // print_shared_memory_info(shmem_tr_info);
+
+         // TODO: this is a hack to extend the shared memory with a known
+         // offset for a temp variable
+         Vector shmem_tr_cache(shmem_tr_info.total_size + residual_tr_size_on_qp *
+                               num_qp);
+
+         Vector direction_tr_e(get_restriction<entity_t>(
+                                  fields[test_space_field_idx],
+                                  element_dof_ordering)->Height());
+
+         derivative_action_tr_l_size[derivative_id] =
+            get_restriction<entity_t>(fields[d_field_idx],
+                                      element_dof_ordering)->Width();
+
+         Vector derivative_action_tr_e(input_e_sz);
+         derivative_action_tr_e = 0.0;
+
+         derivative_action_tr_callbacks[derivative_id].push_back(
+            [
+               // capture by copy:
+               dimension,             // int
+               num_entities,          // int
+               num_trial_dof,         // int
+               num_qp,                // int
+               q1d,                   // int
+               test_vdim,             // int (= output_fop.vdim)
+               test_op_dim,           // int (derived from output_fop)
+               inputs,                // mfem::future::tuple
+               outputs,               // mfem::future::tuple
+               attributes,            // Array<int>
+               ir_weights,            // DeviceTensor
+               use_sum_factorization, // bool
+               input_dtq_maps,        // std::array<DofToQuadMap, num_fields>
+               output_dtq_maps,       // std::array<DofToQuadMap, num_fields>
+               // output_fop,            // class derived from FieldOperator
+               thread_blocks,         // ThreadBlocks
+               shmem_tr_cache,           // Vector (local)
+               shmem_tr_info,            // SharedMemoryInfo
+               // TODO: make this Array<int> a member of the DifferentiableOperator
+               //       and capture it by ref.
+               elem_attributes,       // Array<int>
+
+               input_is_dependent,
+               direction_tr,             // FieldDescriptor
+               direction_tr_e,           // Vector
+               derivative_action_tr_e,   // Vector
+               element_dof_ordering,  // ElementDofOrdering
+               inputs_trial_op_dim,
+               total_trial_op_dim,
+               trial_vdim,
+               input_restriction_transpose,
+               qpdc_idx = cache_index,
+               // capture by ref:
+               &qpdc_mem = derivative_qp_caches_ref
+            ](
+               std::vector<Vector> &f_e, const Vector &dir_tr_l,
+               Vector &derivative_action_tr_l) mutable
+         {
+            restriction<entity_t>(direction_tr, dir_tr_l, direction_tr_e,
+                                  element_dof_ordering);
+            auto ye = Reshape(derivative_action_tr_e.ReadWrite(), num_trial_dof,
+                              trial_vdim, num_entities);
+            auto wrapped_fields_e = wrap_fields(f_e, shmem_tr_info.field_sizes,
+                                                num_entities);
+            auto wrapped_direction_e = Reshape(direction_tr_e.ReadWrite(),
+                                               shmem_tr_info.direction_size,
+                                               num_entities);
+
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
+                                trial_vdim, total_trial_op_dim, num_qp, num_entities);
+
+            auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
+
+            const bool has_attr = attributes.Size() > 0;
+            const auto d_attr = attributes.Read();
+            const auto d_elem_attr = elem_attributes->Read();
+
+            derivative_action_tr_e = 0.0;
+            forall([=] MFEM_HOST_DEVICE (int e, real_t *shmem)
+            {
+               if (has_attr && !d_attr[d_elem_attr[e] - 1]) { return; }
+
+               auto [output_dtq_shmem, input_dtq_shmem, fields_shmem,
+                                       direction_shmem, input_shmem,
+                                       shadow_shmem_, residual_shmem,
+                                       scratch_shmem] =
+                        unpack_shmem(shmem, shmem_tr_info, output_dtq_maps, input_dtq_maps,
+                                     wrapped_fields_e, wrapped_direction_e, num_qp, e);
+               auto &shadow_shmem = shadow_shmem_;
+
+               std::array<bool, num_outputs> all_true{true};
+               map_direction_to_quadrature_data_conditional(
+                  shadow_shmem, direction_shmem, output_dtq_shmem, outputs,
+                  ir_weights, scratch_shmem, all_true, dimension,
+                  use_sum_factorization);
+
+               auto fhat = Reshape(&residual_shmem(0, 0), trial_vdim,
+                                   total_trial_op_dim, num_qp);
+
+               auto qpdce = Reshape(&qpdc(0, 0, 0, 0, 0, e), test_vdim, test_op_dim,
+                                    trial_vdim, total_trial_op_dim, num_qp);
+
+               constexpr bool transpose = true;
+               apply_qpdc(fhat, shadow_shmem, qpdce, itod, q1d, dimension,
+                          use_sum_factorization, transpose);
+
+               auto y = Reshape(&ye(0, 0, e), num_trial_dof, trial_vdim);
+               auto fi_shmem = Reshape(shmem + shmem_tr_info.total_size, trial_vdim,
+                                       total_trial_op_dim, num_qp);
+
+               map_quadrature_data_to_fields_conditional(
+                  y, fhat, inputs, itod, input_dtq_shmem, scratch_shmem, fi_shmem,
+                  input_is_dependent, dimension, use_sum_factorization);
+            }, num_entities, thread_blocks, shmem_tr_info.total_size,
+            shmem_tr_cache.ReadWrite());
+            input_restriction_transpose(derivative_action_tr_e, derivative_action_tr_l);
+         });
+
          assemble_derivative_sparsematrix_callbacks[derivative_id].push_back(
             [
                // capture by copy:
@@ -1190,7 +1359,7 @@ void DifferentiableOperator::AddIntegrator(
                inputs_trial_op_dim,
                Ae_mem,
                output_to_field,
-
+               qpdc_idx = cache_index,
                // capture by ref:
                &qpdc_mem = derivative_qp_caches_ref,
                &fields = fields_ref
@@ -1202,7 +1371,7 @@ void DifferentiableOperator::AddIntegrator(
                                                shmem_info.direction_size,
                                                num_entities);
 
-            auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
                                 trial_vdim, total_trial_op_dim, num_qp, num_entities);
 
             auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1250,7 +1419,10 @@ void DifferentiableOperator::AddIntegrator(
             auto test_fes = *std::get_if<const ParFiniteElementSpace *>
                             (&fields[output_to_field[0]].data);
 
-            A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
+            if (A == nullptr)
+            {
+               A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
+            }
 
             auto tmp = Reshape(Ae_mem.HostReadWrite(), num_test_dof * test_vdim,
                                num_trial_dof * trial_vdim, num_entities);
@@ -1321,14 +1493,15 @@ void DifferentiableOperator::AddIntegrator(
                   A->AddSubMatrix(test_vdofs, trial_vdofs, Aee, 1);
                }
             }
-            A->Finalize();
+            // Don't finalize here since multiple callbacks might contribute to the same matrix
+            // A->Finalize() will be called after all callbacks have contributed
          });
 
          // Create local references for MSVC lambda capture compatibility
          auto& assemble_derivative_sparsematrix_callbacks_ref =
             this->assemble_derivative_sparsematrix_callbacks[derivative_id];
 
-         assemble_derivative_hypreparmatrix_callbacks[derivative_id].push_back(
+         assemble_derivative_hypreparmatrix_callback[derivative_id] =
             [
                input_is_dependent,
                input_to_field,
@@ -1342,6 +1515,7 @@ void DifferentiableOperator::AddIntegrator(
             {
                f(f_e, spmat);
             }
+            spmat->Finalize();
 
             if (spmat == nullptr)
             {
@@ -1395,7 +1569,7 @@ void DifferentiableOperator::AddIntegrator(
                        trial_fes->Dof_TrueDof_Matrix());
             }
             delete spmat;
-         });
+         };
       }, derivative_ids);
    }
 }

fem/dfem/integrate.hpp

@@ -530,4 +530,65 @@ void map_quadrature_data_to_fields(
    }
 }
 
+template <size_t N, typename field_operator_ts>
+MFEM_HOST_DEVICE
+void map_quadrature_data_to_fields_conditional(
+   DeviceTensor<2, real_t> &y,
+   const DeviceTensor<3, real_t> &f,
+   const field_operator_ts &fops,
+   const DeviceTensor<1, const real_t> &op_dims,
+   const std::array<DofToQuadMap, N> &dtqmaps,
+   std::array<DeviceTensor<1>, 6> &scratch_mem,
+   const DeviceTensor<3> &fi_shmem,
+   const std::array<bool, N> &conditions,
+   const int &dimension,
+   const bool &use_sum_factorization)
+{
+   int offset = 0;
+   for_constexpr<N>([&](auto i)
+   {
+      if (conditions[i])
+      {
+         [[maybe_unused]] const auto [K, unused, M] = f.GetShape();
+         const int L = static_cast<int>(op_dims(static_cast<size_t>(i)));
+         auto fi = Reshape(&fi_shmem(0, 0, 0), K, L, M);
+         for (int k = 0; k < K; k++)
+         {
+            for (int l = 0; l < L; l++)
+            {
+               for (int m = 0; m < M; m++)
+               {
+                  fi(k, l, m) = f(k, l + offset, m);
+               }
+            }
+         }
+
+         if (use_sum_factorization)
+         {
+            if (dimension == 1)
+            {
+               map_quadrature_data_to_fields_tensor_impl_1d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else if (dimension == 2)
+            {
+               map_quadrature_data_to_fields_tensor_impl_2d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else if (dimension == 3)
+            {
+               map_quadrature_data_to_fields_tensor_impl_3d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else { MFEM_ABORT_KERNEL("dimension not supported"); }
+         }
+         else
+         {
+            map_quadrature_data_to_fields_impl(y, fi, get<i>(fops), dtqmaps[i]);
+         }
+         offset += L;
+      }
+   });
+}
+
 } // namespace mfem::future

fem/dfem/interpolate.hpp

@@ -505,13 +505,13 @@ void map_field_to_quadrature_data(
    }
 }
 
-template <typename field_operator_ts, size_t num_inputs, size_t num_fields>
+template <typename field_operator_ts, size_t N, size_t M>
 MFEM_HOST_DEVICE inline
 void map_fields_to_quadrature_data(
-   std::array<DeviceTensor<2>, num_inputs> &fields_qp,
-   const std::array<DeviceTensor<1>, num_fields> &fields_e,
-   const std::array<DofToQuadMap, num_inputs> &dtqmaps,
-   const std::array<size_t, num_inputs> &input_to_field,
+   std::array<DeviceTensor<2>, N> &fields_qp,
+   const std::array<DeviceTensor<1>, M> &fields_e,
+   const std::array<DofToQuadMap, N> &dtqmaps,
+   const std::array<size_t, N> &input_to_field,
    const field_operator_ts &fops,
    const DeviceTensor<1, const real_t> &integration_weights,
    const std::array<DeviceTensor<1>, 6> &scratch_mem,
@@ -523,7 +523,7 @@ void map_fields_to_quadrature_data(
    // attached to them and we create a dummy field which is not accessed
    // inside the functions it is passed to.
    const auto dummy_field_weight = DeviceTensor<1>(nullptr, 0);
-   for_constexpr<num_inputs>([&](auto i)
+   for_constexpr<N>([&](auto i)
    {
       const DeviceTensor<1> &field_e =
          (input_to_field[i] == SIZE_MAX) ? dummy_field_weight :
@@ -549,12 +549,7 @@ void map_fields_to_quadrature_data(
                fields_qp[i], dtqmaps[i], field_e, get<i>(fops),
                integration_weights, scratch_mem);
          }
-         else
-         {
-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
-            MFEM_ABORT("unsupported dimension");
-#endif
-         }
+         else { MFEM_ABORT_KERNEL("unsupported dimension"); }
       }
       else
       {
@@ -627,20 +622,20 @@ void map_fields_to_quadrature_data_conditional(
    });
 }
 
-template <size_t num_inputs, typename field_operator_ts>
+template <size_t N, typename field_operator_ts>
 MFEM_HOST_DEVICE
 void map_direction_to_quadrature_data_conditional(
-   std::array<DeviceTensor<2>, num_inputs> &directions_qp,
+   std::array<DeviceTensor<2>, N> &directions_qp,
    const DeviceTensor<1> &direction_e,
-   const std::array<DofToQuadMap, num_inputs> &dtqmaps,
+   const std::array<DofToQuadMap, N> &dtqmaps,
    field_operator_ts fops,
    const DeviceTensor<1, const real_t> &integration_weights,
    const std::array<DeviceTensor<1>, 6> &scratch_mem,
-   const std::array<bool, num_inputs> &conditions,
+   const std::array<bool, N> &conditions,
    const int &dimension,
    const bool &use_sum_factorization)
 {
-   for_constexpr<num_inputs>([&](auto i)
+   for_constexpr<N>([&](auto i)
    {
       if (conditions[i])
       {

fem/dfem/qfunction_apply.hpp

@@ -379,54 +379,102 @@ namespace detail
 /// @param shadow_shmem the shadow shared memory.
 /// @param qpdc the quadrature point data cache holding the resulting
 /// Jacobians on each quadrature point.
-/// @param itod inputs trial operator dimension.
-/// If input is dependent the value corresponds to the spatial dimension, otherwise
-/// a zero indicates non-dependence on the variable.
+/// @param op_dims operator dimensions.
+/// If an operator is dependent, the value corresponds to the spatial dimension.
+/// Otherwise a zero indicates indepence on the variable.
 /// @param q the current quadrature point index.
-template <size_t num_fields>
+/// @param transpose switch to use transpose action.
+template <size_t N>
 MFEM_HOST_DEVICE inline
 void apply_qpdc(
    DeviceTensor<3> &fhat,
-   const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
+   const std::array<DeviceTensor<2>, N> &shadow_shmem,
    const DeviceTensor<5, const real_t> &qpdc,
-   const DeviceTensor<1, const real_t> &itod,
-   const int &q)
+   const DeviceTensor<1, const real_t> &op_dims,
+   const int &q,
+   bool transpose)
 {
+   const size_t num_ops = op_dims.GetShape()[0];
+
    const int test_vdim = qpdc.GetShape()[0];
    const int test_op_dim = qpdc.GetShape()[1];
    const int trial_vdim = qpdc.GetShape()[2];
-   const int num_qp = qpdc.GetShape()[4];
-   const size_t num_inputs = itod.GetShape()[0];
+   const int total_trial_op_dim = qpdc.GetShape()[3];
 
-   for (int i = 0; i < test_vdim; i++)
+   const int num_qp = qpdc.GetShape()[4];
+
+   if (transpose)
    {
-      for (int k = 0; k < test_op_dim; k++)
+      for (int j = 0; j < trial_vdim; j++)
       {
-         real_t sum = 0.0;
-         int m_offset = 0;
-         for (size_t s = 0; s < num_inputs; s++)
+         for (int m = 0; m < total_trial_op_dim; m++)
          {
-            const int trial_op_dim = static_cast<int>(itod(s));
-            if (trial_op_dim == 0)
-            {
-               continue;
-            }
-            const auto d_qp =
-               Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
-            for (int j = 0; j < trial_vdim; j++)
-            {
-               for (int m = 0; m < trial_op_dim; m++)
-               {
-                  sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
-               }
-            }
-            m_offset += trial_op_dim;
+            fhat(j, m, q) = 0.0;
+         }
+      }
+
+      // Since we don't support more than output space right now
+      // shadow_shmem will always be of size 1.
+      constexpr int shadow_idx_tr = 0;
+      auto d_qp = Reshape(&(shadow_shmem[shadow_idx_tr])[0], test_vdim, test_op_dim,
+                          num_qp);
+
+      int m_offset = 0;
+      for (size_t s = 0; s < num_ops; s++)
+      {
+         const int trial_op_dim = static_cast<int>(op_dims(s));
+         if (trial_op_dim == 0) { continue; }
+
+         for (int j = 0; j < trial_vdim; j++)
+         {
+            for (int m = 0; m < trial_op_dim; m++)
+            {
+               real_t sum = 0.0;
+               for (int i = 0; i < test_vdim; i++)
+               {
+                  for (int k = 0; k < test_op_dim; k++)
+                  {
+                     const real_t contrib = qpdc(i, k, j, m + m_offset, q) * d_qp(i, k, q);
+                     sum += contrib;
+                  }
+               }
+               fhat(j, m + m_offset, q) += sum;
+            }
+         }
+         m_offset += trial_op_dim;
+      }
+   }
+   else
+   {
+      for (int i = 0; i < test_vdim; i++)
+      {
+         for (int k = 0; k < test_op_dim; k++)
+         {
+            real_t sum = 0.0;
+            int m_offset = 0;
+            for (size_t s = 0; s < num_ops; s++)
+            {
+               const int trial_op_dim = static_cast<int>(op_dims(s));
+               if (trial_op_dim == 0) { continue; }
+
+               const auto d_qp =
+                  Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
+               for (int j = 0; j < trial_vdim; j++)
+               {
+                  for (int m = 0; m < trial_op_dim; m++)
+                  {
+                     sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
+                  }
+               }
+               m_offset += trial_op_dim;
+            }
+            fhat(i, k, q) = sum;
          }
-         fhat(i, k, q) = sum;
       }
    }
 }
-}
+
+} // namespace detail
 
 /// @brief Apply the quadrature point data cache (qpdc) to a vector
 /// (usually a direction).
@@ -445,16 +493,18 @@ void apply_qpdc(
 /// @param q1d number of quadrature points in 1D.
 /// @param dimension spatial dimension.
 /// @param use_sum_factorization whether to use sum factorization.
-template <size_t num_fields>
+/// @param T switch to use transpose application.
+template <size_t N>
 MFEM_HOST_DEVICE inline
 void apply_qpdc(
    DeviceTensor<3> &fhat,
-   const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
+   const std::array<DeviceTensor<2>, N> &shadow_shmem,
    const DeviceTensor<5, const real_t> &qpdc,
    const DeviceTensor<1, const real_t> &itod,
    const int &q1d,
    const int &dimension,
-   const bool &use_sum_factorization)
+   const bool &use_sum_factorization,
+   const bool T = false)
 {
    if (use_sum_factorization)
    {
@@ -462,7 +512,7 @@ void apply_qpdc(
       {
          MFEM_FOREACH_THREAD_DIRECT(q, x, q1d)
          {
-            detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+            detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
          }
       }
       else if (dimension == 2)
@@ -472,7 +522,7 @@ void apply_qpdc(
             MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
             {
                const int q = qx + q1d * qy;
-               detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+               detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
             }
          }
       }
@@ -485,7 +535,7 @@ void apply_qpdc(
                MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
                {
                   const int q = qx + q1d * (qy + q1d * qz);
-                  detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+                  detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
                }
             }
          }
@@ -500,7 +550,7 @@ void apply_qpdc(
       const int num_qp = qpdc.GetShape()[4];
       MFEM_FOREACH_THREAD_DIRECT(q, x, num_qp)
       {
-         detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+         detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
       }
    }
 }

fem/dfem/qfunction_transform.hpp

@@ -243,6 +243,35 @@ void process_qf_arg(
    }
 }
 
+template <typename T, int n>
+MFEM_HOST_DEVICE inline
+void process_qf_arg(
+   const DeviceTensor<1, T> &u,
+   const DeviceTensor<1, T> &v,
+   tensor<T, n> &arg)
+{
+   for (int i = 0; i < n; i++)
+   {
+      arg(i) = u(i);
+   }
+}
+
+template <typename T, int n, int m>
+MFEM_HOST_DEVICE inline
+void process_qf_arg(
+   const DeviceTensor<1, T> &u,
+   const DeviceTensor<1, T> &v,
+   tensor<T, n, m> &arg)
+{
+   for (int i = 0; i < m; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         arg(j, i) = u((i * n) + j);
+      }
+   }
+}
+
 template <typename arg_type>
 MFEM_HOST_DEVICE inline
 void process_qf_arg(const DeviceTensor<2> &u, arg_type &arg, int qp)
@@ -327,20 +356,4 @@ void process_qf_result(
    }
 }
 
-template <typename T, int n, int m>
-MFEM_HOST_DEVICE inline
-void process_qf_arg(
-   const DeviceTensor<1, T> &u,
-   const DeviceTensor<1, T> &v,
-   tensor<T, n, m> &arg)
-{
-   for (int i = 0; i < m; i++)
-   {
-      for (int j = 0; j < n; j++)
-      {
-         arg(j, i) = u((i * n) + j);
-      }
-   }
-}
-
 } // namespace mfem::future

fem/dfem/util.hpp

@@ -1202,7 +1202,14 @@ std::function<void(const Vector&, Vector&)> get_prolongation_transpose(
    const Operator *P = get_prolongation(f);
    auto PT = [=](const Vector &r_local, Vector &y)
    {
-      P->MultTranspose(r_local, y);
+      if (P)
+      {
+         P->MultTranspose(r_local, y);
+      }
+      else
+      {
+         y = r_local;
+      }
    };
    return PT;
 }
@@ -1580,14 +1587,13 @@ struct SharedMemoryInfo
    std::array<int, 6> temp_sizes;
 };
 
-template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs, typename input_t>
+template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs>
 SharedMemoryInfo<num_fields, num_inputs, num_outputs>
 get_shmem_info(
    const std::array<DofToQuadMap, num_inputs> &input_dtq_maps,
    const std::array<DofToQuadMap, num_outputs> &output_dtq_maps,
    const std::vector<FieldDescriptor> &fields,
    const int &num_entities,
-   const input_t &inputs,
    const int &num_qp,
    const std::vector<int> &input_size_on_qp,
    const int &residual_size_on_qp,

fem/gridfunc.cpp

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

fem/gridfunc.hpp

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

fem/integ/bilininteg_vecdiffusion_pa.hpp

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

fem/integ/bilininteg_vecmass_pa.hpp

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

fem/integ/lininteg_domain_kernels.hpp

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

fem/pgridfunc.cpp

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

fem/pgridfunc.hpp

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

fem/qinterp/det.hpp

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

fem/qinterp/eval.hpp

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

fem/qinterp/eval_hdiv.hpp

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

fem/qinterp/grad.hpp

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

fem/quadinterpolator.cpp

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

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/CMakeLists.txt

@@ -39,6 +39,8 @@ set(UNIT_TESTS_SRCS
   dfem/test_divergence.cpp
   dfem/test_lvector_interface.cpp
   dfem/test_mass.cpp
+  dfem/test_transpose.cpp
+  dfem/test_multiple_integrators.cpp
   general/test_array.cpp
   general/test_scan.cpp
   general/test_arrays_by_name.cpp

tests/unit/dfem/test_mass.cpp

@@ -218,7 +218,7 @@ template <int DIM> void mass_mat_mixed(const char* filename, int p)
 
    SECTION("spmat")
    {
-      SparseMatrix *A;
+      SparseMatrix *A = nullptr;
       ddopdu->Assemble(A);
       TestSameMatrices(*A, blf.SpMat());
       delete A;

tests/unit/dfem/test_multiple_integrators.cpp

@@ -0,0 +1,202 @@
+// 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.
+
+#include "../unit_tests.hpp"
+#include "../linalg/test_same_matrices.hpp"
+#include "mfem.hpp"
+#include "fem/dfem/doperator.hpp"
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+
+#ifdef MFEM_USE_ENZYME
+using dscalar_t = real_t;
+#else
+using mfem::future::dual;
+using dscalar_t = dual<real_t, real_t>;
+#endif
+
+template <int DIM>
+void mult_integ(const char *filename, int p)
+{
+   CAPTURE(filename, DIM, p);
+
+   Mesh smesh(filename);
+   ParMesh pmesh(MPI_COMM_WORLD, smesh);
+   pmesh.EnsureNodes();
+   auto* nodes = static_cast<ParGridFunction*>(pmesh.GetNodes());
+   p = std::max(p, pmesh.GetNodalFESpace()->GetMaxElementOrder());
+   smesh.Clear();
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace fes(&pmesh, &fec);
+
+   ParGridFunction x(&fes), y(&fes), z(&fes);
+   Vector X(fes.GetTrueVSize()), Y(fes.GetTrueVSize()), Z(fes.GetTrueVSize());
+
+   X.Randomize(1);
+   x.SetFromTrueDofs(X);
+
+   const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(), 2 * p);
+
+   Array<int> all_domain_attr;
+   if (pmesh.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(pmesh.attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   ParBilinearForm blf(&fes);
+   blf.AddDomainIntegrator(new DiffusionIntegrator(ir));
+   blf.AddDomainIntegrator(new MassIntegrator(ir));
+   blf.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   blf.Assemble();
+   blf.Mult(x, y);
+   fes.GetProlongationMatrix()->MultTranspose(y, Y);
+
+   const auto mass_qf =
+      [] MFEM_HOST_DEVICE(
+         const dscalar_t &u,
+         const tensor<real_t, DIM, DIM> &J,
+         const real_t &w)
+   {
+      return tuple{u * w * det(J)};
+   };
+
+   const auto diffusion_qf =
+      [] MFEM_HOST_DEVICE(
+         const tensor<dscalar_t, DIM> &dudxi,
+         const tensor<real_t, DIM, DIM> &J,
+         const real_t &w)
+   {
+      return tuple{(dudxi * inv(J)) * transpose(inv(J)) * w * det(J)};
+   };
+
+   static constexpr int U = 0, Coords = 1;
+   const auto sol = std::vector{ FieldDescriptor{ U, &fes } };
+   DifferentiableOperator dop(sol, {{Coords, nodes->ParFESpace()}}, pmesh);
+
+   auto derivatives = std::integer_sequence<size_t, U> {};
+
+   dop.AddDomainIntegrator(diffusion_qf,
+                           tuple{ Gradient<U>{}, Gradient<Coords>{}, Weight{} },
+                           tuple{ Gradient<U>{} },
+                           *ir, all_domain_attr, derivatives);
+
+   dop.AddDomainIntegrator(mass_qf,
+                           tuple{ Value<U>{}, Gradient<Coords>{}, Weight{} },
+                           tuple{ Value<U>{} },
+                           *ir, all_domain_attr, derivatives);
+
+   SECTION("action")
+   {
+      dop.SetParameters({ nodes });
+
+      fes.GetRestrictionMatrix()->Mult(x, X);
+      dop.Mult(X, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   SECTION("linearized action")
+   {
+      auto ddopdu = dop.GetDerivative(U, {&x}, {nodes});
+
+      fes.GetRestrictionMatrix()->Mult(x, X);
+      ddopdu->Mult(X, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   SECTION("linearized assembled SparseMatrix")
+   {
+      auto ddopdu = dop.GetDerivative(U, {&x}, {nodes});
+
+      SparseMatrix *A = nullptr;
+      ddopdu->Assemble(A);
+
+      A->Mult(x, z);
+      fes.GetProlongationMatrix()->MultTranspose(z, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+      delete A;
+   }
+
+   SECTION("linearized assembled HypreParMatrix")
+   {
+      auto ddopdu = dop.GetDerivative(U, {&x}, {nodes});
+
+      HypreParMatrix *A = nullptr;
+      ddopdu->Assemble(A);
+
+      fes.GetRestrictionMatrix()->Mult(x, X);
+      A->Mult(X, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+      delete A;
+   }
+
+}
+
+// no GPU tag to avoid failing 'hypre parallel mat' section
+TEST_CASE("dFEM Multiple Integrators", "[Parallel][dFEM][XXX]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+
+   const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
+
+   SECTION("2d")
+   {
+      const auto filename2d =
+         GENERATE(
+            "../../data/star.mesh",
+            "../../data/star-q3.mesh",
+            "../../data/rt-2d-q3.mesh",
+            "../../data/inline-quad.mesh",
+            "../../data/periodic-square.mesh"
+         );
+      mult_integ<2>(filename2d, p);
+   }
+
+   SECTION("3d")
+   {
+      const auto filename3d =
+         GENERATE(
+            "../../data/fichera.mesh",
+            "../../data/fichera-q3.mesh",
+            "../../data/inline-hex.mesh",
+            "../../data/toroid-hex.mesh",
+            "../../data/periodic-cube.mesh"
+         );
+      mult_integ<3>(filename3d, p);
+   }
+}
+
+#endif // MFEM_USE_MPI

tests/unit/dfem/test_transpose.cpp

@@ -0,0 +1,409 @@
+// 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.
+
+#include "../unit_tests.hpp"
+#include "mfem.hpp"
+
+#include "fem/dfem/doperator.hpp" // TODO: remove before merge
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+
+#ifdef MFEM_USE_ENZYME
+using dscalar_t = real_t;
+#else
+using mfem::future::dual;
+using dscalar_t = dual<real_t, real_t>;
+#endif
+
+template <int DIM>
+void transpose(const char *filename, int p)
+{
+   CAPTURE(filename, DIM, p);
+
+   Mesh serial_mesh(filename);
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   mesh.EnsureNodes();
+   auto* nodes = static_cast<ParGridFunction*>(mesh.GetNodes());
+   p = std::max(p, mesh.GetNodalFESpace()->GetMaxElementOrder());
+
+   Array<int> all_domain_attr;
+   if (mesh.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(mesh.attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace scalar_fes(&mesh, &fec);
+   ParFiniteElementSpace vector_fes(&mesh, &fec, DIM);
+
+   ParGridFunction sgf(&scalar_fes);
+
+   auto f0 = [](const Vector &x)
+   {
+      if constexpr (DIM == 3)
+      {
+         return M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]) * sin(M_PI*x[2]);
+      }
+      return M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]);
+   };
+
+   FunctionCoefficient f0_coeff(f0);
+   sgf.ProjectCoefficient(f0_coeff);
+
+   ParGridFunction vgf(&vector_fes);
+
+   auto gradf1 = [](const Vector &x, Vector &u)
+   {
+      if constexpr (DIM == 3)
+      {
+         u(0) = M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]) * sin(M_PI*x[2]);
+         u(1) = M_PI*sin(M_PI*x[0]) * cos(M_PI*x[1]) * sin(M_PI*x[2]);
+         u(2) = M_PI*sin(M_PI*x[0]) * sin(M_PI*x[1]) * cos(M_PI*x[2]);
+         return;
+      }
+      u(0) = M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]);
+      u(1) = M_PI*sin(M_PI*x[0]) * cos(M_PI*x[1]);
+   };
+
+   VectorFunctionCoefficient gradf1_coeff(DIM, gradf1);
+   vgf.ProjectCoefficient(gradf1_coeff);
+
+   const auto* ir = &IntRules.Get(mesh.GetTypicalElementGeometry(), 2 * p);
+
+   SECTION("Mass Transpose Action")
+   {
+      ParBilinearForm Mblf(&scalar_fes);
+      auto mass_integ = new MassIntegrator;
+      mass_integ->SetIntegrationRule(*ir);
+      Mblf.AddDomainIntegrator(mass_integ);
+      Mblf.Assemble();
+      Mblf.Finalize();
+      auto Mmat = Mblf.ParallelAssemble();
+
+      static constexpr int SCALAR = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{SCALAR, &scalar_fes}};
+      const auto par = std::vector{FieldDescriptor{COORDINATES, nodes->ParFESpace()}};
+      DifferentiableOperator dop(sol, par, mesh);
+      const auto gradient_qf = [] MFEM_HOST_DEVICE(
+                                  const dscalar_t &u,
+                                  const tensor<real_t, DIM, DIM> &J,
+                                  const real_t &w)
+      {
+         return tuple{u * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, SCALAR> {};
+      dop.AddDomainIntegrator(gradient_qf,
+                              tuple{Value<SCALAR>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<SCALAR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      Vector S, T, U;
+      S.SetSize(scalar_fes.GetTrueVSize());
+      T.SetSize(scalar_fes.GetTrueVSize());
+      U.SetSize(scalar_fes.GetTrueVSize());
+
+      sgf.GetTrueDofs(S);
+
+      Mmat->MultTranspose(S, T);
+
+      auto ddop = dop.GetDerivative(SCALAR, {&sgf}, {nodes});
+      ddop->MultTranspose(S, U);
+
+      T -= U;
+      real_t norm_g, norm_l = T.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+
+      delete Mmat;
+   }
+
+   SECTION("Vector Mass Transpose Action")
+   {
+      ParBilinearForm Mvblf(&vector_fes);
+      auto mass_integ = new VectorMassIntegrator;
+      mass_integ->SetIntegrationRule(*ir);
+      Mvblf.AddDomainIntegrator(mass_integ);
+      Mvblf.Assemble();
+      Mvblf.Finalize();
+      auto Mvmat = Mvblf.ParallelAssemble();
+
+      static constexpr int VECTOR = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{VECTOR, &vector_fes}};
+      const auto par = std::vector{FieldDescriptor{COORDINATES, nodes->ParFESpace()}};
+      DifferentiableOperator dop(sol, par, mesh);
+      const auto gradient_qf = [] MFEM_HOST_DEVICE(
+                                  const tensor<dscalar_t, DIM> &u,
+                                  const tensor<real_t, DIM, DIM> &J,
+                                  const real_t &w)
+      {
+         return tuple{u * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, VECTOR> {};
+      dop.AddDomainIntegrator(gradient_qf,
+                              tuple{Value<VECTOR>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<VECTOR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      Vector V, W, Z;
+      V.SetSize(vector_fes.GetTrueVSize());
+      W.SetSize(vector_fes.GetTrueVSize());
+      Z.SetSize(vector_fes.GetTrueVSize());
+
+      vgf.GetTrueDofs(V);
+
+      Mvmat->MultTranspose(V, W);
+
+      auto ddop = dop.GetDerivative(VECTOR, {&vgf}, {nodes});
+      ddop->MultTranspose(V, Z);
+
+      W -= Z;
+      real_t norm_g, norm_l = W.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+
+      delete Mvmat;
+   }
+
+   SECTION("Discrete Gradient Transpose Action")
+   {
+      ParMixedBilinearForm Gblf(&scalar_fes, &vector_fes);
+      auto grad_integ = new GradientIntegrator;
+      grad_integ->SetIntegrationRule(*ir);
+      Gblf.AddDomainIntegrator(grad_integ);
+      Gblf.Assemble();
+      Gblf.Finalize();
+      auto Gmat = Gblf.ParallelAssemble();
+
+      static constexpr int SCALAR = 0, VECTOR = 2, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{SCALAR, &scalar_fes}};
+      const auto par = std::vector
+      {
+         FieldDescriptor{VECTOR, &vector_fes},
+         FieldDescriptor{COORDINATES, nodes->ParFESpace()}
+      };
+      DifferentiableOperator dop(sol, par, mesh);
+      const auto gradient_qf = [] MFEM_HOST_DEVICE(
+                                  const tensor<dscalar_t, DIM> &dudxi,
+                                  const tensor<real_t, DIM, DIM> &J,
+                                  const real_t &w)
+      {
+         const auto dudx = dudxi * inv(J);
+         return tuple{dudx * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, SCALAR> {};
+      dop.AddDomainIntegrator(gradient_qf,
+                              tuple{Gradient<SCALAR>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<VECTOR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({&vgf, nodes});
+
+      Vector S, T, V;
+      S.SetSize(scalar_fes.GetTrueVSize());
+      T.SetSize(scalar_fes.GetTrueVSize());
+      vgf.GetTrueDofs(V);
+
+      Gmat->MultTranspose(V, S);
+
+      auto ddop = dop.GetDerivative(SCALAR, {&sgf}, {&vgf, nodes});
+      ddop->MultTranspose(V, T);
+
+      S -= T;
+      real_t norm_g, norm_l = S.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+
+      delete Gmat;
+   }
+
+   SECTION("Scalar Convection Transpose Action")
+   {
+      auto b_func = [](const Vector &x, Vector &b)
+      {
+         b(0) = cos(x[0] * 2.0 * M_PI);
+         b(1) = 1.0 + cos(x[1] * 2.0 * M_PI);
+         if constexpr (DIM == 3)
+         {
+            b(2) = 2.0 + cos(x[2] * 2.0 * M_PI);
+         }
+      };
+      VectorFunctionCoefficient b_coeff(DIM, b_func);
+
+      ParBilinearForm Gblf(&scalar_fes);
+      auto conv_integ = new ConvectionIntegrator(b_coeff);
+      conv_integ->SetIntegrationRule(*ir);
+      Gblf.AddDomainIntegrator(conv_integ);
+      Gblf.Assemble();
+      Gblf.Finalize();
+      auto Gmat = Gblf.ParallelAssemble();
+
+      static constexpr int SCALAR = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{SCALAR, &scalar_fes}};
+      const auto par = std::vector
+      {
+         FieldDescriptor{COORDINATES, nodes->ParFESpace()}
+      };
+      DifferentiableOperator dop(sol, par, mesh);
+
+      const auto convection_qf =
+         [] MFEM_HOST_DEVICE(
+            const tensor<dscalar_t, DIM> &dudxi,
+            const tensor<real_t, DIM> &x,
+            const tensor<real_t, DIM, DIM> &J,
+            const real_t &w)
+      {
+         const auto dudx = dudxi * inv(J);
+         tensor<dscalar_t, DIM> b{};
+         b(0) = cos(x[0] * 2.0 * M_PI);
+         b(1) = 1.0 + cos(x[1] * 2.0 * M_PI);
+         if constexpr (DIM == 3)
+         {
+            b(2) = 2.0 + cos(x[2] * 2.0 * M_PI);
+         }
+         return tuple{dot(b, dudx) * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, SCALAR> {};
+      dop.AddDomainIntegrator(convection_qf,
+                              tuple{Gradient<SCALAR>{}, Value<COORDINATES>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<SCALAR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      Vector S, T, U;
+      S.SetSize(scalar_fes.GetTrueVSize());
+      T.SetSize(scalar_fes.GetTrueVSize());
+      U.SetSize(scalar_fes.GetTrueVSize());
+      U.Randomize(1);
+
+      {
+         Gmat->MultTranspose(U, S);
+
+         auto ddop = dop.GetDerivative(SCALAR, {&sgf}, {nodes});
+         ddop->MultTranspose(U, T);
+
+         S -= T;
+         real_t norm_g, norm_l = S.Normlinf();
+         MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+         REQUIRE(norm_g == MFEM_Approx(0.0));
+      }
+
+      delete Gmat;
+   }
+
+   SECTION("Nonlinear VectorConvection Transpose Action")
+   {
+      auto b_func = [](const Vector &x, Vector &b)
+      {
+         b(0) = cos(x[0]) * sin(x[0]) * x[1];
+         b(1) = cos(x[1]) * sin(x[1]) * x[0];
+         if constexpr (DIM == 3)
+         {
+            b(2) = cos(x[2]) * sin(x[2]) * x[0];
+         }
+      };
+      VectorFunctionCoefficient b_coeff(DIM, b_func);
+
+      ParGridFunction ugf(&vector_fes);
+      ugf.ProjectCoefficient(b_coeff);
+
+      Vector U(vector_fes.GetTrueVSize());
+      ugf.GetTrueDofs(U);
+
+      ParNonlinearForm nlf(&vector_fes);
+      const auto vcinteg = new VectorConvectionNLFIntegrator();
+      vcinteg->SetIntegrationRule(*ir);
+      nlf.AddDomainIntegrator(vcinteg);
+      HypreParMatrix &Nmat = dynamic_cast<HypreParMatrix&>(nlf.GetGradient(U));
+
+      static constexpr int VELOCITY = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{VELOCITY, &vector_fes}};
+      const auto par = std::vector
+      {
+         FieldDescriptor{COORDINATES, nodes->ParFESpace()}
+      };
+      DifferentiableOperator dop(sol, par, mesh);
+
+      const auto nlconvection_qf =
+         [] MFEM_HOST_DEVICE(
+            const tensor<dscalar_t, DIM> &u,
+            const tensor<dscalar_t, DIM, DIM> &dudxi,
+            const tensor<real_t, DIM, DIM> &J,
+            const real_t &w)
+      {
+         const auto invJ = inv(J);
+         const auto dudx = dudxi * invJ;
+         return tuple{dot(dudx, u) * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, VELOCITY> {};
+      dop.AddDomainIntegrator(nlconvection_qf,
+                              tuple{Value<VELOCITY>{}, Gradient<VELOCITY>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<VELOCITY>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      auto ddop = dop.GetDerivative(VELOCITY, {&ugf}, {nodes});
+
+      Vector S(U.Size()), T(U.Size()), Se(vector_fes.GetVSize());
+
+      Nmat.MultTranspose(U, S);
+      ddop->MultTranspose(U, T);
+
+      S -= T;
+      real_t norm_g, norm_l = S.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+   }
+}
+
+TEST_CASE("dFEM Transpose", "[Parallel][dFEM][XXX]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+
+   const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
+
+   SECTION("2d")
+   {
+      const auto filename2d =
+         GENERATE(
+            "../../data/star.mesh",
+            "../../data/star-q3.mesh",
+            "../../data/rt-2d-q3.mesh",
+            "../../data/inline-quad.mesh"
+         );
+      transpose<2>(filename2d, p);
+   }
+
+   SECTION("3d")
+   {
+      const auto filename3d =
+         GENERATE(
+            "../../data/fichera.mesh",
+            "../../data/fichera-q3.mesh",
+            "../../data/inline-hex.mesh",
+            "../../data/toroid-hex.mesh"
+         );
+      transpose<3>(filename3d, p);
+   }
+}
+
+#endif