diff --git a/cmake/ExternalGitTags.cmake b/cmake/ExternalGitTags.cmake
index 3a058ae25..f536c0de7 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
-  "28e9380ef938a46a6ac49427146fc0ea0cc9b0eb" CACHE STRING
+  "36e1f4a1b98451df0f561349f899e0f5d0e0b6bc" CACHE STRING
   "Git tag for external ButterflyPACK build"
 )
 
@@ -62,7 +62,7 @@ set(EXTERN_HYPRE_GIT_BRANCH
   "Git branch for external HYPRE build"
 )
 set(EXTERN_HYPRE_GIT_TAG
-  "5de2c5ceec9f68fa875f65bb3cfe605504c8229d" CACHE STRING
+  "3caa81955eb8d1b4e35d9b450e27cf6d07b50f6e" CACHE STRING
   "Git tag for external HYPRE build"
 )
 
@@ -76,7 +76,7 @@ set(EXTERN_LIBCEED_GIT_BRANCH
   "Git branch for external libCEED build"
 )
 set(EXTERN_LIBCEED_GIT_TAG
-  "80859cf565d67332b8cfd754474c121c2d43d208" CACHE STRING
+  "83862665467081d1f29e1ddefe74326356b818f0" CACHE STRING
   "Git tag for external libCEED build"
 )
 
@@ -90,7 +90,7 @@ set(EXTERN_LIBXSMM_GIT_BRANCH
   "Git branch for external LIBXSMM build"
 )
 set(EXTERN_LIBXSMM_GIT_TAG
-  "f5cb6f5d32a9527351dfce1b367e32055a1e759b" CACHE STRING
+  "d10cd5d59b03c04d1b0aa2aa44f28a21cd60a062" CACHE STRING
   "Git tag for external LIBXSMM build"
 )
 
@@ -104,7 +104,7 @@ set(EXTERN_MAGMA_GIT_BRANCH
   "Git branch for external MAGMA build"
 )
 set(EXTERN_MAGMA_GIT_TAG
-  "5b7fc4ee3421de71393935751112a939b718c0a4" CACHE STRING
+  "6573c5a752da571b9c5a5c5f7cd68cc24d3c2c6d" CACHE STRING
   "Git tag for external MAGMA build"
 )
 
@@ -132,7 +132,7 @@ set(EXTERN_MFEM_GIT_BRANCH
   "Git branch for external MFEM build"
 )
 set(EXTERN_MFEM_GIT_TAG
-  "e13989a2934e8895bd01f2ebbc8d78102f0c0e8f" CACHE STRING
+  "957ca0d50089795b219e4701312355f875af6019" CACHE STRING
   "Git tag for external MFEM build"
 )
 
@@ -146,7 +146,7 @@ set(EXTERN_MUMPS_GIT_BRANCH
   "Git branch for external MUMPS build"
 )
 set(EXTERN_MUMPS_GIT_TAG
-  "c997fd6236ec157e66ede97c52a296e3930db285" CACHE STRING
+  "8cfe473232bd45a660e247ea909a670ea9186ec3" CACHE STRING
   "Git tag for external MUMPS build"
 )
 
@@ -174,7 +174,7 @@ set(EXTERN_PETSC_GIT_BRANCH
   "Git branch for external PETSc build"
 )
 set(EXTERN_PETSC_GIT_TAG
-  "96e7d5cfaa07a5d3427e87ca13c8795e24a0e757" CACHE STRING
+  "4bbe91384854fec5c16b677958fdcd4ec4bdaf6b" CACHE STRING
   "Git tag for external PETSc build"
 )
 
@@ -188,7 +188,7 @@ set(EXTERN_SCALAPACK_GIT_BRANCH
   "Git branch for external ScaLAPACK build"
 )
 set(EXTERN_SCALAPACK_GIT_TAG
-  "0234af94c6578c53ac4c19f2925eb6e5c4ad6f0f" CACHE STRING
+  "25935e1a7e022ede9fd71bd86dcbaa7a3f1846b7" CACHE STRING
   "Git tag for external ScaLAPACK build"
 )
 
@@ -202,7 +202,7 @@ set(EXTERN_SLEPC_GIT_BRANCH
   "Git branch for external SLEPc build"
 )
 set(EXTERN_SLEPC_GIT_TAG
-  "ffa1d78448240d13360ff0555c46d4e767792d2c" CACHE STRING
+  "e93850835cc2843553d5b759b2277856d10dbc97" CACHE STRING
   "Git tag for external SLEPc build"
 )
 
@@ -216,7 +216,7 @@ set(EXTERN_STRUMPACK_GIT_BRANCH
   "Git branch for external STRUMPACK build"
 )
 set(EXTERN_STRUMPACK_GIT_TAG
-  "cf768ce0977904f2991ea9844e5164dd39a15102" CACHE STRING
+  "b964011b77e1cf6bb087a6e049f61c1c77a8e256" CACHE STRING
   "Git tag for external STRUMPACK build"
 )
 
@@ -230,7 +230,7 @@ set(EXTERN_SUPERLU_GIT_BRANCH
   "Git branch for external SuperLU_DIST build"
 )
 set(EXTERN_SUPERLU_GIT_TAG
-  "ea4d47b206387a1592eea46493519d57cf3984d6" CACHE STRING
+  "a66f013e468524b0f1508447fa2c9e9ff605c3ec" CACHE STRING
   "Git tag for external SuperLU_DIST build"
 )
 
@@ -244,7 +244,7 @@ set(EXTERN_ZFP_GIT_BRANCH
   "Git branch for external ZFP build"
 )
 set(EXTERN_ZFP_GIT_TAG
-  "5c976d8da013988174f931845862b6f94119cade" CACHE STRING
+  "3b6eb55c775e8c22d89c2a584160e3fe81cf9f5a" CACHE STRING
   "Git tag for external ZFP build"
 )
 
diff --git a/extern/patch/mfem/patch_gmsh_parser_performance.diff b/extern/patch/mfem/patch_gmsh_parser_performance.diff
index 716d24e57..9a64aef55 100644
--- a/extern/patch/mfem/patch_gmsh_parser_performance.diff
+++ b/extern/patch/mfem/patch_gmsh_parser_performance.diff
@@ -1,8 +1,8 @@
 diff --git a/mesh/mesh_readers.cpp b/mesh/mesh_readers.cpp
-index c5f01b50d..645b7f167 100644
+index 377cc906e..cf06e5456 100644
 --- a/mesh/mesh_readers.cpp
 +++ b/mesh/mesh_readers.cpp
-@@ -1535,7 +1535,7 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -1548,7 +1548,7 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
     // A map between a serial number of the vertex and its number in the file
     // (there may be gaps in the numbering, and also Gmsh enumerates vertices
     // starting from 1, not 0)
@@ -11,7 +11,7 @@ index c5f01b50d..645b7f167 100644
  
     // A map containing names of physical curves, surfaces, and volumes.
     // The first index is the dimension of the physical manifold, the second
-@@ -1553,9 +1553,10 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -1566,9 +1566,10 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
     // is non-trivial. Note that with these assumptions a 2D mesh parallel to the
     // yz plane will be considered a surface mesh embedded in 3D whereas the same
     // 2D mesh parallel to the xy plane will be considered a 2D mesh.
@@ -25,7 +25,7 @@ index c5f01b50d..645b7f167 100644
  
     // Mesh order
     int mesh_order = 1;
-@@ -1576,7 +1577,6 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -1589,7 +1590,6 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
           getline(input, buff);
           vertices.SetSize(NumOfVertices);
           int serial_number;
@@ -33,7 +33,7 @@ index c5f01b50d..645b7f167 100644
           real_t coord[gmsh_dim];
           for (int ver = 0; ver < NumOfVertices; ++ver)
           {
-@@ -1638,7 +1638,7 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -1651,7 +1651,7 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
  
           // number of nodes for each type of Gmsh elements, type is the index of
           // the array + 1
@@ -42,7 +42,7 @@ index c5f01b50d..645b7f167 100644
           {
              2,   // 2-node line.
              3,   // 3-node triangle.
-@@ -1848,97 +1848,324 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -1861,97 +1861,324 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
  
               This corresponds to the quad9 mapping below.
           */
@@ -443,7 +443,7 @@ index c5f01b50d..645b7f167 100644
              while (n_elem_part < num_of_all_elements)
              {
                 input.read(reinterpret_cast<char*>(header),
-@@ -1946,11 +2173,10 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -1959,11 +2186,10 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                 type_of_element = header[0];
                 n_elem_one_type = header[1];
                 n_tags          = header[2];
@@ -456,7 +456,7 @@ index c5f01b50d..645b7f167 100644
                 for (int el = 0; el < n_elem_one_type; ++el)
                 {
                    input.read(reinterpret_cast<char*>(&data[0]),
-@@ -1968,218 +2194,26 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -1981,218 +2207,26 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                    n_partitions = (n_tags > 2) ? data[dd++] : 0;
                    // we currently just skip the partitions if they exist, and go
                    // directly to vertices describing the mesh element
@@ -680,7 +680,7 @@ index c5f01b50d..645b7f167 100644
                 for (int i = 0; i < n_tags; ++i) { input >> data[i]; }
                 // physical domain - the most important value (to distinguish
                 // materials with different properties)
-@@ -2193,210 +2227,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -2206,210 +2240,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                 // we currently just skip the partitions if they exist, and go
                 // directly to vertices describing the mesh element
                 const int n_elem_nodes = nodes_of_gmsh_element[type_of_element-1];
@@ -896,7 +896,7 @@ index c5f01b50d..645b7f167 100644
              } // el (all elements)
           } // if ASCII
  
-@@ -2520,81 +2363,106 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -2533,81 +2376,106 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
              Nodes_gf.SetSpace(nfes);
              Nodes_gf.MakeOwner(nfec);
  
@@ -1026,7 +1026,7 @@ index c5f01b50d..645b7f167 100644
                       }
                       vm = ho_pyr[el_order];
                       break;
-@@ -2602,62 +2470,47 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -2615,62 +2483,47 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                       MFEM_WARNING("Unsupported Gmsh element type.");
                       break;
                 }
@@ -1121,7 +1121,7 @@ index c5f01b50d..645b7f167 100644
           }
  
           // Suppress warnings (MFEM_CONTRACT_VAR does not work here with nvcc):
-@@ -2795,7 +2648,7 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
+@@ -2808,7 +2661,7 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                 v[j] = v2v[v[j]];
              }
           }
diff --git a/extern/patch/mfem/patch_mesh_prism_vtu_fix.diff b/extern/patch/mfem/patch_mesh_prism_vtu_fix.diff
index 1595722b5..ad3f1314a 100644
--- a/extern/patch/mfem/patch_mesh_prism_vtu_fix.diff
+++ b/extern/patch/mfem/patch_mesh_prism_vtu_fix.diff
@@ -1,8 +1,8 @@
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index 8b32f1eba..e18def619 100644
+index 2c6f144e9..a32ba38e8 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
-@@ -11811,7 +11811,7 @@ void Mesh::PrintVTU(std::ostream &os, int ref, VTKFormat format,
+@@ -11856,7 +11856,7 @@ void Mesh::PrintVTU(std::ostream &os, int ref, VTKFormat format,
              const int *p = VTKGeometry::VertexPermutation[geom];
              for (int k = 0; k < nv; k++, j++)
              {
diff --git a/extern/patch/mfem/patch_mesh_vis_dev.diff b/extern/patch/mfem/patch_mesh_vis_dev.diff
index b4138fe76..73eb3d993 100644
--- a/extern/patch/mfem/patch_mesh_vis_dev.diff
+++ b/extern/patch/mfem/patch_mesh_vis_dev.diff
@@ -1,5 +1,5 @@
 diff --git a/fem/datacollection.cpp b/fem/datacollection.cpp
-index d73d2bb5b..6b3139eb4 100644
+index 1313d2a1e..735f85606 100644
 --- a/fem/datacollection.cpp
 +++ b/fem/datacollection.cpp
 @@ -210,6 +210,11 @@ void DataCollection::Save()
@@ -26,17 +26,18 @@ index d73d2bb5b..6b3139eb4 100644
     if (it != q_field_map.end())
     {
        SaveOneQField(it);
-@@ -765,7 +770,8 @@ ParaViewDataCollection::ParaViewDataCollection(const std::string&
+@@ -765,7 +770,9 @@ ParaViewDataCollection::ParaViewDataCollection(const std::string&
       levels_of_detail(1),
       pv_data_format(VTKFormat::BINARY),
       high_order_output(false),
 -     restart_mode(false)
 +     restart_mode(false),
-+     bdr_output(false)
++     bdr_output(false),
++     length_scale(1.0)
  {
     cycle = 0; // always include a valid cycle index in file names
  
-@@ -910,16 +916,19 @@ void ParaViewDataCollection::Save()
+@@ -910,16 +917,19 @@ void ParaViewDataCollection::Save()
     std::string vtu_prefix = col_path + "/" + GenerateVTUPath() + "/";
  
     // Save the local part of the mesh and grid functions fields to the local
@@ -58,7 +59,7 @@ index d73d2bb5b..6b3139eb4 100644
        const std::string &field_name = qfield.first;
        std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
        qfield.second->SaveVTU(os, pv_data_format, GetCompressionLevel(), field_name);
-@@ -935,7 +944,7 @@ void ParaViewDataCollection::Save()
+@@ -935,7 +945,7 @@ void ParaViewDataCollection::Save()
           std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
           WritePVTUHeader(pvtu_out);
  
@@ -67,8 +68,8 @@ index d73d2bb5b..6b3139eb4 100644
           pvtu_out << "<PPointData>\n";
           for (auto &field_it : field_map)
           {
-@@ -945,7 +954,24 @@ void ParaViewDataCollection::Save()
-                      << "\" NumberOfComponents=\"" << vec_dim << "\" "
+@@ -946,7 +956,24 @@ void ParaViewDataCollection::Save()
+                      << VTKComponentLabels(vec_dim) << " "
                       << "format=\"" << GetDataFormatString() << "\" />\n";
           }
 +         for (auto &field_it : coeff_field_map)
@@ -92,17 +93,17 @@ index d73d2bb5b..6b3139eb4 100644
           // Element attributes
           pvtu_out << "<PCellData>\n";
           pvtu_out << "\t<PDataArray type=\"Int32\" Name=\"" << "attribute"
-@@ -1042,7 +1068,8 @@ void ParaViewDataCollection::SaveDataVTU(std::ostream &os, int ref)
+@@ -1044,7 +1071,8 @@ void ParaViewDataCollection::SaveDataVTU(std::ostream &os, int ref)
     }
     os << " version=\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
     os << "<UnstructuredGrid>\n";
 -   mesh->PrintVTU(os,ref,pv_data_format,high_order_output,GetCompressionLevel());
 +   mesh->PrintVTU(os,ref,pv_data_format,high_order_output,GetCompressionLevel(),
-+                  bdr_output);
++                  bdr_output,length_scale);
  
     // dump out the grid functions as point data
     os << "<PointData >\n";
-@@ -1050,8 +1077,23 @@ void ParaViewDataCollection::SaveDataVTU(std::ostream &os, int ref)
+@@ -1052,8 +1080,23 @@ void ParaViewDataCollection::SaveDataVTU(std::ostream &os, int ref)
     // iterate over all grid functions
     for (FieldMapIterator it=field_map.begin(); it!=field_map.end(); ++it)
     {
@@ -126,15 +127,15 @@ index d73d2bb5b..6b3139eb4 100644
     os << "</PointData>\n";
     // close the mesh
     os << "</Piece>\n"; // close the piece open in the PrintVTU method
-@@ -1073,7 +1115,6 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
-       << " format=\"" << GetDataFormatString() << "\" >" << '\n';
+@@ -1076,7 +1119,6 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
+       << "format=\"" << GetDataFormatString() << "\" >" << '\n';
     if (vec_dim == 1)
     {
 -      // scalar data
        for (int i = 0; i < mesh->GetNE(); i++)
        {
           RefG = GlobGeometryRefiner.Refine(
-@@ -1103,11 +1144,131 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
+@@ -1106,11 +1148,131 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
           }
        }
     }
@@ -144,19 +145,22 @@ index d73d2bb5b..6b3139eb4 100644
 +   }
 +   os << "</DataArray>" << std::endl;
 +}
-+
+ 
+-   if (IsBinaryFormat())
 +void ParaViewDataCollection::SaveCoeffFieldVTU(std::ostream &os, int ref_,
 +                                               const CoeffFieldMapIterator &it)
 +{
 +   RefinedGeometry *RefG;
-+   real_t val;
++   double val;
 +   std::vector<char> buf;
 +   int vec_dim = 1;
 +   os << "<DataArray type=\"" << GetDataTypeString()
 +      << "\" Name=\"" << it->first
 +      << "\" NumberOfComponents=\"" << vec_dim << "\""
 +      << " format=\"" << GetDataFormatString() << "\" >" << '\n';
-+   {
+    {
+-      WriteVTKEncodedCompressed(os,buf.data(),buf.size(),GetCompressionLevel());
+-      os << '\n';
 +      // scalar data
 +      if (!bdr_output)
 +      {
@@ -197,12 +201,11 @@ index d73d2bb5b..6b3139eb4 100644
 +   }
 +   if (pv_data_format != VTKFormat::ASCII)
 +   {
-+      WriteBase64WithSizeAndClear(os, buf, GetCompressionLevel());
++      WriteBase64WithSizeAndClear(os, buf, compression);
 +   }
 +   os << "</DataArray>" << std::endl;
 +}
- 
--   if (IsBinaryFormat())
++
 +void ParaViewDataCollection::SaveVCoeffFieldVTU(std::ostream &os, int ref_,
 +                                                const VCoeffFieldMapIterator &it)
 +{
@@ -214,9 +217,7 @@ index d73d2bb5b..6b3139eb4 100644
 +      << "\" Name=\"" << it->first
 +      << "\" NumberOfComponents=\"" << vec_dim << "\""
 +      << " format=\"" << GetDataFormatString() << "\" >" << '\n';
-    {
--      WriteVTKEncodedCompressed(os,buf.data(),buf.size(),GetCompressionLevel());
--      os << '\n';
++   {
 +      // vector data
 +      if (!bdr_output)
 +      {
@@ -269,7 +270,7 @@ index d73d2bb5b..6b3139eb4 100644
     }
     os << "</DataArray>" << std::endl;
  }
-@@ -1140,6 +1301,11 @@ void ParaViewDataCollection::SetCompression(bool compression_)
+@@ -1143,6 +1305,16 @@ void ParaViewDataCollection::SetCompression(bool compression_)
     compression = compression_;
  }
  
@@ -277,12 +278,17 @@ index d73d2bb5b..6b3139eb4 100644
 +{
 +   bdr_output = bdr_output_;
 +}
++
++void ParaViewDataCollection::SetLengthScale(double length_scale_)
++{
++   length_scale = length_scale_;
++}
 +
  void ParaViewDataCollection::UseRestartMode(bool restart_mode_)
  {
     restart_mode = restart_mode_;
 diff --git a/fem/datacollection.hpp b/fem/datacollection.hpp
-index c0ed950e7..68bb678ba 100644
+index c0ed950e7..9ab719650 100644
 --- a/fem/datacollection.hpp
 +++ b/fem/datacollection.hpp
 @@ -133,6 +133,10 @@ private:
@@ -432,11 +438,12 @@ index c0ed950e7..68bb678ba 100644
  
     /// Load the collection. Not implemented in the base class DataCollection.
     virtual void Load(int cycle_ = 0);
-@@ -512,12 +578,17 @@ private:
+@@ -512,12 +578,18 @@ private:
     VTKFormat pv_data_format;
     bool high_order_output;
     bool restart_mode;
 +   bool bdr_output;
++   double length_scale;
  
  protected:
     void WritePVTUHeader(std::ostream &out);
@@ -450,7 +457,7 @@ index c0ed950e7..68bb678ba 100644
     const char *GetDataFormatString() const;
     const char *GetDataTypeString() const;
     /// @brief If compression is enabled, return the compression level, otherwise
-@@ -531,7 +602,6 @@ protected:
+@@ -531,7 +603,6 @@ protected:
     std::string GeneratePVTUFileName(const std::string &prefix);
     std::string GeneratePVTUPath();
  
@@ -458,85 +465,151 @@ index c0ed950e7..68bb678ba 100644
  public:
     /// Constructor. The collection name is used when saving the data.
     /** If @a mesh_ is NULL, then the mesh can be set later by calling SetMesh().
-@@ -582,6 +652,10 @@ public:
+@@ -582,6 +653,14 @@ public:
     /// by default). Reading high-order data requires ParaView 5.5 or later.
     void SetHighOrderOutput(bool high_order_output_);
  
 +   /// Configures collection to save only fields evaluated on boundaries of
 +   /// the mesh.
 +   void SetBoundaryOutput(bool bdr_output_);
++
++   /// Sets length scale used to scale mesh point coordinates on output.
++   /// The default if unset is 1.0 (unscaled).
++   void SetLengthScale(double length_scale_);
 +
     /// Enable or disable restart mode. If restart is enabled, new writes will
     /// preserve timestep metadata for any solutions prior to the currently
     /// defined time.
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index d59191f90..8b32f1eba 100644
+index 457c4bec1..2c6f144e9 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
-@@ -11662,7 +11662,7 @@ void Mesh::PrintVTU(std::string fname,
+@@ -11698,7 +11698,8 @@ void Mesh::PrintVTU(std::string fname,
                      VTKFormat format,
                      bool high_order_output,
                      int compression_level,
 -                    bool bdr)
-+                    bool bdr_elements)
++                    bool bdr_elements,
++                    double scale_factor)
  {
     int ref = (high_order_output && Nodes)
               ? Nodes->FESpace()->GetMaxElementOrder() : 1;
-@@ -11676,7 +11676,7 @@ void Mesh::PrintVTU(std::string fname,
+@@ -11712,7 +11713,8 @@ void Mesh::PrintVTU(std::string fname,
     }
     os << " byte_order=\"" << VTKByteOrder() << "\">\n";
     os << "<UnstructuredGrid>\n";
 -   PrintVTU(os, ref, format, high_order_output, compression_level, bdr);
-+   PrintVTU(os, ref, format, high_order_output, compression_level, bdr_elements);
++   PrintVTU(os, ref, format, high_order_output, compression_level, bdr_elements,
++            scale_factor);
     os << "</Piece>\n"; // need to close the piece open in the PrintVTU method
     os << "</UnstructuredGrid>\n";
     os << "</VTKFile>" << std::endl;
+@@ -11723,14 +11725,16 @@ void Mesh::PrintVTU(std::string fname,
+ void Mesh::PrintBdrVTU(std::string fname,
+                        VTKFormat format,
+                        bool high_order_output,
+-                       int compression_level)
++                       int compression_level,
++                       double scale_factor)
+ {
+-   PrintVTU(fname, format, high_order_output, compression_level, true);
++   PrintVTU(fname, format, high_order_output, compression_level, true,
++            scale_factor);
+ }
+ 
+ void Mesh::PrintVTU(std::ostream &os, int ref, VTKFormat format,
+                     bool high_order_output, int compression_level,
+-                    bool bdr_elements)
++                    bool bdr_elements, double scale_factor)
+ {
+    RefinedGeometry *RefG;
+    DenseMatrix pmat;
+@@ -11777,6 +11781,11 @@ void Mesh::PrintVTU(std::ostream &os, int ref, VTKFormat format,
+          GetElementTransformation(i)->Transform(RefG->RefPts, pmat);
+       }
+ 
++      if (scale_factor != 1.0)
++      {
++         pmat *= scale_factor;
++      }
++
+       for (int j = 0; j < pmat.Width(); j++)
+       {
+          WriteBinaryOrASCII(os, buf, pmat(0,j), " ", format);
 diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
-index 1114794e0..f328f6c29 100644
+index 297472979..73cf6124f 100644
 --- a/mesh/mesh.hpp
 +++ b/mesh/mesh.hpp
-@@ -2321,7 +2321,7 @@ public:
+@@ -2345,25 +2345,30 @@ public:
+    /** Print the mesh in VTU format. The parameter ref > 0 specifies an element
+        subdivision number (useful for high order fields and curved meshes).
+        If @a bdr_elements is true, then output (only) the boundary elements,
+-       otherwise output only the non-boundary elements. */
++       otherwise output only the non-boundary elements. The parameter
++       @a scale_factor sets a scale factor for the output mesh point
++       coordinates. */
+    void PrintVTU(std::ostream &os,
+                  int ref=1,
+                  VTKFormat format=VTKFormat::ASCII,
+                  bool high_order_output=false,
+                  int compression_level=0,
+-                 bool bdr_elements=false);
++                 bool bdr_elements=false,
++                 double scale_factor=1.0);
+    /** Print the mesh in VTU format with file name fname. */
+    virtual void PrintVTU(std::string fname,
                           VTKFormat format=VTKFormat::ASCII,
                           bool high_order_output=false,
                           int compression_level=0,
 -                         bool bdr=false);
-+                         bool bdr_elements=false);
++                         bool bdr_elements=false,
++                         double scale_factor=1.0);
     /** Print the boundary elements of the mesh in VTU format, and output the
         boundary attributes as a data array (useful for boundary conditions). */
     void PrintBdrVTU(std::string fname,
+                     VTKFormat format=VTKFormat::ASCII,
+                     bool high_order_output=false,
+-                    int compression_level=0);
++                    int compression_level=0,
++                    double scale_factor=1.0);
+ 
+ #ifdef MFEM_USE_NETCDF
+    /// @brief Export a mesh to an Exodus II file.
 diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
-index b703d83e1..46ceb35e2 100644
+index 5f5d9d38d..cc673f78e 100644
 --- a/mesh/pmesh.cpp
 +++ b/mesh/pmesh.cpp
-@@ -6399,7 +6399,7 @@ void ParMesh::PrintVTU(std::string pathname,
+@@ -6398,7 +6398,8 @@ void ParMesh::PrintVTU(std::string pathname,
                         VTKFormat format,
                         bool high_order_output,
                         int compression_level,
 -                       bool bdr)
-+                       bool bdr_elements)
++                       bool bdr_elements,
++                       double scale_factor)
  {
     int pad_digits_rank = 6;
     DataCollection::create_directory(pathname, this, MyRank);
-@@ -6459,7 +6459,8 @@ void ParMesh::PrintVTU(std::string pathname,
+@@ -6458,7 +6459,8 @@ void ParMesh::PrintVTU(std::string pathname,
  
     std::string vtu_fname = pathname + "/" + fname + ".proc"
                             + to_padded_string(MyRank, pad_digits_rank);
 -   Mesh::PrintVTU(vtu_fname, format, high_order_output, compression_level, bdr);
 +   Mesh::PrintVTU(vtu_fname, format, high_order_output, compression_level,
-+                  bdr_elements);
++                  bdr_elements, scale_factor);
  }
  
  int ParMesh::FindPoints(DenseMatrix& point_mat, Array<int>& elem_id,
 diff --git a/mesh/pmesh.hpp b/mesh/pmesh.hpp
-index d85c4c318..e33db8e43 100644
+index c766e37e9..f8a7413c2 100644
 --- a/mesh/pmesh.hpp
 +++ b/mesh/pmesh.hpp
-@@ -731,7 +731,7 @@ public:
+@@ -729,7 +729,8 @@ public:
                   VTKFormat format=VTKFormat::ASCII,
                   bool high_order_output=false,
                   int compression_level=0,
 -                 bool bdr=false) override;
-+                 bool bdr_elements=false) override;
++                 bool bdr_elements=false,
++                 double scale_factor=1.0) override;
  
     /// Parallel version of Mesh::Load().
     void Load(std::istream &input, int generate_edges = 0,
diff --git a/extern/patch/mfem/patch_par_tet_mesh_fix_dev.diff b/extern/patch/mfem/patch_par_tet_mesh_fix_dev.diff
index 08f0b0516..19de8c2ae 100644
--- a/extern/patch/mfem/patch_par_tet_mesh_fix_dev.diff
+++ b/extern/patch/mfem/patch_par_tet_mesh_fix_dev.diff
@@ -13,10 +13,10 @@ index e9dd07217..4129789bb 100644
     virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const { return 0; }
  
 diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
-index e18def619..cc54e78da 100644
+index a32ba38e8..456f9618c 100644
 --- a/mesh/mesh.cpp
 +++ b/mesh/mesh.cpp
-@@ -2152,7 +2152,7 @@ void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
+@@ -2153,7 +2153,7 @@ void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
  
     if (refine)
     {
@@ -25,7 +25,7 @@ index e18def619..cc54e78da 100644
     }
  
     if (generate_edges)
-@@ -2607,82 +2607,113 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
+@@ -2608,82 +2608,113 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
  }
  
  
@@ -177,7 +177,7 @@ index e18def619..cc54e78da 100644
        }
     }
  }
-@@ -3033,9 +3064,7 @@ void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
+@@ -3034,9 +3065,7 @@ void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
  
     if (refine)
     {
@@ -188,7 +188,7 @@ index e18def619..cc54e78da 100644
     }
  
     GetElementToFaceTable();
-@@ -3302,8 +3331,7 @@ void Mesh::Finalize(bool refine, bool fix_orientation)
+@@ -3303,8 +3332,7 @@ void Mesh::Finalize(bool refine, bool fix_orientation)
     // only perform it when Dim == spaceDim.
     if (Dim >= 2 && Dim == spaceDim)
     {
@@ -198,7 +198,7 @@ index e18def619..cc54e78da 100644
        {
           MFEM_VERIFY(faces_info[i].Elem2No < 0 ||
                       faces_info[i].Elem2Inf%2 != 0, "Invalid mesh topology."
-@@ -4003,8 +4031,6 @@ void Mesh::Make2D(int nx, int ny, Element::Type type,
+@@ -4004,8 +4032,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);
        }
@@ -207,7 +207,7 @@ index e18def619..cc54e78da 100644
     }
     else
     {
-@@ -6293,37 +6319,21 @@ static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
+@@ -6295,37 +6321,21 @@ static const char *fixed_or_not[] = { "fixed", "NOT FIXED" };
  
  int Mesh::CheckElementOrientation(bool fix_it)
  {
@@ -250,7 +250,7 @@ index e18def619..cc54e78da 100644
                 switch (GetElementType(i))
                 {
                    case Element::TRIANGLE:
-@@ -6343,88 +6353,41 @@ int Mesh::CheckElementOrientation(bool fix_it)
+@@ -6345,88 +6355,41 @@ int Mesh::CheckElementOrientation(bool fix_it)
           }
        }
     }
@@ -362,7 +362,7 @@ index e18def619..cc54e78da 100644
           }
        }
     }
-@@ -7351,24 +7314,12 @@ void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
+@@ -7384,24 +7347,12 @@ void Mesh::GetBdrPointMatrix(int i,DenseMatrix &pointmat) const
  
     pointmat.SetSize(spaceDim, nv);
     for (k = 0; k < spaceDim; k++)
@@ -389,10 +389,10 @@ index e18def619..cc54e78da 100644
  
  // static method
 diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
-index f328f6c29..c8136833c 100644
+index 73cf6124f..d737768b0 100644
 --- a/mesh/mesh.hpp
 +++ b/mesh/mesh.hpp
-@@ -347,12 +347,12 @@ protected:
+@@ -378,12 +378,12 @@ protected:
     /** Also, initializes #mesh_geoms. */
     void SetMeshGen();
  
@@ -410,7 +410,7 @@ index f328f6c29..c8136833c 100644
  
     // Methods used to prepare and apply permutation of the mesh nodes assuming
 diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
-index a46d9e315..78c224753 100644
+index cc673f78e..0e9274ece 100644
 --- a/mesh/pmesh.cpp
 +++ b/mesh/pmesh.cpp
 @@ -20,6 +20,7 @@
@@ -421,7 +421,7 @@ index a46d9e315..78c224753 100644
  #include <iostream>
  #include <fstream>
  
-@@ -757,8 +758,10 @@ void ParMesh::BuildSharedFaceElems(int ntri_faces, int nquad_faces,
+@@ -762,8 +763,10 @@ void ParMesh::BuildSharedFaceElems(int ntri_faces, int nquad_faces,
              sface_lface[stria_counter] = lface;
              if (meshgen == 1) // Tet-only mesh
              {
@@ -434,7 +434,7 @@ index a46d9e315..78c224753 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.
-@@ -1735,98 +1738,63 @@ void ParMesh::GetSharedTriCommunicator(int ordering,
+@@ -1740,98 +1743,63 @@ void ParMesh::GetSharedTriCommunicator(int ordering,
  
  void ParMesh::MarkTetMeshForRefinement(const DSTable &v_to_v)
  {
@@ -573,7 +573,7 @@ index a46d9e315..78c224753 100644
        }
     }
  
-@@ -1834,13 +1802,16 @@ void ParMesh::MarkTetMeshForRefinement(const DSTable &v_to_v)
+@@ -1839,13 +1807,16 @@ void ParMesh::MarkTetMeshForRefinement(const DSTable &v_to_v)
     {
        if (boundary[i]->GetType() == Element::TRIANGLE)
        {
@@ -592,7 +592,7 @@ index a46d9e315..78c224753 100644
     }
  }
  
-@@ -3244,7 +3215,7 @@ void ParMesh::ReorientTetMesh()
+@@ -3249,7 +3220,7 @@ void ParMesh::ReorientTetMesh()
     Array<int> svert_master_index(svert_lvert);
     for (int i = 0; i < group_svert.Size(); i++)
     {
diff --git a/extern/patch/mumps/patch_build.diff b/extern/patch/mumps/patch_build.diff
index db6145998..6f9bd940c 100644
--- a/extern/patch/mumps/patch_build.diff
+++ b/extern/patch/mumps/patch_build.diff
@@ -1,8 +1,8 @@
 diff --git a/CMakeLists.txt b/CMakeLists.txt
-index 637122a..9306afa 100644
+index 3be37f8..193e8b5 100644
 --- a/CMakeLists.txt
 +++ b/CMakeLists.txt
-@@ -43,10 +43,7 @@ find_package(Threads)
+@@ -48,10 +48,7 @@ find_package(Threads)
  
  if(MUMPS_parallel)
    find_package(MPI COMPONENTS C Fortran REQUIRED)
@@ -14,7 +14,7 @@ index 637122a..9306afa 100644
    set(NUMERIC_LIBS LAPACK::LAPACK)
  
    if(scalapack)
-@@ -94,15 +91,6 @@ endif()
+@@ -99,15 +96,6 @@ endif()
  
  list(APPEND ORDERING_LIBS pord)
  
@@ -539,7 +539,7 @@ index 8647696..0000000
 -
 -mark_as_advanced(LAPACK_LIBRARY LAPACK_INCLUDE_DIR)
 diff --git a/cmake/FindSCALAPACK.cmake b/cmake/FindSCALAPACK.cmake
-index f8c4cdc..718a510 100644
+index 9b2e056..166d05f 100644
 --- a/cmake/FindSCALAPACK.cmake
 +++ b/cmake/FindSCALAPACK.cmake
 @@ -52,8 +52,6 @@ References
@@ -668,12 +668,12 @@ index 3f0997e..14ebbe6 100644
 -find_package(LAPACK REQUIRED COMPONENTS ${LAPACK_VENDOR})
 +find_package(LAPACK REQUIRED)
 diff --git a/cmake/scalapack.cmake b/cmake/scalapack.cmake
-index f4d1a70..7238c21 100644
+index 21503e7..78c9dda 100644
 --- a/cmake/scalapack.cmake
 +++ b/cmake/scalapack.cmake
 @@ -3,27 +3,7 @@ include(GNUInstallDirs)
  
- if(find)
+ if(find AND NOT TARGET SCALAPACK::SCALAPACK)
  
 -if(NOT DEFINED SCALAPACK_VENDOR)
 -  if(LAPACK_VENDOR MATCHES "^MKL")