diff --git a/CHANGELOG.md b/CHANGELOG.md
index d2f51137e..240d76acd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,7 +11,7 @@ The format of this changelog is based on
 [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to
 [Semantic Versioning](https://semver.org/).
 
-## In progress
+## [0.12.0] - 2023-12-21
 
   - Added support for operator partial assembly for high-order finite element spaces based
     on libCEED for mixed and non-tensor product element meshes. This option is disabled by
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 70eff7ce2..992054f01 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,7 +14,7 @@ if(CMAKE_SOURCE_DIR STREQUAL CMAKE_BINARY_DIR)
 endif()
 
 # Initialize the project
-project(palace-superbuild LANGUAGES CXX C VERSION 0.11.2)
+project(palace-superbuild LANGUAGES CXX C VERSION 0.12.0)
 
 # Define build settings and defaults
 set(PALACE_WITH_64BIT_INT OFF CACHE BOOL "Use 64 bit integers")
diff --git a/cmake/ExternalBLASLAPACK.cmake b/cmake/ExternalBLASLAPACK.cmake
index 8ad5dc184..501fa540b 100644
--- a/cmake/ExternalBLASLAPACK.cmake
+++ b/cmake/ExternalBLASLAPACK.cmake
@@ -41,28 +41,29 @@ if(DEFINED ENV{ARMPL_DIR} OR DEFINED ENV{ARMPLROOT} OR DEFINED ENV{ARMPL_ROOT})
     endif()
   else()
     if(PALACE_WITH_OPENMP)
-      set(ARMPL_LIB_SUFFIX "_lp64_mp")
+      set(ARMPL_LIB_SUFFIX "_mp")
     else()
-      set(ARMPL_LIB_SUFFIX "_lp64")
+      set(ARMPL_LIB_SUFFIX "")
     endif()
   endif()
-  find_library(_BLAS_LAPACK_LIBRARIES
-    NAMES armpl${ARMPL_LIB_SUFFIX} armpl
-    PATHS ${ARMPL_DIR}
-    PATH_SUFFIXES lib lib64
-    NO_DEFAULT_PATH
-    REQUIRED
-  )
+  list(APPEND CMAKE_PREFIX_PATH ${ARMPL_DIR})
+  set(BLA_VENDOR "Arm${ARMPL_LIB_SUFFIX}")
+  find_package(BLAS REQUIRED)
+  find_package(LAPACK REQUIRED)
+  set(LAPACK_LIBRARIES "${LAPACK_LIBRARIES};-lm")
+
+  # Locate include directory
   find_path(_BLAS_LAPACK_INCLUDE_DIRS
-    NAMES cblas.h
+    NAMES armpl.h
     PATHS ${ARMPL_DIR}
-    PATH_SUFFIXES include${ARMPL_LIB_SUFFIX} include
+    PATH_SUFFIXES include
     NO_DEFAULT_PATH
     REQUIRED
   )
   message(STATUS "Using BLAS/LAPACK from Arm Performance Libraries (Arm PL)")
 elseif(DEFINED ENV{AOCL_DIR} OR DEFINED ENV{AOCLROOT} OR DEFINED ENV{AOCL_ROOT})
-  # AOCL for x86_64 builds when available
+  # AOCL for x86_64 builds when available (part of CMake's FindBLAS/FindLAPACK as of v3.27
+  # but unnecessarily adds -fopenmp flag)
   if(CMAKE_SYSTEM_PROCESSOR MATCHES "aarch64|arm")
     message(WARNING "AOCL math libraries are not intended for arm64 architecture builds")
   endif()
@@ -85,21 +86,22 @@ elseif(DEFINED ENV{AOCL_DIR} OR DEFINED ENV{AOCLROOT} OR DEFINED ENV{AOCL_ROOT})
   else()
     set(AOCL_LIB_SUFFIX "")
   endif()
-  find_library(BLIS_LIBRARY
+  find_library(BLAS_LIBRARIES
     NAMES blis${AOCL_LIB_SUFFIX} blis
     PATHS ${AOCL_DIR}
     PATH_SUFFIXES lib${AOCL_DIR_SUFFIX} lib lib64
     NO_DEFAULT_PATH
     REQUIRED
   )
-  find_library(FLAME_LIBRARY
+  find_library(LAPACK_LIBRARIES
     NAMES flame FLAME
     PATHS ${AOCL_DIR}
     PATH_SUFFIXES lib${AOCL_DIR_SUFFIX} lib lib64
     NO_DEFAULT_PATH
     REQUIRED
   )
-  set(_BLAS_LAPACK_LIBRARIES "${FLAME_LIBRARY}$<SEMICOLON>${BLIS_LIBRARY}")
+
+  # Locate include directory
   find_path(_BLAS_LAPACK_INCLUDE_DIRS
     NAMES cblas.h
     PATHS ${AOCL_DIR}
@@ -139,9 +141,8 @@ elseif(DEFINED ENV{MKL_DIR} OR DEFINED ENV{MKLROOT} OR DEFINED ENV{MKL_ROOT})
   set(BLA_VENDOR "Intel10${MKL_LIB_SUFFIX}")
   find_package(BLAS REQUIRED)
   find_package(LAPACK REQUIRED)
-  set(_BLAS_LAPACK_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
-  list(REMOVE_DUPLICATES _BLAS_LAPACK_LIBRARIES)
-  string(REPLACE ";" "$<SEMICOLON>" _BLAS_LAPACK_LIBRARIES "${_BLAS_LAPACK_LIBRARIES}")
+
+  # Locate include directories
   find_path(_BLAS_LAPACK_INCLUDE_DIRS
     NAMES mkl_cblas.h
     PATHS ${MKL_DIR}
@@ -165,11 +166,10 @@ else()
   list(APPEND CMAKE_PREFIX_PATH ${OPENBLAS_DIR})
   find_package(BLAS REQUIRED)
   find_package(LAPACK REQUIRED)
-  set(_BLAS_LAPACK_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
-  list(REMOVE_DUPLICATES _BLAS_LAPACK_LIBRARIES)
-  string(REPLACE ";" "$<SEMICOLON>" _BLAS_LAPACK_LIBRARIES "${_BLAS_LAPACK_LIBRARIES}")
+
+  # Locate include directory
   set(_BLAS_LAPACK_DIRS)
-  foreach(LIB IN LISTS _BLAS_LAPACK_LIBRARIES)
+  foreach(LIB IN LISTS LAPACK_LIBRARIES BLAS_LIBRARIES)
     get_filename_component(LIB_DIR ${LIB} DIRECTORY)
     list(APPEND _BLAS_LAPACK_DIRS ${LIB_DIR})
   endforeach()
@@ -184,6 +184,9 @@ else()
 endif()
 
 # Save variables to cache
+set(_BLAS_LAPACK_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
+list(REMOVE_DUPLICATES _BLAS_LAPACK_LIBRARIES)
+string(REPLACE ";" "$<SEMICOLON>" _BLAS_LAPACK_LIBRARIES "${_BLAS_LAPACK_LIBRARIES}")
 set(BLAS_LAPACK_LIBRARIES ${_BLAS_LAPACK_LIBRARIES} CACHE STRING
   "List of library files for BLAS/LAPACK"
 )
diff --git a/cmake/ExternalGitTags.cmake b/cmake/ExternalGitTags.cmake
index 5cf962d59..e29c492f4 100644
--- a/cmake/ExternalGitTags.cmake
+++ b/cmake/ExternalGitTags.cmake
@@ -34,7 +34,7 @@ set(EXTERN_BUTTERFLYPACK_GIT_BRANCH
   "Git branch for external ButterflyPACK build"
 )
 set(EXTERN_BUTTERFLYPACK_GIT_TAG
-  "893a5875e2e3a8f4abd5c7374397c773eb0c9724" CACHE STRING  # 12/20/2023
+  "37fc1a8d09e2a1eed9a6ec6cf9494ac14bc6c700" CACHE STRING  # 12/20/2023
   "Git tag for external ButterflyPACK build"
 )
 
@@ -146,7 +146,7 @@ set(EXTERN_MFEM_GIT_BRANCH
   "Git branch for external MFEM build"
 )
 set(EXTERN_MFEM_GIT_TAG
-  "8a4affef5af84d74c9691f76b49624bdd7199d60" CACHE STRING  # master @ 12/17/2023
+  "c1871afea4dba50b5c4944adc20cfcb17520bd77" CACHE STRING  # master @ 12/21/2023
   "Git tag for external MFEM build"
 )
 
@@ -160,7 +160,7 @@ set(EXTERN_MUMPS_GIT_BRANCH
   "Git branch for external MUMPS build"
 )
 set(EXTERN_MUMPS_GIT_TAG
-  "44fb7e74b5780a3ea1c26a76a55298a2c34279fe" CACHE STRING  # 12/08/2023
+  "88148442304bf9fb82d04d9e8bf3c3dccef3c919" CACHE STRING  # 12/20/2023
   "Git tag for external MUMPS build"
 )
 
@@ -188,7 +188,7 @@ set(EXTERN_PETSC_GIT_BRANCH
   "Git branch for external PETSc build"
 )
 set(EXTERN_PETSC_GIT_TAG
-  "9617a9bff32dac3e651f1a88b0773cde31a476b6" CACHE STRING  # 12/19/2023
+  "cec8e1ee48d2b6b0308890436f9bfafab17cd824" CACHE STRING  # 12/21/2023
   "Git tag for external PETSc build"
 )
 
@@ -216,7 +216,7 @@ set(EXTERN_SLATE_GIT_BRANCH
   "Git branch for external SLATE build"
 )
 set(EXTERN_SLATE_GIT_TAG
-  "f9ac8c71f7dcc918cf3fb4c3955e2dab924ee654" CACHE STRING  # 12/14/2022
+  "aaf28f7410a4832d04d9d6b115eea1e254adf1b0" CACHE STRING  # 12/20/2022
   "Git tag for external SLATE build"
 )
 
diff --git a/cmake/ExternalMFEM.cmake b/cmake/ExternalMFEM.cmake
index 30832649d..072c1fa70 100644
--- a/cmake/ExternalMFEM.cmake
+++ b/cmake/ExternalMFEM.cmake
@@ -359,10 +359,9 @@ message(STATUS "MFEM_OPTIONS: ${MFEM_OPTIONS_PRINT}")
 set(MFEM_PATCH_FILES
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_mfem_device_fixes.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_mesh_vis_dev.diff"
-  "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_par_tet_mesh_fix.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_mesh_partitioner_dev.diff"
+  "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_par_tet_mesh_fix.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_ncmesh_interior_boundary_dev.diff"
-  "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_gslib_leak_fix.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_mesh_const_fix.diff"
 )
 
diff --git a/extern/patch/mfem/patch_gslib_leak_fix.diff b/extern/patch/mfem/patch_gslib_leak_fix.diff
deleted file mode 100644
index 8d9442d59..000000000
--- a/extern/patch/mfem/patch_gslib_leak_fix.diff
+++ /dev/null
@@ -1,12 +0,0 @@
-diff --git a/fem/gslib.cpp b/fem/gslib.cpp
-index 3807694d0..99f12fe35 100644
---- a/fem/gslib.cpp
-+++ b/fem/gslib.cpp
-@@ -67,6 +67,7 @@ FindPointsGSLIB::FindPointsGSLIB()
- 
- FindPointsGSLIB::~FindPointsGSLIB()
- {
-+   comm_free(gsl_comm);
-    delete gsl_comm;
-    delete cr;
-    for (int i = 0; i < 4; i++)
diff --git a/extern/patch/mfem/patch_mesh_const_fix.diff b/extern/patch/mfem/patch_mesh_const_fix.diff
index 45fd7d64d..fe4857a8a 100644
--- a/extern/patch/mfem/patch_mesh_const_fix.diff
+++ b/extern/patch/mfem/patch_mesh_const_fix.diff
@@ -25,7 +25,7 @@ index 8d3955ab4..3d0ccb97f 100644
     ElementTransformation();
  
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index a46277750..627fc17e5 100644
+index a8ff93b24..848992dbe 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
 @@ -83,7 +83,7 @@ void Mesh::GetElementCenter(int i, Vector &center)
@@ -446,7 +446,7 @@ index a46277750..627fc17e5 100644
  Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
  {
     FaceInformation face;
-@@ -10826,7 +10867,7 @@ int Mesh::FindCoarseElement(int i)
+@@ -10829,7 +10870,7 @@ int Mesh::FindCoarseElement(int i)
     return coarse;
  }
  
@@ -456,7 +456,7 @@ index a46277750..627fc17e5 100644
     MFEM_VERIFY(GetLastOperation() == Mesh::REFINE, "");
  
 diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
-index 58e77d388..1972460e2 100644
+index 3eee6c415..3c43cbb97 100644
 --- a/mesh/mesh.hpp
 +++ b/mesh/mesh.hpp
 @@ -239,7 +239,7 @@ protected:
@@ -468,8 +468,8 @@ index 58e77d388..1972460e2 100644
  
     // Nodes are only active for higher order meshes, and share locations with
     // the vertices, plus all the higher- order control points within the
-@@ -436,39 +436,34 @@ protected:
-    void PrintTopo(std::ostream &os, const Array<int> &e_to_k) const;
+@@ -437,39 +437,34 @@ protected:
+                   const std::string &comments = "") const;
  
     /// Used in GetFaceElementTransformations (...)
 -   void GetLocalPtToSegTransformation(IsoparametricTransformation &, int);
@@ -528,7 +528,7 @@ index 58e77d388..1972460e2 100644
  
     /// Returns the orientation of "base" relative to "test"
     /// In other words: GetTriOrientation(test, base) should equal
-@@ -481,7 +476,7 @@ protected:
+@@ -482,7 +477,7 @@ protected:
     static int ComposeTriOrientations(int ori_a_b, int ori_b_c);
  
     /// Returns the orientation of "test" relative to "base"
@@ -537,7 +537,7 @@ index 58e77d388..1972460e2 100644
  
     /// Returns the orientation of "base" relative to "test"
     /// In other words: GetQuadOrientation(test, base) should equal
-@@ -494,7 +489,7 @@ protected:
+@@ -495,7 +490,7 @@ protected:
     static int ComposeQuadOrientations(int ori_a_b, int ori_b_c);
  
     /// Returns the orientation of "test" relative to "base"
@@ -546,7 +546,7 @@ index 58e77d388..1972460e2 100644
  
     static void GetElementArrayEdgeTable(const Array<Element*> &elem_array,
                                          const DSTable &v_to_v,
-@@ -610,8 +605,6 @@ protected:
+@@ -614,8 +609,6 @@ protected:
     void GetElementData(const Array<Element*> &elem_array, int geom,
                         Array<int> &elem_vtx, Array<int> &attr) const;
  
@@ -555,7 +555,7 @@ index 58e77d388..1972460e2 100644
     // Internal helper used in MakeSimplicial (and ParMesh::MakeSimplicial).
     void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal);
  
-@@ -1161,7 +1154,7 @@ public:
+@@ -1165,7 +1158,7 @@ public:
  
         In parallel, the result takes into account elements on all processors.
     */
@@ -564,7 +564,7 @@ index 58e77d388..1972460e2 100644
  
     /// Checks if the mesh has boundary elements
     virtual bool HasBoundaryElements() const { return (NumOfBdrElements > 0); }
-@@ -1253,7 +1246,7 @@ public:
+@@ -1257,7 +1250,7 @@ public:
  
     /// @brief Return pointer to the i'th element object
     ///
@@ -573,7 +573,7 @@ index 58e77d388..1972460e2 100644
     ///
     /// In parallel, @a i is the local element index which is in the
     /// same range mentioned above.
-@@ -1268,7 +1261,7 @@ public:
+@@ -1272,7 +1265,7 @@ public:
  
     /// @brief Return pointer to the i'th boundary element object
     ///
@@ -582,7 +582,7 @@ index 58e77d388..1972460e2 100644
     ///
     /// In parallel, @a i is the local boundary element index which is
     /// in the same range mentioned above.
-@@ -1281,6 +1274,9 @@ public:
+@@ -1285,6 +1278,9 @@ public:
     /// the Element object itself should not be deleted by the caller.
     Element *GetBdrElement(int i) { return boundary[i]; }
  
@@ -592,7 +592,7 @@ index 58e77d388..1972460e2 100644
     const Element *GetFace(int i) const { return faces[i]; }
  
     /// @}
-@@ -1374,6 +1370,8 @@ public:
+@@ -1378,6 +1374,8 @@ public:
  
     double GetElementSize(int i, const Vector &dir);
  
@@ -601,7 +601,7 @@ index 58e77d388..1972460e2 100644
     double GetElementVolume(int i);
  
     void GetElementCenter(int i, Vector &center);
-@@ -1393,12 +1391,12 @@ public:
+@@ -1397,12 +1395,12 @@ public:
        Geometry::Type geom_buf[Geometry::NumGeom];
     public:
        /// Construct a GeometryList of all element geometries in @a mesh.
@@ -616,7 +616,7 @@ index 58e77d388..1972460e2 100644
           : Array<Geometry::Type>(geom_buf, Geometry::NumGeom)
        { mesh.GetGeometries(dim, *this); }
     };
-@@ -1555,7 +1553,8 @@ public:
+@@ -1559,7 +1557,8 @@ public:
  
     /// Builds the transformation defining the i-th element in @a ElTr.
     /// @a ElTr must be allocated in advance and will be owned by the caller.
@@ -626,7 +626,7 @@ index 58e77d388..1972460e2 100644
  
     /// Returns a pointer to the transformation defining the i-th element.
     ///
-@@ -1568,7 +1567,7 @@ public:
+@@ -1572,7 +1571,7 @@ public:
     /// assuming position of the vertices/nodes are given by @a nodes.
     /// @a ElTr must be allocated in advance and will be owned by the caller.
     void GetElementTransformation(int i, const Vector &nodes,
@@ -635,7 +635,7 @@ index 58e77d388..1972460e2 100644
  
     /// Returns a pointer to the transformation defining the i-th boundary
     /// element.
-@@ -1579,11 +1578,18 @@ public:
+@@ -1583,11 +1582,18 @@ public:
  
     /// Builds the transformation defining the i-th boundary element in @a ElTr.
     /// @a ElTr must be allocated in advance and will be owned by the caller.
@@ -656,7 +656,7 @@ index 58e77d388..1972460e2 100644
  
     /** @brief A helper method that constructs a transformation from the
         reference space of a face to the reference space of an element. */
-@@ -1592,17 +1598,11 @@ public:
+@@ -1596,17 +1602,11 @@ public:
         loc_face_orientation. */
     void GetLocalFaceTransformation(int face_type, int elem_type,
                                     IsoparametricTransformation &Transf,
@@ -676,7 +676,7 @@ index 58e77d388..1972460e2 100644
  
     /// Returns a pointer to the transformation defining the given edge element.
     /// @note The returned object is owned by the class and is shared, i.e.,
-@@ -1644,21 +1644,44 @@ public:
+@@ -1648,21 +1648,44 @@ public:
     /// @note The returned object is owned by the class and is shared, i.e.,
     /// calling this function resets pointers obtained from previous calls.
     /// Also, this pointer should NOT be deleted by the caller.
@@ -732,7 +732,7 @@ index 58e77d388..1972460e2 100644
  
     /// @}
  
-@@ -2106,7 +2129,7 @@ public:
+@@ -2110,7 +2133,7 @@ public:
  
     /** Return fine element transformations following a mesh refinement.
         Space uses this to construct a global interpolation matrix. */
@@ -742,7 +742,7 @@ index 58e77d388..1972460e2 100644
     /// Return type of last modification of the mesh.
     Operation GetLastOperation() const { return last_operation; }
 diff --git a/mesh/ncmesh.cpp b/mesh/ncmesh.cpp
-index 5c57adc4a..3ef57cc55 100644
+index 2074d3048..19af5028b 100644
 --- a/mesh/ncmesh.cpp
 +++ b/mesh/ncmesh.cpp
 @@ -4112,7 +4112,7 @@ const NCMesh::PointMatrix& NCMesh::GetGeomIdentity(Geometry::Type geom)
@@ -807,7 +807,7 @@ index 5c57adc4a..3ef57cc55 100644
                             find_node(el, fa.p1),
                             find_node(el, fa.p2),
 diff --git a/mesh/ncmesh.hpp b/mesh/ncmesh.hpp
-index b004cf43e..5874ac644 100644
+index 1e04d16bb..64a56beee 100644
 --- a/mesh/ncmesh.hpp
 +++ b/mesh/ncmesh.hpp
 @@ -337,13 +337,13 @@ public:
@@ -826,7 +826,7 @@ index b004cf43e..5874ac644 100644
  
     /// Free all internal data created by the above three functions.
     void ClearTransforms();
-@@ -1034,15 +1034,15 @@ protected: // implementation
+@@ -1036,15 +1036,15 @@ protected: // implementation
     static const PointMatrix& GetGeomIdentity(Geometry::Type geom);
  
     void GetPointMatrix(Geometry::Type geom, const char* ref_path,
@@ -845,7 +845,7 @@ index b004cf43e..5874ac644 100644
  
     /// state of leaf_elements before Refine(), set by MarkCoarseLevel()
     Array<int> coarse_elements;
-@@ -1067,7 +1067,7 @@ protected: // implementation
+@@ -1069,7 +1069,7 @@ protected: // implementation
  
     int GetEdgeMaster(int node) const;
  
@@ -855,7 +855,7 @@ index b004cf43e..5874ac644 100644
     /**
      * @brief Return the number of splits of this edge that have occurred in the
 diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
-index 5b5657546..e9f791569 100644
+index 5225cfc14..b6424d4e8 100644
 --- a/mesh/pmesh.cpp
 +++ b/mesh/pmesh.cpp
 @@ -1947,65 +1947,6 @@ void ParMesh::GenerateOffsets(int N, HYPRE_BigInt loc_sizes[],
@@ -1250,7 +1250,7 @@ index 5b5657546..e9f791569 100644
  
  int ParMesh::GetNSharedFaces() const
 diff --git a/mesh/pmesh.hpp b/mesh/pmesh.hpp
-index a256e6ad5..45e343f7d 100644
+index d723204f6..b1235ccbf 100644
 --- a/mesh/pmesh.hpp
 +++ b/mesh/pmesh.hpp
 @@ -91,8 +91,6 @@ protected:
diff --git a/extern/patch/mfem/patch_mesh_partitioner_dev.diff b/extern/patch/mfem/patch_mesh_partitioner_dev.diff
index 0f82eab97..81575e3c7 100644
--- a/extern/patch/mfem/patch_mesh_partitioner_dev.diff
+++ b/extern/patch/mfem/patch_mesh_partitioner_dev.diff
@@ -1,5 +1,5 @@
 diff --git a/.gitignore b/.gitignore
-index 739c62d06..6a16c0d9d 100644
+index f7e7265d6..1e8f32ff6 100644
 --- a/.gitignore
 +++ b/.gitignore
 @@ -234,7 +234,7 @@ miniapps/meshing/mobius-strip.mesh
@@ -404,7 +404,7 @@ index 2ed9f4a1b..96373b2d1 100644
  
  ///  C = A * B  (as boolean matrices)
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index 506483992..673765dce 100644
+index 625b049c0..f492d92a1 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
 @@ -19,6 +19,7 @@
@@ -423,7 +423,7 @@ index 506483992..673765dce 100644
  
  // Include the METIS header, if using version 5. If using METIS 4, the needed
  // declarations are inlined below, i.e. no header is needed.
-@@ -1291,7 +1293,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
+@@ -1290,7 +1292,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
     return res;
  }
  
@@ -432,7 +432,7 @@ index 506483992..673765dce 100644
  {
     os << "face topology=";
     switch (info.topology)
-@@ -3150,7 +3152,7 @@ void Mesh::FinalizeTopology(bool generate_bdr)
+@@ -3123,7 +3125,7 @@ void Mesh::FinalizeTopology(bool generate_bdr)
     {
        GetElementToFaceTable();
        GenerateFaces();
@@ -441,7 +441,7 @@ index 506483992..673765dce 100644
        {
           GenerateBoundaryElements();
           GetElementToFaceTable(); // update be_to_face
-@@ -3170,7 +3172,7 @@ void Mesh::FinalizeTopology(bool generate_bdr)
+@@ -3143,7 +3145,7 @@ void Mesh::FinalizeTopology(bool generate_bdr)
        if (Dim == 2)
        {
           GenerateFaces(); // 'Faces' in 2D refers to the edges
@@ -450,7 +450,7 @@ index 506483992..673765dce 100644
           {
              GenerateBoundaryElements();
           }
-@@ -3184,7 +3186,7 @@ void Mesh::FinalizeTopology(bool generate_bdr)
+@@ -3157,7 +3159,7 @@ void Mesh::FinalizeTopology(bool generate_bdr)
     if (Dim == 1)
     {
        GenerateFaces();
@@ -459,7 +459,7 @@ index 506483992..673765dce 100644
        {
           // be_to_face will be set inside GenerateBoundaryElements
           GenerateBoundaryElements();
-@@ -6095,6 +6097,12 @@ const FiniteElementSpace *Mesh::GetNodalFESpace() const
+@@ -6071,6 +6073,12 @@ const FiniteElementSpace *Mesh::GetNodalFESpace() const
  
  void Mesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
  {
@@ -472,16 +472,7 @@ index 506483992..673765dce 100644
     space_dim = (space_dim == -1) ? spaceDim : space_dim;
     FiniteElementCollection* nfec;
     if (discont)
-@@ -11137,7 +11145,7 @@ void Mesh::Printer(std::ostream &os, std::string section_delimiter) const
-    }
- }
- 
--void Mesh::PrintTopo(std::ostream &os,const Array<int> &e_to_k) const
-+void Mesh::PrintTopo(std::ostream &os, const Array<int> &e_to_k) const
- {
-    int i;
-    Array<int> vert;
-@@ -13008,6 +13016,878 @@ void Mesh::GetGeometricParametersFromJacobian(const DenseMatrix &J,
+@@ -13068,6 +13076,878 @@ void Mesh::GetGeometricParametersFromJacobian(const DenseMatrix &J,
  }
  
  
@@ -1361,7 +1352,7 @@ index 506483992..673765dce 100644
                                     int flags, MemoryType d_mt)
  {
 diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
-index 97d3e76bb..1be7e372e 100644
+index 317e4349a..3aebaad29 100644
 --- a/mesh/mesh.hpp
 +++ b/mesh/mesh.hpp
 @@ -28,6 +28,7 @@
@@ -1401,25 +1392,25 @@ index 97d3e76bb..1be7e372e 100644
  
     // Readers for different mesh formats, used in the Load() method.
     // The implementations of these methods are in mesh_readers.cpp.
-@@ -430,7 +433,7 @@ protected:
+@@ -433,7 +436,7 @@ protected:
  
     void UpdateNURBS();
  
--   void PrintTopo(std::ostream &out, const Array<int> &e_to_k) const;
-+   void PrintTopo(std::ostream &os, const Array<int> &e_to_k) const;
+-   void PrintTopo(std::ostream &out, const Array<int> &e_to_k,
++   void PrintTopo(std::ostream &os, const Array<int> &e_to_k,
+                   const std::string &comments = "") const;
  
     /// Used in GetFaceElementTransformations (...)
-    void GetLocalPtToSegTransformation(IsoparametricTransformation &, int);
-@@ -539,7 +542,7 @@ protected:
-    // If NURBS mesh, write NURBS format. If NCMesh, write mfem v1.1 format.
-    // If section_delimiter is empty, write mfem v1.0 format. Otherwise, write
-    // mfem v1.2 format with the given section_delimiter at the end.
+@@ -545,7 +548,7 @@ protected:
+        mfem v1.2 format with the given section_delimiter at the end.
+        If @a comments is non-empty, it will be printed after the first line of
+        the file, and each line should begin with '#'. */
 -   void Printer(std::ostream &out = mfem::out,
 +   void Printer(std::ostream &os = mfem::out,
-                 std::string section_delimiter = "") const;
+                 std::string section_delimiter = "",
+                 const std::string &comments = "") const;
  
-    /// @brief Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz],
-@@ -949,7 +952,7 @@ public:
+@@ -956,7 +959,7 @@ public:
  
     int AddBdrPoint(int v, int attr = 1);
  
@@ -1428,7 +1419,7 @@ index 97d3e76bb..1be7e372e 100644
     /// Finalize the construction of a triangular Mesh.
     void FinalizeTriMesh(int generate_edges = 0, int refine = 0,
                          bool fix_orientation = true);
-@@ -2191,7 +2194,7 @@ public:
+@@ -2202,7 +2205,7 @@ public:
                                 std::ostream &os, int elem_attr = 0) const;
  
     void PrintElementsWithPartitioning (int *partitioning,
@@ -1437,7 +1428,7 @@ index 97d3e76bb..1be7e372e 100644
                                         int interior_faces = 0);
  
     /// Print set of disjoint surfaces:
-@@ -2199,13 +2202,13 @@ public:
+@@ -2210,13 +2213,13 @@ public:
      * If Aface_face(i,j) != 0, print face j as a boundary
      * element with attribute i+1.
      */
@@ -1453,7 +1444,7 @@ index 97d3e76bb..1be7e372e 100644
  
     /** @brief Compute and print mesh characteristics such as number of vertices,
         number of elements, number of boundary elements, minimal and maximal
-@@ -2225,7 +2228,7 @@ public:
+@@ -2236,7 +2239,7 @@ public:
  
  #ifdef MFEM_DEBUG
     /// Output an NCMesh-compatible debug dump.
@@ -1462,7 +1453,7 @@ index 97d3e76bb..1be7e372e 100644
  #endif
  
     /// @}
-@@ -2304,9 +2307,194 @@ public:
+@@ -2315,9 +2318,194 @@ public:
     /// @}
  };
  
@@ -1660,7 +1651,7 @@ index 97d3e76bb..1be7e372e 100644
  
  
  /** @brief Structure for storing mesh geometric factors: coordinates, Jacobians,
-@@ -2315,7 +2503,6 @@ std::ostream &operator<<(std::ostream &out, const Mesh &mesh);
+@@ -2326,7 +2514,6 @@ std::ostream &operator<<(std::ostream &out, const Mesh &mesh);
      Mesh. See Mesh::GetGeometricFactors(). */
  class GeometricFactors
  {
@@ -1668,7 +1659,7 @@ index 97d3e76bb..1be7e372e 100644
  private:
     void Compute(const GridFunction &nodes,
                  MemoryType d_mt = MemoryType::DEFAULT);
-@@ -2363,6 +2550,7 @@ public:
+@@ -2374,6 +2561,7 @@ public:
     Vector detJ;
  };
  
@@ -1676,7 +1667,7 @@ index 97d3e76bb..1be7e372e 100644
  /** @brief Structure for storing face geometric factors: coordinates, Jacobians,
      determinants of the Jacobians, and normal vectors. */
  /** Typically objects of this type are constructed and owned by objects of class
-@@ -2417,6 +2605,7 @@ public:
+@@ -2428,6 +2616,7 @@ public:
     Vector normal;
  };
  
@@ -1684,7 +1675,7 @@ index 97d3e76bb..1be7e372e 100644
  /// Class used to extrude the nodes of a mesh
  class NodeExtrudeCoefficient : public VectorCoefficient
  {
-@@ -2448,8 +2637,12 @@ inline void ShiftRight(int &a, int &b, int &c)
+@@ -2459,8 +2648,12 @@ inline void ShiftRight(int &a, int &b, int &c)
     a = c;  c = b;  b = t;
  }
  
@@ -1699,10 +1690,10 @@ index 97d3e76bb..1be7e372e 100644
  }
  
 diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
-index 72bc3b98a..c43e2ea53 100644
+index d4f6d1abd..8fa83365c 100644
 --- a/mesh/pmesh.cpp
 +++ b/mesh/pmesh.cpp
-@@ -250,6 +250,8 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
+@@ -249,6 +249,8 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
        BuildSharedVertMapping(nsvert, vert_element, vert_global_local);
        delete vert_element;
  
@@ -1711,7 +1702,7 @@ index 72bc3b98a..c43e2ea53 100644
        SetMeshGen();
        meshgen = mesh.meshgen; // copy the global 'meshgen'
     }
-@@ -1526,6 +1528,7 @@ ParMesh ParMesh::MakeSimplicial(ParMesh &orig_mesh)
+@@ -1523,6 +1525,7 @@ ParMesh ParMesh::MakeSimplicial(ParMesh &orig_mesh)
  void ParMesh::Finalize(bool refine, bool fix_orientation)
  {
     const int meshgen_save = meshgen; // Mesh::Finalize() may call SetMeshGen()
@@ -1719,7 +1710,7 @@ index 72bc3b98a..c43e2ea53 100644
  
     Mesh::Finalize(refine, fix_orientation);
  
-@@ -6213,6 +6216,7 @@ void ParMesh::ParPrint(ostream &os) const
+@@ -6249,6 +6252,7 @@ void ParMesh::ParPrint(ostream &os, const std::string &comments) const
     {
        os << "total_shared_faces " << sface_lface.Size() << '\n';
     }
@@ -1728,7 +1719,7 @@ index 72bc3b98a..c43e2ea53 100644
     {
        {
 diff --git a/mesh/pmesh.hpp b/mesh/pmesh.hpp
-index c2ab84434..a256e6ad5 100644
+index af5aa9a7c..d723204f6 100644
 --- a/mesh/pmesh.hpp
 +++ b/mesh/pmesh.hpp
 @@ -108,6 +108,8 @@ protected:
@@ -1741,10 +1732,10 @@ index c2ab84434..a256e6ad5 100644
     void FinalizeParTopo();
  
 diff --git a/mesh/tetrahedron.cpp b/mesh/tetrahedron.cpp
-index 0815fc0a7..c1b0ae6d8 100644
+index c7ebc064b..7c9c8de7a 100644
 --- a/mesh/tetrahedron.cpp
 +++ b/mesh/tetrahedron.cpp
-@@ -55,7 +55,7 @@ void Tetrahedron::Init(int ind1, int ind2, int ind3, int ind4, int attr,
+@@ -53,7 +53,7 @@ void Tetrahedron::Init(int ind1, int ind2, int ind3, int ind4, int attr,
  }
  
  void Tetrahedron::ParseRefinementFlag(int refinement_edges[2], int &type,
@@ -1753,7 +1744,7 @@ index 0815fc0a7..c1b0ae6d8 100644
  {
     int i, f = refinement_flag;
  
-@@ -136,9 +136,10 @@ void Tetrahedron::CreateRefinementFlag(int refinement_edges[2], int type,
+@@ -134,9 +134,10 @@ void Tetrahedron::CreateRefinementFlag(int refinement_edges[2], int type,
     refinement_flag |= refinement_edges[0];
  }
  
@@ -1767,7 +1758,7 @@ index 0815fc0a7..c1b0ae6d8 100644
     switch (face)
     {
 diff --git a/mesh/tetrahedron.hpp b/mesh/tetrahedron.hpp
-index ad018a037..ef8f36eb8 100644
+index c434ae903..b72e5db55 100644
 --- a/mesh/tetrahedron.hpp
 +++ b/mesh/tetrahedron.hpp
 @@ -58,12 +58,13 @@ public:
diff --git a/extern/patch/mfem/patch_mesh_vis_dev.diff b/extern/patch/mfem/patch_mesh_vis_dev.diff
index 4f9d21808..0541362af 100644
--- a/extern/patch/mfem/patch_mesh_vis_dev.diff
+++ b/extern/patch/mfem/patch_mesh_vis_dev.diff
@@ -470,10 +470,10 @@ index c216188af..4a26cd1f8 100644
     /// preserve timestep metadata for any solutions prior to the currently
     /// defined time.
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index 7e57ae2e2..87e446c50 100644
+index 102ee4de6..625b049c0 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
-@@ -11414,7 +11414,7 @@ void Mesh::PrintVTU(std::string fname,
+@@ -11423,7 +11423,7 @@ void Mesh::PrintVTU(std::string fname,
                      VTKFormat format,
                      bool high_order_output,
                      int compression_level,
@@ -482,7 +482,7 @@ index 7e57ae2e2..87e446c50 100644
  {
     int ref = (high_order_output && Nodes)
               ? Nodes->FESpace()->GetMaxElementOrder() : 1;
-@@ -11428,7 +11428,7 @@ void Mesh::PrintVTU(std::string fname,
+@@ -11437,7 +11437,7 @@ void Mesh::PrintVTU(std::string fname,
     }
     os << " byte_order=\"" << VTKByteOrder() << "\">\n";
     os << "<UnstructuredGrid>\n";
@@ -492,10 +492,10 @@ index 7e57ae2e2..87e446c50 100644
     os << "</UnstructuredGrid>\n";
     os << "</VTKFile>" << std::endl;
 diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
-index f396ff91d..40945eb04 100644
+index b82ce4a76..317e4349a 100644
 --- a/mesh/mesh.hpp
 +++ b/mesh/mesh.hpp
-@@ -2178,7 +2178,7 @@ public:
+@@ -2186,7 +2186,7 @@ public:
                           VTKFormat format=VTKFormat::ASCII,
                           bool high_order_output=false,
                           int compression_level=0,
@@ -505,10 +505,10 @@ index f396ff91d..40945eb04 100644
         boundary attributes as a data array (useful for boundary conditions). */
     void PrintBdrVTU(std::string fname,
 diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
-index 9310b7608..f93433b2c 100644
+index 81012851d..d4f6d1abd 100644
 --- a/mesh/pmesh.cpp
 +++ b/mesh/pmesh.cpp
-@@ -6299,7 +6299,7 @@ void ParMesh::PrintVTU(std::string pathname,
+@@ -6303,7 +6303,7 @@ void ParMesh::PrintVTU(std::string pathname,
                         VTKFormat format,
                         bool high_order_output,
                         int compression_level,
@@ -517,7 +517,7 @@ index 9310b7608..f93433b2c 100644
  {
     int pad_digits_rank = 6;
     DataCollection::create_directory(pathname, this, MyRank);
-@@ -6359,7 +6359,8 @@ void ParMesh::PrintVTU(std::string pathname,
+@@ -6363,7 +6363,8 @@ void ParMesh::PrintVTU(std::string pathname,
  
     std::string vtu_fname = pathname + "/" + fname + ".proc"
                             + to_padded_string(MyRank, pad_digits_rank);
@@ -528,10 +528,10 @@ index 9310b7608..f93433b2c 100644
  
  int ParMesh::FindPoints(DenseMatrix& point_mat, Array<int>& elem_id,
 diff --git a/mesh/pmesh.hpp b/mesh/pmesh.hpp
-index 66519301a..c2ab84434 100644
+index 472a0491e..af5aa9a7c 100644
 --- a/mesh/pmesh.hpp
 +++ b/mesh/pmesh.hpp
-@@ -680,7 +680,7 @@ public:
+@@ -690,7 +690,7 @@ public:
                   VTKFormat format=VTKFormat::ASCII,
                   bool high_order_output=false,
                   int compression_level=0,
diff --git a/extern/patch/mfem/patch_ncmesh_interior_boundary_dev.diff b/extern/patch/mfem/patch_ncmesh_interior_boundary_dev.diff
index 8c5173f20..7e8224102 100644
--- a/extern/patch/mfem/patch_ncmesh_interior_boundary_dev.diff
+++ b/extern/patch/mfem/patch_ncmesh_interior_boundary_dev.diff
@@ -1,5 +1,5 @@
 diff --git a/CHANGELOG b/CHANGELOG
-index e76cda016..b01affc9b 100644
+index 30a9e6e6e..b31adac9d 100644
 --- a/CHANGELOG
 +++ b/CHANGELOG
 @@ -15,6 +15,8 @@ Discretization improvements
@@ -12,10 +12,10 @@ index e76cda016..b01affc9b 100644
  - Added functionality for construction of cut-surface and cut-volume
    IntegrationRules through a moment-fitting approach. The cut is specified by
 diff --git a/fem/bilinearform.cpp b/fem/bilinearform.cpp
-index 2fee83d19..e0e7a6b27 100644
+index 194829008..2661ff7f1 100644
 --- a/fem/bilinearform.cpp
 +++ b/fem/bilinearform.cpp
-@@ -1074,7 +1074,8 @@ void BilinearForm::EliminateEssentialBCFromDofs(
+@@ -1077,7 +1077,8 @@ void BilinearForm::EliminateEssentialBCFromDofs(
  void BilinearForm::EliminateEssentialBCFromDofs (const Array<int> &ess_dofs,
                                                   DiagonalPolicy dpolicy)
  {
@@ -25,7 +25,7 @@ index 2fee83d19..e0e7a6b27 100644
  
     for (int i = 0; i < ess_dofs.Size(); i++)
        if (ess_dofs[i] < 0)
-@@ -1086,7 +1087,8 @@ void BilinearForm::EliminateEssentialBCFromDofs (const Array<int> &ess_dofs,
+@@ -1089,7 +1090,8 @@ void BilinearForm::EliminateEssentialBCFromDofs (const Array<int> &ess_dofs,
  void BilinearForm::EliminateEssentialBCFromDofsDiag (const Array<int> &ess_dofs,
                                                       double value)
  {
@@ -1029,7 +1029,7 @@ index a8186c0c8..450cac0ce 100644
  
  extern MFEM_EXPORT class TriLinear3DFiniteElement HexahedronFE;
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index 673765dce..a46277750 100644
+index a037c647c..a8ff93b24 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
 @@ -1416,6 +1416,7 @@ Geometry::Type Mesh::GetFaceGeometry(int Face) const
@@ -1224,10 +1224,10 @@ index 673765dce..a46277750 100644
        {
           case Element::TRIANGLE:
 diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
-index 1be7e372e..58e77d388 100644
+index 0c461fd42..3eee6c415 100644
 --- a/mesh/mesh.hpp
 +++ b/mesh/mesh.hpp
-@@ -597,7 +597,7 @@ protected:
+@@ -601,7 +601,7 @@ protected:
     void Make1D(int n, double sx = 1.0);
  
     /// Internal function used in Mesh::MakeRefined
@@ -1236,7 +1236,7 @@ index 1be7e372e..58e77d388 100644
                       int ref_type);
  
     /// Initialize vertices/elements/boundary/tables from a nonconforming mesh.
-@@ -1788,7 +1788,7 @@ public:
+@@ -1792,7 +1792,7 @@ public:
                            };
  
     /** @brief This structure is used as a human readable output format that
@@ -1246,7 +1246,7 @@ index 1be7e372e..58e77d388 100644
  
         The element indices in this structure don't need further processing,
 diff --git a/mesh/ncmesh.cpp b/mesh/ncmesh.cpp
-index ecb5fb90b..5c57adc4a 100644
+index 5cc8225f7..2074d3048 100644
 --- a/mesh/ncmesh.cpp
 +++ b/mesh/ncmesh.cpp
 @@ -19,6 +19,36 @@
@@ -1784,7 +1784,7 @@ index ecb5fb90b..5c57adc4a 100644
        splits[1] = splits[0];
     }
     else
-@@ -6377,17 +6454,17 @@ void NCMesh::DebugDump(std::ostream &os) const
+@@ -6381,17 +6458,17 @@ void NCMesh::DebugDump(std::ostream &os) const
  
     // dump faces
     os << faces.Size() << "\n";
@@ -1808,7 +1808,7 @@ index ecb5fb90b..5c57adc4a 100644
  
        const int* fv = GI[el.Geom()].faces[lf];
        const int nfv = GI[el.Geom()].nfv[lf];
-@@ -6397,7 +6474,7 @@ void NCMesh::DebugDump(std::ostream &os) const
+@@ -6401,7 +6478,7 @@ void NCMesh::DebugDump(std::ostream &os) const
        {
           os << " " << el.node[fv[i]];
        }
@@ -1818,7 +1818,7 @@ index ecb5fb90b..5c57adc4a 100644
     }
  }
 diff --git a/mesh/ncmesh.hpp b/mesh/ncmesh.hpp
-index 8ac50d342..b004cf43e 100644
+index f5fc20856..1e04d16bb 100644
 --- a/mesh/ncmesh.hpp
 +++ b/mesh/ncmesh.hpp
 @@ -293,6 +293,7 @@ public:
@@ -1846,7 +1846,7 @@ index 8ac50d342..b004cf43e 100644
  
     /// Return element geometry type. @a index is the Mesh element number.
     Geometry::Type GetElementGeometry(int index) const
-@@ -456,7 +459,6 @@ protected: // non-public interface for the Mesh class
+@@ -458,7 +461,6 @@ protected: // non-public interface for the Mesh class
         by calling Mesh::SetCurvature or otherwise setting the Nodes. */
     void MakeTopologyOnly() { coordinates.DeleteAll(); }
  
@@ -1854,7 +1854,7 @@ index 8ac50d342..b004cf43e 100644
  protected: // implementation
  
     int Dim, spaceDim; ///< dimensions of the elements and the vertex coordinates
-@@ -594,7 +596,6 @@ protected: // implementation
+@@ -596,7 +598,6 @@ protected: // implementation
  
     Table element_vertex; ///< leaf-element to vertex table, see FindSetNeighbors
  
@@ -1862,7 +1862,7 @@ index 8ac50d342..b004cf43e 100644
     /// Update the leaf elements indices in leaf_elements
     void UpdateLeafElements();
  
-@@ -712,10 +713,79 @@ protected: // implementation
+@@ -714,10 +715,79 @@ protected: // implementation
  
     mfem::Element* NewMeshElement(int geom) const;
  
@@ -1944,7 +1944,7 @@ index 8ac50d342..b004cf43e 100644
  
     void ForceRefinement(int vn1, int vn2, int vn3, int vn4);
  
-@@ -792,7 +862,6 @@ protected: // implementation
+@@ -794,7 +864,6 @@ protected: // implementation
     virtual void ElementSharesEdge(int elem, int local, int enode) {} // ParNCMesh
     virtual void ElementSharesVertex(int elem, int local, int vnode) {} // ParNCMesh
  
@@ -1952,7 +1952,7 @@ index 8ac50d342..b004cf43e 100644
     // neighbors / element_vertex table
  
     /** Return all vertex-, edge- and face-neighbors of a set of elements.
-@@ -981,9 +1050,7 @@ protected: // implementation
+@@ -983,9 +1052,7 @@ protected: // implementation
     void InitDerefTransforms();
     void SetDerefMatrixCodes(int parent, Array<int> &fine_coarse);
  
@@ -1962,7 +1962,7 @@ index 8ac50d342..b004cf43e 100644
     struct TmpVertex
     {
        bool valid, visited;
-@@ -1002,10 +1069,56 @@ protected: // implementation
+@@ -1004,10 +1071,56 @@ protected: // implementation
  
     void FindFaceNodes(int face, int node[4]);
  
@@ -2019,7 +2019,7 @@ index 8ac50d342..b004cf43e 100644
  
     void CountSplits(int elem, int splits[3]) const;
     void GetLimitRefinements(Array<Refinement> &refinements, int max_level);
-@@ -1042,7 +1155,6 @@ protected: // implementation
+@@ -1044,7 +1157,6 @@ protected: // implementation
     /// Load the deprecated MFEM mesh v1.1 format for backward compatibility.
     void LoadLegacyFormat(std::istream &input, int &curved, int &is_nc);
  
@@ -2028,7 +2028,7 @@ index 8ac50d342..b004cf43e 100644
  
     /// This holds in one place the constants about the geometries we support
 diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
-index c43e2ea53..5b5657546 100644
+index 96e6c43e1..5225cfc14 100644
 --- a/mesh/pmesh.cpp
 +++ b/mesh/pmesh.cpp
 @@ -364,11 +364,8 @@ int ParMesh::BuildLocalVertices(const mfem::Mesh &mesh,
diff --git a/extern/patch/mfem/patch_par_tet_mesh_fix.diff b/extern/patch/mfem/patch_par_tet_mesh_fix.diff
index bcca89d0f..3fca8d363 100644
--- a/extern/patch/mfem/patch_par_tet_mesh_fix.diff
+++ b/extern/patch/mfem/patch_par_tet_mesh_fix.diff
@@ -13,10 +13,10 @@ index c37205eb1..f1b003cae 100644
     virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const { return 0; }
  
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index 87e446c50..506483992 100644
+index f492d92a1..a037c647c 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
-@@ -483,6 +483,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
+@@ -485,6 +485,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
              Nodes->FESpace()->GetTraceElement(elem_id, face_geom);
           MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
                       "Mesh requires nodal Finite Element.");
@@ -24,7 +24,7 @@ index 87e446c50..506483992 100644
           IntegrationRule eir(face_el->GetDof());
           FaceElemTr.Loc1.Transf.ElementNo = elem_id;
           FaceElemTr.Loc1.Transf.mesh = this;
-@@ -2101,7 +2102,7 @@ void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
+@@ -2103,7 +2104,7 @@ void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
  
     if (refine)
     {
@@ -33,7 +33,7 @@ index 87e446c50..506483992 100644
     }
  
     if (generate_edges)
-@@ -2555,82 +2556,110 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
+@@ -2557,82 +2558,110 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
  }
  
  
@@ -180,7 +180,7 @@ index 87e446c50..506483992 100644
        }
     }
  }
-@@ -2980,9 +3009,7 @@ void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
+@@ -2982,9 +3011,7 @@ void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
  
     if (refine)
     {
@@ -191,7 +191,7 @@ index 87e446c50..506483992 100644
     }
  
     GetElementToFaceTable();
-@@ -3249,8 +3276,7 @@ void Mesh::Finalize(bool refine, bool fix_orientation)
+@@ -3251,8 +3278,7 @@ void Mesh::Finalize(bool refine, bool fix_orientation)
     // only perform it when Dim == spaceDim.
     if (Dim >= 2 && Dim == spaceDim)
     {
@@ -201,7 +201,7 @@ index 87e446c50..506483992 100644
        {
           MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
                       faces_info[i].Elem2Inf%2 != 0, "Invalid mesh topology."
-@@ -3950,8 +3976,6 @@ void Mesh::Make2D(int nx, int ny, Element::Type type,
+@@ -3952,8 +3978,6 @@ void Mesh::Make2D(int nx, int ny, Element::Type type,
           boundary[2*nx+j] = new Segment((j+1)*m, j*m, 4);
           boundary[2*nx+ny+j] = new Segment(j*m+nx, (j+1)*m+nx, 2);
        }
@@ -210,7 +210,7 @@ index 87e446c50..506483992 100644
     }
     else
     {
-@@ -6148,37 +6172,21 @@ static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
+@@ -6156,37 +6180,21 @@ static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
  
  int Mesh::CheckElementOrientation(bool fix_it)
  {
@@ -253,7 +253,7 @@ index 87e446c50..506483992 100644
                 switch (GetElementType(i))
                 {
                    case Element::TRIANGLE:
-@@ -6198,88 +6206,41 @@ int Mesh::CheckElementOrientation(bool fix_it)
+@@ -6206,88 +6214,41 @@ int Mesh::CheckElementOrientation(bool fix_it)
           }
        }
     }
@@ -365,7 +365,7 @@ index 87e446c50..506483992 100644
           }
        }
     }
-@@ -7205,24 +7166,12 @@ void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
+@@ -7213,24 +7174,12 @@ void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
  
     pointmat.SetSize(spaceDim, nv);
     for (k = 0; k < spaceDim; k++)
@@ -392,10 +392,10 @@ index 87e446c50..506483992 100644
  
  // static method
 diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
-index 40945eb04..97d3e76bb 100644
+index 3aebaad29..0c461fd42 100644
 --- a/mesh/mesh.hpp
 +++ b/mesh/mesh.hpp
-@@ -332,12 +332,9 @@ protected:
+@@ -335,12 +335,9 @@ protected:
     /** Also, initializes #mesh_geoms. */
     void SetMeshGen();
  
@@ -411,7 +411,7 @@ index 40945eb04..97d3e76bb 100644
  
     // Methods used to prepare and apply permutation of the mesh nodes assuming
 diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
-index f93433b2c..72bc3b98a 100644
+index 8fa83365c..96e6c43e1 100644
 --- a/mesh/pmesh.cpp
 +++ b/mesh/pmesh.cpp
 @@ -20,6 +20,7 @@
@@ -422,7 +422,7 @@ index f93433b2c..72bc3b98a 100644
  #include <iostream>
  #include <fstream>
  
-@@ -762,8 +763,10 @@ void ParMesh::BuildSharedFaceElems(int ntri_faces, int nquad_faces,
+@@ -764,8 +765,10 @@ void ParMesh::BuildSharedFaceElems(int ntri_faces, int nquad_faces,
              sface_lface[stria_counter] = lface;
              if (meshgen == 1) // Tet-only mesh
              {
@@ -435,7 +435,7 @@ index f93433b2c..72bc3b98a 100644
                 // mark the shared face for refinement by reorienting
                 // it according to the refinement flag in the tetrahedron
                 // to which this shared face belongs to.
-@@ -1739,97 +1742,59 @@ void ParMesh::GetSharedTriCommunicator(int ordering,
+@@ -1742,97 +1745,59 @@ void ParMesh::GetSharedTriCommunicator(int ordering,
  
  void ParMesh::MarkTetMeshForRefinement(const DSTable &v_to_v)
  {
@@ -568,7 +568,7 @@ index f93433b2c..72bc3b98a 100644
        }
     }
  
-@@ -1837,13 +1802,16 @@ void ParMesh::MarkTetMeshForRefinement(const DSTable &v_to_v)
+@@ -1840,13 +1805,16 @@ void ParMesh::MarkTetMeshForRefinement(const DSTable &v_to_v)
     {
        if (boundary[i]->GetType() == Element::TRIANGLE)
        {
@@ -588,7 +588,7 @@ index f93433b2c..72bc3b98a 100644
  }
  
 diff --git a/mesh/tetrahedron.cpp b/mesh/tetrahedron.cpp
-index c7ebc064b..0815fc0a7 100644
+index 7c9c8de7a..c1b0ae6d8 100644
 --- a/mesh/tetrahedron.cpp
 +++ b/mesh/tetrahedron.cpp
 @@ -13,6 +13,8 @@
@@ -600,7 +600,7 @@ index c7ebc064b..0815fc0a7 100644
  namespace mfem
  {
  
-@@ -184,19 +186,30 @@ void Tetrahedron::SetVertices(const int *ind)
+@@ -185,19 +187,30 @@ void Tetrahedron::SetVertices(const int *ind)
     }
  }
  
@@ -642,7 +642,7 @@ index c7ebc064b..0815fc0a7 100644
     {
        ind[i] = indices[i];
     }
-@@ -228,13 +241,14 @@ void Tetrahedron::MarkEdge(const DSTable &v_to_v, const int *length)
+@@ -229,13 +242,14 @@ void Tetrahedron::MarkEdge(const DSTable &v_to_v, const int *length)
     // Determine the two longest edges for the other two faces and
     // store them in ind[0] and ind[1]
     ind[0] = 2; ind[1] = 1;
@@ -663,7 +663,7 @@ index c7ebc064b..0815fc0a7 100644
  
     j = 0;
     switch (ind[0])
-@@ -345,4 +359,13 @@ Element *Tetrahedron::Duplicate(Mesh *m) const
+@@ -346,4 +360,13 @@ Element *Tetrahedron::Duplicate(Mesh *m) const
     return tet;
  }
  
@@ -678,10 +678,10 @@ index c7ebc064b..0815fc0a7 100644
 +
  }
 diff --git a/mesh/tetrahedron.hpp b/mesh/tetrahedron.hpp
-index c434ae903..ad018a037 100644
+index b72e5db55..ef8f36eb8 100644
 --- a/mesh/tetrahedron.hpp
 +++ b/mesh/tetrahedron.hpp
-@@ -76,7 +76,9 @@ public:
+@@ -77,7 +77,9 @@ public:
     /** Reorder the vertices so that the longest edge is from vertex 0
         to vertex 1. If called it should be once from the mesh constructor,
         because the order may be used later for setting the edges. **/
diff --git a/palace/CMakeLists.txt b/palace/CMakeLists.txt
index 407224c30..4ae5a7b76 100644
--- a/palace/CMakeLists.txt
+++ b/palace/CMakeLists.txt
@@ -19,7 +19,7 @@ set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
 # Initialize the project
-project(palace LANGUAGES CXX VERSION 0.11.2)
+project(palace LANGUAGES CXX VERSION 0.12.0)
 
 # Define build settings and defaults
 set(PALACE_WITH_OPENMP OFF CACHE BOOL "Use OpenMP for shared-memory parallelism")
diff --git a/test/examples/baseline b/test/examples/baseline
index e20814172..2d17a78c5 100755
--- a/test/examples/baseline
+++ b/test/examples/baseline
@@ -82,10 +82,10 @@ fi
 
 # Configure tests (use defaults set by runtests.jl otherwise)
 if [[ -n "$NUM_PROCS" ]]; then
-	export NUM_PROC_TEST=$NUM_PROCS
+    export NUM_PROC_TEST=$NUM_PROCS
 fi
 if [[ -n "$PALACE" ]]; then
-	export PALACE_TEST=$PALACE
+    export PALACE_TEST=$PALACE
 fi
 
 TEST_DIRS=(