Some suggested refactors.

src/serac/infrastructure/input.hpp

@@ -16,10 +16,15 @@
 #include <variant>
 #include <optional>
 
-#include "mfem.hpp"
 #include "axom/inlet.hpp"
 #include "axom/sidre.hpp"
 
+namespace mfem {
+class Vector;
+class Coefficient;
+class VectorCoefficient;
+}  // namespace mfem
+
 /**
  * @brief The input related helper functions and objects
  *

src/serac/mesh/mesh_utils_base.hpp

@@ -15,10 +15,14 @@
 
 #include <memory>
 #include <variant>
-#include "mfem.hpp"
 
 #include "serac/infrastructure/input.hpp"
 
+namespace mfem {
+class Mesh;
+class ParMesh;
+}  // namespace mfem
+
 namespace serac {
 
 /**

src/serac/numerics/equation_solver.hpp

@@ -16,12 +16,16 @@
 #include <optional>
 #include <variant>
 
-#include "mfem.hpp"
-
 #include "serac/infrastructure/input.hpp"
 #include "serac/numerics/solver_config.hpp"
 #include "serac/numerics/petsc_solvers.hpp"
 
+namespace mfem {
+class Solver;
+class NewtonSolver;
+class vector;
+}  // namespace mfem
+
 namespace serac {
 
 /**

src/serac/numerics/functional/functional.hpp

@@ -444,7 +444,7 @@ class Functional<test(trials...), exec> {
    *  at most one of which may be of the type `differentiate_wrt_this(mfem::Vector)`
    */
   template <uint32_t wrt, typename... T>
-  typename operator_paren_return<wrt>::type operator()(DifferentiateWRT<wrt>, double t, const T&... args)
+  typename operator_paren_return<wrt>::type operator()(DifferentiateWRT<wrt>, double t, const T&... args) const
   {
     const mfem::Vector* input_T[] = {&static_cast<const mfem::Vector&>(args)...};
 
@@ -499,7 +499,7 @@ class Functional<test(trials...), exec> {
 
   /// @overload
   template <typename... T>
-  auto operator()(double t, const T&... args)
+  auto operator()(double t, const T&... args) const
   {
     constexpr int num_differentiated_arguments = (std::is_same_v<T, differentiate_wrt_this> + ...);
     static_assert(num_differentiated_arguments <= 1,
@@ -565,6 +565,13 @@ class Functional<test(trials...), exec> {
       return df_;
     }
 
+    /// @brief syntactic sugar:  df_dx.Mult(dx, df)  <=>  mfem::Vector df = df_dx(dx);
+    const mfem::Vector& operator()(const mfem::Vector& dx) const
+    {
+      form_.ActionOfGradient(dx, df_, which_argument);
+      return df_;
+    }
+
     /// @brief assemble element matrices and form an mfem::HypreParMatrix
     std::unique_ptr<mfem::HypreParMatrix> assemble()
     {

src/serac/numerics/solver_config.hpp

@@ -13,8 +13,6 @@
 #pragma once
 
 #include <variant>
-
-#include "mfem.hpp"
 #include "axom/fmt.hpp"
 
 namespace serac {

src/serac/physics/common.hpp

@@ -11,8 +11,6 @@
  */
 #pragma once
 
-#include "serac/numerics/odes.hpp"
-
 namespace serac {
 
 /**

src/serac/physics/heat_transfer_input.hpp

@@ -14,7 +14,7 @@
 #pragma once
 
 #include "serac/physics/materials/thermal_material_input.hpp"
-#include "serac/physics/common.hpp"
+#include "serac/numerics/odes.hpp"
 
 namespace serac {
 

src/serac/physics/solid_mechanics_input.hpp

@@ -15,6 +15,7 @@
 
 #include "serac/physics/materials/solid_material_input.hpp"
 #include "serac/physics/common.hpp"
+#include "serac/numerics/odes.hpp"
 
 namespace serac {