Fix style

Make unit test that reproduces enzyme error in debug build
Merge pull request #4856 from mfem/phys-range-dim
2026-03-23 08:44:36 -07:00 · 2026-03-21 10:22:20 -07:00 · 2026-03-13 07:44:40 -07:00 · 2026-03-13 07:44:06 -07:00 · 2026-03-11 12:18:25 -04:00 · 2026-03-10 11:08:23 -07:00
45 changed files with 917 additions and 385 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");
   }
 }

@@ -663,58 +663,59 @@ const
   #pragma omp critical (DofToQuad)
 #endif
   {
-      // If the new Dof2Quad is already present, e.g. added in a previous call
-      // or added by another omp thread, return.
+      // Do not run if the new Dof2Quad is already present, e.g. added in a
+      // previous call or added by another omp thread.
      if (DofToQuad::SearchArray(dof2quad_array, ir,
-                                 DofToQuad::LEXICOGRAPHIC_FULL))
-      { return; }
-
-      // Undo the native ordering which is what FiniteElement::GetDofToQuad
-      // returns.
-      auto *d2q_new = new DofToQuad(d2q);
-      d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
-      const int nqpt = ir.GetNPoints();
-
-      const int b_dim = (range_type == VECTOR) ? dim : 1;
-
-      for (int i = 0; i < nqpt; i++)
+                                 DofToQuad::LEXICOGRAPHIC_FULL) == nullptr)
      {
-         for (int d = 0; d < b_dim; d++)
+         // Undo the native ordering which is what FiniteElement::GetDofToQuad
+         // returns.
+         auto *d2q_new = new DofToQuad(d2q);
+         d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
+         const int nqpt = ir.GetNPoints();
+
+         const int b_dim = (range_type == VECTOR) ? dim : 1;
+
+         for (int i = 0; i < nqpt; i++)
         {
-            for (int j = 0; j < dof; j++)
+            for (int d = 0; d < b_dim; d++)
            {
-               const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
-               d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
-               d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
+                  d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
+                  d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               }
            }
         }
-      }

-      const int g_dim = [this]()
-      {
-         switch (deriv_type)
+         const int g_dim = [this]()
         {
-            case GRAD: return dim;
-            case DIV: return 1;
-            case CURL: return cdim;
-            default: return 0;
-         }
-      }();
-
-      for (int i = 0; i < nqpt; i++)
-      {
-         for (int d = 0; d < g_dim; d++)
-         {
-            for (int j = 0; j < dof; j++)
+            switch (deriv_type)
            {
-               const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
-               d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
-               d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               case GRAD: return dim;
+               case DIV: return 1;
+               case CURL: return cdim;
+               default: return 0;
+            }
+         }();
+
+         for (int i = 0; i < nqpt; i++)
+         {
+            for (int d = 0; d < g_dim; d++)
+            {
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
+                  d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
+                  d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               }
            }
         }
+
+         dof2quad_array.Append(d2q_new);
      }

-      dof2quad_array.Append(d2q_new);
   }
 }

@@ -1043,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)
@@ -1063,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
@@ -589,7 +589,7 @@ void H1_TriangleElement::CalcHessian(const IntegrationPoint &ip,
   Vector   shape_x(p + 1),   shape_y(p + 1),   shape_l(p + 1);
   Vector  dshape_x(p + 1),  dshape_y(p + 1),  dshape_l(p + 1);
   Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
-   DenseMatrix ddu(dof, dim);
+   DenseMatrix ddu(dof, (dim*(dim+1))/2);
 #endif

   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
@@ -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,
@@ -282,14 +282,7 @@ int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs_) const
 void FiniteElementSpace::AdjustVDofs(Array<int> &vdofs)
 {
   int n = vdofs.Size(), *vdof = vdofs;
-   for (int i = 0; i < n; i++)
-   {
-      int j;
-      if ((j = vdof[i]) < 0)
-      {
-         vdof[i] = -1-j;
-      }
-   }
+   for (int i = 0; i < n; i++) { vdof[i] = UnsignIndex(vdof[i]); }
 }

 void FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs,
@@ -483,13 +476,14 @@ void FiniteElementSpace::ReorderElementToDofTable()
   for (int k = 0, dof_counter = 0; k < nnz; k++)
   {
      const int sdof = J[k]; // signed dof
-      const int dof = (sdof < 0) ? -1-sdof : sdof;
+      const int dof = UnsignIndex(sdof);
      int new_dof = dof_marker[dof];
      if (new_dof < 0)
      {
         dof_marker[dof] = new_dof = dof_counter++;
      }
-      J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
+      // Preserve the sign of sdof
+      J[k] = (sdof < 0) ? FlipIndexSign(new_dof) : new_dof;
   }
 }

@@ -547,7 +541,7 @@ void MarkDofs(const Array<int> &dofs, Array<int> &mark_array)
 {
   for (auto d : dofs)
   {
-      mark_array[d >= 0 ? d : -1 - d] = -1;
+      mark_array[UnsignIndex(d)] = -1;
   }
 }

@@ -931,7 +925,7 @@ void FiniteElementSpace::AddDependencies(
            if (std::abs(coef) > 1e-12)
            {
               const int mdof = master_dofs[j];
-               if (mdof != sdof && mdof != (-1-sdof))
+               if (mdof != sdof && mdof != FlipIndexSign(sdof))
               {
                  deps.Add(sdof, mdof, coef);
               }
@@ -1024,7 +1018,7 @@ int FiniteElementSpace::GetDegenerateFaceDofs(int index, Array<int> &dofs,
   // FiniteElementSpace::AddDependencies.

   Array<int> edof;
-   int order = GetEdgeDofs(-1 - index, edof, variant);
+   int order = GetEdgeDofs(FlipIndexSign(index), edof, variant);

   int nv = fec->DofForGeometry(Geometry::POINT);
   int ne = fec->DofForGeometry(Geometry::SEGMENT);
@@ -1710,8 +1704,8 @@ SparseMatrix *FiniteElementSpace::RefinementMatrix_main(

         for (int i = 0; i < fine_ldof; i++)
         {
-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);

            if (!mark[m])
            {
@@ -1772,7 +1766,7 @@ SparseMatrix *FiniteElementSpace::VariableOrderRefinementMatrix(
         for (int i = 0; i < fine_ldof; i++)
         {
            const int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int m = UnsignIndex(r);

            if (!mark[m])
            {
@@ -2482,8 +2476,8 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
         {
            if (!std::isfinite(lR(i, 0))) { continue; }

-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);

            if (is_dg || !mark[m])
            {
@@ -3201,7 +3195,7 @@ void FiniteElementSpace::CalcEdgeFaceVarOrders(
            else
            {
               // degenerate face (i.e., edge-face constraint)
-               slave_orders |= edge_orders[-1 - slave.index];
+               slave_orders |= edge_orders[FlipIndexSign(slave.index)];
            }
         }

@@ -1150,7 +1150,7 @@ public:

   /// Helper to return the DOF associated with a sign encoded DOF
   static inline int DecodeDof(int dof)
-   { return (dof >= 0) ? dof : (-1 - dof); }
+   { return UnsignIndex(dof); }

   /// Helper to determine the DOF and sign of a sign encoded DOF
   static inline int DecodeDof(int dof, real_t& sign)
@@ -158,15 +158,16 @@ void LORBase::ConstructLocalDofPermutation(Array<int> &perm_) const
               int i;
               i = dofmap_lor[off_lor + i1 + i2*2];
               int s1 = i < 0 ? -1 : 1;
-               int idof_lor = vdof_lor[absdof(i)];
+               int idof_lor = vdof_lor[UnsignIndex(i)];
               i = dofmap_ho[off_ho + i1*n1 + i2*n2];
               int s2 = i < 0 ? -1 : 1;
-               int idof_ho = vdof_ho[absdof(i)];
+               int idof_ho = vdof_ho[UnsignIndex(i)];
               int s3 = idof_lor < 0 ? -1 : 1;
               int s4 = idof_ho < 0 ? -1 : 1;
               int s = s1*s2*s3*s4;
-               i = absdof(idof_ho);
-               perm_[absdof(idof_lor)] = s < 0 ? -1-absdof(i) : absdof(i);
+               i = UnsignIndex(idof_ho);
+               perm_[UnsignIndex(idof_lor)] = s < 0 ? -1-UnsignIndex(i) :
+                                              UnsignIndex(i);
            }
         }
      };
@@ -232,7 +233,7 @@ void LORBase::ConstructDofPermutation() const
         int j = l_perm[i];
         int s = j < 0 ? -1 : 1;
         int t_i = pfes_lor->GetLocalTDofNumber(i);
-         int t_j = pfes_ho->GetLocalTDofNumber(absdof(j));
+         int t_j = pfes_ho->GetLocalTDofNumber(UnsignIndex(j));
         // Either t_i and t_j both -1, or both non-negative
         if ((t_i < 0 && t_j >=0) || (t_j < 0 && t_i >= 0))
         {
@@ -57,8 +57,6 @@ private:
   /// values (after temporarily changing them for LOR assembly).
   void ResetIntegrationRules(GetIntegratorsFn get_integrators);

-   static inline int absdof(int i) { return i < 0 ? -1-i : i; }
-
 protected:
   enum FESpaceType { H1, ND, RT, L2, INVALID };

@@ -424,7 +424,7 @@ void ParFiniteElementSpace::GetGroupComm(
            {
               if (ind[l] < 0)
               {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                  if (g_ldof_sign)
                  {
                     (*g_ldof_sign)[dofs[l]] = -1;
@@ -462,7 +462,7 @@ void ParFiniteElementSpace::GetGroupComm(
            {
               if (ind[l] < 0)
               {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                  if (g_ldof_sign)
                  {
                     (*g_ldof_sign)[dofs[l]] = -1;
@@ -500,7 +500,7 @@ void ParFiniteElementSpace::GetGroupComm(
            {
               if (ind[l] < 0)
               {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                  if (g_ldof_sign)
                  {
                     (*g_ldof_sign)[dofs[l]] = -1;
@@ -538,16 +538,16 @@ void ParFiniteElementSpace::ApplyLDofSigns(Array<int> &dofs) const
   {
      if (dofs[i] < 0)
      {
-         if (ldof_sign[-1-dofs[i]] < 0)
+         if (ldof_sign[FlipIndexSign(dofs[i])] < 0)
         {
-            dofs[i] = -1-dofs[i];
+            dofs[i] = FlipIndexSign(dofs[i]);
         }
      }
      else
      {
         if (ldof_sign[dofs[i]] < 0)
         {
-            dofs[i] = -1-dofs[i];
+            dofs[i] = FlipIndexSign(dofs[i]);
         }
      }
   }
@@ -699,7 +699,8 @@ void ParFiniteElementSpace::GetSharedEdgeDofs(
      for (int i = 0; i < dofs.Size(); i++)
      {
         const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = di >= 0 ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
      }
   }
 }
@@ -723,7 +724,8 @@ void ParFiniteElementSpace::GetSharedTriangleDofs(
      for (int i = 0; i < dofs.Size(); i++)
      {
         const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = di >= 0 ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
      }
   }
 }
@@ -747,7 +749,8 @@ void ParFiniteElementSpace::GetSharedQuadrilateralDofs(
      for (int i = 0; i < dofs.Size(); i++)
      {
         const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = (di >= 0) ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
      }
   }
 }
@@ -1487,7 +1490,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
         GetElementVDofs(my_elems[i], ldofs);
         for (int j = 0; j < ldofs.Size(); j++)
         {
-            int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
+            int ldof = UnsignIndex(ldofs[j]);

            if (ldof_marker[ldof] != fn)
            {
@@ -1548,7 +1551,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
         GetElementVDofs(my_elems[i], ldofs);
         for (int j = 0; j < ldofs.Size(); j++)
         {
-            int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
+            int ldof = UnsignIndex(ldofs[j]);

            if (ldof_marker[ldof] != fn)
            {
@@ -1573,14 +1576,15 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()

      for (int i = 0; i < num_ldofs; i++)
      {
-         int ldof = (ldofs_fn[i] >= 0 ? ldofs_fn[i] : -1-ldofs_fn[i]);
+         int ldof = UnsignIndex(ldofs_fn[i]);
         ldof_marker[ldof] = i;
      }

      for ( ; j < j_end; j++)
      {
-         int ldof = (send_J[j] >= 0 ? send_J[j] : -1-send_J[j]);
-         send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] : -1-ldof_marker[ldof]);
+         const int ldof = UnsignIndex(send_J[j]);
+         send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] :
+                      FlipIndexSign(ldof_marker[ldof]));
      }
   }

@@ -1672,12 +1676,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
   {
      for (int j_end = face_nbr_ldof.GetI()[fn+1]; j < j_end; j++)
      {
-         int ldof = face_nbr_ldof.GetJ()[j];
-         if (ldof < 0)
-         {
-            ldof = -1-ldof;
-         }
-
+         const int ldof = UnsignIndex(face_nbr_ldof.GetJ()[j]);
         face_nbr_glob_dof_map[j] = dof_face_nbr_offsets[fn] + ldof;
      }
   }
@@ -1721,7 +1720,7 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
   MFEM_ASSERT(Nonconforming() && i >= pmesh->GetNumFaces(), "");
   int el1, el2, inf1, inf2;
   pmesh->GetFaceElements(i, &el1, &el2);
-   el2 = -1 - el2;
+   el2 = FlipIndexSign(el2);
   pmesh->GetFaceInfos(i, &inf1, &inf2);
   MFEM_ASSERT(0 <= el2 && el2 < face_nbr_element_dof.Size(), "");
   const int nd = face_nbr_element_dof.RowSize(el2);
@@ -1737,7 +1736,8 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
   for (int j = 0; j < vdofs.Size(); j++)
   {
      const int ldof = vdofs[j];
-      vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] : -1-vol_vdofs[-1-ldof];
+      vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] :
+                 FlipIndexSign(vol_vdofs[FlipIndexSign(ldof)]);
   }
 }

@@ -2061,8 +2061,8 @@ void ParFiniteElementSpace::GetGhostFaceDofs(const MeshId &face_id,

         for (int j = 0; j < ne; j++)
         {
-            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
-                             /*         */ : (-1 - (first + (-1 - ind[j])));
+            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j]) :
+                             FlipIndexSign(first + FlipIndexSign(ind[j]));
         }
      }
      else
@@ -2072,8 +2072,8 @@ void ParFiniteElementSpace::GetGhostFaceDofs(const MeshId &face_id,
         const int *ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[i]);
         for (int j = 0; j < ne; j++)
         {
-            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
-                             /*         */ : (-1 - (first + (-1 - ind[j])));
+            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j]) :
+                             FlipIndexSign(first + FlipIndexSign(ind[j]));
         }
      }
   }
@@ -2866,7 +2866,7 @@ void NeighborRowMessage::Encode(int rank)

            if (ind && (edof = ind[edof]) < 0)
            {
-               edof = -1 - edof;
+               edof = FlipIndexSign(edof);
               s = -1;
            }

@@ -3067,10 +3067,10 @@ void NeighborRowMessage::Decode(int rank)

         // If edof arrived with a negative index, flip it, and the scaling.
         real_t s = (edof < 0) ? -1.0 : 1.0;
-         edof = (edof < 0) ? -1 - edof : edof;
+         edof = UnsignIndex(edof);
         if (ind && (edof = ind[edof]) < 0)
         {
-            edof = -1 - edof;
+            edof = FlipIndexSign(edof);
            s *= -1.0;
         }

@@ -3121,10 +3121,10 @@ void NeighborRowMessage::Decode(int rank)

            // If edof arrived with a negative index, flip it, and the scaling.
            s = (edof < 0) ? -1.0 : 1.0;
-            edof = (edof < 0) ? -1 - edof : edof;
+            edof = UnsignIndex(edof);
            if (ind && (edof = ind[edof]) < 0)
            {
-               edof = -1 - edof;
+               edof = FlipIndexSign(edof);
               s *= -1.0;
            }

@@ -4405,12 +4405,9 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
         {
            for (int j = 0; j < dofs.Size(); j++)
            {
-               int row = DofToVDof(dofs[j], vd);
-               if (row < 0) { row = -1 - row; }
-
-               int col = DofToVDof(old_dofs[j], vd, old_ndofs);
-               if (col < 0) { col = -1 - col; }
-
+               const int row = UnsignIndex(DofToVDof(dofs[j], vd));
+               const int col = UnsignIndex(DofToVDof(old_dofs[j], vd,
+                                                     old_ndofs));
               i_diag[row] = col;
            }
         }
@@ -4435,9 +4432,7 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
      {
         for (int j = 0; j < dofs.Size(); j++)
         {
-            int row = DofToVDof(dofs[j], vd);
-            if (row < 0) { row = -1 - row; }
-
+            const int row = UnsignIndex(DofToVDof(dofs[j], vd));
            if (i_diag[row] == i_diag[row+1]) // diag row empty?
            {
               i_offd[row] = old_dofs[j + vd * dofs.Size()];
@@ -4546,9 +4541,9 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
   {
      const Embedding &emb = dtrans.embeddings[k];

-      int fine_rank = old_ranks[k];
-      int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
-                        : old_pncmesh->ElementRank(emb.parent);
+      const int fine_rank = old_ranks[k];
+      const int coarse_rank = (emb.parent < 0) ? FlipIndexSign(emb.parent)
+                              : old_pncmesh->ElementRank(emb.parent);

      if (coarse_rank != MyRank && fine_rank == MyRank)
      {
@@ -4636,8 +4631,8 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
            {
               if (!std::isfinite(lR(i, 0))) { continue; }

-               int r = DofToVDof(dofs[i], vd);
-               int m = (r >= 0) ? r : (-1 - r);
+               const int r = DofToVDof(dofs[i], vd);
+               const int m = UnsignIndex(r);

               if (is_dg || !mark[m])
               {
@@ -4686,8 +4681,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
            {
               if (!std::isfinite(lR(i, 0))) { continue; }

-               int r = DofToVDof(dofs[i], vd);
-               int m = (r >= 0) ? r : (-1 - r);
+               const int m = UnsignIndex(DofToVDof(dofs[i], vd));

               if (is_dg || !mark[m])
               {
@@ -271,10 +271,7 @@ inline void QuadratureFunction::GetValues(
   const int s_offset = qspace->Offset(idx);
   const int sl_size = qspace->Offset(idx + 1) - s_offset;
   // Make the values matrix memory an alias of the quadrature function memory
-   Memory<real_t> &values_mem = values.GetMemory();
-   values_mem.Delete();
-   values_mem.MakeAlias(GetMemory(), vdim*s_offset, vdim*sl_size);
-   values.SetSize(vdim, sl_size);
+   values.MakeRef(GetMemory(), vdim*s_offset, vdim, sl_size);
 }

 inline void QuadratureFunction::GetValues(
@@ -844,8 +844,6 @@ void ConformingFaceRestriction::ComputeGatherIndices(
   gather_offsets[0] = 0;
 }

-static inline int absdof(int i) { return i < 0 ? -1-i : i; }
-
 void ConformingFaceRestriction::SetFaceDofsScatterIndices(
   const Mesh::FaceInformation &face,
   const int face_index,
@@ -868,9 +866,9 @@ void ConformingFaceRestriction::SetFaceDofsScatterIndices(
   {
      const int lex_volume_dof = face_map[face_dof];
      const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof]; // signed
-      const int volume_dof = absdof(s_volume_dof);
+      const int volume_dof = UnsignIndex(s_volume_dof);
      const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
-      const int global_dof = absdof(s_global_dof);
+      const int global_dof = UnsignIndex(s_global_dof);
      const int restriction_dof = face_dofs*face_index + face_dof;
      scatter_indices[restriction_dof] = s_global_dof;
      ++gather_offsets[global_dof + 1];
@@ -897,10 +895,10 @@ void ConformingFaceRestriction::SetFaceDofsGatherIndices(
   {
      const int lex_volume_dof = face_map[face_dof];
      const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof];
-      const int volume_dof = absdof(s_volume_dof);
+      const int volume_dof = UnsignIndex(s_volume_dof);
      const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
      const int sgn = (s_global_dof >= 0) ? 1 : -1;
-      const int global_dof = absdof(s_global_dof);
+      const int global_dof = UnsignIndex(s_global_dof);
      const int restriction_dof = face_dofs*face_index + face_dof;
      const int s_restriction_dof = (sgn >= 0) ? restriction_dof : -1 -
                                    restriction_dof;
@@ -4102,8 +4102,11 @@ void TMOP_Integrator::GetSurfaceFittingErrors(const Vector &d_loc,
 #ifdef MFEM_USE_MPI
      // Don't count the overlapping DOFs in parallel.
      // The pfes might be ordered byVDIM, while the loop goes consecutively.
-      const int dof_i = pfes->DofToVDof(i, 0);
-      if (parallel && pfes->GetLocalTDofNumber(dof_i) < 0) { continue; }
+      if (parallel)
+      {
+         const int dof_i = pfes->DofToVDof(i, 0);
+         if (pfes->GetLocalTDofNumber(dof_i) < 0) { continue; }
+      }
 #endif

      dof_cnt++;
@@ -114,10 +114,22 @@ public:
   Array<T> &operator=(const Array<T> &src) { src.Copy(*this); return *this; }

   /// Move assignment operator
+   /** If *this is a non-owning view (e.g., from MakeRef()), the data is copied
+       so that the base is also modified. */
   Array<T> &operator=(Array<T> &&src)
   {
      if (this == &src) { return *this; }
-      Swap(src);  // Swap does not use move assignment!
+      // If *this is a non-owning view (alias), and its capacity is sufficient
+      // to contain src, then copy into *this so that the alias's base memory is
+      // modified.
+      if (!OwnsData() && Capacity() >= src.Size())
+      {
+         *this = src; // Copy assignment.
+      }
+      else
+      {
+         Swap(src); // Swap the pointers only.
+      }
      src.DeleteAll();
      return *this;
   }
@@ -251,6 +263,9 @@ public:
   /// Make this Array a reference to 'master'.
   inline void MakeRef(const Array &master);

+   /// Make this Array a reference to the given sub-Memory of @a base.
+   inline void MakeRef(Memory<T> &base, int offset, int size_);
+
   /// Reset the Array to use the given external Memory @a mem and size @a s.
   /** If @a own_mem is false, the Array will not own any of the pointers of
       @a mem.
@@ -1073,6 +1088,14 @@ inline void Array<T>::MakeRef(const Array &master)
   data.MakeAlias(master.GetMemory(), 0, size);
 }

+template <class T>
+inline void Array<T>::MakeRef(Memory<T> &base, int offset, int size_)
+{
+   data.Delete();
+   size = size_;
+   data.MakeAlias(base, offset, size_);
+}
+
 template <class T>
 inline void Array<T>::NewMemoryAndSize(
   const Memory<T> &mem, int s, bool own_mem)
@@ -113,6 +113,10 @@ void SetGlobalMPI_Comm(MPI_Comm comm);
 /// to suppress the warning.
 const char* GetEnv(const char* name);

+/// Signed indices i -> -1 - i are used as a convention to encode orientation.
+inline MFEM_HOST_DEVICE int FlipIndexSign(int i) { return -1 - i; }
+inline MFEM_HOST_DEVICE int UnsignIndex(int i) { return i < 0 ? -1 - i : i; }
+
 } // namespace mfem

 #endif
@@ -759,7 +759,9 @@ private:
   {
      switch (mt)
      {
-         case MT::HOST_DEBUG: return new MmuHostMemorySpace();
+         case MT::HOST_DEBUG:
+            if (GetEnv("MFEM_MMU_STD")) { return new StdHostMemorySpace(); }
+            return new MmuHostMemorySpace();
 #ifdef MFEM_USE_UMPIRE
         case MT::HOST_UMPIRE:
            return new UmpireHostMemorySpace(
@@ -788,7 +790,9 @@ private:
         case MT::DEVICE_UMPIRE: return new NoDeviceMemorySpace();
         case MT::DEVICE_UMPIRE_2: return new NoDeviceMemorySpace();
 #endif
-         case MT::DEVICE_DEBUG: return new MmuDeviceMemorySpace();
+         case MT::DEVICE_DEBUG:
+            if (GetEnv("MFEM_MMU_STD")) { return new StdDeviceMemorySpace(); }
+            return new MmuDeviceMemorySpace();
         case MT::DEVICE:
         {
 #if defined(MFEM_USE_CUDA)
@@ -82,6 +82,15 @@ public:
      }
   }

+   /// Make the DenseMatrix to reference the given sub-Memory of @a base.
+   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
+       not delete the @a base Memory. */
+   void MakeRef(Memory<real_t> &base, int offset, int h, int w)
+   {
+      data.MakeRef(base, offset, h*w);
+      height = h; width = w;
+   }
+
   /// Change the data array and the size of the DenseMatrix.
   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
       not delete the data array @a d. */
@@ -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)

@@ -493,8 +493,7 @@ void Mesh::GetBdrElementTransformation(int i,
         {
            for (int j = 0; j < n; j++)
            {
-               int idx = vdofs[n*k+j];
-               pm(k,j) = nodes((idx<0)? -1-idx:idx);
+               pm(k,j) = nodes(UnsignIndex(vdofs[n*k+j]));
            }
         }
         ElTr->SetFE(bdr_el);
@@ -1356,7 +1355,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
               face.element[0].conformity = ElementConformity::Coincident;
               face.element[1].conformity = ElementConformity::Coincident;
               face.element[1].location = ElementLocation::FaceNbr;
-               face.element[1].index = -1 - e2;
+               face.element[1].index = FlipIndexSign(e2);
               face.element[1].orientation = inf2%64;
            }
         }
@@ -1379,7 +1378,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
               face.element[1].location = ElementLocation::FaceNbr;
               face.element[0].conformity = ElementConformity::Coincident;
               face.element[1].conformity = ElementConformity::Superset;
-               face.element[1].index = -1 - e2;
+               face.element[1].index = FlipIndexSign(e2);
               face.element[1].orientation = inf2%64;
            }
            face.point_matrix = nc_faces_info[ncface].PointMatrix;
@@ -1405,7 +1404,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
         face.element[1].location = ElementLocation::FaceNbr;
         face.element[0].conformity = ElementConformity::Superset;
         face.element[1].conformity = ElementConformity::Coincident;
-         face.element[1].index = -1 - e2;
+         face.element[1].index = FlipIndexSign(e2);
         face.element[1].orientation = inf2%64;
         face.point_matrix = nc_faces_info[ncface].PointMatrix;
      }
@@ -1438,7 +1437,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
         break;
      case FaceInfoTag::SharedConforming:
         res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
         res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
         res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
         break;
@@ -1448,7 +1447,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
         break;
      case FaceInfoTag::SharedSlaveNonconforming:
         res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
         res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
         res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
         break;
@@ -1456,7 +1455,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
         break;
      case FaceInfoTag::GhostSlave:
         res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
         res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
         res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
         break;
@@ -6514,7 +6513,7 @@ void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_ukv)
         input >> edge_to_ukv[j] >> v[0] >> v[1];
         if (v[0] > v[1])
         {
-            edge_to_ukv[j] = -1 - edge_to_ukv[j];
+            edge_to_ukv[j] = FlipIndexSign(edge_to_ukv[j]);
         }
      }
   }
@@ -6551,9 +6550,6 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
   const int NP = NumOfElements; // number of patches
   const int NPKV = NP * dim;    // number of patch knotvectors
   constexpr int notset = -9999999;
-   // Sign convention
-   auto flipSign = [](int i) { return -1 - i; };
-   auto unSign = [](int i) { return (i < 0) ? -1 - i : i; };
   // Local edge index -> dimension convention
   auto edge_to_dim = [](int i) { return (i < 8) ? ((i & 1) ? 1 : 0) : 2; };

@@ -6569,7 +6565,7 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
      {
         GetElementVertices(i, v);
         // Sign is based on the edge's vertex indices
-         edge_to_ukv[i] = (v[1] > v[0]) ? i : flipSign(i);
+         edge_to_ukv[i] = (v[1] > v[0]) ? i : FlipIndexSign(i);
         ukv_to_rpkv[i] = i;
      }
      return;
@@ -6619,14 +6615,14 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
         // We've set this edge already - link this index to it
         if (edge_to_pkv[edge] != notset)
         {
-            const int pkv_other = unSign(edge_to_pkv[edge]);
+            const int pkv_other = UnsignIndex(edge_to_pkv[edge]);
            unite(pkv, pkv_other);
         }
         else
         {
            GetEdgeVertices(edge, v);
            // Sign is based on the edge's vertex indices
-            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : flipSign(pkv);
+            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : FlipIndexSign(pkv);
         }
      }
   }
@@ -6653,10 +6649,10 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
   edge_to_ukv.SetSize(NumOfEdges);
   for (int i = 0; i < NumOfEdges; i++)
   {
-      const int pkv = unSign(edge_to_pkv[i]);
+      const int pkv = UnsignIndex(edge_to_pkv[i]);
      const int rpkv = pkv_to_rpkv[pkv];
      const int ukv = rpkv_to_ukv[rpkv];
-      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? flipSign(ukv) : ukv;
+      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? FlipIndexSign(ukv) : ukv;
   }

   CorrectPatchTopoOrientations(edge_to_ukv);
@@ -6667,9 +6663,6 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
   const int dim = Dimension(); // Topological (not physical) dimension
   if (dim == 1) { return; }

-   // Sign convention
-   auto flipSign = [](int i) { return -1 - i; };
-
   const Table *face2elem = GetFaceToElementTable();
   Array<int> pfaces, orient;
   Array<int> fe, feo;
@@ -6688,7 +6681,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
         for (auto e : fe)
         {
            const int skv = edge_to_ukv[e];
-            if (skv == kv || flipSign(skv) == kv) { hasKV = true; }
+            if (skv == kv || FlipIndexSign(skv) == kv) { hasKV = true; }
         }
         if (hasKV)
         {
@@ -6718,7 +6711,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
      };
   }

-   Array<int> ukvs((dim==2) ? 4 : 12);
+   Array<int> ukvs((dim == 2) ? 4 : 12);
   Array<int> pe, oe;
   bool initKV = false;

@@ -6732,7 +6725,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
      for (int i = 0; i < pe.Size(); i++)
      {
         ukvs[i] = edge_to_ukv[pe[i]];
-         ukvs[i] = (oe[i] < 0) ? flipSign(ukvs[i]) : ukvs[i];
+         ukvs[i] = (oe[i] < 0) ? FlipIndexSign(ukvs[i]) : ukvs[i];
      }

      // Find the direction with this kv.
@@ -6740,12 +6733,19 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
      for (int d=0; d<dim; ++d) // Loop over directions.
      {
         const int skv = edge_to_ukv[pe[dir_edges[d][0]]];
-         if (skv == kv || flipSign(skv) == kv)
+         if (skv == kv || FlipIndexSign(skv) == kv)
         {
-            thisDir = d;
+            for (auto e : dir_edges[d])
+               if (!edgeSet[pe[e]])
+               {
+                  thisDir = d;
+               }
         }
      }
-      MFEM_VERIFY(thisDir >= 0, "");
+      if (thisDir == -1)
+      {
+         return false;
+      }

      // For this direction, find any edge already set. If no edge is set, we
      // arbitrarily take the first.
@@ -6777,12 +6777,12 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
         }

         const int edge = pe[i];
-         if ((dim == 2 && ukvs[i] != flipSign(ukvs[ref_edge0])) ||
-             (dim == 3 && ukvs[i] == flipSign(ukvs[ref_edge0])))
+         if ((dim == 2 && ukvs[i] != FlipIndexSign(ukvs[ref_edge0])) ||
+             (dim == 3 && ukvs[i] == FlipIndexSign(ukvs[ref_edge0])))
         {
            // Flip the sign of this edge
-            MFEM_VERIFY(!edgeSet[edge], "");
-            edge_to_ukv[edge] = flipSign(edge_to_ukv[edge]);
+            MFEM_ASSERT(!edgeSet[edge], "");
+            edge_to_ukv[edge] = FlipIndexSign(edge_to_ukv[edge]);
         }

         edgeSet[edge] = true;
@@ -6827,10 +6827,11 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
      int unsetDim = -1;
      for (int d=0; d<dim; ++d) // Loop over dimensions.
      {
-         if (!edgeSet[pe[dir_edges[d][0]]])
-         {
-            unsetDim = d;
-         }
+         for (auto e : dir_edges[d])
+            if (!edgeSet[pe[e]])
+            {
+               unsetDim = d;
+            }
      }

      if (unsetDim == -1)
@@ -6839,9 +6840,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
         continue;
      }

-      const int kv_signed = edge_to_ukv[pe[dir_edges[unsetDim][0]]];
-      const int kv = kv_signed < 0 ? flipSign(kv_signed) : kv_signed;
-      MFEM_VERIFY(!edgeSet[pe[dir_edges[unsetDim][0]]], "");
+      const int kv = UnsignIndex(edge_to_ukv[pe[dir_edges[unsetDim][0]]]);

      initKV = false;

@@ -6891,6 +6890,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
      }
   }

+#ifdef MFEM_DEBUG
   bool allSet = true;
   for (auto eset : edgeSet)
   {
@@ -6899,7 +6899,8 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
         allSet = false;
      }
   }
-   MFEM_VERIFY(allSet && unset.size() == 0, "Some edge is not set");
+   MFEM_ASSERT(allSet && unset.size() == 0, "Some edge is not set");
+#endif

   delete face2elem;
 }
@@ -6941,7 +6942,7 @@ void Mesh::LoadNonconformingPatchTopo(std::istream &input,

      if (v[0] > v[1])
      {
-         ukv = -1 - ukv;
+         ukv = FlipIndexSign(ukv);
      }
      edge_to_ukv[j] = ukv;
   }
@@ -12438,11 +12439,7 @@ void Mesh::PrintTopoEdges(std::ostream &os, const Array<int> &e_to_k,
   for (int i = 0; i < NumOfEdges; i++)
   {
      edge_vertex->GetRow(i, vert);
-      int ki = e_to_k[i];
-      if (ki < 0)
-      {
-         ki = -1 - ki;
-      }
+      const int ki = UnsignIndex(e_to_k[i]);

      if (vmap)
      {
@@ -15751,9 +15748,18 @@ Mesh PartitionMPI(int dim, int mpi_cnt, int elem_per_mpi, bool print,
 {
   MFEM_VERIFY(dim > 1, "Not implemented for 1D meshes.");

-   auto factor = [&](int N)
+   // Closest int divisor to the cubit root, going down.
+   auto factor3 = [](int N)
   {
-      for (int i = static_cast<int>(sqrt(N)); i > 0; i--)
+      for (int i = static_cast<int>(round(cbrt(N))); i > 0; i--)
+      {  if (N % i == 0) { return i; } }
+      return 1;
+   };
+
+   // Closest int divisor to the square root, going down.
+   auto factor2 = [](int N)
+   {
+      for (int i = static_cast<int>(round(sqrt(N))); i > 0; i--)
      { if (N % i == 0) { return i; } }
      return 1;
   };
@@ -15777,22 +15783,22 @@ Mesh PartitionMPI(int dim, int mpi_cnt, int elem_per_mpi, bool print,
   int el0_x, el0_y, el0_z;
   if (dim == 2)
   {
-      mpi_x = factor(mpi_cnt);
+      mpi_x = factor2(mpi_cnt);
      mpi_y = mpi_cnt / mpi_x;

      // Switch order for better balance.
-      el0_y = factor(el0);
+      el0_y = factor2(el0);
      el0_x = el0 / el0_y;
   }
   else
   {
-      mpi_x = factor(mpi_cnt);
-      mpi_y = factor(mpi_cnt / mpi_x);
+      mpi_x = factor3(mpi_cnt);
+      mpi_y = factor2(mpi_cnt / mpi_x);
      mpi_z = mpi_cnt / mpi_x / mpi_y;

      // Switch order for better balance.
-      el0_z = factor(el0);
-      el0_y = factor(el0 / el0_z);
+      el0_z = factor3(el0);
+      el0_y = factor2(el0 / el0_z);
      el0_x = el0 / el0_y / el0_z;
   }

@@ -143,7 +143,7 @@ int ThresholdRefiner::ApplyImpl(Mesh &mesh)
   if (num_marked_elements == 0LL) { return STOP; }

   mesh.GeneralRefinement(marked_elements, non_conforming, nc_limit);
-   return CONTINUE + REFINED;
+   return static_cast<int>(CONTINUE) + static_cast<int>(REFINED);
 }

 void ThresholdRefiner::Reset()
@@ -162,7 +162,7 @@ int ThresholdDerefiner::ApplyImpl(Mesh &mesh)
   const Vector &local_err = estimator.GetLocalErrors();
   bool derefs = mesh.DerefineByError(local_err, threshold, nc_limit, op);

-   return derefs ? CONTINUE + DEREFINED : NONE;
+   return derefs ? static_cast<int>(CONTINUE) + static_cast<int>(DEREFINED) : NONE;
 }


@@ -290,7 +290,7 @@ int CoefficientRefiner::PreprocessMesh(Mesh &mesh, int max_it)
   }
   delete l2fes;
   delete gf;
-   return CONTINUE + REFINED;
+   return static_cast<int>(CONTINUE) + static_cast<int>(REFINED);

 }

@@ -310,7 +310,7 @@ int Rebalancer::ApplyImpl(Mesh &mesh)
   if (pmesh && pmesh->Nonconforming())
   {
      pmesh->Rebalance();
-      return CONTINUE + REBALANCED;
+      return static_cast<int>(CONTINUE) + static_cast<int>(REBALANCED);
   }
 #endif
   return NONE;
@@ -3542,7 +3542,7 @@ void NCMesh::TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,

            // create a slave face record with a degenerate point matrix
            face_list.slaves.Append(
-               Slave(-1 - enode.edge_index,
+               Slave(FlipIndexSign(enode.edge_index),
                     eid[0].element, eid[0].local, Geometry::SQUARE));
            Slave &sl = face_list.slaves.Last();

@@ -3589,7 +3589,7 @@ void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
         // non-slave edge is really a (face-)slave itself.
         const MeshId &eid = *eid_and_type.id;
         face_list.slaves.Append(
-            Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
+            Slave(FlipIndexSign(eid.index), eid.element, eid.local, Geometry::TRIANGLE));

         int v0index = nodes[vn0].vert_index;
         int v1index = nodes[vn1].vert_index;
@@ -93,7 +93,7 @@ void NCNURBSExtension::GetMasterEdgeEntities(
      }
      else
      {
-         const int auxEdge = -1 - edge_i;
+         const int auxEdge = FlipIndexSign(edge_i);
         GetAuxEdgeVertices(auxEdge, sverts);
      }

@@ -159,7 +159,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
            {
               if (edge < 0)
               {
-                  sideAuxEdges[s].Append(-1 - edge);
+                  sideAuxEdges[s].Append(FlipIndexSign(edge));
               }
               else
               {
@@ -456,7 +456,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                                       == afverts[j], "");
                        }

-                        ori_f2 = -1 - ori_f2;
+                        ori_f2 = FlipIndexSign(ori_f2);
                     }
                     else
                     {
@@ -468,7 +468,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                     }

                     facePairs.emplace_back(FacePairInfo{fverts[vMinID], f,
-                     SlaveFaceInfo{-1 - afid, ori_f2,
+                     SlaveFaceInfo{FlipIndexSign(afid), ori_f2,
                        {fki(vMinID,0), fki(vMinID,1)},
                        {
                           fki((vMinID + 2) % 4,0) - fki(vMinID,0),
@@ -509,7 +509,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                     auxFaces.push_back(auxFace);

                     facePairs.emplace_back(FacePairInfo{fverts[vMinID], f,
-                     SlaveFaceInfo{-1 - auxFaceId, ori_f,
+                     SlaveFaceInfo{FlipIndexSign(auxFaceId), ori_f,
                        {fki(vMinID,0), fki(vMinID,1)},
                        {
                           fki((vMinID + 2) % 4,0) - fki(vMinID,0),
@@ -622,7 +622,7 @@ void NCNURBSExtension::GetAuxFaceEdges(int auxFace, Array<int> &edges) const
      }
      else // Auxiliary edge
      {
-         edges[i] = -1 - auxv2e.at(edge_v);
+         edges[i] = FlipIndexSign(auxv2e.at(edge_v));
      }
   }
 }
@@ -633,7 +633,7 @@ int OffsetHelper(int i, int j, const Array<int> &a, const Array<int> &b)
 {
   if (i < 0)
   {
-      return b[-1 - i + j];
+      return b[FlipIndexSign(i) + j];
   }
   else if (i + j < a.Size())
   {
@@ -679,7 +679,7 @@ void NCNURBSExtension::GetMasterEdgeDofs(bool dof, int me,
      }
      else // Auxiliary edge
      {
-         GetAuxEdgeVertices(-1 - slaveId, svert);
+         GetAuxEdgeVertices(FlipIndexSign(slaveId), svert);
      }

      bool reverse = false;
@@ -872,7 +872,7 @@ void ReorderArray2D(int i0, int j0, const Array2D<int> &a,
 // Set a quadrilateral vertex index permutation for a given orientation.
 void GetVertexOrdering(int ori, std::array<int, 4> &perm)
 {
-   const int oriAbs = ori < 0 ? -1 - ori : ori;
+   const int oriAbs = UnsignIndex(ori);

   for (int i=0; i<4; ++i)
   {
@@ -1094,7 +1094,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
      if (slaveId < 0)
      {
         // Auxiliary face
-         const int auxFace = -1 - slaveId;
+         const int auxFace = FlipIndexSign(slaveId);

         // Set slave face entity dimensions.
         if (dof)
@@ -1171,7 +1171,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
            }
            else
            {
-               const int auxEdge = -1 - edge;
+               const int auxEdge = FlipIndexSign(edge);
               GetAuxEdgeVertices(auxEdge, evert);
            }
            MFEM_ASSERT(evert[0] == vstart || evert[1] == vstart, "");
@@ -1184,7 +1184,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
         // dimensions of the master face, by using ori.
         int e1 = -1, e2 = -1;
         {
-            const int aori = ori < 0 ? -1 - ori : ori;
+            const int aori = UnsignIndex(ori);
            if (aori % 2 == 0)
            {
               e1 = 0;
@@ -1416,14 +1416,15 @@ void NCNURBSExtension::ProcessVertexToKnot2D(const VertexToKnotSpan &v2k,
         {
            // Create a new auxiliary edge
            auxv2e[childPair] = auxEdges.size();
-            auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ?
-                                                parentEdge : -1 - parentEdge,
+            auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ? parentEdge :
+                                                FlipIndexSign(parentEdge),
            {childPair.first, childPair.second},
            {newParentEdge ? 0 : prevKI, ks}});
         }
      }

-      const int childEdge = childPairTopo ? v2e[childPair] : -1 - auxv2e[childPair];
+      const int childEdge = childPairTopo ? v2e[childPair] :
+                            FlipIndexSign(auxv2e[childPair]);

      // Check whether this is the final vertex in this parent edge. Note that
      // the logic for comparing (pv[0],pv[1]) to the next parents assumes the
@@ -1460,14 +1461,15 @@ void NCNURBSExtension::ProcessVertexToKnot2D(const VertexToKnotSpan &v2k,

               // -1 denotes `ne` at endpoint
               auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ?
-                                                   -1 - parentEdge : parentEdge,
+                                                   FlipIndexSign(parentEdge) :
+                                                   parentEdge,
               {finalChildPair.first, finalChildPair.second},
               {ks, -1}});
            }
         }

         const int finalChildEdge = finalChildPairTopo ? v2e[finalChildPair] :
-                                    -1 - auxv2e[finalChildPair];
+                                    FlipIndexSign(auxv2e[finalChildPair]);
         edgePairs.emplace_back(-1, -1, finalChildEdge, parentEdge);
      }

@@ -1805,7 +1807,7 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
                  auxFaces.push_back(auxFace);
                  facePairs.emplace_back(
                     FacePairInfo{cv[0], parentFace,
-                                  SlaveFaceInfo{-1 - auxv2f[childPair],
+                                  SlaveFaceInfo{FlipIndexSign(auxv2f[childPair]),
                                                0, {i0, j0}, {d0, d1}}});
               }
            }
@@ -2111,7 +2113,7 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
                     auxv2e[childPair] = auxEdges.size();
                     auxEdges.emplace_back(AuxiliaryEdge{pv0 < pv1 ?
                                                         parentEdge :
-                                                         -1 - parentEdge,
+                                                         FlipIndexSign(parentEdge),
                     {childPair.first, childPair.second},
                     {knotIndex0, knotIndex1}});
                  }
@@ -2131,7 +2133,8 @@ void NCNURBSExtension::ProcessVertexToKnot3D(

                  const EdgePairInfo ep_e((e_idx == n_d - de) ? -1 : tv,
                                          (e_idx == n_d - de) ? -1 : tvki,
-                                          -1 - auxv2e[childPair], parentEdge);
+                                          FlipIndexSign(auxv2e[childPair]),
+                                          parentEdge);

                  const bool unset = !edgePairs[edgePairOS[parentEdge] + e_idx].isSet;
                  if (unset)
@@ -2226,7 +2229,7 @@ void NCNURBSExtension::GetAuxFaceToPatchTable(Array2D<int> &auxface2patch)
               if (s < 0)
               {
                  // Auxiliary face.
-                  const int aux = -1 - s;
+                  const int aux = FlipIndexSign(s);
                  if (auxface2patch(aux, 0) >= 0)
                  {
                     if (auxface2patch(aux, 1) != -1) { consistent = false; }
@@ -2316,7 +2319,7 @@ void NCNURBSExtension::UpdateAuxiliaryKnotSpans(const Array<int> &rf)
   for (auto auxEdge : auxEdges)
   {
      const int p = auxEdge.parent;
-      const int parent = p < 0 ? -1 - p : p;
+      const int parent = UnsignIndex(p);
      const int kv = KnotInd(parent);
      for (int i=0; i<2; ++i)
      {
@@ -2382,14 +2385,8 @@ int NCNURBSExtension::AuxiliaryEdgeNE(int aux_edge)
   const int signedParentEdge = auxEdges[aux_edge].parent;
   const int ki0 = auxEdges[aux_edge].ksi[0];
   const int ki1raw = auxEdges[aux_edge].ksi[1];
-   int ki1 = ki1raw;
-   if (ki1raw == -1)
-   {
-      const bool rev = signedParentEdge < 0;
-      const int parentEdge = rev ? -1 - signedParentEdge : signedParentEdge;
-      ki1 = KnotVec(parentEdge)->GetNE();
-   }
-
+   const int ki1 = ki1raw == -1 ? KnotVec(UnsignIndex(signedParentEdge))->GetNE()
+                   : ki1raw;
   return ki1 - ki0;
 }

@@ -2403,7 +2400,7 @@ void NCNURBSExtension::SlaveEdgeToParent(int se, int parent,
   Array<int> sev(2);
   if (se < 0) // Auxiliary edge
   {
-      for (int i=0; i<2; ++i) { sev[i] = auxEdges[-1 - se].v[i]; }
+      for (int i=0; i<2; ++i) { sev[i] = auxEdges[FlipIndexSign(se)].v[i]; }
   }
   else
   {
@@ -2459,7 +2456,7 @@ void NCNURBSExtension::GetMasterEdgePieceOffsets(int mid, Array<int> &os)
      }
      else
      {
-         nes = AuxiliaryEdgeNE(-1 - s);
+         nes = AuxiliaryEdgeNE(FlipIndexSign(s));
      }

      os[i+1] = os[i] + nes;
@@ -2565,7 +2562,7 @@ int NCNURBSExtension::SetPatchFactors(int p)
               }
               else // Aux edge
               {
-                  const int aux_edge = -1 - s;
+                  const int aux_edge = FlipIndexSign(s);
                  if (auxef[aux_edge].Size() == 0)
                  {
                     auxef[aux_edge].SetSize(AuxiliaryEdgeNE(aux_edge));
@@ -2611,7 +2608,7 @@ int NCNURBSExtension::SetPatchFactors(int p)
   }

   MFEM_VERIFY(consistent, "");
-   return partialChange ? -1 - dirSet : dirSet;
+   return partialChange ? FlipIndexSign(dirSet) : dirSet;
 }

 void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
@@ -2715,7 +2712,7 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
               if (s < 0)
               {
                  // Auxiliary face.
-                  const int aux = -1 - s;
+                  const int aux = FlipIndexSign(s);
                  for (int i=0; i<2; ++i)
                  {
                     const int patch = auxface2patch(aux, i);
@@ -2763,7 +2760,7 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)

         const int dirSetSigned = SetPatchFactors(p);
         const bool partialChange = dirSetSigned < 0;
-         const int dirSet = partialChange ? -1 - dirSetSigned : dirSetSigned;
+         const int dirSet = UnsignIndex(dirSetSigned);
         const bool changed = (patchState[p] != dirSet) || partialChange;
         patchState[p] = dirSet;

@@ -2806,8 +2803,8 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
            {
               const int dirSetSigned_i = SetPatchFactors(i);
               const bool partialChange_i = dirSetSigned_i < 0;
-               const int dirSet_i = partialChange_i ? -1 - dirSetSigned_i :
-                                    dirSetSigned_i;
+               const int dirSet_i = partialChange_i ?
+                                    FlipIndexSign(dirSetSigned_i) : dirSetSigned_i;
               const bool changed_i = (patchState[i] != dirSet_i) ||
                                      partialChange_i;
               patchState[p] = dirSet_i;
@@ -3027,7 +3024,7 @@ int GetFaceOrientation(const Mesh *mesh, const int face,

   // Check whether ordering is reversed.
   const bool rev = verts[(s + 1) % 4] != fverts[1];
-   if (rev) { s = -1 - s; } // Reversed order is encoded by the sign.
+   if (rev) { s = FlipIndexSign(s); } // Reversed order is encoded by the sign.
   return s;
 }

@@ -3040,7 +3037,7 @@ int GetFaceOrientation(const Mesh *mesh, const int face,
 // see GetFaceOrientation.
 bool Reorder2D(int ori, std::array<int, 2> &s0)
 {
-   const int shift = ori < 0 ? -1 - ori : ori;
+   const int shift = UnsignIndex(ori);

   // Shift is an F1 index in the counter-clockwise ordering of 4 quad vertices.
   // Now find the (i,j) indices of this index, with i,j in {0,1}.
@@ -3064,7 +3061,7 @@ void GetInverseShiftedDimensions2D(int signedShift, int sm, int sn, int &m,
                                   int &n)
 {
   const bool rev = (signedShift < 0);
-   const int shift = rev ? -1 - signedShift : signedShift;
+   const int shift = UnsignIndex(signedShift);
   MFEM_ASSERT(0 <= shift && shift < 4, "");

   // We consider 8 cases for the possible values of rev and shift.
@@ -3136,7 +3133,7 @@ void GetShiftedGridPoints2D(int m, int n, int i, int j, int signedShift,
                            int& sm, int& sn, int& si, int& sj)
 {
   const bool rev = (signedShift < 0);
-   const int shift = rev ? -1 - signedShift : signedShift;
+   const int shift = UnsignIndex(signedShift);
   MFEM_ASSERT(0 <= shift && shift < 4, "");

   // (0,0) <= (i,j) < (m,n) are old indices, and old vertex [shift] maps
@@ -3798,8 +3795,7 @@ void NCNURBSExtension::GenerateOffsets()
      const int signedParentEdge = auxEdges[e].parent;
      const int ki0 = auxEdges[e].ksi[0];
      const int ki1raw = auxEdges[e].ksi[1];
-      const bool rev = signedParentEdge < 0;
-      const int parentEdge = rev ? -1 - signedParentEdge : signedParentEdge;
+      const int parentEdge = UnsignIndex(signedParentEdge);
      const int masterNE = KnotVec(parentEdge)->GetNE();
      const int ki1 = ki1raw == -1 ? masterNE : ki1raw;
      const int perEdgeCP = GetNCPperEdge(KnotVec(e));
@@ -3723,10 +3723,7 @@ bool NURBSExtension::CheckPatches()
      for (int i = 0; i < edges.Size(); i++)
      {
         edges[i] = edge_to_ukv[edges[i]];
-         if (oedge[i] < 0)
-         {
-            edges[i] = -1 - edges[i];
-         }
+         if (oedge[i] < 0) { edges[i] = FlipIndexSign(edges[i]); }
      }

      // In 2d - opposite edges must be same knotvector with opposite sign.
@@ -3740,7 +3737,7 @@ bool NURBSExtension::CheckPatches()
      // {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7} for Geometry::CUBE in 3D
      // See fem/geom.cpp for these definitions.
      if ((dim == 2 &&
-           (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
+           (edges[0] != FlipIndexSign(edges[2]) || edges[1] != FlipIndexSign(edges[3]))) ||

          (dim == 3 &&
           (edges[0] != edges[2] || edges[0] != edges[4] ||
@@ -3769,7 +3766,7 @@ void NURBSExtension::CheckBdrPatches()
         edges[i] = edge_to_ukv[edges[i]];
         if (oedge[i] < 0)
         {
-            edges[i] = -1 - edges[i];
+            edges[i] = FlipIndexSign(edges[i]);
         }
      }

@@ -4766,14 +4763,13 @@ void NURBSExtension::GenerateBdrElementDofTable()
   SetPatchToBdrElements();

   int *dof = bel_dof->GetJ();
-   int ndof = bel_dof->Size_of_connections();
+   const int ndof = bel_dof->Size_of_connections();
   for (int i = 0; i < ndof; i++)
   {
-      int idx = dof[i];
+      const int idx = dof[i];
      if (idx < 0)
      {
-         dof[i] = -1 - (activeDof[-1-idx] - 1);
-         dof[i] = -activeDof[-1-idx];
+         dof[i] = -activeDof[FlipIndexSign(idx)];
      }
      else
      {
@@ -4858,12 +4854,12 @@ void NURBSExtension::Generate2DBdrElementDofTable()
                  for (int ii = 0; ii <= ord0; ii++)
                  {
                     conn.to = DofMap(p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]);
-                     if (s == -1) { conn.to = -1 -conn.to; }
+                     if (s == -1) { conn.to = FlipIndexSign(conn.to); }
                     bel_dof_list.Append(conn);
                  }
               }
               bel_to_patch[lbe] = b;
-               bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
+               bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : FlipIndexSign(i);
               lbe++;
            }
            gbe++;
@@ -4936,14 +4932,14 @@ void NURBSExtension::Generate3DBdrElementDofTable()
                           {
                              const int ii_ = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
                              conn.to = DofMap(p2g(ii_, jj_));
-                              if (s == -1) { conn.to = -1 -conn.to; }
+                              if (s == -1) { conn.to = FlipIndexSign(conn.to); }
                              bel_dof_list.Append(conn);
                           }
                        }
                     }
                     bel_to_patch[lbe] = b;
-                     bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
-                     bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
+                     bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : FlipIndexSign(i);
+                     bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : FlipIndexSign(j);
                     lbe++;
                  }
                  gbe++;
@@ -220,7 +220,7 @@ public:
       @a u.
       The main purpose of this function is its use in FindInterpolant.
       Use GetBotella instead for each shape function separately, perhaps in
-       conjuction with GetSpan and GetRefPoint.*/
+       conjunction with GetSpan and GetRefPoint.*/
   MFEM_DEPRECATED void FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const;

   /** @brief Global curve interpolation through the points @a x (overwritten).
@@ -1398,8 +1398,7 @@ inline const real_t &NURBSPatch::operator()(int i, int j, int k, int l) const

 inline int NURBSExtension::KnotInd(int edge) const
 {
-   const int kv = edge_to_ukv[edge];
-   return kv >= 0 ? kv : -1 - kv;
+   return UnsignIndex(edge_to_ukv[edge]);
 }

 inline int NURBSExtension::KnotSign(int edge) const
@@ -1429,7 +1428,7 @@ const
   else
   {
      *okv = -oedge;
-      return knotVectors[-1-kv];
+      return knotVectors[FlipIndexSign(kv)];
   }
 }

@@ -400,7 +400,7 @@ void ParNCMesh::MakeSharedList(const NCList &list, NCList &shared)
         }
         else // special case: prism edge-face constraint
         {
-            if (entity_owner[1][-1-si] != MyRank)
+            if (entity_owner[1][FlipIndexSign(si)] != MyRank)
            {
               master_flag |= 0x2;
            }
@@ -571,9 +571,10 @@ void ParNCMesh::CalculatePMatrixGroups()
      ranks.SetSize(0);
      for (int j = master_face.slaves_begin; j < master_face.slaves_end; j++)
      {
-         int si = face_list.slaves[j].index;
-         int owner = (si >= 0) ? entity_owner[2][si] // standard face dependency
-                     /*     */ : entity_owner[1][-1 - si]; // prism edge-face dep
+         const int si = face_list.slaves[j].index;
+         const int owner =
+            (si >= 0) ? entity_owner[2][si] : // standard face dependency
+            entity_owner[1][FlipIndexSign(si)]; // prism edge-face dep
         ranks.Append(groups[owner][0]);
      }
      ranks.Sort();
@@ -1181,7 +1182,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
      if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }

      Mesh::FaceInfo &fi = pmesh.faces_info[cf.index];
-      fi.Elem2No = -1 - fnbr_index[e[0]->index - NElements];
+      fi.Elem2No = FlipIndexSign(fnbr_index[e[0]->index - NElements]);

      if (Dim == 3)
      {
@@ -1270,7 +1271,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
               // In other words, side 1 IS the side that generated the face.
            }
            MFEM_ASSERT(fi.Elem2No >= NElements, "");
-            fi.Elem2No = -1 - fnbr_index[fi.Elem2No - NElements];
+            fi.Elem2No = FlipIndexSign(fnbr_index[fi.Elem2No - NElements]);

            const DenseMatrix* pm = full_list.point_matrices[sf.geom][sf.matrix];
            if (!sloc && Dim == 3)
@@ -2286,7 +2287,7 @@ void ParNCMesh::Derefine(const Array<int> &derefs)
      if (element_type[index] == 0)
      {
         // this coarse element will get pruned, encode who owns it now
-         index = -1 - elements[coarse[i]].rank;
+         index = FlipIndexSign(elements[coarse[i]].rank);
      }
      transforms.embeddings[i].parent = index;
   }
@@ -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)
@@ -264,7 +264,7 @@ int main(int argc, char *argv[])
   //    parallel degrees of freedom, with offsets given by array true_offset.
   FiniteElementCollection *fec = nullptr;
   NURBSExtension *NURBSext = nullptr;
-   if (mesh->NURBSext)
+   if (pmesh->NURBSext)
   {
      NURBSext = new NURBSExtension(pmesh->NURBSext, order);
      fec = new NURBSFECollection(order);
@@ -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
+
@@ -39,6 +39,7 @@ set(UNIT_TESTS_SRCS
  dfem/test_divergence.cpp
  dfem/test_lvector_interface.cpp
  dfem/test_mass.cpp
+  dfem/test_tensor_derivatives.cpp
  general/test_array.cpp
  general/test_scan.cpp
  general/test_arrays_by_name.cpp
@@ -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

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

@@ -0,0 +1,48 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#ifdef MFEM_USE_ENZYME
+
+using namespace mfem::future;
+
+tensor<double, 3, 3> Shift(
+    const tensor<double, 3, 3>& X, const tensor<double, 3, 3>& p)
+{
+  return X + p;
+}
+
+TEST_CASE("Enzyme derivatives of tensors", "[tensor][enzyme]")
+{
+    tensor<double, 3, 3> x{{{ 1.0, -1.0,  0.0},
+                            {-1.0,  2.0, -1.0},
+                            { 0.0, -1.0,  1.0}}};
+    tensor p = 3*IdentityMatrix<3>();
+    tensor<double, 3, 3> x_dot{{{1.0, 0.0, 0.0},
+                                {0.0, 0.0, 0.0},
+                                {0.0, 0.0, 0.0}}};
+
+    SECTION("Tensor function with constant parameter")
+    {
+        auto y_dot = __enzyme_fwddiff<tensor<double, 3, 3>>((void*)Shift, enzyme_dup, x, x_dot, enzyme_const, p);
+
+        // correct answer is simply y_dot = x_dot
+        for (int i = 0; i < 3; i++) {
+            for (int j = 0; j < 3; j++) {
+                CHECK(y_dot[i][j] == MFEM_Approx(x_dot[i][j]));        
+            }
+        }
+    }
+}
+
+#endif // MFEM_USE_ENZYME
@@ -151,6 +151,8 @@ TEST_CASE("Collocated Derivative Kernels", "[QuadratureInterpolator]")
      auto L = GENERATE(QVectorLayout::byNODES, QVectorLayout::byVDIM);
      auto P = GENERATE(true, false);

+      CAPTURE(L, P);
+
      const int nd = maps.ndof;
      const int nq = maps.nqpt;

@@ -163,7 +165,15 @@ TEST_CASE("Collocated Derivative Kernels", "[QuadratureInterpolator]")
      CGK::Run(dim, L, P, vdim, nd, nelem, maps.G.Read(), geom->J.Read(),
               evec_values.Read(), col_der.Write(), sdim, vdim, nd);

+      const real_t max_norm = qp_der.Normlinf();
+
      qp_der -= col_der;
-      REQUIRE(qp_der.Normlinf() == MFEM_Approx(0.0, 1e-10, 1e-10));
+
+      const real_t abs_err = qp_der.Normlinf();
+      const real_t rel_err = max_norm > 0_r ?
+                             abs_err/max_norm :
+                             abs_err > 0_r ? mfem::infinity() : 0_r;
+      CAPTURE(rel_err, max_norm);
+      CHECK(rel_err <= 1e-13);
   }
 }
@@ -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))
@@ -125,6 +125,59 @@ TEST_CASE("Array stl-interactions", "[Array]")
   }
 }

+
+TEST_CASE("Array move assignment to view (MakeRef)", "[Array]")
+{
+   constexpr int n = 5;
+
+   // Helper function that returns an Array by value (rvalue)
+   auto make_sequence_array = [](int n)
+   {
+      Array<int> arr(n);
+      for (int i = 0; i < n; i++)
+      {
+         arr[i] = i + 1;  // 1, 2, 3, ...
+      }
+      return arr;
+   };
+
+   // Create backing storage
+   Array<int> backing1(n), backing2(n);
+   backing1 = 0.0; backing2 = 0.0;
+
+   // Create a view into the backing storage
+   Array<int> view1, view2;
+   view1.MakeRef(backing1);
+   view2.MakeRef(backing2);
+
+   auto seq_array = make_sequence_array(n);
+   view1 = seq_array; // copy assign
+   view2 = std::move(seq_array); // move assign
+
+   CHECK(seq_array.Size() == 0); // seq_array is invalidated by the move
+
+   // Both assignments should be semantically equivalent.
+   for (int i = 0; i < n; i++)
+   {
+      CHECK(backing1[i] == backing2[i]);
+   }
+
+   // Create backing storage as raw array
+   int backing3[n] = {0, 0, 0, 0, 0};
+
+   // Create a view into the backing storage
+   Array<int> view3;
+   view3.MakeRef(backing3, n);
+
+   // Assign from rvalue
+   view3 = std::move(view2);
+
+   for (int i = 0; i < n; i++)
+   {
+      CHECK(backing3[i] == backing1[i]);
+   }
+}
+
 TEST_CASE("Array delete at indices", "[Array],[GPU]")
 {
   for (int use_dev = 0; use_dev < 2; use_dev++)
@@ -961,3 +961,38 @@ TEST_CASE("NNLS", "[DenseMatrix]")
 }

 #endif // if MFEM_USE_LAPACK
+
+TEST_CASE("DenseTensor slice copy and move assign equivalency",
+          "[DenseMatrix][DenseTensor]")
+{
+   auto fill_matrix = [](int n)
+   {
+      DenseMatrix M(n, n);
+      M = 0.0;
+      for (int i = 0; i < n; i++)
+         for (int j = 0; j < n; j++)
+         {
+            M(i, j) = i+j*n+1;
+         }
+      return M;
+   };
+
+   constexpr int n = 3;
+   constexpr int k = 2;
+
+   DenseTensor tensor1(n, n, k), tensor2(n, n, k);
+   tensor1 = 0.0; tensor2 = 0.0;
+
+   DenseMatrix temp = fill_matrix(n);
+   tensor1(0) = temp; // copy assign
+   tensor2(0) = std::move(temp); // move assign
+
+   // Check that the tensor was actually updated
+   for (int i = 0; i < n; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         CHECK(tensor1(i, j, 0) == tensor2(i, j, 0));
+      }
+   }
+}
@@ -222,13 +222,13 @@ TEST_CASE("Greville, Botella and Demko points", "[NURBS]")
      mfem::out<<"Knotvector : "; kvp.Print(mfem::out);

      kvp.GetGreville(xi);
-      mfem::out<<"Greville points : "; xi.Print(std::cout,999);
+      mfem::out<<"Greville points : "; xi.Print(mfem::out,999);

      kvp.GetBotella(xi);
-      mfem::out<<"Botella  points : "; xi.Print(std::cout,999);
+      mfem::out<<"Botella  points : "; xi.Print(mfem::out,999);

      kvp.GetDemko(xi);
-      mfem::out<<"Demko    points : "; xi.Print(std::cout,999);
+      mfem::out<<"Demko    points : "; xi.Print(mfem::out,999);
   }

   KnotVector kv(3, Vector({0.0, 0.3, 0.3, 0.3, 0.6, 1.0}));
Author	SHA1	Message	Date
Brandon Talamini	6fd649f1c2	Fix style	2026-03-23 08:44:36 -07:00
Brandon Talamini	02a66e39b6	Make unit test that reproduces enzyme error in debug build	2026-03-21 10:22:20 -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	82abd48bba	Merge pull request #5080 from mfem/cmake-config CMake config.mk for CUDA and HIP	2026-03-11 12:18:25 -04: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
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
Tzanio Kolev	5ce6e90ceb	Merge pull request #5253 from mfem/bugfix-nurbs-orientation NURBS orientation bug fix	2026-03-01 11:46:34 -08:00
Tzanio Kolev	03ba184adb	Merge pull request #4882 from mfem/debug-mem-silent Add option to run debug memory backend without issuing errors	2026-03-01 11:46:15 -08:00
Tzanio Kolev	cfa87477da	Merge pull request #5248 from lindsayad/dont-do-math-with-enums Don't do arithmetic with enums	2026-03-01 11:45:47 -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
dylan-copeland	939932bc68	Merge branch 'master' of github.com:mfem/mfem into bugfix-nurbs-orientation	2026-02-27 13:59:37 -08:00
Tzanio Kolev	d3ae34710c	Merge pull request #5195 from mfem/flip-index Index sign functions	2026-02-27 13:32:22 -08:00
Dylan Copeland	9b652996b2	Remove an unnecessary assertion.	2026-02-27 11:48:01 -08:00
Dylan Copeland	5cfbbe5fad	Bug fix for NURBS orientation, when the same knotvector index is used on all edges of some patches.	2026-02-27 11:24:36 -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
Tzanio Kolev	71909cd5e3	Merge pull request #5203 from mfem/bugfix/unit-test-cout Changing std::cout to mfem::out in unit test	2026-02-26 05:50:05 -08:00
Alex Lindsay	a5a3169064	Don't do arithmetic with enums Else with gcc 13.3 with `-std=c++20` I get warnings	2026-02-25 13:35:10 -07: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
Dylan Copeland	afded067a7	Replace absdof.	2026-02-24 11:17:54 -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
John Camier	629e93afd9	Merge branch 'master' into debug-mem-silent	2026-02-24 09:20:59 -08:00
Tzanio Kolev	3e277808a9	Merge pull request #5210 from mfem/hughcars/array-move-assignment-bugfix Fix Array non-owning move assignment	2026-02-24 08:41:48 -08:00
Tzanio Kolev	864fb1ce9e	Merge pull request #5217 from mfem/quadrature-function-fix Fix a bug in `QuadratureFunction::GetValues` that returns a `DenseMatrix` view	2026-02-24 08:40:53 -08:00
Tzanio Kolev	ac7415cc69	Merge branch 'master' into bugfix/unit-test-cout	2026-02-24 08:38:11 -08:00
Veselin Dobrev	38030d4395	Merge pull request #5222 from mfem/bugfix/chapman39/rm-return-from-omp Remove return from openmp section	2026-02-24 08:18:53 -08:00
Alex Tyler Chapman	2c96dc6a1f	Merge branch 'master' into bugfix/chapman39/rm-return-from-omp	2026-02-23 11:53:58 -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
John Camier	ecbc7bf8c2	Merge branch 'master' into bugfix/unit-test-cout	2026-02-23 08:00:52 -08:00
John Camier	5b5a21edac	Merge branch 'master' into debug-mem-silent	2026-02-22 12:28:39 -08:00
Tzanio Kolev	1f84ba036e	Merge pull request #5179 from izaid/pfes Fixed NULL pointer segfault in GetSurfaceFittingErrors	2026-02-21 12:52:33 -08:00
Tzanio Kolev	76d0312309	Merge branch 'master' into hughcars/array-move-assignment-bugfix	2026-02-21 12:52:11 -08:00
Tzanio Kolev	c7ed339260	Merge pull request #5228 from mfem/mesh-3d-part-fix Improve 3d mesh partitions	2026-02-21 12:43:49 -08:00
Tzanio Kolev	f1b3a33fb2	Merge pull request #5234 from mfem/few-small-fixes A few small fixes	2026-02-21 12:43:02 -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
Veselin Dobrev	e3dedbbd5b	Increase a relative tolerance due to failures on some machines	2026-02-18 10:27:44 -08:00
John Camier	261f3805b8	Merge branch 'master' into debug-mem-silent	2026-02-17 10:59:40 -08:00
John Camier	b386b2d6b6	Merge branch 'master' into bugfix/chapman39/rm-return-from-omp	2026-02-17 08:28:34 -08:00
Tzanio Kolev	6e98055eb7	Merge branch 'master' into quadrature-function-fix	2026-02-17 08:28:12 -08:00
Tzanio Kolev	2b14134496	Merge branch 'master' into najlkin/fix-cmplx-assign	2026-02-17 08:27:01 -08:00
Veselin Dobrev	75012728db	Formatting: re-wrap comment to 80 chars/line.	2026-02-14 13:10:30 -08:00
Dylan Copeland	1f094244f8	Typo	2026-02-13 10:37:58 -08:00
Dylan Copeland	cc16ddadbf	Merge branch 'master' of github.com:mfem/mfem into flip-index	2026-02-13 10:35:32 -08:00
Veselin Dobrev	8acd5cd3a2	Fix a size bug in thread-safe mode in H1_TriangleElement::CalcHessian. Fix use-after-delete bug in nurbs_ex10p.cpp. Use relative tolerance in the "Collocated Derivative Kernels" unit test to resolve failures in some setups with the original absolute tolerance.	2026-02-13 10:29:22 -08:00
chapman39	1f5bc1c3d8	wording	2026-02-13 10:03:00 -08:00
chapman39	3588d47ec1	Merge remote-tracking branch 'origin/master' into bugfix/chapman39/rm-return-from-omp	2026-02-13 09:59:31 -08:00
Alex Tyler Chapman	f0fe5b0ec0	Merge branch 'master' into bugfix/chapman39/rm-return-from-omp	2026-02-11 15:14:48 -08:00
chapman39	a92983051a	style	2026-02-11 15:14:23 -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
Vladimir Z Tomov	6bb6745c0e	improve 3d partitions	2026-02-10 16:23:38 -08:00
chapman39	32fb4bf244	remove return from openmp section	2026-02-05 14:56:40 -08:00
Veselin Dobrev	2d021685de	Fix a bug in QuadratureFunction::GetValues when getting a reference to the data for one mesh element as a DenseMatrix. Added new methods: - Array<T>::MakeRef(Memory<T> &base, int offset, int size_) - DenseMatrix::MakeRef(Memory<real_t> &base, int offset, int h, int w)	2026-02-02 08:42:58 -08:00
Hugh Carson	d6ea262498	Address PR feedback - Hoist src deletion - Remove unneeded explicit cast given src is named	2026-01-30 14:01:59 -05:00
Andrew Ho	26e9057f02	revert change, updated comment to why libdl gets special treatment	2026-01-30 07:27:04 -08:00
Dylan Copeland	6b147fd9ff	Merge branch 'master' of github.com:mfem/mfem into flip-index	2026-01-29 15:45:52 -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
Hugh Carson	a0981cb363	Refactor test for showing equivalency of copy and move assignment	2026-01-27 14:15:08 -05:00
Hugh Carson	894de992da	Move assignment of non-owned arrays must fallback to copy assignment	2026-01-27 14:01:23 -05:00
Tzanio Kolev	942249395b	Merge branch 'master' into debug-mem-silent	2026-01-24 11:40:50 -08:00
Stowell, Mark L.	a3f6d5b971	Changing std::cout to mfem::out	2026-01-22 15:37:51 -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	5b065ad7f2	Merge branch 'master' into debug-mem-silent	2026-01-20 17:59:07 -08:00
Dylan Copeland	0161ad9d92	Minor simplifications.	2026-01-14 12:22:47 -08:00
Dylan Copeland	3589479481	Merge branch 'master' of github.com:mfem/mfem into flip-index	2026-01-14 12:10:30 -08:00
Dylan Copeland	e9acfeccda	Using new sign functions.	2026-01-14 12:10:07 -08:00
John Camier	caa973d6a0	Merge branch 'master' into cmake-config	2026-01-13 08:13:47 -08:00
izaid	5cd3ec521b	applied astyle	2026-01-08 23:19:16 +00:00
Tzanio Kolev	9488637956	Merge branch 'master' into pfes	2026-01-07 14:20:49 -08:00
Andrew Ho	ad40704e20	Merge branch 'master' into cmake-config	2026-01-06 11:55:40 -08:00
Dylan Copeland	6f280d81b5	Introducing new functions for flipping index signs.	2026-01-05 18:46:25 -08:00
izaid	90353c437e	Fixed bug in GetSurfaceFittingErrors	2026-01-01 23:09:09 +01:00
John Camier	fbb50af208	Merge branch 'master' into debug-mem-silent	2025-12-30 12:50:33 -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
John Camier	ce8cd01cfd	Merge branch 'master' into debug-mem-silent	2025-10-28 10:20:51 -07:00
John Camier	a60baf8ce6	Merge branch 'master' into debug-mem-silent	2025-10-25 11:28:07 -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
John Camier	96bba18449	Merge branch 'master' into debug-mem-silent	2025-10-20 09:36:28 -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
John Camier	f9cce3ab62	Merge branch 'master' into debug-mem-silent	2025-09-29 09:25:14 -07:00
John Camier	292700bb52	Merge branch 'master' into debug-mem-silent	2025-09-25 15:58:22 -07:00
John Camier	fd2f0df34f	Merge branch 'master' into debug-mem-silent	2025-09-21 20:13:41 -07:00
John Camier	4ac41a6427	Merge branch 'master' into debug-mem-silent	2025-09-21 10:27:05 -07:00
Tzanio Kolev	58826d64c9	Merge branch 'master' into debug-mem-silent	2025-08-31 15:32:36 -07:00
Victor DeCaria	a126203ccd	implement suggested change	2025-08-21 15:08:57 -06:00
victor-decaria-nnl	6656a7ef72	Merge branch 'master' into debug-mem-silent	2025-08-20 15:50:34 -04:00
Mark L. Stowell	83fd119b95	Merge branch 'master' into phys-range-dim	2025-07-01 10:16:58 -07:00
John Camier	8804df317d	Merge branch 'master' into debug-mem-silent	2025-07-01 09:51:06 -07:00
John Camier	95408b0fae	Merge branch 'master' into debug-mem-silent	2025-06-22 20:25:51 -07:00
John Camier	79819a5563	Merge branch 'master' into debug-mem-silent	2025-06-18 10:22:04 -07:00
John Camier	d57fc7c0d9	Merge branch 'master' into debug-mem-silent	2025-06-16 08:33:43 -07:00
John Camier	9fb590d79d	Merge branch 'master' into debug-mem-silent	2025-06-13 08:48:05 -07:00
Victor DeCaria	cb4ca9228f	add option to disable protection with debug backend	2025-06-05 09:34:25 -06: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