Merge branch 'master' into gf-conv-dev

# Conflicts:
#	fem/fe.cpp

examples/ex1p_vo.cpp

@@ -1,348 +0,0 @@
-//                       MFEM Example 1 - Parallel Version
-//
-// Compile with: make ex1p
-//
-// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
-//               mpirun -np 4 ex1p -m ../data/star.mesh
-//               mpirun -np 4 ex1p -m ../data/star-mixed.mesh
-//               mpirun -np 4 ex1p -m ../data/escher.mesh
-//               mpirun -np 4 ex1p -m ../data/fichera.mesh
-//               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
-//               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
-//               mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
-//               mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
-//               mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
-//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
-//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
-//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
-//               mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
-//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
-//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
-//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
-//               mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
-//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
-//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
-//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
-//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
-//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
-//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
-//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
-//
-// Device sample runs:
-//               mpirun -np 4 ex1p -pa -d cuda
-//               mpirun -np 4 ex1p -pa -d occa-cuda
-//               mpirun -np 4 ex1p -pa -d raja-omp
-//               mpirun -np 4 ex1p -pa -d ceed-cpu
-//               mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
-//             * mpirun -np 4 ex1p -pa -d ceed-cuda
-//             * mpirun -np 4 ex1p -pa -d ceed-hip
-//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
-//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
-//
-// Description:  This example code demonstrates the use of MFEM to define a
-//               simple finite element discretization of the Laplace problem
-//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
-//               Specifically, we discretize using a FE space of the specified
-//               order, or if order < 1 using an isoparametric/isogeometric
-//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
-//               NURBS mesh, etc.)
-//
-//               The example highlights the use of mesh refinement, finite
-//               element grid functions, as well as linear and bilinear forms
-//               corresponding to the left-hand side and right-hand side of the
-//               discrete linear system. We also cover the explicit elimination
-//               of essential boundary conditions, static condensation, and the
-//               optional connection to the GLVis tool for visualization.
-
-#include "mfem.hpp"
-#include "linalg/vector_operator.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-class CoordCoefficient : public Coefficient
-{
-private:
-   int d;
-
-   mutable Vector x;
-
-public:
-   CoordCoefficient(int d) : d(d), x(3) {}
-
-   double Eval(ElementTransformation &T, const IntegrationPoint &ip)
-   {
-      if (d == -1) { return 1.0; }
-
-      T.Transform(ip, x);
-      return x[d];
-   }
-};
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   MPI_Session mpi;
-   int num_procs = mpi.WorldSize();
-   int myid = mpi.WorldRank();
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int order = 1;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = true;
-   bool algebraic_ceed = false;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) or -1 for"
-                  " isoparametric space.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-#ifdef MFEM_USE_CEED
-   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
-                  "-no-a", "--no-algebraic",
-                  "Use algebraic Ceed solver");
-#endif
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 3. Enable hardware devices such as GPUs, and programming models such as
-   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
-   Device device(device_config);
-   if (myid == 0) { device.Print(); }
-
-   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-
-   // 5. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement. We choose
-   //    'ref_levels' to be the largest number that gives a final mesh with no
-   //    more than 10,000 elements.
-   {
-      int ref_levels =
-         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-   {
-      int par_ref_levels = 2;
-      for (int l = 0; l < par_ref_levels; l++)
-      {
-         pmesh.UniformRefinement();
-      }
-   }
-
-   // 7. Define a parallel finite element space on the parallel mesh. Here we
-   //    use continuous Lagrange finite elements of the specified order. If
-   //    order < 1, we instead use an isoparametric/isogeometric space.
-   FiniteElementCollection *fec;
-   bool delete_fec;
-   if (order > 0)
-   {
-      fec = new H1_FECollection(order, dim);
-      delete_fec = true;
-   }
-   else if (pmesh.GetNodes())
-   {
-      fec = pmesh.GetNodes()->OwnFEC();
-      delete_fec = false;
-      if (myid == 0)
-      {
-         cout << "Using isoparametric FEs: " << fec->Name() << endl;
-      }
-   }
-   else
-   {
-      fec = new H1_FECollection(order = 1, dim);
-      delete_fec = true;
-   }
-   ParFiniteElementSpace fespace(&pmesh, fec);
-   HYPRE_BigInt size = fespace.GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of finite element unknowns: " << size << endl;
-   }
-
-   // 8. Determine the list of true (i.e. parallel conforming) essential
-   //    boundary dofs. In this example, the boundary conditions are defined
-   //    by marking all the boundary attributes from the mesh as essential
-   //    (Dirichlet) and converting them to a list of true dofs.
-   Array<int> ess_tdof_list;
-   if (pmesh.bdr_attributes.Size())
-   {
-      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
-      ess_bdr = 1;
-      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   // 9. Set up the parallel linear form b(.) which corresponds to the
-   //    right-hand side of the FEM linear system, which in this case is
-   //    (1,phi_i) where phi_i are the basis functions in fespace.
-   ParLinearForm b(&fespace);
-   ConstantCoefficient one(1.0);
-   b.AddDomainIntegrator(new DomainLFIntegrator(one));
-   b.Assemble();
-
-   // 10. Define the solution vector x as a parallel finite element grid
-   //     function corresponding to fespace. Initialize x with initial guess of
-   //     zero, which satisfies the boundary conditions.
-   ParGridFunction x(&fespace);
-   x = 0.0;
-
-   ParVectorOperator vo(MPI_COMM_WORLD, myid, fespace.TrueVSize(), dim + 1);
-   {
-      for (int d=0; d <= dim; d++)
-      {
-         ParLinearForm bd(&fespace);
-         CoordCoefficient dCoef(d - 1);
-         bd.AddDomainIntegrator(new DomainLFIntegrator(dCoef));
-         bd.Assemble();
-
-         Vector *dv = new Vector(fespace.TrueVSize());
-         bd.ParallelAssemble(*dv);
-
-         vo.SetVector(d, dv, 1.0, true);
-      }
-   }
-
-   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
-   //     corresponding to the Laplacian operator -Delta, by adding the
-   //     Diffusion domain integrator.
-   ParBilinearForm a(&fespace);
-   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   a.AddDomainIntegrator(new DiffusionIntegrator(one));
-
-   // 12. Assemble the parallel bilinear form and the corresponding linear
-   //     system, applying any necessary transformations such as: parallel
-   //     assembly, eliminating boundary conditions, applying conforming
-   //     constraints for non-conforming AMR, static condensation, etc.
-   if (static_cond) { a.EnableStaticCondensation(); }
-   a.Assemble();
-
-   OperatorPtr A;
-   Vector B, X;
-   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
-
-   // 13. Solve the linear system A X = B.
-   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
-   //     * With partial assembly, use Jacobi smoothing, for now.
-   Solver *prec = NULL;
-   if (pa)
-   {
-      if (UsesTensorBasis(fespace))
-      {
-         if (algebraic_ceed)
-         {
-            prec = new ceed::AlgebraicSolver(a, ess_tdof_list);
-         }
-         else
-         {
-            prec = new OperatorJacobiSmoother(a, ess_tdof_list);
-         }
-      }
-   }
-   else
-   {
-      prec = new HypreBoomerAMG;
-   }
-   CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
-   cg.SetPrintLevel(1);
-   if (prec) { cg.SetPreconditioner(*prec); }
-   cg.SetOperator(*A);
-   cg.Mult(B, X);
-   delete prec;
-
-   {
-      Vector com((myid == 0) ? dim+1 : 0);
-      vo.Mult(X, com);
-
-      if (myid == 0)
-      {
-         cout << "Mass: " << com[0] << endl;
-         cout << "Center of mass: (";
-         for (int d=1; d<=dim; d++)
-         {
-            cout << com[d]/com[0];
-            if (d < dim) { cout << " ,"; }
-         }
-         cout << ")" << endl;
-      }
-   }
-
-   // 14. Recover the parallel grid function corresponding to X. This is the
-   //     local finite element solution on each processor.
-   a.RecoverFEMSolution(X, b, x);
-
-   // 15. Save the refined mesh and the solution in parallel. This output can
-   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
-   {
-      ostringstream mesh_name, sol_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "sol." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh.Print(mesh_ofs);
-
-      ofstream sol_ofs(sol_name.str().c_str());
-      sol_ofs.precision(8);
-      x.Save(sol_ofs);
-   }
-
-   // 16. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock(vishost, visport);
-      sol_sock << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock.precision(8);
-      sol_sock << "solution\n" << pmesh << x << flush;
-   }
-
-   // 17. Free the used memory.
-   if (delete_fec)
-   {
-      delete fec;
-   }
-
-   return 0;
-}

examples/ex_sconv_p.cpp

@@ -0,0 +1,503 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double freq = 0.5, kappa;
+static int dim;
+
+double u_func(const Vector &);
+
+enum SCA_TYPE {INVALID_SCA_TYPE = -1,
+               H1_TYPE = 0,
+               L2_TYPE,
+               L2I_TYPE,
+               NUM_SCA_TYPES
+              };
+enum CONV_TYPE {INVALID_CONV_TYPE = -1,
+                PROJECTION = 0,
+                INTERPOLATION_OP,
+                SOLVE,
+                SOLVE_W_DBC,
+                NUM_CONV_TYPES
+               };
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
+                                   FiniteElementCollection &fec);
+string GetTypeName(SCA_TYPE type);
+string GetConvTypeName(CONV_TYPE type);
+string GetConvTypeShortName(CONV_TYPE type);
+
+void Projection(const ParGridFunction &v0, ParGridFunction &v1);
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1);
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   MPI_Session mpi(argc, argv);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order0 = 1;
+   int order1 = 1;
+   int type0 = 0;
+   int type1 = 1;
+   int conv_type = -1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order0, "-o0", "--initial-order",
+                  "Finite element order (polynomial degree) "
+                  "for initial field.");
+   args.AddOption(&order1, "-o1", "--final-order",
+                  "Finite element order (polynomial degree) "
+                  "for final field.");
+   args.AddOption(&type0, "-t0", "--initial-type",
+                  "Set the basis type for the initial field: "
+                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
+   args.AddOption(&type1, "-t1", "--final-type",
+                  "Set the basis type for the final field: "
+                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
+   args.AddOption(&conv_type, "-c", "--conversion-type",
+                  "Set the conversion scheme: "
+                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
+                  "3-Least Squares with BC, -1 loop over all.");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root()) { args.PrintUsage(cout); }
+      return 1;
+   }
+   if (mpi.Root()) { args.PrintOptions(cout); }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (mpi.Root()) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution (1 time by
+   //    default, or specified on the command line with -rp). Once the parallel
+   //    mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   FunctionCoefficient uCoef(u_func);
+
+   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
+
+   if (mpi.Root())
+   {
+      cout << "L2 Errors:" << endl;
+   }
+   int t0a = (type0 == -1) ? 0 : type0;
+   int t0b = (type0 == -1) ? NUM_SCA_TYPES : (type0+1);
+   for (int t0 = t0a; t0 < t0b; t0++)
+   {
+      FiniteElementCollection *fec0 = GetFECollection((SCA_TYPE)t0, order0);
+      ParFiniteElementSpace   *fes0 = GetFESpace((SCA_TYPE)t0, pmesh, *fec0);
+
+      ParGridFunction x0(fes0);
+      x0.ProjectCoefficient(uCoef);
+
+      double err0 = x0.ComputeL2Error(uCoef);
+      if (mpi.Root())
+      {
+         cout << "Initial " << GetTypeName((SCA_TYPE)t0)
+              << ": \t\t" << err0 << endl;
+      }
+
+      // nn. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")";
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock0(vishost, visport);
+         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
+                   << pmesh.GetMyRank() << '\n';
+         sol_sock0.precision(8);
+         sol_sock0 << "solution\n" << pmesh << x0
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
+                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
+                   << flush;
+      }
+
+      int t1a = (type1 == -1) ? 0 : type1;
+      int t1b = (type1 == -1) ? NUM_SCA_TYPES : (type1+1);
+      for (int t1 = t1a; t1 < t1b; t1++)
+      {
+         FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, order1);
+         ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, pmesh, *fec1);
+
+         ParGridFunction y1(fes1);
+
+         if (mpi.Root())
+         {
+            cout << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")"
+                 << " -> "
+                 << GetTypeName((SCA_TYPE)t1) << "(" << order1 << ")"
+                 << ":" << endl;
+         }
+
+         int c01a = (conv_type == -1) ? 0 : conv_type;
+         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
+         for (int c01 = c01a; c01 < c01b; c01++)
+         {
+            string cmnt = "";
+
+            switch ((CONV_TYPE)c01)
+            {
+               case PROJECTION:
+                  Projection(x0, y1);
+                  break;
+               case INTERPOLATION_OP:
+                  cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
+                         "(should match projection)" : "(not expected to succeed)";
+                  InterpolationOp(x0, y1);
+                  break;
+               case SOLVE:
+                  LeastSquares((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1);
+                  break;
+               case SOLVE_W_DBC:
+                  LeastSquaresBC((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1, uCoef);
+                  break;
+               default:
+                  y1 = 0.0;
+            }
+
+            double err1 = y1.ComputeL2Error(uCoef);
+            cout << GetConvTypeName((CONV_TYPE)c01)
+                 << "\t\t" << err1 << "\t" << cmnt << endl;
+
+            if (visualization)
+            {
+               ostringstream oss;
+               oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")" << " --"
+                   << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
+                   << GetTypeName((SCA_TYPE)t1)<< "(" << order1 << ")";
+               char vishost[] = "localhost";
+               int  visport   = 19916;
+               socketstream sol_sock1(vishost, visport);
+               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
+                         << pmesh.GetMyRank() << '\n';
+               sol_sock1.precision(8);
+               sol_sock1 << "solution\n" << pmesh << y1
+                         << "window_title '" << oss.str() << "'"
+                         << "window_geometry "
+                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
+                                  Ws * (t0 - t0a))
+                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
+                         << " " << Ww << " " << Wh
+                         << flush;
+            }
+         }
+         if (mpi.Root())
+         {
+            cout << endl;
+         }
+
+         delete fes1;
+         delete fec1;
+      }
+
+      delete fes0;
+      delete fec0;
+
+      if (t0 < t0b - 1)
+      {
+         char c;
+         if (mpi.Root())
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      if (mpi.Root())
+      {
+         cout << endl;
+      }
+   }
+
+   return 0;
+}
+
+double u_func(const Vector &x)
+{
+   double kx = kappa * x[0];
+   double ky = kappa * x[1];
+   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
+
+   // Add the gradient of a scalar function
+   return cos(kx) * cos(ky) * cos(kz);
+}
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return new H1_FECollection(p, dim);
+      case L2_TYPE:
+         return new L2_FECollection(p-1, dim);
+      case L2I_TYPE:
+         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
+                                    FiniteElement::INTEGRAL);
+      default:
+         return NULL;
+   }
+}
+
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
+                                   ParMesh &pmesh,
+                                   FiniteElementCollection &fec)
+{
+   return new ParFiniteElementSpace(&pmesh, &fec);
+}
+
+string GetTypeName(SCA_TYPE type)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return "    H1";
+      case L2_TYPE:
+         return "    L2";
+      case L2I_TYPE:
+         return "    L2I";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Projection            ";
+      case INTERPOLATION_OP:
+         return "Interpolation Operator";
+      case SOLVE:
+         return "Least Squares         ";
+      case SOLVE_W_DBC:
+         return "Least Squares with BC ";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeShortName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Proj";
+      case INTERPOLATION_OP:
+         return "Interp";
+      case SOLVE:
+         return "LS";
+      case SOLVE_W_DBC:
+         return "LSwBC";
+      default:
+         return "--";
+   }
+}
+
+/** Perform a naive projection from one scalar field to another.
+
+   This scheme simply evaluates v0 at the interpolation points of v1.
+
+   If v0 has reduced continuity compared to v1 this can produce
+   results that depend on the order in which the elements are
+   traversed.
+
+   Suitable conversions:
+   H1 -> L2
+   H1 -> DG (same as L2)
+ */
+void Projection(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   GridFunctionCoefficient v0Coef(&v0);
+   v1.ProjectCoefficient(v0Coef);
+}
+
+/** In theory this interpolation scheme should be equivalent to projection.
+
+    Building an interpolastion matrix could lead to computational
+    efficiency compared to simple projection if the operator will be
+    used several times.
+
+    Unfortunately this is broken for several combinations of source
+    and target fields.
+*/
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
+   op.AddDomainInterpolator(new IdentityInterpolator);
+   op.Assemble();
+   op.Finalize();
+
+   op.Mult(v0, v1);
+}
+
+/** Compute a least-squares best fit using the target basis functions.
+
+    This scheme is more difficult to setup and more computationally
+    expensive but the results can be significantly better than simple
+    projections.
+*/
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   HypreParMatrix * M = m.ParallelAssemble();
+
+   HypreDiagScale diag(*M);
+   HyprePCG pcg(*M);
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   Vector B, X;
+   b.ParallelAssemble(B);
+
+   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+
+   delete M;
+}
+
+/** Compute a least-squares best fit with boundary conditions.
+
+    This scheme is virtually identical to the previous one but it
+    makes use of boundary values, when available, to improve the
+    accuracy.  This scheme can produce significantly better results
+    when the normal derivative of the field is large near the
+    boundary.  This is particularly true when the field is
+    under-resolved near the boundary.
+*/
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
+   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
+      ess_bdr = 1;
+      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   if (t1 == H1_TYPE)
+   {
+      v1.ProjectBdrCoefficient(c, ess_bdr);
+   }
+
+   OperatorPtr M;
+   Vector B, X;
+   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
+
+   HypreDiagScale diag(*M.As<HypreParMatrix>());
+   HyprePCG pcg(*M.As<HypreParMatrix>());
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+}

examples/ex_vconv_p.cpp

@@ -0,0 +1,586 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double freq = 0.5, kappa;
+static int dim;
+
+void u_func(const Vector &, Vector &);
+
+enum VEC_TYPE {INVALID_VEC_TYPE = -1,
+               H1V_TYPE = 0,
+               ND_TYPE,
+               RT_TYPE,
+               L2V_TYPE,
+               NUM_VEC_TYPES
+              };
+enum CONV_TYPE {INVALID_CONV_TYPE = -1,
+                PROJECTION = 0,
+                INTERPOLATION_OP,
+                SOLVE,
+                SOLVE_W_DBC,
+                NUM_CONV_TYPES
+               };
+
+FiniteElementCollection * GetFECollection(VEC_TYPE type, int p);
+ParFiniteElementSpace * GetFESpace(VEC_TYPE type, ParMesh &pmesh,
+                                   FiniteElementCollection &fec);
+string GetTypeName(VEC_TYPE type);
+string GetConvTypeName(CONV_TYPE type);
+
+void Projection(const ParGridFunction &v0, ParGridFunction &v1);
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
+void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
+                  VEC_TYPE t1, ParGridFunction &v1);
+void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
+                    VEC_TYPE t1, ParGridFunction &v1,
+                    VectorCoefficient &vc);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   MPI_Session mpi(argc, argv);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order0 = 1;
+   int order1 = 1;
+   int type0 = 0;
+   int type1 = 1;
+   int conv_type = -1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order0, "-o0", "--initial-order",
+                  "Finite element order (polynomial degree) "
+                  "for initial field.");
+   args.AddOption(&order1, "-o1", "--final-order",
+                  "Finite element order (polynomial degree) "
+                  "for final field.");
+   args.AddOption(&type0, "-t0", "--initial-type",
+                  "Set the basis type for the initial field: "
+                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
+   args.AddOption(&type1, "-t1", "--final-type",
+                  "Set the basis type for the final field: "
+                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
+   args.AddOption(&conv_type, "-c", "--conversion-type",
+                  "Set the conversion scheme: "
+                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
+                  "3-Least Squares with BC, -1 loop over all.");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root()) { args.PrintUsage(cout); }
+      return 1;
+   }
+   if (mpi.Root()) { args.PrintOptions(cout); }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (mpi.Root()) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution (1 time by
+   //    default, or specified on the command line with -rp). Once the parallel
+   //    mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   VectorFunctionCoefficient uCoef(dim, u_func);
+
+   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
+
+   if (mpi.Root())
+   {
+      cout << "L2 Errors:" << endl;
+   }
+   int t0a = (type0 == -1) ? 0 : type0;
+   int t0b = (type0 == -1) ? NUM_VEC_TYPES : (type0+1);
+   for (int t0 = t0a; t0 < t0b; t0++)
+   {
+      FiniteElementCollection *fec0 = GetFECollection((VEC_TYPE)t0, order0);
+      ParFiniteElementSpace   *fes0 = GetFESpace((VEC_TYPE)t0, pmesh, *fec0);
+
+      ParGridFunction x0(fes0);
+      x0.ProjectCoefficient(uCoef);
+
+      double err0 = x0.ComputeL2Error(uCoef);
+      if (mpi.Root())
+      {
+         cout << "Initial " << GetTypeName((VEC_TYPE)t0)
+              << ": \t\t" << err0 << endl;
+      }
+
+      // nn. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << GetTypeName((VEC_TYPE)t0);
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock0(vishost, visport);
+         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
+                   << pmesh.GetMyRank() << '\n';
+         sol_sock0.precision(8);
+         sol_sock0 << "solution\n" << pmesh << x0
+                   << "keys vvv "
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
+                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
+                   << flush;
+      }
+
+      int t1a = (type1 == -1) ? 0 : type1;
+      int t1b = (type1 == -1) ? NUM_VEC_TYPES : (type1+1);
+      for (int t1 = t1a; t1 < t1b; t1++)
+      {
+         FiniteElementCollection *fec1 = GetFECollection((VEC_TYPE)t1, order1);
+         ParFiniteElementSpace   *fes1 = GetFESpace((VEC_TYPE)t1, pmesh, *fec1);
+
+         ParGridFunction x1(fes1);
+
+         if (mpi.Root())
+         {
+            cout << GetTypeName((VEC_TYPE)t0) << " -> "
+                 << GetTypeName((VEC_TYPE)t1) << ":" << endl;
+         }
+
+         int c01a = (conv_type == -1) ? 0 : conv_type;
+         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
+         for (int c01 = c01a; c01 < c01b; c01++)
+         {
+            switch ((CONV_TYPE)c01)
+            {
+               case PROJECTION:
+                  Projection(x0, x1);
+                  break;
+               case INTERPOLATION_OP:
+                  // InterpolationOp(x0, x1);
+                  x1 = 0.0;
+                  break;
+               case SOLVE:
+                  LeastSquares((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1);
+                  break;
+               case SOLVE_W_DBC:
+                  LeastSquaresBC((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1, uCoef);
+                  break;
+               default:
+                  x1 = 0.0;
+            }
+
+            double err1 = x1.ComputeL2Error(uCoef);
+            cout << GetConvTypeName((CONV_TYPE)c01)
+                 << "\t\t" << err1 << endl;
+
+            if (visualization)
+            {
+               ostringstream oss;
+               oss << GetTypeName((VEC_TYPE)t0) << " --" << c01 << "--> "
+                   << GetTypeName((VEC_TYPE)t1);
+               char vishost[] = "localhost";
+               int  visport   = 19916;
+               socketstream sol_sock1(vishost, visport);
+               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
+                         << pmesh.GetMyRank() << '\n';
+               sol_sock1.precision(8);
+               sol_sock1 << "solution\n" << pmesh << x1
+                         << "keys vvv "
+                         << "window_title '" << oss.str() << "'"
+                         << "window_geometry "
+                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
+                                  Ws * (t0 - t0a))
+                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
+                         << " " << Ww << " " << Wh
+                         << flush;
+            }
+         }
+         if (mpi.Root())
+         {
+            cout << endl;
+         }
+
+         delete fes1;
+         delete fec1;
+      }
+
+      delete fes0;
+      delete fec0;
+
+      if (t0 < t0b - 1)
+      {
+         char c;
+         if (mpi.Root())
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      if (mpi.Root())
+      {
+         cout << endl;
+      }
+   }
+
+   return 0;
+}
+
+void u_func(const Vector &x, Vector &u)
+{
+   u.SetSize(dim);
+
+   double kx = kappa * x[0];
+   double ky = kappa * x[1];
+   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
+
+   // Add the gradient of a scalar function
+   u(0) = sin(kx) * cos(ky);
+   u(1) = cos(kx) * sin(ky);
+   if (dim == 3)
+   {
+      u(0) *= cos(kz);
+      u(1) *= cos(kz);
+      u(2) = cos(kx) * cos(ky) * sin(kz);
+   }
+
+   // Add the curl of a vector function
+   u(0) -= cos(kx) * sin(ky);
+   u(1) += sin(kx) * cos(ky);
+   if (dim == 3)
+   {
+      u(0) += cos(kx) * sin(kz);
+      u(1) -= cos(ky) * sin(kz);
+      u(2) += (sin(ky) - sin(kx)) * cos(kz);
+   }
+}
+
+FiniteElementCollection * GetFECollection(VEC_TYPE type, int p)
+{
+   switch (type)
+   {
+      case H1V_TYPE:
+         return new H1_FECollection(p, dim);
+      case ND_TYPE:
+         return new ND_FECollection(p, dim);
+      case RT_TYPE:
+         return new RT_FECollection(p-1, dim);
+      case L2V_TYPE:
+         return new L2_FECollection(p-1, dim);
+      default:
+         return NULL;
+   }
+}
+
+ParFiniteElementSpace * GetFESpace(VEC_TYPE type,
+                                   ParMesh &pmesh,
+                                   FiniteElementCollection &fec)
+{
+   switch (type)
+   {
+      case H1V_TYPE:
+      case L2V_TYPE:
+         return new ParFiniteElementSpace(&pmesh, &fec, dim);
+      case ND_TYPE:
+      case RT_TYPE:
+         return new ParFiniteElementSpace(&pmesh, &fec);
+      default:
+         return NULL;
+   }
+}
+
+string GetTypeName(VEC_TYPE type)
+{
+   switch (type)
+   {
+      case H1V_TYPE:
+         return "    H1V";
+      case ND_TYPE:
+         return "H(Curl)";
+      case RT_TYPE:
+         return " H(Div)";
+      case L2V_TYPE:
+         return "    L2V";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Projection            ";
+      case INTERPOLATION_OP:
+         return "Interpolation Operator";
+      case SOLVE:
+         return "Least Squares         ";
+      case SOLVE_W_DBC:
+         return "Least Squares with BC ";
+      default:
+         return "--";
+   }
+}
+
+/** Perform a naive projection from one vector field to another.
+
+   This scheme simply evaluates v0 at the interpolation points of v1.
+
+   If v0 has reduced continuity compared to v1 this can produce
+   results that depend on the order in which the elements are
+   traversed.
+
+   Suitable conversions:
+   H1V     -> H(Curl), H(Div), or L2V
+   H(Curl) -> L2V
+   H(Div)  -> L2V
+
+ */
+void Projection(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   VectorGridFunctionCoefficient v0Coef(&v0);
+   v1.ProjectCoefficient(v0Coef);
+}
+
+/** In theory this interpolation scheme should be equivalent to projection.
+
+    Building an interpolastion matrix could lead to computational
+    efficiency compared to simple projection if the operator will be
+    used several times.
+
+    Unfortunately this is broken for several combinations of source
+    and target fields.
+*/
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
+   op.AddDomainInterpolator(new IdentityInterpolator);
+   op.Assemble();
+   op.Finalize();
+
+   op.Mult(v0, v1);
+}
+
+/** Compute a least-squares best fit using the target basis functions.
+
+    This scheme is more difficult to setup and more computationally
+    expensive but the results can be significantly better than simple
+    projections.
+*/
+void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
+                  VEC_TYPE t1, ParGridFunction &v1)
+{
+   bool trans = false;
+   ParFiniteElementSpace *fes0, *fes1;
+   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
+       (t1 == ND_TYPE || t1 == RT_TYPE))
+   {
+      fes0 = v1.ParFESpace();
+      fes1 = v0.ParFESpace();
+      trans = true;
+   }
+   else
+   {
+      fes0 = v0.ParFESpace();
+      fes1 = v1.ParFESpace();
+   }
+   ParMixedBilinearForm op(fes0, fes1);
+   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      op.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      op.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   if (trans)
+   {
+      op.MultTranspose(v0, b);
+   }
+   else
+   {
+      op.Mult(v0, b);
+   }
+
+   ParBilinearForm m(v1.ParFESpace());
+   if (t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      m.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      m.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   m.Assemble();
+   m.Finalize();
+
+   HypreParMatrix * M = m.ParallelAssemble();
+
+   HypreDiagScale diag(*M);
+   HyprePCG pcg(*M);
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   Vector B, X;
+   b.ParallelAssemble(B);
+
+   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+
+   delete M;
+}
+
+/** Compute a least-squares best fit with boundary conditions.
+
+    This scheme is virtually identical to the previous one but it
+    makes use of boundary values, when available, to improve the
+    accuracy.  This scheme can produce significantly better results
+    when the normal derivative of the field is large near the
+    boundary.  This is particularly true when the field is
+    under-resolved near the boundary.
+*/
+void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
+                    VEC_TYPE t1, ParGridFunction &v1,
+                    VectorCoefficient &vc)
+{
+   bool trans = false;
+   ParFiniteElementSpace *fes0, *fes1;
+   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
+       (t1 == ND_TYPE || t1 == RT_TYPE))
+   {
+      fes0 = v1.ParFESpace();
+      fes1 = v0.ParFESpace();
+      trans = true;
+   }
+   else
+   {
+      fes0 = v0.ParFESpace();
+      fes1 = v1.ParFESpace();
+   }
+   ParMixedBilinearForm op(fes0, fes1);
+   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      op.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      op.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   if (trans)
+   {
+      op.MultTranspose(v0, b);
+   }
+   else
+   {
+      op.Mult(v0, b);
+   }
+
+   ParBilinearForm m(v1.ParFESpace());
+   if (t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      m.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      m.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   m.Assemble();
+   m.Finalize();
+
+   Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
+   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
+      ess_bdr = 1;
+      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   if (t1 == H1V_TYPE)
+   {
+      v1.ProjectBdrCoefficient(vc, ess_bdr);
+   }
+   if (t1 == ND_TYPE)
+   {
+      v1.ProjectBdrCoefficientTangent(vc, ess_bdr);
+   }
+   else if (t1 == RT_TYPE)
+   {
+      v1.ProjectBdrCoefficientNormal(vc, ess_bdr);
+   }
+
+   OperatorPtr M;
+   Vector B, X;
+   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
+
+   HypreDiagScale diag(*M.As<HypreParMatrix>());
+   HyprePCG pcg(*M.As<HypreParMatrix>());
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+}

linalg/vector_operator.cpp

@@ -1,80 +0,0 @@
-// Copyright (c) 2010-2022, 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 "vector_operator.hpp"
-
-namespace mfem
-{
-
-#ifdef MFEM_USE_MPI
-
-ParVectorOperator::ParVectorOperator(MPI_Comm comm,
-                                     int myid,
-                                     int local_vec_size,
-                                     int num_vecs)
-   : Operator((myid == 0) ? num_vecs : 0, local_vec_size),
-     comm(comm),
-     myid(myid),
-     vecs(num_vecs),
-     coefs(num_vecs),
-     owns(num_vecs)
-{
-   vecs = NULL;
-   coefs = 1.0;
-   owns = false;
-}
-
-ParVectorOperator::~ParVectorOperator()
-{
-   for (int i=0; i < vecs.Size(); i++)
-   {
-      if (owns[i]) { delete vecs[i]; }
-      vecs[i] = NULL;
-   }
-}
-
-void ParVectorOperator::SetVector(int idx, Vector *vec,
-                                  double c, bool own_vec)
-{
-   MFEM_VERIFY(idx >= 0 && idx < vecs.Size(),
-               "ParVectorOperator: Index out of range");
-
-   vecs[idx] = vec;
-   coefs[idx] = c;
-   owns[idx] = own_vec;
-}
-
-void ParVectorOperator::Mult(const Vector &x, Vector &y) const
-{
-   for (int i=0; i<vecs.Size(); i++)
-   {
-      double vo = coefs[i] * (*vecs[i] * x);
-      double vi = 0.0;
-      MPI_Reduce(&vo, &vi, 1, MPI_DOUBLE, MPI_SUM, 0, comm);
-      if (myid == 0) { y[i] = vi; }
-   }
-}
-
-/// Action of the transpose operator: `y=A^t(x)`.
-void ParVectorOperator::MultTranspose(const Vector &x, Vector &y) const
-{
-   y = 0.0;
-   for (int i=0; i<vecs.Size(); i++)
-   {
-      double xi = (myid == 0) ? x[i] : 0.0;
-      MPI_Bcast(&xi, 1, MPI_DOUBLE, 0, comm);
-      y.Add(xi * coefs[i], *vecs[i]);
-   }
-}
-
-}
-
-#endif // MFEM_USE_MPI

linalg/vector_operator.hpp

@@ -1,55 +0,0 @@
-// 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_VECTOR_OPERATOR
-#define MFEM_VECTOR_OPERATOR
-
-#include "operator.hpp"
-#include "vector.hpp"
-
-namespace mfem
-{
-
-#ifdef MFEM_USE_MPI
-
-class ParVectorOperator : public Operator
-{
-private:
-   MPI_Comm comm;
-   int myid;
-
-   Array<Vector*> vecs;
-   Array<double> coefs;
-   Array<bool> owns;
-
-public:
-   ParVectorOperator(MPI_Comm comm,
-                     int myid,
-                     int local_vec_size,
-                     int num_vecs);
-
-   ~ParVectorOperator();
-
-   void SetVector(int idx, Vector *vec,
-                  double c = 1.0, bool own_vec = false);
-
-   /// Operator application: `y=A(x)`.
-   void Mult(const Vector &x, Vector &y) const;
-
-   /// Action of the transpose operator: `y=A^t(x)`.
-   void MultTranspose(const Vector &x, Vector &y) const;
-};
-
-#endif // MFEM_USE_MPI
-
-} // namespace mfem
-
-#endif // MFEM_VECTOR_OPERATOR

miniapps/tools/gf-conv.cpp

@@ -0,0 +1,577 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static int dim;
+
+enum SCA_TYPE {INVALID_SCA_TYPE = -1,
+               H1_TYPE = 0,
+               L2_TYPE,
+               L2I_TYPE,
+               NUM_SCA_TYPES
+              };
+enum CONV_TYPE {INVALID_CONV_TYPE = -1,
+                PROJECTION = 0,
+                INTERPOLATION_OP,
+                SOLVE,
+                SOLVE_W_DBC,
+                NUM_CONV_TYPES
+               };
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
+                                   FiniteElementCollection &fec);
+void parseFieldNames(const char * field_name_c_str,
+                     vector<string> &field_names);
+
+string GetTypeName(SCA_TYPE type);
+string GetTypeShortName(SCA_TYPE type);
+string GetConvTypeName(CONV_TYPE type);
+string GetConvTypeShortName(CONV_TYPE type);
+
+void Projection(const ParGridFunction &v0, ParGridFunction &v1);
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1);
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c);
+
+int main(int argc, char *argv[])
+{
+#ifdef MFEM_USE_MPI
+   Mpi::Init();
+   if (!Mpi::Root()) { mfem::out.Disable(); mfem::err.Disable(); }
+   Hypre::Init();
+#endif
+
+   // Parse command-line options.
+   const char *coll_name = NULL;
+   int cycle = 0;
+
+   const char *field_name_c_str = "ALL";
+
+   Array<int> orders;
+   Array<int> types;
+   Array<int> conv_types;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&coll_name, "-r", "--root-file",
+                  "Set the VisIt data collection root file prefix.", true);
+   args.AddOption(&cycle, "-c", "--cycle", "Set the cycle index to read.");
+   args.AddOption(&field_name_c_str, "-fn", "--field-names",
+                  "List of field names to get values from.");
+   args.AddOption(&orders, "-o", "--final-order",
+                  "Finite element orders for each final field "
+                  "(an array of integers for multiple fields).");
+   args.AddOption(&types, "-t", "--final-type",
+                  "Set the basis type for the final fields: "
+                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
+   args.AddOption(&conv_types, "-ct", "--conversion-type",
+                  "Set the conversion schemes: "
+                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
+                  "3-Least Squares with BC, -1 loop over all.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+#ifdef MFEM_USE_MPI
+   VisItDataCollection dc(MPI_COMM_WORLD, coll_name);
+#else
+   VisItDataCollection dc(coll_name);
+#endif
+
+   dc.Load(cycle);
+
+   if (dc.Error() != DataCollection::NO_ERROR)
+   {
+      mfem::out << "Error loading VisIt data collection: " << coll_name << endl;
+      return 1;
+   }
+
+   dim = dc.GetMesh()->Dimension();
+   int spaceDim = dc.GetMesh()->SpaceDimension();
+
+   mfem::out << endl;
+   mfem::out << "Collection Name:    " << dc.GetCollectionName() << endl;
+   mfem::out << "Manifold Dimension: " << dim << endl;
+   mfem::out << "Space Dimension:    " << spaceDim << endl;
+   mfem::out << "Cycle:              " << dc.GetCycle() << endl;
+   mfem::out << "Time:               " << dc.GetTime() << endl;
+   mfem::out << "Time Step:          " << dc.GetTimeStep() << endl;
+   mfem::out << endl;
+
+   typedef DataCollection::FieldMapType fields_t;
+   const fields_t &fields = dc.GetFieldMap();
+   // Print the names of all fields.
+   mfem::out << "fields: [ ";
+   for (fields_t::const_iterator it = fields.begin(); it != fields.end(); ++it)
+   {
+      if (it != fields.begin()) { mfem::out << ", "; }
+      mfem::out << it->first;
+   }
+   mfem::out << " ]" << endl;
+
+   // Parsing desired field names
+   vector<string> field_names;
+   parseFieldNames(field_name_c_str, field_names);
+
+   if (field_names.size() == 1)
+   {
+      if (field_names[0] == "ALL")
+      {
+         fields_t::const_iterator it = fields.begin();
+         field_names[0] = it->first; it++;
+         for ( ; it != fields.end(); ++it)
+         {
+            field_names.push_back(it->first);
+         }
+      }
+   }
+
+   if (orders.Size() < field_names.size())
+   {
+      int size = orders.Size();
+      int order = (size > 0) ? orders[0] : 1;
+
+      orders.SetSize(field_names.size());
+      for (int i=size; i < field_names.size(); i++)
+      {
+         orders[i] = order;
+      }
+   }
+
+   if (types.Size() < field_names.size())
+   {
+      int size = types.Size();
+      int type = (size > 0) ? types[0] : 0;
+
+      types.SetSize(field_names.size());
+      for (int i=size; i < field_names.size(); i++)
+      {
+         types[i] = type;
+      }
+   }
+
+   if (conv_types.Size() < field_names.size())
+   {
+      int size = conv_types.Size();
+      int type = (size > 0) ? conv_types[0] : 0;
+
+      conv_types.SetSize(field_names.size());
+      for (int i=size; i < field_names.size(); i++)
+      {
+         conv_types[i] = type;
+      }
+   }
+
+   // Print field names to be extracted
+   mfem::out << "Extracting fields: ";
+   for (int i=0; i < field_names.size(); i++)
+   {
+      mfem::out << " \"" << field_names[i] << "\"";
+   }
+   mfem::out << endl;
+
+#ifdef MFEM_USE_MPI
+   ParMesh *mesh = dynamic_cast<ParMesh*>(dc.GetMesh());
+#else
+   Mesh *mesh = dc.GetMesh();
+#endif
+   if (mesh == NULL)
+   {
+      mfem::out << "Problem with mesh\n";
+      return 1;
+   }
+
+   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
+
+   // Loop over all requested fields.
+   for (int i=0; i < field_names.size(); i++)
+   {
+#ifdef MFEM_USE_MPI
+      ParGridFunction *x0 = dc.GetParField(field_names[i]);
+#else
+      GridFunction *x0 = dc.GetField(field_names[i]);
+#endif
+      if (x0 == NULL)
+      {
+         mfem::out << "Problem with x0 for field \"" << field_names[i] << "\"\n";
+         continue;
+      }
+
+      int t0 = 0;
+
+      // nn. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << field_names[i];
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock0(vishost, visport);
+#ifdef MFEM_USE_MPI
+         sol_sock0 << "parallel " << mesh->GetNRanks() << ' '
+                   << mesh->GetMyRank() << '\n';
+#endif
+         sol_sock0.precision(8);
+         sol_sock0 << "solution\n" << *mesh << *x0
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << Ws * (t0) << " " << Ws * (t0) << " "
+                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
+                   << flush;
+      }
+
+      int t1 = types[i];
+      FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, orders[i]);
+      ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, *mesh, *fec1);
+
+      ParGridFunction *y1 = new ParGridFunction(fes1);
+
+      mfem::out << GetTypeName((SCA_TYPE)t1) << "(" << orders[i] << ")"
+                << ":" << endl;
+
+      int c01 = conv_types[i];
+      string cmnt = "";
+
+      switch ((CONV_TYPE)c01)
+      {
+         case PROJECTION:
+            Projection(*x0, *y1);
+            break;
+         case INTERPOLATION_OP:
+            cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
+                   "(should match projection)" : "(not expected to succeed)";
+            InterpolationOp(*x0, *y1);
+            break;
+         case SOLVE:
+            LeastSquares((SCA_TYPE)t0, *x0, (SCA_TYPE)t1, *y1);
+            break;
+         default:
+            *y1 = 0.0;
+      }
+
+      {
+         ostringstream oss;
+         oss << field_names[i] << "_" << GetConvTypeShortName((CONV_TYPE)c01)
+             << "_" << GetTypeShortName((SCA_TYPE)t1) << "_o" << orders[i];
+
+         dc.RegisterField(oss.str(), y1);
+      }
+
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
+             << GetTypeName((SCA_TYPE)t1)<< "(" << orders[i] << ")";
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock1(vishost, visport);
+#ifdef MFEM_USE_MPI
+         sol_sock1 << "parallel " << mesh->GetNRanks() << ' '
+                   << mesh->GetMyRank() << '\n';
+#endif
+         sol_sock1.precision(8);
+         sol_sock1 << "solution\n" << *mesh << y1
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << (int)((Ww + Fw) * (1.5 + c01) +
+                            Ws * (t0))
+                   << " " << (Wh + Fh) * (t1) + Ws * (t0)
+                   << " " << Ww << " " << Wh
+                   << flush;
+      }
+      mfem::out << endl;
+
+      // delete fes1;
+      // delete fec1;
+   }
+   dc.Save();
+
+   return 0;
+}
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return new H1_FECollection(p, dim);
+      case L2_TYPE:
+         return new L2_FECollection(p-1, dim);
+      case L2I_TYPE:
+         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
+                                    FiniteElement::INTEGRAL);
+      default:
+         return NULL;
+   }
+}
+
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
+                                   ParMesh &pmesh,
+                                   FiniteElementCollection &fec)
+{
+   return new ParFiniteElementSpace(&pmesh, &fec);
+}
+
+string GetTypeName(SCA_TYPE type)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return "    H1";
+      case L2_TYPE:
+         return "    L2";
+      case L2I_TYPE:
+         return "    L2I";
+      default:
+         return "--";
+   }
+}
+
+string GetTypeShortName(SCA_TYPE type)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return "H1";
+      case L2_TYPE:
+         return "L2";
+      case L2I_TYPE:
+         return "L2I";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Projection            ";
+      case INTERPOLATION_OP:
+         return "Interpolation Operator";
+      case SOLVE:
+         return "Least Squares         ";
+      case SOLVE_W_DBC:
+         return "Least Squares with BC ";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeShortName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Proj";
+      case INTERPOLATION_OP:
+         return "Interp";
+      case SOLVE:
+         return "LS";
+      case SOLVE_W_DBC:
+         return "LSwBC";
+      default:
+         return "--";
+   }
+}
+
+void parseFieldNames(const char * field_name_c_str, vector<string> &field_names)
+{
+   string field_name_str(field_name_c_str);
+   string field_name;
+
+   for (string::iterator it=field_name_str.begin();
+        it!=field_name_str.end(); it++)
+   {
+      if (*it == '\\')
+      {
+         it++;
+         field_name.push_back(*it);
+      }
+      else if (*it == ' ')
+      {
+         if (!field_name.empty())
+         {
+            field_names.push_back(field_name);
+         }
+         field_name.clear();
+      }
+      else if (it == field_name_str.end() - 1)
+      {
+         field_name.push_back(*it);
+         field_names.push_back(field_name);
+      }
+      else
+      {
+         field_name.push_back(*it);
+      }
+   }
+   if (field_names.size() == 0)
+   {
+      field_names.push_back("ALL");
+   }
+}
+
+/** Perform a naive projection from one scalar field to another.
+
+   This scheme simply evaluates v0 at the interpolation points of v1.
+
+   If v0 has reduced continuity compared to v1 this can produce
+   results that depend on the order in which the elements are
+   traversed.
+
+   Suitable conversions:
+   H1 -> L2
+   H1 -> DG (same as L2)
+ */
+void Projection(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   GridFunctionCoefficient v0Coef(&v0);
+   v1.ProjectCoefficient(v0Coef);
+}
+
+/** In theory this interpolation scheme should be equivalent to projection.
+
+    Building an interpolastion matrix could lead to computational
+    efficiency compared to simple projection if the operator will be
+    used several times.
+
+    Unfortunately this is broken for several combinations of source
+    and target fields.
+*/
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
+   op.AddDomainInterpolator(new IdentityInterpolator);
+   op.Assemble();
+   op.Finalize();
+
+   op.Mult(v0, v1);
+}
+
+/** Compute a least-squares best fit using the target basis functions.
+
+    This scheme is more difficult to setup and more computationally
+    expensive but the results can be significantly better than simple
+    projections.
+*/
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   HypreParMatrix * M = m.ParallelAssemble();
+
+   HypreDiagScale diag(*M);
+   HyprePCG pcg(*M);
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   Vector B, X;
+   b.ParallelAssemble(B);
+
+   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+
+   delete M;
+}
+
+/** Compute a least-squares best fit with boundary conditions.
+
+    This scheme is virtually identical to the previous one but it
+    makes use of boundary values, when available, to improve the
+    accuracy.  This scheme can produce significantly better results
+    when the normal derivative of the field is large near the
+    boundary.  This is particularly true when the field is
+    under-resolved near the boundary.
+*/
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
+   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
+      ess_bdr = 1;
+      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   if (t1 == H1_TYPE)
+   {
+      v1.ProjectBdrCoefficient(c, ess_bdr);
+   }
+
+   OperatorPtr M;
+   Vector B, X;
+   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
+
+   HypreDiagScale diag(*M.As<HypreParMatrix>());
+   HyprePCG pcg(*M.As<HypreParMatrix>());
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+}