mod

2026-04-01 23:42:36 -07:00 · 2026-04-01 23:39:06 -07:00 · 2026-04-01 23:33:36 -07:00 · 2026-04-01 23:25:07 -07:00 · 2026-04-01 23:20:44 -07:00 · 2026-03-31 10:10:51 -07:00
34 changed files with 3900 additions and 1465 deletions
@@ -652,6 +652,8 @@ foreach(TPL IN LISTS MFEM_TPLS)
  endif()
 endforeach(TPL)

+# reverse to remove the first instance of entries in TPL_LIBRARIES
+# so later duplicates are kept (for dependency ordering)
 list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
 list(REVERSE TPL_LIBRARIES)
@@ -1015,5 +1017,7 @@ install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
 # Create 'config.mk' from 'config.mk.in' for the build and install locations and
 # define install rules for 'config.mk' and 'test.mk'
 #-------------------------------------------------------------------------------
-
+if (MFEM_USE_CUDA OR MFEM_USE_HIP)
+  option(MFEM_EXPORT_GPU_CONFIG "Export config.mk for GPU-enabled downstream packages" ON)
+endif()
 mfem_export_mk_files()
@@ -701,7 +701,6 @@ endfunction(mfem_find_library)
 # Extract compile and link options needed by the given target.
 #
 function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
-
  if (NOT TARGET ${Target})
    return()
  endif()
@@ -799,7 +798,12 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
          # message(STATUS "Lib = ${Lib}")
          # Filter-out generator expressions
          if (NOT ("${Lib}" MATCHES "^\\$"))
-            list(APPEND LinkOpts "${Lib}")
+            if(NOT ("${Lib}" STREQUAL "dl"))
+              list(APPEND LinkOpts "${Lib}")
+            else()
+              # for some reason libdl doesn't include the "-l"
+              list(APPEND LinkOpts "-ldl")
+            endif()
          endif()
        else()
          mfem_get_target_options(${Lib} COpts LOpts)
@@ -888,9 +892,18 @@ function(mfem_export_mk_files)
      set(${var} NO)
    endif()
  endforeach()
-  # TODO: Add support for MFEM_USE_CUDA=YES
-  set(MFEM_CXX ${CMAKE_CXX_COMPILER})
-  set(MFEM_HOST_CXX ${MFEM_CXX})
+  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
+    set(MFEM_CXX ${CMAKE_CUDA_COMPILER})
+    if(MFEM_CUDA_COMPILER_IS_NVCC)
+      set(MFEM_HOST_CXX ${CMAKE_CUDA_HOST_COMPILER})
+    else()
+      set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
+    endif()
+  else()
+    # mfem doesn't use enable_language(HIP)
+    set(MFEM_CXX ${CMAKE_CXX_COMPILER})
+    set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
+  endif()
  set(MFEM_CPPFLAGS "")
  get_target_property(cxx_std mfem CXX_STANDARD)
  # For now, we ignore the setting of the CXX_EXTENSIONS property. If this
@@ -900,6 +913,50 @@ function(mfem_export_mk_files)
  string(STRIP
         "${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
         MFEM_CXXFLAGS)
+  if(MFEM_EXPORT_GPU_CONFIG)
+    if (MFEM_USE_CUDA)
+      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} ${CMAKE_CUDA_FLAGS}")
+      if (MFEM_CUDA_COMPILER_IS_NVCC)
+        set(MFEM_CXXFLAGS "-x=cu ${MFEM_CXXFLAGS} -ccbin ${CMAKE_CXX_COMPILER} --forward-unknown-to-host-compiler")
+        # The following intentionally hides CUDA deprecation warnings
+        foreach(ENTRY IN LISTS CUDAToolkit_INCLUDE_DIRS)
+          set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -isystem ${ENTRY}")
+        endforeach()
+        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
+          # architecture flags not part of CMAKE_CUDA_FLAGS
+          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
+            set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -arch=${CMAKE_CUDA_ARCHITECTURES}")
+          else()
+            foreach (ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
+              set(MFEM_CXXFLAGS
+                "${MFEM_CXXFLAGS} -gencode arch=compute_${ENTRY},code=sm_${ENTRY}")
+            endforeach()
+          endif()
+        endif()
+      else()
+        set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xcuda --cuda-path=${CUDAToolkit_LIBRARY_ROOT}")
+        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
+          # architecture flags not part of CMAKE_CUDA_FLAGS
+          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
+            # TODO: not supported
+          else()
+            foreach(ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
+              set(MFEM_CXXFLAGS "-cuda-gpu-arch=sm_${ENTRY} ${MFEM_CXXFLAGS}")
+            endforeach()
+          endif()
+        endif()
+      endif()
+    elseif (MFEM_USE_HIP)
+      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xhip")
+      foreach(ENTRY IN LISTS CMAKE_HIP_ARCHITECTURES)
+        set(MFEM_CXXFLAGS "--offload-arch=${ENTRY} ${MFEM_CXXFLAGS}")
+      endforeach()
+    endif()
+  endif()
  set(MFEM_TPLFLAGS "")
  foreach(dir ${TPL_INCLUDE_DIRS})
    set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
@@ -930,6 +987,9 @@ function(mfem_export_mk_files)
    set(MFEM_SHARED NO)
    set(MFEM_STATIC YES)
  endif()
+  if (MFEM_USE_CUDA)
+    set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} -lcudart")
+  endif()
  set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
  set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
  # For the next 4 variables, these are the values for the build-tree version of
@@ -938,8 +998,15 @@ function(mfem_export_mk_files)
  set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
  set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
  set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
-  # TODO: CUDA/HIP support:
-  set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
+  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
+    if (MFEM_CUDA_COMPILER_IS_NVCC)
+      set(MFEM_XLINKER "-Xlinker=")
+    else()
+      set(MFEM_XLINKER "${CMAKE_CUDA_LINKER_WRAPPER_FLAG}")
+    endif()
+  else()
+    set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
+  endif()
  set(MFEM_MPIEXEC ${MPIEXEC})
  if (NOT MFEM_MPIEXEC)
    set(MFEM_MPIEXEC "mpirun")
@@ -987,16 +1054,21 @@ function(mfem_export_mk_files)
    # handle interfaces (e.g., SCOREC::apf)
    if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*" OR "${lib}" MATCHES "ParMoonolith::.*")
    elseif (TARGET "${lib}")
-      mfem_get_target_options(${lib} CompileOpts LinkOpts)
+      mfem_get_target_options(${lib} CompileOpts2 LinkOpts2)
+      # remove generator expressions
+      string(GENEX_STRIP "${CompileOpts2}" CompileOpts)
+      string(GENEX_STRIP "${LinkOpts2}" LinkOpts)
      # Removing duplicates may lead to issues:
      # list(REMOVE_DUPLICATES CompileOpts)
      # list(REMOVE_DUPLICATES LinkOpts)
-      string(REPLACE ";" " " COpts "${CompileOpts}")
-      string(REPLACE ";" " " LOpts "${LinkOpts}")
-      # message(STATUS "${lib}[COpts]: '${COpts}'")
-      # message(STATUS "${lib}[LOpts]: '${LOpts}'")
-      set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpts}")
-      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpts}")
+      # message(WARNING "${lib}[LinkOpts]: ${LinkOpts}")      
+      # message(WARNING "${lib}[CompileOpts]: ${CompileOpts}")
+      foreach(LOpt IN LISTS LinkOpts)
+        set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpt}")
+      endforeach()
+      foreach(COpt IN LISTS CompileOpts)
+        set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpt}")
+      endforeach()
      # message(FATAL_ERROR "***** interface lib found ... exiting *****")
      # handle static and shared libs
    elseif ("${suffix}" STREQUAL "${CMAKE_SHARED_LIBRARY_SUFFIX}")
@@ -1004,7 +1076,7 @@ function(mfem_export_mk_files)
      get_filename_component(fullLibName ${lib} NAME_WE)
      string(REGEX REPLACE "^lib" "" libname ${fullLibName})
      set(MFEM_EXT_LIBS
-          "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
+        "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
    else()
      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${lib}")
    endif()
@@ -1013,7 +1085,7 @@ function(mfem_export_mk_files)
  # Create the build-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config.mk")
+    "${PROJECT_BINARY_DIR}/config/config.mk" @ONLY)
  # Copy 'test.mk' from the source-tree to the build-tree
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/test.mk"
@@ -1031,7 +1103,7 @@ function(mfem_export_mk_files)
  # Create the install-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config-install.mk")
+    "${PROJECT_BINARY_DIR}/config/config-install.mk" @ONLY)

  # Install rules for 'config.mk' and 'test.mk'
  install(FILES ${PROJECT_SOURCE_DIR}/config/test.mk
@@ -82,6 +82,25 @@ public:
   /// underlying #fes
   int VectorDim() const;

+   /// Copy assignment. Only the data of the base class Vector is copied.
+   /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
+       have the same size.
+
+       @note Defining this method overwrites the implicitly defined copy
+       assignment operator. */
+   ComplexGridFunction &operator=(const ComplexGridFunction &rhs)
+   { return operator=((const Vector &)rhs); }
+
+   /// Copy the data from @a v.
+   /** The size of @a v must be equal to double of the size of the associated
+       FiniteElementSpace #fes. */
+   ComplexGridFunction &operator=(const Vector &v)
+   {
+      MFEM_ASSERT(fes && v.Size() == 2*fes->GetVSize(), "");
+      Vector::operator=(v);
+      return *this;
+   }
+
   /// Assign constant values to the ComplexGridFunction data.
   ComplexGridFunction &operator=(const std::complex<real_t> & value)
   { *gfr = value.real(); *gfi = value.imag(); return *this; }
@@ -90,8 +90,8 @@ void map_quadrature_data_to_fields_impl(
   }
   else
   {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        " for this field descriptor");
   }
 }

@@ -169,8 +169,9 @@ void map_quadrature_data_to_fields_tensor_impl_1d(
   }
   else
   {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor with sum factorization on tensor product elements");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        "for this field descriptor with sum factorization on"
+                        " tensor product elements");
   }
 }

@@ -306,8 +307,9 @@ void map_quadrature_data_to_fields_tensor_impl_2d(
   }
   else
   {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor with sum factorization on tensor product elements");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        " for this field descriptor with sum factorization on"
+                        " tensor product elements");
   }
 }

@@ -492,8 +494,9 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
   }
   else
   {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor with sum factorization on tensor product elements");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        " for this field descriptor with sum factorization on"
+                        " tensor product elements");
   }
 }

@@ -82,7 +82,7 @@ protected:
 };

 /// @brief Uniform parameter space
-class UniformParameterSpace final : public ParameterSpace
+class UniformParameterSpace : public ParameterSpace
 {
 public:
   /// @brief Constructor for a uniform parameter space
@@ -1044,9 +1044,50 @@ void VectorFiniteElement::SetDerivMembers()
   switch (map_type)
   {
      case H_DIV:
-         deriv_type = DIV;
-         deriv_range_type = SCALAR;
-         deriv_map_type = INTEGRAL;
+         switch (dim)
+         {
+            case 3: // div: 3D H_DIV -> 3D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            case 2: // div: 2D H_DIV -> 2D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_DIV_R2D:
+         switch (dim)
+         {
+            case 2: // div: 2D H_DIV_R2D -> 2D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            case 1: // div: 1D H_DIV_R2D -> 1D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_DIV_R1D:
+         switch (dim)
+         {
+            case 1: // div: 1D H_DIV_R1D -> 1D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
         break;
      case H_CURL:
         switch (dim)
@@ -1064,13 +1105,49 @@ void VectorFiniteElement::SetDerivMembers()
               break;
            case 1:
               deriv_type = NONE;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
+               deriv_range_type = UNKNOWN_RANGE_TYPE;
+               deriv_map_type = UNKNOWN_MAP_TYPE;
               break;
            default:
               MFEM_ABORT("Invalid dimension, Dim = " << dim);
         }
         break;
+      case H_CURL_R2D:
+         switch (dim)
+         {
+            case 2:
+               // curl: 2D H_CURL_R2D -> H_DIV_R2D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R2D;
+               break;
+            case 1:
+               // curl: 1D H_CURL_R2D -> H_DIV_R2D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R2D;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_CURL_R1D:
+         switch (dim)
+         {
+            case 1:
+               // curl: 1D H_CURL_R1D -> H_DIV_R1D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R1D;
+               break;
+            case 0:
+               deriv_type = NONE;
+               deriv_range_type = UNKNOWN_RANGE_TYPE;
+               deriv_map_type = UNKNOWN_MAP_TYPE;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
      default:
         MFEM_ABORT("Invalid MapType = " << map_type);
   }
@@ -295,10 +295,20 @@ public:
                          $ u(x) = (1/w) \hat u(\hat x) $ */
      H_DIV,     /**< For vector fields; preserves surface integrals of the
                          normal component $ u(x) = (J/w) \hat u(\hat x) $ */
-      H_CURL     /**< For vector fields; preserves line integrals of the
+      H_CURL,    /**< For vector fields; preserves line integrals of the
                          tangential component
                          $ u(x) = J^{-t} \hat u(\hat x) $ (square J),
                          $ u(x) = J(J^t J)^{-1} \hat u(\hat x) $ (general J) */
+      H_DIV_R2D, /**< For 3-component vector fields in 2D; equivalent to a
+                          direct sum of an H_DIV basis and an INTEGRAL basis */
+      H_CURL_R2D,/**< For 3-component vector fields in 2D; equivalent to a
+                          direct sum of an H_CURL basis and a VALUE basis */
+      H_DIV_R1D, /**< For 3-component vector fields in 1D; equivalent to a
+                          direct sum of a VALUE basis and a pair of INTEGRAL
+                          bases */
+      H_CURL_R1D /**< For 3-component vector fields in 1D; equivalent to a
+                          direct sum of an INTEGRAL basis and a pair of VALUE
+                          bases */
   };

   /** @brief Enumeration for DerivType: defines which derivative method
@@ -330,12 +340,28 @@ public:
   int GetDim() const { return dim; }

   /** @brief Returns the vector dimension for vector-valued finite elements,
-       which is also the dimension of the interpolation operation. */
+       which is also the dimension of the interpolation operation and the
+       width of the DenseMatrix argument in
+       CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape). */
   int GetRangeDim() const { return vdim; }

-   /// Returns the dimension of the curl for vector-valued finite elements.
+   /** @brief Returns the vector dimension, in physical space, for
+       vector-valued finite elements, which is also the width of the
+       DenseMatrix argument in
+       CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape). */
+   int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
+
+   /** Returns the dimension of the curl for vector-valued finite elements,
+       which is also the width of the DenseMatrix argument in
+       CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape). */
   int GetCurlDim() const { return cdim; }

+   /** Returns the dimension, in physical space, of the curl for vector-valued
+       finite elements, which is also the width of the DenseMatrix argument in
+       CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape).
+   */
+   int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
+
   /// Returns the Geometry::Type of the reference element.
   Geometry::Type GetGeomType() const { return geom_type; }

@@ -990,6 +1016,8 @@ protected:
 public:
   VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M,
                       int F = FunctionSpace::Pk);
+
+   int GetPhysRangeDim(int space_dim) const { return space_dim; }
 };

 /// @brief Class for computing 1D special polynomials and their associated basis
@@ -2531,7 +2531,7 @@ void ND_FuentesPyramidElement::calcCurlBasis(const int p,

 ND_R1D_PointElement::ND_R1D_PointElement(int p)
   : VectorFiniteElement(1, Geometry::POINT, 2, p,
-                         H_CURL, FunctionSpace::Pk)
+                         H_CURL_R1D, FunctionSpace::Pk)
 {
   // VectorFiniteElement::SetDerivMembers doesn't support 0D H_CURL elements
   // so we mimic a 1D element and then correct the dimension here.
@@ -2562,7 +2562,7 @@ ND_R1D_SegmentElement::ND_R1D_SegmentElement(const int p,
                                             const int cb_type,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 2, p,
-                         H_CURL, FunctionSpace::Pk),
+                         H_CURL_R1D, FunctionSpace::Pk),
     dof2tk(dof),
     cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
     obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -2839,7 +2839,7 @@ ND_R2D_SegmentElement::ND_R2D_SegmentElement(const int p,
                                             const int cb_type,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, 2 * p + 1, p,
-                         H_CURL, FunctionSpace::Pk),
+                         H_CURL_R2D, FunctionSpace::Pk),
     dof2tk(dof),
     cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
     obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -3023,7 +3023,7 @@ void ND_R2D_SegmentElement::Project(VectorCoefficient &vc,
 ND_R2D_FiniteElement::ND_R2D_FiniteElement(int p, Geometry::Type G, int Do,
                                           const real_t *tk_fe)
   : VectorFiniteElement(2, G, Do, p,
-                         H_CURL, FunctionSpace::Pk),
+                         H_CURL_R2D, FunctionSpace::Pk),
     tk(tk_fe),
     dof_map(dof),
     dof2tk(dof)
@@ -663,6 +663,9 @@ public:
                         const int cb_type = BasisType::GaussLobatto,
                         const int ob_type = BasisType::GaussLegendre);

+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { return 1; }
+
   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;

@@ -705,6 +708,9 @@ private:
                           DenseMatrix &I) const;

 public:
+   int GetPhysRangeDim(int space_dim) const { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 3; }
+
   using FiniteElement::CalcVShape;
   using FiniteElement::CalcPhysCurlShape;

@@ -2006,7 +2006,7 @@ RT_R1D_SegmentElement::RT_R1D_SegmentElement(const int p,
                                             const int cb_type,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 4, p + 1,
-                         H_DIV, FunctionSpace::Pk),
+                         H_DIV_R1D, FunctionSpace::Pk),
     dof2nk(dof),
     cbasis1d(poly1d.GetBasis(p + 1, VerifyClosed(cb_type))),
     obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
@@ -2281,7 +2281,7 @@ const real_t RT_R2D_SegmentElement::nk[2] = { 0.,1.};
 RT_R2D_SegmentElement::RT_R2D_SegmentElement(const int p,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, p + 1, p + 1,
-                         H_DIV, FunctionSpace::Pk),
+                         H_DIV_R2D, FunctionSpace::Pk),
     dof2nk(dof),
     obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
 {
@@ -2392,7 +2392,7 @@ void RT_R2D_SegmentElement::LocalInterpolation(const VectorFiniteElement &cfe,
 RT_R2D_FiniteElement::RT_R2D_FiniteElement(int p, Geometry::Type G, int Do,
                                           const real_t *nk_fe)
   : VectorFiniteElement(2, G, Do, p + 1,
-                         H_DIV, FunctionSpace::Pk),
+                         H_DIV_R2D, FunctionSpace::Pk),
     nk(nk_fe),
     dof_map(dof),
     dof2nk(dof)
@@ -510,6 +510,9 @@ public:
   RT_R2D_SegmentElement(const int p,
                         const int ob_type = BasisType::GaussLegendre);

+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { return 0; }
+
   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;

@@ -547,6 +550,9 @@ private:
                           DenseMatrix &I) const;

 public:
+   int GetPhysRangeDim(int space_dim) const { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 0; }
+
   using FiniteElement::CalcVShape;

   void CalcVShape(ElementTransformation &Trans,
@@ -321,12 +321,17 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
   const int vd = vdim;
   const bool t = byvdim;
   const int threshold = ndofs;
-   const int nsdofs = pfes.GetFaceNbrVSize();
+   const int nsdofs = pfes.GetFaceNbrVSize() / vd;
   auto d_indices1 = scatter_indices1.Read();
   auto d_indices2 = scatter_indices2.Read();
   auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
-   auto d_x_shared = Reshape(face_nbr_data.Read(),
-                             t?vd:nsdofs, t?nsdofs:vd);
+   const int ne_shared = nsdofs / elem_dofs;
+   const int nedof = elem_dofs;
+   // Note: the shape of face_nbr_data, as determined by
+   // ParFiniteElementSpace::ExchangeFaceNbrData, is (elem_dofs, vdim,
+   // ne_shared), independent of the ordering (byNODES or byVDIM) of the finite
+   // element space.
+   auto d_x_shared = Reshape(face_nbr_data.Read(), elem_dofs, vd, ne_shared);
   auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
   mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
   {
@@ -346,8 +351,9 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
         }
         else if (idx2>=threshold) // shared boundary
         {
-            d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
-                                              t?(idx2-threshold):c);
+            const int e_shared = (idx2 - threshold) / nedof;
+            const int i_shared = (idx2 - threshold) % nedof;
+            d_y(dof, c, 1, face) = d_x_shared(i_shared,c,e_shared);
         }
         else // true boundary
         {
@@ -1398,20 +1398,17 @@ void L2FaceRestriction::PermuteAndSetSharedFaceDofsScatterIndices2(
   const int dim = fes.GetMesh()->Dimension();
   const int dof1d = fes.GetTypicalFE()->GetOrder()+1;
   fes.GetTypicalFE()->GetFaceMap(face_id2, face_map);
-   Array<int> face_nbr_dofs;
-   const ParFiniteElementSpace &pfes =
-      static_cast<const ParFiniteElementSpace&>(this->fes);
-   pfes.GetFaceNbrElementVDofs(elem_index, face_nbr_dofs);

   for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
   {
      const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
                                               orientation, dof1d, face_dof_elem1);
      const int volume_dof_elem2 = face_map[face_dof_elem2];
-      const int global_dof_elem2 = face_nbr_dofs[volume_dof_elem2];
+      // Encode the volume DOF index and element index
+      const int global_dof_elem2 = elem_index*elem_dofs + volume_dof_elem2;
      const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
      // Trick to differentiate dof location inter/shared
-      scatter_indices2[restriction_dof_elem2] = ndofs+global_dof_elem2;
+      scatter_indices2[restriction_dof_elem2] = ndofs + global_dof_elem2;
   }
 #endif
 }
@@ -22,14 +22,6 @@
 //#define _WIN32
 //#define _aligned_malloc(s,a) malloc(s)

-#ifdef NVTX_DEBUG_HPP
-#undef NVTX_COLOR
-#define NVTX_COLOR ::nvtx::kGold
-#include NVTX_DEBUG_HPP
-#else
-#define dbg(...)
-#endif
-
 #ifndef _WIN32
 #include <unistd.h>
 #include <signal.h>
@@ -768,12 +760,7 @@ private:
      switch (mt)
      {
         case MT::HOST_DEBUG:
-            if (GetEnv("MFEM_MMU_STD"))
-            {
-               dbg("Using STD memory space for debug device!");
-               return new StdHostMemorySpace();
-            }
-            dbg("Using MMU memory space for debug device!");
+            if (GetEnv("MFEM_MMU_STD")) { return new StdHostMemorySpace(); }
            return new MmuHostMemorySpace();
 #ifdef MFEM_USE_UMPIRE
         case MT::HOST_UMPIRE:
@@ -284,7 +284,7 @@ void RK4Solver::Init(TimeDependentOperator &f_)
   int n = f->Width();
   y.SetSize(n, mem_type);
   k.SetSize(n, mem_type);
-   z.SetSize(n, mem_type);
+   z.SetSize(n, mem_type);   
 }

 void RK4Solver::Step(Vector &x, real_t &t, real_t &dt)
@@ -316,6 +316,67 @@ void RK4Solver::Step(Vector &x, real_t &t, real_t &dt)
   t += dt;
 }

+void RK4Solver::SetSolution(const Vector& x, real_t ts)
+{
+   hist.clear(); //we need only the last time step at time ts
+   StepRecord rec{ts,x};
+   hist.push_back(rec);
+}
+
+void RK4Solver::AdjointStep(Vector &lambda, real_t &t, real_t &dt)
+{
+   MFEM_VERIFY(!hist.empty(), "No forward history.");
+   MFEM_VERIFY(adj_mode == AdjointMode::Discrete, "Wrong adjoint mode.");
+
+   //Restore the solution history, i.e., perform a forward step.
+   auto rec = hist.back();
+   hist.pop_back();
+   const real_t tn = rec.t;
+   const real_t h  = dt;
+   const Vector &xn = rec.x;
+   // Recompute forward stage states and stage RHS:
+   f->SetTime(tn);
+   f->Mult(xn, k1); // k1
+   add(xn, h/2.0, k1, y2); // y2
+
+   f->SetTime(tn+h/2.0);
+   f->Mult(y2, k); // k2
+   add(xn, h/2.0, k, y3); // y3
+
+   f->Mult(y3, k); // k3
+   add(xn, h, k, y4); // y4
+
+
+   //set the initial values for a1,a2,a3,a4
+   a1.Set(h/6.0,lambda);
+   a2.Set(h/3.0,lambda);
+   a3.Set(h/3.0,lambda);
+   a4.Set(h/6.0,lambda);
+   
+   //stage 4
+   f->SetTime(tn+h);
+   f->JacobianMultTranspose(y4, a4, k);
+   lambda.Add(1.0,k);
+   a3.Add(h,k);
+
+   //stage 3
+   f->SetTime(tn+h/2.0);
+   f->JacobianMultTranspose(y3, a3, k);
+   lambda.Add(1.0,k);
+   a2.Add(h/2.0,k);
+
+   //stage 2
+   f->SetTime(tn+h/2.0);
+   f->JacobianMultTranspose(y2, a2, k);
+   lambda.Add(1.0,k);
+   a1.Add(h/2.0,k);
+
+   //stage 1
+   f->SetTime(tn);
+   f->JacobianMultTranspose(xn, a1, k);
+   lambda.Add(1.0,k);
+}
+
 ExplicitRKSolver::ExplicitRKSolver(int s_, const real_t *a_, const real_t *b_,
                                   const real_t *c_)
 {
@@ -190,6 +190,78 @@ public:
      while (t < tf) { Step(x, t, dt); }
   }

+   /// @brief The enum defines the supported adjoint modes.
+   enum class AdjointMode { None, Continuous, Discrete };
+
+   AdjointMode adj_mode= AdjointMode::None;
+
+   /// Sets the solution at time ts. The solution will be used by the adjoint step
+   /// to recompute the internal stages/steps between ts and the time t used in the 
+   /// AdjointStep operator.
+   /// SetSolution(x,ts); and then calling AdjointStep(lambda, t, dt), will
+   /// first restore the solution/parametrization between  ts and t (if required by 
+   /// the time integration scheme) and then execute the adjoint step. 
+   virtual void SetSolution(const Vector& x, real_t ts){} 
+
+   
+   /// Return primal state x(t) for t in the recorded time interval.
+   /// Default false => solver does not provide primal interpolation access.
+   virtual bool GetForwardSolution(real_t t, Vector &x) const 
+   { 
+      return false; 
+   }
+
+   /// Perform one backward step for the adjoint variable lambda.
+   /// dt > 0 requests stepping from t to t - dt.
+   /// The solver may adjust dt (adaptive / grid snapping) and returns dt used.
+   /** @param[in,out] lambda   Approximate solution.
+       @param[in,out] t   Time associated with the approximate solution @a x.
+       @param[in,out] dt  Time step size.
+   */
+   virtual void AdjointStep(Vector &lambda, real_t &t, real_t &dt)
+   {
+      MFEM_ABORT("Adjoint stepping is not supported by this ODESolver.");
+   }
+
+   /// Perform adjoint integration from time @a t [in] to time @a tf [in].
+   /** @param[in,out] lambda   Approximate adjoint solution.
+       @param[in,out] t        Time associated with the approximate solution @a lambda.
+       @param[in,out] dt  Time step size.
+       @param[in]     t0  Requested final time. (t0<t)
+   */
+   virtual void AdjointRun(Vector &lambda, real_t &t, real_t &dt, real_t t0)
+   {
+      MFEM_VERIFY(t0 <= t, "AdjointRun expects t0 <= t.");
+      while (t > t0) { 
+         real_t d = std::min(dt, t - t0); 
+         AdjointStep(lambda, t, d); 
+         dt = d; 
+      }
+   }
+
+   /// Enable adjoint functionality and tell the solver 
+   /// what to record during forward stepping. 
+   /// Default: abort (adjoint mode is not supported).
+   virtual void EnableAdjoint(AdjointMode mod)
+   {
+      if(AdjointMode::None!=mod)
+      {
+         MFEM_ABORT("Adjoint is not supported by the ODESolver.");
+      }
+   }
+
+   /// Return true of the ODESolver supports the AdjointMode
+   /// and false if not.
+   virtual bool SupportsAdjoint(AdjointMode mod) const
+   { 
+      if(AdjointMode::None==mod)
+      {
+         return true;
+      }
+      //else
+      return false;
+   }
+
   /// Returns how many State vectors the ODE requires
   virtual int GetStateSize() { return 0; };

@@ -309,10 +381,73 @@ class RK4Solver : public ODESolver
 private:
   Vector y, k, z;

+   //storage for the adjoint solver
+   /// @brief stores the state x at time t
+   struct StepRecord{
+      real_t t;
+      Vector x;
+   };
+
+   /// @brief stores several steps
+   std::vector<StepRecord> hist;
+
+   Vector k1,y2,y3,y4; //RK4 stages
+   Vector a1,a2,a3,a4; //adjoint vectors
+
 public:
   void Init(TimeDependentOperator &f_) override;

   void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   /// Return true of the ODESolver supports the AdjointMode and false
+   /// if not.
+   virtual bool SupportsAdjoint(AdjointMode mod) const override
+   { 
+      if(AdjointMode::Discrete==mod)
+      {
+         return true;
+      }
+      //else
+      return false;
+   }
+
+   /// Enable adjoint functionality and tell the solver 
+   /// what to record during forward stepping. 
+   /// Default: abort (adjoint mode is not supported).
+   virtual void EnableAdjoint(AdjointMode mod) override
+   {
+      if(AdjointMode::Continuous==mod)
+      {
+         MFEM_ABORT("Continuous adjoint not supported by the RK4Solver.");
+      }
+      adj_mode=mod;
+
+      // allocate the storage memeory
+      if(AdjointMode::Discrete==mod){
+         int n = f->Width();
+         k1.SetSize(n, mem_type);
+         y2.SetSize(n, mem_type);
+         y3.SetSize(n, mem_type);
+         y4.SetSize(n, mem_type);
+
+         a1.SetSize(n, mem_type);
+         a2.SetSize(n, mem_type);
+         a3.SetSize(n, mem_type);
+         a4.SetSize(n, mem_type);
+      }
+   }
+
+
+   void AdjointStep(Vector &lam, real_t &t, real_t &dt) override;
+
+   //Two posible ways to restore the histroy
+   //1) Use Step(Vector &x, real_t &t, real_t &dt) in discrete adjoint mode and record the internal states
+   //2) Use additional Method SetState(const Vector &x, real_t ts);
+   //   Calling AdjointStep(lam,t,dt) should restore the history between ts and t, and then execute the 
+   //   adjoint step in reverse to time t-dt with the recorded discrete states between ts and t. 
+   // Here we will use SetSolution
+   virtual void SetSolution(const Vector& x, real_t ts) override;
+
 };


@@ -640,6 +640,15 @@ public:
   virtual int SUNMassMult(const Vector &x, Vector &v);

   virtual ~TimeDependentOperator() { }
+
+   /// y = (df/dx(x,t))^T * w
+   /// Default uses GetGradient(x).MultTranspose(w,y).
+   virtual void JacobianMultTranspose(const Vector &x,
+                                      const Vector &w,
+                                      Vector &y) const
+   {
+      GetGradient(x).MultTranspose(w, y);
+   }
 };


@@ -38,6 +38,13 @@
 #if PETSC_VERSION_LT(3,19,0)
 #define PETSC_SUCCESS 0
 #endif
+#if PETSC_VERSION_LT(3,23,0)
+#define PetscContainerSetCtxDestroy(A,B) PetscContainerSetUserDestroy(A,B)
+typedef PetscErrorCode (PetscCtxDestroyFn)(void**);
+#endif
+#if PETSC_VERSION_LT(3,24,0)
+typedef PetscErrorCode KSPMonitorFn(KSP,PetscInt,PetscReal,void*);
+#endif

 #include <fstream>
 #include <iomanip>
@@ -77,13 +84,17 @@ static PetscErrorCode __mfem_mat_shell_apply_transpose(Mat,Vec,Vec);
 static PetscErrorCode __mfem_mat_shell_destroy(Mat);
 static PetscErrorCode __mfem_mat_shell_copy(Mat,Mat,MatStructure);
 #if PETSC_VERSION_LT(3,23,0)
-static PetscErrorCode __mfem_array_container_destroy(void*);
-static PetscErrorCode __mfem_matarray_container_destroy(void *);
-#else
-static PetscErrorCode __mfem_array_container_destroy(void**);
-static PetscErrorCode __mfem_matarray_container_destroy(void**);
+typedef void *PetscCtxRt;
+#elif PETSC_VERSION_LT(3,25,0)
+typedef void **PetscCtxRt;
 #endif
+static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt);
+static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt);
+#if PETSC_VERSION_LT(3,23,0)
 static PetscErrorCode __mfem_monitor_ctx_destroy(void**);
+#else
+static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt);
+#endif

 // auxiliary functions
 static PetscErrorCode Convert_Array_IS(MPI_Comm,bool,const mfem::Array<int>*,
@@ -1317,11 +1328,7 @@ BlockDiagonalConstructor(MPI_Comm comm,

         ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
         ierr = PetscContainerSetPointer(c,ptrs[i]); CCHKERRQ(comm,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-         ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
         ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
         CCHKERRQ(comm,ierr);
         ierr = PetscObjectCompose((PetscObject)A,names[i],(PetscObject)c);
         CCHKERRQ(comm,ierr);
@@ -1648,11 +1655,7 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
         PetscContainer c;
         ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
         ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-         ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
-#else
         ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
-#endif
         PCHKERRQ(c,ierr);
         ierr = PetscObjectCompose((PetscObject)(*A),"_MatIS_PtAP_l2l",(PetscObject)c);
         PCHKERRQ((*A),ierr);
@@ -1748,11 +1751,7 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,

            ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); PCHKERRQ(B,ierr);
            ierr = PetscContainerSetPointer(c,ptrs[i]); PCHKERRQ(B,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-            ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
            ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
            PCHKERRQ(B,ierr);
            ierr = PetscObjectCompose((PetscObject)(B),names[i],(PetscObject)c);
            PCHKERRQ(B,ierr);
@@ -2198,11 +2197,7 @@ PetscParMatrix * RAP(PetscParMatrix *Rt, PetscParMatrix *A, PetscParMatrix *P)
         ierr = PetscContainerCreate(PetscObjectComm((PetscObject)B),&c);
         PCHKERRQ(B,ierr);
         ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-         ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
-#else
         ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
-#endif
         PCHKERRQ(c,ierr);
         ierr = PetscObjectCompose((PetscObject)B,"_MatIS_PtAP_l2l",(PetscObject)c);
         PCHKERRQ(B,ierr);
@@ -2485,7 +2480,6 @@ void PetscSolver::SetMaxIter(int max_iter)

 void PetscSolver::SetPrintLevel(int plev)
 {
-   typedef PetscErrorCode (*myPetscFunc)(void**);
   PetscViewerAndFormat *vf = NULL;
   PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm(obj));

@@ -2498,7 +2492,6 @@ void PetscSolver::SetPrintLevel(int plev)
   {
      // there are many other options, see the function KSPSetFromOptions() in
      // src/ksp/ksp/interface/itcl.c
-      typedef PetscErrorCode (*myMonitor)(KSP,PetscInt,PetscReal,void*);
      KSP ksp = (KSP)obj;
      if (plev >= 0)
      {
@@ -2507,29 +2500,29 @@ void PetscSolver::SetPrintLevel(int plev)
      if (plev == 1)
      {
 #if PETSC_VERSION_LT(3,15,0)
-         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorDefault,vf,
+         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorDefault,vf,
 #else
-         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorResidual,vf,
+         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorResidual,vf,
 #endif
-                              (myPetscFunc)PetscViewerAndFormatDestroy);
+                              (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
         PCHKERRQ(ksp,ierr);
      }
      else if (plev > 1)
      {
         ierr = KSPSetComputeSingularValues(ksp,PETSC_TRUE); PCHKERRQ(ksp,ierr);
-         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorSingularValue,vf,
-                              (myPetscFunc)PetscViewerAndFormatDestroy);
+         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorSingularValue,vf,
+                              (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
         PCHKERRQ(ksp,ierr);
         if (plev > 2)
         {
            ierr = PetscViewerAndFormatCreate(viewer,PETSC_VIEWER_DEFAULT,&vf);
            PCHKERRQ(viewer,ierr);
 #if PETSC_VERSION_LT(3,15,0)
-            ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidualNorm,vf,
+            ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidualNorm,vf,
 #else
-            ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidual,vf,
+            ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidual,vf,
 #endif
-                                 (myPetscFunc)PetscViewerAndFormatDestroy);
+                                 (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
            PCHKERRQ(ksp,ierr);
         }
      }
@@ -2545,7 +2538,7 @@ void PetscSolver::SetPrintLevel(int plev)
      if (plev > 0)
      {
         ierr = SNESMonitorSet(snes,(myMonitor)SNESMonitorDefault,vf,
-                               (myPetscFunc)PetscViewerAndFormatDestroy);
+                               (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
         PCHKERRQ(snes,ierr);
      }
   }
@@ -5329,21 +5322,27 @@ static PetscErrorCode __mfem_pc_shell_destroy(PC pc)
   PetscFunctionReturn(PETSC_SUCCESS);
 }

+static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt ptr)
+{
+   PetscErrorCode ierr;
+
+   PetscFunctionBeginUser;
 #if PETSC_VERSION_LT(3,23,0)
-
-static PetscErrorCode __mfem_array_container_destroy(void *ptr)
-{
-   PetscErrorCode ierr;
-
-   PetscFunctionBeginUser;
   ierr = PetscFree(ptr); CHKERRQ(ierr);
+#else
+   ierr = PetscFree(*(void**)ptr); CHKERRQ(ierr);
+#endif
   PetscFunctionReturn(PETSC_SUCCESS);
 }

-static PetscErrorCode __mfem_matarray_container_destroy(void *ptr)
+static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt ptr)
 {
+#if PETSC_VERSION_LT(3,23,0)
   mfem::Array<Mat> *a = (mfem::Array<Mat>*)ptr;
-   PetscErrorCode   ierr;
+#else
+   mfem::Array<Mat> *a = *(mfem::Array<Mat>**)ptr;
+#endif
+   PetscErrorCode ierr;

   PetscFunctionBeginUser;
   for (int i=0; i<a->Size(); i++)
@@ -5356,41 +5355,16 @@ static PetscErrorCode __mfem_matarray_container_destroy(void *ptr)
   PetscFunctionReturn(PETSC_SUCCESS);
 }

+#if PETSC_VERSION_LT(3,23,0)
+static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
 #else
-
-static PetscErrorCode __mfem_array_container_destroy(void **ptr)
+static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt ctx)
+#endif
 {
   PetscErrorCode ierr;

   PetscFunctionBeginUser;
-   ierr = PetscFree(*ptr); CHKERRQ(ierr);
-   PetscFunctionReturn(PETSC_SUCCESS);
-}
-
-static PetscErrorCode __mfem_matarray_container_destroy(void **ptr)
-{
-   mfem::Array<Mat> *a = (mfem::Array<Mat>*)*ptr;
-   PetscErrorCode   ierr;
-
-   PetscFunctionBeginUser;
-   for (int i=0; i<a->Size(); i++)
-   {
-      Mat M = (*a)[i];
-      MPI_Comm comm = PetscObjectComm((PetscObject)M);
-      ierr = MatDestroy(&M); CCHKERRQ(comm,ierr);
-   }
-   delete a;
-   PetscFunctionReturn(PETSC_SUCCESS);
-}
-
-#endif
-
-static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
-{
-   PetscErrorCode  ierr;
-
-   PetscFunctionBeginUser;
-   ierr = PetscFree(*ctx); CHKERRQ(ierr);
+   ierr = PetscFree(*(void**)ctx); CHKERRQ(ierr);
   PetscFunctionReturn(PETSC_SUCCESS);
 }

@@ -5635,11 +5609,7 @@ static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix* hA,Mat* pA)

      ierr = PetscContainerCreate(comm,&c); CHKERRQ(ierr);
      ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
-#if PETSC_VERSION_LT(3,23,0)
-      ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
      ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
      CHKERRQ(ierr);
      ierr = PetscObjectCompose((PetscObject)(*pA),names[i],(PetscObject)c);
      CHKERRQ(ierr);
@@ -5733,11 +5703,7 @@ static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix* hA,Mat* pA)

      ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); CHKERRQ(ierr);
      ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
-#if PETSC_VERSION_LT(3,23,0)
-      ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
      ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
      CHKERRQ(ierr);
      ierr = PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
      CHKERRQ(ierr);
@@ -126,11 +126,11 @@ EXAMPLE_TEST_DIRS := examples
 MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
 toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
 hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
- fluids/navier fluids/schrodinger-flow
+ fluids/navier fluids/schrodinger-flow plasma
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys shifted dpg diag-smoothers fluids/navier)
+ toys shifted dpg diag-smoothers fluids/navier plasma)

 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)

@@ -37,6 +37,7 @@ add_subdirectory(multidomain)
 add_subdirectory(nurbs)
 add_subdirectory(parelag)
 add_subdirectory(performance)
+add_subdirectory(plasma)
 add_subdirectory(shifted)
 add_subdirectory(solvers)
 add_subdirectory(spde)
@@ -26,5 +26,10 @@ if(MFEM_USE_MPI AND MFEM_USE_DOUBLE)
 					 MAIN mtop_test_iso_elasticity.cpp
 					 ${PAR_MTOP_COMMON_FILES}
 					 LIBRARIES mfem)
+	add_mfem_miniapp(tst_rk4_adj
+				 	MAIN tst_rk4_adj.cpp
+					${PAR_MTOP_COMMON_FILES}
+                                         LIBRARIES mfem)
+

 endif (MFEM_USE_MPI AND MFEM_USE_DOUBLE)
@@ -31,6 +31,11 @@ add_mfem_miniapp(mtop_test_adjoint_le
                      ${EX_COMMON_FILES}
                      LIBRARIES mfem)

+add_mfem_miniapp(mtop_test_adjoint_op
+                      MAIN test_adjoint_op.cpp
+                      ${EX_COMMON_FILES}
+                      LIBRARIES mfem)
+

 endif (MFEM_USE_MPI)

@@ -1,329 +1,289 @@
-// 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.
 #ifndef LINEAR_ELASTICITY_HPP
 #define LINEAR_ELASTICITY_HPP
-
 #include "mfem.hpp"

-///////////////////////////////////////////////////////////////////////////////
-class LinearElasticityTimeDependentOperator: public mfem::TimeDependentOperator
+
+class LinearElasticityTimeDependentOperator : public mfem::TimeDependentOperator
 {
 public:

-   LinearElasticityTimeDependentOperator(mfem::ParMesh &mesh_, int vorder =1);
+LinearElasticityTimeDependentOperator(mfem::ParMesh &mesh_, int vorder =1);

-   virtual ~LinearElasticityTimeDependentOperator() { }
+virtual ~LinearElasticityTimeDependentOperator() override
+{
+}

-   void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+void MultInvMass(const mfem::Vector &x, mfem::Vector &y) const;

-   virtual void AdjointMult(const mfem::Vector &x, mfem::Vector &y) const;
+virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override;

-   void ImplicitSolve(const mfem::real_t dt,
-                      const mfem::Vector &x,
-                      mfem::Vector &k) override;
+void JacobianMultTranspose(const mfem::Vector &x,
+                                      const mfem::Vector &w,
+                                      mfem::Vector &y) const override;

-   // sets the objective integraand which will be integrated with the state of the system
-   void SetObjective(std::shared_ptr<mfem::Operator> op_);
+virtual void ImplicitSolve(const mfem::real_t dt,
+                             const mfem::Vector &x,
+                             mfem::Vector &k) override;

-   // Assemble the explicit operators
-   // must be called after setting all material coefficients
-   // and before time stepping
-   void AssembleExplicit();

-   void SetDensity(mfem::Coefficient &rho)
-   {
-      density = std::make_shared<mfem::CoefficientVector>(*qs,
-                                                          mfem::CoefficientStorage::FULL);
-      cdensity = &rho;
-      density->Project(rho);
-   }
+// Assemble the explicit operators
+// must be called after setting all material coefficients
+// and before time stepping
+void AssembleExplicit();

-   // l1, m1 are the Lame parameters for material 1
-   // l2, m2 are the Lame parameters for material 2
-   void SetElasticityCoefficients(mfem::Coefficient& l1_,
-                                  mfem::Coefficient& m1_,
-                                  mfem::Coefficient& l2_,
-                                  mfem::Coefficient& m2_)
-   {
-      l1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      l2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      m1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      m2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
+void SetDensity(mfem::Coefficient &rho)
+{
+    density = std::make_shared<mfem::CoefficientVector>(*qs, mfem::CoefficientStorage::FULL); 
+    cdensity = &rho;
+    density->Project(rho);
+}

-      l1->Project(l1_);
-      m1->Project(m1_);
-      l2->Project(l2_);
-      m2->Project(m2_);
-   }

-   // rho1 and rho2 are the density coefficients for material 1 and 2
-   void SetDensityMaterialCoefficients(mfem::Coefficient& rho1_,
-                                       mfem::Coefficient& rho2_)
-   {
-      dens1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      dens2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
+// l1, m1 are the Lame parameters for material 1
+// l2, m2 are the Lame parameters for material 2
+void SetElasticityCoefficients(mfem::Coefficient& l1_,
+                               mfem::Coefficient& m1_,
+                               mfem::Coefficient& l2_,
+                               mfem::Coefficient& m2_)
+{
+    l1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 
+    l2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 
+    m1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 
+    m2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));

-      cdens1 = &rho1_;
-      cdens2 = &rho2_;
+    l1->Project(l1_);
+    m1->Project(m1_);
+    l2->Project(l2_);
+    m2->Project(m2_);
+}

-      dens1->Project(rho1_);
-      dens2->Project(rho2_);
-   }

-   // c1 and c2 are velocity proportional damping coefficients for material 1 and 2
-   // both of them are considered to be time dependent
-   void SetDampingMaterialCoefficients(mfem::Coefficient& c1_,
-                                       mfem::Coefficient& c2_)
-   {
-      cm1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      cm2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
+// rho1 and rho2 are the density coefficients for material 1 and 2
+void SetDensityMaterialCoefficients(mfem::Coefficient& rho1_,
+                                    mfem::Coefficient& rho2_)
+{
+    dens1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 
+    dens2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 

-      cm1->Project(c1_);
-      cm2->Project(c2_);
-   }
+    cdens1 = &rho1_;
+    cdens2 = &rho2_;

-   // dl1, dm1 are the strain velocity proportional damping coefficients
-   // for material 1 and dl2, dm2 for material 2
-   // all of them are considered to be time dependent
-   void SetDampingMaterialCoefficients(mfem::Coefficient& dl1_,
-                                       mfem::Coefficient& dm1_,
-                                       mfem::Coefficient& dl2_,
-                                       mfem::Coefficient& dm2_)
-   {
-      dl1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      dl2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      dm1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
-      dm2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));
+    dens1->Project(rho1_);
+    dens2->Project(rho2_);
+}

-      dl1->Project(dl1_);
-      dm1->Project(dm1_);
-      dl2->Project(dl2_);
-      dm2->Project(dm2_);
-   }
+// c1 and c2 are velocity proportional damping coefficients for material 1 and 2
+// both of them are considered to be time dependent
+void SetDampingMaterialCoefficients(mfem::Coefficient& c1_,
+                                    mfem::Coefficient& c2_)
+{
+    cm1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));     
+    cm2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));     

-   mfem::ParGridFunction& GetDisplacement() { return displ; }
-   mfem::ParGridFunction& GetVelocity() { return veloc; }
-   mfem::Array<int>& GetTrueBlockOffsets() { return block_true_offsets;}
+    cm1->Project(c1_);
+    cm2->Project(c2_);
+}

-   mfem::Vector& GetState() {return sol;}
+// dl1, dm1 are the strain velocity proportional damping coefficients
+// for material 1 and dl2, dm2 for material 2
+// all of them are considered to be time dependent 
+void SetDampingMaterialCoefficients(mfem::Coefficient& dl1_,
+                                    mfem::Coefficient& dm1_,
+                                    mfem::Coefficient& dl2_,
+                                    mfem::Coefficient& dm2_)
+{
+    dl1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 
+    dl2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 
+    dm1.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL)); 
+    dm2.reset(new mfem::CoefficientVector(*qs, mfem::CoefficientStorage::FULL));    

-   const mfem::ParFiniteElementSpace* GetFESpace() { return fespace.get();}
+    dl1->Project(dl1_);
+    dm1->Project(dm1_);
+    dl2->Project(dl2_);
+    dm2->Project(dm2_);
+}             

-   void SetZeroBdr(int bdr_attr)
-   {
-      zero_bdrs.insert(bdr_attr);
-   }
+mfem::ParGridFunction& GetDisplacement() { return displ; }
+mfem::ParGridFunction& GetVelocity() { return veloc; }
+mfem::Array<int>& GetTrueBlockOffsets(){ return block_true_offsets;}

-   void SetBdrLoad(int attr)
-   {
-      bdr_loads_markers.insert(attr);
-   }
+mfem::Vector& GetState(){return sol;}

-   void SetVolForce(mfem::real_t period, mfem::real_t amplitude, mfem::real_t rad,
-                    mfem::real_t xc=0.0, mfem::real_t yc=0.0, mfem::real_t zc=0.0,
-                    mfem::real_t L=5.0, mfem::real_t t0=0.0, mfem::real_t n=2.0)
-   {
-      // copy data to the host
-      mfem::real_t* pvol_force_mem = vol_force_mem.HostReadWrite();
+mfem::ParFiniteElementSpace* GetFESpace(){ return fespace.get();}

-      pvol_force_mem[1] = period;
-      pvol_force_mem[2] = amplitude;
-      pvol_force_mem[3] = rad;
-      pvol_force_mem[4] = xc;
-      pvol_force_mem[5] = yc;
-      pvol_force_mem[6] = zc;
+void SetZeroBdr(int bdr_attr)
+{
+    zero_bdrs.insert(bdr_attr);
+}

-      pvol_force_mem[7] = L;
-      pvol_force_mem[8] = t0;
-      pvol_force_mem[9] = n;
+void SetBdrLoad(int attr)
+{
+    bdr_loads_markers.insert(attr);
+}

-      // copy data to the device
-      vol_force_mem.Read();
-   }
+void SetVolForce(mfem::real_t period, mfem::real_t amplitude, mfem::real_t rad,
+                     mfem::real_t xc=0.0, mfem::real_t yc=0.0, mfem::real_t zc=0.0,
+                     mfem::real_t L=5.0, mfem::real_t t0=0.0, mfem::real_t n=2.0)
+{
+    // copy data to the host 
+    mfem::real_t* pvol_force_mem = vol_force_mem.HostReadWrite();

-   void AddState(mfem::real_t t, mfem::Vector& state)
-   {
-      int cind=adjoint_data.ind % adjoint_data.max_states;
-      adjoint_data.states[cind]=state;
-      adjoint_data.times[cind]=t;
-      adjoint_data.ind++;
-   }
+    pvol_force_mem[1] = period;
+    pvol_force_mem[2] = amplitude;
+    pvol_force_mem[3] = rad;
+    pvol_force_mem[4] = xc;
+    pvol_force_mem[5] = yc;
+    pvol_force_mem[6] = zc;
+
+    pvol_force_mem[7] = L;
+    pvol_force_mem[8] = t0;
+    pvol_force_mem[9] = n;
+
+    // copy data to the device
+    vol_force_mem.Read();
+}

 private:
-   mfem::ParMesh &mesh;
-   int order;
+mfem::ParMesh &mesh;
+int order;

-   std::unique_ptr<mfem::FiniteElementCollection>  fec;
-   std::unique_ptr<mfem::ParFiniteElementSpace>  fespace;
-   int dim;
-   int space_dim;
+std::unique_ptr<mfem::FiniteElementCollection>  fec;
+std::unique_ptr<mfem::ParFiniteElementSpace>  fespace;
+int dim;
+int space_dim;

-   int myrank;
+int myrank;

-   mfem::ParGridFunction *nodes;
-   mfem::ParFiniteElementSpace *mfes;
-   mfem::Array<int> domain_attributes;
-   const mfem::IntegrationRule *ir;
+mfem::ParGridFunction *nodes;
+mfem::ParFiniteElementSpace *mfes;
+mfem::Array<int> domain_attributes;
+const mfem::IntegrationRule *ir;

-   mutable mfem::ParGridFunction displ;
-   mutable mfem::ParGridFunction veloc;
-   mutable mfem::ParGridFunction accel;
+mutable mfem::ParGridFunction displ;
+mutable mfem::ParGridFunction veloc;
+mutable mfem::ParGridFunction accel;

-   mutable mfem::BlockVector sol;
-   mutable mfem::BlockVector rhs;
-   mutable mfem::BlockVector tmp;
-   mutable mfem::Vector  res;
+mutable mfem::BlockVector sol;
+mutable mfem::BlockVector rhs;
+mutable mfem::BlockVector tmp;
+mutable mfem::Vector  res;

-   mfem::Array<int> block_true_offsets;
+mfem::Array<int> block_true_offsets;

-   std::unique_ptr<mfem::future::UniformParameterSpace> ups;
-   std::unique_ptr<mfem::QuadratureSpace> qs;
-   std::unique_ptr<mfem::FaceQuadratureSpace> fqs;
+std::unique_ptr<mfem::future::UniformParameterSpace> ups;
+std::unique_ptr<mfem::QuadratureSpace> qs;
+std::unique_ptr<mfem::FaceQuadratureSpace> fqs;

-   // linear elasticty coefficients in dFEM form
-   // l1, m1 - material 1
-   // l2, m2 - material 2
-   std::unique_ptr<mfem::CoefficientVector> l1, l2;
-   std::unique_ptr<mfem::CoefficientVector> m1, m2;
+// linear elasticty coefficients in dFEM form
+// l1, m1 - material 1
+// l2, m2 - material 2
+std::unique_ptr<mfem::CoefficientVector> l1, l2;
+std::unique_ptr<mfem::CoefficientVector> m1, m2;

-   // linear elasticity damping coefficients in dFEM form
-   // strain proportional damping for material 1 and 2
-   std::unique_ptr<mfem::CoefficientVector> dl1, dl2;
-   std::unique_ptr<mfem::CoefficientVector> dm1, dm2;
+// linear elasticity damping coefficients in dFEM form
+// strain proportional damping for material 1 and 2
+std::unique_ptr<mfem::CoefficientVector> dl1, dl2;
+std::unique_ptr<mfem::CoefficientVector> dm1, dm2;

-   // damping mass coefficients in dFEM form
-   // (velocity proportional damping) for material 1 and 2
-   std::unique_ptr<mfem::CoefficientVector> cm1, cm2;
+// damping mass coefficients in dFEM form 
+// (velocity proportional damping) for material 1 and 2
+std::unique_ptr<mfem::CoefficientVector> cm1, cm2;

-   // density coefficients in dFEM form for material 1 and 2
-   std::unique_ptr<mfem::CoefficientVector> dens1, dens2;
-   mfem::Coefficient *cdens1, *cdens2;
+// density coefficients in dFEM form for material 1 and 2
+std::unique_ptr<mfem::CoefficientVector> dens1, dens2;
+mfem::Coefficient *cdens1, *cdens2;

-   // density coefficient for topology optimization
-   std::shared_ptr<mfem::CoefficientVector> density;
-   mfem::Coefficient *cdensity;
+// density coefficient for topology optimization
+std::shared_ptr<mfem::CoefficientVector> density;
+mfem::Coefficient *cdensity;

-   static constexpr int FDispl = 0; //grid function displacement
-   static constexpr int FVeloc = 1; //grid function velocity
-   // elasticity Coefficient Vectors
-   static constexpr int Lambda1 = 2, Lambda2 = 3, Mu1 = 4, Mu2 = 5;
-   // damping Coefficient Vectors
-   static constexpr int DLambda1 = 6, DMu1 = 7, DLambda2 = 8, DMu2 = 9;
-   static constexpr int CMass1 = 10, CMass2 = 11; // damping mass coeff vectors
-   static constexpr int Dens1 = 12, Dens2 = 13; // density Coefficient Vectors
-   // density for topology optimization
-   static constexpr int Density = 14; // coefficient vector
-   static constexpr int Coords = 15; // coordinates grid function

-   // DFEM forward related definitions
-   std::unique_ptr<mfem::future::DifferentiableOperator> dfem_forward_op;
-   std::unique_ptr<mfem::future::DifferentiableOperator> dfem_mass_op;
-   std::unique_ptr<mfem::future::DifferentiableOperator> dfem_damp_op;
-   std::unique_ptr<mfem::future::DifferentiableOperator> dfem_vol_force_op;

-   std::unique_ptr<mfem::HypreParMatrix> M_lor;
+static constexpr int FDispl = 0; //grid function displacement
+static constexpr int FVeloc = 1; //grid function velocity
+// elasticity Coefficient Vectors
+static constexpr int Lambda1 = 2, Lambda2 = 3, Mu1 = 4, Mu2 = 5;
+// damping Coefficient Vectors
+static constexpr int DLambda1 = 6, DMu1 = 7, DLambda2 = 8, DMu2 = 9;
+static constexpr int CMass1 = 10, CMass2 = 11; // damping mass coeff vectors
+static constexpr int Dens1 = 12, Dens2 = 13; // density Coefficient Vectors   
+// density for topology optimization
+static constexpr int Density = 14; // coefficient vector
+static constexpr int Coords = 15; // coordinates grid function

-   std::unique_ptr<mfem::CGSolver> cg;
-   std::unique_ptr<mfem::HypreBoomerAMG> amg;
+// DFEM forward related definitions
+std::unique_ptr<mfem::future::DifferentiableOperator> dfem_forward_op;
+std::unique_ptr<mfem::future::DifferentiableOperator> dfem_mass_op;
+std::unique_ptr<mfem::future::DifferentiableOperator> dfem_damp_op;
+std::unique_ptr<mfem::future::DifferentiableOperator> dfem_vol_force_op;

-   // zero bdrs
-   std::set<int> zero_bdrs;
+std::unique_ptr<mfem::HypreParMatrix> M_lor;

-   // time dependent memory vector for dynamic force
-   // the force is applied on all boundary attributes in bdr_loads_markers
-   std::set<int> bdr_loads_markers;
-   mutable mfem::Vector bdr_force_mem; // [0] - time, [1] - period, [2] - amplitude
+std::unique_ptr<mfem::CGSolver> cg;
+std::unique_ptr<mfem::HypreBoomerAMG> amg;

-   // volumetric force parameters
-   // A*sin(2*pi*t/T)*cos^n (pi (t-t_0)/L)
-   // [0] - time, [1] - period, [2] - amplitude, [3] - radius
-   // [4],[5],[6] - point coordinates of the center of the force application
-   // [7] - L total train length - could be proportional to the period [1]
-   // [8] - t_0  center of the train
-   // [9] - n the envelope power
-   mutable mfem::Vector vol_force_mem;
+// zero bdrs 
+std::set<int> zero_bdrs;

-   // zero bdr dofs - constructed during the corrsponding Assemble calls
-   mfem::Array<int> ess_tdof_list;
+// time dependent memory vector for dynamic force
+// the force is applied on all boundary attributes in bdr_loads_markers
+std::set<int> bdr_loads_markers;
+mutable mfem::Vector bdr_force_mem; // [0] - time, [1] - period, [2] - amplitude

-   //objective/constraints integrand
-   //obj->Mult(x,y)
-   //takes state vector s  and returns y which consists of multiple objectives/constraints
-   std::shared_ptr<mfem::Operator> obj;
+// volumetric force parameters
+// A*sin(2*pi*t/T)*cos^n (pi (t-t_0)/L)
+// [0] - time, [1] - period, [2] - amplitude, [3] - radius
+// [4],[5],[6] - point coordinates of the center of the force application
+// [7] - L total train length - could be proportional to the period [1]
+// [8] - t_0  center of the train
+// [9] - n the envelope power
+mutable mfem::Vector vol_force_mem; 
+
+// zero bdr dofs - constructed during the corrsponding Assemble calls
+mfem::Array<int> ess_tdof_list;

-   struct
-   {
-      mfem::real_t Tfinal;
-      mfem::Vector states[2];
-      mfem::real_t times[2];
-      int ind=0;
-      int max_states=2;
-   } adjoint_data;
 };


-///////////////////////////////////////////////////////////////////////////////
+
+
 class ExampleObjectiveIntegrand: public mfem::Operator
 {
 public:

-   ExampleObjectiveIntegrand(mfem::ParFiniteElementSpace* fes_,
-                             std::shared_ptr<mfem::Coefficient> objc);
+    ExampleObjectiveIntegrand(mfem::ParFiniteElementSpace* fes_, 
+                                std::shared_ptr<mfem::Coefficient> objc={});

-   void SetCoefficients( std::shared_ptr<mfem::Coefficient> objc);

-   //evaluates the QoIs y[1] for a given state x[2 x fes_->GetTrueVSize()]
-   virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+    void SetCoefficients( std::shared_ptr<mfem::Coefficient> objc);

-   mfem::real_t EvalScalar(const mfem::Vector &x) const
-   {
-      mfem::Vector y(1);
-      Mult(x, y);
-      return y[0];
-   }
+    //evaluates the QoIs y[1] for a given state x[2 x fes_->GetTrueVSize()]
+    virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override;

-   void EvalGradient(const mfem::Vector &x, mfem::Vector &grad) const;
+    mfem::real_t EvalScalar(const mfem::Vector &x) const
+    {
+        mfem::Vector y(1);
+        Mult(x, y);
+        return y[0];
+    }
+
+    void EvalGradient(const mfem::Vector &x, mfem::Vector &grad) const;
+    

 private:
-   std::shared_ptr<mfem::Coefficient> co;
-   mfem::ParFiniteElementSpace* fes;
-   mutable mfem::ParGridFunction disp;
-   mutable mfem::ParGridFunction velo;
+    std::shared_ptr<mfem::Coefficient> co;
+    mfem::ParFiniteElementSpace* fes;
+    mutable mfem::ParGridFunction disp;
+    mutable mfem::ParGridFunction velo;

-   mfem::Operator* grad;
+    mfem::Operator* grad;
+    
+    mfem::Array<int> block_true_offsets;

-   static constexpr int FDispl = 0; //grid function displacement
-   static constexpr int FVeloc = 1; //grid function velocity
-   static constexpr int Density = 14; // coefficient vector
-   static constexpr int Coords = 15; // coordinates grid function
+    mutable mfem::Vector res; 

-   // DFEM related definitions (3 objectives)
-   std::unique_ptr<mfem::future::DifferentiableOperator> obj;
+    std::unique_ptr<mfem::ParBilinearForm> mass;

-   // density coefficient for computing the objective function
-   std::shared_ptr<mfem::CoefficientVector> density;
-
-   mfem::Array<int> block_true_offsets;
-
-   //uniform parameter space
-   std::unique_ptr<mfem::future::UniformParameterSpace> ups;
-   //quadrature space for the coefficient
-   std::unique_ptr<mfem::QuadratureSpace> qs;
-
-   mfem::ParGridFunction *nodes;
-   mfem::ParFiniteElementSpace *mfes;
-   mfem::Array<int> domain_attributes;
-   const mfem::IntegrationRule *ir;
-
-   mutable mfem::Vector res;
 };

 #endif // LINEAR_ELASTICITY_HPP
@@ -0,0 +1,545 @@
+#include "linear_elasticity.hpp"
+
+#include "../chpt/dynamic_checkpointing.hpp"
+#include "../chpt/fixed_slot_checkpoint_storage.hpp"
+
+#include <cmath>
+#include <iomanip>
+#include <cstring>
+
+using namespace mfem;
+using namespace std;
+
+constexpr auto MESH_TRI = MFEM_SOURCE_DIR "/miniapps/mtop/examples/dyn_hex2d_tri.msh";
+constexpr auto MESH_QUAD = MFEM_SOURCE_DIR "/miniapps/mtop/examples/dyn_hex2d_quad.msh";
+
+struct State
+{
+   mfem::real_t time = 0.0; //time of the state
+   mfem::real_t dt=0.0;
+   mfem::real_t obj  = 0.0; //accumulated objective 
+   mfem::BlockVector v; //state of the system
+};
+
+// Snapshot = *view* (non-owning) used only during Store() packing and Read() callback
+struct StateSnapshotView
+{
+   mfem::real_t time = 0.0;
+   mfem::real_t dt=0.0;
+   mfem::real_t obj  = 0.0;
+
+   // Points to n*sizeof(real_t) bytes:
+   // - during Store(): points to current State::v data (host)
+   // - during Read(): points into the storage slot bytes (valid only during callback)
+   const unsigned char *v_bytes = nullptr;
+};
+
+class StateSnapshotViewPacker
+{
+public:
+   explicit StateSnapshotViewPacker(int n) : n_(n)
+   {
+      MFEM_VERIFY(n_ > 0, "StateSnapshotViewPacker: n must be > 0.");
+   }
+
+   std::size_t SlotBytes() const
+   {
+      return (std::size_t)(3 + n_) * sizeof(mfem::real_t);
+   }
+
+   void Pack(const StateSnapshotView &s, void *dst) const
+   {
+      MFEM_VERIFY(dst != nullptr, "Pack: dst is null.");
+      MFEM_VERIFY(s.v_bytes != nullptr, "Pack: snapshot v_bytes is null.");
+
+      unsigned char *b = static_cast<unsigned char*>(dst);
+
+      std::memcpy(b + 0*sizeof(mfem::real_t), &s.time, sizeof(mfem::real_t));
+      std::memcpy(b + 1*sizeof(mfem::real_t), &s.dt,   sizeof(mfem::real_t));
+      std::memcpy(b + 2*sizeof(mfem::real_t), &s.obj,  sizeof(mfem::real_t));
+
+      std::memcpy(b + 3*sizeof(mfem::real_t),
+                  s.v_bytes,
+                  (std::size_t)n_ * sizeof(mfem::real_t));
+   }
+
+   // Important: Unpack returns a *view* pointing into src bytes (no allocation).
+   void Unpack(const void *src, StateSnapshotView &out) const
+   {
+      MFEM_VERIFY(src != nullptr, "Unpack: src is null.");
+
+      const unsigned char *b = static_cast<const unsigned char*>(src);
+
+      std::memcpy(&out.time, b + 0*sizeof(mfem::real_t), sizeof(mfem::real_t));
+      std::memcpy(&out.dt,   b + 1*sizeof(mfem::real_t), sizeof(mfem::real_t));
+      std::memcpy(&out.obj,  b + 2*sizeof(mfem::real_t), sizeof(mfem::real_t));
+
+      out.v_bytes = b + 3*sizeof(mfem::real_t);
+   }
+
+   int VectorSize() const { return n_; }
+
+private:
+   int n_ = 0;
+};
+
+struct AdjState
+{
+   mfem::real_t time;
+   mfem::real_t obj;
+   mfem::BlockVector adj;
+};
+
+
+///////////////////////////////////////////////////////////////////////////////
+/// \brief The IsoElasticyLambdaCoeff class converts E modulus of elasticity
+/// and Poisson's ratio to Lame's lambda coefficient
+class IsoElasticyLambdaCoeff : public mfem::Coefficient
+{
+   mfem::Coefficient *E, *nu;
+
+public:
+   /// Constructor - takes as inputs E modulus and Poisson's ratio
+   IsoElasticyLambdaCoeff(mfem::Coefficient *E,
+                          mfem::Coefficient *nu):
+      E(E), nu(nu) { }
+
+   /// Evaluates the Lame's lambda coefficient
+   real_t Eval(mfem::ElementTransformation &T,
+               const mfem::IntegrationPoint &ip) override
+   {
+      const real_t EE = E->Eval(T, ip);
+      const real_t nn = nu->Eval(T, ip);
+      constexpr auto Lambda = [](const real_t E, const real_t nu)
+      {
+         return E * nu / (1.0 + nu) / (1.0 - 2.0 * nu);
+      };
+      return Lambda(EE, nn);
+   }
+};
+
+
+
+///////////////////////////////////////////////////////////////////////////////
+/// \brief The IsoElasticySchearCoeff class converts E modulus of elasticity
+/// and Poisson's ratio to Shear coefficient
+///
+class IsoElasticySchearCoeff : public mfem::Coefficient
+{
+   mfem::Coefficient *E, *nu;
+
+public:
+   /// Constructor - takes as inputs E modulus and Poisson's ratio
+   IsoElasticySchearCoeff(mfem::Coefficient *E_, mfem::Coefficient *nu_):
+      E(E_), nu(nu_) { }
+
+   /// Evaluates the shear coefficient coefficient
+   real_t Eval(mfem::ElementTransformation &T,
+               const mfem::IntegrationPoint &ip) override
+   {
+      const real_t EE = E->Eval(T, ip);
+      const real_t nn = nu->Eval(T, ip);
+      constexpr auto Schear = [](const real_t E, const real_t nu)
+      {
+         return E / (2.0 * (1.0 + nu));
+      };
+      return Schear(EE, nn);
+   }
+};
+
+
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = MESH_QUAD;
+   const char *device_config = "cpu";
+   int order = 3;
+   bool mesh_tri = false;
+   bool mesh_quad = false;
+   int par_ref_levels = 1;
+   bool paraview = true;
+   bool visualization = true;
+   int ode_solver_type = 4;
+   real_t Tfinal = 0.07;
+   real_t dt = 0.005;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&mesh_tri, "-tri", "--triangular", "-no-tri",
+                  "--no-triangular", "Enable or not triangular mesh.");
+   args.AddOption(&mesh_quad, "-quad", "--quadrilateral", "-no-quad",
+                  "--no-quadrilateral", "Enable or not quadrilateral mesh.");
+   args.AddOption(&par_ref_levels, "-prl", "--par-ref-levels",
+                  "Number of parallel mesh refinement levels.");
+   args.AddOption(&paraview, "-pa", "--paraview", "-no-pa",
+                  "--no-paraview", "Enable or not Paraview output.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization", "Enable or not visualization.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  ODESolver::Types.c_str());
+
+   args.AddOption(&Tfinal, "-T",   
+                    "--tfinal",   "Terminate when accumulated time reaches Tfinal.");
+   args.AddOption(&dt,    "-dt", "--dt",  "Time step.");   
+
+   args.ParseCheck();
+
+   // Enable hardware devices such as GPUs, and programming models such as
+   // CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+   
+   // Read the (serial) mesh from the given mesh file on all processors.  We
+   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   // and volume meshes with the same code.
+   Mesh mesh(mesh_tri ? MESH_TRI : mesh_quad ? MESH_QUAD : mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+   
+   // Refine the serial mesh on all processors to increase the resolution. In
+   // this example we do 'ref_levels' of uniform refinement. We choose
+   // 'ref_levels' to be the largest number that gives a final mesh with no
+   // more than 1000 elements.
+   {
+      const int ref_levels =
+         (int)floor(log(1000. / mesh.GetNE()) / log(2.) / dim);
+      for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
+   }
+   if (Mpi::Root())
+   {
+      std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+   }
+
+   // Define a parallel mesh by a partitioning of the serial mesh. Refine
+   // this mesh further in parallel to increase the resolution. Once the
+   // parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   for (int l = 0; l < par_ref_levels; l++) { pmesh.UniformRefinement(); }   
+
+   // Allocate the time dependent linear elasticity operator
+   LinearElasticityTimeDependentOperator lin_elasticity_op(pmesh, order);
+
+   // Set the material coefficients
+   ConstantCoefficient rho_coef(0.5); // density coefficient for topology optimization
+   
+   // Set elasticity coefficients for material 1 and 2
+   ConstantCoefficient E1(0.1);
+   ConstantCoefficient E2(1.0);
+   ConstantCoefficient nu1(0.3);
+   ConstantCoefficient nu2(0.3);
+
+   //Lame coefficients
+   IsoElasticyLambdaCoeff lambda1(&E1, &nu1);
+   IsoElasticySchearCoeff mu1(&E1, &nu1);
+   IsoElasticyLambdaCoeff lambda2(&E2, &nu2);
+   IsoElasticySchearCoeff mu2(&E2, &nu2);
+
+   // Set density coefficients for material 1 and 2
+   ConstantCoefficient dens1_coef(0.5);
+   ConstantCoefficient dens2_coef(1.0);
+
+   // Set damping coefficients
+   ProductCoefficient cm1_coef(0.02, dens1_coef);
+   ProductCoefficient cm2_coef(0.02, dens2_coef);
+
+   ProductCoefficient cl1_coef(0.01, lambda1);
+   ProductCoefficient cmu1_coef(0.01, mu1);
+   ProductCoefficient cl2_coef(0.01, lambda2);
+   ProductCoefficient cmu2_coef(0.01, mu2);
+
+   lin_elasticity_op.SetElasticityCoefficients(lambda1, mu1, lambda2, mu2);
+   
+   lin_elasticity_op.SetDensityMaterialCoefficients(dens1_coef, dens2_coef);
+   
+   lin_elasticity_op.SetDampingMaterialCoefficients(cm1_coef, cm2_coef);
+   lin_elasticity_op.SetDampingMaterialCoefficients(cl1_coef, cmu1_coef,
+                                                    cl2_coef, cmu2_coef);
+
+   lin_elasticity_op.SetDensity(rho_coef);
+   
+   //set bottom bdr to zero (both the velocities and the displacements)
+   lin_elasticity_op.SetZeroBdr(1);
+
+   lin_elasticity_op.SetVolForce(1.0 /*period*/, 1.0 /*amplitude*/, 0.2 /*radius*/,
+                                 0.0 /*x center*/ , 0.0 /*y center*/, 0.0 /*z center*/,
+                                 5.0 /* train length*/, 2.5 /*center of the train*/, 2.0 /*power*/);
+
+   lin_elasticity_op.AssembleExplicit();
+
+   // test mult explicit
+   {
+      BlockVector tst; tst.Update(lin_elasticity_op.GetTrueBlockOffsets()); 
+      tst=0.0; //tst.Randomize();
+      tst.UseDevice(true); tst.Read();
+
+      BlockVector grd; grd.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      grd=0.0;
+      lin_elasticity_op.Mult(tst,grd);
+      lin_elasticity_op.GetVelocity().SetFromTrueDofs(grd.GetBlock(1));
+   }   
+
+   std::shared_ptr<ExampleObjectiveIntegrand> eobj=
+      std::make_shared<ExampleObjectiveIntegrand>(lin_elasticity_op.GetFESpace());
+   //set the objective for the integration process                                                   
+
+   ParaViewDataCollection paraview_dc("isoel", &pmesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.RegisterField("disp", &(lin_elasticity_op.GetDisplacement()));
+   paraview_dc.RegisterField("velo", &(lin_elasticity_op.GetVelocity()));
+
+   // 4. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   unique_ptr<ODESolver> ode_solver = ODESolver::Select(ode_solver_type);
+
+   lin_elasticity_op.SetTime(0.0);
+   ode_solver->Init(lin_elasticity_op);
+
+   //random perturbation of the initial conditions
+   BlockVector p; p.Update(lin_elasticity_op.GetTrueBlockOffsets()); p=0.0;
+   Vector tmpv(p.GetBlock(0)); tmpv=0.0;
+
+   //Forward computations
+   {
+      int s=10; //number of snapshots to be stored by the checkpointing process
+
+      //define the packer object
+      int n=lin_elasticity_op.GetState().Size();
+      StateSnapshotViewPacker packer(n);
+
+      // storage stores StateSnapshotView snapshots using fixed-size slots
+      using Storage=mfem::FixedSlotMemoryCheckpointStorage<StateSnapshotView, StateSnapshotViewPacker>;
+      Storage storage(s, packer);
+
+      // Snapshot type is StateSnapshotView
+      using Checkpointing = mfem::DynamicCheckpointing<StateSnapshotView, Storage>;
+      Checkpointing ckpt(s, storage);
+
+      // Returns view of the State and avoids data transfer
+      auto make_snapshot = [&](const State &u) -> StateSnapshotView
+      {
+         MFEM_VERIFY(u.v.Size() == n, "make_snapshot: State.v size changed!");
+
+         // Ensure host access if MFEM device is in use:
+         const mfem::real_t *vh = u.v.HostRead();
+
+         StateSnapshotView snap;
+         snap.time = u.time;
+         snap.obj  = u.obj;
+         snap.v_bytes = reinterpret_cast<const unsigned char*>(vh);
+         return snap;
+      };
+
+       //Transfers data from the snaphot view to the State u_out.
+      auto restore_snapshot = [&](const StateSnapshotView &snap, State &u_out)
+      {
+         u_out.time = snap.time;
+         u_out.obj  = snap.obj;
+
+         if (u_out.v.Size() != n) { u_out.v.SetSize(n); }
+
+         mfem::real_t *vh = u_out.v.HostWrite();
+         std::memcpy(vh,
+                     snap.v_bytes,
+                     (std::size_t)n * sizeof(mfem::real_t));
+
+         //make sure that the date is on the device            
+         u_out.v.Read(true);
+      };
+
+      using Step = mfem::DynamicCheckpointing<StateSnapshotView, Storage>::Step;
+
+
+      //execute one integration step
+      auto primal_step = [&](State &u_st, Step i)
+      {
+         //begin with curent state u_st
+         real_t t=u_st.time;
+         real_t ldt=dt;
+         real_t obj=u_st.obj;
+
+         //make sure the integration does not overjump Tfinal
+         if((t+ldt)>Tfinal){
+            ldt=Tfinal-t;
+         }
+
+         //advance u_st
+         ode_solver->Step(u_st.v,t,ldt);
+         //update objective 
+         obj=eobj->EvalScalar(u_st.v);
+
+         //return updated u_st
+         u_st.dt=t-u_st.time;
+         u_st.time=t;
+         u_st.obj=obj;
+      };
+
+
+      auto adjoint_step = [&](AdjState &adj_st, const State &u_st, Step i)
+      {
+         real_t t=u_st.time;
+         real_t ldt=u_st.dt;
+         real_t obj=u_st.obj;
+
+         ode_solver->EnableAdjoint(mfem::ODESolver::AdjointMode::Discrete);
+         ode_solver->SetSolution(u_st.v,t);
+         ode_solver->AdjointStep(adj_st.adj,t,ldt);
+         ode_solver->EnableAdjoint(mfem::ODESolver::AdjointMode::None);
+         adj_st.time=t;
+         adj_st.obj=obj;
+      };
+
+      BlockVector x0; x0.Update(lin_elasticity_op.GetTrueBlockOffsets()); 
+      x0=0.0;
+
+   //optimization loop   
+   for(int ii=0;ii<1;ii++){
+
+      ParaViewDataCollection paraview_dc("frw", &pmesh);
+      paraview_dc.SetPrefixPath("ParaView");
+      paraview_dc.SetLevelsOfDetail(order);
+      paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      paraview_dc.SetHighOrderOutput(true);
+      paraview_dc.RegisterField("disp", &(lin_elasticity_op.GetDisplacement()));
+      paraview_dc.RegisterField("velo", &(lin_elasticity_op.GetVelocity()));
+
+
+      State u;
+      //u.v.SetSize(lin_elasticity_op.GetState().Size());
+      u.v.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      u.obj=0.0;
+      u.time=0.0;
+      u.dt=0.0;
+
+      //set initial state to 0
+      u.v=0.0;
+      u.v.Set(1.0,x0);
+
+      // Forward sweep (unknown number of steps) 
+      real_t t = 0.0;
+      Step i = 0;
+
+      paraview_dc.SetCycle(0);
+      paraview_dc.SetTime(t);
+      lin_elasticity_op.GetVelocity().SetFromTrueDofs(u.v.GetBlock(1));
+      lin_elasticity_op.GetDisplacement().SetFromTrueDofs(u.v.GetBlock(0));
+      paraview_dc.Save();
+
+
+      while(t<Tfinal)
+      {
+         ckpt.ForwardStep(i,u, primal_step, make_snapshot);
+         t=u.time;
+         ++i;
+
+         if((i%5)==0){
+            paraview_dc.SetCycle(i+1);
+            paraview_dc.SetTime(t);
+            lin_elasticity_op.GetVelocity().SetFromTrueDofs(u.v.GetBlock(1));
+            lin_elasticity_op.GetDisplacement().SetFromTrueDofs(u.v.GetBlock(0));
+            paraview_dc.Save();
+         }
+      }
+
+      //Do one more time step without checkpointing
+      //advance u_st
+      ode_solver->Step(u.v,t,dt);
+      //update objective 
+      real_t obj=eobj->EvalScalar(u.v);
+
+      {
+         paraview_dc.SetCycle(i+1);
+         paraview_dc.SetTime(t);
+         lin_elasticity_op.GetVelocity().SetFromTrueDofs(u.v.GetBlock(1));
+         lin_elasticity_op.GetDisplacement().SetFromTrueDofs(u.v.GetBlock(0));
+         paraview_dc.Save();
+      }
+
+      if(Mpi::Root())
+      {
+         std::cout<<"Total number of steps i="<<i<<" t="<<t<<" obj="<<obj<<std::endl;
+         std::cout<<"Start adjoint steps!"<<std::endl;
+      }
+
+      ParGridFunction adisp(lin_elasticity_op.GetDisplacement());
+      ParGridFunction avelo(lin_elasticity_op.GetVelocity());
+      ParaViewDataCollection paraview_ac("adj", &pmesh);
+      paraview_ac.SetPrefixPath("ParaView");
+      paraview_ac.SetLevelsOfDetail(order);
+      paraview_ac.SetDataFormat(VTKFormat::BINARY);
+      paraview_ac.SetHighOrderOutput(true);
+      paraview_ac.RegisterField("adisp", &(adisp));
+      paraview_ac.RegisterField("avelo", &(avelo));
+
+      // define the adjoint state
+      AdjState adj_st; 
+      //adj_st.adj.SetSize(lin_elasticity_op.GetState().Size()); 
+      adj_st.adj.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      adj_st.adj=0.0;
+      eobj->EvalGradient(u.v,adj_st.adj);
+      adj_st.obj=obj;
+      adj_st.time=t;
+
+      paraview_ac.SetCycle(i);
+      paraview_ac.SetTime(t);
+      lin_elasticity_op.MultInvMass(adj_st.adj.GetBlock(0),tmpv);
+      adisp.SetFromTrueDofs(tmpv);
+      lin_elasticity_op.MultInvMass(adj_st.adj.GetBlock(1),tmpv);
+      avelo.SetFromTrueDofs(tmpv);
+      paraview_ac.Save();
+
+      
+      //step backward
+      const Step m=i;
+      
+      //tmp state for stepping backward
+      State u_wrk;
+      //u_wrk.v.SetSize(lin_elasticity_op.GetState().Size());
+      u_wrk.v.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      u_wrk.obj=0.0;
+      u_wrk.time=0.0;
+      u_wrk.dt=0.0;
+
+      for (Step i = m - 1; i >= 0; --i)
+      {
+         ckpt.BackwardStep(i,adj_st, u_wrk,
+                           primal_step, adjoint_step,
+                           make_snapshot, restore_snapshot);
+
+         if((i%5)==0){
+            paraview_ac.SetCycle(i);
+            paraview_ac.SetTime(adj_st.time);
+            //adisp.SetFromTrueDofs(adj_st.adj.GetBlock(0));
+            //avelo.SetFromTrueDofs(adj_st.adj.GetBlock(1));
+            lin_elasticity_op.MultInvMass(adj_st.adj.GetBlock(0),tmpv);
+            adisp.SetFromTrueDofs(tmpv);
+            lin_elasticity_op.MultInvMass(adj_st.adj.GetBlock(1),tmpv);
+            avelo.SetFromTrueDofs(tmpv);
+            paraview_ac.Save();
+         }
+
+         if (i == 0) { break; }
+      }
+
+      lin_elasticity_op.MultInvMass(adj_st.adj.GetBlock(0),tmpv);
+      x0.GetBlock(0).Add(-0.1,tmpv);
+      lin_elasticity_op.MultInvMass(adj_st.adj.GetBlock(1),tmpv);
+      x0.GetBlock(1).Add(-0.1,tmpv);
+   }
+
+   }
+
+   return EXIT_SUCCESS;
+}
@@ -1,26 +1,17 @@
 #include "linear_elasticity.hpp"

-// #include "../chpt/dynamic_checkpointing.hpp"
-// #include "../chpt/fixed_slot_checkpoint_storage.hpp"
+#include "../chpt/dynamic_checkpointing.hpp"
+#include "../chpt/fixed_slot_checkpoint_storage.hpp"

 #include <cmath>
+#include <iomanip>
 #include <cstring>

 using namespace mfem;
 using namespace std;

-#ifdef NVTX_DEBUG_HPP
-#undef NVTX_COLOR
-#define NVTX_COLOR ::nvtx::kCyan
-#include NVTX_DEBUG_HPP
-#else
-#define dbg(...)
-#endif
-
-constexpr auto MESH_TRI = MFEM_SOURCE_DIR
-                          "/miniapps/mtop/examples/dyn_hex2d_tri.msh";
-constexpr auto MESH_QUAD = MFEM_SOURCE_DIR
-                           "/miniapps/mtop/examples/dyn_hex2d_quad.msh";
+constexpr auto MESH_TRI = MFEM_SOURCE_DIR "/miniapps/mtop/examples/dyn_hex2d_tri.msh";
+constexpr auto MESH_QUAD = MFEM_SOURCE_DIR "/miniapps/mtop/examples/dyn_hex2d_quad.msh";

 struct State
 {
@@ -157,8 +148,6 @@ public:

 int main(int argc, char *argv[])
 {
-   dbg();
-
   // Initialize MPI and HYPRE.
   Mpi::Init();
   Hypre::Init();
@@ -170,8 +159,7 @@ int main(int argc, char *argv[])
   bool mesh_tri = false;
   bool mesh_quad = false;
   int par_ref_levels = 1;
-   int max_steps = 1000;
-   bool paraview = false;
+   bool paraview = true;
   bool visualization = true;
   int ode_solver_type = 4;

@@ -188,29 +176,26 @@ int main(int argc, char *argv[])
                  "--no-quadrilateral", "Enable or not quadrilateral mesh.");
   args.AddOption(&par_ref_levels, "-prl", "--par-ref-levels",
                  "Number of parallel mesh refinement levels.");
-   args.AddOption(&max_steps, "-ms", "--max-steps",
-                  "Maximum number of time steps.");
   args.AddOption(&paraview, "-pa", "--paraview", "-no-pa",
                  "--no-paraview", "Enable or not Paraview output.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization", "Enable or not visualization.");
-   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                  ODESolver::Types.c_str());

   args.ParseCheck();
-   dbg("mesh_file: {}", mesh_file);

   // Enable hardware devices such as GPUs, and programming models such as
   // CUDA, OCCA, RAJA and OpenMP based on command line options.
   Device device(device_config);
   if (Mpi::Root()) { device.Print(); }
-
+   
   // Read the (serial) mesh from the given mesh file on all processors.  We
   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
   // and volume meshes with the same code.
   Mesh mesh(mesh_tri ? MESH_TRI : mesh_quad ? MESH_QUAD : mesh_file, 1, 1);
   const int dim = mesh.Dimension();
-
+   
   // Refine the serial mesh on all processors to increase the resolution. In
   // this example we do 'ref_levels' of uniform refinement. We choose
   // 'ref_levels' to be the largest number that gives a final mesh with no
@@ -223,7 +208,6 @@ int main(int argc, char *argv[])
   if (Mpi::Root())
   {
      std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
-      dbg("Number of elements: {}", mesh.GetNE());
   }

   // Define a parallel mesh by a partitioning of the serial mesh. Refine
@@ -231,15 +215,14 @@ int main(int argc, char *argv[])
   // parallel mesh is defined, the serial mesh can be deleted.
   ParMesh pmesh(MPI_COMM_WORLD, mesh);
   mesh.Clear();
-   for (int l = 0; l < par_ref_levels; l++) { pmesh.UniformRefinement(); }
+   for (int l = 0; l < par_ref_levels; l++) { pmesh.UniformRefinement(); }   

   // Allocate the time dependent linear elasticity operator
   LinearElasticityTimeDependentOperator lin_elasticity_op(pmesh, order);

   // Set the material coefficients
-   ConstantCoefficient rho_coef(
-      0.5); // density coefficient for topology optimization
-
+   ConstantCoefficient rho_coef(0.5); // density coefficient for topology optimization
+   
   // Set elasticity coefficients for material 1 and 2
   ConstantCoefficient E1(0.1);
   ConstantCoefficient E2(1.0);
@@ -266,27 +249,27 @@ int main(int argc, char *argv[])
   ProductCoefficient cmu2_coef(0.01, mu2);

   lin_elasticity_op.SetElasticityCoefficients(lambda1, mu1, lambda2, mu2);
-
+   
   lin_elasticity_op.SetDensityMaterialCoefficients(dens1_coef, dens2_coef);
-
+   
   lin_elasticity_op.SetDampingMaterialCoefficients(cm1_coef, cm2_coef);
   lin_elasticity_op.SetDampingMaterialCoefficients(cl1_coef, cmu1_coef,
                                                    cl2_coef, cmu2_coef);

   lin_elasticity_op.SetDensity(rho_coef);
-
+   
   //set bottom bdr to zero (both the velocities and the displacements)
   lin_elasticity_op.SetZeroBdr(1);

   lin_elasticity_op.SetVolForce(1.0 /*period*/, 1.0 /*amplitude*/, 0.2 /*radius*/,
-                                 0.0 /*x center*/, 0.0 /*y center*/, 0.0 /*z center*/,
+                                 0.0 /*x center*/ , 0.0 /*y center*/, 0.0 /*z center*/,
                                 5.0 /* train length*/, 2.5 /*center of the train*/, 2.0 /*power*/);

   lin_elasticity_op.AssembleExplicit();

   // test mult explicit
   {
-      BlockVector tst; tst.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      BlockVector tst; tst.Update(lin_elasticity_op.GetTrueBlockOffsets()); 
      tst=0.0; //tst.Randomize();
      tst.UseDevice(true); tst.Read();

@@ -294,7 +277,7 @@ int main(int argc, char *argv[])
      grd=0.0;
      lin_elasticity_op.Mult(tst,grd);
      lin_elasticity_op.GetVelocity().SetFromTrueDofs(grd.GetBlock(1));
-   }
+   }   

   ParaViewDataCollection paraview_dc("isoel", &pmesh);
   paraview_dc.SetPrefixPath("ParaView");
@@ -380,7 +363,7 @@ int main(int argc, char *argv[])
      };


-      // Initial condition
+      // Initial condition 
      StateCheckPoint spt; spt.obj=-1.0; spt.time=-1.0; spt.state=(lin_elasticity_op.GetState());

      Step i=0;
@@ -398,7 +381,7 @@ int main(int argc, char *argv[])
         std::cout<<" Total number of steps="<<m<<std::endl;
      }

-      AdjState ast; ast.obj=1.0; ast.time=spt.time;
+      AdjState ast; ast.obj=1.0; ast.time=spt.time; 
      ast.adj=(lin_elasticity_op.GetState());
      ast.grd=(lin_elasticity_op.GetState());

@@ -409,19 +392,77 @@ int main(int argc, char *argv[])
         }

         ckpt.BackwardStep(j, ast, spt, primal_step, adjoint_step, make_snapshot, restore_snapshot);
-
+         
         if (j == 0) { break; }
      }
-
+      


   }
      */

+   //Test Jacobian transpose   
+   {
+      BlockVector x; x.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      BlockVector w; w.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      BlockVector d; d.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      BlockVector g; g.Update(lin_elasticity_op.GetTrueBlockOffsets());
+      x.Randomize();
+
+      real_t t = 0.0;
+
+      ExampleObjectiveIntegrand* obj=new ExampleObjectiveIntegrand(lin_elasticity_op.GetFESpace());
+
+      lin_elasticity_op.Mult(x,d);
+      real_t oo=obj->EvalScalar(d);
+      if(Mpi::Root())
+      {
+         std::cout<<"t="<<t<<" oo="<<oo<<std::endl;
+      }
+
+      obj->EvalGradient(d,w);
+      lin_elasticity_op.JacobianMultTranspose(x,w,g); //the gradient is in g
+
+      //FD check
+      {
+         BlockVector rnd; rnd.Update(lin_elasticity_op.GetTrueBlockOffsets());
+         rnd.Randomize();
+         real_t sca=1.0;
+
+         real_t ipr=mfem::InnerProduct(pmesh.GetComm(),rnd,g);
+
+         for(int i=0;i<20;i++){
+            w.Set(sca,rnd);
+            w.Add(1.0,x);
+            lin_elasticity_op.Mult(w,d);
+            real_t co=obj->EvalScalar(d);
+
+            w.Set(-sca,rnd);
+            w.Add(1.0,x);
+            lin_elasticity_op.Mult(w,d);
+            real_t mo=obj->EvalScalar(d);
+
+
+            if(Mpi::Root())
+            {
+               std::cout<<"s="<<sca<<" o="<<oo<<" c="<<co<<" gr="
+                                 <<ipr<<" fd="<<(co-oo)/sca
+                                 <<" cd="<<(co-mo)/(2.0*sca)<<std::endl;
+            }
+            sca=sca/2;
+         }   
+
+
+      }
+
+      delete obj;
+
+   }



   //test time integration
+   /*
   {
      real_t t = 0.0;

@@ -429,6 +470,16 @@ int main(int argc, char *argv[])
      //set initial conditions at time t
      tsol=0.0;

+      ExampleObjectiveIntegrand* obj=new ExampleObjectiveIntegrand(lin_elasticity_op.GetFESpace());
+
+      real_t oo=obj->EvalScalar(tsol);
+      if(Mpi::Root())
+      {
+         std::cout<<"t="<<t<<" oo="<<oo<<std::endl;
+      }
+
+      delete obj;
+
      // 4. Define the ODE solver used for time integration. Several explicit
      //    Runge-Kutta methods are available.
      unique_ptr<ODESolver> ode_solver = ODESolver::Select(ode_solver_type);
@@ -445,47 +496,34 @@ int main(int argc, char *argv[])
      real_t dt_real = 0.005;
      //ode_solver->Run(tsol, t, dt_real, 1.0);

-      for (int i=0; i< std::min(max_steps, 6000); i++)
-      {
-         dbg("#{}", i);
+      for(int i=0;i<2000; i++){
         ode_solver->Step(tsol, t, dt_real);

         if (Mpi::Root())
         {
-            std::cout << "t: " << t << std::endl;
+            std::cout << "t: " << t << std::endl; 
         }

-         if (paraview && (i%5)==0)
-         {
-            dbg("Paraview output at step {} time {}", i, t);
+         if((i%5)==0){
            paraview_dc.SetCycle(i+1);
            paraview_dc.SetTime(t);
            lin_elasticity_op.GetVelocity().SetFromTrueDofs(tsol.GetBlock(1));
            lin_elasticity_op.GetDisplacement().SetFromTrueDofs(tsol.GetBlock(0));
            paraview_dc.Save();
         }
-
-         if (socketstream glvis; visualization && (i%10)==0 &&
-             ((glvis.open("localhost", 19916)),
-              (glvis.is_open() ? glvis.precision(8) : 0),
-              (glvis.is_open() ? (glvis << "mgjR" << std::endl, 0) : 0),
-              (glvis.is_open())))
-         {
-            dbg("GLVis output at step {} time {}", i, t);
-            glvis << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
-
-            lin_elasticity_op.GetVelocity().SetFromTrueDofs(tsol.GetBlock(1));
-            ParGridFunction velocity(lin_elasticity_op.GetVelocity());
-
-            // lin_elasticity_op.GetDisplacement().SetFromTrueDofs(tsol.GetBlock(0));
-            // ParGridFunction displacement(lin_elasticity_op.GetDisplacement());
-
-            glvis << "solution\n" << pmesh << velocity << std::flush;
-            glvis << "keys g\n" << std::flush; // pause
-         }
      }

+      obj=new ExampleObjectiveIntegrand(lin_elasticity_op.GetFESpace());
+      oo=obj->EvalScalar(tsol);
+      if(Mpi::Root())
+      {
+         std::cout<<"t="<<t<<" oo="<<oo<<std::endl;
+      }
+
+      delete obj;
+
   }
+   */   

   /*
   if (paraview)
@@ -0,0 +1,538 @@
+#include "mfem.hpp"
+
+
+class InterpCoeff: public mfem::Coefficient
+{
+public:
+    InterpCoeff(std::shared_ptr<mfem::Coefficient> c1_,
+                std::shared_ptr<mfem::Coefficient> c2_,
+                std::shared_ptr<mfem::Coefficient> rho_)
+                : c1(c1_), c2(c2_), rho(rho_)
+                {
+
+                }
+
+    virtual mfem::real_t 	Eval (mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
+    {
+        mfem::real_t v1=c1->Eval(T,ip);
+        mfem::real_t v2=c2->Eval(T,ip);
+        mfem::real_t r=rho->Eval(T,ip);
+
+        return v1*r+v2*(1.0-r);
+    }
+
+private:
+    std::shared_ptr<mfem::Coefficient> c1;
+    std::shared_ptr<mfem::Coefficient> c2;
+    std::shared_ptr<mfem::Coefficient> rho;
+};
+
+
+namespace mfem{
+
+class PAExplicitDiffusionOperator : public TimeDependentOperator
+{
+private:
+   ParFiniteElementSpace &pfes;
+   Array<int> ess_tdof_list;
+
+   // Partially assembled bilinear forms (no assembled sparse matrices)
+   ParBilinearForm k_form; // diffusion stiffness
+
+   // Minv(i) = 1 / diag(M)(i) on TRUE dofs
+   Vector Minv;
+   mutable Vector tmpv;
+
+public:
+   /// ess_bdr: boundary attribute marker (size = pmesh->bdr_attributes.Max()), 1 -> essential (Dirichlet)
+   /// kappa: diffusivity coefficient
+   PAExplicitDiffusionOperator(ParFiniteElementSpace &pfes_,
+                               const Array<int> &ess_tdof_list_,
+                               Coefficient &kappa)
+      : TimeDependentOperator(pfes_.GetTrueVSize(), 0.0, TimeDependentOperator::EXPLICIT),
+        pfes(pfes_),
+        k_form(&pfes_),
+        Minv(height),
+        tmpv(height)
+   {
+      // Essential TRUE dofs
+      // pfes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      ess_tdof_list=ess_tdof_list_;
+
+      // --- Mass operator (PA) ---
+      auto mifi=new MassIntegrator();
+      int order=pfes.GetOrder(0);
+      IntegrationRules gll_rules(0, Quadrature1D::GaussLobatto);
+      const IntegrationRule &ir_ni = gll_rules.Get(pfes.GetParMesh()->GetTypicalElementGeometry(),
+                                                        2 * order - 1);
+      mifi->SetIntRule(&ir_ni);
+
+      ParBilinearForm m_form(&pfes); // mass
+      m_form.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      m_form.AddDomainIntegrator(mifi);
+      m_form.Assemble();
+
+      // Diagonal on TRUE dofs (works with PA in ParBilinearForm)
+      Vector Mdiag(height);
+      m_form.AssembleDiagonal(Mdiag);
+
+      Minv = Mdiag;
+      Minv.Reciprocal();
+
+      // --- Diffusion operator (PA) ---
+      k_form.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      k_form.AddDomainIntegrator(new DiffusionIntegrator(kappa));
+      k_form.Assemble();
+   }
+
+   /// Compute du_dt = -Minv .* (K u)
+   void Mult(const Vector &u, Vector &du_dt) const override
+   {
+      du_dt = 0.0;
+
+      // du_dt += -1 * (P^T K_local P) u  (true-dof action; PA supported)
+      k_form.TrueAddMult(u, du_dt, -1.0);
+
+      // Apply diagonal inverse mass
+      du_dt *= Minv;
+
+      // Strongly enforce Dirichlet: derivative is zero on essential tdofs
+      if (ess_tdof_list.Size())
+      {
+         du_dt.SetSubVector(ess_tdof_list, 0.0);
+      }
+   }
+
+   const Array<int> &GetEssentialTrueDofs() const { return ess_tdof_list; }
+
+   /// y = (df/dx(x,t))^T * w
+   void JacobianMultTranspose(const Vector &x,
+                                      const Vector &w,
+                                      Vector &y) const override
+   {
+      y=0.0;
+      tmpv.Set(1.0,w);
+      if (ess_tdof_list.Size()){
+        tmpv.SetSubVector(ess_tdof_list, 0.0);}
+      // Apply diagonal inverse mass
+      tmpv*=Minv;
+      // y += -1 * (P^T K_local P) u  (true-dof action; PA supported)
+      k_form.TrueAddMult(tmpv, y, -1.0); 
+   }
+};    
+
+};
+
+class DiffusionTDOP:public mfem::TimeDependentOperator
+{
+public:
+    DiffusionTDOP(mfem::ParMesh* mesh_, int vorder =1, int dorder=1):order(vorder)
+    {
+        mesh=mesh_;
+        int dim=mesh->Dimension();
+        fec.reset(new mfem::H1_FECollection(vorder,dim));
+        fes.reset(new mfem::ParFiniteElementSpace(mesh,fec.get()));
+
+        dfec.reset(new mfem::H1_FECollection(dorder,dim));
+        dfes.reset(new mfem::ParFiniteElementSpace(mesh,dfec.get()));
+
+        gfdens.SetSpace(dfes.get());
+        
+
+        //the state vector of the TDOP consists of [solution, density, obj]
+        siz_u=fes->GetTrueVSize();  //solution
+        siz_d=dfes->GetTrueVSize(); //design
+        siz_f=1; //objective;
+
+        block_true_offsets.SetSize(4);
+        block_true_offsets[0] = 0;
+        block_true_offsets[1] = siz_u;
+        block_true_offsets[2] = siz_d;
+        block_true_offsets[3] = siz_f;
+        block_true_offsets.PartialSum();
+
+        //set the width and the height of the operator
+        this->width=  block_true_offsets[3];
+        this->height= block_true_offsets[3];
+    }
+
+    const mfem::ParFiniteElementSpace* GetStateFES(){ return fes.get(); }
+    const mfem::ParFiniteElementSpace* GetDesignFEM(){ return dfes.get(); }
+
+    virtual 
+    ~DiffusionTDOP()
+    {
+
+    }
+
+    void SetLoad(int ind, std::shared_ptr<mfem::Coefficient> l)
+    {
+        load[ind]=l;
+    }
+
+    void SetZeroBC(int ind)
+    {
+        zero_bc.insert(ind);
+    }
+
+    // free the allocated bilinear forms, matrices and solvers
+    void Reset()
+    {
+
+    }
+
+    /// Asembles the load, the mass and the stiffness matrices for
+    /// a give true vector tvdens with densities varying from 0 to 1
+    /// and a given time t.
+    void Assemble(const mfem::Vector& tvdens, mfem::real_t t)
+    {
+        
+        gfdens.SetFromTrueDofs(tvdens);
+        gfc.reset(new mfem::GridFunctionCoefficient(&gfdens));
+
+        if(nullptr==kbf.get())
+        {
+            // set constrained dofs
+            SetEssTDofs(ess_tdofv);
+            //allocate the system matrices and solvers
+            //allocate diagonal mass matrix and its inverse
+            {
+                mfem::IntegrationRules gll_rules(0, mfem::Quadrature1D::GaussLobatto);
+                const mfem::IntegrationRule &ir_ni = gll_rules.Get(mesh->GetTypicalElementGeometry(),
+                                                        2 * order - 1);
+
+                mfem::ParBilinearForm mform(fes.get());
+                mform.SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+                cs.reset(new InterpCoeff(cs1,cs2,gfc));
+                auto *mblfi = new mfem::MassIntegrator(*cs);
+                mblfi->SetIntRule(&ir_ni);
+                mform.AddDomainIntegrator(mblfi);
+
+                mfem::Vector diag(fes->GetTrueVSize());
+                mform.AssembleDiagonal(diag);
+                minv=diag;
+                minv.Reciprocal();
+                // Ensure essential dofs stay fixed: zero inverse there (extra safety)
+                if (ess_tdofv.Size())
+                {
+                    minv.SetSubVector(ess_tdofv, 0.0);
+                }
+            }
+
+            //allocate the stiffness matrix
+            {
+                kap.reset(new InterpCoeff(kap1,kap2,gfc));
+                kbf.reset(new mfem::ParBilinearForm(fes.get()));
+                kbf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+                kbf->AddDomainIntegrator(new mfem::DiffusionIntegrator(*kap));
+                kbf->Assemble();
+            }
+        }
+
+        //allocate the RHS
+
+
+    }
+
+    virtual
+    void Mult(const mfem::Vector& x, mfem::Vector& y) const override
+    {
+        mfem::BlockVector bx(const_cast<mfem::Vector&>(x), block_true_offsets);
+        mfem::BlockVector by(y, block_true_offsets);
+        //compute the time derivatives
+        //compute the time derivative of the objective 
+    }
+
+
+    void SetEssTDofs(mfem::Array<int>& ess_dofs)
+    {
+        ess_dofs.DeleteAll();
+
+        mfem::Array<int> ess_bdr(mesh->bdr_attributes.Max());
+        ess_bdr=0;
+        for (auto it=zero_bc.begin(); it!=zero_bc.end(); ++it)
+        {
+            int attr = *it;
+            ess_bdr[attr-1] = 1;
+        }
+
+        fes->GetEssentialTrueDofs(ess_bdr,ess_dofs);
+    }
+
+    /// Sets Dirichlet dofs to zero
+    void SetEssTDofs(mfem::Vector& v) const
+    {
+        mfem::Array<int> loc_tdofs;
+        mfem::Array<int> ess_bdr(mesh->bdr_attributes.Max());
+        ess_bdr=0;
+        for (auto it=zero_bc.begin(); it!=zero_bc.end(); ++it)
+        {
+            int attr = *it;
+            ess_bdr[attr-1] = 1;
+        }
+
+        fes->GetEssentialTrueDofs(ess_bdr, loc_tdofs);
+
+        for (int j=0; j<loc_tdofs.Size(); j++)
+        {
+            v[loc_tdofs[j]]=0.0;
+        }
+    }
+
+private:
+
+int order;
+mfem::ParMesh* mesh;
+std::unique_ptr<mfem::FiniteElementCollection> fec;
+std::unique_ptr<mfem::ParFiniteElementSpace> fes;
+
+std::unique_ptr<mfem::FiniteElementCollection> dfec;
+std::unique_ptr<mfem::ParFiniteElementSpace> dfes;
+
+std::unique_ptr<mfem::HypreBoomerAMG> prec;
+std::unique_ptr<mfem::CGSolver> ls;
+
+std::shared_ptr<mfem::Coefficient> kap1, kap2;
+std::shared_ptr<mfem::Coefficient> cs1, cs2;
+std::shared_ptr<InterpCoeff> cs,kap;
+std::shared_ptr<mfem::GridFunctionCoefficient> gfc;
+mfem::ParGridFunction gfdens;
+
+std::map<int, std::shared_ptr<mfem::Coefficient>> load;
+std::set<int> zero_bc;
+// holds the constrained DOFs
+mfem::Array<int> ess_tdofv;
+
+mfem::Array<int> block_true_offsets;
+int siz_u;
+int siz_d;
+int siz_f;
+
+std::unique_ptr<mfem::ParBilinearForm> kbf;
+mfem::Vector minv;
+
+};
+
+using namespace std;
+using namespace mfem;
+
+constexpr auto MESH_TRI = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_tri.mesh";
+constexpr auto MESH_QUAD = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_quad.mesh";
+
+int main(int argc, char *argv[])
+{
+    // Initialize MPI and HYPRE.
+   Mpi::Init();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = MESH_QUAD;
+   const char *device_config = "cpu";
+   int order = 2;
+   bool pa = false;
+   bool dfem = false;
+   bool mesh_tri = false;
+   bool mesh_quad = false;
+   int par_ref_levels = 1;
+   bool paraview = false;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&dfem, "-dfem", "--dFEM", "-no-dfem", "--no-dFEM",
+                  "Enable or not dFEM.");
+   args.AddOption(&mesh_tri, "-tri", "--triangular", "-no-tri",
+                  "--no-triangular", "Enable or not triangular mesh.");
+   args.AddOption(&mesh_quad, "-quad", "--quadrilateral", "-no-quad",
+                  "--no-quadrilateral", "Enable or not quadrilateral mesh.");
+   args.AddOption(&par_ref_levels, "-prl", "--par-ref-levels",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Enable or not Paraview visualization");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+   MFEM_VERIFY(!(pa && dfem), "pa and dfem cannot be both set");
+
+   // Enable hardware devices such as GPUs, and programming models such as
+   // CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   // Read the (serial) mesh from the given mesh file on all processors.  We
+   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   // and volume meshes with the same code.
+   Mesh mesh(mesh_tri ? MESH_TRI : mesh_quad ? MESH_QUAD : mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+
+   // Refine the serial mesh on all processors to increase the resolution. In
+   // this example we do 'ref_levels' of uniform refinement. We choose
+   // 'ref_levels' to be the largest number that gives a final mesh with no
+   // more than 1000 elements.
+   {
+      const int ref_levels =
+         (int)floor(log(1000. / mesh.GetNE()) / log(2.) / dim);
+      for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
+   }
+   if (Mpi::Root())
+   {
+      std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+   }
+
+   // Define a parallel mesh by a partitioning of the serial mesh. Refine
+   // this mesh further in parallel to increase the resolution. Once the
+   // parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   for (int l = 0; l < par_ref_levels; l++) { pmesh.UniformRefinement(); }
+
+   mfem::H1_FECollection fec(order, pmesh.Dimension());
+   mfem::ParFiniteElementSpace pfes(&pmesh, &fec);
+
+   mfem::Array<int> ess_tdof_list;
+   {
+       // Dirichlet boundary marker (example: all boundary attributes)
+       mfem::Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+       ess_bdr = 0;
+       ess_bdr[9]=1.0;
+       pfes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+
+   mfem::ConstantCoefficient kappa(1.0);
+
+   mfem::PAExplicitDiffusionOperator oper(pfes, ess_tdof_list, kappa);
+
+   // Initial condition as a ParGridFunction
+   mfem::ParGridFunction u_gf(&pfes);
+   u_gf = 0.0; 
+   // replace with projection
+   {
+        mfem::FunctionCoefficient fc([](const Vector &x) -> real_t
+        {   return std::sin(4.0*x[0]*M_PI)*std::sin(4.0*x[1]*M_PI); });
+        //project
+        u_gf.ProjectCoefficient(fc);
+   }
+
+
+   // True dof vector state
+   mfem::Vector u;
+   u_gf.GetTrueDofs(u);
+
+   // Enforce homogeneous Dirichlet initially
+   u.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
+   u_gf.SetFromTrueDofs(u);
+
+   // set paraview output
+   ParaViewDataCollection paraview_dc("tdiff", &pmesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.RegisterField("temp", &u_gf);
+
+   // Pick an explicit solver
+   mfem::RK4Solver ode;
+   ode.Init(oper);
+
+   double t = 0.0;
+   double dt = 1e-5; // diffusion CFL ~ O(h^2); must be small
+
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetTime(t);
+   paraview_dc.Save();
+
+
+   /*
+   for (int ti = 0; ti < 10; ti++)
+   {
+      ode.Step(u, t, dt);
+
+      // Keep Dirichlet dofs pinned (avoid drift)
+      u.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
+
+      u_gf.SetFromTrueDofs(u);
+      paraview_dc.SetCycle(ti+1);
+      paraview_dc.SetTime(t);
+      paraview_dc.Save();
+   }
+   */
+
+
+   //test the RK4 adjoint
+   {
+        Vector x;
+        {
+            mfem::FunctionCoefficient fc([](const Vector &x) -> real_t
+            {   return std::sin(4.0*x[0]*M_PI)*std::sin(4.0*x[1]*M_PI); });
+            //project
+            u_gf.ProjectCoefficient(fc);
+        }
+        u_gf.GetTrueDofs(x);
+        // Enforce homogeneous Dirichlet initially
+        x.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
+        u_gf.SetFromTrueDofs(x);
+        u.Set(1.0,x);
+
+        Vector vone(x); vone=1.0;
+
+        //compute one step
+        t=0.0;
+        ode.Step(x, t, dt);
+        //compute objective 
+        real_t obj=mfem::InnerProduct(pmesh.GetComm(),x,vone);
+
+        //compute the adjoint
+        mfem::Vector lam(u.Size()); lam=0.0;
+        {
+            ode.EnableAdjoint(mfem::ODESolver::AdjointMode::Discrete);
+            ode.SetSolution(u,0.0); //set the solution at t=0.0;
+            lam.Set(1.0,vone);
+            ode.AdjointStep(lam,t,dt); //on exit lam is the adjoint at t-dt
+        }
+
+        mfem::Vector rnd(u.Size()); rnd.Randomize();
+        rnd.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
+
+        real_t iprod=mfem::InnerProduct(pmesh.GetComm(),rnd,lam);
+
+        real_t sca=1.0;
+        for(int i=0;i<20;i++){
+            x.Set(sca,rnd);
+            x.Add(1.0,u);
+            t=0.0;
+            ode.Step(x, t, dt);
+            real_t pobj=mfem::InnerProduct(pmesh.GetComm(),x,vone);
+
+            x.Set(-sca,rnd);
+            x.Add(1.0,u);
+            t=0.0;
+            ode.Step(x, t, dt);
+            real_t mobj=mfem::InnerProduct(pmesh.GetComm(),x,vone);
+
+            if(mfem::Mpi::Root())
+            {
+                std::cout<<" scale="<<sca<<" "<<" o="<<obj<<" p="<<pobj
+                            <<" do="<<(pobj-mobj)/(2.0*sca)
+                            <<" oo="<<(obj-mobj)/sca
+                            <<" to="<<iprod<<std::endl;
+            }
+
+            sca=sca/10.0;
+        }
+
+   }
+
+
+
+   return EXIT_SUCCESS;
+
+}
@@ -0,0 +1,25 @@
+# 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.
+
+if (MFEM_USE_MPI)
+  list(APPEND PLASMA_COMMON_SOURCES)
+
+  list(APPEND PLASMA_COMMON_HEADERS
+  plasma.hpp)
+
+  convert_filenames_to_full_paths(PLASMA_COMMON_SOURCES)
+  convert_filenames_to_full_paths(PLASMA_COMMON_HEADERS)
+
+  set(PLASMA_COMMON_FILES
+      EXTRA_SOURCES ${PLASMA_COMMON_SOURCES}
+      EXTRA_HEADERS ${PLASMA_COMMON_HEADERS})
+
+endif()
@@ -0,0 +1,85 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS = 
+PAR_MINIAPPS = 
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all lib-common clean clean-build clean-exec
+.PRECIOUS: %.o
+
+COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
+
+# If MFEM_SHARED is set, add the ../common rpath
+COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+COMMON_O=
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(MINIAPPS)
+
+# Rules for building the miniapps
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_O) $(COMMON_LIB) \
+        $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+# Rule for building lib-common
+lib-common:
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Specific execution options
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+
@@ -0,0 +1,62 @@
+// 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.
+
+#ifndef MFEM_PLASMA_HPP
+#define MFEM_PLASMA_HPP
+
+#include <cmath>
+#include <complex>
+
+namespace mfem
+{
+
+namespace plasma
+{
+
+// Physical Constants
+
+// Permittivity of Free Space (units F/m)
+static const real_t epsilon0_ = 8.8541878176e-12;
+
+// Permeability of Free Space (units H/m)
+static const real_t mu0_ = 4.0e-7 * M_PI;
+
+// Speed of light in Free Space (units m/s)
+static const real_t c0_ = 1.0 / sqrt(epsilon0_ * mu0_);
+
+// Impedance of Free Space (units Ohm)
+static const real_t Z0_ = sqrt(mu0_ / epsilon0_);
+
+static const real_t q_     = 1.602176634e-19; // Elementary charge in coulombs
+static const real_t eV_    = 1.602176634e-19; // 1 eV in Joules
+static const real_t amu_   = 1.660539040e-27; // Atomic mass unit in kilograms
+static const real_t me_kg_ = 9.10938356e-31;  // Mass of electron in kilograms
+static const real_t me_u_  = 5.4857990907e-4; // Mass of electron in a.m.u
+
+/**
+   Returns the cyclotron frequency in radians/second
+   m is the mass in a.m.u
+   q is the charge in units of elementary electric charge
+   B is the magnetic field magnitude in tesla
+ */
+inline real_t cyclotronFrequency(real_t B, real_t m, real_t q)
+{
+   return fabs(q * q_ * B / (m * amu_));
+}
+
+typedef std::complex<real_t> complex_t;
+
+} // namespace plasma
+
+} // namespace mfem
+
+#endif // MFEM_PLASMA_HPP
+
@@ -295,8 +295,13 @@ namespace Catch {
 // Otherwise all supported compilers support COUNTER macro,
 // but user still might want to turn it off
 #if ( !defined(__JETBRAINS_IDE__) || __JETBRAINS_IDE__ >= 20170300L )
+#if ( !(defined(__clang__) && __clang_major__ >= 22 ) )
+// don't use __COUNTER__ if compiling with clang 22+ to avoid compiler warning
+// https://github.com/llvm/llvm-project/pull/162662
+// TODO: can enable if building with C2y
    #define CATCH_INTERNAL_CONFIG_COUNTER
 #endif
+#endif

 ////////////////////////////////////////////////////////////////////////////////

@@ -117,3 +117,52 @@ TEST_CASE("Vector FE Face Restriction", "[FaceRestriction]")
   gf2 -= gf;
   REQUIRE(gf2.Normlinf() == MFEM_Approx(0.0));
 }
+
+#ifdef MFEM_USE_MPI
+
+TEST_CASE("L2 Face Restriction", "[FaceRestriction][Parallel]")
+{
+   const int dim = GENERATE(2, 3);
+   constexpr int nx = 3;
+   constexpr int order = 2;
+   constexpr int vdim = 2;
+   const Ordering::Type ordering = GENERATE(Ordering::byNODES, Ordering::byVDIM);
+
+   Mesh serial_mesh = MakeCartesianMesh(nx, dim);
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+
+   L2_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace fes(&mesh, &fec, vdim, ordering);
+
+   auto *R = fes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,
+                                    FaceType::Interior);
+
+   Vector vals({1.0, 2.0});
+   VectorConstantCoefficient coeff(vals);
+
+   ParGridFunction gf(&fes);
+   gf.ProjectCoefficient(coeff);
+
+   Vector face_vec(R->Height());
+   R->Mult(gf, face_vec);
+
+   const int nf = mesh.GetNFbyType(FaceType::Interior);
+   const int face_dofs = fes.GetTypicalTraceElement()->GetDof();
+   auto h_face_vec = Reshape(face_vec.HostRead(), face_dofs, vdim, 2, nf);
+
+   for (int f = 0; f < nf; ++f)
+   {
+      for (int m = 0; m < 2; ++m)
+      {
+         for (int c = 0; c < vdim; ++c)
+         {
+            for (int i = 0; i < face_dofs; ++i)
+            {
+               REQUIRE(h_face_vec(i, c, m, f) == vals[c]);
+            }
+         }
+      }
+   }
+}
+
+#endif
@@ -281,8 +281,10 @@ TEST_CASE("Nedelec Segment Finite Element",
            REQUIRE( fe.GetRangeType()      == (int) FiniteElement::VECTOR  );
            REQUIRE( fe.GetMapType()        == (int) FiniteElement::H_CURL  );
            REQUIRE( fe.GetDerivType()      == (int) FiniteElement::NONE    );
-            REQUIRE( fe.GetDerivRangeType() == (int) FiniteElement::SCALAR  );
-            REQUIRE( fe.GetDerivMapType()   == (int) FiniteElement::INTEGRAL);
+            REQUIRE( fe.GetDerivRangeType() ==
+                     (int) FiniteElement::UNKNOWN_RANGE_TYPE);
+            REQUIRE( fe.GetDerivMapType()   ==
+                     (int) FiniteElement::UNKNOWN_MAP_TYPE);
         }
      }
      SECTION("Sizes for p = " + std::to_string(p))
Author	SHA1	Message	Date
blaz	68f9353cd5	mod	2026-04-01 23:42:36 -07:00
blaz	2bb5f45bb6	mod	2026-04-01 23:39:06 -07:00
blaz	cf8743dd5e	mod	2026-04-01 23:33:36 -07:00
blaz	f5b9c0afeb	mod	2026-04-01 23:25:07 -07:00
blaz	cf5447d879	mod	2026-04-01 23:20:44 -07:00
bslazarov@github.com	30558cd938	mod	2026-03-31 10:10:51 -07:00
blaz	77511860ef	mod	2026-03-26 23:04:50 -07:00
blaz	d017504e5a	mod	2026-03-26 22:58:50 -07:00
blaz	88fa1d60cf	mod	2026-03-26 22:52:21 -07:00
blaz	da83d4961a	opt	2026-03-26 22:44:05 -07:00
blaz	db1e90d02c	mod	2026-03-26 21:08:01 -07:00
blaz	010600b18a	mod	2026-03-26 20:46:37 -07:00
blaz	613d1a3bcb	mod	2026-03-26 20:11:24 -07:00
blaz	e07760dfc2	tst grad	2026-03-26 19:20:43 -07:00
blaz	acf49b667b	explicit time adjoint example	2026-03-26 10:54:59 -07:00
bslazarov@github.com	09ba64bb29	fix	2026-03-24 21:53:59 -07:00
blaz	bb31355d93	fix	2026-03-24 19:35:50 -07:00
blaz	a483e50abf	ch	2026-03-24 19:19:18 -07:00
blaz	2f55b418ed	assembly level	2026-03-24 19:02:49 -07:00
blaz	313d74fa1f	objective added	2026-03-24 18:55:30 -07:00
blaz	829905c8fc	Merge branch 'ode_adjoints' into checkpointing_examples	2026-03-23 21:07:56 -07:00
blaz	272cbe594d	Merge branch 'master' into ode_adjoints	2026-03-23 21:07:21 -07:00
Tzanio Kolev	55e42eeefe	Merge pull request #5238 from mfem/face-nbr-restr-vdim-bugfix Fix bug in ParL2FaceRestriction with vdim > 1	2026-03-22 10:24:53 -07:00
blaz	e58de8335a	changes	2026-03-16 10:22:32 -07:00
blaz	a327f612be	tst file	2026-03-14 10:32:33 -07:00
Tzanio Kolev	fa89c5e98c	Merge pull request #4856 from mfem/phys-range-dim Range and curl dimension in physical space	2026-03-13 07:44:40 -07:00
Tzanio Kolev	0980bda63b	Merge pull request #5215 from balay/barry/update-for-petsc-v3.25-PetscCtx Update to change in PETSc API (in v3.25) for PetscCtx and PetscCtxRt	2026-03-13 07:44:06 -07:00
Will Pazner	a1758e51e5	Merge remote-tracking branch 'origin/master' into face-nbr-restr-vdim-bugfix	2026-03-12 18:02:19 -07:00
Will Pazner	ccf84aab7c	Use constexpr in unit test	2026-03-12 18:01:40 -07:00
Will Pazner	82abd48bba	Merge pull request #5080 from mfem/cmake-config CMake config.mk for CUDA and HIP	2026-03-11 12:18:25 -04:00
blaz	fad600cbd0	RK4 adjoint	2026-03-10 22:43:49 -07:00
Tzanio Kolev	d4374a9d5f	Merge branch 'master' into cmake-config	2026-03-10 11:08:23 -07:00
Tzanio Kolev	dcd3a25730	Merge branch 'master' into barry/update-for-petsc-v3.25-PetscCtx	2026-03-10 11:01:46 -07:00
blaz	ec1ae6eaf5	initial not-ready version for RK4 adjoint	2026-03-10 10:52:05 -07:00
Andrew Ho	0f99528c62	Merge branch 'master' into phys-range-dim	2026-03-05 12:45:42 -08:00
Veselin Dobrev	ddfd74e899	Merge pull request #5255 from mfem/catch-tests fix clang compiler warning for __COUNTER__	2026-03-05 12:39:29 -08:00
Mark L. Stowell	0248720eeb	Merge branch 'master' into phys-range-dim	2026-03-05 09:35:14 -08:00
Andrew Ho	feded39641	Merge branch 'master' into catch-tests	2026-03-05 09:28:31 -08:00
Tzanio Kolev	09128b9a5d	Merge pull request #5240 from mfem/bugfix/chapan39/use-mfem-abort-kernel-in-device dfem integrate: use mfem abort kernel in device code	2026-03-04 09:55:23 -08:00
Tzanio Kolev	68383b462b	Merge pull request #5231 from mfem/plasma-dir-dev Plasma Miniapp Directory	2026-03-04 09:54:48 -08:00
Tzanio Kolev	24d5609585	Merge pull request #5212 from mfem/najlkin/fix-cmplx-assign [BUG] Complex grid function copy assignment	2026-03-04 09:54:25 -08:00
Andrew Ho	670a3f9a45	comment on why TPL_LIBRARIES is reversed twice	2026-03-03 11:43:02 -08:00
Andrew Ho	7baae02d65	Merge remote-tracking branch 'base/cmake-config' into cmake-config	2026-03-02 16:31:18 -08:00
Andrew Ho	728a0f313b	move cudart to MFEM_EXT_LIBS	2026-03-02 16:30:31 -08:00
Andrew Ho	1bb624e2a8	fix clang compiler warning for __COUNTER__	2026-03-02 14:05:28 -08:00
Tzanio Kolev	ee7ccd6464	Merge branch 'master' into plasma-dir-dev	2026-03-02 11:56:51 -08:00
Stowell, Mark L.	a3ae5a6f01	Changing copyright date to pass CI checks	2026-03-02 09:22:29 -08:00
Andrew HoandNuno Nobre	9243d00549	Update config/cmake/modules/MfemCmakeUtilities.cmake Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2026-02-28 14:53:22 -08:00
Andrew HoandNuno Nobre	4fe3db5a5f	Update config/cmake/modules/MfemCmakeUtilities.cmake Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2026-02-28 14:53:15 -08:00
Andrew Ho	55bb710cba	fixed wrong dir being marked as system	2026-02-27 14:02:24 -08:00
Andrew HoandNuno Nobre	7ad6939454	Update config/cmake/modules/MfemCmakeUtilities.cmake Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2026-02-27 07:06:16 -08:00
Barry Smith Nuno Nobre Satish Balay	60cc94e5a1	Update to use PetscCtxRt from (3,25,0), and cleanup duplicate code Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk> Co-authored-by: Satish Balay <balay@mcs.anl.gov>	2026-02-26 11:57:27 -06:00
Satish Balay	9122ac1839	update KSPMonitorFn usage for < (3,24,0)	2026-02-26 11:57:22 -06:00
Satish Balay	864186117d	update PetscCtxDestroyFn usage for < (3,23,0)	2026-02-26 11:56:15 -06:00
Ketan Mittal	35de169fd0	Merge branch 'master' into plasma-dir-dev	2026-02-26 09:42:51 -08:00
Andrew Ho	0a3184ab31	MFEM_EXPORT_GPU_CONFIG should export CPU config.mk when set to off	2026-02-24 11:56:38 -08:00
Mark L. Stowell	a438e09caf	Merge branch 'master' into plasma-dir-dev	2026-02-24 11:02:04 -08:00
Alex Tyler Chapman	ea03a86df2	Merge branch 'master' into bugfix/chapan39/use-mfem-abort-kernel-in-device	2026-02-24 10:32:30 -08:00
chapman39	6ef7a9e6fb	80 chars/ line	2026-02-24 10:32:18 -08:00
Alex Tyler Chapman	db7dd30d32	Merge branch 'master' into bugfix/chapan39/use-mfem-abort-kernel-in-device	2026-02-23 09:33:23 -08:00
chapman39	a1fe3a19b1	dfem integrate: use mfem abort kernel in device code	2026-02-19 17:45:11 -08:00
Satish Balay	213ccd7a4e	rework PetscContainerSetCtxDestroy() usage for < (3,23,0)	2026-02-18 15:51:42 -06:00
Will Pazner	8e78471fdf	Add comment about the shape of FaceNbrData	2026-02-18 09:13:16 -08:00
Will Pazner	c0f8501950	Add unit test for parallel L2 face restriction with vdim > 1	2026-02-18 09:13:03 -08:00
Will Pazner	c31510289f	Fix bug in ParL2FaceRestriction with vdim > 1 The layout of the FaceNbrData vector was not handled properly	2026-02-17 21:07:17 -08:00
Tzanio Kolev	2b14134496	Merge branch 'master' into najlkin/fix-cmplx-assign	2026-02-17 08:27:01 -08:00
Stowell, Mark L.	794a5fbfc2	Adding miniapps/plasma subdirectory to build system	2026-02-11 11:56:06 -08:00
Stowell, Mark L.	746a62f017	Adding plasma miniapp directory	2026-02-11 11:53:01 -08:00
Andrew Ho	26e9057f02	revert change, updated comment to why libdl gets special treatment	2026-01-30 07:27:04 -08:00
Jan Nikl	16dfa11f27	Minor docstring correction.	2026-01-28 15:45:20 -08:00
Jan Nikl	c7774e3c1c	Fixed complex grid function copy assignment.	2026-01-28 15:34:55 -08:00
Andrew Ho	194f3d8140	suggestions from Veselin	2026-01-21 11:45:20 -08:00
Mark L. Stowell	3f9b44a9cd	Merge branch 'master' into phys-range-dim	2026-01-21 10:38:56 -08:00
John Camier	caa973d6a0	Merge branch 'master' into cmake-config	2026-01-13 08:13:47 -08:00
Andrew Ho	ad40704e20	Merge branch 'master' into cmake-config	2026-01-06 11:55:40 -08:00
Andrew Ho	d3470c07c9	Merge branch 'master' into cmake-config	2025-12-16 12:05:27 -08:00
Andrew Ho	06a15cb7a9	missed one old unsetting of shared_link_flag	2025-12-01 17:02:05 -08:00
Andrew Ho	d19ff6c676	Merge branch 'master' into cmake-config	2025-12-01 12:33:24 -08:00
Andrew Ho	d85fbc6504	review suggestions	2025-11-25 14:35:38 -08:00
Andrew Ho	29346a87b6	Merge branch 'master' into cmake-config	2025-11-25 14:31:21 -05:00
Andrew Ho	3464f7a004	Merge branch 'master' into cmake-config	2025-10-28 11:08:38 -07:00
Andrew Ho	7de48e47ad	Merge branch 'master' into cmake-config	2025-10-24 09:24:50 -07:00
Andrew Ho	70814c640b	fixes for hip	2025-10-20 14:01:50 -07:00
Andrew Ho	e9d3ae80f7	remove debug printout	2025-10-20 13:30:20 -07:00
Andrew Ho	c8efc23c12	seems to be building external laghos now	2025-10-20 13:26:04 -07:00
Andrew Ho	f26eb33252	Merge remote-tracking branch 'base/cmake-gpu' into cmake-config	2025-10-20 10:25:58 -07:00
Andrew Ho	05e622f837	improving config.mk file generated by cmake to work with hip/cuda Still need to export compiler flags	2025-10-20 08:33:17 -07:00
Mark L. Stowell	e9f84b033f	Merge branch 'master' into phys-range-dim	2025-10-15 06:49:33 -07:00
Mark L. Stowell	83fd119b95	Merge branch 'master' into phys-range-dim	2025-07-01 10:16:58 -07:00
Mark L. Stowell	4aecb86d71	Merge branch 'master' into phys-range-dim	2025-05-19 17:51:08 -07:00
Stowell, Mark L.	776a4c1815	Updating unit tests	2025-05-11 14:40:05 -07:00
Stowell, Mark L.	c870d7dc1c	Using new MapType entries and implementing new GetPhys*Dim methods	2025-05-11 14:39:49 -07:00
Stowell, Mark L.	8519889074	Adding new MapType entries for R2D and R1D classes	2025-05-11 14:38:50 -07:00