Merge branch 'master' into gf-conv-dev

Improvements of Par(Mixed)BilinearForm and Par(Block)NonlinearForm

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);
+}

fem/nonlinearform.cpp

@@ -436,7 +436,7 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
    // In parallel, the result is in 'py' which is an alias for 'aux2'.
 }
 
-Operator &NonlinearForm::GetGradient(const Vector &x) const
+Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
 {
    if (ext)
    {
@@ -644,6 +644,8 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
       }
    }
 
+   if (!finalize) { return *Grad; }
+
    if (!Grad->Finalized())
    {
       Grad->Finalize(skip_zeros);
@@ -1203,7 +1205,14 @@ const BlockVector &BlockNonlinearForm::Prolongate(const BlockVector &bx) const
       aux1.Update(block_offsets);
       for (int s = 0; s < fes.Size(); s++)
       {
-         P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
+         if (P[s])
+         {
+            P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
+         }
+         else
+         {
+            aux1.GetBlock(s) = bx.GetBlock(s);
+         }
       }
       return aux1;
    }
@@ -1232,11 +1241,16 @@ void BlockNonlinearForm::Mult(const Vector &x, Vector &y) const
       {
          cP[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
       }
+      else if (needs_prolongation)
+      {
+         by.GetBlock(s) = pby.GetBlock(s);
+      }
       by.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
    }
 }
 
-void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
+void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
+                                                bool finalize) const
 {
    const int skip_zeros = 0;
    Array<Array<int> *> vdofs(fes.Size());
@@ -1490,7 +1504,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (!Grads(0,0)->Finalized())
+   if (finalize && !Grads(0,0)->Finalized())
    {
       for (int i=0; i<fes.Size(); ++i)
       {
@@ -1529,7 +1543,23 @@ Operator &BlockNonlinearForm::GetGradient(const Vector &x) const
          for (int s2 = 0; s2 < fes.Size(); ++s2)
          {
             delete cGrads(s1, s2);
-            cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
+            if (cP[s1] && cP[s2])
+            {
+               cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
+            }
+            else if (cP[s1])
+            {
+               cGrads(s1, s2) = TransposeMult(*cP[s1], *Grads(s1, s2));
+            }
+            else if (cP[s2])
+            {
+               cGrads(s1, s2) = mfem::Mult(*Grads(s1, s2), *cP[s2]);
+            }
+            else
+            {
+               cGrads(s1, s2) = NULL;
+               continue;
+            }
             mGrads(s1, s2) = cGrads(s1, s2);
          }
       }

fem/nonlinearform.hpp

@@ -217,7 +217,12 @@ public:
        In general, @a x may have non-homogeneous essential boundary values.
 
        The state @a x must be a true-dof vector. */
-   Operator &GetGradient(const Vector &x) const override;
+   Operator &GetGradient(const Vector &x) const override { return GetGradient(x, true); }
+
+   /** @brief Compute the gradient Operator of the NonlinearForm corresponding
+       to the state @a x with optional finalization and elimintaion. */
+   /** @see GetGradient(const Vector &) */
+   Operator &GetGradient(const Vector &x, bool finalize) const;
 
    /// Update the NonlinearForm to propagate updates of the associated FE space.
    /** After calling this method, the essential boundary conditions need to be
@@ -308,7 +313,7 @@ protected:
    void MultBlocked(const BlockVector &bx, BlockVector &by) const;
 
    /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx) const;
+   void ComputeGradientBlocked(const BlockVector &bx, bool finalize = true) const;
 
 public:
    /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().

fem/pbilinearform.cpp

@@ -151,6 +151,15 @@ void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
    }
 }
 
+HypreParMatrix *ParBilinearForm::ParallelAssembleInternalMatrix()
+{
+   if (p_mat.Ptr() == NULL)
+   {
+      ParallelAssemble(p_mat, mat);
+   }
+   return p_mat.As<HypreParMatrix>();
+}
+
 void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
 {
    A.Clear();
@@ -333,6 +342,15 @@ void ParBilinearForm
    A.EliminateRowsCols(dof_list, X, B);
 }
 
+void ParBilinearForm::ParallelEliminateEssentialBC(
+   const Array<int> &bdr_attr_is_ess, const HypreParVector &X, HypreParVector &B)
+{
+   Array<int> dof_list;
+   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
+
+   p_mat.As<HypreParMatrix>()->EliminateRowsCols(dof_list, X, B);
+}
+
 HypreParMatrix *ParBilinearForm::
 ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                              HypreParMatrix &A) const
@@ -344,6 +362,26 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
    return A.EliminateRowsCols(dof_list);
 }
 
+void ParBilinearForm::ParallelEliminateEssentialBC(const Array<int>
+                                                   &bdr_attr_is_ess)
+{
+   Array<int> tdofs_list;
+   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, tdofs_list);
+
+   ParallelEliminateTDofs(tdofs_list);
+}
+
+void ParBilinearForm::ParallelEliminateTDofs(const Array<int> &tdofs_list)
+{
+   p_mat_e.EliminateRowsCols(p_mat, tdofs_list);
+}
+
+void ParBilinearForm::ParallelEliminateTDofsInRHS(
+   const Array<int> &tdofs_list, const Vector &x, Vector &b)
+{
+   p_mat.EliminateBC(p_mat_e, tdofs_list, x, b);
+}
+
 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const real_t a)
 const
 {
@@ -485,7 +523,7 @@ void ParBilinearForm::FormLinearSystem(
       HypreParVector true_X(pfes), true_B(pfes);
       P.MultTranspose(b, true_B);
       R.Mult(x, true_X);
-      p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
+      ParallelEliminateTDofsInRHS(ess_tdof_list, true_X, true_B);
       R.MultTranspose(true_B, b);
       hybridization->ReduceRHS(true_B, B);
       X.SetSize(B.Size());
@@ -498,17 +536,11 @@ void ParBilinearForm::FormLinearSystem(
       B.SetSize(X.Size());
       P.MultTranspose(b, B);
       R.Mult(x, X);
-      p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
+      ParallelEliminateTDofsInRHS(ess_tdof_list, X, B);
       if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
    }
 }
 
-void ParBilinearForm::EliminateVDofsInRHS(
-   const Array<int> &vdofs, const Vector &x, Vector &b)
-{
-   p_mat.EliminateBC(p_mat_e, vdofs, x, b);
-}
-
 void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
                                        OperatorHandle &A)
 {
@@ -553,7 +585,7 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
          mat = NULL;
          delete mat_e;
          mat_e = NULL;
-         p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
+         ParallelEliminateTDofs(ess_tdof_list);
       }
       if (hybridization)
       {
@@ -615,36 +647,180 @@ void ParBilinearForm::Update(FiniteElementSpace *nfes)
    p_mat_e.Clear();
 }
 
-
-HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
+void ParMixedBilinearForm::pAllocMat()
 {
-   // construct the block-diagonal matrix A
-   HypreParMatrix *A =
-      new HypreParMatrix(trial_pfes->GetComm(),
-                         test_pfes->GlobalVSize(),
-                         trial_pfes->GlobalVSize(),
-                         test_pfes->GetDofOffsets(),
-                         trial_pfes->GetDofOffsets(),
-                         mat);
+   const int trial_nbr_size = trial_pfes->GetFaceNbrVSize();
+   const int test_nbr_size = test_pfes->GetFaceNbrVSize();
 
-   HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
-                             trial_pfes->Dof_TrueDof_Matrix());
-
-   delete A;
-
-   return rap;
+   if (keep_nbr_block)
+   {
+      mat = new SparseMatrix(height + test_nbr_size, width + trial_nbr_size);
+   }
+   else
+   {
+      mat = new SparseMatrix(height, width + trial_nbr_size);
+   }
 }
 
-void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
+void ParMixedBilinearForm::AssembleSharedFaces(int skip_zeros)
 {
-   // construct the rectangular block-diagonal matrix dA
-   OperatorHandle dA(A.Type());
-   dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
-                               test_pfes->GlobalVSize(),
-                               trial_pfes->GlobalVSize(),
-                               test_pfes->GetDofOffsets(),
-                               trial_pfes->GetDofOffsets(),
-                               mat);
+   ParMesh *pmesh = trial_pfes->GetParMesh();
+   FaceElementTransformations *T;
+   Array<int> tr_vdofs1, tr_vdofs2, tr_vdofs_all;
+   Array<int> te_vdofs1, te_vdofs2, te_vdofs_all;
+   DenseMatrix elemmat;
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+      trial_pfes->GetElementVDofs(T->Elem1No, tr_vdofs1);
+      test_pfes->GetElementVDofs(T->Elem1No, te_vdofs1);
+      trial_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, tr_vdofs2);
+      test_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, te_vdofs2);
+
+      tr_vdofs1.Copy(tr_vdofs_all);
+      for (int j = 0; j < tr_vdofs2.Size(); j++)
+      {
+         if (tr_vdofs2[j] >= 0)
+         {
+            tr_vdofs2[j] += width;
+         }
+         else
+         {
+            tr_vdofs2[j] -= width;
+         }
+      }
+      tr_vdofs_all.Append(tr_vdofs2);
+
+      if (keep_nbr_block)
+      {
+         te_vdofs1.Copy(te_vdofs_all);
+         for (int j = 0; j < te_vdofs2.Size(); j++)
+         {
+            if (te_vdofs2[j] >= 0)
+            {
+               te_vdofs2[j] += height;
+            }
+            else
+            {
+               te_vdofs2[j] -= height;
+            }
+         }
+         te_vdofs_all.Append(te_vdofs2);
+      }
+
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k]->
+         AssembleFaceMatrix(*trial_pfes->GetFE(T->Elem1No),
+                            *test_pfes->GetFE(T->Elem1No),
+                            *trial_pfes->GetFaceNbrFE(Elem2NbrNo),
+                            *test_pfes->GetFaceNbrFE(Elem2NbrNo),
+                            *T, elemmat);
+         if (keep_nbr_block)
+         {
+            mat->AddSubMatrix(te_vdofs_all, tr_vdofs_all, elemmat, skip_zeros);
+         }
+         else
+         {
+            mat->AddSubMatrix(te_vdofs1, tr_vdofs_all, elemmat, skip_zeros);
+         }
+      }
+   }
+}
+
+void ParMixedBilinearForm::Assemble(int skip_zeros)
+{
+   if (interior_face_integs.Size())
+   {
+      trial_pfes->ExchangeFaceNbrData();
+      test_pfes->ExchangeFaceNbrData();
+      if (!ext && mat == NULL)
+      {
+         pAllocMat();
+      }
+   }
+
+   MixedBilinearForm::Assemble(skip_zeros);
+
+   if (!ext && interior_face_integs.Size() > 0)
+   {
+      AssembleSharedFaces(skip_zeros);
+   }
+}
+
+HypreParMatrix *ParMixedBilinearForm::ParallelAssembleInternalMatrix()
+{
+   if (p_mat.Ptr() == NULL)
+   {
+      ParallelAssemble(p_mat, mat);
+   }
+   return p_mat.As<HypreParMatrix>();
+}
+
+HypreParMatrix *ParMixedBilinearForm::ParallelAssemble(SparseMatrix *m)
+{
+   OperatorHandle Mh(Operator::Hypre_ParCSR);
+   ParallelAssemble(Mh, m);
+   Mh.SetOperatorOwner(false);
+   return Mh.As<HypreParMatrix>();
+}
+
+void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A,
+                                            SparseMatrix *A_local)
+{
+   A.Clear();
+
+   if (A_local == NULL) { return; }
+   MFEM_VERIFY(A_local->Finalized(), "the local matrix must be finalized");
+
+   OperatorHandle dA(A.Type()), hdA;
+
+   if (interior_face_integs.Size() == 0)
+   {
+      // construct the rectangular block-diagonal matrix dA
+      dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
+                                  test_pfes->GlobalVSize(),
+                                  trial_pfes->GlobalVSize(),
+                                  test_pfes->GetDofOffsets(),
+                                  trial_pfes->GetDofOffsets(),
+                                  A_local);
+   }
+   else
+   {
+      // handle the case when 'a' contains off-diagonal
+      const int lvrows = test_pfes->GetVSize();
+      const int lvcols = trial_pfes->GetVSize();
+      const HYPRE_BigInt *face_nbr_glob_lcol = trial_pfes->GetFaceNbrGlobalDofMap();
+      const HYPRE_BigInt lcol_offset = trial_pfes->GetMyDofOffset();
+
+      Array<HYPRE_BigInt> glob_J(A_local->NumNonZeroElems());
+      const int *J = A_local->GetJ();
+      for (int i = 0; i < glob_J.Size(); i++)
+      {
+         if (J[i] < lvcols)
+         {
+            glob_J[i] = J[i] + lcol_offset;
+         }
+         else
+         {
+            glob_J[i] = face_nbr_glob_lcol[J[i] - lvcols];
+         }
+      }
+
+      // TODO - construct dA directly in the A format
+      hdA.Reset(
+         new HypreParMatrix(trial_pfes->GetComm(), lvrows, test_pfes->GlobalVSize(),
+                            trial_pfes->GlobalVSize(), A_local->GetI(), glob_J,
+                            A_local->GetData(), test_pfes->GetDofOffsets(),
+                            trial_pfes->GetDofOffsets()));
+      // - hdA owns the new HypreParMatrix
+      // - the above constructor copies all input arrays
+      glob_J.DeleteAll();
+      dA.ConvertFrom(hdA);
+   }
 
    OperatorHandle P_test(A.Type()), P_trial(A.Type());
 
@@ -670,6 +846,44 @@ void ParMixedBilinearForm::TrueAddMult(const Vector &x, Vector &y,
    test_pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Yaux, 1.0, y);
 }
 
+void ParMixedBilinearForm::ParallelEliminateTrialEssentialBC(
+   const Array<int> &bdr_attr_is_ess)
+{
+   Array<int> trial_tdof_list;
+   trial_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, trial_tdof_list);
+
+   ParallelEliminateTrialTDofs(trial_tdof_list);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTrialTDofs(
+   const Array<int> &trial_tdof_list)
+{
+   HypreParMatrix *temp = p_mat.As<HypreParMatrix>()->EliminateCols(
+                             trial_tdof_list);
+   p_mat_e.Reset(temp, true);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTrialTDofsInRHS(
+   const Array<int> &trial_tdof_list, const Vector &x, Vector &b)
+{
+   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, x, 1.0, b);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTestEssentialBC(
+   const Array<int> &bdr_attr_is_ess)
+{
+   Array<int> test_tdof_list;
+   test_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, test_tdof_list);
+
+   ParallelEliminateTestTDofs(test_tdof_list);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTestTDofs(
+   const Array<int> &test_tdof_list)
+{
+   p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
+}
+
 void ParMixedBilinearForm::FormRectangularSystemMatrix(
    const Array<int>
    &trial_tdof_list,
@@ -690,10 +904,8 @@ void ParMixedBilinearForm::FormRectangularSystemMatrix(
       mat = NULL;
       delete mat_e;
       mat_e = NULL;
-      HypreParMatrix *temp =
-         p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
-      p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
-      p_mat_e.Reset(temp, true);
+      ParallelEliminateTrialTDofs(trial_tdof_list);
+      ParallelEliminateTestTDofs(test_tdof_list);
    }
 
    A = p_mat;
@@ -723,7 +935,7 @@ void ParMixedBilinearForm::FormRectangularLinearSystem(
    test_P->MultTranspose(b, B);
    trial_R->Mult(x, X);
 
-   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
+   ParallelEliminateTrialTDofsInRHS(trial_tdof_list, X, B);
    B.SetSubVector(test_tdof_list, 0.0);
 }
 

fem/pbilinearform.hpp

@@ -73,7 +73,7 @@ public:
    /** When set to true and the ParBilinearForm has interior face integrators,
        the local SparseMatrix will include the rows (in addition to the columns)
        corresponding to face-neighbor dofs. The default behavior is to disregard
-       those rows. Must be called before the first Assemble call. */
+       those rows. Must be called before the first Assemble() call. */
    void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
 
    /** @brief Set the operator type id for the parallel matrix/operator when
@@ -101,6 +101,14 @@ public:
        diagonal for this case. */
    void AssembleDiagonal(Vector &diag) const override;
 
+   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   /** The returned matrix is the internal one, owned by the form. It is not
+       reassembled if it has been already constructed. If FormSystemMatrix()
+       has been called before, it is the system matrix with eliminated
+       essential DOFs, otherwise the parallel matrix is assembled here without
+       the elimination process. */
+   HypreParMatrix *ParallelAssembleInternalMatrix();
+
    /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
    /** The returned matrix has to be deleted by the caller. */
    HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
@@ -146,6 +154,13 @@ public:
                                      const HypreParVector &X,
                                      HypreParVector &B) const;
 
+   /// Eliminate essential boundary DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. */
+   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
+                                     const HypreParVector &X,
+                                     HypreParVector &B);
+
    /// Eliminate essential boundary DOFs from a parallel assembled matrix @a A.
    /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
        the essential part of the boundary. The eliminated part is stored in a
@@ -157,6 +172,12 @@ public:
    HypreParMatrix *ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                                                 HypreParMatrix &A) const;
 
+   /// Eliminate essential boundary DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. This method relies on
+       ParallelEliminateTDofs(const Array<int> &), see it for details. */
+   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess);
+
    /// Eliminate essential true DOFs from a parallel assembled matrix @a A.
    /** Given a list of essential true dofs and the parallel assembled matrix
        @a A, eliminate the true dofs from the matrix, storing the eliminated
@@ -169,6 +190,28 @@ public:
                                           HypreParMatrix &A) const
    { return A.EliminateRowsCols(tdofs_list); }
 
+   /// Eliminate essential true DOFs from the parallel system matrix.
+   /** Given a list of essential true dofs, eliminate the true dofs from
+       the parallel assembled system matrix, storing the eliminated part
+       internally. This method works in conjunction with
+       ParallelEliminateTDofsInRHS() and allows elimination of boundary
+       conditions in multiple right-hand sides. */
+   void ParallelEliminateTDofs(const Array<int> &tdofs_list);
+
+   /** @brief Use the stored eliminated part of the parallel system matrix for
+       elimination of boundary conditions in the r.h.s. */
+   /** Given a list of essential true dofs, eliminate the true dofs from the
+       right-hand side @a b using the solution vector @a x and the previously
+       stored eliminated part of the parallel assembled system matrix produced
+       by ParallelEliminateTDofs(const Array<int> &). */
+   void ParallelEliminateTDofsInRHS(const Array<int> &tdofs, const Vector &x,
+                                    Vector &b);
+
+   /// @deprecated Use ParallelEliminateTDofsInRHS() instead.
+   MFEM_DEPRECATED void EliminateVDofsInRHS(const Array<int> &vdofs,
+                                            const Vector &x, Vector &b)
+   { ParallelEliminateTDofsInRHS(vdofs, x, b); }
+
    /** @brief Compute @a y += @a a (P^t A P) @a x, where @a x and @a y are
        vectors on the true dofs. */
    void TrueAddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
@@ -238,8 +281,6 @@ public:
 
    void Update(FiniteElementSpace *nfes = NULL) override;
 
-   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
-
    virtual ~ParBilinearForm() { }
 };
 
@@ -257,6 +298,13 @@ protected:
    /// Matrix and eliminated matrix
    OperatorHandle p_mat, p_mat_e;
 
+   bool keep_nbr_block;
+
+   // Allocate mat - called when (mat == NULL && fbfi.Size() > 0)
+   void pAllocMat();
+
+   void AssembleSharedFaces(int skip_zeros = 1);
+
 private:
    /// Copy construction is not supported; body is undefined.
    ParMixedBilinearForm(const ParMixedBilinearForm &);
@@ -276,6 +324,7 @@ public:
    {
       trial_pfes = trial_fes;
       test_pfes  = test_fes;
+      keep_nbr_block = false;
    }
 
    /** @brief Create a ParMixedBilinearForm on the given FiniteElementSpace%s
@@ -295,15 +344,89 @@ public:
    {
       trial_pfes = trial_fes;
       test_pfes  = test_fes;
+      keep_nbr_block = false;
    }
 
+   /** When set to true and the ParMixedBilinearForm has interior face
+       integrators, the local SparseMatrix will include the rows (in addition
+       to the columns) corresponding to face-neighbor dofs. The default
+       behavior is to disregard those rows. Must be called before the first
+       Assemble() call. */
+   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
    /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
-   HypreParMatrix *ParallelAssemble();
+   /** The returned matrix is the internal one, owned by the form. It is not
+       reassembled if it has been already constructed. If
+       FormRectangularSystemMatrix() has been called before, it is the system
+       matrix with eliminated essential DOFs, otherwise the parallel matrix is
+       assembled here without the elimination process. */
+   HypreParMatrix *ParallelAssembleInternalMatrix();
+
+   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
+
+   /** @brief Returns the eliminated matrix assembled on the true dofs, i.e.
+       P_test^t A_local P_trial. */
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssembleElim() { return ParallelAssemble(mat_e); }
+
+   /** @brief Return the matrix @a m assembled on the true dofs, i.e. P_test^t
+       A_local P_trial. */
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssemble(SparseMatrix *m);
 
    /** @brief Returns the matrix assembled on the true dofs, i.e.
        @a A = P_test^t A_local P_trial, in the format (type id) specified by
        @a A. */
-   void ParallelAssemble(OperatorHandle &A);
+   void ParallelAssemble(OperatorHandle &A) { ParallelAssemble(A, mat); }
+
+   /** Returns the eliminated matrix assembled on the true dofs, i.e.
+       @a A_elim = P^t A_elim_local P in the format (type id) specified by @a A.
+    */
+   void ParallelAssembleElim(OperatorHandle &A_elim)
+   { ParallelAssemble(A_elim, mat_e); }
+
+   /** Returns the matrix @a A_local assembled on the true dofs, i.e.
+       @a A = P_test^t A_local P_trial in the format (type id) specified by
+       @a A. */
+   void ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local);
+
+   /// Eliminate essential boundary trial DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. This method relies on
+       ParallelEliminateTrialTDofs(const Array<int> &), see it for details. */
+   void ParallelEliminateTrialEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// Eliminate essential trial true DOFs from the parallel system matrix.
+   /** Given a list of essential trial true dofs, eliminate the trial true dofs
+       from the parallel assembled system matrix, storing the eliminated part
+       internally. This method works in conjunction with
+       ParallelEliminateTrialTDofsInRHS() and allows elimination of boundary
+       conditions in multiple right-hand sides. */
+   void ParallelEliminateTrialTDofs(const Array<int> &trial_tdof_list);
+
+   /** @brief Use the stored eliminated part of the parallel system matrix for
+       elimination of boundary conditions in the r.h.s. */
+   /** Given a list of essential trial true dofs, eliminate the trial true dofs
+       from the right-hand side @a B using the solution vector @a X and the
+       previously stored eliminated part of the parallel assembled system
+       matrix produced by ParallelEliminateTrialTDofs(const Array<int> &). */
+   void ParallelEliminateTrialTDofsInRHS(const Array<int> &trial_tdof_list,
+                                         const Vector &X, Vector &B);
+
+   /// Eliminate essential boundary test DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. */
+   void ParallelEliminateTestEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// Eliminate essential test true DOFs from the parallel system matrix.
+   /** Given a list of essential test true dofs, eliminate the test true dofs
+       from the parallel assembled system matrix. */
+   void ParallelEliminateTestTDofs(const Array<int> &test_tdof_list);
 
    using MixedBilinearForm::FormRectangularSystemMatrix;
    using MixedBilinearForm::FormRectangularLinearSystem;

fem/pnonlinearform.cpp

@@ -105,6 +105,59 @@ const SparseMatrix &ParNonlinearForm::GetLocalGradient(const Vector &x) const
    return *Grad;
 }
 
+void ParNonlinearForm::GradientSharedFaces(const Vector &x,
+                                           int skip_zeros) const
+{
+   ParFiniteElementSpace *pfes = ParFESpace();
+   ParMesh *pmesh = pfes->GetParMesh();
+   FaceElementTransformations *T;
+   Array<int> vdofs1, vdofs2, vdofs_all;
+   DenseMatrix elemmat;
+   Vector el_x, nbr_x, face_x;
+   const Vector &px = Prolongate(x);
+
+   ParGridFunction pgf(pfes, const_cast<Vector&>(px), 0);
+   pgf.ExchangeFaceNbrData();
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+
+      pfes->GetElementVDofs(T->Elem1No, vdofs1);
+      pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
+      face_x.SetSize(vdofs1.Size() + vdofs2.Size());
+
+      el_x.MakeRef(face_x, 0, vdofs1.Size());
+      pgf.GetSubVector(vdofs1, el_x);
+
+      nbr_x.MakeRef(face_x, vdofs1.Size(), vdofs2.Size());
+      pgf.FaceNbrData().GetSubVector(vdofs2, nbr_x);
+
+      vdofs1.Copy(vdofs_all);
+      for (int j = 0; j < vdofs2.Size(); j++)
+      {
+         if (vdofs2[j] >= 0)
+         {
+            vdofs2[j] += height;
+         }
+         else
+         {
+            vdofs2[j] -= height;
+         }
+      }
+      vdofs_all.Append(vdofs2);
+      for (int k = 0; k < fnfi.Size(); k++)
+      {
+         fnfi[k]->AssembleFaceGrad(*pfes->GetFE(T->Elem1No),
+                                   *pfes->GetFaceNbrFE(Elem2NbrNo),
+                                   *T, face_x, elemmat);
+         Grad->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
+      }
+   }
+}
+
 Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 {
    if (NonlinearForm::ext) { return NonlinearForm::GetGradient(x); }
@@ -112,19 +165,61 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
    ParFiniteElementSpace *pfes = ParFESpace();
 
    pGrad.Clear();
+   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type()), hdA;
 
-   NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
-
-   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
-
-   if (fnfi.Size() == 0)
+   if (fnfi.Size())
    {
-      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
-                             pfes->GetDofOffsets(), Grad);
+      const int skip_zeros = 0;
+
+      pfes->ExchangeFaceNbrData();
+      if (Grad == NULL)
+      {
+         int nbr_size = pfes->GetFaceNbrVSize();
+         Grad = new SparseMatrix(pfes->GetVSize(), pfes->GetVSize() + nbr_size);
+      }
+
+      NonlinearForm::GetGradient(x, false); // (re)assemble Grad, no b.c.
+
+      GradientSharedFaces(x, skip_zeros);
+
+      Grad->Finalize(skip_zeros);
+
+      // handle the case when 'a' contains off-diagonal
+      int lvsize = pfes->GetVSize();
+      const HYPRE_BigInt *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
+      HYPRE_BigInt ldof_offset = pfes->GetMyDofOffset();
+
+      Array<HYPRE_BigInt> glob_J(Grad->NumNonZeroElems());
+      int *J = Grad->GetJ();
+      for (int i = 0; i < glob_J.Size(); i++)
+      {
+         if (J[i] < lvsize)
+         {
+            glob_J[i] = J[i] + ldof_offset;
+         }
+         else
+         {
+            glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
+         }
+      }
+
+      // TODO - construct dA directly in the A format
+      hdA.Reset(
+         new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
+                            pfes->GlobalVSize(), Grad->GetI(), glob_J,
+                            Grad->GetData(), pfes->GetDofOffsets(),
+                            pfes->GetDofOffsets()));
+      // - hdA owns the new HypreParMatrix
+      // - the above constructor copies all input arrays
+      glob_J.DeleteAll();
+      dA.ConvertFrom(hdA);
    }
    else
    {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+      NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
+
+      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
+                             pfes->GetDofOffsets(), Grad);
    }
 
    // RAP the local gradient dA.
@@ -271,7 +366,70 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
 
    if (fnfi.Size() > 0)
    {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+      // Terms over shared interior faces in parallel.
+      ParMesh *pmesh = ParFESpace(0)->GetParMesh();
+      FaceElementTransformations *tr;
+
+      Array<Array<int> *>vdofs(fes.Size());
+      Array<Array<int> *>vdofs2(fes.Size());
+      Array<Vector *> el_x(fes.Size());
+      Array<const Vector *> el_x_const(fes.Size());
+      Array<Vector *> el_y(fes.Size());
+      Array<const FiniteElement *> fe(fes.Size());
+      Array<const FiniteElement *> fe2(fes.Size());
+      Array<ParGridFunction *> pgfs(fes.Size());
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         el_x_const[s] = el_x[s] = new Vector();
+         el_y[s] = new Vector();
+         vdofs[s] = new Array<int>;
+         vdofs2[s] = new Array<int>;
+         pgfs[s] = new ParGridFunction(const_cast<ParFiniteElementSpace*>(ParFESpace(s)),
+                                       xs.GetBlock(s));
+         pgfs[s]->ExchangeFaceNbrData();
+      }
+
+      const int n_shared_faces = pmesh->GetNSharedFaces();
+      for (int i = 0; i < n_shared_faces; i++)
+      {
+         tr = pmesh->GetSharedFaceTransformations(i, true);
+         int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
+
+         for (int s=0; s<fes.Size(); ++s)
+         {
+            const ParFiniteElementSpace *pfes = ParFESpace(s);
+            fe[s] = pfes->GetFE(tr->Elem1No);
+            fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
+
+            pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
+            pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
+
+            el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
+            xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
+            pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
+                                                el_x[s]->GetData() + vdofs[s]->Size());
+         }
+
+         for (int k = 0; k < fnfi.Size(); ++k)
+         {
+            fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
+
+            for (int s=0; s<fes.Size(); ++s)
+            {
+               if (el_y[s]->Size() == 0) { continue; }
+               ys.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
+            }
+         }
+      }
+
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         delete pgfs[s];
+         delete vdofs2[s];
+         delete vdofs[s];
+         delete el_y[s];
+         delete el_x[s];
+      }
    }
 
    for (int s=0; s<fes.Size(); ++s)
@@ -328,6 +486,106 @@ void ParBlockNonlinearForm::SetGradientType(Operator::Type tid)
    }
 }
 
+void ParBlockNonlinearForm::GradientSharedFaces(const BlockVector &xs,
+                                                int skip_zeros) const
+{
+   // Terms over shared interior faces in parallel.
+   ParMesh *pmesh = ParFESpace(0)->GetParMesh();
+   FaceElementTransformations *tr;
+
+   Array<Array<int> *>vdofs(fes.Size());
+   Array<Array<int> *>vdofs2(fes.Size());
+   Array<Array<int> *>vdofs_all(fes.Size());
+   Array<Vector *> el_x(fes.Size());
+   Array<const Vector *> el_x_const(fes.Size());
+   Array2D<DenseMatrix *> elmats(fes.Size(), fes.Size());
+   Array<const FiniteElement *> fe(fes.Size());
+   Array<const FiniteElement *> fe2(fes.Size());
+   Array<ParGridFunction *> pgfs(fes.Size());
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      el_x_const[s1] = el_x[s1] = new Vector();
+      vdofs[s1] = new Array<int>;
+      vdofs2[s1] = new Array<int>;
+      vdofs_all[s1] = new Array<int>;
+      pgfs[s1] = new ParGridFunction(
+         const_cast<ParFiniteElementSpace*>(ParFESpace(s1)),
+         const_cast<Vector&>(xs.GetBlock(s1)));
+      pgfs[s1]->ExchangeFaceNbrData();
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         elmats(s1,s2) = new DenseMatrix();
+      }
+   }
+
+   const int n_shared_faces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < n_shared_faces; i++)
+   {
+      tr = pmesh->GetSharedFaceTransformations(i, true);
+      int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
+
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         const ParFiniteElementSpace *pfes = ParFESpace(s);
+         fe[s] = pfes->GetFE(tr->Elem1No);
+         fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
+
+         pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
+         pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
+
+         el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
+         xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
+         pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
+                                             el_x[s]->GetData() + vdofs[s]->Size());
+
+         vdofs[s]->Copy(*vdofs_all[s]);
+
+         const int lvsize = pfes->GetVSize();
+         for (int j = 0; j < vdofs2[s]->Size(); j++)
+         {
+            if ((*vdofs2[s])[j] >= 0)
+            {
+               (*vdofs2[s])[j] += lvsize;
+            }
+            else
+            {
+               (*vdofs2[s])[j] -= lvsize;
+            }
+         }
+         vdofs_all[s]->Append(*(vdofs2[s]));
+      }
+
+      for (int k = 0; k < fnfi.Size(); ++k)
+      {
+         fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
+
+         for (int s1=0; s1<fes.Size(); ++s1)
+         {
+            for (int s2=0; s2<fes.Size(); ++s2)
+            {
+               if (elmats(s1,s2)->Height() == 0) { continue; }
+               Grads(s1,s2)->AddSubMatrix(*vdofs[s1], *vdofs_all[s2],
+                                          *elmats(s1,s2), skip_zeros);
+            }
+         }
+      }
+   }
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      delete pgfs[s1];
+      delete vdofs_all[s1];
+      delete vdofs2[s1];
+      delete vdofs[s1];
+      delete el_x[s1];
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         delete elmats(s1,s2);
+      }
+   }
+}
+
 BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
 {
    if (pBlockGrad == NULL)
@@ -347,49 +605,155 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
       }
    }
 
-   GetLocalGradient(x); // gradients are stored in 'Grads'
+   // xs_true is not modified, so const_cast is okay
+   xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
+   xs.Update(block_offsets);
+
+   for (int s=0; s<fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->Mult(
+         xs_true.GetBlock(s), xs.GetBlock(s));
+   }
 
    if (fnfi.Size() > 0)
    {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
+      const int skip_zeros = 0;
 
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
+      for (int s=0; s<fes.Size(); ++s)
       {
-         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                        Ph(phBlockGrad(s1,s2)->Type()),
-                        Rh(phBlockGrad(s1,s2)->Type());
+         const_cast<ParFiniteElementSpace*>(pfes[s])->ExchangeFaceNbrData();
+      }
 
-         if (s1 == s2)
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
          {
-            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-            OperatorHandle Ae;
-            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            if (Grads(s1,s2) == NULL)
+            {
+               int nbr_size = pfes[s2]->GetFaceNbrVSize();
+               Grads(s1,s2) = new SparseMatrix(pfes[s1]->GetVSize(),
+                                               pfes[s2]->GetVSize() + nbr_size);
+            }
          }
-         else
+      }
+
+      // (re)assemble Grad without b.c. into 'Grads'
+      BlockNonlinearForm::ComputeGradientBlocked(xs, false);
+
+      GradientSharedFaces(xs, skip_zeros);
+
+      // finalize the gradients
+      for (int s1=0; s1<fes.Size(); ++s1)
+         for (int s2=0; s2<fes.Size(); ++s2)
          {
-            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                        pfes[s1]->GlobalVSize(),
-                                        pfes[s2]->GlobalVSize(),
-                                        pfes[s1]->GetDofOffsets(),
-                                        pfes[s2]->GetDofOffsets(),
-                                        Grads(s1,s2));
-            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            Grads(s1,s2)->Finalize(skip_zeros);
          }
 
-         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
+         {
+            OperatorHandle hdA;
+            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                           Ph(phBlockGrad(s1,s2)->Type()),
+                           Rh(phBlockGrad(s1,s2)->Type());
+
+            // handle the case when 'a' contains off-diagonal
+            int lvsize = pfes[s2]->GetVSize();
+            const HYPRE_BigInt *face_nbr_glob_ldof =
+               const_cast<ParFiniteElementSpace*>(pfes[s2])->GetFaceNbrGlobalDofMap();
+            HYPRE_BigInt ldof_offset = pfes[s2]->GetMyDofOffset();
+
+            Array<HYPRE_BigInt> glob_J(Grads(s1,s2)->NumNonZeroElems());
+            int *J = Grads(s1,s2)->GetJ();
+            for (int i = 0; i < glob_J.Size(); i++)
+            {
+               if (J[i] < lvsize)
+               {
+                  glob_J[i] = J[i] + ldof_offset;
+               }
+               else
+               {
+                  glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
+               }
+            }
+
+            // TODO - construct dA directly in the A format
+            hdA.Reset(
+               new HypreParMatrix(pfes[s2]->GetComm(), pfes[s1]->GetVSize(),
+                                  pfes[s1]->GlobalVSize(), pfes[s2]->GlobalVSize(),
+                                  Grads(s1,s2)->GetI(), glob_J, Grads(s1,s2)->GetData(),
+                                  pfes[s1]->GetDofOffsets(), pfes[s2]->GetDofOffsets()));
+            // - hdA owns the new HypreParMatrix
+            // - the above constructor copies all input arrays
+            glob_J.DeleteAll();
+            dA.ConvertFrom(hdA);
+
+            if (s1 == s2)
+            {
+               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
+
+               OperatorHandle Ae;
+               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            }
+            else
+            {
+               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            }
+
+            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+         }
+      }
+   }
+   else
+   {
+      // (re)assemble Grad without b.c. into 'Grads'
+      BlockNonlinearForm::ComputeGradientBlocked(xs);
+
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
+         {
+            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                           Ph(phBlockGrad(s1,s2)->Type()),
+                           Rh(phBlockGrad(s1,s2)->Type());
+
+            if (s1 == s2)
+            {
+               dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
+                                      pfes[s1]->GetDofOffsets(), Grads(s1,s1));
+               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
+
+               OperatorHandle Ae;
+               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            }
+            else
+            {
+               dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
+                                           pfes[s1]->GlobalVSize(),
+                                           pfes[s2]->GlobalVSize(),
+                                           pfes[s1]->GetDofOffsets(),
+                                           pfes[s2]->GetDofOffsets(),
+                                           Grads(s1,s2));
+               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            }
+
+            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+         }
       }
    }
 

fem/pnonlinearform.hpp

@@ -29,6 +29,8 @@ protected:
    mutable ParGridFunction X, Y;
    mutable OperatorHandle pGrad;
 
+   void GradientSharedFaces(const Vector &x, int skip_zeros = 1) const;
+
 public:
    ParNonlinearForm(ParFiniteElementSpace *pf);
 
@@ -81,6 +83,8 @@ protected:
    mutable Array2D<OperatorHandle *> phBlockGrad;
    mutable BlockOperator *pBlockGrad;
 
+   void GradientSharedFaces(const BlockVector &xs, int skip_zeros) const;
+
 public:
    /// Computes the energy of the system
    real_t GetEnergy(const Vector &x) const override;

linalg/hypre.cpp

@@ -561,7 +561,8 @@ void CopyMemory(Memory<T> &src, Memory<T> &dst, MemoryClass dst_mc,
     this function. In particular, @a dst should be empty or deleted before
     calling this function. */
 template <typename SrcT, typename DstT>
-void CopyConvertMemory(Memory<SrcT> &src, MemoryClass dst_mc, Memory<DstT> &dst)
+void CopyConvertMemory(const Memory<SrcT> &src, MemoryClass dst_mc,
+                       Memory<DstT> &dst)
 {
    auto capacity = src.Capacity();
    dst.New(capacity, GetMemoryType(dst_mc));
@@ -842,8 +843,8 @@ static int GetPartitioningArraySize(MPI_Comm comm)
 ///
 /// Both @a row and @a col are partitioning arrays, whose length is returned by
 /// GetPartitioningArraySize(), see @ref hypre_partitioning_descr.
-static bool RowAndColStartsAreEqual(MPI_Comm comm, HYPRE_BigInt *rows,
-                                    HYPRE_BigInt *cols)
+static bool RowAndColStartsAreEqual(MPI_Comm comm, const HYPRE_BigInt *rows,
+                                    const HYPRE_BigInt *cols)
 {
    const int part_size = GetPartitioningArraySize(comm);
    bool are_equal = true;
@@ -1131,7 +1132,7 @@ HypreParMatrix::HypreParMatrix(
 HypreParMatrix::HypreParMatrix(MPI_Comm comm,
                                HYPRE_BigInt *row_starts,
                                HYPRE_BigInt *col_starts,
-                               SparseMatrix *sm_a)
+                               const SparseMatrix *sm_a)
 {
    MFEM_ASSERT(sm_a != NULL, "invalid input");
    MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
@@ -1145,7 +1146,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,
 
    hypre_CSRMatrixSetDataOwner(csr_a,0);
    MemoryIJData mem_a;
-   CopyCSR(sm_a, mem_a, csr_a, false);
+   CopyCSR(const_cast<SparseMatrix*>(sm_a), mem_a, csr_a, false);
    hypre_CSRMatrixSetRownnz(csr_a);
 
    // NOTE: this call creates a matrix on host even when device support is
@@ -1307,10 +1308,11 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows,
                                HYPRE_BigInt glob_nrows,
                                HYPRE_BigInt glob_ncols,
-                               int *I, HYPRE_BigInt *J,
-                               real_t *data,
-                               HYPRE_BigInt *rows,
-                               HYPRE_BigInt *cols)
+                               const int *I,
+                               const HYPRE_BigInt *J,
+                               const real_t *data,
+                               const HYPRE_BigInt *rows,
+                               const HYPRE_BigInt *cols)
 {
    Init();
 

linalg/hypre.hpp

@@ -565,7 +565,7 @@ public:
        partitioning arrays @a row_starts and @a col_starts. */
    HypreParMatrix(MPI_Comm comm, HYPRE_BigInt *row_starts,
                   HYPRE_BigInt *col_starts,
-                  SparseMatrix *a); // constructor with 4 arguments, v2
+                  const SparseMatrix *a); // constructor with 4 arguments, v2
 
    /// Creates boolean block-diagonal rectangular parallel matrix.
    /** The new HypreParMatrix does not take ownership of any of the input
@@ -594,9 +594,9 @@ public:
        arrays (so they can be deleted). See @ref hypre_partitioning_descr "here"
        for a description of the partitioning arrays @a rows and @a cols. */
    HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_BigInt glob_nrows,
-                  HYPRE_BigInt glob_ncols, int *I, HYPRE_BigInt *J,
-                  real_t *data, HYPRE_BigInt *rows,
-                  HYPRE_BigInt *cols); // constructor with 9 arguments
+                  HYPRE_BigInt glob_ncols, const int *I, const HYPRE_BigInt *J,
+                  const real_t *data, const HYPRE_BigInt *rows,
+                  const HYPRE_BigInt *cols); // constructor with 9 arguments
 
    /** @brief Copy constructor for a ParCSR matrix which creates a deep copy of
        structure and data from @a P. */

miniapps/navier/navier_cht.cpp

@@ -461,7 +461,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 
    Kmat.Mult(u, z);
    z.Neg(); // z = -z
-   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
+   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
 
    M_solver.Mult(z, du_dt);
    du_dt.Print();
@@ -483,7 +483,7 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
    MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
    Kmat.Mult(u, z);
    z.Neg();
-   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
+   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
 
    T_solver.Mult(z, du_dt);
    du_dt.SetSubVector(ess_tdof_list, 0.0);

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);
+}