Modify ex26 with anisotropic diffusion and ILU smoothing

Variable order transfer operator

.gitignore

@@ -259,6 +259,8 @@ miniapps/performance/sol.*
 
 miniapps/shifted/distance
 miniapps/shifted/ParaViewDistance
+miniapps/shifted/extrapolate
+miniapps/shifted/ParaViewExtrapolate
 miniapps/shifted/diffusion
 miniapps/shifted/diffusion.mesh
 miniapps/shifted/diffusion.gf

CHANGELOG

@@ -21,6 +21,9 @@ Version 4.3.1 (development)
 - More explicit and consistent formating of the output of iterative solvers
   with the new IterativeSolver::PrintLevel options. See linalg/solvers.hpp.
 
+- Added a miniapp for PDE-based extrapolation of finite element functions. See
+  miniapps/shifted/extrapolate.cpp.
+
 - Added support for automatic differentiation. Users can select between native
   implementation and external library implementation during configuration. One
   parallel and two serial examples are implemented in the miniapps/autodiff/
@@ -95,6 +98,9 @@ Version 4.3.1 (development)
 - The HPC versions of ex1 and ex1p (in miniapps/performance) now support
   runtime selection of either 2D or 3D meshes.
 
+- Added ParaView visualization of `QuadratureFunction` fields, through both
+  `QuadratureFunction::SaveVTU` and `ParaViewDataCollection::RegisterQField`.
+
 
 Version 4.3, released on July 29, 2021
 ======================================

examples/ex26.cpp

@@ -34,6 +34,43 @@
 using namespace std;
 using namespace mfem;
 
+MatrixConstantCoefficient AnisotropicCoefficient(int dim, double anisotropy)
+{
+   DenseMatrix coeff(dim, dim);
+   coeff = 0.0;
+   coeff(0,0) = anisotropy;
+   for (int d = 1; d < dim; ++d)
+   {
+      coeff(d, d) = 1.0;
+   }
+   return coeff;
+}
+
+class SymmetricILUSmoother : public Solver
+{
+   BlockILU ilu;
+   double alpha;
+
+public:
+   SymmetricILUSmoother(Operator &op, double alpha_)
+      : ilu(op), alpha(alpha_)
+   { }
+
+   void Mult(const Vector &b, Vector &x) const
+   {
+      ilu.Mult(b, x);
+      x *= alpha;
+   }
+
+   void MultTranspose(const Vector &b, Vector &x) const
+   {
+      ilu.Mult(b, x);
+      x *= alpha;
+   }
+
+   void SetOperator(const Operator &op) { }
+};
+
 // Class for constructing a multigrid preconditioner for the diffusion operator.
 // This example multigrid preconditioner class demonstrates the creation of the
 // diffusion bilinear forms and operators using partial assembly for all spaces
@@ -43,13 +80,18 @@ using namespace mfem;
 class DiffusionMultigrid : public GeometricMultigrid
 {
 private:
-   ConstantCoefficient one;
+   MatrixConstantCoefficient coeff;
+   bool use_ilu;
 
 public:
    // Constructs a diffusion multigrid for the given FiniteElementSpaceHierarchy
    // and the array of essential boundaries
-   DiffusionMultigrid(FiniteElementSpaceHierarchy& fespaces, Array<int>& ess_bdr)
-      : GeometricMultigrid(fespaces), one(1.0)
+   DiffusionMultigrid(FiniteElementSpaceHierarchy& fespaces, Array<int>& ess_bdr,
+                      double anisotropy, bool use_ilu_)
+      : GeometricMultigrid(fespaces),
+        coeff(AnisotropicCoefficient(fespaces.GetFinestFESpace().GetMesh()->Dimension(),
+                                     anisotropy)),
+        use_ilu(use_ilu_)
    {
       ConstructCoarseOperatorAndSolver(fespaces.GetFESpaceAtLevel(0), ess_bdr);
 
@@ -63,8 +105,7 @@ private:
    void ConstructBilinearForm(FiniteElementSpace& fespace, Array<int>& ess_bdr)
    {
       BilinearForm* form = new BilinearForm(&fespace);
-      form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
-      form->AddDomainIntegrator(new DiffusionIntegrator(one));
+      form->AddDomainIntegrator(new DiffusionIntegrator(coeff));
       form->Assemble();
       bfs.Append(form);
 
@@ -78,18 +119,13 @@ private:
       ConstructBilinearForm(coarse_fespace, ess_bdr);
 
       OperatorPtr opr;
-      opr.SetType(Operator::ANY_TYPE);
+      opr.SetType(Operator::MFEM_SPARSEMAT);
       bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), opr);
       opr.SetOperatorOwner(false);
 
-      CGSolver* pcg = new CGSolver();
-      pcg->SetPrintLevel(-1);
-      pcg->SetMaxIter(200);
-      pcg->SetRelTol(sqrt(1e-4));
-      pcg->SetAbsTol(0.0);
-      pcg->SetOperator(*opr.Ptr());
+      UMFPackSolver *coarse_solver = new UMFPackSolver(*opr.As<SparseMatrix>());
 
-      AddLevel(opr.Ptr(), pcg, true, true);
+      AddLevel(opr.Ptr(), coarse_solver, false, true);
    }
 
    void ConstructOperatorAndSmoother(FiniteElementSpace& fespace,
@@ -98,16 +134,25 @@ private:
       ConstructBilinearForm(fespace, ess_bdr);
 
       OperatorPtr opr;
-      opr.SetType(Operator::ANY_TYPE);
+      opr.SetType(Operator::MFEM_SPARSEMAT);
       bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), opr);
       opr.SetOperatorOwner(false);
 
-      Vector diag(fespace.GetTrueVSize());
-      bfs.Last()->AssembleDiagonal(diag);
+      Solver *smoother;
 
-      Solver* smoother = new OperatorChebyshevSmoother(*opr, diag,
-                                                       *essentialTrueDofs.Last(), 2);
-      AddLevel(opr.Ptr(), smoother, true, true);
+      if (use_ilu)
+      {
+         smoother = new SymmetricILUSmoother(*opr, 0.5);
+      }
+      else
+      {
+         Vector diag(fespace.GetTrueVSize());
+         bfs.Last()->AssembleDiagonal(diag);
+         smoother = new OperatorChebyshevSmoother(
+            *opr, diag, *essentialTrueDofs.Last(), 2);
+      }
+
+      AddLevel(opr.Ptr(), smoother, false, true);
    }
 };
 
@@ -120,6 +165,8 @@ int main(int argc, char *argv[])
    int order_refinements = 2;
    const char *device_config = "cpu";
    bool visualization = true;
+   double anisotropy = 1.0;
+   bool use_ilu = false;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -128,8 +175,11 @@ int main(int argc, char *argv[])
                   "Number of geometric refinements done prior to order refinements.");
    args.AddOption(&order_refinements, "-or", "--order-refinements",
                   "Number of order refinements. Finest level in the hierarchy has order 2^{or}.");
+   args.AddOption(&anisotropy, "-a", "--anisotropy", "Anisotropy coefficient.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
+   args.AddOption(&use_ilu, "-i", "--use-ilu", "-no-i", "--no-ilu",
+                  "Use ILU smoothing?");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -210,7 +260,7 @@ int main(int argc, char *argv[])
    Array<int> ess_bdr(mesh->bdr_attributes.Max());
    ess_bdr = 1;
 
-   DiffusionMultigrid M(fespaces, ess_bdr);
+   DiffusionMultigrid M(fespaces, ess_bdr, anisotropy, use_ilu);
    M.SetCycleType(Multigrid::CycleType::VCYCLE, 1, 1);
 
    OperatorPtr A;

fem/datacollection.cpp

@@ -11,6 +11,7 @@
 
 #include "fem.hpp"
 #include "../mesh/nurbs.hpp"
+#include "../mesh/vtk.hpp"
 #include "../general/binaryio.hpp"
 #include "../general/text.hpp"
 #include "picojson.h"
@@ -787,53 +788,44 @@ void ParaViewDataCollection::Load(int )
 
 std::string ParaViewDataCollection::GenerateCollectionPath()
 {
-   std::string out = "";
-   out = prefix_path + DataCollection::GetCollectionName();
-   return out;
+   return prefix_path + DataCollection::GetCollectionName();
 }
 
 std::string ParaViewDataCollection::GeneratePVTUPath()
 {
-   std::string out = "Cycle" + to_padded_string(cycle,pad_digits_cycle);
-   return out;
+   return "Cycle" + to_padded_string(cycle,pad_digits_cycle);
 }
 
 std::string ParaViewDataCollection::GenerateVTUPath()
 {
-   std::string out = GeneratePVTUPath();
-   return out;
+   return GeneratePVTUPath();
 }
 
 std::string ParaViewDataCollection::GeneratePVDFileName()
 {
-   std::string out = GetCollectionName()+".pvd";
-   return out;
+   return GetCollectionName() + ".pvd";
 }
 
-std::string ParaViewDataCollection::GeneratePVTUFileName()
+std::string ParaViewDataCollection::GeneratePVTUFileName(
+   const std::string &prefix)
 {
-   std::string out = "data.pvtu";
-   return out;
+   return prefix + ".pvtu";
 }
 
-std::string ParaViewDataCollection::GenerateVTUFileName()
+std::string ParaViewDataCollection::GenerateVTUFileName(
+   const std::string &prefix, int rank)
 {
-   std::string out = "proc" + to_padded_string(myid,pad_digits_rank)+".vtu";
-   return out;
-}
-std::string ParaViewDataCollection::GenerateVTUFileName(int crank)
-{
-   std::string out = "proc" + to_padded_string(crank,pad_digits_rank)+".vtu";
-   return out;
+   return prefix + to_padded_string(rank, pad_digits_rank) + ".vtu";
 }
 
 void ParaViewDataCollection::Save()
 {
    // add a new collection to the PDV file
 
+   std::string col_path = GenerateCollectionPath();
    // check if the directories are created
    {
-      std::string path = GenerateCollectionPath()+"/"+GenerateVTUPath();
+      std::string path = col_path + "/" + GenerateVTUPath();
       int err = create_directory(path, mesh, myid);
       if (err)
       {
@@ -850,8 +842,7 @@ void ParaViewDataCollection::Save()
 
    if (myid == 0 && !pvd_stream.is_open())
    {
-      std::string dpath=GenerateCollectionPath();
-      std::string pvdname=dpath+"/"+GeneratePVDFileName();
+      std::string pvdname = col_path + "/" + GeneratePVDFileName();
 
       bool write_header = true;
       std::ifstream pvd_in;
@@ -915,80 +906,87 @@ void ParaViewDataCollection::Save()
       }
    }
 
-   // define the vtu file
+   std::string vtu_prefix = col_path + "/" + GenerateVTUPath() + "/";
+
+   // Save the local part of the mesh and grid functions fields to the local
+   // VTU file
    {
-      std::string fname = GenerateCollectionPath()+"/"+GenerateVTUPath()+"/"
-                          +GenerateVTUFileName();
-      std::fstream out(fname, std::ios::out);
+      std::ofstream out(vtu_prefix + GenerateVTUFileName("proc", myid));
       out.precision(precision);
-      SaveDataVTU(out,levels_of_detail);
-      out.close();
+      SaveDataVTU(out, levels_of_detail);
    }
 
-   // define the pvtu file only on process 0
-   if (myid==0)
+   // Save the local part of the quadrature function fields
+   for (const auto &qfield : q_field_map)
    {
-      std::string fname = GenerateCollectionPath()+"/"+GeneratePVTUPath()+"/"
-                          +GeneratePVTUFileName();
-      std::fstream out(fname, std::ios::out);
+      const std::string &field_name = qfield.first;
+      std::ofstream out(vtu_prefix + GenerateVTUFileName(field_name, myid));
+      qfield.second->SaveVTU(out, pv_data_format, compression);
+   }
 
-      out << "<?xml version=\"1.0\"?>\n";
-      out << "<VTKFile type=\"PUnstructuredGrid\"";
-      out << " version =\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
-      out << "<PUnstructuredGrid GhostLevel=\"0\">\n";
-
-      out << "<PPoints>\n";
-      out << "\t<PDataArray type=\"" << GetDataTypeString() << "\" ";
-      out << " Name=\"Points\" NumberOfComponents=\"3\""
-          << " format=\"" << GetDataFormatString() << "\"/>\n";
-      out << "</PPoints>\n";
-
-      out << "<PCells>\n";
-      out << "\t<PDataArray type=\"Int32\" ";
-      out << " Name=\"connectivity\" NumberOfComponents=\"1\""
-          << " format=\"" << GetDataFormatString() << "\"/>\n";
-      out << "\t<PDataArray type=\"Int32\" ";
-      out << " Name=\"offsets\"      NumberOfComponents=\"1\""
-          << " format=\"" << GetDataFormatString() << "\"/>\n";
-      out << "\t<PDataArray type=\"UInt8\" ";
-      out << " Name=\"types\"        NumberOfComponents=\"1\""
-          << " format=\"" << GetDataFormatString() << "\"/>\n";
-      out << "</PCells>\n";
-
-      out << "<PPointData>\n";
-      for (FieldMapIterator it=field_map.begin(); it!=field_map.end(); ++it)
+   // MPI rank 0 also creates a "PVTU" file that points to all of the separately
+   // written VTU files.
+   // This file path is then appended to the PVD file.
+   if (myid == 0)
+   {
+      // Create the main PVTU file
       {
-         int vec_dim=it->second->VectorDim();
-         out << "<PDataArray type=\"" << GetDataTypeString()
-             << "\" Name=\"" << it->first
-             << "\" NumberOfComponents=\"" << vec_dim << "\" "
-             << "format=\"" << GetDataFormatString() << "\" />\n";
+         std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
+         WritePVTUHeader(pvtu_out);
+
+         // Grid function fields
+         pvtu_out << "<PPointData>\n";
+         for (auto &field_it : field_map)
+         {
+            int vec_dim = field_it.second->VectorDim();
+            pvtu_out << "<PDataArray type=\"" << GetDataTypeString()
+                     << "\" Name=\"" << field_it.first
+                     << "\" NumberOfComponents=\"" << vec_dim << "\" "
+                     << "format=\"" << GetDataFormatString() << "\" />\n";
+         }
+         pvtu_out << "</PPointData>\n";
+         // Element attributes
+         pvtu_out << "<PCellData>\n";
+         pvtu_out << "\t<PDataArray type=\"Int32\" Name=\"" << "attribute"
+                  << "\" NumberOfComponents=\"1\""
+                  << " format=\"" << GetDataFormatString() << "\"/>\n";
+         pvtu_out << "</PCellData>\n";
+
+         WritePVTUFooter(pvtu_out, "proc");
       }
-      out << "</PPointData>\n";
 
-      // CELL DATA
-      out << "<PCellData>\n";
-      out << "\t<PDataArray type=\"Int32\" Name=\"" << "attribute"
-          << "\" NumberOfComponents=\"1\""
-          << " format=\"" << GetDataFormatString() << "\"/>\n";
-      out << "</PCellData>\n";
+      // Add the latest PVTU to the PVD
+      pvd_stream << "<DataSet timestep=\"" << GetTime()
+                 << "\" group=\"\" part=\"" << 0 << "\" file=\""
+                 << GeneratePVTUPath() + "/" + GeneratePVTUFileName("data")
+                 << "\" name=\"mesh\"/>\n";
 
-      for (int ii=0; ii<num_procs; ii++)
+      // Create PVTU files for each quadrature field and add them to the PVD
+      // file
+      for (auto &q_field : q_field_map)
       {
-         // this one is generated without the path
-         std::string nfname=GenerateVTUFileName(ii);
-         out << "<Piece Source=\"" << nfname << "\"/>\n";
-      }
-      out << "</PUnstructuredGrid>\n";
-      out << "</VTKFile>\n";
-      out.close();
+         const std::string &q_field_name = q_field.first;
+         std::string q_fname = GeneratePVTUPath() + "/"
+                               + GeneratePVTUFileName(q_field_name);
 
-      fname = GeneratePVTUPath()+"/"+GeneratePVTUFileName();
-      // add the pvtu file to the pvd_stream
-      pvd_stream << "<DataSet timestep=\"" << GetTime();  // GetCycle();
-      pvd_stream << "\" group=\"\" part=\"" << 0 << "\" file=\"";
-      pvd_stream << fname << "\"/>\n";
+         std::ofstream pvtu_out(col_path + "/" + q_fname);
+         WritePVTUHeader(pvtu_out);
+         int vec_dim = q_field.second->GetVDim();
+         pvtu_out << "<PPointData>\n";
+         pvtu_out << "<PDataArray type=\"" << GetDataTypeString()
+                  << "\" Name=\"" << q_field_name
+                  << "\" NumberOfComponents=\"" << vec_dim << "\" "
+                  << "format=\"" << GetDataFormatString() << "\" />\n";
+         pvtu_out << "</PPointData>\n";
+         WritePVTUFooter(pvtu_out, q_field_name);
+
+         pvd_stream << "<DataSet timestep=\"" << GetTime()
+                    << "\" group=\"\" part=\"" << 0 << "\" file=\""
+                    << q_fname << "\" name=\"" << q_field_name << "\"/>\n";
+      }
       pvd_stream.flush();
+      // Move the insertion point before the closing collection tag, so that
+      // the PVD file is valid even when writing incrementally.
       std::fstream::pos_type pos = pvd_stream.tellp();
       pvd_stream << "</Collection>\n";
       pvd_stream << "</VTKFile>" << std::endl;
@@ -996,6 +994,44 @@ void ParaViewDataCollection::Save()
    }
 }
 
+void ParaViewDataCollection::WritePVTUHeader(std::ostream &out)
+{
+   out << "<?xml version=\"1.0\"?>\n";
+   out << "<VTKFile type=\"PUnstructuredGrid\"";
+   out << " version =\"0.1\" byte_order=\"" << VTKByteOrder() << "\">\n";
+   out << "<PUnstructuredGrid GhostLevel=\"0\">\n";
+
+   out << "<PPoints>\n";
+   out << "\t<PDataArray type=\"" << GetDataTypeString() << "\" ";
+   out << " Name=\"Points\" NumberOfComponents=\"3\""
+       << " format=\"" << GetDataFormatString() << "\"/>\n";
+   out << "</PPoints>\n";
+
+   out << "<PCells>\n";
+   out << "\t<PDataArray type=\"Int32\" ";
+   out << " Name=\"connectivity\" NumberOfComponents=\"1\""
+       << " format=\"" << GetDataFormatString() << "\"/>\n";
+   out << "\t<PDataArray type=\"Int32\" ";
+   out << " Name=\"offsets\"      NumberOfComponents=\"1\""
+       << " format=\"" << GetDataFormatString() << "\"/>\n";
+   out << "\t<PDataArray type=\"UInt8\" ";
+   out << " Name=\"types\"        NumberOfComponents=\"1\""
+       << " format=\"" << GetDataFormatString() << "\"/>\n";
+   out << "</PCells>\n";
+}
+
+void ParaViewDataCollection::WritePVTUFooter(std::ostream &out,
+                                             const std::string &vtu_prefix)
+{
+   for (int ii=0; ii<num_procs; ii++)
+   {
+      std::string vtu_filename = GenerateVTUFileName(vtu_prefix, ii);
+      out << "<Piece Source=\"" << vtu_filename << "\"/>\n";
+   }
+   out << "</PUnstructuredGrid>\n";
+   out << "</VTKFile>\n";
+}
+
 void ParaViewDataCollection::SaveDataVTU(std::ostream &out, int ref)
 {
    out << "<VTKFile type=\"UnstructuredGrid\"";
@@ -1015,16 +1051,6 @@ void ParaViewDataCollection::SaveDataVTU(std::ostream &out, int ref)
    {
       SaveGFieldVTU(out,ref,it);
    }
-   // iterate over all quadrature functions
-   // if the Quadrature functions are dumped as cell data
-   // the cycle should be moved before the grid functions
-   // and the PrintVTU CellData section should be open in the mesh dump
-   for (QFieldMapIterator it=q_field_map.begin(); it!=q_field_map.end(); ++it)
-   {
-      // save the quadrature functions
-      // this one is not implemented yet
-      SaveQFieldVTU(out,ref,it);
-   }
    out << "</PointData>\n";
    // close the mesh
    out << "</Piece>\n"; // close the piece open in the PrintVTU method
@@ -1032,27 +1058,21 @@ void ParaViewDataCollection::SaveDataVTU(std::ostream &out, int ref)
    out << "</VTKFile>" << std::endl;
 }
 
-void ParaViewDataCollection::SaveQFieldVTU(std::ostream &out, int ref,
-                                           const QFieldMapIterator& it )
-{
-   MFEM_WARNING("SaveQFieldVTU is not currently implemented - field name:"<<it->second);
-}
-
 void ParaViewDataCollection::SaveGFieldVTU(std::ostream &out, int ref_,
-                                           const FieldMapIterator& it)
+                                           const FieldMapIterator &it)
 {
    RefinedGeometry *RefG;
    Vector val;
    DenseMatrix vval, pmat;
    std::vector<char> buf;
    int vec_dim = it->second->VectorDim();
+   out << "<DataArray type=\"" << GetDataTypeString()
+       << "\" Name=\"" << it->first;
+   out << "\" NumberOfComponents=\"" << vec_dim << "\""
+       << " format=\"" << GetDataFormatString() << "\" >" << '\n';
    if (vec_dim == 1)
    {
       // scalar data
-      out << "<DataArray type=\"" << GetDataTypeString()
-          << "\" Name=\"" << it->first;
-      out << "\" NumberOfComponents=\"1\" format=\""
-          << GetDataFormatString() << "\" >\n";
       for (int i = 0; i < mesh->GetNE(); i++)
       {
          RefG = GlobGeometryRefiner.Refine(
@@ -1060,51 +1080,23 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &out, int ref_,
          it->second->GetValues(i, RefG->RefPts, val, pmat);
          for (int j = 0; j < val.Size(); j++)
          {
-            if (pv_data_format == VTKFormat::ASCII)
-            {
-               out << ZeroSubnormal(val(j)) << '\n';
-            }
-            else if (pv_data_format == VTKFormat::BINARY)
-            {
-               bin_io::AppendBytes(buf, val(j));
-            }
-            else
-            {
-               bin_io::AppendBytes<float>(buf, float(val(j)));
-            }
+            WriteBinaryOrASCII(out, buf, val(j), "\n", pv_data_format);
          }
       }
    }
    else
    {
       // vector data
-      out << "<DataArray type=\"" << GetDataTypeString()
-          << "\" Name=\"" << it->first;
-      out << "\" NumberOfComponents=\"" << vec_dim << "\""
-          << " format=\"" << GetDataFormatString() << "\" >" << '\n';
       for (int i = 0; i < mesh->GetNE(); i++)
       {
          RefG = GlobGeometryRefiner.Refine(
                    mesh->GetElementBaseGeometry(i), ref_, 1);
-
          it->second->GetVectorValues(i, RefG->RefPts, vval, pmat);
-
          for (int jj = 0; jj < vval.Width(); jj++)
          {
             for (int ii = 0; ii < vval.Height(); ii++)
             {
-               if (pv_data_format == VTKFormat::ASCII)
-               {
-                  out << ZeroSubnormal(vval(ii,jj)) << ' ';
-               }
-               else if (pv_data_format == VTKFormat::BINARY)
-               {
-                  bin_io::AppendBytes(buf, vval(ii,jj));
-               }
-               else
-               {
-                  bin_io::AppendBytes<float>(buf, float(vval(ii,jj)));
-               }
+               WriteBinaryOrASCII(out, buf, vval(ii,jj), " ", pv_data_format);
             }
             if (pv_data_format == VTKFormat::ASCII) { out << '\n'; }
          }

fem/datacollection.hpp

@@ -491,19 +491,20 @@ private:
    bool restart_mode;
 
 protected:
+   void WritePVTUHeader(std::ostream &out);
+   void WritePVTUFooter(std::ostream &out, const std::string &vtu_prefix);
    void SaveDataVTU(std::ostream &out, int ref);
    void SaveGFieldVTU(std::ostream& out, int ref_, const FieldMapIterator& it);
-   void SaveQFieldVTU(std::ostream &out, int ref, const QFieldMapIterator& it);
    const char *GetDataFormatString() const;
    const char *GetDataTypeString() const;
 
-   std::string  GenerateCollectionPath();
-   std::string  GenerateVTUFileName();
-   std::string  GenerateVTUFileName(int rank);
-   std::string  GenerateVTUPath();
-   std::string  GeneratePVDFileName();
-   std::string  GeneratePVTUFileName();
-   std::string  GeneratePVTUPath();
+   std::string GenerateCollectionPath();
+   std::string GenerateVTUFileName(const std::string &prefix, int rank);
+   std::string GenerateVTUPath();
+   std::string GeneratePVDFileName();
+   std::string GeneratePVTUFileName(const std::string &prefix);
+   std::string GeneratePVTUPath();
+
 
 public:
    /// Constructor. The collection name is used when saving the data.

fem/gridfunc.cpp

@@ -3948,6 +3948,129 @@ std::ostream &operator<<(std::ostream &out, const QuadratureFunction &qf)
    return out;
 }
 
+void QuadratureFunction::SaveVTU(std::ostream &out, VTKFormat format,
+                                 int compression_level) const
+{
+   out << R"(<VTKFile type="UnstructuredGrid" version="0.1")";
+   if (compression_level != 0)
+   {
+      out << R"( compressor="vtkZLibDataCompressor")";
+   }
+   out << " byte_order=\"" << VTKByteOrder() << "\">\n";
+   out << "<UnstructuredGrid>\n";
+
+   const char *fmt_str = (format == VTKFormat::ASCII) ? "ascii" : "binary";
+   const char *type_str = (format != VTKFormat::BINARY32) ? "Float64" : "Float32";
+   std::vector<char> buf;
+
+   int np = qspace->GetSize();
+   int ne = qspace->GetNE();
+   int sdim = qspace->GetMesh()->SpaceDimension();
+
+   // For quadrature functions, each point is a vertex cell, so number of cells
+   // is equal to number of points
+   out << "<Piece NumberOfPoints=\"" << np
+       << "\" NumberOfCells=\"" << np << "\">\n";
+
+   // print out the points
+   out << "<Points>\n";
+   out << "<DataArray type=\"" << type_str
+       << "\" NumberOfComponents=\"3\" format=\"" << fmt_str << "\">\n";
+
+   Vector pt(sdim);
+   for (int i = 0; i < ne; i++)
+   {
+      ElementTransformation &T = *qspace->GetMesh()->GetElementTransformation(i);
+      const IntegrationRule &ir = GetElementIntRule(i);
+      for (int j = 0; j < ir.Size(); j++)
+      {
+         T.Transform(ir[j], pt);
+         WriteBinaryOrASCII(out, buf, pt[0], " ", format);
+         if (sdim > 1) { WriteBinaryOrASCII(out, buf, pt[1], " ", format); }
+         else { WriteBinaryOrASCII(out, buf, 0.0, " ", format); }
+         if (sdim > 2) { WriteBinaryOrASCII(out, buf, pt[2], "", format); }
+         else { WriteBinaryOrASCII(out, buf, 0.0, "", format); }
+         if (format == VTKFormat::ASCII) { out << '\n'; }
+      }
+   }
+   if (format != VTKFormat::ASCII)
+   {
+      WriteBase64WithSizeAndClear(out, buf, compression_level);
+   }
+   out << "</DataArray>\n";
+   out << "</Points>\n";
+
+   // Write cells (each cell is just a vertex)
+   out << "<Cells>\n";
+   // Connectivity
+   out << R"(<DataArray type="Int32" Name="connectivity" format=")"
+       << fmt_str << "\">\n";
+
+   for (int i=0; i<np; ++i) { WriteBinaryOrASCII(out, buf, i, "\n", format); }
+   if (format != VTKFormat::ASCII)
+   {
+      WriteBase64WithSizeAndClear(out, buf, compression_level);
+   }
+   out << "</DataArray>\n";
+   // Offsets
+   out << R"(<DataArray type="Int32" Name="offsets" format=")"
+       << fmt_str << "\">\n";
+   for (int i=0; i<np; ++i) { WriteBinaryOrASCII(out, buf, i, "\n", format); }
+   if (format != VTKFormat::ASCII)
+   {
+      WriteBase64WithSizeAndClear(out, buf, compression_level);
+   }
+   out << "</DataArray>\n";
+   // Types
+   out << R"(<DataArray type="UInt8" Name="types" format=")"
+       << fmt_str << "\">\n";
+   for (int i = 0; i < np; i++)
+   {
+      uint8_t vtk_cell_type = VTKGeometry::POINT;
+      WriteBinaryOrASCII(out, buf, vtk_cell_type, "\n", format);
+   }
+   if (format != VTKFormat::ASCII)
+   {
+      WriteBase64WithSizeAndClear(out, buf, compression_level);
+   }
+   out << "</DataArray>\n";
+   out << "</Cells>\n";
+
+   out << "<PointData>\n";
+   out << "<DataArray type=\"" << type_str << "\" Name=\"u\" format=\""
+       << fmt_str << "\" NumberOfComponents=\"" << vdim << "\">\n";
+   for (int i = 0; i < ne; i++)
+   {
+      DenseMatrix vals;
+      GetElementValues(i, vals);
+      for (int j = 0; j < vals.Size(); ++j)
+      {
+         for (int vd = 0; vd < vdim; ++vd)
+         {
+            WriteBinaryOrASCII(out, buf, vals(vd, j), " ", format);
+         }
+         if (format == VTKFormat::ASCII) { out << '\n'; }
+      }
+   }
+   if (format != VTKFormat::ASCII)
+   {
+      WriteBase64WithSizeAndClear(out, buf, compression_level);
+   }
+   out << "</DataArray>\n";
+   out << "</PointData>\n";
+
+   out << "</Piece>\n";
+   out << "</UnstructuredGrid>\n";
+   out << "</VTKFile>" << std::endl;
+}
+
+void QuadratureFunction::SaveVTU(const std::string &filename, VTKFormat format,
+                                 int compression_level) const
+{
+   std::ofstream f(filename + ".vtu");
+   SaveVTU(f, format, compression_level);
+}
+
 
 double ZZErrorEstimator(BilinearFormIntegrator &blfi,
                         GridFunction &u,

fem/gridfunc.hpp

@@ -902,6 +902,22 @@ public:
 
    /// Write the QuadratureFunction to the stream @a out.
    void Save(std::ostream &out) const;
+
+   /// @brief Write the QuadratureFunction to @a out in VTU (ParaView) format.
+   ///
+   /// The data will be uncompressed if @a compression_level is zero, or if the
+   /// format is VTKFormat::ASCII. Otherwise, zlib compression will be used for
+   /// binary data.
+   void SaveVTU(std::ostream &out, VTKFormat format=VTKFormat::ASCII,
+                int compression_level=0) const;
+
+   /// @brief Save the QuadratureFunction to a VTU (ParaView) file.
+   ///
+   /// The extension ".vtu" will be appended to @a filename.
+   /// @sa SaveVTU(std::ostream &out, VTKFormat format=VTKFormat::ASCII,
+   ///             int compression_level=0)
+   void SaveVTU(const std::string &filename, VTKFormat format=VTKFormat::ASCII,
+                int compression_level=0) const;
 };
 
 /// Overload operator<< for std::ostream and QuadratureFunction.

fem/transfer.cpp

@@ -903,7 +903,8 @@ TransferOperator::TransferOperator(const FiniteElementSpace& lFESpace_,
                                    const FiniteElementSpace& hFESpace_)
    : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize())
 {
-   if (lFESpace_.FEColl() == hFESpace_.FEColl())
+   bool isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
+   if (lFESpace_.FEColl() == hFESpace_.FEColl() && !isvar_order)
    {
       OperatorPtr P(Operator::ANY_TYPE);
       hFESpace_.GetTransferOperator(lFESpace_, P);
@@ -912,8 +913,11 @@ TransferOperator::TransferOperator(const FiniteElementSpace& lFESpace_,
    }
    else if (lFESpace_.GetMesh()->GetNE() > 0
             && hFESpace_.GetMesh()->GetNE() > 0
+            && lFESpace_.GetVDim() == 1
+            && hFESpace_.GetVDim() == 1
             && dynamic_cast<const TensorBasisElement*>(lFESpace_.GetFE(0))
             && dynamic_cast<const TensorBasisElement*>(hFESpace_.GetFE(0))
+            && !isvar_order
             && (hFESpace_.FEColl()->GetContType() ==
                 mfem::FiniteElementCollection::CONTINUOUS ||
                 hFESpace_.FEColl()->GetContType() ==
@@ -945,6 +949,7 @@ PRefinementTransferOperator::PRefinementTransferOperator(
    : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
      hFESpace(hFESpace_)
 {
+   isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
 }
 
 PRefinementTransferOperator::~PRefinementTransferOperator() {}
@@ -969,7 +974,7 @@ void PRefinementTransferOperator::Mult(const Vector& x, Vector& y) const
       DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
 
       const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
-      if (geom != cached_geom)
+      if (geom != cached_geom || isvar_order)
       {
          h_fe = hFESpace.GetFE(i);
          l_fe = lFESpace.GetFE(i);
@@ -1026,7 +1031,7 @@ void PRefinementTransferOperator::MultTranspose(const Vector& x,
       DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
 
       const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
-      if (geom != cached_geom)
+      if (geom != cached_geom || isvar_order)
       {
          h_fe = hFESpace.GetFE(i);
          l_fe = lFESpace.GetFE(i);
@@ -1424,20 +1429,36 @@ void TensorProductPRefinementTransferOperator::MultTranspose(const Vector& x,
    elem_restrict_lex_l->MultTranspose(localL, y);
 }
 
-#ifdef MFEM_USE_MPI
-TrueTransferOperator::TrueTransferOperator(const
-                                           ParFiniteElementSpace& lFESpace_,
-                                           const ParFiniteElementSpace& hFESpace_)
+
+TrueTransferOperator::TrueTransferOperator(const FiniteElementSpace& lFESpace_,
+                                           const FiniteElementSpace& hFESpace_)
    : Operator(hFESpace_.GetTrueVSize(), lFESpace_.GetTrueVSize()),
      lFESpace(lFESpace_),
      hFESpace(hFESpace_)
 {
    localTransferOperator = new TransferOperator(lFESpace_, hFESpace_);
 
-   tmpL.SetSize(lFESpace_.GetVSize());
-   tmpH.SetSize(hFESpace_.GetVSize());
+   P = lFESpace.GetProlongationMatrix();
+   R = hFESpace.IsVariableOrder() ? hFESpace.GetHpRestrictionMatrix() :
+       hFESpace.GetRestrictionMatrix();
 
-   hFESpace.GetRestrictionMatrix()->BuildTranspose();
+   // P and R can be both null
+   // P can be null and R not null
+   // If P is not null it is assumed that R is not null as well
+   if (P) { MFEM_VERIFY(R, "Both P and R have to be not NULL") }
+
+   if (P)
+   {
+      tmpL.SetSize(lFESpace_.GetVSize());
+      tmpH.SetSize(hFESpace_.GetVSize());
+      R->BuildTranspose();
+   }
+   // P can be null and R not null
+   else if (R)
+   {
+      tmpH.SetSize(hFESpace_.GetVSize());
+      R->BuildTranspose();
+   }
 }
 
 TrueTransferOperator::~TrueTransferOperator()
@@ -1447,17 +1468,40 @@ TrueTransferOperator::~TrueTransferOperator()
 
 void TrueTransferOperator::Mult(const Vector& x, Vector& y) const
 {
-   lFESpace.GetProlongationMatrix()->Mult(x, tmpL);
-   localTransferOperator->Mult(tmpL, tmpH);
-   hFESpace.GetRestrictionMatrix()->Mult(tmpH, y);
+   if (P)
+   {
+      P->Mult(x, tmpL);
+      localTransferOperator->Mult(tmpL, tmpH);
+      R->Mult(tmpH, y);
+   }
+   else if (R)
+   {
+      localTransferOperator->Mult(x, tmpH);
+      R->Mult(tmpH, y);
+   }
+   else
+   {
+      localTransferOperator->Mult(x, y);
+   }
 }
 
 void TrueTransferOperator::MultTranspose(const Vector& x, Vector& y) const
 {
-   hFESpace.GetRestrictionMatrix()->MultTranspose(x, tmpH);
-   localTransferOperator->MultTranspose(tmpH, tmpL);
-   lFESpace.GetProlongationMatrix()->MultTranspose(tmpL, y);
+   if (P)
+   {
+      R->MultTranspose(x, tmpH);
+      localTransferOperator->MultTranspose(tmpH, tmpL);
+      P->MultTranspose(tmpL, y);
+   }
+   else if (R)
+   {
+      R->MultTranspose(x, tmpH);
+      localTransferOperator->MultTranspose(tmpH, y);
+   }
+   else
+   {
+      localTransferOperator->MultTranspose(x, y);
+   }
 }
-#endif
 
 } // namespace mfem

fem/transfer.hpp

@@ -387,6 +387,7 @@ class PRefinementTransferOperator : public Operator
 private:
    const FiniteElementSpace& lFESpace;
    const FiniteElementSpace& hFESpace;
+   bool isvar_order;
 
 public:
    /// @brief Constructs a transfer operator from \p lFESpace to \p hFESpace
@@ -452,14 +453,15 @@ public:
    virtual void MultTranspose(const Vector& x, Vector& y) const override;
 };
 
-#ifdef MFEM_USE_MPI
 /// @brief Matrix-free transfer operator between finite element spaces working
 /// on true degrees of freedom
 class TrueTransferOperator : public Operator
 {
 private:
-   const ParFiniteElementSpace& lFESpace;
-   const ParFiniteElementSpace& hFESpace;
+   const FiniteElementSpace& lFESpace;
+   const FiniteElementSpace& hFESpace;
+   const Operator * P = nullptr;
+   const SparseMatrix * R = nullptr;
    TransferOperator* localTransferOperator;
    mutable Vector tmpL;
    mutable Vector tmpH;
@@ -467,8 +469,8 @@ private:
 public:
    /// @brief Constructs a transfer operator working on true degrees of freedom
    /// from \p lFESpace to \p hFESpace
-   TrueTransferOperator(const ParFiniteElementSpace& lFESpace_,
-                        const ParFiniteElementSpace& hFESpace_);
+   TrueTransferOperator(const FiniteElementSpace& lFESpace_,
+                        const FiniteElementSpace& hFESpace_);
 
    /// Destructor
    ~TrueTransferOperator();
@@ -484,7 +486,6 @@ public:
        the true dof vector \p y corresponding to the coarse space. */
    virtual void MultTranspose(const Vector& x, Vector& y) const override;
 };
-#endif
 
 } // namespace mfem
 

general/binaryio.hpp

@@ -58,16 +58,19 @@ void AppendBytes(std::vector<char> &vec, const T &val)
    vec.insert(vec.end(), ptr, ptr + sizeof(T));
 }
 
-/// Given a buffer @a buf of length @a nbytes, encode the data in base-64
-/// format, and write the encoded data to the output stream @a out.
+/// @brief Given a buffer @a bytes of length @a nbytes, encode the data in
+/// base-64 format, and write the encoded data to the output stream @a out.
 void WriteBase64(std::ostream &out, const void *bytes, size_t nbytes);
 
-/// Decode @a len base-64 encoded characters in the buffer @a src, and store the
-/// resulting decoded data in @a buf. @a buf will be resized as needed.
+/// @brief Decode @a len base-64 encoded characters in the buffer @a src, and
+/// store the resulting decoded data in @a buf. @a buf will be resized as
+/// needed.
 void DecodeBase64(const char *src, size_t len, std::vector<char> &buf);
 
-/// Return the number of characters needed to encode @a nbytes in base-64. This
-/// is equal to 4*nbytes/3, rounded up to the nearest multiple of 4.
+/// @brief Return the number of characters needed to encode @a nbytes in
+/// base-64.
+///
+/// This is equal to 4*nbytes/3, rounded up to the nearest multiple of 4.
 size_t NumBase64Chars(size_t nbytes);
 
 } // namespace mfem::bin_io

mesh/mesh.cpp

@@ -10152,61 +10152,6 @@ void Mesh::PrintBdrVTU(std::string fname,
    PrintVTU(fname, format, high_order_output, compression_level, true);
 }
 
-template <typename T>
-void WriteBinaryOrASCII(std::ostream &out, std::vector<char> &buf, const T &val,
-                        const char *suffix, VTKFormat format)
-{
-   if (format == VTKFormat::ASCII) { out << val << suffix; }
-   else { bin_io::AppendBytes(buf, val); }
-}
-
-// Ensure ASCII output of uint8_t to stream is integer rather than character
-template <>
-void WriteBinaryOrASCII<uint8_t>(std::ostream &out, std::vector<char> &buf,
-                                 const uint8_t &val, const char *suffix,
-                                 VTKFormat format)
-{
-   if (format == VTKFormat::ASCII) { out << static_cast<int>(val) << suffix; }
-   else { bin_io::AppendBytes(buf, val); }
-}
-
-template <>
-void WriteBinaryOrASCII<double>(std::ostream &out, std::vector<char> &buf,
-                                const double &val, const char *suffix,
-                                VTKFormat format)
-{
-   if (format == VTKFormat::BINARY32)
-   {
-      bin_io::AppendBytes<float>(buf, float(val));
-   }
-   else if (format == VTKFormat::BINARY)
-   {
-      bin_io::AppendBytes(buf, val);
-   }
-   else
-   {
-      out << val << suffix;
-   }
-}
-
-template <>
-void WriteBinaryOrASCII<float>(std::ostream &out, std::vector<char> &buf,
-                               const float &val, const char *suffix,
-                               VTKFormat format)
-{
-   if (format == VTKFormat::BINARY) { bin_io::AppendBytes<double>(buf, val); }
-   else if (format == VTKFormat::BINARY32) { bin_io::AppendBytes(buf, val); }
-   else { out << val << suffix; }
-}
-
-void WriteBase64WithSizeAndClear(std::ostream &out, std::vector<char> &buf,
-                                 int compression_level)
-{
-   WriteVTKEncodedCompressed(out, buf.data(), buf.size(), compression_level);
-   out << '\n';
-   buf.clear();
-}
-
 void Mesh::PrintVTU(std::ostream &out, int ref, VTKFormat format,
                     bool high_order_output, int compression_level,
                     bool bdr_elements)

mesh/vtk.cpp

@@ -600,4 +600,51 @@ const char *VTKByteOrder()
 
 }
 
+// Ensure ASCII output of uint8_t to stream is integer rather than character
+template <>
+void WriteBinaryOrASCII<uint8_t>(std::ostream &out, std::vector<char> &buf,
+                                 const uint8_t &val, const char *suffix,
+                                 VTKFormat format)
+{
+   if (format == VTKFormat::ASCII) { out << static_cast<int>(val) << suffix; }
+   else { bin_io::AppendBytes(buf, val); }
+}
+
+template <>
+void WriteBinaryOrASCII<double>(std::ostream &out, std::vector<char> &buf,
+                                const double &val, const char *suffix,
+                                VTKFormat format)
+{
+   if (format == VTKFormat::BINARY32)
+   {
+      bin_io::AppendBytes<float>(buf, float(val));
+   }
+   else if (format == VTKFormat::BINARY)
+   {
+      bin_io::AppendBytes(buf, val);
+   }
+   else
+   {
+      out << ZeroSubnormal(val) << suffix;
+   }
+}
+
+template <>
+void WriteBinaryOrASCII<float>(std::ostream &out, std::vector<char> &buf,
+                               const float &val, const char *suffix,
+                               VTKFormat format)
+{
+   if (format == VTKFormat::BINARY) { bin_io::AppendBytes<double>(buf, val); }
+   else if (format == VTKFormat::BINARY32) { bin_io::AppendBytes(buf, val); }
+   else { out << ZeroSubnormal(val) << suffix; }
+}
+
+void WriteBase64WithSizeAndClear(std::ostream &out, std::vector<char> &buf,
+                                 int compression_level)
+{
+   WriteVTKEncodedCompressed(out, buf.data(), buf.size(), compression_level);
+   out << '\n';
+   buf.clear();
+}
+
 } // namespace mfem

mesh/vtk.hpp

@@ -13,23 +13,34 @@
 #define MFEM_VTK
 
 #include "../fem/geom.hpp"
+#include "../general/binaryio.hpp"
 
 namespace mfem
 {
 
 // Helpers for reading and writing VTK format
 
-// VTK element types defined at: https://git.io/JvZLm
+/// @brief Helper class for converting between MFEM and VTK geometry types.
+///
+/// Note: The VTK element types defined are at: https://git.io/JvZLm
 struct VTKGeometry
 {
+   /// @name VTK geometry types
+   ///@{
    static const int POINT = 1;
+
+   /// @name Low-order (linear, straight-sided) VTK geometric types
+   ///@{
    static const int SEGMENT = 3;
    static const int TRIANGLE = 5;
    static const int SQUARE = 9;
    static const int TETRAHEDRON = 10;
    static const int CUBE = 12;
    static const int PRISM = 13;
+   ///@}
 
+   /// @name Legacy quadratic VTK geometric types
+   ///@{
    static const int QUADRATIC_SEGMENT = 21;
    static const int QUADRATIC_TRIANGLE = 22;
    static const int BIQUADRATIC_SQUARE = 28;
@@ -37,49 +48,138 @@ struct VTKGeometry
    static const int TRIQUADRATIC_CUBE = 29;
    static const int QUADRATIC_PRISM = 26;
    static const int BIQUADRATIC_QUADRATIC_PRISM = 32;
+   ///@}
 
+   /// @name Arbitrary-order VTK geometric types
+   ///@{
    static const int LAGRANGE_SEGMENT = 68;
    static const int LAGRANGE_TRIANGLE = 69;
    static const int LAGRANGE_SQUARE = 70;
    static const int LAGRANGE_TETRAHEDRON = 71;
    static const int LAGRANGE_CUBE = 72;
    static const int LAGRANGE_PRISM = 73;
+   ///@}
+   ///@}
 
+   /// Permutation from MFEM's prism ordering to VTK's prism ordering.
    static const int PrismMap[6];
+
+   /// @brief Permutation from MFEM's vertex ordering to VTK's vertex ordering.
+   /// @note If the MFEM and VTK orderings are the same, the vertex permutation
+   /// will be NULL.
    static const int *VertexPermutation[Geometry::NUM_GEOMETRIES];
 
+   /// Map from MFEM's Geometry::Type to linear VTK geometries.
    static const int Map[Geometry::NUM_GEOMETRIES];
+   /// Map from MFEM's Geometry::Type to legacy quadratic VTK geometries/
    static const int QuadraticMap[Geometry::NUM_GEOMETRIES];
+   /// Map from MFEM's Geometry::Type to arbitrary-order Lagrange VTK geometries
    static const int HighOrderMap[Geometry::NUM_GEOMETRIES];
 
+   /// Given a VTK geometry type, return the corresponding MFEM Geometry::Type.
    static Geometry::Type GetMFEMGeometry(int vtk_geom);
+   /// @brief Does the given VTK geometry type describe an arbitrary-order
+   /// Lagrange element?
    static bool IsLagrange(int vtk_geom);
+   /// @brief Does the given VTK geometry type describe a legacy quadratic
+   /// element?
    static bool IsQuadratic(int vtk_geom);
+   /// @brief For the given VTK geometry type and number of points, return the
+   /// order of the element.
    static int GetOrder(int vtk_geom, int npoints);
 };
 
+/// Data array format for VTK and VTU files.
 enum class VTKFormat
 {
+   /// Data arrays will be written in ASCII format.
    ASCII,
+   /// Data arrays will be written in binary format. Floating point numbers will
+   /// be be output with 64 bits of precision.
    BINARY,
+   /// Data arrays will be written in binary format. Floating point numbers will
+   /// be be output with 32 bits of precision.
    BINARY32
 };
 
-/// Create the VTK element connectivity array for a given element geometry and
-/// refinement level. Converts node numbers from MFEM to VTK ordering.
+/// @brief Create the VTK element connectivity array for a given element
+/// geometry and refinement level.
+///
+/// The output array @a con will be such that, for the @a ith VTK node index,
+/// con[i] will contain the index of the corresponding node in MFEM ordering.
 void CreateVTKElementConnectivity(Array<int> &con, Geometry::Type geom,
                                   int ref);
 
-/// Outputs encoded binary data in the format needed by VTK. The binary data
-/// will be base 64 encoded, and compressed if @a compression_level is not
-/// zero. The proper header will be prepended to the data.
+/// @brief Outputs encoded binary data in the base 64 format needed by VTK.
+///
+/// The binary data will be base 64 encoded, and compressed if @a
+/// compression_level is not zero. The proper header will be prepended to the
+/// data.
 void WriteVTKEncodedCompressed(std::ostream &out, const void *bytes,
                                uint32_t nbytes, int compression_level);
 
+/// @brief Return the VTK node index of the barycentric point @a b in a
+/// triangle with refinement level @a ref.
+///
+/// The barycentric index @a b has three components, satisfying b[0] + b[1] +
+/// b[2] == ref.
 int BarycentricToVTKTriangle(int *b, int ref);
 
+/// Determine the byte order and return either "BigEndian" or "LittleEndian"
 const char *VTKByteOrder();
 
+/// @brief Write either ASCII data to the stream or binary data to the buffer
+/// depending on the given format.
+///
+/// If @a format is VTK::ASCII, write the canonical ASCII representation of @a
+/// val to the output stream. Subnormal floating point numbers are rounded to
+/// zero. Otherwise, append its raw binary data to the byte buffer @a buf.
+///
+/// Note that there are specializations for @a uint8_t (to write as a numeric
+/// value rather than a character), and for @a float and @a double values to use
+/// the precision specified by @a format.
+template <typename T>
+void WriteBinaryOrASCII(std::ostream &out, std::vector<char> &buf, const T &val,
+                        const char *suffix, VTKFormat format)
+{
+   if (format == VTKFormat::ASCII) { out << val << suffix; }
+   else { bin_io::AppendBytes(buf, val); }
+}
+
+/// @brief Specialization of @ref WriteBinaryOrASCII for @a uint8_t to ensure
+/// ASCII output is numeric (rather than interpreting @a val as a character.)
+template <>
+void WriteBinaryOrASCII<uint8_t>(std::ostream &out, std::vector<char> &buf,
+                                 const uint8_t &val, const char *suffix,
+                                 VTKFormat format);
+
+/// @brief Specialization of @ref WriteBinaryOrASCII for @a double.
+///
+/// If @a format is equal to VTKFormat::BINARY32, @a val is converted to a @a
+/// float and written as 32 bits. Subnormals are rounded to zero in ASCII
+/// output.
+template <>
+void WriteBinaryOrASCII<double>(std::ostream &out, std::vector<char> &buf,
+                                const double &val, const char *suffix,
+                                VTKFormat format);
+
+/// @brief Specialization of @ref WriteBinaryOrASCII<T> for @a float.
+///
+/// If @a format is equal to VTKFormat::BINARY, @a val is converted to a @a
+/// double and written as 64 bits. Subnormals are rounded to zero in ASCII
+/// output.
+template <>
+void WriteBinaryOrASCII<float>(std::ostream &out, std::vector<char> &buf,
+                               const float &val, const char *suffix,
+                               VTKFormat format);
+
+/// @brief Encode in base 64 (and potentially compress) the given data, write it
+/// to the output stream (with a header) and clear the buffer.
+///
+/// @sa WriteVTKEncodedCompressed.
+void WriteBase64WithSizeAndClear(std::ostream &out, std::vector<char> &buf,
+                                 int compression_level);
+
 } // namespace mfem
 
 #endif

miniapps/common/pfem_extras.cpp

@@ -310,7 +310,7 @@ void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
 }
 
 void VisualizeField(socketstream &sock, const char *vishost, int visport,
-                    ParGridFunction &gf, const char *title,
+                    const ParGridFunction &gf, const char *title,
                     int x, int y, int w, int h, const char *keys, bool vec)
 {
    ParMesh &pmesh = *gf.ParFESpace()->GetParMesh();

miniapps/common/pfem_extras.hpp

@@ -197,7 +197,7 @@ void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
 /// specified host and port. Set the visualization window title, and optionally,
 /// its geometry.
 void VisualizeField(socketstream &sock, const char *vishost, int visport,
-                    ParGridFunction &gf, const char *title,
+                    const ParGridFunction &gf, const char *title,
                     int x = 0, int y = 0, int w = 400, int h = 400,
                     const char *keys = NULL, bool vec = false);
 

miniapps/shifted/CMakeLists.txt

@@ -13,12 +13,14 @@ if (MFEM_USE_MPI)
   list(APPEND DIST_COMMON_SOURCES
         dist_solver.cpp
         sbm_solver.cpp
-        marking.cpp)
+        marking.cpp
+        extrapolator.cpp)
   list(APPEND DIST_COMMON_HEADERS
         dist_solver.hpp
         sbm_solver.hpp
         sbm_aux.hpp
-        marking.hpp)
+        marking.hpp
+        extrapolator.hpp)
 
   convert_filenames_to_full_paths(DIST_COMMON_SOURCES)
   convert_filenames_to_full_paths(DIST_COMMON_HEADERS)
@@ -37,6 +39,11 @@ if (MFEM_USE_MPI)
     ${DIST_COMMON_FILES}
     LIBRARIES mfem mfem-common)
 
+  add_mfem_miniapp(extrapolate
+    MAIN extrapolate.cpp
+    ${DIST_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+
   if (MFEM_ENABLE_TESTING)
     add_test(NAME shifted_distance_np${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}

miniapps/shifted/extrapolate.cpp

@@ -0,0 +1,218 @@
+// Copyright (c) 2010-2021, 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.
+//
+//            ------------------------------------------------
+//            Extrapolation Miniapp: PDE-based extrapolation
+//            ------------------------------------------------
+//
+// This miniapp extrapolates a finite element function from a set of elements
+// (known values) to the rest of the domain. The set of elements that contains
+// the known values is specified by the positive values of a level set
+// Coefficient. The known values are not modified. The miniapp supports two
+// PDE-based approaches [1, 2], both of which rely on solving a sequence of
+// advection problems in the direction of the unknown parts of the domain.
+// The extrapolation can be constant (1st order), linear (2nd order), or
+// quadratic (3rd order). These formal orders hold for a limited band around
+// the zero level set, see the given references for more info.
+//
+// [1] Aslam, "A Partial Differential Equation Approach to Multidimensional
+// Extrapolation", JCP 193(1), 2004.
+// [2] Bochkov, Gibou, "PDE-Based Multidimensional Extrapolation of Scalar
+// Fields over Interfaces with Kinks and High Curvatures", SISC 42(4), 2020.
+//
+// Compile with: make extrapolate
+//
+// Sample runs:
+//     mpirun -np 4 extrapolate -m "../../data/inline-segment.mesh" -rs 6 -ed 2
+//     mpirun -np 4 extrapolate -rs 5 -p 0 -ed 2
+//     mpirun -np 4 extrapolate -rs 5 -p 1 -ed 2
+//     mpirun -np 4 extrapolate -rs 5 -p 1 -et 1 -ed 1 -dg 1
+//     mpirun -np 4 extrapolate -m "../../data/inline-hex.mesh" -ed 1 -rs 1
+//     mpirun -np 4 extrapolate -m "../../data/inline-hex.mesh" -p 1 -ed 1 -rs 1
+
+#include "extrapolator.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int problem = 0;
+
+double domainLS(const Vector &coord)
+{
+   // Map from [0,1] to [-1,1].
+   const int dim = coord.Size();
+   const double x = coord(0)*2.0 - 1.0,
+                y = (dim > 1) ? coord(1)*2.0 - 1.0 : 0.0,
+                z = (dim > 2) ? coord(2)*2.0 - 1.0 : 0.0;
+
+   switch (problem)
+   {
+      case 0:
+      {
+         // Sphere.
+         return 0.75 - sqrt(x*x + y*y + z*z + 1e-12);
+      }
+      case 1:
+      {
+         // Star.
+         MFEM_VERIFY(dim > 1, "Problem 1 is not applicable to 1D.");
+
+         return 0.60 - sqrt(x*x + y*y + z*z + 1e-12) +
+                0.25 * (y*y*y*y*y + 5.0*x*x*x*x*y - 10.0*x*x*y*y*y) /
+                pow(x*x + y*y + z*z + 1e-12, 2.5) *
+                std::cos(0.5*M_PI * z / 0.6);
+      }
+      default: MFEM_ABORT("Bad option for --problem!"); return 0.0;
+   }
+}
+
+double solution0(const Vector &coord)
+{
+   // Map from [0,1] to [-1,1].
+   const int dim = coord.Size();
+   const double x = coord(0)*2.0 - 1.0 + 0.25,
+                y = (dim > 1) ? coord(1)*2.0 - 1.0 : 0.0,
+                z = (dim > 2) ? coord(2)*2.0 - 1.0 : 0.0;
+
+   return std::cos(M_PI * x) * std::cos(M_PI * y) * std::cos(M_PI * z);
+}
+
+void PrintNorm(int myid, Vector &v, std::string text)
+{
+   double norm = v.Norml1();
+   MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   if (myid == 0)
+   {
+      std::cout << std::setprecision(12) << std::fixed
+                << text << norm << std::endl;
+   }
+}
+
+void PrintIntegral(int myid, ParGridFunction &g, std::string text)
+{
+   ConstantCoefficient zero(0.0);
+   double norm = g.ComputeL1Error(zero);
+   if (myid == 0)
+   {
+      std::cout << std::setprecision(12) << std::fixed
+                << text << norm << std::endl;
+   }
+}
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI.
+   MPI_Session mpi;
+   int myid = mpi.WorldRank();
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int rs_levels = 2;
+   Extrapolator::XtrapType ex_type   = Extrapolator::ASLAM;
+   AdvectionOper::AdvectionMode dg_mode = AdvectionOper::HO;
+   int ex_degree = 1;
+   int order = 2;
+   double distance = 0.35;
+   bool vis_on = true;
+   int vis_steps_cnt = 50;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption((int*)&ex_type, "-et", "--extrap-type",
+                  "Extrapolation type: Aslam (0) or Bochkov (1).");
+   args.AddOption((int*)&dg_mode, "-dg", "--dg-mode",
+                  "DG advection mode: 0 - Standard High-Order,\n\t"
+                  "                   1 - Low-Order Upwind Diffusion.");
+   args.AddOption(&ex_degree, "-ed", "--extrap-degree",
+                  "Extrapolation degree: 0/1/2 for constant/linear/quadratic.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&distance, "-d", "--distance",
+                  "Extrapolation distance.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "0 - 2D circle,\n\t"
+                  "1 - 2D star");
+   args.AddOption(&vis_on, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&vis_steps_cnt, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0) { args.PrintUsage(cout); }
+      return 1;
+   }
+   if (myid == 0) { args.PrintOptions(cout); }
+
+   // Refine the mesh and distribute.
+   Mesh mesh(mesh_file, 1, 1);
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   const int dim = pmesh.Dimension();
+
+   // Input function.
+   L2_FECollection fec_L2(order, dim);
+   ParFiniteElementSpace pfes_L2(&pmesh, &fec_L2);
+   ParGridFunction u(&pfes_L2);
+   FunctionCoefficient u0_coeff(solution0);
+   u.ProjectCoefficient(u0_coeff);
+
+   // Extrapolate.
+   Extrapolator xtrap;
+   xtrap.xtrap_type     = ex_type;
+   xtrap.advection_mode = dg_mode;
+   xtrap.xtrap_degree   = ex_degree;
+   xtrap.visualization  = vis_on;
+   xtrap.vis_steps      = vis_steps_cnt;
+   FunctionCoefficient ls_coeff(domainLS);
+   ParGridFunction ux(&pfes_L2);
+   xtrap.Extrapolate(ls_coeff, u, distance, ux);
+
+   PrintNorm(myid, ux, "Solution l1 norm: ");
+   PrintIntegral(myid, ux, "Solution L1 norm: ");
+
+   GridFunctionCoefficient u_exact_coeff(&u);
+   double err_L1 = ux.ComputeL1Error(u_exact_coeff),
+          err_L2 = ux.ComputeL2Error(u_exact_coeff);
+   if (myid == 0)
+   {
+      std::cout << "Global L1 error: " << err_L1 << std::endl
+                << "Global L2 error: " << err_L2 << std::endl;
+   }
+   double loc_error_L1, loc_error_L2, loc_error_LI;
+   xtrap.ComputeLocalErrors(ls_coeff, u, ux,
+                            loc_error_L1, loc_error_L2, loc_error_LI);
+   if (myid == 0)
+   {
+      std::cout << "Local  L1 error: " << loc_error_L1 << std::endl
+                << "Local  L2 error: " << loc_error_L2 << std::endl
+                << "Local  Li error: " << loc_error_LI << std::endl;
+   }
+
+   // ParaView output.
+   ParGridFunction ls_gf(&pfes_L2);
+   ls_gf.ProjectCoefficient(ls_coeff);
+   ParaViewDataCollection dacol("ParaViewExtrapolate", &pmesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.RegisterField("Level Set Function", &ls_gf);
+   dacol.RegisterField("Extrapolated Solution", &ux);
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+   return 0;
+}

miniapps/shifted/extrapolator.cpp

@@ -0,0 +1,565 @@
+// Copyright (c) 2010-2021, 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 "extrapolator.hpp"
+#include "../common/mfem-common.hpp"
+#include "marking.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+const char vishost[] = "localhost";
+const int  visport   = 19916;
+int wsize            = 350; // glvis window size
+
+AdvectionOper::AdvectionOper(Array<bool> &zones, ParBilinearForm &Mbf,
+                             ParBilinearForm &Kbf, const Vector &rhs)
+   : TimeDependentOperator(Mbf.Size()),
+     active_zones(zones),
+     M(Mbf), K(Kbf), K_mat(NULL), b(rhs),
+     lo_solver(NULL), lumpedM(NULL)
+{
+   K_mat = K.ParallelAssemble(&K.SpMat());
+
+   ParBilinearForm M_Lump(M.ParFESpace());
+   lumpedM = new Vector;
+   M_Lump.AddDomainIntegrator(new LumpedIntegrator(new MassIntegrator));
+   M_Lump.Assemble();
+   M_Lump.Finalize();
+   M_Lump.SpMat().GetDiag(*lumpedM);
+   lo_solver = new DiscreteUpwindLOSolver(*M.ParFESpace(),
+                                          K.SpMat(), *lumpedM);
+}
+
+AdvectionOper::~AdvectionOper()
+{
+   delete lo_solver;
+   delete lumpedM;
+   delete K_mat;
+}
+
+void AdvectionOper::Mult(const Vector &x, Vector &dx) const
+{
+   ParFiniteElementSpace &pfes = *M.ParFESpace();
+   const int NE = pfes.GetNE();
+   const int nd = pfes.GetFE(0)->GetDof();
+   Array<int> dofs(nd);
+
+   if (adv_mode == LO)
+   {
+      lo_solver->CalcLOSolution(x, b, dx);
+      for (int k = 0; k < NE; k++)
+      {
+         pfes.GetElementDofs(k, dofs);
+         if (active_zones[k] == false)
+         {
+            dx.SetSubVector(dofs, 0.0);
+            continue;
+         }
+      }
+      return;
+   }
+
+   MFEM_VERIFY(adv_mode == HO, "Wrong input for avection mode (-dg).");
+
+   Vector rhs(x.Size());
+   K_mat->Mult(x, rhs);
+   rhs += b;
+
+   DenseMatrix M_loc(nd);
+   DenseMatrixInverse M_loc_inv(&M_loc);
+   Vector rhs_loc(nd), dx_loc(nd);
+   for (int k = 0; k < NE; k++)
+   {
+      pfes.GetElementDofs(k, dofs);
+
+      if (active_zones[k] == false)
+      {
+         dx.SetSubVector(dofs, 0.0);
+         continue;
+      }
+
+      rhs.GetSubVector(dofs, rhs_loc);
+      M.SpMat().GetSubMatrix(dofs, dofs, M_loc);
+      M_loc_inv.Factor();
+      M_loc_inv.Mult(rhs_loc, dx_loc);
+      dx.SetSubVector(dofs, dx_loc);
+   }
+}
+
+void AdvectionOper::ComputeElementsMinMax(const ParGridFunction &gf,
+                                          Vector &el_min, Vector &el_max) const
+{
+   ParFiniteElementSpace &pfes = *gf.ParFESpace();
+   const int NE = pfes.GetNE(), ndof = pfes.GetFE(0)->GetDof();
+   for (int k = 0; k < NE; k++)
+   {
+      el_min(k) = numeric_limits<double>::infinity();
+      el_max(k) = -numeric_limits<double>::infinity();
+
+      for (int i = 0; i < ndof; i++)
+      {
+         el_min(k) = min(el_min(k), gf(k*ndof + i));
+         el_max(k) = max(el_max(k), gf(k*ndof + i));
+      }
+   }
+}
+
+void AdvectionOper::ComputeBounds(const ParFiniteElementSpace &pfes,
+                                  const Vector &el_min, const Vector &el_max,
+                                  Vector &dof_min, Vector &dof_max) const
+{
+   ParMesh *pmesh = pfes.GetParMesh();
+   L2_FECollection fec_bounds(0, pmesh->Dimension());
+   ParFiniteElementSpace pfes_bounds(pmesh, &fec_bounds);
+   ParGridFunction el_min_gf(&pfes_bounds), el_max_gf(&pfes_bounds);
+   const int NE = pmesh->GetNE(), ndofs = dof_min.Size() / NE;
+
+   el_min_gf = el_min;
+   el_max_gf = el_max;
+
+   el_min_gf.ExchangeFaceNbrData(); el_max_gf.ExchangeFaceNbrData();
+   const Vector &min_nbr = el_min_gf.FaceNbrData();
+   const Vector &max_nbr = el_max_gf.FaceNbrData();
+   const Table &el_to_el = pmesh->ElementToElementTable();
+   Array<int> face_nbr_el;
+   for (int k = 0; k < NE; k++)
+   {
+      double k_min = el_min_gf(k), k_max = el_max_gf(k);
+
+      el_to_el.GetRow(k, face_nbr_el);
+      for (int n = 0; n < face_nbr_el.Size(); n++)
+      {
+         if (face_nbr_el[n] < NE)
+         {
+            // Local neighbor.
+            k_min = std::min(k_min, el_min_gf(face_nbr_el[n]));
+            k_max = std::max(k_max, el_max_gf(face_nbr_el[n]));
+         }
+         else
+         {
+            // MPI face neighbor.
+            k_min = std::min(k_min, min_nbr(face_nbr_el[n] - NE));
+            k_max = std::max(k_max, max_nbr(face_nbr_el[n] - NE));
+         }
+      }
+
+      for (int j = 0; j < ndofs; j++)
+      {
+         dof_min(k*ndofs + j) = k_min;
+         dof_max(k*ndofs + j) = k_max;
+      }
+   }
+}
+
+void Extrapolator::Extrapolate(Coefficient &level_set,
+                               const ParGridFunction &input,
+                               const double time_period,
+                               ParGridFunction &xtrap)
+{
+   ParMesh &pmesh = *input.ParFESpace()->GetParMesh();
+   const int order = input.ParFESpace()->GetOrder(0),
+             dim   = pmesh.Dimension(), NE = pmesh.GetNE();
+
+   // Get a ParGridFunction and mark elements.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace pfes_H1(&pmesh, &fec);
+   ParGridFunction ls_gf(&pfes_H1);
+   ls_gf.ProjectCoefficient(level_set);
+   if (visualization)
+   {
+      socketstream sock1, sock2;
+      common::VisualizeField(sock1, vishost, visport, ls_gf,
+                             "Domain level set", 0, 0, wsize, wsize,
+                             "rRjlmm********A");
+      common::VisualizeField(sock2, vishost, visport, input,
+                             "Input u", 0, wsize+60, wsize, wsize,
+                             "rRjlmm********A");
+      MPI_Barrier(pmesh.GetComm());
+   }
+   // Mark elements.
+   Array<int> elem_marker;
+   ShiftedFaceMarker marker(pmesh, pfes_H1, false);
+   ls_gf.ExchangeFaceNbrData();
+   marker.MarkElements(ls_gf, elem_marker);
+
+   // The active zones are where we extrapolate (where the PDE is solved).
+   Array<bool> active_zones(NE);
+   for (int k = 0; k < NE; k++)
+   {
+      // Extrapolation is done in zones that are CUT or OUTSIDE.
+      active_zones[k] =
+         (elem_marker[k] == ShiftedFaceMarker::INSIDE) ? false : true;
+   }
+
+   // Setup a VectorCoefficient for n = - grad_ls / |grad_ls|.
+   // The sign makes it point out of the known region.
+   // The coefficient must be continuous to have well-defined transport.
+   LevelSetNormalGradCoeff ls_n_coeff_L2(ls_gf);
+   ParFiniteElementSpace pfes_H1_vec(&pmesh, &fec, dim);
+   ParGridFunction lsn_gf(&pfes_H1_vec);
+   ls_gf.ExchangeFaceNbrData();
+   lsn_gf.ProjectDiscCoefficient(ls_n_coeff_L2, GridFunction::ARITHMETIC);
+   VectorGridFunctionCoefficient ls_n_coeff(&lsn_gf);
+
+   // Initial solution.
+   // Trim to the known values (only elements inside the known region).
+   Array<int> dofs;
+   L2_FECollection fec_L2(order, dim);
+   ParFiniteElementSpace pfes_L2(&pmesh, &fec_L2);
+   ParGridFunction u(&pfes_L2), vis_marking(&pfes_L2);
+   u.ProjectGridFunction(input);
+   for (int k = 0; k < NE; k++)
+   {
+      pfes_L2.GetElementDofs(k, dofs);
+      if (elem_marker[k] != ShiftedFaceMarker::INSIDE)
+      { u.SetSubVector(dofs, 0.0); }
+      vis_marking.SetSubVector(dofs, elem_marker[k]);
+   }
+   if (visualization)
+   {
+      socketstream sock1, sock2;
+      common::VisualizeField(sock1, vishost, visport, u,
+                             "Fixed (known) u values", wsize, 0,
+                             wsize, wsize, "rRjlmm********A");
+      common::VisualizeField(sock2, vishost, visport, vis_marking,
+                             "Element markings", 0, 2*wsize+60,
+                             wsize, wsize, "rRjlmm********A");
+   }
+
+   // Normal derivative function.
+   ParGridFunction n_grad_u(&pfes_L2);
+   NormalGradCoeff n_grad_u_coeff(u, ls_n_coeff);
+   n_grad_u.ProjectCoefficient(n_grad_u_coeff);
+   if (visualization && xtrap_degree >= 1)
+   {
+      socketstream sock;
+      common::VisualizeField(sock, vishost, visport, n_grad_u,
+                             "n.grad(u)", 2*wsize, 0, wsize, wsize,
+                             "rRjlmm********A");
+   }
+
+   // 2nd normal derivative function.
+   ParGridFunction n_grad_n_grad_u(&pfes_L2);
+   NormalGradCoeff n_grad_n_grad_u_coeff(n_grad_u, ls_n_coeff);
+   n_grad_n_grad_u.ProjectCoefficient(n_grad_n_grad_u_coeff);
+   if (visualization && xtrap_degree == 2)
+   {
+      socketstream sock;
+      common::VisualizeField(sock, vishost, visport, n_grad_n_grad_u,
+                             "n.grad(n.grad(u))", 3*wsize, 0, wsize, wsize,
+                             "rRjmm********A");
+   }
+
+   ParBilinearForm lhs_bf(&pfes_L2), rhs_bf(&pfes_L2);
+   lhs_bf.AddDomainIntegrator(new MassIntegrator);
+   const double alpha = -1.0;
+   rhs_bf.AddDomainIntegrator(new ConvectionIntegrator(ls_n_coeff, alpha));
+   auto trace_i = new NonconservativeDGTraceIntegrator(ls_n_coeff, alpha);
+   rhs_bf.AddInteriorFaceIntegrator(trace_i);
+   rhs_bf.KeepNbrBlock(true);
+
+   ls_gf.ExchangeFaceNbrData();
+   lhs_bf.Assemble();
+   lhs_bf.Finalize();
+   rhs_bf.Assemble(0);
+   rhs_bf.Finalize(0);
+
+   // Compute a CFL time step.
+   double h_min = std::numeric_limits<double>::infinity();
+   for (int k = 0; k < NE; k++)
+   {
+      h_min = std::min(h_min, pmesh.GetElementSize(k));
+   }
+   MPI_Allreduce(MPI_IN_PLACE, &h_min, 1, MPI_DOUBLE, MPI_MIN, pmesh.GetComm());
+   // The propagation speed is 1.
+   double dt = 0.25 * h_min / order / 1.0;
+   double half_dt = 0.5 * dt;
+   if (advection_mode == AdvectionOper::LO)
+   {
+      dt = half_dt;
+   }
+
+   // Time loops.
+   Vector rhs(pfes_L2.GetVSize());
+   AdvectionOper adv_oper(active_zones, lhs_bf, rhs_bf, rhs);
+   adv_oper.adv_mode = advection_mode;
+   RK2Solver ode_solver(1.0);
+   ode_solver.Init(adv_oper);
+
+   if (xtrap_degree == 0)
+   {
+      // Constant extrapolation of u (always LO).
+      rhs = 0.0;
+      adv_oper.adv_mode = AdvectionOper::LO;
+      TimeLoop(u, ode_solver, time_period, half_dt,
+               wsize, "Extrap const u -- LO");
+      xtrap.ProjectGridFunction(u);
+      return;
+   }
+
+   std::string mode_text = "HO";
+   if (advection_mode == AdvectionOper::LO)  { mode_text = "LO"; }
+
+   MFEM_VERIFY(xtrap_degree == 1 || xtrap_degree == 2, "Wrong order input.");
+   if (xtrap_type == ASLAM)
+   {
+      if (xtrap_degree == 1)
+      {
+         // Constant extrapolation of [n.grad_u] (always LO).
+         rhs = 0.0;
+         adv_oper.adv_mode = AdvectionOper::LO;
+         TimeLoop(n_grad_u, ode_solver, time_period, half_dt,
+                  2*wsize, "Extrap const n.grad(u) -- Aslam -- LO");
+
+         adv_oper.adv_mode = advection_mode;
+
+         // Linear extrapolation of u.
+         lhs_bf.Mult(n_grad_u, rhs);
+         TimeLoop(u, ode_solver, time_period, dt,
+                  wsize, "Extrap linear u -- Aslam -- " + mode_text);
+      }
+
+      if (xtrap_degree == 2)
+      {
+         // Constant extrapolation of [n.grad(n.grad(u))] (always LO).
+         rhs = 0.0;
+         adv_oper.adv_mode = AdvectionOper::LO;
+         TimeLoop(n_grad_n_grad_u, ode_solver, time_period, half_dt,
+                  3*wsize, "Extrap const n.grad(n.grad(u)) -- Aslam -- LO");
+
+         adv_oper.adv_mode = advection_mode;
+
+         // Linear extrapolation of [n.grad_u].
+         lhs_bf.Mult(n_grad_n_grad_u, rhs);
+         TimeLoop(n_grad_u, ode_solver, time_period, dt,
+                  2*wsize, "Extrap linear n.grad(u) -- Aslam -- " + mode_text);
+
+         // Quadratic extrapolation of u.
+         lhs_bf.Mult(n_grad_u, rhs);
+         TimeLoop(u, ode_solver, time_period, dt,
+                  wsize, "Extrap quadratic u -- Aslam -- " + mode_text);
+      }
+   }
+   else if (xtrap_type == BOCHKOV)
+   {
+      if (xtrap_degree == 1)
+      {
+         // Constant extrapolation of all grad(u) components (always LO).
+         rhs = 0.0;
+         adv_oper.adv_mode = AdvectionOper::LO;
+         ParGridFunction grad_u_0(&pfes_L2), grad_u_1(&pfes_L2);
+         GradComponentCoeff grad_u_0_coeff(u, 0), grad_u_1_coeff(u, 1);
+         grad_u_0.ProjectCoefficient(grad_u_0_coeff);
+         grad_u_1.ProjectCoefficient(grad_u_1_coeff);
+         TimeLoop(grad_u_0, ode_solver, time_period, half_dt,
+                  2*wsize, "Extrap const du_dx -- Bochkov -- LO");
+         TimeLoop(grad_u_1, ode_solver, time_period, half_dt,
+                  3*wsize, "Extrap const du_dy -- Bochkov -- LO");
+
+         adv_oper.adv_mode = advection_mode;
+
+         // Linear extrapolation of u.
+         ParLinearForm rhs_lf(&pfes_L2);
+         NormalGradComponentCoeff grad_u_n(grad_u_0, grad_u_1, ls_n_coeff);
+         rhs_lf.AddDomainIntegrator(new DomainLFIntegrator(grad_u_n));
+         rhs_lf.Assemble();
+         rhs = rhs_lf;
+         TimeLoop(u, ode_solver, time_period, dt,
+                  wsize, "Extrap linear u -- Bochkov -- " + mode_text);
+      }
+
+      if (xtrap_degree == 2)
+      {
+         MFEM_ABORT("Quadratic Bochkov method is not implemented.");
+      }
+   }
+   else { MFEM_ABORT("Wrong input for extrapolation type (-et)."); }
+
+   xtrap.ProjectGridFunction(u);
+}
+
+// Errors in cut elements.
+void Extrapolator::ComputeLocalErrors(Coefficient &level_set,
+                                      const ParGridFunction &exact,
+                                      const ParGridFunction &xtrap,
+                                      double &err_L1, double &err_L2,
+                                      double &err_LI)
+{
+   ParMesh &pmesh = *exact.ParFESpace()->GetParMesh();
+   const int order = exact.ParFESpace()->GetOrder(0),
+             dim   = pmesh.Dimension(), NE = pmesh.GetNE();
+
+   // Get a ParGridFunction and mark elements.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace pfes_H1(&pmesh, &fec);
+   ParGridFunction ls_gf(&pfes_H1);
+   ls_gf.ProjectCoefficient(level_set);
+   // Mark elements.
+   Array<int> elem_marker;
+   ShiftedFaceMarker marker(pmesh, pfes_H1, false);
+   ls_gf.ExchangeFaceNbrData();
+   marker.MarkElements(ls_gf, elem_marker);
+
+   Vector errors_L1(NE), errors_L2(NE), errors_LI(NE);
+   GridFunctionCoefficient exact_coeff(&exact);
+
+   xtrap.ComputeElementL1Errors(exact_coeff, errors_L1);
+   xtrap.ComputeElementL2Errors(exact_coeff, errors_L2);
+   xtrap.ComputeElementMaxErrors(exact_coeff, errors_LI);
+   err_L1 = 0.0, err_L2 = 0.0, err_LI = 0.0;
+   double cut_volume = 0.0;
+   for (int k = 0; k < NE; k++)
+   {
+      if (elem_marker[k] == ShiftedFaceMarker::CUT)
+      {
+         err_L1 += errors_L1(k);
+         err_L2 += errors_L2(k);
+         err_LI = std::max(err_LI, errors_LI(k));
+         cut_volume += pmesh.GetElementVolume(k);
+      }
+   }
+   MPI_Comm comm = pmesh.GetComm();
+   MPI_Allreduce(MPI_IN_PLACE, &err_L1, 1, MPI_DOUBLE, MPI_SUM, comm);
+   MPI_Allreduce(MPI_IN_PLACE, &err_L2, 1, MPI_DOUBLE, MPI_SUM, comm);
+   MPI_Allreduce(MPI_IN_PLACE, &err_LI, 1, MPI_DOUBLE, MPI_MAX, comm);
+   MPI_Allreduce(MPI_IN_PLACE, &cut_volume, 1, MPI_DOUBLE, MPI_SUM, comm);
+   err_L1 /= cut_volume;
+   err_L2 /= cut_volume;
+}
+
+void Extrapolator::TimeLoop(ParGridFunction &sltn, ODESolver &ode_solver,
+                            double t_final, double dt,
+                            int vis_x_pos, std::string vis_name)
+{
+   socketstream sock;
+
+   const int myid  = sltn.ParFESpace()->GetMyRank();
+   bool done = false;
+   double t = 0.0;
+   for (int ti = 0; !done;)
+   {
+      double dt_real = min(dt, t_final - t);
+      ode_solver.Step(sltn, t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+      if (done || ti % vis_steps == 0)
+      {
+         if (myid == 0)
+         {
+            cout << vis_name+" / time step: " << ti << ", time: " << t << endl;
+         }
+         if (visualization)
+         {
+            common::VisualizeField(sock, vishost, visport, sltn,
+                                   vis_name.c_str(), vis_x_pos, wsize+60,
+                                   wsize, wsize, "rRjlmm********A");
+            MPI_Barrier(sltn.ParFESpace()->GetComm());
+         }
+      }
+   }
+}
+
+
+
+DiscreteUpwindLOSolver::DiscreteUpwindLOSolver(ParFiniteElementSpace &space,
+                                               const SparseMatrix &adv,
+                                               const Vector &Mlump)
+   : pfes(space), K(adv), D(adv), K_smap(), M_lumped(Mlump)
+{
+   // Assuming it is finalized.
+   const int *I = K.GetI(), *J = K.GetJ(), n = K.Size();
+   K_smap.SetSize(I[n]);
+   for (int row = 0, j = 0; row < n; row++)
+   {
+      for (int end = I[row+1]; j < end; j++)
+      {
+         int col = J[j];
+         // Find the offset, _j, of the (col,row) entry and store it in smap[j].
+         for (int _j = I[col], _end = I[col+1]; true; _j++)
+         {
+            MFEM_VERIFY(_j != _end, "Can't find the symmetric entry!");
+
+            if (J[_j] == row) { K_smap[j] = _j; break; }
+         }
+      }
+   }
+
+   ComputeDiscreteUpwindMatrix();
+}
+
+void DiscreteUpwindLOSolver::CalcLOSolution(const Vector &u, const Vector &rhs,
+                                            Vector &du) const
+{
+   ParGridFunction u_gf(&pfes);
+   u_gf = u;
+   ApplyDiscreteUpwindMatrix(u_gf, du);
+
+   const int s = du.Size();
+   for (int i = 0; i < s; i++)
+   {
+      du(i) = (du(i) + rhs(i)) / M_lumped(i);
+   }
+}
+
+void DiscreteUpwindLOSolver::ComputeDiscreteUpwindMatrix() const
+{
+   const int *I = K.HostReadI(), *J = K.HostReadJ(), n = K.Size();
+
+   const double *K_data = K.HostReadData();
+
+   double *D_data = D.HostReadWriteData();
+   D.HostReadWriteI(); D.HostReadWriteJ();
+
+   for (int i = 0, k = 0; i < n; i++)
+   {
+      double rowsum = 0.;
+      for (int end = I[i+1]; k < end; k++)
+      {
+         int j = J[k];
+         double kij = K_data[k];
+         double kji = K_data[K_smap[k]];
+         double dij = fmax(fmax(0.0,-kij),-kji);
+         D_data[k] = kij + dij;
+         D_data[K_smap[k]] = kji + dij;
+         if (i != j) { rowsum += dij; }
+      }
+      D(i,i) = K(i,i) - rowsum;
+   }
+}
+
+void DiscreteUpwindLOSolver::ApplyDiscreteUpwindMatrix(ParGridFunction &u,
+                                                       Vector &du) const
+{
+   const int s = u.Size();
+   const int *I = D.HostReadI(), *J = D.HostReadJ();
+   const double *D_data = D.HostReadData();
+
+   u.ExchangeFaceNbrData();
+   const Vector &u_np = u.FaceNbrData();
+
+   for (int i = 0; i < s; i++)
+   {
+      du(i) = 0.0;
+      for (int k = I[i]; k < I[i + 1]; k++)
+      {
+         int j = J[k];
+         double u_j  = (j < s) ? u(j) : u_np[j - s];
+         double d_ij = D_data[k];
+         du(i) += d_ij * u_j;
+      }
+   }
+}
+
+}

miniapps/shifted/extrapolator.hpp

@@ -0,0 +1,191 @@
+// Copyright (c) 2010-2021, 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_EXTRAPOLATOR_HPP
+#define MFEM_EXTRAPOLATOR_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+class DiscreteUpwindLOSolver;
+class FluxBasedFCT;
+
+class AdvectionOper : public TimeDependentOperator
+{
+private:
+   Array<bool> &active_zones;
+   ParBilinearForm &M, &K;
+   HypreParMatrix *K_mat;
+   const Vector &b;
+
+   DiscreteUpwindLOSolver *lo_solver;
+   Vector *lumpedM;
+
+   void ComputeElementsMinMax(const ParGridFunction &gf,
+                              Vector &el_min, Vector &el_max) const;
+   void ComputeBounds(const ParFiniteElementSpace &pfes,
+                      const Vector &el_min, const Vector &el_max,
+                      Vector &dof_min, Vector &dof_max) const;
+   void ZeroOutInactiveZones(Vector &dx);
+
+public:
+   // HO is standard FE advection solve; LO is upwind diffusion.
+   enum AdvectionMode {HO, LO} adv_mode = AdvectionOper::HO;
+
+   AdvectionOper(Array<bool> &zones, ParBilinearForm &Mbf,
+                 ParBilinearForm &Kbf, const Vector &rhs);
+
+   ~AdvectionOper();
+
+   virtual void Mult(const Vector &x, Vector &dx) const;
+};
+
+// Extrapolates through DG advection based on:
+// [1] Aslam, "A Partial Differential Equation Approach to Multidimensional
+// Extrapolation", JCP 193(1), 2004.
+// [2] Bochkov, Gibou, "PDE-Based Multidimensional Extrapolation of Scalar
+// Fields over Interfaces with Kinks and High Curvatures", SISC 42(4), 2020.
+class Extrapolator
+{
+public:
+   enum XtrapType {ASLAM, BOCHKOV} xtrap_type = ASLAM;
+   AdvectionOper::AdvectionMode advection_mode = AdvectionOper::HO;
+   int xtrap_degree   = 1;
+   bool visualization = false;
+   int vis_steps      = 5;
+
+   Extrapolator() { }
+
+   // The known values taken from elements where level_set > 0, and extrapolated
+   // to all other elements. The known values are not changed.
+   void Extrapolate(Coefficient &level_set, const ParGridFunction &input,
+                    const double time_period, ParGridFunction &xtrap);
+
+   // Errors in cut elements, given an exact solution.
+   void ComputeLocalErrors(Coefficient &level_set, const ParGridFunction &exact,
+                           const ParGridFunction &xtrap,
+                           double &err_L1, double &err_L2, double &err_LI);
+
+private:
+   void TimeLoop(ParGridFunction &sltn, ODESolver &ode_solver, double t_final,
+                 double dt, int vis_x_pos, std::string vis_name);
+};
+
+class LevelSetNormalGradCoeff : public VectorCoefficient
+{
+private:
+   const ParGridFunction &ls_gf;
+
+public:
+   LevelSetNormalGradCoeff(const ParGridFunction &ls) :
+      VectorCoefficient(ls.ParFESpace()->GetMesh()->Dimension()), ls_gf(ls) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      Vector grad_ls(vdim), n(vdim);
+      ls_gf.GetGradient(T, grad_ls);
+      const double norm_grad = grad_ls.Norml2();
+      V = grad_ls;
+      if (norm_grad > 0.0) { V /= norm_grad; }
+
+      // Since positive level set values correspond to the known region, we
+      // transport into the opposite direction of the gradient.
+      V *= -1;
+   }
+};
+
+class GradComponentCoeff : public Coefficient
+{
+private:
+   const ParGridFunction &u_gf;
+   int comp;
+
+public:
+   GradComponentCoeff(const ParGridFunction &u, int c) : u_gf(u), comp(c) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      Vector grad_u(T.GetDimension());
+      u_gf.GetGradient(T, grad_u);
+      return grad_u(comp);
+   }
+};
+
+class NormalGradCoeff : public Coefficient
+{
+private:
+   const ParGridFunction &u_gf;
+   VectorCoefficient &n_coeff;
+
+public:
+   NormalGradCoeff(const ParGridFunction &u, VectorCoefficient &n)
+      : u_gf(u), n_coeff(n) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      const int dim = T.GetDimension();
+      Vector n(dim), grad_u(dim);
+      n_coeff.Eval(n, T, ip);
+      u_gf.GetGradient(T, grad_u);
+      return n * grad_u;
+   }
+};
+
+class NormalGradComponentCoeff : public Coefficient
+{
+private:
+   const ParGridFunction &du_dx, &du_dy;
+   VectorCoefficient &n_coeff;
+
+public:
+   NormalGradComponentCoeff(const ParGridFunction &dx,
+                            const ParGridFunction &dy, VectorCoefficient &n)
+      : du_dx(dx), du_dy(dy), n_coeff(n) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      const int dim = T.GetDimension();
+      Vector n(dim), grad_u(dim);
+      n_coeff.Eval(n, T, ip);
+      grad_u(0) = du_dx.GetValue(T, ip);
+      grad_u(1) = du_dy.GetValue(T, ip);
+      return n * grad_u;
+   }
+};
+
+class DiscreteUpwindLOSolver
+{
+public:
+   DiscreteUpwindLOSolver(ParFiniteElementSpace &space, const SparseMatrix &adv,
+                          const Vector &Mlump);
+
+   void CalcLOSolution(const Vector &u, const Vector &rhs, Vector &du) const;
+
+   Array<int> &GetKmap() { return K_smap; }
+
+protected:
+   ParFiniteElementSpace &pfes;
+   const SparseMatrix &K;
+   mutable SparseMatrix D;
+
+   Array<int> K_smap;
+   const Vector &M_lumped;
+
+   void ComputeDiscreteUpwindMatrix() const;
+   void ApplyDiscreteUpwindMatrix(ParGridFunction &u, Vector &du) const;
+};
+
+} // namespace mfem
+
+#endif

miniapps/shifted/makefile

@@ -29,8 +29,10 @@ DIFFUSION_SRC = diffusion.cpp dist_solver.cpp sbm_solver.cpp marking.cpp
 DIFFUSION_OBJ = $(DIFFUSION_SRC:.cpp=.o)
 DISTANCE_SRC = distance.cpp dist_solver.cpp
 DISTANCE_OBJ = $(DISTANCE_SRC:.cpp=.o)
+EXTRAPOLATE_SRC = extrapolate.cpp extrapolator.cpp marking.cpp
+EXTRAPOLATE_OBJ = $(EXTRAPOLATE_SRC:.cpp=.o)
 
-PAR_MINIAPPS = distance diffusion
+PAR_MINIAPPS = distance diffusion extrapolate
 
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS =
@@ -65,6 +67,9 @@ distance: $(DISTANCE_OBJ)
 diffusion: $(DIFFUSION_OBJ)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(DIFFUSION_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
 
+extrapolate: $(EXTRAPOLATE_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(EXTRAPOLATE_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
 # Rule for building lib-common
 lib-common:
 	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
@@ -89,9 +94,9 @@ $(MFEM_LIB_FILE):
 clean: clean-build clean-exec
 
 clean-build:
-	rm -f *.o *~ distance diffusion
+	rm -f *.o *~ distance diffusion extrapolate
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
 	@rm -f diffusion.mesh diffusion.gf
-	@rm -rf ParaViewDistance ParaViewDiffusion
+	@rm -rf ParaViewDistance ParaViewDiffusion ParaViewExtrapolate

miniapps/shifted/marking.hpp

@@ -51,7 +51,9 @@ public:
         include_cut_cell(include_cut_cell_), initial_marking_done(false),
         level_set_index(0) { }
 
-   /// Mark all the elements in the mesh using the @a SBElementType
+   /// Mark all the elements in the mesh using the @a SBElementType.
+   /// A point is considered inside when the level set function is positive.
+   /// Assumes the ExchangeFaceNbrData() has been called for pmesh, ls_func.
    void MarkElements(const ParGridFunction &ls_func, Array<int> &elem_marker);
 
    /// List dofs associated with the surrogate boundary.

tests/unit/fem/test_transfer.cpp

@@ -14,16 +14,54 @@
 
 using namespace mfem;
 
-int dimension;
-double coeff(const Vector& x)
+int RandomPRefinement(FiniteElementSpace & fes)
 {
+   Mesh * mesh = fes.GetMesh();
+   int maxorder = 0;
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      const int order = fes.GetElementOrder(i);
+      maxorder = std::max(maxorder,order);
+      if ((double) rand() / RAND_MAX < 0.5)
+      {
+         fes.SetElementOrder(i,order+1);
+         maxorder = std::max(maxorder,order+1);
+      }
+   }
+   fes.Update(false);
+   return maxorder;
+}
+
+
+int dimension;
+int coeff_order;
+double coeff(const Vector& X)
+{
+   double x = X[0];
+   double y = X[1];
+   double z = 0.;
    if (dimension == 2)
    {
-      return 1.1 * x[0] + 2.0 * x[1];
+      if (coeff_order == 1)
+      {
+         return 1.1 * x + 2.0 * y;
+      }
+      else
+      {
+         return (1.-x)*x*(1.-y)*y;
+      }
    }
    else
    {
-      return 1.1 * x[0] + 2.0 * x[1] + 3.0 * x[2];
+      z = X[2];
+      if (coeff_order == 1)
+      {
+         return 1.1 * x + 2.0 * y + 3.0 * z;
+      }
+      else
+      {
+         return (1.-x)*x*(1.-y)*y*(1.-z)*z;
+      }
    }
 }
 
@@ -37,7 +75,6 @@ void vectorcoeff(const Vector& x, Vector& y)
    }
 }
 
-
 TEST_CASE("transfer")
 {
    for (int vectorspace = 0; vectorspace <= 3; ++vectorspace)
@@ -82,26 +119,26 @@ TEST_CASE("transfer")
                         }
                         mesh = Mesh::MakeCartesian3D(ne, ne, ne, type, 1.0, 1.0, 1.0);
                      }
-                     FiniteElementCollection* c_h1_fec = nullptr;
-                     FiniteElementCollection* f_h1_fec = nullptr;
+                     FiniteElementCollection* c_fec = nullptr;
+                     FiniteElementCollection* f_fec = nullptr;
 
                      if (vectorspace < 2)
                      {
-                        c_h1_fec = new H1_FECollection(order, dimension);
-                        f_h1_fec = (geometric == 1) ? c_h1_fec : new
-                                   H1_FECollection(fineOrder, dimension);
+                        c_fec = new H1_FECollection(order, dimension);
+                        f_fec = (geometric == 1) ? c_fec : new
+                                H1_FECollection(fineOrder, dimension);
                      }
                      else if (vectorspace == 2)
                      {
-                        c_h1_fec = new ND_FECollection(order+1, dimension);
-                        f_h1_fec = (geometric == 1) ? c_h1_fec : new
-                                   ND_FECollection(fineOrder, dimension);
+                        c_fec = new ND_FECollection(order+1, dimension);
+                        f_fec = (geometric == 1) ? c_fec : new
+                                ND_FECollection(fineOrder, dimension);
                      }
                      else
                      {
-                        c_h1_fec = new RT_FECollection(order, dimension);
-                        f_h1_fec = (geometric == 1) ? c_h1_fec : new
-                                   RT_FECollection(fineOrder, dimension);
+                        c_fec = new RT_FECollection(order, dimension);
+                        f_fec = (geometric == 1) ? c_fec : new
+                                RT_FECollection(fineOrder, dimension);
                      }
 
                      Mesh fineMesh(mesh);
@@ -117,34 +154,34 @@ TEST_CASE("transfer")
                         spaceDimension = dimension;
                      }
 
-                     FiniteElementSpace* c_h1_fespace =
-                        new FiniteElementSpace(&mesh, c_h1_fec, spaceDimension);
-                     FiniteElementSpace* f_h1_fespace =
-                        new FiniteElementSpace(&fineMesh, f_h1_fec,spaceDimension);
+                     FiniteElementSpace* c_fespace =
+                        new FiniteElementSpace(&mesh, c_fec, spaceDimension);
+                     FiniteElementSpace* f_fespace =
+                        new FiniteElementSpace(&fineMesh, f_fec,spaceDimension);
 
 
                      Operator* referenceOperator = nullptr;
 
                      if (geometric == 0)
                      {
-                        referenceOperator = new PRefinementTransferOperator(*c_h1_fespace,
-                                                                            *f_h1_fespace);
+                        referenceOperator = new PRefinementTransferOperator(*c_fespace,
+                                                                            *f_fespace);
                      }
                      else
                      {
                         OperatorPtr P(Operator::ANY_TYPE);
-                        f_h1_fespace->GetTransferOperator(*c_h1_fespace, P);
+                        f_fespace->GetTransferOperator(*c_fespace, P);
                         P.SetOperatorOwner(false);
                         referenceOperator = P.Ptr();
                      }
 
-                     TransferOperator testTransferOperator(*c_h1_fespace, *f_h1_fespace);
-                     GridFunction X(c_h1_fespace);
-                     GridFunction X_cmp(c_h1_fespace);
-                     GridFunction Y_exact(f_h1_fespace);
-                     GridFunction Y_std(f_h1_fespace);
-                     GridFunction Y_test(f_h1_fespace);
-
+                     TransferOperator testTransferOperator(*c_fespace, *f_fespace);
+                     GridFunction X(c_fespace);
+                     GridFunction X_cmp(c_fespace);
+                     GridFunction Y_exact(f_fespace);
+                     GridFunction Y_std(f_fespace);
+                     GridFunction Y_test(f_fespace);
+                     coeff_order = 1;
                      if (vectorspace == 0)
                      {
                         FunctionCoefficient funcCoeff(&coeff);
@@ -166,31 +203,25 @@ TEST_CASE("transfer")
                      Y_std -= Y_exact;
                      REQUIRE(Y_std.Norml2() < 1e-12 * Y_exact.Norml2());
 
-                     if (vectorspace == 0)
-                     {
-                        testTransferOperator.Mult(X, Y_test);
+                     testTransferOperator.Mult(X, Y_test);
 
-                        Y_test -= Y_exact;
-                        REQUIRE(Y_test.Norml2() < 1e-12 * Y_exact.Norml2());
-                     }
+                     Y_test -= Y_exact;
+                     REQUIRE(Y_test.Norml2() < 1e-12 * Y_exact.Norml2());
 
-                     if (vectorspace == 0)
-                     {
-                        referenceOperator->MultTranspose(Y_exact, X);
-                        testTransferOperator.MultTranspose(Y_exact, X_cmp);
+                     referenceOperator->MultTranspose(Y_exact, X);
+                     testTransferOperator.MultTranspose(Y_exact, X_cmp);
 
-                        X -= X_cmp;
-                        REQUIRE(X.Norml2() < 1e-12 * X_cmp.Norml2());
-                     }
+                     X -= X_cmp;
+                     REQUIRE(X.Norml2() < 1e-12 * X_cmp.Norml2());
 
                      delete referenceOperator;
-                     delete f_h1_fespace;
-                     delete c_h1_fespace;
+                     delete f_fespace;
+                     delete c_fespace;
                      if (geometric == 0)
                      {
-                        delete f_h1_fec;
+                        delete f_fec;
                      }
-                     delete c_h1_fec;
+                     delete c_fec;
                   }
                }
             }
@@ -199,6 +230,250 @@ TEST_CASE("transfer")
    }
 }
 
+
+TEST_CASE("variable_order_transfer")
+{
+   for (int vectorspace = 0; vectorspace <= 3; ++vectorspace)
+   {
+      for (dimension = 2; dimension <= 3; ++dimension)
+      {
+         for (int ne = 1; ne <= 3; ++ne)
+         {
+            for (int order = 1; order <= 4; order *= 2)
+            {
+               std::cout << "Testing variable order transfer:\n"
+                         << "  Vectorspace:    " << vectorspace << "\n"
+                         << "  Dimension:      " << dimension << "\n"
+                         << "  Elements:       " << std::pow(ne, dimension) << "\n"
+                         << "  Coarse order:   " << order << "\n";
+
+               Mesh mesh;
+               if (dimension == 2)
+               {
+                  Element::Type type = Element::QUADRILATERAL;
+                  mesh = Mesh::MakeCartesian2D(ne, ne, type, 1, 1.0, 1.0);
+               }
+               else
+               {
+                  Element::Type type = Element::HEXAHEDRON;
+                  mesh = Mesh::MakeCartesian3D(ne, ne, ne, type, 1.0, 1.0, 1.0);
+               }
+               FiniteElementCollection* c_fec = nullptr;
+               FiniteElementCollection* f_fec = nullptr;
+               if (vectorspace < 2)
+               {
+                  c_fec = new H1_FECollection(order, dimension);
+                  f_fec = new H1_FECollection(order, dimension);
+               }
+               else if (vectorspace == 2)
+               {
+                  c_fec = new ND_FECollection(order+1, dimension);
+                  f_fec = new ND_FECollection(order+1, dimension);
+               }
+               else
+               {
+                  c_fec = new RT_FECollection(order, dimension);
+                  f_fec = new RT_FECollection(order, dimension);
+               }
+
+               mesh.EnsureNCMesh();
+               mesh.RandomRefinement(0.5);
+
+               int spaceDimension = 1;
+
+               if (vectorspace == 1)
+               {
+                  spaceDimension = dimension;
+               }
+
+               FiniteElementSpace* c_fespace =
+                  new FiniteElementSpace(&mesh, c_fec, spaceDimension);
+               FiniteElementSpace* f_fespace =
+                  new FiniteElementSpace(&mesh, f_fec,spaceDimension);
+               int maxorder = RandomPRefinement(*f_fespace);
+
+               std::cout  << "  Max fine order: " << maxorder << "\n";
+
+               Operator* referenceOperator = nullptr;
+
+               referenceOperator = new PRefinementTransferOperator(*c_fespace,
+                                                                   *f_fespace);
+
+               TransferOperator testTransferOperator(*c_fespace, *f_fespace);
+               GridFunction X(c_fespace); X = 0.;
+               GridFunction X_cmp(c_fespace); X_cmp = 0.;
+               GridFunction Y_exact(f_fespace); Y_exact = 0.;
+               GridFunction Y_std(f_fespace); Y_std = 0.;
+               GridFunction Y_test(f_fespace); Y_test = 0.;
+               coeff_order = std::min(2,order);
+               if (vectorspace == 0)
+               {
+                  FunctionCoefficient funcCoeff(&coeff);
+                  X.ProjectCoefficient(funcCoeff);
+                  Y_exact.ProjectCoefficient(funcCoeff);
+               }
+               else
+               {
+                  VectorFunctionCoefficient funcCoeff(dimension, &vectorcoeff);
+                  X.ProjectCoefficient(funcCoeff);
+                  Y_exact.ProjectCoefficient(funcCoeff);
+               }
+
+               Y_std = 0.0;
+               Y_test = 0.0;
+
+               referenceOperator->Mult(X, Y_std);
+               Y_std -= Y_exact;
+               REQUIRE(Y_std.Norml2() < 1e-12 * Y_exact.Norml2());
+
+               testTransferOperator.Mult(X, Y_test);
+
+               Y_test -= Y_exact;
+               REQUIRE(Y_test.Norml2() < 1e-12 * Y_exact.Norml2());
+
+               referenceOperator->MultTranspose(Y_exact, X);
+               testTransferOperator.MultTranspose(Y_exact, X_cmp);
+
+               X -= X_cmp;
+               REQUIRE(X.Norml2() < 1e-12 * X_cmp.Norml2());
+
+               delete referenceOperator;
+               delete f_fespace;
+               delete c_fespace;
+               delete f_fec;
+               delete c_fec;
+            }
+         }
+      }
+   }
+}
+
+TEST_CASE("variable_order_true_transfer")
+{
+   for (int vectorspace = 0; vectorspace <= 1; ++vectorspace)
+   {
+      for (dimension = 2; dimension <= 3; ++dimension)
+      {
+         for (int order = 1; order <= 3; order++)
+         {
+            std::cout << "Testing variable order true transfer:\n"
+                      << "  Vectorspace:    " << vectorspace << "\n"
+                      << "  Dimension:      " << dimension << "\n"
+                      << "  Coarse order:   " << order << "\n";
+
+            Mesh mesh;
+            if (dimension == 2)
+            {
+               Element::Type type = Element::QUADRILATERAL;
+               mesh = Mesh::MakeCartesian2D(3, 3, type, 1, 1.0, 1.0);
+            }
+            else
+            {
+               Element::Type type = Element::HEXAHEDRON;
+               mesh = Mesh::MakeCartesian3D(3, 3, 3, type, 1.0, 1.0, 1.0);
+            }
+            FiniteElementCollection* c_fec = nullptr;
+            FiniteElementCollection* f_fec = nullptr;
+            c_fec = new H1_FECollection(order, dimension);
+            f_fec = new H1_FECollection(order, dimension);
+            mesh.EnsureNCMesh();
+            mesh.RandomRefinement(0.5);
+            int spaceDimension = 1;
+
+            if (vectorspace == 1)
+            {
+               spaceDimension = dimension;
+            }
+
+            FiniteElementSpace* c_fespace =
+               new FiniteElementSpace(&mesh, c_fec, spaceDimension);
+            FiniteElementSpace* f_fespace =
+               new FiniteElementSpace(&mesh, f_fec,spaceDimension);
+
+            int maxorder = RandomPRefinement(*f_fespace);
+
+            std::cout  << "  Max fine order: " << maxorder << "\n";
+
+            const SparseMatrix * Rc = c_fespace->GetRestrictionMatrix();
+            TrueTransferOperator T(*c_fespace, *f_fespace);
+            GridFunction xc(c_fespace);
+            Vector Xc(c_fespace->GetTrueVSize());
+            Vector Diff(c_fespace->GetTrueVSize());
+            Vector Yc(c_fespace->GetTrueVSize());
+            Vector Xf(f_fespace->GetTrueVSize());
+            Vector Yf(f_fespace->GetTrueVSize());
+
+            coeff_order = 2;
+
+            if (vectorspace == 0)
+            {
+               FunctionCoefficient funcCoeff(&coeff);
+               xc.ProjectCoefficient(funcCoeff);
+            }
+            else
+            {
+               VectorFunctionCoefficient funcCoeff(dimension, &vectorcoeff);
+               xc.ProjectCoefficient(funcCoeff);
+            }
+            if (Rc)
+            {
+               Rc->Mult(xc,Xc);
+            }
+            else
+            {
+               Xc.MakeRef(xc,0);
+            }
+            T.Mult(Xc, Xf);
+
+            BilinearFormIntegrator * massc=nullptr;
+            BilinearFormIntegrator * massf=nullptr;
+
+            switch (vectorspace)
+            {
+               case 0:
+                  massc = new MassIntegrator;
+                  massf = new MassIntegrator;
+                  break;
+               default:
+                  massc = new VectorMassIntegrator;
+                  massf = new VectorMassIntegrator;
+                  break;
+            }
+
+            BilinearForm mc(c_fespace);
+            mc.AddDomainIntegrator(massc);
+            mc.Assemble();
+            SparseMatrix Mc;
+            Array<int> empty;
+            mc.FormSystemMatrix(empty, Mc);
+
+            BilinearForm mf(f_fespace);
+            mf.AddDomainIntegrator(massf);
+            mf.Assemble();
+            SparseMatrix Mf;
+            mf.FormSystemMatrix(empty, Mf);
+
+            Mf.Mult(Xf,Yf);
+
+            T.MultTranspose(Yf,Yc);
+
+            GSSmoother M(Mc);
+            Diff = 0.;
+            PCG(Mc, M, Yc, Diff, 0, 500, 1e-24, 0.0);
+
+            Diff -= Xc;
+
+            REQUIRE(Diff.Norml2() < 1e-10);
+            delete f_fespace;
+            delete c_fespace;
+            delete f_fec;
+            delete c_fec;
+         }
+      }
+   }
+}
+
+
 #ifdef MFEM_USE_MPI
 
 TEST_CASE("partransfer", "[Parallel]")
@@ -230,6 +505,7 @@ TEST_CASE("partransfer", "[Parallel]")
                                << "  Fine order:   " << fineOrder << "\n"
                                << "  Geometric:    " << geometric << "\n";
                   }
+                  coeff_order = 1;
 
                   Mesh mesh;
                   if (dimension == 2)
@@ -310,7 +586,6 @@ TEST_CASE("partransfer", "[Parallel]")
                   ParGridFunction X(c_h1_fespace);
                   ParGridFunction Y_exact(f_h1_fespace);
                   ParGridFunction Y(f_h1_fespace);
-
                   FunctionCoefficient funcCoeff(&coeff);
                   X.ProjectCoefficient(funcCoeff);
                   Y_exact.ProjectCoefficient(funcCoeff);