updting io benchmark

Boundary integral order bug [bdr-int-order-fix]

.gitignore

@@ -231,6 +231,14 @@ miniapps/gslib/field-diff
 miniapps/gslib/findpts
 miniapps/gslib/pfindpts
 
+miniapps/navier/navier_mms
+miniapps/navier/navier_kovasznay
+miniapps/navier/navier_tgv
+miniapps/navier/navier_shear
+miniapps/navier/navier_3dfoc
+miniapps/navier/tgv_out*.txt
+miniapps/navier/*_output
+
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests

CHANGELOG

@@ -52,6 +52,12 @@ New and updated examples and miniapps
 - Added a new Example 26/26p to demonstrate the construction of a matrix-free
   geometric and p-multigrid preconditioner for the Laplace problem.
 
+- Added a new example, Example 27/27p, to demonstrate the enforcement of
+  various boundary conditions with the Laplace operator.  The example shows the
+  procedures for applying Dirichlet, Neumann (both homogeneous and
+  inhomogeneous), Robin, and periodic boundary conditions with either H1 or DG
+  discretizations.
+
 - Added a simple meshing miniapp, Twist, which demonstrates MFEM's strategy of
   stitching together opposite surfaces of a mesh to create a topologically
   periodic mesh.

doc/CodeDocumentation.conf.in

@@ -774,6 +774,7 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                          @MFEM_SOURCE_DIR@/miniapps/electromagnetics \
                          @MFEM_SOURCE_DIR@/miniapps/gslib \
                          @MFEM_SOURCE_DIR@/miniapps/meshing \
+                         @MFEM_SOURCE_DIR@/miniapps/navier \
                          @MFEM_SOURCE_DIR@/miniapps/nurbs \
                          @MFEM_SOURCE_DIR@/miniapps/performance \
                          @MFEM_SOURCE_DIR@/miniapps/tools \

doc/CodeDocumentation.dox

@@ -156,6 +156,7 @@ namespace mfem {
  * - <a class="el" href="lor-transfer_8cpp_source.html">LOR Transfer</a>: map functions between high-order and low-order refined spaces
  * - <a class="el" href="findpts_8cpp_source.html">Find Points</a>: evaluate grid function in physical space, <a class="el" href="findpts_8cpp_source.html">serial</a> and <a class="el" href="pfindpts_8cpp_source.html">parallel</a> versions
  * - <a class="el" href="field-diff_8cpp_source.html">Field Diff</a>: compare grid functions on different meshes
+ * - <a class="el" href="classmfem_1_1navier_1_1NavierSolver.html">Navier</a>: solve the transient incompressible Navier-Stokes equations
  * - <a class="el" href="miniapps_2performance_2ex1_8cpp_source.html">HPC Example 1</a>: high-performance nodal H1 FEM for the Laplace problem
  * - <a class="el" href="miniapps_2performance_2ex1p_8cpp_source.html">HPC Example 1p</a>: high-performance parallel nodal H1 FEM for the Laplace problem
  *

examples/CMakeLists.txt

@@ -33,6 +33,7 @@ list(APPEND ALL_EXE_SRCS
   ex24.cpp
   ex25.cpp
   ex26.cpp
+  ex27.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -62,6 +63,9 @@ if (MFEM_USE_MPI)
     ex24p.cpp
     ex25p.cpp
     ex26p.cpp
+    ex27p.cpp
+    pa_oper.cpp
+    io_benchmark.cpp
     )
 endif()
 
@@ -81,6 +85,8 @@ foreach(SRC_FILE ${ALL_EXE_SRCS})
     list(APPEND THIS_TEST_OPTIONS "-tf" "5")
   elseif(${TEST_NAME} MATCHES "ex15p*")
     list(APPEND THIS_TEST_OPTIONS "-e" "1")
+  elseif(${TEST_NAME} MATCHES "ex27p*")
+    list(APPEND THIS_TEST_OPTIONS "-dg")
   endif()
 
   if (NOT (${TEST_NAME} MATCHES ".*p$"))

examples/ex1p.cpp

@@ -70,6 +70,7 @@ int main(int argc, char *argv[])
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   int nfiles = 1;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -86,6 +87,7 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&nfiles, "-nf", "--num-files", "Number of files to write.");
    args.Parse();
    if (!args.Good())
    {
@@ -158,7 +160,7 @@ int main(int argc, char *argv[])
    {
       fec = new H1_FECollection(order = 1, dim);
    }
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec, 1, 0);
    HYPRE_Int size = fespace->GlobalTrueVSize();
    if (myid == 0)
    {
@@ -237,20 +239,52 @@ int main(int argc, char *argv[])
    //     local finite element solution on each processor.
    a->RecoverFEMSolution(X, *b, x);
 
+   std::string filename("nranks_");
+   filename += to_string(num_procs);
+   filename += ".gf";
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      x.Save(filename.c_str(), nfiles);
+      double t2 = MPI_Wtime();
+
+      if (myid == 0)
+      {
+         err << "elapsed write time: " << t2 - t1 << endl;
+      }
+   }
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      ParGridFunction new_x(fespace, filename.c_str());
+      double t2 = MPI_Wtime();
+      if (myid == 0)
+      {
+         err << "elapsed read time: " << t2 - t1 << endl;
+      }
+      // new_x -= x;
+      // out << "GF difference: " << new_x.Norml1() << endl;
+   }
+
    // 15. Save the refined mesh and the solution in parallel. This output can
    //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
    {
       ostringstream mesh_name, sol_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "sol." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
+      //mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << num_procs << setfill('0') << setw(6) << myid;
 
+      //ofstream mesh_ofs(mesh_name.str().c_str());
+      //mesh_ofs.precision(8);
+      //pmesh->Print(mesh_ofs);
+      double t1 = MPI_Wtime();
       ofstream sol_ofs(sol_name.str().c_str());
       sol_ofs.precision(8);
       x.Save(sol_ofs);
+      double t2 = MPI_Wtime();
+      if (myid == 0)
+      {
+         err << t2 - t1 << endl;
+      }
    }
 
    // 16. Send the solution by socket to a GLVis server.

examples/ex27.cpp

@@ -0,0 +1,736 @@
+//                       MFEM Example 27 - Serial Version
+//
+// Compile with: make ex27
+//
+// Sample runs:  ex27
+//               ex27 -dg
+//               ex27 -dg -dbc 8 -nbc -2
+//               ex27 -rbc-a 1 -rbc-b 8
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 0 with a variety of boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order using a continuous or discontinuous space.  We then
+//               apply Dirichlet, Neumann (both homogeneous and inhomogeneous),
+//               Robin, and Periodic boundary conditions on different portions
+//               of a predefined mesh.
+//
+//               The predefined mesh consists of a rectangle with two
+//               holes removed (see below).  The narrow ends of the
+//               mesh are connected to form a Periodic boundary
+//               condition.  The lower edge (tagged with attribute 1)
+//               receives an inhomogeneous Neumann boundary condition.
+//               A Robin boundary condition is applied to upper edge
+//               (attribute 2).  The circular hole on the left
+//               (attribute 3) enforces a Dirichlet boundary
+//               condition.  Finally, a natural boundary condition, or
+//               homogeneous Neumann BC, is applied to the circular
+//               hole on the right (attribute 4).
+//
+//                  Attribute 3    ^ y  Attribute 2
+//                        \        |      /
+//                     +-----------+-----------+
+//                     |    \_     |     _     |
+//                     |    / \    |    / \    |
+//                  <--+---+---+---+---+---+---+--> x
+//                     |    \_/    |    \_/    |
+//                     |           |      \    |
+//                     +-----------+-----------+       (hole radii are
+//                          /      |        \            adjustable)
+//                  Attribute 1    v    Attribute 4
+//
+//               The boundary conditions are defined as (where u is
+//               the solution field):
+//                  Dirichlet: u = d
+//                  Neumann:   n.Grad(u) = g
+//                  Robin:     n.Grad(u) + a u = b
+//
+//               The user can adjust the values of 'd', 'g', 'a', and
+//               'b' with command line options.
+//
+//               This example highlights the differing implementations of
+//               boundary conditions with continuous and discontinuous Galerkin
+//               formulations of the Laplace problem.
+//
+//               We recommend viewing examples 1 and 14 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double a_ = 0.2;
+
+// Normal to hole with boundary attribute 4
+void n4Vec(const Vector &x, Vector &n) { n = x; n[0] -= 0.5; n /= -n.Norml2(); }
+
+Mesh * GenerateSerialMesh(int ref);
+
+// Compute the average value of alpha*n.Grad(sol) + beta*sol over the boundary
+// attributes marked in bdr_marker.  Also computes the L2 norm of
+// alpha*n.Grad(sol) + beta*sol - gamma  over the same boundary.
+double IntegrateBC(const GridFunction &sol, const Array<int> &bdr_marker,
+                   double alpha, double beta, double gamma,
+                   double &err);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int ser_ref_levels = 2;
+   int order = 1;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool h1 = true;
+   bool visualization = true;
+
+   double mat_val = 1.0;
+   double dbc_val = 0.0;
+   double nbc_val = 1.0;
+   double rbc_a_val = 1.0; // du/dn + a * u = b
+   double rbc_b_val = 1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&h1, "-h1", "--continuous", "-dg", "--discontinuous",
+                  "Select continuous \"H1\" or discontinuous \"DG\" basis.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&mat_val, "-mat", "--material-value",
+                  "Constant value for material coefficient "
+                  "in the Laplace operator.");
+   args.AddOption(&dbc_val, "-dbc", "--dirichlet-value",
+                  "Constant value for Dirichlet Boundary Condition.");
+   args.AddOption(&nbc_val, "-nbc", "--neumann-value",
+                  "Constant value for Neumann Boundary Condition.");
+   args.AddOption(&rbc_a_val, "-rbc-a", "--robin-a-value",
+                  "Constant 'a' value for Robin Boundary Condition: "
+                  "du/dn + a * u = b.");
+   args.AddOption(&rbc_b_val, "-rbc-b", "--robin-b-value",
+                  "Constant 'b' value for Robin Boundary Condition: "
+                  "du/dn + a * u = b.");
+   args.AddOption(&a_, "-a", "--radius",
+                  "Radius of holes in the mesh.");
+   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;
+   }
+   if (kappa < 0 && !h1)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   args.PrintOptions(mfem::out);
+
+   if (a_ < 0.01)
+   {
+      mfem::out << "Hole radius too small, resetting to 0.01.\n";
+      a_ = 0.01;
+   }
+   if (a_ > 0.49)
+   {
+      mfem::out << "Hole radius too large, resetting to 0.49.\n";
+      a_ = 0.49;
+   }
+
+   // 2. Construct the (serial) mesh and refine it if requested.
+   Mesh *mesh = GenerateSerialMesh(ser_ref_levels);
+   int dim = mesh->Dimension();
+
+   // 3. Define a finite element space on the serial mesh.  Here we
+   //    use either continuous Lagrange finite elements or discontinuous
+   //    Galerkin finite elements of the specified order.
+   FiniteElementCollection *fec =
+      h1 ? (FiniteElementCollection*)new H1_FECollection(order, dim) :
+      (FiniteElementCollection*)new DG_FECollection(order, dim);
+   FiniteElementSpace fespace(mesh, fec);
+   int size = fespace.GetTrueVSize();
+   mfem::out << "Number of finite element unknowns: " << size << endl;
+
+   // 4. Create "marker arrays" to define the portions of the boundary
+   //    associated with each type of boundary condition.  These arrays
+   //    have an entry corresponding to each boundary attribute.
+   //    Placing a '1' in entry i marks attribute i+1 as being
+   //    active, '0' is inactive.
+   Array<int> nbc_bdr(mesh->bdr_attributes.Max());
+   Array<int> rbc_bdr(mesh->bdr_attributes.Max());
+   Array<int> dbc_bdr(mesh->bdr_attributes.Max());
+
+   nbc_bdr = 0; nbc_bdr[0] = 1;
+   rbc_bdr = 0; rbc_bdr[1] = 1;
+   dbc_bdr = 0; dbc_bdr[2] = 1;
+
+   Array<int> ess_tdof_list(0);
+   if (h1 && mesh->bdr_attributes.Size())
+   {
+      // For a continuous basis the linear system must be modifed to enforce
+      // an essential (Dirichlet) boundary condition.  In the DG case this is
+      // not necessary as the boundary condition will only be enforced weakly.
+      fespace.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list);
+   }
+
+   // 5. Setup the various coefficients needed for the Laplace operator and
+   //    the various boundary conditions.  In general these coefficients could
+   //    be functions of position but here we use only constants.
+   ConstantCoefficient matCoef(mat_val);
+   ConstantCoefficient dbcCoef(dbc_val);
+   ConstantCoefficient nbcCoef(nbc_val);
+   ConstantCoefficient rbcACoef(rbc_a_val);
+   ConstantCoefficient rbcBCoef(rbc_b_val);
+
+   // Since the n.Grad(u) terms arise by integrating -Div(m Grad(u)) by parts
+   // we must introduce the coefficient 'm' into the boundary conditions.
+   // Therefore, in the case of the Neumann BC, we actually enforce
+   // m n.Grad(u) = m g rather than simply n.Grad(u) = g.
+   ProductCoefficient m_nbcCoef(matCoef, nbcCoef);
+   ProductCoefficient m_rbcACoef(matCoef, rbcACoef);
+   ProductCoefficient m_rbcBCoef(matCoef, rbcBCoef);
+
+   // 6. Define the solution vector u as a finite element grid function
+   //    corresponding to fespace. Initialize u with initial guess of zero.
+   GridFunction u(&fespace);
+   u = 0.0;
+
+   // 7. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //    domain integrator.
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(matCoef));
+   if (h1)
+   {
+      // Add a Mass integrator on the Robin boundary
+      a.AddBoundaryIntegrator(new MassIntegrator(m_rbcACoef), rbc_bdr);
+   }
+   else
+   {
+      // Add the interfacial portion of the Lapalce operator
+      a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(matCoef,
+                                                            sigma, kappa));
+
+      // Counteract the n.Grad(u) term on the Dirichlet portion of the boundary
+      a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(matCoef, sigma, kappa),
+                             dbc_bdr);
+
+      // Augment the n.Grad(u) term with a*u on the Robin portion of boundary
+      a.AddBdrFaceIntegrator(new BoundaryMassIntegrator(m_rbcACoef),
+                             rbc_bdr);
+   }
+   a.Assemble();
+
+   // 8. Assemble the linear form for the right hand side vector.
+   LinearForm b(&fespace);
+
+   if (h1)
+   {
+      // Set the Dirchlet values in the solution vector
+      u.ProjectBdrCoefficient(dbcCoef, dbc_bdr);
+
+      // Add the desired value for n.Grad(u) on the Neumann boundary
+      b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_nbcCoef), nbc_bdr);
+
+      // Add the desired value for n.Grad(u) + a*u on the Robin boundary
+      b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_rbcBCoef), rbc_bdr);
+   }
+   else
+   {
+      // Add the desired value for the Dirchlet boundary
+      b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(dbcCoef, matCoef,
+                                                         sigma, kappa),
+                             dbc_bdr);
+
+      // Add the desired value for n.Grad(u) on the Neumann boundary
+      b.AddBdrFaceIntegrator(new BoundaryLFIntegrator(m_nbcCoef),
+                             nbc_bdr);
+
+      // Add the desired value for n.Grad(u) + a*u on the Robin boundary
+      b.AddBdrFaceIntegrator(new BoundaryLFIntegrator(m_rbcBCoef),
+                             rbc_bdr);
+   }
+   b.Assemble();
+
+   // 9. Construct the linear system.
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
+   //     solve the system AX=B with PCG in the symmetric case, and GMRES in the
+   //     non-symmetric one.
+   {
+      GSSmoother M((SparseMatrix&)(*A));
+      if (sigma == -1.0)
+      {
+         PCG(*A, M, B, X, 1, 500, 1e-12, 0.0);
+      }
+      else
+      {
+         GMRES(*A, M, B, X, 1, 500, 10, 1e-12, 0.0);
+      }
+   }
+#else
+   // 11. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
+   //     system.
+   UMFPackSolver umf_solver;
+   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   umf_solver.SetOperator(*A);
+   umf_solver.Mult(B, X);
+#endif
+
+   // 12. Recover the grid function corresponding to U. This is the
+   //     local finite element solution.
+   a.RecoverFEMSolution(X, b, u);
+
+   // 13. Build a mass matrix to help solve for n.Grad(u) where 'n' is
+   //     a surface normal.
+   BilinearForm m(&fespace);
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+
+   ess_tdof_list.SetSize(0);
+   OperatorPtr M;
+   m.FormSystemMatrix(ess_tdof_list, M);
+
+   // 14. Compute the various boundary integrals.
+   mfem::out << endl
+             << "Verifying boundary conditions" << endl
+             << "=============================" << endl;
+   {
+      // Integrate the solution on the Dirichlet boundary and compare
+      // to the expected value.
+      double err, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, err);
+
+      bool hom_dbc = (dbc_val == 0.0);
+      err /=  hom_dbc ? 1.0 : fabs(dbc_val);
+      mfem::out << "Average of solution on Gamma_dbc:\t"
+                << avg << ", \t"
+                << (hom_dbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+   {
+      // Integrate n.Grad(u) on the inhomogeneous Neumann boundary and
+      // compare to the expected value.
+      double err, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, err);
+
+      bool hom_nbc = (nbc_val == 0.0);
+      err /=  hom_nbc ? 1.0 : fabs(nbc_val);
+      mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
+                << avg << ", \t"
+                << (hom_nbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+   {
+      // Integrate n.Grad(u) on the homogeneous Neumann boundary and compare
+      // to the expected value of zero.
+      Array<int> nbc0_bdr(mesh->bdr_attributes.Max());
+      nbc0_bdr = 0;
+      nbc0_bdr[3] = 1;
+
+      double err, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, err);
+
+      bool hom_nbc = true;
+      mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
+                << avg << ", \t"
+                << (hom_nbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+   {
+      // Integrate n.Grad(u) + a * u on the Robin boundary and compare to
+      // the expected value.
+      double err, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, err);
+
+      bool hom_rbc = (rbc_b_val == 0.0);
+      err /=  hom_rbc ? 1.0 : fabs(rbc_b_val);
+      mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
+                << avg << ", \t"
+                << (hom_rbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+
+   // 15. Save the refined mesh and the solution. This output can be viewed
+   //     later using GLVis: "glvis -m refined.mesh -g sol.gf".
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      u.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      string title_str = h1 ? "H1" : "DG";
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << u
+               << "window_title '" << title_str << " Solution'"
+               << " keys 'mmc'" << flush;
+   }
+
+   // 17. Free the used memory.
+   delete fec;
+   delete mesh;
+
+   return 0;
+}
+
+void quad_trans(double u, double v, double &x, double &y, bool log = false)
+{
+   double a = a_; // Radius of disc
+
+   double d = 4.0 * a * (M_SQRT2 - 2.0 * a) * (1.0 - 2.0 * v);
+
+   double v0 = (1.0 + M_SQRT2) * (M_SQRT2 * a - 2.0 * v) *
+               ((4.0 - 3 * M_SQRT2) * a +
+                (8.0 * (M_SQRT2 - 1.0) * a - 2.0) * v) / d;
+
+   double r = 2.0 * ((M_SQRT2 - 1.0) * a * a * (1.0 - 4.0 *v) +
+                     2.0 * (1.0 + M_SQRT2 *
+                            (1.0 + 2.0 * (2.0 * a - M_SQRT2 - 1.0) * a)) * v * v
+                    ) / d;
+
+   double t = asin(v / r) * u / v;
+   if (log)
+   {
+      mfem::out << "u, v, r, v0, t "
+                << u << " " << v << " " << r << " " << v0 << " " << t
+                << endl;
+   }
+   x = r * sin(t);
+   y = r * cos(t) - v0;
+}
+
+void trans(const Vector &u, Vector &x)
+{
+   double tol = 1e-4;
+
+   if (u[1] > 0.5 - tol || u[1] < -0.5 + tol)
+   {
+      x = u;
+      return;
+   }
+   if (u[0] > 1.0 - tol || u[0] < -1.0 + tol || fabs(u[0]) < tol)
+   {
+      x = u;
+      return;
+   }
+
+   if (u[0] > 0.0)
+   {
+      if (u[1] > fabs(u[0] - 0.5))
+      {
+         quad_trans(u[0] - 0.5, u[1], x[0], x[1]);
+         x[0] += 0.5;
+         return;
+      }
+      if (u[1] < -fabs(u[0] - 0.5))
+      {
+         quad_trans(u[0] - 0.5, -u[1], x[0], x[1]);
+         x[0] += 0.5;
+         x[1] *= -1.0;
+         return;
+      }
+      if (u[0] - 0.5 > fabs(u[1]))
+      {
+         quad_trans(u[1], u[0] - 0.5, x[1], x[0]);
+         x[0] += 0.5;
+         return;
+      }
+      if (u[0] - 0.5 < -fabs(u[1]))
+      {
+         quad_trans(u[1], 0.5 - u[0], x[1], x[0]);
+         x[0] *= -1.0;
+         x[0] += 0.5;
+         return;
+      }
+   }
+   else
+   {
+      if (u[1] > fabs(u[0] + 0.5))
+      {
+         quad_trans(u[0] + 0.5, u[1], x[0], x[1]);
+         x[0] -= 0.5;
+         return;
+      }
+      if (u[1] < -fabs(u[0] + 0.5))
+      {
+         quad_trans(u[0] + 0.5, -u[1], x[0], x[1]);
+         x[0] -= 0.5;
+         x[1] *= -1.0;
+         return;
+      }
+      if (u[0] + 0.5 > fabs(u[1]))
+      {
+         quad_trans(u[1], u[0] + 0.5, x[1], x[0]);
+         x[0] -= 0.5;
+         return;
+      }
+      if (u[0] + 0.5 < -fabs(u[1]))
+      {
+         quad_trans(u[1], -0.5 - u[0], x[1], x[0]);
+         x[0] *= -1.0;
+         x[0] -= 0.5;
+         return;
+      }
+   }
+   x = u;
+}
+
+Mesh * GenerateSerialMesh(int ref)
+{
+   Mesh * mesh = new Mesh(2, 29, 16, 24, 2);
+
+   int vi[4];
+
+   for (int i=0; i<2; i++)
+   {
+      int o = 13 * i;
+      vi[0] = o + 0; vi[1] = o + 3; vi[2] = o + 4; vi[3] = o + 1;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 1; vi[1] = o + 4; vi[2] = o + 5; vi[3] = o + 2;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 5; vi[1] = o + 8; vi[2] = o + 9; vi[3] = o + 2;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 8; vi[1] = o + 12; vi[2] = o + 15; vi[3] = o + 9;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 11; vi[1] = o + 14; vi[2] = o + 15; vi[3] = o + 12;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 10; vi[1] = o + 13; vi[2] = o + 14; vi[3] = o + 11;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 6; vi[1] = o + 13; vi[2] = o + 10; vi[3] = o + 7;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 0; vi[1] = o + 6; vi[2] = o + 7; vi[3] = o + 3;
+      mesh->AddQuad(vi);
+   }
+
+   vi[0] =  0; vi[1] =  6; mesh->AddBdrSegment(vi, 1);
+   vi[0] =  6; vi[1] = 13; mesh->AddBdrSegment(vi, 1);
+   vi[0] = 13; vi[1] = 19; mesh->AddBdrSegment(vi, 1);
+   vi[0] = 19; vi[1] = 26; mesh->AddBdrSegment(vi, 1);
+
+   vi[0] = 28; vi[1] = 22; mesh->AddBdrSegment(vi, 2);
+   vi[0] = 22; vi[1] = 15; mesh->AddBdrSegment(vi, 2);
+   vi[0] = 15; vi[1] =  9; mesh->AddBdrSegment(vi, 2);
+   vi[0] =  9; vi[1] =  2; mesh->AddBdrSegment(vi, 2);
+
+   for (int i=0; i<2; i++)
+   {
+      int o = 13 * i;
+      vi[0] = o +  7; vi[1] = o +  3; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o + 10; vi[1] = o +  7; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o + 11; vi[1] = o + 10; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o + 12; vi[1] = o + 11; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  8; vi[1] = o + 12; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  5; vi[1] = o +  8; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  4; vi[1] = o +  5; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  3; vi[1] = o +  4; mesh->AddBdrSegment(vi, 3 + i);
+   }
+
+   double d[2];
+   double a = a_ / M_SQRT2;
+
+   d[0] = -1.0; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] = -1.0; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] = -1.0; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] = -0.5 - a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] = -0.5 - a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] = -0.5 - a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] = -0.5; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] = -0.5; d[1] =   -a; mesh->AddVertex(d);
+   d[0] = -0.5; d[1] =    a; mesh->AddVertex(d);
+   d[0] = -0.5; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] = -0.5 + a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] = -0.5 + a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] = -0.5 + a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] =  0.0; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] =  0.0; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  0.0; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] =  0.5 - a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] =  0.5 - a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  0.5 - a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] =  0.5; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] =  0.5; d[1] =   -a; mesh->AddVertex(d);
+   d[0] =  0.5; d[1] =    a; mesh->AddVertex(d);
+   d[0] =  0.5; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] =  0.5 + a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] =  0.5 + a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  0.5 + a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] =  1.0; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] =  1.0; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  1.0; d[1] =  0.5; mesh->AddVertex(d);
+
+   mesh->FinalizeTopology();
+
+   mesh->SetCurvature(1, true);
+
+   // Stitch the ends of the stack together
+   {
+      Array<int> v2v(mesh->GetNV());
+      for (int i = 0; i < v2v.Size() - 3; i++)
+      {
+         v2v[i] = i;
+      }
+      // identify vertices on the narrow ends of the rectangle
+      v2v[v2v.Size() - 3] = 0;
+      v2v[v2v.Size() - 2] = 1;
+      v2v[v2v.Size() - 1] = 2;
+
+      // renumber elements
+      for (int i = 0; i < mesh->GetNE(); i++)
+      {
+         Element *el = mesh->GetElement(i);
+         int *v = el->GetVertices();
+         int nv = el->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            v[j] = v2v[v[j]];
+         }
+      }
+      // renumber boundary elements
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         Element *el = mesh->GetBdrElement(i);
+         int *v = el->GetVertices();
+         int nv = el->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            v[j] = v2v[v[j]];
+         }
+      }
+      mesh->RemoveUnusedVertices();
+      mesh->RemoveInternalBoundaries();
+   }
+   mesh->SetCurvature(3, true);
+
+   for (int l = 0; l < ref; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   mesh->Transform(trans);
+
+   return mesh;
+}
+
+double IntegrateBC(const GridFunction &x, const Array<int> &bdr,
+                   double alpha, double beta, double gamma,
+                   double &err)
+{
+   double nrm = 0.0;
+   double avg = 0.0;
+   err = 0.0;
+
+   const bool a_is_zero = alpha == 0.0;
+   const bool b_is_zero = beta == 0.0;
+
+   const FiniteElementSpace &fes = *x.FESpace();
+   MFEM_ASSERT(fes.GetVDim() == 1, "");
+   Mesh &mesh = *fes.GetMesh();
+   Vector shape, loc_dofs, w_nor;
+   DenseMatrix dshape;
+   Array<int> dof_ids;
+   for (int i = 0; i < mesh.GetNBE(); i++)
+   {
+      if (bdr[mesh.GetBdrAttribute(i)-1] == 0) { continue; }
+
+      FaceElementTransformations *FTr = mesh.GetBdrFaceTransformations(i);
+      if (FTr == nullptr) { continue; }
+
+      const FiniteElement &fe = *fes.GetFE(FTr->Elem1No);
+      MFEM_ASSERT(fe.GetMapType() == FiniteElement::VALUE, "");
+      const int int_order = 2*fe.GetOrder() + 3;
+      const IntegrationRule &ir = IntRules.Get(FTr->FaceGeom, int_order);
+
+      fes.GetElementDofs(FTr->Elem1No, dof_ids);
+      x.GetSubVector(dof_ids, loc_dofs);
+      if (!a_is_zero)
+      {
+         const int sdim = FTr->Face->GetSpaceDim();
+         w_nor.SetSize(sdim);
+         dshape.SetSize(fe.GetDof(), sdim);
+      }
+      if (!b_is_zero)
+      {
+         shape.SetSize(fe.GetDof());
+      }
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+         IntegrationPoint eip;
+         FTr->Loc1.Transform(ip, eip);
+         FTr->Face->SetIntPoint(&ip);
+         double face_weight = FTr->Face->Weight();
+         double val = 0.0;
+         if (!a_is_zero)
+         {
+            FTr->Elem1->SetIntPoint(&eip);
+            fe.CalcPhysDShape(*FTr->Elem1, dshape);
+            CalcOrtho(FTr->Face->Jacobian(), w_nor);
+            val += alpha * dshape.InnerProduct(w_nor, loc_dofs) / face_weight;
+         }
+         if (!b_is_zero)
+         {
+            fe.CalcShape(eip, shape);
+            val += beta * (shape * loc_dofs);
+         }
+
+         // Measure the length of the boundary
+         nrm += ip.weight * face_weight;
+
+         // Integrate alpha * n.Grad(x) + beta * x
+         avg += val * ip.weight * face_weight;
+
+         // Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
+         val -= gamma;
+         err += (val*val) * ip.weight * face_weight;
+      }
+   }
+
+   // Normalize by the length of the boundary
+   if (std::abs(nrm) > 0.0)
+   {
+      err /= nrm;
+      avg /= nrm;
+   }
+
+   // Compute l2 norm of the error in the boundary condition
+   // (negative quadrature weights may produce negative 'err')
+   err = (err >= 0.0) ? sqrt(err) : -sqrt(-err);
+
+   // Return the average value of alpha * n.Grad(x) + beta * x
+   return avg;
+}

examples/ex27p.cpp

@@ -0,0 +1,773 @@
+//                       MFEM Example 27 - Parallel Version
+//
+// Compile with: make ex27p
+//
+// Sample runs:  mpirun -np 4 ex27p
+//               mpirun -np 4 ex27p -dg
+//               mpirun -np 4 ex27p -dg -dbc 8 -nbc -2
+//               mpirun -np 4 ex27p -rbc-a 1 -rbc-b 8
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 0 with a variety of boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order using a continuous or discontinuous space.  We then
+//               apply Dirichlet, Neumann (both homogeneous and inhomogeneous),
+//               Robin, and Periodic boundary conditions on different portions
+//               of a predefined mesh.
+//
+//               The predefined mesh consists of a rectangle with two
+//               holes removed (see below).  The narrow ends of the
+//               mesh are connected to form a Periodic boundary
+//               condition.  The lower edge (tagged with attribute 1)
+//               receives an inhomogeneous Neumann boundary condition.
+//               A Robin boundary condition is applied to upper edge
+//               (attribute 2).  The circular hole on the left
+//               (attribute 3) enforces a Dirichlet boundary
+//               condition.  Finally, a natural boundary condition, or
+//               homogeneous Neumann BC, is applied to the circular
+//               hole on the right (attribute 4).
+//
+//                  Attribute 3    ^ y  Attribute 2
+//                        \        |      /
+//                     +-----------+-----------+
+//                     |    \_     |     _     |
+//                     |    / \    |    / \    |
+//                  <--+---+---+---+---+---+---+--> x
+//                     |    \_/    |    \_/    |
+//                     |           |      \    |
+//                     +-----------+-----------+       (hole radii are
+//                          /      |        \            adjustable)
+//                  Attribute 1    v    Attribute 4
+//
+//               The boundary conditions are defined as (where u is
+//               the solution field):
+//                  Dirichlet: u = d
+//                  Neumann:   n.Grad(u) = g
+//                  Robin:     n.Grad(u) + a u = b
+//
+//               The user can adjust the values of 'd', 'g', 'a', and
+//               'b' with command line options.
+//
+//               This example highlights the differing implementations of
+//               boundary conditions with continuous and discontinuous Galerkin
+//               formulations of the Laplace problem.
+//
+//               We recommend viewing examples 1 and 14 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double a_ = 0.2;
+
+// Normal to hole with boundary attribute 4
+void n4Vec(const Vector &x, Vector &n) { n = x; n[0] -= 0.5; n /= -n.Norml2(); }
+
+Mesh * GenerateSerialMesh(int ref);
+
+// Compute the average value of alpha*n.Grad(sol) + beta*sol over the boundary
+// attributes marked in bdr_marker.  Also computes the L2 norm of
+// alpha*n.Grad(sol) + beta*sol - gamma  over the same boundary.
+double IntegrateBC(const ParGridFunction &sol, const Array<int> &bdr_marker,
+                   double alpha, double beta, double gamma,
+                   double &err);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   MPI_Session mpi;
+   if (!mpi.Root()) { mfem::out.Disable(); mfem::err.Disable(); }
+
+   // 2. Parse command-line options.
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   int order = 1;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool h1 = true;
+   bool visualization = true;
+
+   double mat_val = 1.0;
+   double dbc_val = 0.0;
+   double nbc_val = 1.0;
+   double rbc_a_val = 1.0; // du/dn + a * u = b
+   double rbc_b_val = 1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&h1, "-h1", "--continuous", "-dg", "--discontinuous",
+                  "Select continuous \"H1\" or discontinuous \"DG\" basis.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   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(&mat_val, "-mat", "--material-value",
+                  "Constant value for material coefficient "
+                  "in the Laplace operator.");
+   args.AddOption(&dbc_val, "-dbc", "--dirichlet-value",
+                  "Constant value for Dirichlet Boundary Condition.");
+   args.AddOption(&nbc_val, "-nbc", "--neumann-value",
+                  "Constant value for Neumann Boundary Condition.");
+   args.AddOption(&rbc_a_val, "-rbc-a", "--robin-a-value",
+                  "Constant 'a' value for Robin Boundary Condition: "
+                  "du/dn + a * u = b.");
+   args.AddOption(&rbc_b_val, "-rbc-b", "--robin-b-value",
+                  "Constant 'b' value for Robin Boundary Condition: "
+                  "du/dn + a * u = b.");
+   args.AddOption(&a_, "-a", "--radius",
+                  "Radius of holes in the mesh.");
+   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;
+   }
+   if (kappa < 0 && !h1)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   args.PrintOptions(mfem::out);
+
+   if (a_ < 0.01)
+   {
+      mfem::out << "Hole radius too small, resetting to 0.01.\n";
+      a_ = 0.01;
+   }
+   if (a_ > 0.49)
+   {
+      mfem::out << "Hole radius too large, resetting to 0.49.\n";
+      a_ = 0.49;
+   }
+
+   // 3. Construct the (serial) mesh and refine it if requested.
+   Mesh *mesh = GenerateSerialMesh(ser_ref_levels);
+   int dim = mesh->Dimension();
+
+   // 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int l = 0; l < par_ref_levels; l++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   // 5. Define a parallel finite element space on the parallel mesh. Here we
+   //    use either continuous Lagrange finite elements or discontinuous
+   //    Galerkin finite elements of the specified order.
+   FiniteElementCollection *fec =
+      h1 ? (FiniteElementCollection*)new H1_FECollection(order, dim) :
+      (FiniteElementCollection*)new DG_FECollection(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, fec);
+   HYPRE_Int size = fespace.GlobalTrueVSize();
+   mfem::out << "Number of finite element unknowns: " << size << endl;
+
+   // 6. Create "marker arrays" to define the portions of the boundary
+   //    associated with each type of boundary condition.  These arrays
+   //    have an entry corresponding to each boundary attribute.
+   //    Placing a '1' in entry i marks attribute i+1 as being
+   //    active, '0' is inactive.
+   Array<int> nbc_bdr(pmesh.bdr_attributes.Max());
+   Array<int> rbc_bdr(pmesh.bdr_attributes.Max());
+   Array<int> dbc_bdr(pmesh.bdr_attributes.Max());
+
+   nbc_bdr = 0; nbc_bdr[0] = 1;
+   rbc_bdr = 0; rbc_bdr[1] = 1;
+   dbc_bdr = 0; dbc_bdr[2] = 1;
+
+   Array<int> ess_tdof_list(0);
+   if (h1 && pmesh.bdr_attributes.Size())
+   {
+      // For a continuous basis the linear system must be modifed to enforce
+      // an essential (Dirichlet) boundary condition.  In the DG case this is
+      // not necessary as the boundary condition will only be enforced weakly.
+      fespace.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list);
+   }
+
+   // 7. Setup the various coefficients needed for the Laplace operator and
+   //    the various boundary conditions.  In general these coefficients could
+   //    be functions of position but here we use only constants.
+   ConstantCoefficient matCoef(mat_val);
+   ConstantCoefficient dbcCoef(dbc_val);
+   ConstantCoefficient nbcCoef(nbc_val);
+   ConstantCoefficient rbcACoef(rbc_a_val);
+   ConstantCoefficient rbcBCoef(rbc_b_val);
+
+   // Since the n.Grad(u) terms arise by integrating -Div(m Grad(u)) by parts
+   // we must introduce the coefficient 'm' into the boundary conditions.
+   // Therefore, in the case of the Neumann BC, we actually enforce
+   // m n.Grad(u) = m g rather than simply n.Grad(u) = g.
+   ProductCoefficient m_nbcCoef(matCoef, nbcCoef);
+   ProductCoefficient m_rbcACoef(matCoef, rbcACoef);
+   ProductCoefficient m_rbcBCoef(matCoef, rbcBCoef);
+
+   // 8. Define the solution vector u as a parallel finite element grid function
+   //    corresponding to fespace. Initialize u with initial guess of zero.
+   ParGridFunction u(&fespace);
+   u = 0.0;
+
+   // 9. Set up the parallel bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //    domain integrator.
+   ParBilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(matCoef));
+   if (h1)
+   {
+      // Add a Mass integrator on the Robin boundary
+      a.AddBoundaryIntegrator(new MassIntegrator(m_rbcACoef), rbc_bdr);
+   }
+   else
+   {
+      // Add the interfacial portion of the Lapalce operator
+      a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(matCoef,
+                                                            sigma, kappa));
+
+      // Counteract the n.Grad(u) term on the Dirichlet portion of the boundary
+      a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(matCoef, sigma, kappa),
+                             dbc_bdr);
+
+      // Augment the n.Grad(u) term with a*u on the Robin portion of boundary
+      a.AddBdrFaceIntegrator(new BoundaryMassIntegrator(m_rbcACoef),
+                             rbc_bdr);
+   }
+   a.Assemble();
+
+   // 10. Assemble the parallel linear form for the right hand side vector.
+   ParLinearForm b(&fespace);
+
+   if (h1)
+   {
+      // Set the Dirchlet values in the solution vector
+      u.ProjectBdrCoefficient(dbcCoef, dbc_bdr);
+
+      // Add the desired value for n.Grad(u) on the Neumann boundary
+      b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_nbcCoef), nbc_bdr);
+
+      // Add the desired value for n.Grad(u) + a*u on the Robin boundary
+      b.AddBoundaryIntegrator(new BoundaryLFIntegrator(m_rbcBCoef), rbc_bdr);
+   }
+   else
+   {
+      // Add the desired value for the Dirchlet boundary
+      b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(dbcCoef, matCoef,
+                                                         sigma, kappa),
+                             dbc_bdr);
+
+      // Add the desired value for n.Grad(u) on the Neumann boundary
+      b.AddBdrFaceIntegrator(new BoundaryLFIntegrator(m_nbcCoef),
+                             nbc_bdr);
+
+      // Add the desired value for n.Grad(u) + a*u on the Robin boundary
+      b.AddBdrFaceIntegrator(new BoundaryLFIntegrator(m_rbcBCoef),
+                             rbc_bdr);
+   }
+   b.Assemble();
+
+   // 11. Construct the linear system.
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
+
+   // 12. Solve the linear system A X = B.
+   HypreSolver *amg = new HypreBoomerAMG;
+   if (h1 || sigma == -1.0)
+   {
+      HyprePCG pcg(MPI_COMM_WORLD);
+      pcg.SetTol(1e-12);
+      pcg.SetMaxIter(200);
+      pcg.SetPrintLevel(2);
+      pcg.SetPreconditioner(*amg);
+      pcg.SetOperator(*A);
+      pcg.Mult(B, X);
+   }
+   else
+   {
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetAbsTol(0.0);
+      gmres.SetRelTol(1e-12);
+      gmres.SetMaxIter(200);
+      gmres.SetKDim(10);
+      gmres.SetPrintLevel(1);
+      gmres.SetPreconditioner(*amg);
+      gmres.SetOperator(*A);
+      gmres.Mult(B, X);
+   }
+   delete amg;
+
+   // 13. Recover the parallel grid function corresponding to U. This is the
+   //     local finite element solution on each processor.
+   a.RecoverFEMSolution(X, b, u);
+
+   // 14. Build a mass matrix to help solve for n.Grad(u) where 'n' is
+   //     a surface normal.
+   ParBilinearForm m(&fespace);
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+
+   ess_tdof_list.SetSize(0);
+   OperatorPtr M;
+   m.FormSystemMatrix(ess_tdof_list, M);
+
+   // 15. Compute the various boundary integrals.
+   mfem::out << endl
+             << "Verifying boundary conditions" << endl
+             << "=============================" << endl;
+   {
+      // Integrate the solution on the Dirichlet boundary and compare
+      // to the expected value.
+      double err, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, err);
+
+      bool hom_dbc = (dbc_val == 0.0);
+      err /=  hom_dbc ? 1.0 : fabs(dbc_val);
+      mfem::out << "Average of solution on Gamma_dbc:\t"
+                << avg << ", \t"
+                << (hom_dbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+   {
+      // Integrate n.Grad(u) on the inhomogeneous Neumann boundary and
+      // compare to the expected value.
+      double err, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, err);
+
+      bool hom_nbc = (nbc_val == 0.0);
+      err /=  hom_nbc ? 1.0 : fabs(nbc_val);
+      mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
+                << avg << ", \t"
+                << (hom_nbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+   {
+      // Integrate n.Grad(u) on the homogeneous Neumann boundary and compare
+      // to the expected value of zero.
+      Array<int> nbc0_bdr(pmesh.bdr_attributes.Max());
+      nbc0_bdr = 0;
+      nbc0_bdr[3] = 1;
+
+      double err, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, err);
+
+      bool hom_nbc = true;
+      mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
+                << avg << ", \t"
+                << (hom_nbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+   {
+      // Integrate n.Grad(u) + a * u on the Robin boundary and compare to
+      // the expected value.
+      double err, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, err);
+
+      bool hom_rbc = (rbc_b_val == 0.0);
+      err /=  hom_rbc ? 1.0 : fabs(rbc_b_val);
+      mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
+                << avg << ", \t"
+                << (hom_rbc ? "absolute" : "relative")
+                << " error " << err << endl;
+   }
+
+   // 16. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << mpi.WorldRank();
+      sol_name << "sol." << setfill('0') << setw(6) << mpi.WorldRank();
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh.Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      u.Save(sol_ofs);
+   }
+
+   // 17. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      string title_str = h1 ? "H1" : "DG";
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << mpi.WorldSize()
+               << " " << mpi.WorldRank() << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << u
+               << "window_title '" << title_str << " Solution'"
+               << " keys 'mmc'" << flush;
+   }
+
+   // 18. Free the used memory.
+   delete fec;
+
+   return 0;
+}
+
+void quad_trans(double u, double v, double &x, double &y, bool log = false)
+{
+   double a = a_; // Radius of disc
+
+   double d = 4.0 * a * (M_SQRT2 - 2.0 * a) * (1.0 - 2.0 * v);
+
+   double v0 = (1.0 + M_SQRT2) * (M_SQRT2 * a - 2.0 * v) *
+               ((4.0 - 3 * M_SQRT2) * a +
+                (8.0 * (M_SQRT2 - 1.0) * a - 2.0) * v) / d;
+
+   double r = 2.0 * ((M_SQRT2 - 1.0) * a * a * (1.0 - 4.0 *v) +
+                     2.0 * (1.0 + M_SQRT2 *
+                            (1.0 + 2.0 * (2.0 * a - M_SQRT2 - 1.0) * a)) * v * v
+                    ) / d;
+
+   double t = asin(v / r) * u / v;
+   if (log)
+   {
+      mfem::out << "u, v, r, v0, t "
+                << u << " " << v << " " << r << " " << v0 << " " << t
+                << endl;
+   }
+   x = r * sin(t);
+   y = r * cos(t) - v0;
+}
+
+void trans(const Vector &u, Vector &x)
+{
+   double tol = 1e-4;
+
+   if (u[1] > 0.5 - tol || u[1] < -0.5 + tol)
+   {
+      x = u;
+      return;
+   }
+   if (u[0] > 1.0 - tol || u[0] < -1.0 + tol || fabs(u[0]) < tol)
+   {
+      x = u;
+      return;
+   }
+
+   if (u[0] > 0.0)
+   {
+      if (u[1] > fabs(u[0] - 0.5))
+      {
+         quad_trans(u[0] - 0.5, u[1], x[0], x[1]);
+         x[0] += 0.5;
+         return;
+      }
+      if (u[1] < -fabs(u[0] - 0.5))
+      {
+         quad_trans(u[0] - 0.5, -u[1], x[0], x[1]);
+         x[0] += 0.5;
+         x[1] *= -1.0;
+         return;
+      }
+      if (u[0] - 0.5 > fabs(u[1]))
+      {
+         quad_trans(u[1], u[0] - 0.5, x[1], x[0]);
+         x[0] += 0.5;
+         return;
+      }
+      if (u[0] - 0.5 < -fabs(u[1]))
+      {
+         quad_trans(u[1], 0.5 - u[0], x[1], x[0]);
+         x[0] *= -1.0;
+         x[0] += 0.5;
+         return;
+      }
+   }
+   else
+   {
+      if (u[1] > fabs(u[0] + 0.5))
+      {
+         quad_trans(u[0] + 0.5, u[1], x[0], x[1]);
+         x[0] -= 0.5;
+         return;
+      }
+      if (u[1] < -fabs(u[0] + 0.5))
+      {
+         quad_trans(u[0] + 0.5, -u[1], x[0], x[1]);
+         x[0] -= 0.5;
+         x[1] *= -1.0;
+         return;
+      }
+      if (u[0] + 0.5 > fabs(u[1]))
+      {
+         quad_trans(u[1], u[0] + 0.5, x[1], x[0]);
+         x[0] -= 0.5;
+         return;
+      }
+      if (u[0] + 0.5 < -fabs(u[1]))
+      {
+         quad_trans(u[1], -0.5 - u[0], x[1], x[0]);
+         x[0] *= -1.0;
+         x[0] -= 0.5;
+         return;
+      }
+   }
+   x = u;
+}
+
+Mesh * GenerateSerialMesh(int ref)
+{
+   Mesh * mesh = new Mesh(2, 29, 16, 24, 2);
+
+   int vi[4];
+
+   for (int i=0; i<2; i++)
+   {
+      int o = 13 * i;
+      vi[0] = o + 0; vi[1] = o + 3; vi[2] = o + 4; vi[3] = o + 1;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 1; vi[1] = o + 4; vi[2] = o + 5; vi[3] = o + 2;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 5; vi[1] = o + 8; vi[2] = o + 9; vi[3] = o + 2;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 8; vi[1] = o + 12; vi[2] = o + 15; vi[3] = o + 9;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 11; vi[1] = o + 14; vi[2] = o + 15; vi[3] = o + 12;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 10; vi[1] = o + 13; vi[2] = o + 14; vi[3] = o + 11;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 6; vi[1] = o + 13; vi[2] = o + 10; vi[3] = o + 7;
+      mesh->AddQuad(vi);
+
+      vi[0] = o + 0; vi[1] = o + 6; vi[2] = o + 7; vi[3] = o + 3;
+      mesh->AddQuad(vi);
+   }
+
+   vi[0] =  0; vi[1] =  6; mesh->AddBdrSegment(vi, 1);
+   vi[0] =  6; vi[1] = 13; mesh->AddBdrSegment(vi, 1);
+   vi[0] = 13; vi[1] = 19; mesh->AddBdrSegment(vi, 1);
+   vi[0] = 19; vi[1] = 26; mesh->AddBdrSegment(vi, 1);
+
+   vi[0] = 28; vi[1] = 22; mesh->AddBdrSegment(vi, 2);
+   vi[0] = 22; vi[1] = 15; mesh->AddBdrSegment(vi, 2);
+   vi[0] = 15; vi[1] =  9; mesh->AddBdrSegment(vi, 2);
+   vi[0] =  9; vi[1] =  2; mesh->AddBdrSegment(vi, 2);
+
+   for (int i=0; i<2; i++)
+   {
+      int o = 13 * i;
+      vi[0] = o +  7; vi[1] = o +  3; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o + 10; vi[1] = o +  7; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o + 11; vi[1] = o + 10; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o + 12; vi[1] = o + 11; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  8; vi[1] = o + 12; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  5; vi[1] = o +  8; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  4; vi[1] = o +  5; mesh->AddBdrSegment(vi, 3 + i);
+      vi[0] = o +  3; vi[1] = o +  4; mesh->AddBdrSegment(vi, 3 + i);
+   }
+
+   double d[2];
+   double a = a_ / M_SQRT2;
+
+   d[0] = -1.0; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] = -1.0; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] = -1.0; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] = -0.5 - a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] = -0.5 - a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] = -0.5 - a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] = -0.5; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] = -0.5; d[1] =   -a; mesh->AddVertex(d);
+   d[0] = -0.5; d[1] =    a; mesh->AddVertex(d);
+   d[0] = -0.5; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] = -0.5 + a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] = -0.5 + a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] = -0.5 + a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] =  0.0; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] =  0.0; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  0.0; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] =  0.5 - a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] =  0.5 - a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  0.5 - a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] =  0.5; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] =  0.5; d[1] =   -a; mesh->AddVertex(d);
+   d[0] =  0.5; d[1] =    a; mesh->AddVertex(d);
+   d[0] =  0.5; d[1] =  0.5; mesh->AddVertex(d);
+
+   d[0] =  0.5 + a; d[1] =   -a; mesh->AddVertex(d);
+   d[0] =  0.5 + a; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  0.5 + a; d[1] =    a; mesh->AddVertex(d);
+
+   d[0] =  1.0; d[1] = -0.5; mesh->AddVertex(d);
+   d[0] =  1.0; d[1] =  0.0; mesh->AddVertex(d);
+   d[0] =  1.0; d[1] =  0.5; mesh->AddVertex(d);
+
+   mesh->FinalizeTopology();
+
+   mesh->SetCurvature(1, true);
+
+   // Stitch the ends of the stack together
+   {
+      Array<int> v2v(mesh->GetNV());
+      for (int i = 0; i < v2v.Size() - 3; i++)
+      {
+         v2v[i] = i;
+      }
+      // identify vertices on the narrow ends of the rectangle
+      v2v[v2v.Size() - 3] = 0;
+      v2v[v2v.Size() - 2] = 1;
+      v2v[v2v.Size() - 1] = 2;
+
+      // renumber elements
+      for (int i = 0; i < mesh->GetNE(); i++)
+      {
+         Element *el = mesh->GetElement(i);
+         int *v = el->GetVertices();
+         int nv = el->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            v[j] = v2v[v[j]];
+         }
+      }
+      // renumber boundary elements
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         Element *el = mesh->GetBdrElement(i);
+         int *v = el->GetVertices();
+         int nv = el->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            v[j] = v2v[v[j]];
+         }
+      }
+      mesh->RemoveUnusedVertices();
+      mesh->RemoveInternalBoundaries();
+   }
+   mesh->SetCurvature(3, true);
+
+   for (int l = 0; l < ref; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   mesh->Transform(trans);
+
+   return mesh;
+}
+
+double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
+                   double alpha, double beta, double gamma,
+                   double &glb_err)
+{
+   double loc_vals[3];
+   double &nrm = loc_vals[0];
+   double &avg = loc_vals[1];
+   double &err = loc_vals[2];
+
+   nrm = 0.0;
+   avg = 0.0;
+   err = 0.0;
+
+   const bool a_is_zero = alpha == 0.0;
+   const bool b_is_zero = beta == 0.0;
+
+   const ParFiniteElementSpace &fes = *x.ParFESpace();
+   MFEM_ASSERT(fes.GetVDim() == 1, "");
+   ParMesh &mesh = *fes.GetParMesh();
+   Vector shape, loc_dofs, w_nor;
+   DenseMatrix dshape;
+   Array<int> dof_ids;
+   for (int i = 0; i < mesh.GetNBE(); i++)
+   {
+      if (bdr[mesh.GetBdrAttribute(i)-1] == 0) { continue; }
+
+      FaceElementTransformations *FTr = mesh.GetBdrFaceTransformations(i);
+      if (FTr == nullptr) { continue; }
+
+      const FiniteElement &fe = *fes.GetFE(FTr->Elem1No);
+      MFEM_ASSERT(fe.GetMapType() == FiniteElement::VALUE, "");
+      const int int_order = 2*fe.GetOrder() + 3;
+      const IntegrationRule &ir = IntRules.Get(FTr->FaceGeom, int_order);
+
+      fes.GetElementDofs(FTr->Elem1No, dof_ids);
+      x.GetSubVector(dof_ids, loc_dofs);
+      if (!a_is_zero)
+      {
+         const int sdim = FTr->Face->GetSpaceDim();
+         w_nor.SetSize(sdim);
+         dshape.SetSize(fe.GetDof(), sdim);
+      }
+      if (!b_is_zero)
+      {
+         shape.SetSize(fe.GetDof());
+      }
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+         IntegrationPoint eip;
+         FTr->Loc1.Transform(ip, eip);
+         FTr->Face->SetIntPoint(&ip);
+         double face_weight = FTr->Face->Weight();
+         double val = 0.0;
+         if (!a_is_zero)
+         {
+            FTr->Elem1->SetIntPoint(&eip);
+            fe.CalcPhysDShape(*FTr->Elem1, dshape);
+            CalcOrtho(FTr->Face->Jacobian(), w_nor);
+            val += alpha * dshape.InnerProduct(w_nor, loc_dofs) / face_weight;
+         }
+         if (!b_is_zero)
+         {
+            fe.CalcShape(eip, shape);
+            val += beta * (shape * loc_dofs);
+         }
+
+         // Measure the length of the boundary
+         nrm += ip.weight * face_weight;
+
+         // Integrate alpha * n.Grad(x) + beta * x
+         avg += val * ip.weight * face_weight;
+
+         // Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
+         val -= gamma;
+         err += (val*val) * ip.weight * face_weight;
+      }
+   }
+
+   double glb_vals[3];
+   MPI_Allreduce(loc_vals, glb_vals, 3, MPI_DOUBLE, MPI_SUM, fes.GetComm());
+
+   double glb_nrm = glb_vals[0];
+   double glb_avg = glb_vals[1];
+   glb_err = glb_vals[2];
+
+   // Normalize by the length of the boundary
+   if (std::abs(glb_nrm) > 0.0)
+   {
+      glb_err /= glb_nrm;
+      glb_avg /= glb_nrm;
+   }
+
+   // Compute l2 norm of the error in the boundary condition
+   // (negative quadrature weights may produce negative 'err')
+   glb_err = (glb_err >= 0.0) ? sqrt(glb_err) : -sqrt(-glb_err);
+
+   // Return the average value of alpha * n.Grad(x) + beta * x
+   return glb_avg;
+}

examples/ex2p.cpp

@@ -274,6 +274,13 @@ int main(int argc, char *argv[])
       pmesh->SetNodalFESpace(fespace);
    }
 
+   {
+      x.Save("ex2p.gf", 1);
+      ParGridFunction new_x(fespace, "ex2p.gf");
+      new_x -= x;
+      out << "GF difference: " << new_x.Norml1() << endl;
+   }
+
    // 16. Save in parallel the displaced mesh and the inverted solution (which
    //     gives the backward displacements to the original grid). This output
    //     can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".

examples/io_benchmark.cpp

@@ -0,0 +1,292 @@
+//                       MFEM Example 1 - Parallel Version
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex1p -m ../data/star.mesh
+//               mpirun -np 4 ex1p -m ../data/star-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/escher.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
+//
+// Device sample runs:
+//               mpirun -np 4 ex1p -pa -d cuda
+//               mpirun -np 4 ex1p -pa -d occa-cuda
+//               mpirun -np 4 ex1p -pa -d raja-omp
+//               mpirun -np 4 ex1p -pa -d ceed-cpu
+//               mpirun -np 4 ex1p -pa -d ceed-cuda
+//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mpi.h"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   // const char *mesh_file = "../data/star.mesh";
+   const char *mesh_file = "../data/square-disc.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = false;
+   int nfiles = 1;
+   // const char *out_file = "0_0.gf";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&nfiles, "-nf", "--num-files", "Number of files to write.");
+   // args.AddOption(&out_file, "-o", "--outfile",
+   //                "Name of file to write.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (pmesh->GetNodes())
+   {
+      fec = pmesh->GetNodes()->OwnFEC();
+      if (myid == 0)
+      {
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec, 1, 0);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm *b = new ParLinearForm(fespace);
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x with initial guess of zero,
+   //     which satisfies the boundary conditions.
+   ParGridFunction x(fespace);
+   x = 0.0;
+
+   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //     domain integrator.
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   // 13. Solve the linear system A X = B.
+   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
+   //     * With partial assembly, use Jacobi smoothing, for now.
+   Solver *prec = NULL;
+   if (pa)
+   {
+      if (UsesTensorBasis(*fespace))
+      {
+         prec = new OperatorJacobiSmoother(*a, ess_tdof_list);
+      }
+   }
+   else
+   {
+      prec = new HypreBoomerAMG;
+   }
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   if (prec) { cg.SetPreconditioner(*prec); }
+   cg.SetOperator(*A);
+   cg.Mult(B, X);
+   delete prec;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   std::string filename = to_string(num_procs) + "_" + to_string(nfiles) + "_";
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      x.Save(filename.c_str(), nfiles);
+      double t2 = MPI_Wtime();
+
+      double write_time = t2 - t1;
+      double average_write_time;
+      MPI_Reduce(&write_time, &average_write_time, 1,
+               MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+      if (myid == 0)
+      {
+         std::cout << "Average write time: " << average_write_time / num_procs << " for "
+                   << nfiles << " files and "  << num_procs << " ranks\n";
+      }
+   }
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      ParGridFunction temp_gf(fespace, filename.c_str());
+      double t2 = MPI_Wtime();
+
+      double read_time = t2 - t1;
+      double average_read_time;
+      MPI_Reduce(&read_time, &average_read_time, 1,
+                 MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+      if (myid == 0)
+      {
+         std::cout << "Average read time: " << average_read_time / num_procs << " for "
+                   << nfiles << " files and "  << num_procs << " ranks\n";
+      }
+   }
+
+   // 17. Free the used memory.
+   delete a;
+   delete b;
+   delete fespace;
+   if (order > 0) { delete fec; }
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

examples/makefile

@@ -22,9 +22,10 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26
+  ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p ex25p ex26p
+  ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p ex25p\
+  ex26p ex27p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -103,6 +104,10 @@ ex15-test-seq: ex15
 	@$(call mfem-test,$<,, Serial example,-e 1)
 ex15p-test-par: ex15p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-e 1)
+ex27-test-seq: ex27
+	@$(call mfem-test,$<,, Serial example,-dg)
+ex27p-test-par: ex27p
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-dg)
 # Testing: optional tests
 ifeq ($(MFEM_USE_STRUMPACK),YES)
 ex11p-test-strumpack: ex11p

examples/pa_oper.cpp

@@ -0,0 +1,906 @@
+//                                MFEM Example 9
+//
+// Compile with: make serial_nogpu
+//
+// Description:  This code solves the time-dependent advection-diffusion
+//               equation:
+//               \frac(\partial u}{\partial t}
+//                                 = \mathbf{a} \cdot \Nabla u - \nu \Nabla^2 u
+//               where a is a given advection velocity, \nu is the diffusion
+//               parameter, and u0(x) = u(0,x) is a given initial condition.
+//
+//               The demonstrates explicit time marching with H1 elements of
+//               arbitrary order. Periodic boundary conditions are used through
+//               periodic meshes. GLVis can be used for visualization of a
+//               time-evolving solution.
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+#include "mfem.hpp"
+#include "mpi.h"
+
+using namespace std;
+using namespace mfem;
+
+/** A time-dependent operator for the right-hand side of the ODE. The weak
+    form of du/dt = -a.grad(u) + nu Delta(u) is M du/dt = K u + b, where M and
+    K are the mass and advection-diffusion matrices, and b describes the flow
+    on the boundary. This can be written as a general ODE,
+    du/dt = M^{-1} (K u + b), and this class is used to evaluate the right-hand
+    side. */
+class AdvectionDiffusionEvolution : public mfem::TimeDependentOperator
+{
+public:
+   /// \param[in] M - bilinear form for mass matrix
+   /// \param[in] K - bilinear form for stiffness matrix
+   /// \param[in] b - load vector
+   AdvectionDiffusionEvolution(mfem::BilinearForm &M, mfem::BilinearForm &K,
+                               const mfem::Vector &b);
+
+   /// Perform the action of the operator: y = k = f(x, t), where k solves
+   /// Compute k = M^-1(Kx + l)
+   void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   /// Solve the implicit equation: k = f(x + dt k, t), for the unknown k at
+   /// the current time t. 
+   void ImplicitSolve(const double dt, const mfem::Vector &x,
+                      mfem::Vector &k) override;
+
+   virtual ~AdvectionDiffusionEvolution();
+
+private:
+   mfem::BilinearForm &M, &K;
+   const mfem::Vector &b;
+   /// solver for inverting mass matrix for explicit time-marching
+   std::unique_ptr<mfem::Solver> M_prec;
+   mfem::CGSolver M_solver;
+   /// solver for implicit time-marching
+   mfem::GSSmoother prec;
+   mfem::GMRESSolver linear_solver;
+   mfem::NewtonSolver newton;
+
+   mutable mfem::Vector z;
+
+   /// pointer-to-implementation idiom
+   /// Hides implementation details of this operator
+   class SystemOperator;
+   /// Operator that combines the linear spatial discretization with
+   /// the load vector into one operator used for implicit solves
+   std::unique_ptr<SystemOperator> combined_oper;
+
+   /// sets the state and dt for the combined operator
+   /// \param[in] dt - time increment
+   /// \param[in] x - the current state
+   void setOperParameters(double dt, const mfem::Vector *x);
+
+};
+
+class PAJacobianOperator : public mfem::Operator
+{
+public:
+   PAJacobianOperator(mfem::ParBilinearForm &_mass,
+                      mfem::ParBilinearForm &_stiff);
+
+   /// Compute r = J@k = M@k + dt*K@k
+   /// \param[in] k - dx/dt
+   /// \param[out] r - J@k = M@k + dt*K@k
+   void Mult(const mfem::Vector &k, mfem::Vector &r) const override;
+      
+   /// Set current dt values - needed to compute action of Jacobian.
+   void setParameters(double dt);
+
+private:
+   mfem::ParBilinearForm &mass;
+   mfem::ParBilinearForm &stiff;
+
+   double dt;
+};
+
+class ParSystemOperator : public mfem::Operator
+{
+public:
+   /// Nonlinear operator of the form that combines the mass, res, stiff,
+   /// and load elements for implicit/explicit ODE integration
+   /// \param[in] ess_bdr - array of boundaries attributes marked essential
+   /// \param[in] mass - bilinear form for mass matrix (not owned)
+   /// \param[in] res - nonlinear residual operator (not owned)
+   /// \param[in] stiff - bilinear form for stiffness matrix (not owned)
+   /// \param[in] load - load vector (not owned)
+   /// \param[in] a - used to move the spatial residual to the rhs
+   ParSystemOperator(mfem::ParBilinearForm &_mass,
+                     mfem::ParBilinearForm &_stiff);
+
+   /// Compute r = M@k + K@(x+dt*k)
+   /// (with `@` denoting matrix-vector multiplication)
+   /// \param[in] k - dx/dt 
+   /// \param[out] r - the residual
+   /// \note the signs on each operator must be accounted for elsewhere
+   void Mult(const mfem::Vector &k, mfem::Vector &r) const override;
+
+   /// Compute J = M + dt * K
+   /// \param[in] k - dx/dt 
+   mfem::Operator &GetGradient(const mfem::Vector &k) const override;
+
+   /// Set current dt and x values - needed to compute action and Jacobian.
+   void setParameters(double _dt, const mfem::Vector *_x);
+
+   ~ParSystemOperator();
+
+private:
+   mfem::ParBilinearForm &mass;
+   mfem::ParBilinearForm &stiff;
+   mutable mfem::HypreParMatrix *jacobian, *stiff_jacobian;
+
+   double dt;
+   const mfem::Vector *x;
+
+   mutable mfem::Vector work, work2;
+
+   std::unique_ptr<PAJacobianOperator> pa_jac;
+
+};
+
+/** A time-dependent operator for the right-hand side of the ODE. The weak
+    form of du/dt = -a.grad(u) + nu Delta(u) is M du/dt = K u + b, where M and
+    K are the mass and advection-diffusion matrices, and b describes the flow
+    on the boundary. This can be written as a general ODE,
+    du/dt = M^{-1} (K u + b), and this class is used to evaluate the right-hand
+    side. */
+class ParAdvectionDiffusionEvolution : public mfem::TimeDependentOperator
+{
+public:
+   /// \param[in] M - parallel bilinear form for mass matrix
+   /// \param[in] K - parallel bilinear form for stiffness matrix
+   ParAdvectionDiffusionEvolution(mfem::ParBilinearForm &M,
+                                  mfem::ParBilinearForm &K);
+
+   /// Perform the action of the operator: y = k = f(x, t), where k solves
+   /// Compute k = M^-1(Kx + l)
+   void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   /// Solve the implicit equation: k = f(x + dt k, t), for the unknown k at
+   /// the current time t. 
+   void ImplicitSolve(const double dt, const mfem::Vector &x,
+                      mfem::Vector &k) override;
+
+   virtual ~ParAdvectionDiffusionEvolution();
+
+private:
+   mfem::OperatorHandle M_;
+   mfem::ParBilinearForm &M, &K;
+   /// solver for inverting mass matrix for explicit time-marching
+   std::unique_ptr<mfem::Solver> M_prec;
+   mfem::CGSolver M_solver;
+   /// solver for implicit time-marching
+   mfem::Solver *prec;
+   mfem::GMRESSolver linear_solver;
+   mfem::NewtonSolver newton;
+
+   mfem::Vector diag;
+   mutable mfem::Vector z, work, work2;
+
+   /// pointer-to-implementation idiom
+   /// Hides implementation details of this operator
+   /// Operator that combines the linear spatial discretization with
+   /// the load vector into one operator used for implicit solves
+   std::unique_ptr<ParSystemOperator> combined_oper;
+
+   /// sets the state and dt for the combined operator
+   /// \param[in] dt - time increment
+   /// \param[in] x - the current state
+   void setOperParameters(double dt, const mfem::Vector *x);
+
+};
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &X, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &X);
+
+// Inflow boundary condition
+double inflow_function(const Vector &X, const double t);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   problem = 3;
+   const char *mesh_file = "../data/periodic-square.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order = 3;
+   const char *device_config = "cpu";
+   int ode_solver_type = 22;
+   double t_final = 3 * 2*M_PI;
+   double dt = 0.01;
+   bool glvis = false;
+   bool paraview = false;
+   int vis_steps = 5;
+
+   double nu_val = 0.001;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   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(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler,\n\t"
+                  "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&glvis, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&nu_val, "-nu", "--nu-value",
+                  "Value for \nu, the parameter that controls diffusion.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      std::cout << "Num ranks: " << num_procs << "\n";
+      args.PrintOptions(cout);
+   }
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 3. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle geometrically periodic meshes in this code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 5. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter. If the mesh is of NURBS type, we convert it
+   //    to a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+   mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 7. Define the finite element space of the given
+   //    polynomial order on the refined mesh.
+   H1_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
+
+   HYPRE_Int global_vSize = fes->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << global_vSize << endl;
+   }
+
+   // 8. Set up and assemble the bilinear and linear forms corresponding to the
+   //    CG discretization.
+   /// negative to move the diffusion terms to the right side
+   ConstantCoefficient nu(-nu_val);
+   ConstantCoefficient one(1.0);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   ParBilinearForm *m_pa = new ParBilinearForm(fes);
+   ParBilinearForm *k_pa = new ParBilinearForm(fes);
+   m_pa->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   k_pa->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+
+   /// create mass matrix
+   m_pa->AddDomainIntegrator(new MassIntegrator(one));
+   /// add advection terms to stiffness matrix
+   k_pa->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   /// add diffusion terms to stiffness matrix
+   k_pa->AddDomainIntegrator(new DiffusionIntegrator(nu));
+
+   m_pa->Assemble();
+   int skip_zeros = 0;
+   k_pa->Assemble(skip_zeros);
+   m_pa->Finalize();
+   k_pa->Finalize(skip_zeros);
+
+   ParBilinearForm *m = new ParBilinearForm(fes);
+   ParBilinearForm *k = new ParBilinearForm(fes);
+   /// create mass matrix
+   m->AddDomainIntegrator(new MassIntegrator);
+   /// add advection terms to stiffness matrix
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   /// add diffusion terms to stiffness matrix
+   k->AddDomainIntegrator(new DiffusionIntegrator(nu));
+
+   m->Assemble();
+   k->Assemble(skip_zeros);
+   m->Finalize();
+   k->Finalize(skip_zeros);
+
+
+   ParGridFunction *u = new ParGridFunction(fes);
+   u->UseDevice(true);
+   u->ProjectCoefficient(u0);
+
+
+   HypreParVector *U = u->GetTrueDofs();
+
+   ParSystemOperator pso(*m, *k);
+   ParSystemOperator pso_pa(*m_pa, *k_pa);
+
+   pso.setParameters(dt, U);
+   pso_pa.setParameters(dt, U);
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_r(U->Size());
+   double t1 = MPI_Wtime();
+   pso.Mult(*U, pso_r);
+   double t2 = MPI_Wtime();
+   double fa_mult_time = t2 - t1;
+   double average_fa_mult_time;
+   MPI_Reduce(&fa_mult_time, &average_fa_mult_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "FA Mult time: " << average_fa_mult_time / num_procs << endl;
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_pa_r(U->Size());
+   double t3 = MPI_Wtime();
+   pso_pa.Mult(*U, pso_pa_r);
+   double t4 = MPI_Wtime();
+   double pa_mult_time = t4 - t3;
+   double average_pa_mult_time;
+   MPI_Reduce(&pa_mult_time, &average_pa_mult_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "FA Mult time: " << average_pa_mult_time / num_procs << endl;
+   
+   double local_mult_speedup = (t2-t1) / (t4-t3);
+   double global_mult_speedup;
+   MPI_Reduce(&local_mult_speedup, &global_mult_speedup, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+
+   if (myid == 0)
+      std::cout << "PA mult speedup: " << global_mult_speedup / num_procs << endl;
+
+   mfem::Vector diff_r(pso_pa_r);
+   diff_r -= pso_r;
+   // std::cout << "r diff: " << diff_r.Norml2() << std::endl;
+
+   mfem::Operator &pso_jac = pso.GetGradient(*U);
+   mfem::Operator &pso_pa_jac = pso_pa.GetGradient(*U);
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_jac_r(U->Size());
+   double t5 = MPI_Wtime();
+   pso_jac.Mult(*U, pso_jac_r);
+   double t6 = MPI_Wtime();
+   double fa_jac_mult_time = t6-t5;
+   double average_fa_jac_time;
+   MPI_Reduce(&fa_jac_mult_time, &average_fa_jac_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "FA Jac Mult time: " << average_fa_jac_time / num_procs << endl;
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_pa_jac_r(U->Size());
+   double t7 = MPI_Wtime();
+   pso_pa_jac.Mult(*U, pso_pa_jac_r);
+   double t8 = MPI_Wtime();
+   double pa_jac_mult_time = t8-t7;
+   double average_pa_jac_time;
+   MPI_Reduce(&pa_jac_mult_time, &average_pa_jac_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "PA Jac Mult time: " << average_pa_jac_time / num_procs << endl;
+
+   double local_jac_speedup = (t6-t5) / (t8-t7);
+   double global_jac_speedup;
+   MPI_Reduce(&local_jac_speedup, &global_jac_speedup, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "PA Jac mult speedup: " << global_jac_speedup / num_procs << endl;
+
+   // 13. Free the used memory.
+   delete U;
+   delete u;
+   delete k;
+   delete m;
+   delete fes;
+   delete pmesh;
+
+   MPI_Finalize();
+   return 0;
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 3:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 0:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x, const double t)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}
+
+class AdvectionDiffusionEvolution::SystemOperator : public mfem::Operator
+{
+public:
+   /// Nonlinear operator of the form that combines the mass, res, stiff,
+   /// and load elements for implicit/explicit ODE integration
+   /// \param[in] mass - bilinear form for mass matrix (not owned)
+   /// \param[in] res - nonlinear residual operator (not owned)
+   /// \param[in] stiff - bilinear form for stiffness matrix (not owned)
+   /// \param[in] load - load vector (not owned)
+   /// \param[in] a - used to move the spatial residual to the rhs
+   SystemOperator(BilinearForm &_mass, BilinearForm &_stiff,
+                  const mfem::Vector &b)
+      : Operator(_mass.Height()), mass(_mass), stiff(_stiff),
+        load(b), Jacobian(NULL), dt(0.0), x(NULL), work(height)
+      { }
+
+   /// Compute r = M@k + K@(x+dt*k) + l
+   /// (with `@` denoting matrix-vector multiplication)
+   /// \param[in] k - dx/dt 
+   /// \param[out] r - the residual
+   /// \note the signs on each operator must be accounted for elsewhere
+   void Mult(const mfem::Vector &k, mfem::Vector &r) const override
+   {
+      /// work = x+dt*k = x+dt*dx/dt = x+dx
+      add(1.0, *x, dt, k, work);
+      r = 0.0;
+      stiff.AddMult(work, r);
+      r += load;
+      mass.AddMult(k, r, -1.0);
+   }
+
+   /// Compute J = M + dt * K
+   /// \param[in] k - dx/dt 
+   mfem::Operator &GetGradient(const mfem::Vector &k) const override
+   {
+      delete Jacobian;
+      Jacobian = Add(-1.0, mass.SpMat(), dt, stiff.SpMat());
+      return *Jacobian;
+   }
+
+   /// Set current dt and x values - needed to compute action and Jacobian.
+   void setParameters(double _dt, const mfem::Vector *_x)
+   {
+      dt = _dt;
+      x = _x;
+   };
+
+   ~SystemOperator() {delete Jacobian;};
+
+private:
+   BilinearForm &mass;
+   BilinearForm &stiff;
+   const mfem::Vector &load;
+   mutable mfem::SparseMatrix *Jacobian;
+
+   double dt;
+   const mfem::Vector *x;
+
+   mutable mfem::Vector work, work2;
+
+};
+
+AdvectionDiffusionEvolution::AdvectionDiffusionEvolution(
+      BilinearForm &_M, BilinearForm &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()), M(_M), K(_K), b(_b),
+     z(_M.Height())
+{
+   bool pa = M.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   Array<int> ess_tdof_list;
+   if (pa)
+   {
+      M_prec.reset(new OperatorJacobiSmoother(M, ess_tdof_list));
+      M_solver.SetOperator(M);
+   }
+   else
+   {
+      M_prec.reset(new DSmoother(M.SpMat()));
+      M_solver.SetOperator(M.SpMat());
+   }
+
+   combined_oper.reset(new SystemOperator(_M, _K, _b));
+
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+
+   linear_solver.iterative_mode = true;
+   linear_solver.SetRelTol(1e-12);
+   linear_solver.SetAbsTol(0.0);
+   linear_solver.SetMaxIter(100);
+   linear_solver.SetPrintLevel(0);
+   linear_solver.SetPreconditioner(prec);
+
+   newton.iterative_mode = false;
+   newton.SetRelTol(1e-9);
+   newton.SetAbsTol(0.0);
+   newton.SetMaxIter(100);
+   newton.SetPrintLevel(-1);
+   newton.SetSolver(linear_solver);
+   newton.SetOperator(*combined_oper);
+}
+
+void AdvectionDiffusionEvolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void AdvectionDiffusionEvolution::ImplicitSolve(const double dt,
+                                                const Vector &x,
+                                                Vector &k)
+{
+   setOperParameters(dt, &x);
+   Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
+   newton.Mult(zero, k);
+   MFEM_VERIFY(newton.GetConverged(), "Newton solver did not converge!");
+}
+
+void AdvectionDiffusionEvolution::setOperParameters(double dt,
+                                                    const mfem::Vector *x)
+{
+   combined_oper->setParameters(dt, x);
+}
+
+AdvectionDiffusionEvolution::~AdvectionDiffusionEvolution() {}
+
+
+PAJacobianOperator::PAJacobianOperator(ParBilinearForm &_mass, ParBilinearForm &_stiff)
+ : Operator(_mass.ParFESpace()->GetTrueVSize()), mass(_mass), stiff(_stiff),
+   dt(0.0) { }
+
+
+void PAJacobianOperator::Mult(const mfem::Vector &k, mfem::Vector &r) const
+{
+   r.UseDevice(true);
+   r = 0.0;
+   stiff.TrueAddMult(k, r, dt);
+   mass.TrueAddMult(k, r, -1.0);
+}
+   
+void PAJacobianOperator::setParameters(const double _dt)
+{
+   dt = _dt;
+};
+
+ParSystemOperator::ParSystemOperator(ParBilinearForm &_mass, ParBilinearForm &_stiff)
+   : Operator(_mass.ParFESpace()->GetTrueVSize()), mass(_mass), stiff(_stiff),
+   jacobian(NULL), stiff_jacobian(NULL), dt(0.0), x(NULL),
+   work(height)
+{
+   pa_jac.reset(new PAJacobianOperator(mass, stiff));
+}
+
+/// Compute r = M@k + K@(x+dt*k)
+/// (with `@` denoting matrix-vector multiplication)
+/// \param[in] k - dx/dt 
+/// \param[out] r - the residual
+/// \note the signs on each operator must be accounted for elsewhere
+void ParSystemOperator::Mult(const mfem::Vector &k, mfem::Vector &r) const
+{
+   r = 0.0;
+   work.UseDevice(true);
+   work = 0.0;
+   /// work = x+dt*k = x+dt*dx/dt = x+dx
+   if (x)
+   {
+      add(1.0, *x, dt, k, work);
+   }
+
+   stiff.TrueAddMult(work, r);
+   mass.TrueAddMult(k, r, -1.0);
+}
+
+/// Compute J = M + dt * K
+/// \param[in] k - dx/dt 
+mfem::Operator &ParSystemOperator::GetGradient(const mfem::Vector &k) const
+{
+   bool mass_pa = mass.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   bool stiff_pa = stiff.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   
+   if (mass_pa && stiff_pa)
+   {
+      return *pa_jac.get();
+   }
+   else
+   {
+      delete stiff_jacobian;
+      delete jacobian;
+      jacobian = mass.ParallelAssemble();
+      *jacobian *= -1.0; //alpha;
+      stiff_jacobian = stiff.ParallelAssemble();
+      jacobian->Add(dt, *stiff_jacobian);
+      return *jacobian;
+   }
+}
+
+/// Set current dt and x values - needed to compute action and Jacobian.
+void ParSystemOperator::setParameters(const double _dt, const mfem::Vector *_x)
+{
+   dt = _dt;
+   x = _x;
+   pa_jac->setParameters(_dt);
+};
+
+ParSystemOperator::~ParSystemOperator()
+{
+   delete jacobian;
+   delete stiff_jacobian;
+};
+
+ParAdvectionDiffusionEvolution::ParAdvectionDiffusionEvolution(
+      ParBilinearForm &_M, ParBilinearForm &_K)
+   : TimeDependentOperator(_M.ParFESpace()->GetTrueVSize()), M(_M), K(_K), z(_M.Height())
+{
+   bool mass_pa = M.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   bool stiff_pa = K.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+
+   Array<int> ess_tdof_list;
+   M_solver = CGSolver(MPI_COMM_WORLD);
+   if (mass_pa)
+   {
+      M_prec.reset(new OperatorJacobiSmoother(M, ess_tdof_list));
+      M_solver.SetOperator(M);
+   }
+   else
+   {
+      M_.Reset(_M.ParallelAssemble(), true);
+
+      // M_prec.reset(new HypreSmoother());
+      // M_solver.SetOperator(M.As<HypreParMatrix>());
+      HypreParMatrix &M_mat = *M_.As<HypreParMatrix>();
+      // HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
+      M_prec.reset(new HypreSmoother(M_mat, HypreSmoother::Jacobi));
+   }
+
+   combined_oper.reset(new ParSystemOperator(_M, _K));
+
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+
+   if (mass_pa && stiff_pa)
+   {
+      diag.UseDevice(true);
+      diag.SetSize(M.ParFESpace()->GetTrueVSize());
+      diag = 0.0;
+      work.UseDevice(true);
+      work2.UseDevice(true);
+      work.SetSize(M.ParFESpace()->GetTrueVSize());
+      work2.SetSize(M.ParFESpace()->GetTrueVSize());
+      work = 0.0;
+      work2 = 0.0;
+      M.AssembleDiagonal(work);
+
+      ParBilinearForm k(M.ParFESpace());
+      ConstantCoefficient nu(-0.01);
+      k.AddDomainIntegrator(new mfem::DiffusionIntegrator(nu));
+      k.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      k.Assemble(0);
+      k.Finalize(0);
+      k.AssembleDiagonal(work2);
+
+      double dt = 0.1;
+      add(-1.0, work, dt, work2, diag);
+
+      prec = new OperatorChebyshevSmoother(combined_oper.get(), diag,
+                                           ess_tdof_list, 5,
+                                           M.ParFESpace()->GetComm());
+   }
+   else
+   {
+      prec = new HypreSmoother();
+   }
+   
+   linear_solver = GMRESSolver(MPI_COMM_WORLD);
+   linear_solver.iterative_mode = true;
+   linear_solver.SetRelTol(1e-12);
+   linear_solver.SetAbsTol(0.0);
+   linear_solver.SetMaxIter(2000);
+   linear_solver.SetPrintLevel(0);
+   linear_solver.SetPreconditioner(*prec);
+   linear_solver.SetKDim(2000);
+
+   newton.iterative_mode = true;
+   newton.SetRelTol(1e-9);
+   newton.SetAbsTol(0.0);
+   newton.SetMaxIter(10);
+   newton.SetPrintLevel(-1);
+   newton.SetSolver(linear_solver);
+   newton.SetOperator(*combined_oper);
+}
+
+void ParAdvectionDiffusionEvolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   M_solver.Mult(z, y);
+}
+
+void ParAdvectionDiffusionEvolution::ImplicitSolve(const double dt,
+                                                   const Vector &x,
+                                                   Vector &k)
+{
+   setOperParameters(dt, &x);
+   Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
+   newton.Mult(zero, k);
+   MFEM_VERIFY(newton.GetConverged(), "Newton solver did not converge!");
+}
+
+void ParAdvectionDiffusionEvolution::setOperParameters(const double dt,
+                                                       const mfem::Vector *x)
+{
+   combined_oper->setParameters(dt, x);
+}
+
+ParAdvectionDiffusionEvolution::~ParAdvectionDiffusionEvolution() {delete prec;}

fem/bilinearform.cpp

@@ -467,8 +467,17 @@ void BilinearForm::Assemble(int skip_zeros)
          const FiniteElement &be = *fes->GetBE(i);
          fes -> GetBdrElementVDofs (i, vdofs);
          eltrans = fes -> GetBdrElementTransformation (i);
-         bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
-         for (int k = 1; k < bbfi.Size(); k++)
+         int k = 0;
+         for (; k < bbfi.Size(); k++)
+         {
+            if (bbfi_marker[k] &&
+                (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
+
+            bbfi[k]->AssembleElementMatrix(be, *eltrans, elmat);
+            k++;
+            break;
+         }
+         for (; k < bbfi.Size(); k++)
          {
             if (bbfi_marker[k] &&
                 (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }

fem/bilininteg.hpp

@@ -453,6 +453,16 @@ public:
                                        ElementTransformation &Trans,
                                        DenseMatrix &elmat);
 
+   /// Support for use in BilinearForm. Can be used only when appropriate.
+   /** Appropriate use cases are classes derived from
+       MixedScalarVectorIntegrator where the trial and test spaces can be the
+       same. Examples of such classes are: MixedVectorDivergenceIntegrator,
+       MixedScalarWeakDivergenceIntegrator, etc. */
+   virtual void AssembleElementMatrix(const FiniteElement &fe,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat)
+   { AssembleElementMatrix2(fe, fe, Trans, elmat); }
+
 protected:
 
    MixedScalarVectorIntegrator(VectorCoefficient &vq, bool _transpose = false,

fem/coefficient.cpp

@@ -174,7 +174,24 @@ void VectorGridFunctionCoefficient::SetGridFunction(GridFunction *gf)
 void VectorGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
                                          const IntegrationPoint &ip)
 {
-   GridFunc->GetVectorValue(T.ElementNo, ip, V);
+   Mesh *mesh = GridFunc->FESpace()->GetMesh();
+   if (mesh->Dimension() == T.GetDimension())
+   {
+      GridFunc->GetVectorValue(T.ElementNo, ip, V);
+   }
+   else // Assuming T is a boundary element transformation
+   {
+      int el_id, el_info;
+      mesh->GetBdrElementAdjacentElement(T.ElementNo, el_id, el_info);
+      IntegrationPointTransformation loc_T;
+      mesh->GetLocalFaceTransformation(mesh->GetBdrElementType(T.ElementNo),
+                                       mesh->GetElementType(el_id),
+                                       loc_T.Transf,
+                                       el_info);
+      IntegrationPoint eip;
+      loc_T.Transform(ip, eip);
+      GridFunc->GetVectorValue(el_id, eip, V);
+   }
 }
 
 void VectorGridFunctionCoefficient::Eval(

fem/lininteg.cpp

@@ -93,6 +93,37 @@ void BoundaryLFIntegrator::AssembleRHSElementVect(
    }
 }
 
+void BoundaryLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
+{
+   int dof = el.GetDof();
+
+   shape.SetSize(dof);        // vector of size dof
+   elvect.SetSize(dof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int intorder = oa * el.GetOrder() + ob;    // <------ user control
+      ir = &IntRules.Get(Tr.FaceGeom, intorder); // of integration order
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      IntegrationPoint eip;
+      Tr.Loc1.Transform(ip, eip);
+
+      Tr.Face->SetIntPoint (&ip);
+      double val = Tr.Face->Weight() * ip.weight * Q.Eval(*Tr.Face, ip);
+
+      el.CalcShape(eip, shape);
+
+      add(elvect, val, shape, elvect);
+   }
+}
+
 void BoundaryNormalLFIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {

fem/lininteg.hpp

@@ -126,7 +126,8 @@ class BoundaryLFIntegrator : public LinearFormIntegrator
    Coefficient &Q;
    int oa, ob;
 public:
-   /// Constructs a boundary integrator with a given Coefficient QG
+   /** @brief Constructs a boundary integrator with a given Coefficient @a QG.
+       Integration order will be @a a * basis_order + @a b. */
    BoundaryLFIntegrator(Coefficient &QG, int a = 1, int b = 1)
       : Q(QG), oa(a), ob(b) { }
 
@@ -135,8 +136,9 @@ public:
    virtual void AssembleRHSElementVect(const FiniteElement &el,
                                        ElementTransformation &Tr,
                                        Vector &elvect);
-
-   using LinearFormIntegrator::AssembleRHSElementVect;
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       FaceElementTransformations &Tr,
+                                       Vector &elvect);
 };
 
 /// Class for boundary integration \f$ L(v) = (g \cdot n, v) \f$

fem/pbilinearform.cpp

@@ -283,7 +283,14 @@ const
    }
 
    X.Distribute(&x);
-   mat->Mult(X, Y);
+   if (ext)
+   {
+      ext->Mult(X, Y);
+   }
+   else
+   {
+      mat->Mult(X, Y);
+   }
    pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Y, 1.0, y);
 }
 

fem/pfespace.hpp

@@ -240,7 +240,7 @@ public:
    int GetNRanks() const { return NRanks; }
    int GetMyRank() const { return MyRank; }
 
-   inline ParMesh *GetParMesh() { return pmesh; }
+   inline ParMesh *GetParMesh() const { return pmesh; }
 
    int GetDofSign(int i)
    { return NURBSext || Nonconforming() ? 1 : ldof_sign[VDofToDof(i)]; }

fem/pgridfunc.cpp

@@ -16,6 +16,7 @@
 #include "fem.hpp"
 #include <iostream>
 #include <limits>
+#include <string>
 #include "../general/forall.hpp"
 using namespace std;
 
@@ -78,6 +79,229 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
    fes = pfes;
 }
 
+ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf,
+                                 const char *_filename)
+   : GridFunction(pf), pfes(pf)
+{
+   MPI_Comm fes_comm;
+   int fes_rank, n_fes_ranks;
+   fes_comm = pfes->GetComm();
+   MPI_Comm_size(fes_comm, &n_fes_ranks);
+   MPI_Comm_rank(fes_comm, &fes_rank);
+
+   std::string filename(_filename);
+   std::string file_prefix;
+   std::string file_ext;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         file_prefix = (filename.substr(0, i));
+         file_ext = (filename.substr(i, filename.length() - i));
+      }
+   }
+
+   int nfiles = 1;
+   if (fes_rank == 0)
+   {
+      int n_rfes_ranks;
+      int tmp[2];
+      std::string mpi_filename;
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         mpi_filename = file_prefix + to_string(0) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + to_string(0);
+      }
+      
+      MPI_File fh;
+      MPI_File_open(MPI_COMM_SELF, mpi_filename.c_str(), MPI_MODE_RDONLY,
+                                                      MPI_INFO_NULL, &fh);
+      MPI_File_read_at(fh, 0, tmp, 2, MPI_INT, MPI_STATUS_IGNORE);
+      MPI_File_close(&fh);
+      
+      n_rfes_ranks = tmp[0];
+      nfiles = tmp[1];
+
+      MFEM_ASSERT(n_fes_ranks == n_rfes_ranks,
+         "ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf,"
+               " const char *_filename):\n"
+               "\tThe number of MPI ranks used to save the GridFunction is\n"
+               "\tnot the same as the number used to load it!");
+   }
+   MPI_Bcast(&nfiles, 1, MPI_INT, 0, fes_comm);
+
+   int color = fes_rank * nfiles / n_fes_ranks;
+   
+   MPI_Comm file_comm;
+   MPI_Comm_split(fes_comm, color, fes_rank, &file_comm);
+
+   int file_rank, n_file_ranks;
+   MPI_Comm_size(file_comm, &n_file_ranks);
+   MPI_Comm_rank(file_comm, &file_rank);
+
+   std::string mpi_filename;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos) {
+         mpi_filename = file_prefix + std::to_string(color) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + std::to_string(color);
+      }
+   }
+
+   MPI_File fh;
+   MPI_File_open(file_comm, mpi_filename.c_str(), MPI_MODE_RDONLY,
+                                                  MPI_INFO_NULL, &fh);
+
+   int *dof_counts = new int[5*n_file_ranks];
+   int **nv = new int*[n_file_ranks];
+   int **nvdofs = new int*[n_file_ranks];
+   int **nedofs = new int*[n_file_ranks];
+   int **nfdofs = new int*[n_file_ranks];
+   int **nrdofs = new int*[n_file_ranks];
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      nv[i]     = &dof_counts[i*5+0];
+      nvdofs[i] = &dof_counts[i*5+1];
+      nedofs[i] = &dof_counts[i*5+2];
+      nfdofs[i] = &dof_counts[i*5+3];
+      nrdofs[i] = &dof_counts[i*5+4];
+   }
+
+   *nv[file_rank] = pfes->GetVSize();
+   *nvdofs[file_rank] = pfes->GetNVDofs();
+   *nedofs[file_rank] = pfes->GetNEDofs();
+   *nfdofs[file_rank] = pfes->GetNFDofs();
+
+   int vdim = pfes->GetVDim();
+   *nrdofs[file_rank] = *nv[file_rank] / vdim - *nvdofs[file_rank] -
+                                       *nedofs[file_rank] - *nfdofs[file_rank];
+
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &dof_counts[0], 5,
+                 MPI_INT, file_comm);
+
+   double *data_ = HostWrite();
+
+   MPI_Offset header_offset = 0;
+   header_offset += 2 * sizeof(int);
+   MPI_Offset v_offset, e_offset, f_offset, r_offset;
+
+   int total_vdofs = 0, total_edofs = 0, total_fdofs = 0, total_rdofs = 0;
+   int total_scalar_dofs = 0;
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      total_vdofs += *nvdofs[i];
+      total_edofs += *nedofs[i];
+      total_fdofs += *nfdofs[i];
+      total_rdofs += *nrdofs[i];
+      total_scalar_dofs += *nv[i];
+   }
+
+   total_scalar_dofs /= vdim;
+
+   if (pfes->GetOrdering() == Ordering::byNODES)
+   {
+      for (int d = 0; d < vdim; ++d)
+      {
+         int v_data_offset = 0 + *nv[file_rank] * d / vdim ;
+         int e_data_offset = v_data_offset + *nvdofs[file_rank];
+         int f_data_offset = e_data_offset + *nedofs[file_rank];
+         int r_data_offset = f_data_offset + *nfdofs[file_rank];
+
+         v_offset = header_offset;
+         e_offset = header_offset;
+         f_offset = header_offset;
+         r_offset = header_offset;
+         
+         v_offset += total_scalar_dofs * d * sizeof(double);
+         e_offset += (total_vdofs + total_scalar_dofs * d) * sizeof(double);
+         f_offset += (total_vdofs + total_edofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+         r_offset += (total_vdofs + total_edofs + total_fdofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+
+
+         for (int i = 0; i < file_rank; ++i)
+         {
+            v_offset += *nvdofs[i] * sizeof(double);
+            e_offset += *nedofs[i] * sizeof(double);
+            f_offset += *nfdofs[i] * sizeof(double);
+            r_offset += *nrdofs[i] * sizeof(double);
+         }
+
+         MPI_File_read_at_all(fh, v_offset, &data_[v_data_offset],
+                              *nvdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, e_offset, &data_[e_data_offset],
+                              *nedofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, f_offset, &data_[f_data_offset],
+                              *nfdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, r_offset, &data_[r_data_offset],
+                              *nrdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+      }
+   }
+   else
+   {
+      v_offset = header_offset;
+      e_offset = v_offset + total_vdofs * vdim * sizeof(double);
+      f_offset = e_offset + total_edofs * vdim * sizeof(double);
+      r_offset = f_offset + total_fdofs * vdim * sizeof(double);
+
+      for (int i = 0; i < file_rank; ++i)
+      {
+         v_offset += *nvdofs[i] * sizeof(double) * vdim;
+         e_offset += *nedofs[i] * sizeof(double) * vdim;
+         f_offset += *nfdofs[i] * sizeof(double) * vdim;
+         r_offset += *nrdofs[i] * sizeof(double) * vdim;
+      }
+
+      int v_data_offset = 0;
+      int e_data_offset = v_data_offset + *nvdofs[file_rank] * vdim;
+      int f_data_offset = e_data_offset + *nedofs[file_rank] * vdim;
+      int r_data_offset = f_data_offset + *nfdofs[file_rank] * vdim;
+      
+      MPI_File_read_at_all(fh, v_offset, &data_[v_data_offset],
+                           *nvdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, e_offset, &data_[e_data_offset],
+                           *nedofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, f_offset, &data_[f_data_offset],
+                           *nfdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, r_offset, &data_[r_data_offset],
+                           *nrdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+   }
+
+   MPI_File_close(&fh);
+   MPI_Comm_free(&file_comm);
+
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   delete[] dof_counts;
+   delete[] nv;
+   delete[] nvdofs;
+   delete[] nedofs;
+   delete[] nfdofs;
+   delete[] nrdofs;
+}
+
+
 void ParGridFunction::Update()
 {
    face_nbr_data.Destroy();
@@ -518,6 +742,201 @@ void ParGridFunction::Save(adios2stream &out,
 }
 #endif
 
+void ParGridFunction::Save(const char *_filename, const int nfiles)
+{
+   MPI_Comm fes_comm;
+   int fes_rank, n_fes_ranks;
+   fes_comm = pfes->GetComm();
+
+   MPI_Comm_size(fes_comm, &n_fes_ranks);
+   MPI_Comm_rank(fes_comm, &fes_rank);
+
+   int color = fes_rank * nfiles / n_fes_ranks;
+   
+   MPI_Comm file_comm;
+   MPI_Comm_split(fes_comm, color, fes_rank, &file_comm);
+
+   int file_rank, n_file_ranks;
+   MPI_Comm_size(file_comm, &n_file_ranks);
+   MPI_Comm_rank(file_comm, &file_rank);
+
+   std::string filename(_filename);
+   std::string file_prefix;
+   std::string file_ext;
+   std::string mpi_filename;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         file_prefix = (filename.substr(0, i));
+         file_ext = (filename.substr(i, filename.length() - i));
+         mpi_filename = file_prefix + std::to_string(color) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + std::to_string(color);
+      }
+   }
+
+   MPI_File fh;
+   MPI_File_open(file_comm, mpi_filename.c_str(), MPI_MODE_CREATE |
+                                                  MPI_MODE_WRONLY,
+                                                  MPI_INFO_NULL, &fh);
+
+   int *dof_counts = new int[5*n_file_ranks];
+   int **nv = new int*[n_file_ranks];
+   int **nvdofs = new int*[n_file_ranks];
+   int **nedofs = new int*[n_file_ranks];
+   int **nfdofs = new int*[n_file_ranks];
+   int **nrdofs = new int*[n_file_ranks];
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      nv[i]     = &dof_counts[i*5+0];
+      nvdofs[i] = &dof_counts[i*5+1];
+      nedofs[i] = &dof_counts[i*5+2];
+      nfdofs[i] = &dof_counts[i*5+3];
+      nrdofs[i] = &dof_counts[i*5+4];
+   }
+
+   *nv[file_rank] = pfes->GetVSize();
+   *nvdofs[file_rank] = pfes->GetNVDofs();
+   *nedofs[file_rank] = pfes->GetNEDofs();
+   *nfdofs[file_rank] = pfes->GetNFDofs();
+
+   int vdim = pfes->GetVDim();
+   *nrdofs[file_rank] = *nv[file_rank] / vdim - *nvdofs[file_rank] -
+                                       *nedofs[file_rank] - *nfdofs[file_rank];
+
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &dof_counts[0], 5,
+                 MPI_INT, file_comm);
+
+   double *data_  = const_cast<double*>(HostRead());
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   MPI_Offset header_offset = 0;
+
+   if (file_rank == 0)
+   {
+      int tmp[] = {n_fes_ranks, nfiles};
+      MPI_File_write_at(fh, header_offset, &tmp, 2, MPI_INT,
+                        MPI_STATUS_IGNORE);
+   }
+
+   header_offset += 2 * sizeof(int);
+
+   MPI_Offset v_offset, e_offset, f_offset, r_offset;
+
+   int total_vdofs = 0, total_edofs = 0, total_fdofs = 0, total_rdofs = 0;
+   int total_scalar_dofs = 0;
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      total_vdofs += *nvdofs[i];
+      total_edofs += *nedofs[i];
+      total_fdofs += *nfdofs[i];
+      total_rdofs += *nrdofs[i];
+      total_scalar_dofs += *nv[i];
+   }
+
+   total_scalar_dofs /= vdim;
+
+   if (pfes->GetOrdering() == Ordering::byNODES)
+   {
+      for (int d = 0; d < vdim; ++d)
+      {
+         int v_data_offset = 0 + *nv[file_rank] * d / vdim ;
+         int e_data_offset = v_data_offset + *nvdofs[file_rank];
+         int f_data_offset = e_data_offset + *nedofs[file_rank];
+         int r_data_offset = f_data_offset + *nfdofs[file_rank];
+
+         v_offset = header_offset;
+         e_offset = header_offset;
+         f_offset = header_offset;
+         r_offset = header_offset;
+         
+         v_offset += total_scalar_dofs * d * sizeof(double);
+         e_offset += (total_vdofs + total_scalar_dofs * d) * sizeof(double);
+         f_offset += (total_vdofs + total_edofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+         r_offset += (total_vdofs + total_edofs + total_fdofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+
+         for (int i = 0; i < file_rank; ++i)
+         {
+            v_offset += *nvdofs[i] * sizeof(double);
+            e_offset += *nedofs[i] * sizeof(double);
+            f_offset += *nfdofs[i] * sizeof(double);
+            r_offset += *nrdofs[i] * sizeof(double);
+         }
+
+         MPI_File_write_at_all(fh, v_offset, &data_[v_data_offset],
+                               *nvdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, e_offset, &data_[e_data_offset],
+                               *nedofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, f_offset, &data_[f_data_offset],
+                               *nfdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, r_offset, &data_[r_data_offset],
+                               *nrdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+      }
+   }
+   else
+   {
+      v_offset = header_offset;
+      e_offset = v_offset + total_vdofs * vdim * sizeof(double);
+      f_offset = e_offset + total_edofs * vdim * sizeof(double);
+      r_offset = f_offset + total_fdofs * vdim * sizeof(double);
+
+      for (int i = 0; i < file_rank; ++i)
+      {
+         v_offset += *nvdofs[i] * sizeof(double) * vdim;
+         e_offset += *nedofs[i] * sizeof(double) * vdim;
+         f_offset += *nfdofs[i] * sizeof(double) * vdim;
+         r_offset += *nrdofs[i] * sizeof(double) * vdim;
+      }
+
+      int v_data_offset = 0;
+      int e_data_offset = v_data_offset + *nvdofs[file_rank] * vdim;
+      int f_data_offset = e_data_offset + *nedofs[file_rank] * vdim;
+      int r_data_offset = f_data_offset + *nfdofs[file_rank] * vdim;
+      
+      MPI_File_write_at_all(fh, v_offset, &data_[v_data_offset],
+                            *nvdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, e_offset, &data_[e_data_offset],
+                            *nedofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, f_offset, &data_[f_data_offset],
+                            *nfdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, r_offset, &data_[r_data_offset],
+                            *nrdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+   }
+
+   MPI_File_close(&fh);
+   MPI_Comm_free(&file_comm);
+   
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   delete[] dof_counts;
+   delete[] nv;
+   delete[] nvdofs;
+   delete[] nedofs;
+   delete[] nfdofs;
+   delete[] nrdofs;
+}
+
 void ParGridFunction::SaveAsOne(std::ostream &out)
 {
    int i, p;

fem/pgridfunc.hpp

@@ -83,6 +83,13 @@ public:
        constructed. The new ParGridFunction assumes ownership of both. */
    ParGridFunction(ParMesh *pmesh, std::istream &input);
 
+   /// Construct a ParGridFunction by loading a ParGridFunction saved using
+   /// ParGridFunction::Save(char *filename, int nfiles).
+   /** The parallel space @a *pf and the space used by the GridFunction saved
+       in @a *filename should match. The number of ranks used when loading the
+       ParGridFunction must be the same as when it was saved. */
+   ParGridFunction(ParFiniteElementSpace *pf, const char *filename);
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use ParFiniteElementSpace%s
        that have the same size.
@@ -324,6 +331,20 @@ public:
       const adios2stream::data_type type = adios2stream::data_type::point_data) const;
 #endif
 
+   /** Save the local grid functions to n number of files, where each file will
+       contain the grid functions from potentially multiple ranks. This is
+       similar to the syncIO approach from "Fu, Jing, et al. 'Scalable parallel
+       I/O alternatives for massively parallel partitioned solver systems.'
+       2010 IEEE International Symposium on Parallel & Distributed Processing,
+       Workshops and Phd Forum (IPDPSW). IEEE, 2010." 
+       @param[in] filename - filename for output files with extension
+       @param[in] nfiles - number of files to write using MPI-IO
+       @note - takes into account the signs of the local dofs.
+       @note - writes a binary file without the FESpace header; the saved file
+       should only be loaded by the accompanying constructor:
+       ParGridFunction(ParFiniteElementSpace *pf, const char *filename) */
+   void Save(const char *filename, const int nfiles = 1);
+
    /// Merge the local grid functions
    void SaveAsOne(std::ostream &out = mfem::out);
 

makefile

@@ -117,7 +117,7 @@ EXAMPLE_SUBDIRS = sundials petsc pumi hiop ginkgo
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
-MINIAPP_SUBDIRS = common electromagnetics meshing performance tools toys nurbs gslib
+MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools toys nurbs gslib
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics tools toys)

miniapps/CMakeLists.txt

@@ -16,6 +16,7 @@ set(MINIAPPS_COMMON_SOURCES "")
 set(MINIAPPS_COMMON_HEADERS "")
 add_subdirectory(common)
 add_subdirectory(electromagnetics)
+add_subdirectory(navier)
 add_subdirectory(meshing)
 add_subdirectory(performance)
 add_subdirectory(tools)

miniapps/navier/CMakeLists.txt

@@ -0,0 +1,45 @@
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+if (MFEM_USE_MPI)
+  list(APPEND NAVIER_COMMON_FILES
+      navier_solver.cpp
+      navier_solver.hpp
+      ortho_solver.cpp
+      ortho_solver.hpp)
+
+  convert_filenames_to_full_paths(NAVIER_COMMON_FILES)
+
+  add_mfem_miniapp(navier_mms
+    MAIN navier_mms.cpp
+    ${NAVIER_COMMON_FILES}
+    LIBRARIES mfem)
+
+  add_mfem_miniapp(navier_kovasznay
+    MAIN navier_kovasznay.cpp
+    ${NAVIER_COMMON_FILES}
+    LIBRARIES mfem)
+
+  add_mfem_miniapp(navier_tgv
+    MAIN navier_tgv.cpp
+    ${NAVIER_COMMON_FILES}
+    LIBRARIES mfem)
+
+  add_mfem_miniapp(navier_shear
+    MAIN navier_shear.cpp
+    ${NAVIER_COMMON_FILES}
+    LIBRARIES mfem)
+
+  add_mfem_miniapp(navier_3dfoc
+    MAIN navier_3dfoc.cpp
+    ${NAVIER_COMMON_FILES}
+    LIBRARIES mfem)
+endif ()

miniapps/navier/makefile

@@ -0,0 +1,77 @@
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+# SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/navier/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+NAVIER_COMMON_SRC = navier_solver.cpp ortho_solver.cpp
+NAVIER_COMMON_OBJ = $(NAVIER_COMMON_SRC:.cpp=.o)
+
+PAR_MINIAPPS = navier_mms navier_kovasznay navier_tgv navier_shear navier_3dfoc
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS =
+else
+   MINIAPPS = $(PAR_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+%: %.cpp $(NAVIER_COMMON_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(NAVIER_COMMON_OBJ) $(MFEM_LIBS)
+
+%.o: %.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+all: $(MINIAPPS)
+
+PAR_TEST_EXECS = navier_mms navier_kovasznay navier_tgv
+MFEM_TESTS = TEST_EXECS
+include $(MFEM_TEST_MK)
+
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+
+navier_mms-test-par: navier_mms
+	@$(call mfem-test,$<, $(RUN_MPI), Navier,-cr)
+
+navier_kovasznay-test-par: navier_kovasznay
+	@$(call mfem-test,$<, $(RUN_MPI), Navier,-cr)
+
+navier_tgv-test-par: navier_tgv
+	@$(call mfem-test,$<, $(RUN_MPI), Navier,-cr)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf tgv_out*.txt *_output

miniapps/navier/navier_3dfoc.cpp

@@ -0,0 +1,137 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// 3d flow over a cylinder benchmark example
+
+#include "navier_solver.hpp"
+#include <fstream>
+
+using namespace mfem;
+using namespace navier;
+
+struct s_NavierContext
+{
+   int order = 4;
+   double kin_vis = 0.001;
+   double t_final = 8.0;
+   double dt = 1e-3;
+} ctx;
+
+void vel(const Vector &x, double t, Vector &u)
+{
+   double xi = x(0);
+   double yi = x(1);
+   double zi = x(2);
+
+   double U = 2.25;
+
+   if (xi <= 1e-8)
+   {
+      u(0) = 16.0 * U * yi * zi * sin(M_PI * t / 8.0) * (0.41 - yi)
+             * (0.41 - zi) / pow(0.41, 4.0);
+   }
+   else
+   {
+      u(0) = 0.0;
+   }
+   u(1) = 0.0;
+   u(2) = 0.0;
+}
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   int serial_refinements = 0;
+
+   Mesh *mesh = new Mesh("box-cylinder.mesh");
+
+   for (int i = 0; i < serial_refinements; ++i)
+   {
+      mesh->UniformRefinement();
+   }
+
+   if (mpi.Root())
+   {
+      std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
+   }
+
+   auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Create the flow solver.
+   NavierSolver flowsolver(pmesh, ctx.order, ctx.kin_vis);
+   flowsolver.EnablePA(true);
+
+   // Set the initial condition.
+   ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
+   VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
+   u_ic->ProjectCoefficient(u_excoeff);
+
+   // Add Dirichlet boundary conditions to velocity space restricted to
+   // selected attributes on the mesh.
+   Array<int> attr(pmesh->bdr_attributes.Max());
+   // Inlet is attribute 1.
+   attr[0] = 1;
+   // Walls is attribute 3.
+   attr[2] = 1;
+   flowsolver.AddVelDirichletBC(vel, attr);
+
+   double t = 0.0;
+   double dt = ctx.dt;
+   double t_final = ctx.t_final;
+   bool last_step = false;
+
+   flowsolver.Setup(dt);
+
+   ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
+   ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
+
+   ParaViewDataCollection pvdc("3dfoc", pmesh);
+   pvdc.SetDataFormat(VTKFormat::BINARY32);
+   pvdc.SetHighOrderOutput(true);
+   pvdc.SetLevelsOfDetail(ctx.order);
+   pvdc.SetCycle(0);
+   pvdc.SetTime(t);
+   pvdc.RegisterField("velocity", u_gf);
+   pvdc.RegisterField("pressure", p_gf);
+   pvdc.Save();
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (t + dt >= t_final - dt / 2)
+      {
+         last_step = true;
+      }
+
+      flowsolver.Step(t, dt, step);
+
+      if (step % 10 == 0)
+      {
+         pvdc.SetCycle(step);
+         pvdc.SetTime(t);
+         pvdc.Save();
+      }
+
+      if (mpi.Root())
+      {
+         printf("%11s %11s\n", "Time", "dt");
+         printf("%.5E %.5E\n", t, dt);
+         fflush(stdout);
+      }
+   }
+
+   flowsolver.PrintTimingData();
+
+   delete pmesh;
+
+   return 0;
+}

miniapps/navier/navier_kovasznay.cpp

@@ -0,0 +1,264 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+// Navier Kovasznay example
+//
+// Solve for the steady Kovasznay flow at Re = 40 defined by
+//
+// u = [1 - exp(L * x) * cos(2 * pi * y),
+//      L / (2 * pi) * exp(L * x) * sin(2 * pi * y)],
+//
+// p = 1/2 * (1 - exp(2 * L * x)),
+//
+// with L = Re/2 - sqrt(Re^2/4 + 4 * pi^2).
+//
+// The problem domain is set up like this
+//
+//            +-------------+
+//            |             |
+//            |             |
+//            |             |
+//            |             |
+//  Inflow -> |             | -> Outflow
+//            |             |
+//            |             |
+//            |             |
+//            |             |
+//            |             |
+//            +-------------+
+//
+// and Dirichlet boundary conditions are applied for the velocity on every
+// boundary. The problem, although steady state, is time integrated up to the
+// final time and the solution is compared with the known exact solution.
+
+#include "navier_solver.hpp"
+#include <fstream>
+
+using namespace mfem;
+using namespace navier;
+
+struct s_NavierContext
+{
+   int ser_ref_levels = 1;
+   int order = 6;
+   double kinvis = 1.0 / 40.0;
+   double t_final = 10 * 0.001;
+   double dt = 0.001;
+   double reference_pressure = 0.0;
+   double reynolds = 1.0 / kinvis;
+   double lam = 0.5 * reynolds
+                - sqrt(0.25 * reynolds * reynolds + 4.0 * M_PI * M_PI);
+   bool pa = true;
+   bool ni = false;
+   bool visualization = false;
+   bool checkres = false;
+} ctx;
+
+void vel_kovasznay(const Vector &x, double t, Vector &u)
+{
+   double xi = x(0);
+   double yi = x(1);
+
+   u(0) = 1.0 - exp(ctx.lam * xi) * cos(2.0 * M_PI * yi);
+   u(1) = ctx.lam / (2.0 * M_PI) * exp(ctx.lam * xi) * sin(2.0 * M_PI * yi);
+}
+
+double pres_kovasznay(const Vector &x, double t)
+{
+   double xi = x(0);
+
+   return 0.5 * (1.0 - exp(2.0 * ctx.lam * xi)) + ctx.reference_pressure;
+}
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ctx.ser_ref_levels,
+                  "-rs",
+                  "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&ctx.order,
+                  "-o",
+                  "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
+   args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
+   args.AddOption(&ctx.pa,
+                  "-pa",
+                  "--enable-pa",
+                  "-no-pa",
+                  "--disable-pa",
+                  "Enable partial assembly.");
+   args.AddOption(&ctx.ni,
+                  "-ni",
+                  "--enable-ni",
+                  "-no-ni",
+                  "--disable-ni",
+                  "Enable numerical integration rules.");
+   args.AddOption(&ctx.visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(
+      &ctx.checkres,
+      "-cr",
+      "--checkresult",
+      "-no-cr",
+      "--no-checkresult",
+      "Enable or disable checking of the result. Returns -1 on failure.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root())
+      {
+         args.PrintUsage(mfem::out);
+      }
+      return 1;
+   }
+   if (mpi.Root())
+   {
+      args.PrintOptions(mfem::out);
+   }
+
+   Mesh *mesh = new Mesh(2, 4, Element::QUADRILATERAL, false, 1.5, 2.0);
+
+   mesh->EnsureNodes();
+   GridFunction *nodes = mesh->GetNodes();
+   *nodes -= 0.5;
+
+   for (int i = 0; i < ctx.ser_ref_levels; ++i)
+   {
+      mesh->UniformRefinement();
+   }
+
+   if (mpi.Root())
+   {
+      std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
+   }
+
+   auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Create the flow solver.
+   NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
+   flowsolver.EnablePA(ctx.pa);
+   flowsolver.EnableNI(ctx.ni);
+
+   // Set the initial condition.
+   ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
+   VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_kovasznay);
+   u_ic->ProjectCoefficient(u_excoeff);
+
+   FunctionCoefficient p_excoeff(pres_kovasznay);
+
+   // Add Dirichlet boundary conditions to velocity space restricted to
+   // selected attributes on the mesh.
+   Array<int> attr(pmesh->bdr_attributes.Max());
+   attr = 1;
+   flowsolver.AddVelDirichletBC(vel_kovasznay, attr);
+
+   double t = 0.0;
+   double dt = ctx.dt;
+   double t_final = ctx.t_final;
+   bool last_step = false;
+
+   flowsolver.Setup(dt);
+
+   double err_u = 0.0;
+   double err_p = 0.0;
+   ParGridFunction *u_gf = nullptr;
+   ParGridFunction *p_gf = nullptr;
+
+   ParGridFunction p_ex_gf(flowsolver.GetCurrentPressure()->ParFESpace());
+   GridFunctionCoefficient p_ex_gf_coeff(&p_ex_gf);
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (t + dt >= t_final - dt / 2)
+      {
+         last_step = true;
+      }
+
+      flowsolver.Step(t, dt, step);
+
+      // Compare against exact solution of velocity and pressure.
+      u_gf = flowsolver.GetCurrentVelocity();
+      p_gf = flowsolver.GetCurrentPressure();
+
+      u_excoeff.SetTime(t);
+      p_excoeff.SetTime(t);
+
+      // Remove mean value from exact pressure solution.
+      p_ex_gf.ProjectCoefficient(p_excoeff);
+      flowsolver.MeanZero(p_ex_gf);
+
+      err_u = u_gf->ComputeL2Error(u_excoeff);
+      err_p = p_gf->ComputeL2Error(p_ex_gf_coeff);
+
+      double cfl = flowsolver.ComputeCFL(*u_gf, dt);
+
+      if (mpi.Root())
+      {
+         printf("%5s %8s %8s %8s %11s %11s\n",
+                "Order",
+                "CFL",
+                "Time",
+                "dt",
+                "err_u",
+                "err_p");
+         printf("%5.2d %8.2E %.2E %.2E %.5E %.5E err\n",
+                ctx.order,
+                cfl,
+                t,
+                dt,
+                err_u,
+                err_p);
+         fflush(stdout);
+      }
+   }
+
+   if (ctx.visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "parallel " << mpi.WorldSize() << " " << mpi.WorldRank()
+               << "\n";
+      sol_sock << "solution\n" << *pmesh << *u_ic << std::flush;
+   }
+
+   flowsolver.PrintTimingData();
+
+   // Test if the result for the test run is as expected.
+   if (ctx.checkres)
+   {
+      double tol_u = 1e-6;
+      double tol_p = 1e-5;
+      if (err_u > tol_u || err_p > tol_p)
+      {
+         if (mpi.Root())
+         {
+            mfem::out << "Result has a larger error than expected."
+                      << std::endl;
+         }
+         return -1;
+      }
+   }
+
+   delete pmesh;
+
+   return 0;
+}

miniapps/navier/navier_mms.cpp

@@ -0,0 +1,248 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+// Navier MMS example
+//
+// A manufactured solution is defined as
+//
+// u = [pi * sin(t) * sin(pi * x)^2 * sin(2 * pi * y),
+//      -(pi * sin(t) * sin(2 * pi * x)) * sin(pi * y)^2].
+//
+// p = cos(pi * x) * sin(t) * sin(pi * y)
+//
+// The solution is used to compute the symbolic forcing term (right hand side),
+// of the equation. Then the numerical solution is computed and compared to the
+// exact manufactured solution to determine the error.
+
+#include "navier_solver.hpp"
+#include <fstream>
+
+using namespace mfem;
+using namespace navier;
+
+struct s_NavierContext
+{
+   int ser_ref_levels = 1;
+   int order = 5;
+   double kinvis = 1.0;
+   double t_final = 10 * 0.25e-4;
+   double dt = 0.25e-4;
+   bool pa = true;
+   bool ni = false;
+   bool visualization = false;
+   bool checkres = false;
+} ctx;
+
+void vel(const Vector &x, double t, Vector &u)
+{
+   double xi = x(0);
+   double yi = x(1);
+
+   u(0) = M_PI * sin(t) * pow(sin(M_PI * xi), 2.0) * sin(2.0 * M_PI * yi);
+   u(1) = -(M_PI * sin(t) * sin(2.0 * M_PI * xi) * pow(sin(M_PI * yi), 2.0));
+}
+
+double p(const Vector &x, double t)
+{
+   double xi = x(0);
+   double yi = x(1);
+
+   return cos(M_PI * xi) * sin(t) * sin(M_PI * yi);
+}
+
+void accel(const Vector &x, double t, Vector &u)
+{
+   double xi = x(0);
+   double yi = x(1);
+
+   u(0) = M_PI * sin(t) * sin(M_PI * xi) * sin(M_PI * yi)
+          * (-1.0
+             + 2.0 * pow(M_PI, 2.0) * sin(t) * sin(M_PI * xi)
+             * sin(2.0 * M_PI * xi) * sin(M_PI * yi))
+          + M_PI
+          * (2.0 * ctx.kinvis * pow(M_PI, 2.0)
+             * (1.0 - 2.0 * cos(2.0 * M_PI * xi)) * sin(t)
+             + cos(t) * pow(sin(M_PI * xi), 2.0))
+          * sin(2.0 * M_PI * yi);
+
+   u(1) = M_PI * cos(M_PI * yi) * sin(t)
+          * (cos(M_PI * xi)
+             + 2.0 * ctx.kinvis * pow(M_PI, 2.0) * cos(M_PI * yi)
+             * sin(2.0 * M_PI * xi))
+          - M_PI * (cos(t) + 6.0 * ctx.kinvis * pow(M_PI, 2.0) * sin(t))
+          * sin(2.0 * M_PI * xi) * pow(sin(M_PI * yi), 2.0)
+          + 4.0 * pow(M_PI, 3.0) * cos(M_PI * yi) * pow(sin(t), 2.0)
+          * pow(sin(M_PI * xi), 2.0) * pow(sin(M_PI * yi), 3.0);
+}
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ctx.ser_ref_levels,
+                  "-rs",
+                  "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&ctx.order,
+                  "-o",
+                  "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
+   args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
+   args.AddOption(&ctx.pa,
+                  "-pa",
+                  "--enable-pa",
+                  "-no-pa",
+                  "--disable-pa",
+                  "Enable partial assembly.");
+   args.AddOption(&ctx.ni,
+                  "-ni",
+                  "--enable-ni",
+                  "-no-ni",
+                  "--disable-ni",
+                  "Enable numerical integration rules.");
+   args.AddOption(&ctx.visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(
+      &ctx.checkres,
+      "-cr",
+      "--checkresult",
+      "-no-cr",
+      "--no-checkresult",
+      "Enable or disable checking of the result. Returns -1 on failure.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root())
+      {
+         args.PrintUsage(mfem::out);
+      }
+      return 1;
+   }
+   if (mpi.Root())
+   {
+      args.PrintOptions(mfem::out);
+   }
+
+   Mesh *mesh = new Mesh("../../data/inline-quad.mesh");
+   mesh->EnsureNodes();
+   GridFunction *nodes = mesh->GetNodes();
+   *nodes *= 2.0;
+   *nodes -= 1.0;
+
+   for (int i = 0; i < ctx.ser_ref_levels; ++i)
+   {
+      mesh->UniformRefinement();
+   }
+
+   if (mpi.Root())
+   {
+      std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
+   }
+
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Create the flow solver.
+   NavierSolver naviersolver(pmesh, ctx.order, ctx.kinvis);
+   naviersolver.EnablePA(ctx.pa);
+   naviersolver.EnableNI(ctx.ni);
+
+   // Set the initial condition.
+   ParGridFunction *u_ic = naviersolver.GetCurrentVelocity();
+   VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
+   u_ic->ProjectCoefficient(u_excoeff);
+
+   FunctionCoefficient p_excoeff(p);
+
+   // Add Dirichlet boundary conditions to velocity space restricted to
+   // selected attributes on the mesh.
+   Array<int> attr(pmesh->bdr_attributes.Max());
+   attr = 1;
+   naviersolver.AddVelDirichletBC(vel, attr);
+
+   Array<int> domain_attr(pmesh->attributes.Max());
+   domain_attr = 1;
+   naviersolver.AddAccelTerm(accel, domain_attr);
+
+   double t = 0.0;
+   double dt = ctx.dt;
+   double t_final = ctx.t_final;
+   bool last_step = false;
+
+   naviersolver.Setup(dt);
+
+   double err_u = 0.0;
+   double err_p = 0.0;
+   ParGridFunction *u_gf = nullptr;
+   ParGridFunction *p_gf = nullptr;
+   u_gf = naviersolver.GetCurrentVelocity();
+   p_gf = naviersolver.GetCurrentPressure();
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (t + dt >= t_final - dt / 2)
+      {
+         last_step = true;
+      }
+
+      naviersolver.Step(t, dt, step);
+
+      // Compare against exact solution of velocity and pressure.
+      u_excoeff.SetTime(t);
+      p_excoeff.SetTime(t);
+      err_u = u_gf->ComputeL2Error(u_excoeff);
+      err_p = p_gf->ComputeL2Error(p_excoeff);
+
+      if (mpi.Root())
+      {
+         printf("%11s %11s %11s %11s\n", "Time", "dt", "err_u", "err_p");
+         printf("%.5E %.5E %.5E %.5E err\n", t, dt, err_u, err_p);
+         fflush(stdout);
+      }
+   }
+
+   if (ctx.visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << mpi.WorldSize() << " " << mpi.WorldRank()
+               << "\n";
+      sol_sock << "solution\n" << *pmesh << *u_ic << std::flush;
+   }
+
+   naviersolver.PrintTimingData();
+
+   // Test if the result for the test run is as expected.
+   if (ctx.checkres)
+   {
+      double tol = 1e-3;
+      if (err_u > tol || err_p > tol)
+      {
+         if (mpi.Root())
+         {
+            mfem::out << "Result has a larger error than expected."
+                      << std::endl;
+         }
+         return -1;
+      }
+   }
+
+   delete pmesh;
+
+   return 0;
+}

miniapps/navier/navier_shear.cpp

@@ -0,0 +1,153 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+// Navier double shear layer example
+//
+// Solve the double shear problem in the following configuration.
+//
+//       +-------------------+
+//       |                   |
+//       |      u0 = ua      |
+//       |                   |
+//  -------------------------------- y = 0.5
+//       |                   |
+//       |      u0 = ub      |
+//       |                   |
+//       +-------------------+
+//
+// The initial condition u0 is chosen to be a varying velocity in the y
+// direction. It includes a perturbation at x = 0.5 which leads to an
+// instability and the dynamics of the flow. The boundary conditions are fully
+// periodic.
+
+#include "navier_solver.hpp"
+#include <fstream>
+
+using namespace mfem;
+using namespace navier;
+
+struct s_NavierContext
+{
+   int order = 6;
+   double kinvis = 1.0 / 100000.0;
+   double t_final = 10 * 1e-3;
+   double dt = 1e-3;
+} ctx;
+
+void vel_shear_ic(const Vector &x, double t, Vector &u)
+{
+   double xi = x(0);
+   double yi = x(1);
+
+   double rho = 30.0;
+   double delta = 0.05;
+
+   if (yi <= 0.5)
+   {
+      u(0) = tanh(rho * (yi - 0.25));
+   }
+   else
+   {
+      u(0) = tanh(rho * (0.75 - yi));
+   }
+
+   u(1) = delta * sin(2.0 * M_PI * xi);
+}
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   int serial_refinements = 2;
+
+   Mesh *mesh = new Mesh("../../data/periodic-square.mesh");
+   mesh->EnsureNodes();
+   GridFunction *nodes = mesh->GetNodes();
+   *nodes -= -1.0;
+   *nodes /= 2.0;
+
+   for (int i = 0; i < serial_refinements; ++i)
+   {
+      mesh->UniformRefinement();
+   }
+
+   if (mpi.Root())
+   {
+      std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
+   }
+
+   auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Create the flow solver.
+   NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
+   flowsolver.EnablePA(true);
+
+   // Set the initial condition.
+   ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
+   VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_shear_ic);
+   u_ic->ProjectCoefficient(u_excoeff);
+
+   double t = 0.0;
+   double dt = ctx.dt;
+   double t_final = ctx.t_final;
+   bool last_step = false;
+
+   flowsolver.Setup(dt);
+
+   ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
+   ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
+
+   ParGridFunction w_gf(*u_gf);
+   flowsolver.ComputeCurl2D(*u_gf, w_gf);
+
+   ParaViewDataCollection pvdc("shear_output", pmesh);
+   pvdc.SetDataFormat(VTKFormat::BINARY32);
+   pvdc.SetHighOrderOutput(true);
+   pvdc.SetLevelsOfDetail(ctx.order);
+   pvdc.SetCycle(0);
+   pvdc.SetTime(t);
+   pvdc.RegisterField("velocity", u_gf);
+   pvdc.RegisterField("pressure", p_gf);
+   pvdc.RegisterField("vorticity", &w_gf);
+   pvdc.Save();
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (t + dt >= t_final - dt / 2)
+      {
+         last_step = true;
+      }
+
+      flowsolver.Step(t, dt, step);
+
+      if (step % 10 == 0)
+      {
+         flowsolver.ComputeCurl2D(*u_gf, w_gf);
+         pvdc.SetCycle(step);
+         pvdc.SetTime(t);
+         pvdc.Save();
+      }
+
+      if (mpi.Root())
+      {
+         printf("%11s %11s\n", "Time", "dt");
+         printf("%.5E %.5E\n", t, dt);
+         fflush(stdout);
+      }
+   }
+
+   flowsolver.PrintTimingData();
+
+   delete pmesh;
+
+   return 0;
+}

miniapps/navier/navier_solver.hpp

@@ -0,0 +1,388 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_NAVIER_SOLVER_HPP
+#define MFEM_NAVIER_SOLVER_HPP
+
+#define NAVIER_VERSION 0.1
+
+#include "mfem.hpp"
+#include "ortho_solver.hpp"
+
+namespace mfem
+{
+namespace navier
+{
+using VecFuncT = void(const Vector &x, double t, Vector &u);
+using ScalarFuncT = double(const Vector &x, double t);
+
+/// Container for a Dirichlet boundary condition of the velocity field.
+class VelDirichletBC_T
+{
+public:
+   VelDirichletBC_T(Array<int> attr, VectorCoefficient *coeff)
+      : attr(attr), coeff(coeff)
+   {}
+
+   ~VelDirichletBC_T() { delete coeff; }
+
+   Array<int> attr;
+   VectorCoefficient *coeff;
+};
+
+/// Container for a Dirichlet boundary condition of the pressure field.
+class PresDirichletBC_T
+{
+public:
+   PresDirichletBC_T(Array<int> attr, Coefficient *coeff)
+      : attr(attr), coeff(coeff)
+   {}
+
+   ~PresDirichletBC_T() { delete coeff; }
+
+   Array<int> attr;
+   Coefficient *coeff;
+};
+
+/// Container for an acceleration term.
+class AccelTerm_T
+{
+public:
+   AccelTerm_T(Array<int> attr, VectorCoefficient *coeff)
+      : attr(attr), coeff(coeff)
+   {}
+
+   ~AccelTerm_T() { delete coeff; }
+
+   Array<int> attr;
+   VectorCoefficient *coeff;
+};
+
+/// Transient incompressible Navier Stokes solver in a split scheme formulation.
+/**
+ * This implementation of a transient incompressible Navier Stokes solver uses
+ * the non-dimensionalized formulation. The coupled momentum and
+ * incompressibilty equations are decoupled using the split scheme described in
+ * [1]. This leads to three solving steps.
+ *
+ * 1. An extrapolation step for all nonlinear terms which are treated
+ *    explicitly. This step avoids a fully coupled nonlinear solve and only
+ *    requires a solve of the mass matrix in velocity space \f$M_v^{-1}\f$. On
+ *    the other hand this introduces a CFL stability condition on the maximum
+ *    timestep.
+ *
+ * 2. A Poisson solve \f$S_p^{-1}\f$.
+ *
+ * 3. A Helmholtz like solve \f$(M_v - \partial t K_v)^{-1}\f$.
+ *
+ * The numerical solver setup for each step are as follows.
+ *
+ * \f$M_v^{-1}\f$ is solved using CG with Jacobi as preconditioner.
+ *
+ * \f$S_p^{-1}\f$ is solved using CG with AMG applied to the low order refined
+ * (LOR) assembled pressure poisson matrix. To avoid assembling a matrix for
+ * preconditioning, one can use p-MG as an alternative (NYI).
+ *
+ * \f$(M_v - \partial t K_v)^{-1}\f$ due to the CFL condition we expect the time
+ * step to be small. Therefore this is solved using CG with Jacobi as
+ * preconditioner. For large time steps a preconditioner like AMG or p-MG should
+ * be used (NYI).
+ *
+ * Statements marked with NYI mean this feature is planned but Not Yet
+ * Implemented.
+ *
+ * A detailed description is available in [1] section 4.2. The algorithm is
+ * originated from [2].
+ *
+ * [1] Michael Franco, Jean-Sylvain Camier, Julian Andrej, Will Pazner (2020)
+ * High-order matrix-free incompressible flow solvers with GPU acceleration and
+ * low-order refined preconditioners (https://arxiv.org/abs/1910.03032)
+ *
+ * [2] A. G. Tomboulides, J. C. Y. Lee & S. A. Orszag (1997) Numerical
+ * Simulation of Low Mach Number Reactive Flows
+ */
+class NavierSolver
+{
+public:
+   /// Initialize data structures, set FE space order and kinematic viscosity.
+   /**
+    * The ParMesh @a mesh can be a linear or curved parallel mesh. The @a order
+    * of the finite element spaces is this algorithm is of equal order
+    * \f$(P_N)^d P_N\f$ for velocity and pressure respectively. This means the
+    * pressure is in discretized in the same space (just scalar instead of a
+    * vector space) as the velocity.
+    *
+    * Kinematic viscosity (dimensionless) is set using @a kin_vis and
+    * automatically converted to the Reynolds number. If you want to set the
+    * Reynolds number directly, you can provide the inverse.
+    */
+   NavierSolver(ParMesh *mesh, int order, double kin_vis);
+
+   /// Initialize forms, solvers and preconditioners.
+   void Setup(double dt);
+
+   /// Compute solution at the next time step t+dt.
+   void Step(double &time, double dt, int cur_step);
+
+   /// Return a pointer to the current velocity ParGridFunction.
+   ParGridFunction *GetCurrentVelocity() { return &un_gf; }
+
+   /// Return a pointer to the current pressure ParGridFunction.
+   ParGridFunction *GetCurrentPressure() { return &pn_gf; }
+
+   /// Add a Dirichlet boundary condition to the velocity field.
+   void AddVelDirichletBC(VectorCoefficient *coeff, Array<int> &attr);
+
+   void AddVelDirichletBC(VecFuncT *f, Array<int> &attr);
+
+   /// Add a Dirichlet boundary condition to the pressure field.
+   void AddPresDirichletBC(Coefficient *coeff, Array<int> &attr);
+
+   void AddPresDirichletBC(ScalarFuncT *f, Array<int> &attr);
+
+   /// Add an accelaration term to the RHS of the equation.
+   /**
+    * The VecFuncT @a f is evaluated at the current time t and extrapolated
+    * together with the nonlinear parts of the Navier Stokes equation.
+    */
+   void AddAccelTerm(VectorCoefficient *coeff, Array<int> &attr);
+
+   void AddAccelTerm(VecFuncT *f, Array<int> &attr);
+
+   /// Enable partial assembly for every operator.
+   void EnablePA(bool pa) { partial_assembly = pa; }
+
+   /// Enable numerical integration rules. This means collocated quadrature at
+   /// the nodal points.
+   void EnableNI(bool ni) { numerical_integ = ni; }
+
+   /// Print timing summary of the solving routine.
+   /**
+    * The summary shows the timing in seconds in the first row of
+    *
+    * 1. SETUP: Time spent for the setup of all forms, solvers and
+    *    preconditioners.
+    * 2. STEP: Time spent computing a full time step. It includes allthree
+    *    solves.
+    * 3. EXTRAP: Time spent for extrapolation of all forcing and nonlinear
+    *    terms.
+    * 4. CURLCURL: Time spent for computing the curl curl term in the pressure
+    *    Poisson equation (see references for detailed explanation).
+    * 5. PSOLVE: Time spent in the pressure Poisson solve.
+    * 6. HSOLVE: Time spent in the Helmholtz solve.
+    *
+    * The second row shows a proportion of a column relative to the whole
+    * time step.
+    */
+   void PrintTimingData();
+
+   ~NavierSolver();
+
+   /// Compute \f$\nabla \times \nabla \times u\f$ for \f$u \in (H^1)^2\f$.
+   void ComputeCurl2D(ParGridFunction &u,
+                      ParGridFunction &cu,
+                      bool assume_scalar = false);
+
+   /// Compute \f$\nabla \times \nabla \times u\f$ for \f$u \in (H^1)^3\f$.
+   void ComputeCurl3D(ParGridFunction &u, ParGridFunction &cu);
+
+   /// Remove mean from a Vector.
+   /**
+    * Modify the Vector @a v by subtracting its mean using
+    * \f$v = v - \frac{\sum_i^N v_i}{N} \f$
+    */
+   void Orthogonalize(Vector &v);
+
+   /// Remove the mean from a ParGridFunction.
+   /**
+    * Modify the ParGridFunction @a v by subtracting its mean using
+    * \f$ v = v - \int_\Omega \frac{v}{vol(\Omega)} dx \f$.
+    */
+   void MeanZero(ParGridFunction &v);
+
+   /// Compute CFL
+   double ComputeCFL(ParGridFunction &u, double dt);
+
+protected:
+   /// Print informations about the Navier version.
+   void PrintInfo();
+
+   /// Update the EXTk/BDF time integration coefficient.
+   /**
+    * Depending on which time step the computation is in, the EXTk/BDF time
+    * integration coefficients have to be set accordingly. This allows
+    * bootstrapping with a BDF scheme of order 1 and increasing the order each
+    * following time step, up to order 3.
+    */
+   void SetTimeIntegrationCoefficients(int step);
+
+   /// Eliminate essential BCs in an Operator and apply to RHS.
+   void EliminateRHS(Operator &A,
+                     ConstrainedOperator &constrainedA,
+                     const Array<int> &ess_tdof_list,
+                     Vector &x,
+                     Vector &b,
+                     Vector &X,
+                     Vector &B,
+                     int copy_interior = 0);
+
+   /// Enable/disable debug output.
+   bool debug = false;
+
+   /// Enable/disable verbose output.
+   bool verbose = true;
+
+   /// Enable/disable partial assembly of forms.
+   bool partial_assembly = false;
+
+   /// Enable/disable numerical integration rules of forms.
+   bool numerical_integ = false;
+
+   /// The parallel mesh.
+   ParMesh *pmesh = nullptr;
+
+   /// The order of the velocity and pressure space.
+   int order;
+
+   /// Kinematic viscosity (dimensionless).
+   double kin_vis;
+
+   /// Velocity \f$H^1\f$ finite element collection.
+   FiniteElementCollection *vfec = nullptr;
+
+   /// Pressure \f$H^1\f$ finite element collection.
+   FiniteElementCollection *pfec = nullptr;
+
+   /// Velocity \f$(H^1)^d\f$ finite element space.
+   ParFiniteElementSpace *vfes = nullptr;
+
+   /// Pressure \f$H^1\f$ finite element space.
+   ParFiniteElementSpace *pfes = nullptr;
+
+   ParNonlinearForm *N = nullptr;
+
+   ParBilinearForm *Mv_form = nullptr;
+
+   ParBilinearForm *Sp_form = nullptr;
+
+   ParMixedBilinearForm *D_form = nullptr;
+
+   ParMixedBilinearForm *G_form = nullptr;
+
+   ParBilinearForm *H_form = nullptr;
+
+   VectorGridFunctionCoefficient *FText_gfcoeff = nullptr;
+
+   ParLinearForm *FText_bdr_form = nullptr;
+
+   ParLinearForm *f_form = nullptr;
+
+   ParLinearForm *g_bdr_form = nullptr;
+
+   /// Linear form to compute the mass matrix in various subroutines.
+   ParLinearForm *mass_lf = nullptr;
+   ConstantCoefficient onecoeff;
+   double volume = 0.0;
+
+   ConstantCoefficient nlcoeff;
+   ConstantCoefficient Sp_coeff;
+   ConstantCoefficient H_lincoeff;
+   ConstantCoefficient H_bdfcoeff;
+
+   OperatorHandle Mv;
+   OperatorHandle Sp;
+   OperatorHandle D;
+   OperatorHandle G;
+   OperatorHandle H;
+
+   Solver *MvInvPC = nullptr;
+   CGSolver *MvInv = nullptr;
+
+   HypreBoomerAMG *SpInvPC = nullptr;
+   OrthoSolver *SpInvOrthoPC = nullptr;
+   CGSolver *SpInv = nullptr;
+
+   Solver *HInvPC = nullptr;
+   CGSolver *HInv = nullptr;
+
+   Vector fn, un, unm1, unm2, Nun, Nunm1, Nunm2, Fext, FText, Lext, resu;
+   Vector tmp1;
+
+   Vector pn, resp, FText_bdr, g_bdr;
+
+   ParGridFunction un_gf, curlu_gf, curlcurlu_gf, Lext_gf, FText_gf, resu_gf;
+
+   ParGridFunction pn_gf, resp_gf;
+
+   // All essential attributes.
+   Array<int> vel_ess_attr;
+   Array<int> pres_ess_attr;
+
+   // All essential true dofs.
+   Array<int> vel_ess_tdof;
+   Array<int> pres_ess_tdof;
+
+   // Bookkeeping for velocity dirichlet bcs.
+   std::vector<VelDirichletBC_T> vel_dbcs;
+
+   // Bookkeeping for pressure dirichlet bcs.
+   std::vector<PresDirichletBC_T> pres_dbcs;
+
+   // Bookkeeping for acceleration (forcing) terms.
+   std::vector<AccelTerm_T> accel_terms;
+
+   int cur_step = 0;
+
+   // BDFk/EXTk coefficients.
+   double bd0 = 0.0;
+   double bd1 = 0.0;
+   double bd2 = 0.0;
+   double bd3 = 0.0;
+   double ab1 = 0.0;
+   double ab2 = 0.0;
+   double ab3 = 0.0;
+
+   // Timers.
+   StopWatch sw_setup, sw_step, sw_extrap, sw_curlcurl, sw_spsolve, sw_hsolve;
+
+   // Print levels.
+   int pl_mvsolve = 0;
+   int pl_spsolve = 0;
+   int pl_hsolve = 0;
+   int pl_amg = 0;
+
+   // Relative tolerances.
+   double rtol_spsolve = 1e-6;
+   double rtol_hsolve = 1e-8;
+
+   // Iteration counts.
+   int iter_mvsolve = 0, iter_spsolve = 0, iter_hsolve = 0;
+
+   // Residuals.
+   double res_mvsolve = 0.0, res_spsolve = 0.0, res_hsolve = 0.0;
+
+   // LOR related.
+   ParMesh *pmesh_lor = nullptr;
+   FiniteElementCollection *pfec_lor = nullptr;
+   ParFiniteElementSpace *pfes_lor = nullptr;
+   InterpolationGridTransfer *vgt = nullptr, *pgt = nullptr;
+
+   ParBilinearForm *Mv_form_lor = nullptr;
+   ParBilinearForm *Sp_form_lor = nullptr;
+   ParBilinearForm *H_form_lor = nullptr;
+
+   OperatorHandle Mv_lor;
+   OperatorHandle Sp_lor;
+   OperatorHandle H_lor;
+};
+} // namespace navier
+} // namespace mfem
+#endif

miniapps/navier/navier_tgv.cpp

@@ -0,0 +1,405 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// 3d Taylor-Green vortex benchmark example at Re=1600
+// Unsteady flow of a decaying vortex is computed and compared against a known,
+// analytical solution.
+
+#include "navier_solver.hpp"
+#include <fstream>
+
+using namespace mfem;
+using namespace navier;
+
+struct s_NavierContext
+{
+   int element_subdivisions = 1;
+   int order = 4;
+   double kinvis = 1.0 / 1600.0;
+   double t_final = 10 * 1e-3;
+   double dt = 1e-3;
+   bool pa = true;
+   bool ni = false;
+   bool visualization = false;
+   bool checkres = false;
+} ctx;
+
+void vel_tgv(const Vector &x, double t, Vector &u)
+{
+   double xi = x(0);
+   double yi = x(1);
+   double zi = x(2);
+
+   u(0) = sin(xi) * cos(yi) * cos(zi);
+   u(1) = -cos(xi) * sin(yi) * cos(zi);
+   u(2) = 0.0;
+}
+
+class QuantitiesOfInterest
+{
+public:
+   QuantitiesOfInterest(ParMesh *pmesh)
+   {
+      H1_FECollection h1fec(1);
+      ParFiniteElementSpace h1fes(pmesh, &h1fec);
+
+      onecoeff.constant = 1.0;
+      mass_lf = new ParLinearForm(&h1fes);
+      mass_lf->AddDomainIntegrator(new DomainLFIntegrator(onecoeff));
+      mass_lf->Assemble();
+
+      ParGridFunction one_gf(&h1fes);
+      one_gf.ProjectCoefficient(onecoeff);
+
+      volume = mass_lf->operator()(one_gf);
+   };
+
+   double ComputeKineticEnergy(ParGridFunction &v)
+   {
+      Vector velx, vely, velz;
+      double integ = 0.0;
+      const FiniteElement *fe;
+      ElementTransformation *T;
+      FiniteElementSpace *fes = v.FESpace();
+
+      for (int i = 0; i < fes->GetNE(); i++)
+      {
+         fe = fes->GetFE(i);
+         double intorder = 2 * fe->GetOrder();
+         const IntegrationRule *ir = &(
+                                        IntRules.Get(fe->GetGeomType(), intorder));
+
+         v.GetValues(i, *ir, velx, 1);
+         v.GetValues(i, *ir, vely, 2);
+         v.GetValues(i, *ir, velz, 3);
+
+         T = fes->GetElementTransformation(i);
+         for (int j = 0; j < ir->GetNPoints(); j++)
+         {
+            const IntegrationPoint &ip = ir->IntPoint(j);
+            T->SetIntPoint(&ip);
+
+            double vel2 = velx(j) * velx(j) + vely(j) * vely(j)
+                          + velz(j) * velz(j);
+
+            integ += ip.weight * T->Weight() * vel2;
+         }
+      }
+
+      double global_integral = 0.0;
+      MPI_Allreduce(&integ,
+                    &global_integral,
+                    1,
+                    MPI_DOUBLE,
+                    MPI_SUM,
+                    MPI_COMM_WORLD);
+
+      return 0.5 * global_integral / volume;
+   };
+
+   ~QuantitiesOfInterest() { delete mass_lf; };
+
+private:
+   ConstantCoefficient onecoeff;
+   ParLinearForm *mass_lf;
+   double volume;
+};
+
+template<typename T>
+T sq(T x)
+{
+   return x * x;
+}
+
+// Computes Q = 0.5*(tr(\nabla u)^2 - tr(\nabla u \cdot \nabla u))
+void ComputeQCriterion(ParGridFunction &u, ParGridFunction &q)
+{
+   FiniteElementSpace *v_fes = u.FESpace();
+   FiniteElementSpace *fes = q.FESpace();
+
+   // AccumulateAndCountZones
+   Array<int> zones_per_vdof;
+   zones_per_vdof.SetSize(fes->GetVSize());
+   zones_per_vdof = 0;
+
+   q = 0.0;
+
+   // Local interpolation
+   int elndofs;
+   Array<int> v_dofs, dofs;
+   Vector vals;
+   Vector loc_data;
+   int vdim = v_fes->GetVDim();
+   DenseMatrix grad_hat;
+   DenseMatrix dshape;
+   DenseMatrix grad;
+
+   for (int e = 0; e < fes->GetNE(); ++e)
+   {
+      fes->GetElementVDofs(e, dofs);
+      v_fes->GetElementVDofs(e, v_dofs);
+      u.GetSubVector(v_dofs, loc_data);
+      vals.SetSize(dofs.Size());
+      ElementTransformation *tr = fes->GetElementTransformation(e);
+      const FiniteElement *el = fes->GetFE(e);
+      elndofs = el->GetDof();
+      int dim = el->GetDim();
+      dshape.SetSize(elndofs, dim);
+
+      for (int dof = 0; dof < elndofs; ++dof)
+      {
+         // Project
+         const IntegrationPoint &ip = el->GetNodes().IntPoint(dof);
+         tr->SetIntPoint(&ip);
+
+         // Eval
+         // GetVectorGradientHat
+         el->CalcDShape(tr->GetIntPoint(), dshape);
+         grad_hat.SetSize(vdim, dim);
+         DenseMatrix loc_data_mat(loc_data.GetData(), elndofs, vdim);
+         MultAtB(loc_data_mat, dshape, grad_hat);
+
+         const DenseMatrix &Jinv = tr->InverseJacobian();
+         grad.SetSize(grad_hat.Height(), Jinv.Width());
+         Mult(grad_hat, Jinv, grad);
+
+         double q_val = 0.5 * (sq(grad(0, 0)) + sq(grad(1, 1)) + sq(grad(2, 2)))
+                        + grad(0, 1) * grad(1, 0) + grad(0, 2) * grad(2, 0)
+                        + grad(1, 2) * grad(2, 1);
+
+         vals(dof) = q_val;
+      }
+
+      // Accumulate values in all dofs, count the zones.
+      for (int j = 0; j < dofs.Size(); j++)
+      {
+         int ldof = dofs[j];
+         q(ldof) += vals[j];
+         zones_per_vdof[ldof]++;
+      }
+   }
+
+   // Communication
+
+   // Count the zones globally.
+   GroupCommunicator &gcomm = q.ParFESpace()->GroupComm();
+   gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
+   gcomm.Bcast(zones_per_vdof);
+
+   // Accumulate for all vdofs.
+   gcomm.Reduce<double>(q.GetData(), GroupCommunicator::Sum);
+   gcomm.Bcast<double>(q.GetData());
+
+   // Compute means
+   for (int i = 0; i < q.Size(); i++)
+   {
+      const int nz = zones_per_vdof[i];
+      if (nz)
+      {
+         q(i) /= nz;
+      }
+   }
+}
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ctx.element_subdivisions,
+                  "-es",
+                  "--element-subdivisions",
+                  "Number of 1d uniform subdivisions for each element.");
+   args.AddOption(&ctx.order,
+                  "-o",
+                  "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
+   args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
+   args.AddOption(&ctx.pa,
+                  "-pa",
+                  "--enable-pa",
+                  "-no-pa",
+                  "--disable-pa",
+                  "Enable partial assembly.");
+   args.AddOption(&ctx.ni,
+                  "-ni",
+                  "--enable-ni",
+                  "-no-ni",
+                  "--disable-ni",
+                  "Enable numerical integration rules.");
+   args.AddOption(&ctx.visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(
+      &ctx.checkres,
+      "-cr",
+      "--checkresult",
+      "-no-cr",
+      "--no-checkresult",
+      "Enable or disable checking of the result. Returns -1 on failure.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root())
+      {
+         args.PrintUsage(mfem::out);
+      }
+      return 1;
+   }
+   if (mpi.Root())
+   {
+      args.PrintOptions(mfem::out);
+   }
+
+   Mesh *orig_mesh = new Mesh("../../data/periodic-cube.mesh");
+   Mesh *mesh = new Mesh(orig_mesh,
+                         ctx.element_subdivisions,
+                         BasisType::ClosedUniform);
+   delete orig_mesh;
+
+   mesh->EnsureNodes();
+   GridFunction *nodes = mesh->GetNodes();
+   *nodes *= M_PI;
+
+   int nel = mesh->GetNE();
+   if (mpi.Root())
+   {
+      mfem::out << "Number of elements: " << nel << std::endl;
+   }
+
+   auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Create the flow solver.
+   NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
+   flowsolver.EnablePA(ctx.pa);
+   flowsolver.EnableNI(ctx.ni);
+
+   // Set the initial condition.
+   ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
+   VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_tgv);
+   u_ic->ProjectCoefficient(u_excoeff);
+
+   double t = 0.0;
+   double dt = ctx.dt;
+   double t_final = ctx.t_final;
+   bool last_step = false;
+
+   flowsolver.Setup(dt);
+
+   ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
+   ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
+
+   ParGridFunction w_gf(*u_gf);
+   ParGridFunction q_gf(*p_gf);
+   flowsolver.ComputeCurl3D(*u_gf, w_gf);
+   ComputeQCriterion(*u_gf, q_gf);
+
+   QuantitiesOfInterest kin_energy(pmesh);
+
+   ParaViewDataCollection pvdc("shear_output", pmesh);
+   pvdc.SetDataFormat(VTKFormat::BINARY32);
+   pvdc.SetHighOrderOutput(true);
+   pvdc.SetLevelsOfDetail(ctx.order);
+   pvdc.SetCycle(0);
+   pvdc.SetTime(t);
+   pvdc.RegisterField("velocity", u_gf);
+   pvdc.RegisterField("pressure", p_gf);
+   pvdc.RegisterField("vorticity", &w_gf);
+   pvdc.RegisterField("qcriterion", &q_gf);
+   pvdc.Save();
+
+   double u_inf_loc = u_gf->Normlinf();
+   double p_inf_loc = p_gf->Normlinf();
+   double u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
+   double p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
+   double ke = kin_energy.ComputeKineticEnergy(*u_gf);
+
+   std::string fname = "tgv_out_p_" + std::to_string(ctx.order) + ".txt";
+   FILE *f;
+
+   if (mpi.Root())
+   {
+      int nel1d = std::round(pow(nel, 1.0 / 3.0));
+      int ngridpts = p_gf->ParFESpace()->GlobalVSize();
+      printf("%11s %11s %11s %11s %11s\n", "Time", "dt", "u_inf", "p_inf", "ke");
+      printf("%.5E %.5E %.5E %.5E %.5E\n", t, dt, u_inf, p_inf, ke);
+
+      f = fopen(fname.c_str(), "w");
+      fprintf(f, "3D Taylor Green Vortex\n");
+      fprintf(f, "order = %d\n", ctx.order);
+      fprintf(f, "grid = %d x %d x %d\n", nel1d, nel1d, nel1d);
+      fprintf(f, "dofs per component = %d\n", ngridpts);
+      fprintf(f, "=================================================\n");
+      fprintf(f, "        time                   kinetic energy\n");
+      fprintf(f, "%20.16e     %20.16e\n", t, ke);
+      fflush(f);
+      fflush(stdout);
+   }
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (t + dt >= t_final - dt / 2)
+      {
+         last_step = true;
+      }
+
+      flowsolver.Step(t, dt, step);
+
+      if ((step + 1) % 100 == 0 || last_step)
+      {
+         flowsolver.ComputeCurl3D(*u_gf, w_gf);
+         ComputeQCriterion(*u_gf, q_gf);
+         pvdc.SetCycle(step);
+         pvdc.SetTime(t);
+         pvdc.Save();
+      }
+
+      double u_inf_loc = u_gf->Normlinf();
+      double p_inf_loc = p_gf->Normlinf();
+      double u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
+      double p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
+      double ke = kin_energy.ComputeKineticEnergy(*u_gf);
+      if (mpi.Root())
+      {
+         printf("%.5E %.5E %.5E %.5E %.5E\n", t, dt, u_inf, p_inf, ke);
+         fprintf(f, "%20.16e     %20.16e\n", t, ke);
+         fflush(f);
+         fflush(stdout);
+      }
+   }
+
+   flowsolver.PrintTimingData();
+
+   // Test if the result for the test run is as expected.
+   if (ctx.checkres)
+   {
+      double tol = 1e-5;
+      double ke_expected = 1.25e-1;
+      if (fabs(ke - ke_expected) > tol)
+      {
+         if (mpi.Root())
+         {
+            mfem::out << "Result has a larger error than expected."
+                      << std::endl;
+         }
+         return -1;
+      }
+   }
+
+   delete pmesh;
+
+   return 0;
+}

miniapps/navier/ortho_solver.cpp

@@ -0,0 +1,52 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ortho_solver.hpp"
+
+using namespace mfem;
+using namespace navier;
+
+OrthoSolver::OrthoSolver() : Solver(0, true) {}
+
+void OrthoSolver::SetOperator(const Operator &op)
+{
+   oper = &op;
+}
+
+void OrthoSolver::Mult(const Vector &b, Vector &x) const
+{
+   // Orthoganlize input
+   Orthogonalize(b, b_ortho);
+
+   // Apply operator
+   oper->Mult(b_ortho, x);
+
+   // Orthoganlize output
+   Orthogonalize(x, x);
+}
+
+void OrthoSolver::Orthogonalize(const Vector &v, Vector &v_ortho) const
+{
+   double loc_sum = v.Sum();
+   double global_sum = 0.0;
+   int loc_size = v.Size();
+   int global_size = 0;
+
+   MPI_Allreduce(&loc_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   MPI_Allreduce(&loc_size, &global_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+   double ratio = global_sum / static_cast<double>(global_size);
+   v_ortho.SetSize(v.Size());
+   for (int i = 0; i < v_ortho.Size(); ++i)
+   {
+      v_ortho(i) = v(i) - ratio;
+   }
+}

miniapps/navier/ortho_solver.hpp

@@ -0,0 +1,47 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_NAVIER_ORTHO_SOLVER_HPP
+#define MFEM_NAVIER_ORTHO_SOLVER_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+namespace navier
+{
+/// Solver wrapper which orthogonalizes the input and output vector
+/**
+ * OrthoSolver wraps an existing Operator and orthogonalizes the input vector
+ * before passing it to the Mult method of the Operator. This is a convenience
+ * implementation to handle e.g. a Poisson problem with pure Neumann boundary
+ * conditions, where this procedure removes the Nullspace.
+ */
+class OrthoSolver : public Solver
+{
+public:
+   OrthoSolver();
+
+   virtual void SetOperator(const Operator &op);
+
+   void Mult(const Vector &b, Vector &x) const;
+
+private:
+   const Operator *oper = nullptr;
+
+   mutable Vector b_ortho;
+
+   void Orthogonalize(const Vector &v, Vector &v_ortho) const;
+};
+} // namespace navier
+} // namespace mfem
+
+#endif

tests/unit/fem/test_bilinearform.cpp

@@ -0,0 +1,67 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "catch.hpp"
+
+#include <iostream>
+
+using namespace mfem;
+
+TEST_CASE("Test order of boundary integrators",
+          "[BilinearForm]")
+{
+   // Create a simple mesh
+   int dim = 2, nx = 2, ny = 2, order = 2;
+   Element::Type e_type = Element::QUADRILATERAL;
+   Mesh mesh(nx, ny, e_type);
+
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   SECTION("Order of restricted boundary integrators")
+   {
+      ConstantCoefficient one(1.0);
+      ConstantCoefficient two(2.0);
+      ConstantCoefficient three(3.0);
+      ConstantCoefficient four(4.0);
+
+      Array<int> bdr1(4); bdr1 = 0; bdr1[0] = 1;
+      Array<int> bdr2(4); bdr2 = 0; bdr2[1] = 1;
+      Array<int> bdr3(4); bdr3 = 0; bdr3[2] = 1;
+      Array<int> bdr4(4); bdr4 = 0; bdr4[3] = 1;
+
+      BilinearForm a1234(&fes);
+      a1234.AddBoundaryIntegrator(new MassIntegrator(one), bdr1);
+      a1234.AddBoundaryIntegrator(new MassIntegrator(two), bdr2);
+      a1234.AddBoundaryIntegrator(new MassIntegrator(three), bdr3);
+      a1234.AddBoundaryIntegrator(new MassIntegrator(four), bdr4);
+      a1234.Assemble(0);
+      a1234.Finalize(0);
+
+      BilinearForm a4321(&fes);
+      a4321.AddBoundaryIntegrator(new MassIntegrator(four), bdr4);
+      a4321.AddBoundaryIntegrator(new MassIntegrator(three), bdr3);
+      a4321.AddBoundaryIntegrator(new MassIntegrator(two), bdr2);
+      a4321.AddBoundaryIntegrator(new MassIntegrator(one), bdr1);
+      a4321.Assemble(0);
+      a4321.Finalize(0);
+
+      const SparseMatrix &A1234 = a1234.SpMat();
+      const SparseMatrix &A4321 = a4321.SpMat();
+
+      SparseMatrix *D = Add(1.0, A1234, -1.0, A4321);
+
+      REQUIRE(D->MaxNorm() == Approx(0.0));
+
+      delete D;
+   }
+}