Silenced tee.

Range and curl dimension in physical space

CMakeLists.txt

@@ -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()

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -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

examples/ex1p.cpp

@@ -66,6 +66,26 @@
 using namespace std;
 using namespace mfem;
 
+class VisServer
+{
+   const char *path{};
+   int port{-1};
+
+   struct Deleter { void operator()(FILE *); };
+   std::unique_ptr<FILE, Deleter> fglvis;
+
+   void Open();
+public:
+   VisServer(const char *path_, int port_ = 19916)
+      : path(path_), port(port_) { Open(); }
+
+   explicit VisServer() = default;
+
+   inline const char *GetPath() const { return path; }
+   inline int GetPort() const { return port; }
+   inline bool IsOpen() const { return (fglvis != nullptr); }
+};
+
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI and HYPRE.
@@ -82,6 +102,7 @@ int main(int argc, char *argv[])
    bool fa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   const char *visserver = "";
    bool algebraic_ceed = false;
 
    OptionsParser args(argc, argv);
@@ -106,6 +127,8 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&visserver, "-vs", "--visualization-server",
+                  "Path to GLVis binary to start own server.");
    args.Parse();
    if (!args.Good())
    {
@@ -301,6 +324,14 @@ int main(int argc, char *argv[])
    {
       char vishost[] = "localhost";
       int  visport   = 19916;
+
+      VisServer server;
+      if (strlen(visserver) > 0)
+      {
+         server = VisServer(visserver, visport);
+         visport = server.GetPort();
+      }
+
       socketstream sol_sock(vishost, visport);
       sol_sock << "parallel " << num_procs << " " << myid << "\n";
       sol_sock.precision(8);
@@ -315,3 +346,66 @@ int main(int argc, char *argv[])
 
    return 0;
 }
+
+void VisServer::Open()
+{
+   if (!path || strlen(path) <= 0) { return; }
+
+#ifdef MFEM_USE_MPI
+   if (Mpi::Root())
+#endif // MFEM_USE_MPI
+   {
+      // Open a log file
+      FILE *ftmp = popen("mktemp --tmpdir glvis-log.XXXXXX", "r");
+      constexpr size_t ssize = 256;
+      char stmp[ssize];
+      fgets(stmp, ssize, ftmp);
+      pclose(ftmp);
+      int len = strlen(stmp);
+      if (len > 0 && stmp[len-1] == '\n') { stmp[--len] = '\0'; }
+      std::cout << "Starting GLVis log in: " << stmp << std::endl;
+
+      // Start the server
+      std::stringstream ss;
+      ss << "trap '' SIGPIPE && " << path <<
+         " -no-pr -no-ex 2>&1 | tee \"" << stmp <<
+         "\" 2> /dev/null | grep -m 1 ^GLVIS_SERVER_PORT";
+      fglvis.reset(popen(ss.str().c_str(), "r"));
+
+      // Capture the port number
+      char line[ssize];
+      bool captured = false;
+      while (fgets(line, ssize, fglvis.get()))
+      {
+         if (strncmp(line, "GLVIS_SERVER_PORT=", 18) == 0)
+         {
+            sscanf(line, "GLVIS_SERVER_PORT=%d", &port);
+            captured = true;
+            break;
+         }
+      }
+      if (captured)
+      {
+         std::cout << "Started GLVis server at port: " << port << std::endl;
+      }
+      else
+      {
+         std::cerr << "GLVis server did not start normally." << std::endl;
+      }
+   }
+#ifdef MFEM_USE_MPI
+   MPI_Bcast(&port, 1, MPI_INT, 0, MPI_COMM_WORLD);
+#endif // MFEM_USE_MPI
+}
+
+void VisServer::Deleter::operator()(FILE *f)
+{
+   if (f != nullptr)
+   {
+      int ierr = pclose(f);
+      if (ierr != 0)
+      {
+         std::cerr << "GLVis server pclose() returns: " << ierr << std::endl;
+      }
+   }
+}

fem/complex_fem.hpp

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

fem/dfem/integrate.hpp

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

fem/fe/fe_base.cpp

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

fem/fe/fe_base.hpp

@@ -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

fem/fe/fe_nd.cpp

@@ -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)

fem/fe/fe_nd.hpp

@@ -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;
 

fem/fe/fe_rt.cpp

@@ -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)

fem/fe/fe_rt.hpp

@@ -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,

linalg/petsc.cpp

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

makefile

@@ -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)
 

miniapps/CMakeLists.txt

@@ -35,6 +35,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)

miniapps/plasma/CMakeLists.txt

@@ -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()

miniapps/plasma/makefile

@@ -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:
+

miniapps/plasma/plasma.hpp

@@ -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
+

tests/unit/catch.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
 
 ////////////////////////////////////////////////////////////////////////////////
 

tests/unit/fem/test_fe.cpp

@@ -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))