remove print

add derivative call

commit to check compatibility on CI

dual type for non enzyme ci

change dependent variable

add dual transformation

back to non generic lambda for nvcc

check constexpr mechanism with CI

testing out lambda-less approach

cleanup

Revert "testing out lambda-less approach"

This reverts commit 6808fe0693c2589abf3c0c929e0a10b9a5050d4d.

# Conflicts:
#	tests/unit/dfem/test_functional.cpp

simpler approach

fem/dfem/doperator.hpp

@@ -367,9 +367,16 @@ private:
    mutable Vector residual_e;
 
    std::function<void(Vector &, Vector &)> prolongation_transpose;
+
+   // TODO: This can probably be removed, as it is only used in
+   // the callbacks so they can copy it during the capture.
    std::function<void(Vector &, Vector &)> output_restriction_transpose;
+
    restriction_callback_t restriction_callback;
 
+   std::map<size_t, std::function<void(const Vector &, Vector &)>>
+                                                                daction_prolongation_transpose;
+
    std::map<size_t, size_t> assembled_vector_sizes;
 
    bool use_tensor_product_structure = true;
@@ -512,7 +519,7 @@ void DifferentiableOperator::AddDomainIntegrator(
    const int num_entities = GetNumEntities<entity_t>(mesh);
    const int num_qp = integration_rule.GetNPoints();
 
-   if constexpr (is_sum_fop<decltype(output_fop)>::value)
+   if constexpr (is_sum_fop<std::remove_cv_t<decltype(output_fop)>>::value)
    {
       residual_l.SetSize(1);
       height = 1;
@@ -546,8 +553,19 @@ void DifferentiableOperator::AddDomainIntegrator(
 
    const int test_vdim = output_fop.vdim;
    const int test_op_dim = output_fop.size_on_qp / output_fop.vdim;
-   const int num_test_dof =
-      num_entities ? (output_e_size / output_fop.vdim / num_entities) : 0;
+
+   int num_test_dof = 0;
+   if (num_entities)
+   {
+      if constexpr (is_sum_fop<std::decay_t<decltype(output_fop)>>::value)
+      {
+         num_test_dof = 1;
+      }
+      else
+      {
+         num_test_dof = output_e_size / output_fop.vdim / num_entities;
+      }
+   }
 
    auto ir_weights = Reshape(integration_rule.GetWeights().Read(), num_qp);
 
@@ -665,11 +683,16 @@ void DifferentiableOperator::AddDomainIntegrator(
    if constexpr (derivative_ids_t::size() != 0)
    {
       // Create the action of the derivatives
-      for_constexpr([&, &or_transpose =
-                        this->output_restriction_transpose](const std::size_t derivative_id)
+      for_constexpr([&,
+                     &or_transpose = this->output_restriction_transpose,
+                     &dapr_transpose = this->daction_prolongation_transpose]
+                    (const std::size_t derivative_id)
       {
          const size_t d_field_idx = FindIdx(derivative_id, fields);
          const auto direction = fields[d_field_idx];
+
+         dapr_transpose[derivative_id] = get_generic_prolongation_transpose(direction);
+
          const int da_size_on_qp =
             GetSizeOnQP<entity_t>(output_fop, fields[test_space_field_idx]);
 

fem/dfem/integrate.hpp

@@ -225,6 +225,17 @@ void map_quadrature_data_to_fields_tensor_impl_2d(
          MFEM_SYNC_THREAD;
       }
    }
+   else if constexpr (is_sum_fop<std::decay_t<output_t>>::value)
+   {
+      const auto [q1d, unused, d1d] = B.GetShape();
+      auto fqp = Reshape(&f(0, 0, 0), q1d * q1d);
+      auto yqp = Reshape(&y(0, 0), output.size_on_qp);
+      for (int i = 0; i < q1d * q1d; i++)
+      {
+         yqp(0) += fqp(i);
+      }
+      MFEM_SYNC_THREAD;
+   }
    else
    {
       MFEM_ABORT("quadrature data mapping to field is not implemented for"
@@ -411,6 +422,17 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
          MFEM_SYNC_THREAD;
       }
    }
+   else if constexpr (is_sum_fop<std::decay_t<output_t>>::value)
+   {
+      const auto [q1d, unused, d1d] = B.GetShape();
+      auto fqp = Reshape(&f(0, 0, 0), q1d * q1d * q1d);
+      auto yqp = Reshape(&y(0, 0), output.size_on_qp);
+      for (int i = 0; i < q1d * q1d * q1d; i++)
+      {
+         yqp(0) += fqp(i);
+      }
+      MFEM_SYNC_THREAD;
+   }
    else
    {
       MFEM_ABORT("quadrature data mapping to field is not implemented for"

fem/dfem/qfunction_transform.hpp

@@ -181,6 +181,14 @@ void process_derivative_from_native_dual(
    }
 }
 
+template <typename T>
+MFEM_HOST_DEVICE inline
+void process_derivative_from_native_dual(
+   DeviceTensor<1, T> &r,
+   const dual<T, T> &x)
+{
+   r(0) = x.gradient;
+}
 
 template <typename T0, typename T1>
 MFEM_HOST_DEVICE inline

fem/dfem/util.hpp

@@ -892,6 +892,37 @@ int GetDimension(const FieldDescriptor &f)
    }, f.data);
 }
 
+/// @brief Get the number of elements from a field descriptor.
+///
+/// @param f the field descriptor.
+/// @tparam entity_t the entity type (see Entity).
+/// @returns the number of elements from the field descriptor.
+template <typename entity_t>
+int GetNumElements(const FieldDescriptor &f)
+{
+   return std::visit([](auto && arg)
+   {
+      using T = std::decay_t<decltype(arg)>;
+      if constexpr (std::is_same_v<T, const FiniteElementSpace *> ||
+                    std::is_same_v<T, const ParFiniteElementSpace *>)
+      {
+         if constexpr (std::is_same_v<entity_t, Entity::Element>)
+         {
+            return arg->GetMesh()->GetNE();
+         }
+      }
+      else if constexpr (std::is_same_v<T, const ParameterSpace *>)
+      {
+         // TODO: Implement GetNumElements for ParameterSpace
+         return 0;
+      }
+      else
+      {
+         static_assert(dfem::always_false<T>, "can't use GetDimension on type");
+      }
+      return 0; // Unreachable, but avoids compiler warning
+   }, f.data);
+}
 
 /// @brief Get the prolongation operator for a field descriptor.
 ///
@@ -989,7 +1020,9 @@ get_restriction_transpose(
       {
          v_l = v_e;
       };
-      return std::make_tuple(RT, 1);
+      // The size is always 1 * number of elements for sum operators, e.g.
+      // each element is reduced to a single value on the E-vector level.
+      return std::make_tuple(RT, GetNumElements<entity_t>(f));
    }
    else
    {
@@ -1076,6 +1109,18 @@ void prolongation(const std::vector<FieldDescriptor> fields,
    }
 }
 
+inline
+std::function<void(const Vector&, Vector&)> get_generic_prolongation_transpose(
+   const FieldDescriptor &f)
+{
+   const Operator *P = get_prolongation(f);
+   auto PT = [=](const Vector &r_local, Vector &y)
+   {
+      P->MultTranspose(r_local, y);
+   };
+   return PT;
+}
+
 /// @brief Get a transpose prolongation callback for a field descriptor.
 ///
 /// In the special case of a one field operator, the transpose prolongation
@@ -1110,11 +1155,8 @@ std::function<void(const Vector&, Vector&)> get_prolongation_transpose(
       };
       return PT;
    }
-   const Operator *P = get_prolongation(f);
-   auto PT = [=](const Vector &r_local, Vector &y)
-   {
-      P->MultTranspose(r_local, y);
-   };
+
+   auto PT = get_generic_prolongation_transpose(f);
    return PT;
 }
 

tests/unit/dfem/test_functional.cpp

@@ -0,0 +1,134 @@
+// 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"
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+
+namespace dfem_pa_kernels
+{
+
+template <int DIM>
+void dfem_functional(const char *filename, int p, const int r)
+{
+   CAPTURE(filename, DIM, p, r);
+
+   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();
+
+   Array<int> all_domain_attr;
+   if (pmesh.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(pmesh.attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace fes(&pmesh, &fec);
+   const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(), 2 * p + r);
+
+   ParGridFunction x(&fes), y(&fes), z(&fes);
+   Vector X(fes.GetTrueVSize()), Y(fes.GetTrueVSize()), Z(fes.GetTrueVSize());
+
+   X = 1.0;
+   x.SetFromTrueDofs(X);
+
+   ConstantCoefficient one(1.0);
+   ParBilinearForm blf(&fes);
+   blf.AddDomainIntegrator(new MassIntegrator(one, ir));
+   blf.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   blf.Assemble();
+   blf.Mult(x, y);
+   fes.GetProlongationMatrix()->MultTranspose(y, Y);
+
+   real_t sum_g, sum_l = Y.Sum();
+   MPI_Allreduce(&sum_l, &sum_g, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+
+   static constexpr int U = 0, Coords = 1;
+   const auto sol = std::vector{ FieldDescriptor{ U, &fes } };
+   DifferentiableOperator dop(sol, {{Coords, nodes->ParFESpace()}}, pmesh);
+
+   const auto functional_qf =
+      [] MFEM_HOST_DEVICE(
+#ifdef MFEM_USE_ENZYME
+         const real_t &u,
+#else
+         const dual<real_t, real_t> &u,
+#endif
+         const tensor<real_t, DIM, DIM> &J,
+         const real_t &w)
+   {
+      return tuple{u * w * det(J)};
+   };
+
+   auto derivatives = std::integer_sequence<size_t, U> {};
+   dop.AddDomainIntegrator(functional_qf,
+                           tuple{ Value<U>{}, Gradient<Coords>{}, Weight{} },
+                           tuple{ Sum<U>{} },
+                           *ir, all_domain_attr, derivatives);
+   dop.SetParameters({ nodes });
+
+   fes.GetRestrictionMatrix()->Mult(x, X);
+   Vector sum(1);
+   dop.Mult(X, sum);
+
+   REQUIRE(MFEM_Approx(0.0) == std::abs(sum_g - sum(0)));
+
+   auto dRdu = dop.GetDerivative(U, {&x}, {nodes});
+   Vector d(1);
+   dRdu->Mult(X, d);
+   REQUIRE(MFEM_Approx(0.0) == std::abs(d(0) - sum(0)));
+
+   MPI_Barrier(MPI_COMM_WORLD);
+}
+
+TEST_CASE("DFEM Functional", "[Parallel][DFEM][Functional]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+
+   const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
+   const auto r = !all_tests ? 1 : GENERATE(0, 1, 2, 3);
+
+   SECTION("2D p=" + std::to_string(p) + " r=" + std::to_string(r))
+   {
+      const auto filename =
+         GENERATE("../../data/star.mesh",
+                  "../../data/star-q3.mesh",
+                  "../../data/rt-2d-q3.mesh",
+                  "../../data/inline-quad.mesh",
+                  "../../data/periodic-square.mesh");
+      dfem_functional<2>(filename, p, r);
+   }
+
+   SECTION("3D p=" + std::to_string(p) + " r=" + std::to_string(r))
+   {
+      const auto filename =
+         GENERATE("../../data/fichera.mesh",
+                  "../../data/fichera-q3.mesh",
+                  "../../data/inline-hex.mesh",
+                  "../../data/toroid-hex.mesh",
+                  "../../data/periodic-cube.mesh");
+      dfem_functional<3>(filename, p, r);
+   }
+}
+
+} // namespace dfem_pa_kernels
+
+#endif // MFEM_USE_MPI