make style

Disable Build Deploy Container workflow in forks

.github/workflows/build-container.yaml

@@ -20,6 +20,7 @@ on:
 
 jobs:
   build:
+    if: github.repository == 'mfem/mfem' # Don't run in forks
     permissions:
       packages: write
     strategy:

CMakeLists.txt

@@ -556,9 +556,9 @@ endif()
 message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")
 target_include_directories(mfem
   PUBLIC
-    ${TPL_INCLUDE_DIRS}
     $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
-    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>)
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+    ${TPL_INCLUDE_DIRS})
 set_target_properties(mfem PROPERTIES VERSION "${mfem_VERSION}")
 set_target_properties(mfem PROPERTIES SOVERSION "${mfem_VERSION}")
 

examples/ex33.cpp

@@ -3,34 +3,63 @@
 // Compile with: make ex33
 //
 // Sample runs:  ex33 -m ../data/square-disc.mesh -alpha 0.33 -o 2
+//               ex33 -m ../data/square-disc.mesh -alpha 4.5 -o 3
+//               ex33 -m ../data/star.mesh -alpha 1.4 -o 3
 //               ex33 -m ../data/star.mesh -alpha 0.99 -o 3
 //               ex33 -m ../data/inline-quad.mesh -alpha 0.5 -o 3
+//               ex33 -m ../data/amr-quad.mesh -alpha 1.5 -o 3
 //               ex33 -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3
+//               ex33 -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4
 //               ex33 -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
+//               ex33 -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5
+//
+// Verification runs:
+//    ex33 -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2
+//    ex33 -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2
+//    ex33 -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
+//    ex33 -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1
+//
+//  Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i)
+//        for all alpha.
 //
 // Description:
 //
 //  In this example we solve the following fractional PDE with MFEM:
 //
-//    ( - Δ )^α u = f  in Ω,      u = 0  on ∂Ω,      0 < α < 1,
+//    ( - Δ )^α u = f  in Ω,      u = 0  on ∂Ω,      0 < α,
 //
-//  To solve this FPDE, we rely on a rational approximation [2] of the normal
-//  linear operator A^{-α}, where A = - Δ (with associated homogeneous
-//  boundary conditions). Namely, we first approximate the operator
+//  To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer
+//  N is given by floor(α). By doing so, we obtain
 //
-//    A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1},      d_0 = 0,   d_i > 0,
+//    ( - Δ )^(α-N) u = ( - Δ )^(-N) f  in Ω,      u = 0  on ∂Ω,      0 < α.
+//
+//  We first compute the right hand side by solving the integer order PDE
+//
+//   ( - Δ )^N g = f  in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1
+//
+//  The remaining FPDE is then given by
+//
+//  ( - Δ )^(α-N) u = g  in Ω,      u = 0  on ∂Ω.
+//
+//  We rely on a rational approximation [2] of the normal linear operator
+//  A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions)
+//  and (a-N) in (0,1). We approximate the operator
+//
+//    A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1},      d_0 = 0,   d_i > 0,
 //
 //  where I is the L2-identity operator and the coefficients c_i and d_i
 //  are generated offline to a prescribed accuracy in a pre-processing step.
 //  We use the triple-A algorithm [1] to generate the rational approximation
-//  that this partial fractional expansion derives from. We then solve N+1
+//  that this partial fractional expansion derives from. We then solve M+1
 //  independent integer-order PDEs,
 //
-//    A u_i + d_i u_i = c_i f  in Ω,      u_i = 0  on ∂Ω,      i=0,...,N,
+//    A u_i + d_i u_i = c_i g  in Ω,      u_i = 0  on ∂Ω,      i=0,...,M,
 //
 //  using MFEM and sum u_i to arrive at an approximate solution of the FPDE
 //
-//    u ≈ Σ_{i=0}^N u_i.
+//    u ≈ Σ_{i=0}^M u_i.
+//
+//  (If alpha is an integer, we stop after the first PDE was solved.)
 //
 // References:
 //
@@ -47,6 +76,8 @@
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>
+#include <math.h>
+#include <string>
 
 #include "ex33.hpp"
 
@@ -59,8 +90,9 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    int num_refs = 3;
-   bool visualization = true;
    double alpha = 0.5;
+   bool visualization = true;
+   bool verification = false;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -75,6 +107,9 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&verification, "-ver", "--verification", "-no-ver",
+                  "--no-verification",
+                  "Use sinusoidal function (f) for analytic comparison.");
    args.Parse();
    if (!args.Good())
    {
@@ -84,9 +119,31 @@ int main(int argc, char *argv[])
    args.PrintOptions(cout);
 
    Array<double> coeffs, poles;
+   int progress_steps = 1;
 
-   // 2. Compute the coefficients that define the integer-order PDEs.
-   ComputePartialFractionApproximation(alpha,coeffs,poles);
+   // 2. Compute the rational expansion coefficients that define the
+   //    integer-order PDEs.
+   const int power_of_laplace = floor(alpha);
+   double exponent_to_approximate = alpha - power_of_laplace;
+   bool integer_order = false;
+   // Check if alpha is an integer or not.
+   if (abs(exponent_to_approximate) > 1e-12)
+   {
+      mfem::out << "Approximating the fractional exponent "
+                << exponent_to_approximate
+                << endl;
+      ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
+                                          poles);
+
+      // If the example is build without LAPACK, the exponent_to_approximate
+      // might be modified by the function call above.
+      alpha = exponent_to_approximate + power_of_laplace;
+   }
+   else
+   {
+      integer_order = true;
+      mfem::out << "Treating integer order PDE." << endl;
+   }
 
    // 3. Read the mesh from the given mesh file.
    Mesh mesh(mesh_file, 1, 1);
@@ -99,8 +156,8 @@ int main(int argc, char *argv[])
    }
 
    // 5. Define a finite element space on the mesh.
-   FiniteElementCollection *fec = new H1_FECollection(order, dim);
-   FiniteElementSpace fespace(&mesh, fec);
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
    cout << "Number of finite element unknowns: "
         << fespace.GetTrueVSize() << endl;
 
@@ -114,79 +171,234 @@ int main(int argc, char *argv[])
    }
 
    // 7. Define diffusion coefficient, load, and solution GridFunction.
-   ConstantCoefficient f(1.0);
+   auto func = [&alpha](const Vector &x)
+   {
+      double val = 1.0;
+      for (int i=0; i<x.Size(); i++)
+      {
+         val *= sin(M_PI*x(i));
+      }
+      return pow(x.Size()*pow(M_PI,2), alpha) * val;
+   };
+   FunctionCoefficient f(func);
    ConstantCoefficient one(1.0);
    GridFunction u(&fespace);
-   u = 0.;
+   GridFunction x(&fespace);
+   GridFunction g(&fespace);
+   u = 0.0;
+   x = 0.0;
+   g = 0.0;
 
    // 8. Prepare for visualization.
    char vishost[] = "localhost";
    int  visport   = 19916;
-   socketstream xout, uout;
-   ostringstream oss_x, oss_u;
-   if (visualization)
+
+   // 9. Set up the linear form b(.) for integer-order PDE solves.
+   LinearForm b(&fespace);
+   if (verification)
    {
-      xout.open(vishost, visport);
-      xout.precision(8);
-      uout.open(vishost, visport);
-      uout.precision(8);
+      // This statement is only relevant for the verification of the code. It
+      // uses a different f such that an analytic solution is known and easy
+      // to compare with the numerical one. The FPDE becomes:
+      // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
+      // -> u(x,y) = sin(\pi x) sin(\pi y)
+      b.AddDomainIntegrator(new DomainLFIntegrator(f));
    }
-
-   for (int i = 0; i < coeffs.Size(); i++)
+   else
    {
-      // 9. Set up the linear form b(.) for integer-order PDE solve.
-      LinearForm b(&fespace);
-      ProductCoefficient cf(coeffs[i], f);
-      b.AddDomainIntegrator(new DomainLFIntegrator(cf));
-      b.Assemble();
+      b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   }
+   b.Assemble();
 
-      // 10. Define GridFunction for integer-order PDE solve.
-      GridFunction x(&fespace);
-      x = 0.0;
+   // ------------------------------------------------------------------------
+   // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
+   // ------------------------------------------------------------------------
 
-      // 11. Set up the bilinear form a(.,.) for integer-order PDE solve.
-      BilinearForm a(&fespace);
-      a.AddDomainIntegrator(new DiffusionIntegrator(one));
-      ConstantCoefficient c2(-poles[i]);
-      a.AddDomainIntegrator(new MassIntegrator(c2));
-      a.Assemble();
+   if (power_of_laplace > 0)
+   {
+      // 10.1 Compute Stiffnes Matrix
+      BilinearForm k(&fespace);
+      k.AddDomainIntegrator(new DiffusionIntegrator(one));
+      k.Assemble();
 
-      // 12. Assemble the bilinear form and the corresponding linear system.
-      OperatorPtr A;
+      // 10.2 Compute Mass Matrix
+      BilinearForm m(&fespace);
+      m.AddDomainIntegrator(new MassIntegrator(one));
+      m.Assemble();
+      SparseMatrix mass;
+      Array<int> empty;
+      m.FormSystemMatrix(empty, mass);
+
+      // 10.3 Form the system of equations
       Vector B, X;
-      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+      OperatorPtr Op;
+      k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
+      GSSmoother M((SparseMatrix&)(*Op));
 
-      // 13. Solve the linear system A X = B.
-      GSSmoother M((SparseMatrix&)(*A));
-
-      mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
-                << " u = " << coeffs[i] << " f " << endl;
-      PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
-
-      // 14. Recover the solution as a finite element grid function.
-      a.RecoverFEMSolution(X, b, x);
-
-      // 15. Accumulate integer-order PDE solutions.
-      u+=x;
-
-      // 16. Send the solutions by socket to a GLVis server.
-      if (visualization)
+      mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
+                << " ( f ) " << endl;
+      for (int i = 0; i < power_of_laplace; i++)
       {
-         oss_x.str(""); oss_x.clear();
-         oss_x << "Solution of PDE -Δ u + " << -poles[i]
-               << " u = " << coeffs[i] << " f";
-         xout << "solution\n" << mesh << x
-              << "window_title '" << oss_x.str() << "'" << flush;
+         // 10.4 Solve the linear system Op X = B (N times).
+         PCG(*Op, M, B, X, 3, 300, 1e-12, 0.0);
 
-         oss_u.str(""); oss_u.clear();
-         oss_u << "Solution of fractional PDE -Δ^" << alpha
-               << " u = f";
-         uout << "solution\n" << mesh << u
-              << "window_title '" << oss_u.str() << "'" << flush;
+         // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
+         if (i == power_of_laplace - 1)
+         {
+            // Needed for visualization and solution verification.
+            k.RecoverFEMSolution(X, b, g);
+            if (integer_order && verification)
+            {
+               // For an integer order PDE, g is also our solution u.
+               u+=g;
+            }
+            if (visualization)
+            {
+               socketstream fout;
+               ostringstream oss_f;
+               fout.open(vishost, visport);
+               fout.precision(8);
+               oss_f.str(""); oss_f.clear();
+               oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
+                     << power_of_laplace
+                     << " g = f";
+               fout << "solution\n" << mesh << g
+                    << "window_title '" << oss_f.str() << "'" << flush;
+            }
+         }
+
+         // 10.6 Prepare for next iteration (primal / dual space)
+         mass.Mult(X, B);
+         X.SetSubVectorComplement(ess_tdof_list,0.0);
+      }
+
+      // 10.7 Extract solution for the next step. The b now corresponds to the
+      //      function g in the PDE.
+      const SparseMatrix * R = fespace.GetRestrictionMatrix();
+      if (R)
+      {
+         R->MultTranspose(B,b);
+      }
+      else
+      {
+         b = B;
       }
    }
 
-   // 17. Free the used memory.
-   delete fec;
+   // ------------------------------------------------------------------------
+   // 11. Solve the fractional PDE by solving M integer order PDEs and adding
+   //     up the solutions.
+   // ------------------------------------------------------------------------
+   if (!integer_order)
+   {
+      // Setup visualization.
+      socketstream xout, uout;
+      ostringstream oss_x, oss_u;
+      if (visualization)
+      {
+         xout.open(vishost, visport);
+         xout.precision(8);
+         uout.open(vishost, visport);
+         uout.precision(8);
+      }
+      // Iterate over all expansion coefficient that contribute to the
+      // solution.
+      for (int i = 0; i < coeffs.Size(); i++)
+      {
+         mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
+                   << " u = " << coeffs[i] << " g " << endl;
+
+
+         // 11.1 Reset GridFunction for integer-order PDE solve.
+         x = 0.0;
+
+         // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
+         BilinearForm a(&fespace);
+         a.AddDomainIntegrator(new DiffusionIntegrator(one));
+         ConstantCoefficient d_i(-poles[i]);
+         a.AddDomainIntegrator(new MassIntegrator(d_i));
+         a.Assemble();
+
+         // 11.3 Assemble the bilinear form and the corresponding linear system.
+         OperatorPtr A;
+         Vector B, X;
+         a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+         // 11.4 Solve the linear system A X = B.
+         GSSmoother M((SparseMatrix&)(*A));
+
+         PCG(*A, M, B, X, 3, 300, 1e-12, 0.0);
+
+         // 11.5 Recover the solution as a finite element grid function.
+         a.RecoverFEMSolution(X, b, x);
+
+         // 11.6 Accumulate integer-order PDE solutions.
+         x *= coeffs[i];
+         u += x;
+
+         // 11.7 Send fractional PDE solution to a GLVis server.
+         if (visualization)
+         {
+            oss_x.str(""); oss_x.clear();
+            oss_x << "Step " << progress_steps
+                  << ": Solution of PDE -Δ u + " << -poles[i]
+                  << " u = " << coeffs[i] << " g";
+            xout << "solution\n" << mesh << x
+                 << "window_title '" << oss_x.str() << "'" << flush;
+
+            oss_u.str(""); oss_u.clear();
+            oss_u << "Step " << progress_steps + 1
+                  << ": Solution of fractional PDE (-Δ)^" << alpha
+                  << " u = f";
+            uout << "solution\n" << mesh << u
+                 << "window_title '" << oss_u.str() << "'"
+                 << flush;
+         }
+      }
+   }
+
+   // ------------------------------------------------------------------------
+   // 12. (optional) Verify the solution.
+   // ------------------------------------------------------------------------
+   if (verification)
+   {
+      auto solution = [] (const Vector &x)
+      {
+         double val = 1.0;
+         for (int i=0; i<x.Size(); i++)
+         {
+            val *= sin(M_PI*x(i));
+         }
+         return val;
+      };
+      FunctionCoefficient sol(solution);
+      double l2_error = u.ComputeL2Error(sol);
+
+      string analytic_solution,expected_mesh;
+      switch (dim)
+      {
+         case 1:
+            analytic_solution = "sin(π x)";
+            expected_mesh = "inline_segment.mesh";
+            break;
+         case 2:
+            analytic_solution = "sin(π x) sin(π y)";
+            expected_mesh = "inline_quad.mesh";
+            break;
+         default:
+            analytic_solution = "sin(π x) sin(π y) sin(π z)";
+            expected_mesh = "inline_hex.mesh";
+            break;
+      }
+
+      mfem::out << "\n" << string(80,'=')
+                << "\n\nSolution Verification in "<< dim << "D \n\n"
+                << "Analytic solution : " << analytic_solution << "\n"
+                << "Expected mesh     : " << expected_mesh <<"\n"
+                << "Your mesh         : " << mesh_file << "\n"
+                << "L2 error          : " << l2_error << "\n\n"
+                << string(80,'=') << endl;
+   }
+
    return 0;
 }

examples/ex33.hpp

@@ -32,6 +32,7 @@
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>
+#include <string>
 
 using namespace std;
 using namespace mfem;
@@ -249,6 +250,13 @@ void PartialFractionExpansion(double scale, Array<double> & poles,
    coeffs.SetSize(psize);
    coeffs = scale;
 
+   // Note: C p(z)/q(z) = Σ_i c_i / (z - p_i) results in an system of equations
+   // where the N unknowns are the coefficients c_i. After multiplying the
+   // system with q(z), the coefficients c_i can be computed analytically by
+   // choosing N values for z. Choosing z_j = = p_j diagonalizes the system and
+   // one can obtain an analytic form for the c_i coefficients. The result is
+   // implemented in the code block below.
+
    for (int i=0; i<psize; i++)
    {
       double tmp_numer=1.0;
@@ -305,9 +313,12 @@ void ComputePartialFractionApproximation(double & alpha,
    if (print_warning)
    {
       mfem::out
-            << "\nMFEM is compiled without LAPACK.\nUsing precomputed values for PartialFractionApproximation. \n"
-            << "Only alpha = 0.33, 0.5, and 0.99 are available.\nThe default is alpha = 0.5."
-            << std::endl;
+            << "\n" << string(80, '=')
+            << "\nMFEM is compiled without LAPACK."
+            << "\nUsing precomputed values for PartialFractionApproximation."
+            << "\nOnly alpha = 0.33, 0.5, and 0.99 are available."
+            << "\nThe default is alpha = 0.5.\n" << string(80, '=') << "\n"
+            << endl;
    }
    const double eps = std::numeric_limits<double>::epsilon();
 
@@ -351,7 +362,7 @@ void ComputePartialFractionApproximation(double & alpha,
 
    if (print_warning)
    {
-      mfem::out << "Using precomputed values for alpha = "
+      mfem::out << "=> Using precomputed values for alpha = "
                 << alpha << "\n" << std::endl;
    }
 

examples/ex33p.cpp

@@ -3,34 +3,63 @@
 // Compile with: make ex33p
 //
 // Sample runs:  mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 0.33 -o 2
+//               mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 4.5 -o 3
+//               mpirun -np 4 ex33p -m ../data/star.mesh -alpha 1.4 -o 3
 //               mpirun -np 4 ex33p -m ../data/star.mesh -alpha 0.99 -o 3
 //               mpirun -np 4 ex33p -m ../data/inline-quad.mesh -alpha 0.5 -o 3
-//               mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3
+//               mpirun -np 4 ex33p -m ../data/amr-quad.mesh -alpha 1.5 -o 3
+//               mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 -r 2
+//               mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4
 //               mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
+//               mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5
+//
+// Verification runs:
+//    mpirun -np 4 ex33p -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2
+//    mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2
+//    mpirun -np 4 ex33p -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
+//    mpirun -np 4 ex33p -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1
+
+//  Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i)
+//        for all alpha.
 //
 // Description:
 //
 //  In this example we solve the following fractional PDE with MFEM:
 //
-//    ( - Δ )^α u = f  in Ω,      u = 0  on ∂Ω,      0 < α < 1,
+//    ( - Δ )^α u = f  in Ω,      u = 0  on ∂Ω,      0 < α,
 //
-//  To solve this FPDE, we rely on a rational approximation [2] of the normal
-//  linear operator A^{-α}, where A = - Δ (with associated homogeneous
-//  boundary conditions). Namely, we first approximate the operator
+//  To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer
+//  N is given by floor(α). By doing so, we obtain
 //
-//    A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1},      d_0 = 0,   d_i > 0,
+//    ( - Δ )^(α-N) u = ( - Δ )^(-N) f  in Ω,      u = 0  on ∂Ω,      0 < α.
+//
+//  We first compute the right hand side by solving the integer order PDE
+//
+//   ( - Δ )^N g = f  in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1
+//
+//  The remaining FPDE is then given by
+//
+//  ( - Δ )^(α-N) u = g  in Ω,      u = 0  on ∂Ω.
+//
+//  We rely on a rational approximation [2] of the normal linear operator
+//  A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions)
+//  and (a-N) in (0,1). We approximate the operator
+//
+//    A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1},      d_0 = 0,   d_i > 0,
 //
 //  where I is the L2-identity operator and the coefficients c_i and d_i
 //  are generated offline to a prescribed accuracy in a pre-processing step.
 //  We use the triple-A algorithm [1] to generate the rational approximation
-//  that this partial fractional expansion derives from. We then solve N+1
+//  that this partial fractional expansion derives from. We then solve M+1
 //  independent integer-order PDEs,
 //
-//    A u_i + d_i u_i = c_i f  in Ω,      u_i = 0  on ∂Ω,      i=0,...,N,
+//    A u_i + d_i u_i = c_i g  in Ω,      u_i = 0  on ∂Ω,      i=0,...,M,
 //
 //  using MFEM and sum u_i to arrive at an approximate solution of the FPDE
 //
-//    u ≈ Σ_{i=0}^N u_i.
+//    u ≈ Σ_{i=0}^M u_i.
+//
+//  (If alpha is an integer, we stop after the first PDE was solved.)
 //
 // References:
 //
@@ -47,6 +76,8 @@
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>
+#include <math.h>
+#include <string>
 
 #include "ex33.hpp"
 
@@ -65,9 +96,9 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    int num_refs = 3;
-   bool visualization = true;
-   bool visualize_x = false;
    double alpha = 0.5;
+   bool visualization = true;
+   bool verification = false;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -79,12 +110,12 @@ int main(int argc, char *argv[])
                   "Number of uniform refinements");
    args.AddOption(&alpha, "-alpha", "--alpha",
                   "Fractional exponent");
-   args.AddOption(&visualize_x, "-vis_x", "--visualize_x", "-no-vis_x",
-                  "--no-visualization_x",
-                  "Enable or disable GLVis visualization of each integer-order PDE solution.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
-                  "Enable or disable GLVis visualization of the fractional PDE solution.");
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&verification, "-ver", "--verification", "-no-ver",
+                  "--no-verification",
+                  "Use sinusoidal function (f) for analytic comparison.");
    args.Parse();
    if (!args.Good())
    {
@@ -97,61 +128,51 @@ int main(int argc, char *argv[])
    }
 
    Array<double> coeffs, poles;
+   int progress_steps = 1;
 
-   // 2. Compute the coefficients that define the integer-order PDEs.
-   ComputePartialFractionApproximation(alpha,coeffs,poles);
-
-   int num_par_solves;
-   int max_par_solves = max(1,num_procs/2);
-   for (num_par_solves=max_par_solves; num_par_solves>0; num_par_solves--)
+   // 2. Compute the rational expansion coefficients that define the
+   //    integer-order PDEs.
+   const int power_of_laplace = floor(alpha);
+   double exponent_to_approximate = alpha - power_of_laplace;
+   bool integer_order = false;
+   // Check if alpha is an integer or not.
+   if (abs(exponent_to_approximate) > 1e-12)
    {
-      if (num_procs%num_par_solves==0 && num_par_solves<coeffs.Size())
+      if (Mpi::Root())
       {
-         break;
+         mfem::out << "Approximating the fractional exponent "
+                   << exponent_to_approximate
+                   << endl;
+      }
+      ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
+                                          poles);
+
+      // If the example is build without LAPACK, the exponent_to_approximate
+      // might be modified by the function call above.
+      alpha = exponent_to_approximate + power_of_laplace;
+   }
+   else
+   {
+      integer_order = true;
+      if (Mpi::Root())
+      {
+         mfem::out << "Treating integer order PDE." << endl;
       }
    }
-   if (num_par_solves == 1) {num_par_solves = num_procs;}
 
-   int solver_ranks = num_procs/num_par_solves;
-
-   // 3. Split the MPI communicator:
-   //    row_comm is used for parallel partition of the mesh
-   //    col_comm is used for independent integer-order solves
-   int row_color = myid / solver_ranks; // Determine color based on row
-   int col_color = myid % solver_ranks; // Determine color based on col
-
-   MPI_Comm row_comm, col_comm;
-   MPI_Comm_split(MPI_COMM_WORLD, row_color, myid, &row_comm);
-   MPI_Comm_split(MPI_COMM_WORLD, col_color, myid, &col_comm);
-
-   int row_rank, row_size, col_rank, col_size;
-   MPI_Comm_rank(row_comm, &row_rank);
-   MPI_Comm_size(row_comm, &row_size);
-   MPI_Comm_rank(col_comm, &col_rank);
-   MPI_Comm_size(col_comm, &col_size);
-
-   if (Mpi::Root())
-   {
-      mfem::out << "\nTotal number of MPI ranks = " << num_procs << endl;
-      mfem::out << "Number of independent parallel solves = " << col_size << endl;
-      mfem::out << "Number of MPI ranks within each solve = " << row_size
-                <<"\n" << endl;
-   }
-
-   // 4. Read the mesh from the given mesh file.
+   // 3. Read the mesh from the given mesh file.
    Mesh mesh(mesh_file, 1, 1);
    int dim = mesh.Dimension();
 
-   // 5. Refine the mesh to increase the resolution.
+   // 4. Refine the mesh to increase the resolution.
    for (int i = 0; i < num_refs; i++)
    {
       mesh.UniformRefinement();
    }
-
-   ParMesh pmesh(row_comm, mesh);
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
    mesh.Clear();
 
-   // 6. Define a finite element space on the mesh.
+   // 5. Define a finite element space on the mesh.
    H1_FECollection fec(order, dim);
    ParFiniteElementSpace fespace(&pmesh, &fec);
    if (Mpi::Root())
@@ -160,7 +181,7 @@ int main(int argc, char *argv[])
            << fespace.GetTrueVSize() << endl;
    }
 
-   // 7. Determine the list of true (i.e. conforming) essential boundary dofs.
+   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
    Array<int> ess_tdof_list;
    if (pmesh.bdr_attributes.Size())
    {
@@ -169,120 +190,250 @@ int main(int argc, char *argv[])
       fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
 
-   // 8. Define diffusion coefficient, load, and solution GridFunction.
-   ConstantCoefficient f(1.0);
+   // 7. Define diffusion coefficient, load, and solution GridFunction.
+   auto func = [&alpha](const Vector &x)
+   {
+      double val = 1.0;
+      for (int i=0; i<x.Size(); i++)
+      {
+         val *= sin(M_PI*x(i));
+      }
+      return pow(x.Size()*pow(M_PI,2), alpha) * val;
+   };
+   FunctionCoefficient f(func);
    ConstantCoefficient one(1.0);
    ParGridFunction u(&fespace);
    ParGridFunction x(&fespace);
+   ParGridFunction g(&fespace);
    u = 0.0;
+   x = 0.0;
+   g = 0.0;
+
+   // 8. Prepare for visualization.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
 
    // 9. Set up the linear form b(.) for integer-order PDE solves.
    ParLinearForm b(&fespace);
-   b.AddDomainIntegrator(new DomainLFIntegrator(f));
+   if (verification)
+   {
+      // This statement is only relevant for the verification of the code. It
+      // uses a different f such that an analytic solution is known and easy
+      // to compare with the numerical one. The FPDE becomes:
+      // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
+      // -> u(x,y) = sin(\pi x) sin(\pi y)
+      b.AddDomainIntegrator(new DomainLFIntegrator(f));
+   }
+   else
+   {
+      b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   }
    b.Assemble();
 
-   int my_coeff_size = max(coeffs.Size()/col_size,1);
-   int ibeg = col_rank*my_coeff_size;
-   if (ibeg + 2*my_coeff_size > coeffs.Size())
+   // ------------------------------------------------------------------------
+   // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
+   // ------------------------------------------------------------------------
+
+   if (power_of_laplace > 0)
    {
-      my_coeff_size = coeffs.Size()-col_rank*my_coeff_size;
-   }
-   else if (ibeg > coeffs.Size() - 1)
-   {
-      my_coeff_size = 0;
-   }
+      // 10.1 Compute Stiffnes Matrix
+      ParBilinearForm k(&fespace);
+      k.AddDomainIntegrator(new DiffusionIntegrator(one));
+      k.Assemble();
 
-   int iend = ibeg+my_coeff_size;
+      // 10.2 Compute Mass Matrix
+      ParBilinearForm m(&fespace);
+      m.AddDomainIntegrator(new MassIntegrator(one));
+      m.Assemble();
+      HypreParMatrix mass;
+      Array<int> empty;
+      m.FormSystemMatrix(empty, mass);
 
-
-   for (int i = ibeg; i < iend; i++)
-   {
-      // 10. Reset GridFunction for integer-order PDE solve.
-      x = 0.0;
-
-      // 11. Set up the bilinear form a(.,.) for integer-order PDE solve.
-      ParBilinearForm a(&fespace);
-      a.AddDomainIntegrator(new DiffusionIntegrator(one));
-      ConstantCoefficient d_i(-poles[i]);
-      a.AddDomainIntegrator(new MassIntegrator(d_i));
-      a.Assemble();
-
-      // 12. Assemble the bilinear form and the corresponding linear system.
-      OperatorPtr A;
+      // 10.3 Form the system of equations
       Vector B, X;
-      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
-
-      // 13. Solve the linear system A X = B.
-      HypreBoomerAMG * prec = new HypreBoomerAMG;
-      prec->SetPrintLevel(-1);
-
-      int print_level = (col_rank==0) ? 3 : 0;
-      if (Mpi::Root())
-      {
-         mfem::out << "\nMPI rank " <<  myid
-                   << ": Solving PDE -Δ u + " << -poles[i]
-                   << " u = " << coeffs[i] << " f " << endl;
-      }
-      CGSolver cg(row_comm);
+      OperatorPtr Op;
+      k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
+      HypreBoomerAMG prec;
+      prec.SetPrintLevel(-1);
+      CGSolver cg(MPI_COMM_WORLD);
       cg.SetRelTol(1e-12);
       cg.SetMaxIter(2000);
-      cg.SetPrintLevel(print_level);
-      cg.SetPreconditioner(*prec);
-      cg.SetOperator(*A);
-      cg.Mult(B, X);
-      delete prec;
+      cg.SetPrintLevel(3);
+      cg.SetPreconditioner(prec);
+      cg.SetOperator(*Op);
 
-      // 14. Recover the solution as a finite element grid function.
-      a.RecoverFEMSolution(X, b, x);
-
-      // 15. Accumulate integer-order PDE solutions.
-      x *= coeffs[i];
-      u += x;
-
-      // 16. Send integer-order PDE solutions to a GLVis server.
-      if (visualize_x)
+      if (Mpi::Root())
       {
-         if (col_rank > 0 && i < iend-1)
+         mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
+                   << " ( f ) " << endl;
+      }
+      for (int i = 0; i < power_of_laplace; i++)
+      {
+         // 10.4 Solve the linear system Op X = B (N times).
+         cg.Mult(B, X);
+         // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
+         if (i == power_of_laplace - 1)
          {
-            MPI_Status status;
-            MPI_Recv(nullptr,0,MPI_INT, col_rank-1,0,col_comm,&status);
+            // Needed for visualization and solution verification.
+            k.RecoverFEMSolution(X, b, g);
+            if (integer_order && verification)
+            {
+               // For an integer order PDE, g is also our solution u.
+               u+=g;
+            }
+            if (visualization)
+            {
+               socketstream fout;
+               ostringstream oss_f;
+               fout.open(vishost, visport);
+               fout.precision(8);
+               oss_f.str(""); oss_f.clear();
+               oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
+                     << power_of_laplace
+                     << " g = f";
+               fout << "parallel " << num_procs << " " << myid << "\n"
+                    << "solution\n" << pmesh << g
+                    << "window_title '" << oss_f.str() << "'" << flush;
+            }
          }
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream xout(vishost, visport);
+
+         // 10.6 Prepare for next iteration (primal / dual space)
+         mass.Mult(X, B);
+         X.SetSubVectorComplement(ess_tdof_list,0.0);
+      }
+
+      // 10.7 Extract solution for the next step. The b now corresponds to the
+      //      function g in the PDE.
+      const SparseMatrix* rm = fespace.GetRestrictionMatrix();
+      rm->MultTranspose(B, b);
+   }
+
+   // ------------------------------------------------------------------------
+   // 11. Solve the fractional PDE by solving M integer order PDEs and adding
+   //     up the solutions.
+   // ------------------------------------------------------------------------
+   if (!integer_order)
+   {
+      // Setup visualization.
+      socketstream xout, uout;
+      ostringstream oss_x, oss_u;
+      if (visualization)
+      {
+         xout.open(vishost, visport);
          xout.precision(8);
-         ostringstream oss;
-         oss << "Solution of PDE -Δ u + " << -poles[i]
-             << " u = " << coeffs[i] << " f" ;
-         xout << "parallel " << row_size << " " << row_rank << "\n";
-         xout << "solution\n" << pmesh << x
-              << "window_title '" << oss.str() << "'" << flush;
-         if (col_rank < col_size-1)
+         uout.open(vishost, visport);
+         uout.precision(8);
+      }
+      // Iterate over all expansion coefficient that contribute to the
+      // solution.
+      for (int i = 0; i < coeffs.Size(); i++)
+      {
+         if (Mpi::Root())
          {
-            MPI_Send(nullptr,0,MPI_INT,col_rank+1,0,col_comm);
+            mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
+                      << " u = " << coeffs[i] << " g " << endl;
+         }
+
+         // 11.1 Reset GridFunction for integer-order PDE solve.
+         x = 0.0;
+
+         // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
+         ParBilinearForm a(&fespace);
+         a.AddDomainIntegrator(new DiffusionIntegrator(one));
+         ConstantCoefficient d_i(-poles[i]);
+         a.AddDomainIntegrator(new MassIntegrator(d_i));
+         a.Assemble();
+
+         // 11.3 Assemble the bilinear form and the corresponding linear system.
+         OperatorPtr A;
+         Vector B, X;
+         a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+         // 11.4 Solve the linear system A X = B.
+         HypreBoomerAMG prec;
+         prec.SetPrintLevel(-1);
+
+         CGSolver cg(MPI_COMM_WORLD);
+         cg.SetRelTol(1e-12);
+         cg.SetMaxIter(2000);
+         cg.SetPrintLevel(3);
+         cg.SetPreconditioner(prec);
+         cg.SetOperator(*A);
+         cg.Mult(B, X);
+
+         // 11.5 Recover the solution as a finite element grid function.
+         a.RecoverFEMSolution(X, b, x);
+
+         // 11.6 Accumulate integer-order PDE solutions.
+         x *= coeffs[i];
+         u += x;
+
+         // 11.7 Send fractional PDE solution to a GLVis server.
+         if (visualization)
+         {
+            oss_x.str(""); oss_x.clear();
+            oss_x << "Step " << progress_steps
+                  << ": Solution of PDE -Δ u + " << -poles[i]
+                  << " u = " << coeffs[i] << " g";
+            xout << "parallel " << num_procs << " " << myid << "\n"
+                 << "solution\n" << pmesh << x
+                 << "window_title '" << oss_x.str() << "'" << flush;
+
+            oss_u.str(""); oss_u.clear();
+            oss_u << "Step " << progress_steps + 1
+                  << ": Solution of fractional PDE (-Δ)^" << alpha
+                  << " u = f";
+            uout << "parallel " << num_procs << " " << myid << "\n"
+                 << "solution\n" << pmesh << u
+                 << "window_title '" << oss_u.str() << "'"
+                 << flush;
          }
       }
    }
 
-   // 17. Accumulate for the fractional PDE solution
-   MPI_Allreduce(MPI_IN_PLACE, u.GetData(), u.Size(),
-                 MPI_DOUBLE, MPI_SUM,col_comm);
-
-   // 18. Send fractional PDE solution to a GLVis server.
-   if (visualization)
+   // ------------------------------------------------------------------------
+   // 12. (optional) Verify the solution.
+   // ------------------------------------------------------------------------
+   if (verification)
    {
-      if (col_rank == 0)
+      auto solution = [] (const Vector &x)
       {
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream uout(vishost, visport);
-         uout.precision(8);
-         ostringstream oss;
-         oss << "Solution of fractional PDE -Δ^" << alpha
-             << " u = f" ;
-         uout << "parallel " << row_size << " " << row_rank << "\n";
-         uout << "solution\n" << pmesh << u
-              << "window_title '" << oss.str() << "'" << flush;
+         double val = 1.0;
+         for (int i=0; i<x.Size(); i++)
+         {
+            val *= sin(M_PI*x(i));
+         }
+         return val;
+      };
+      FunctionCoefficient sol(solution);
+      double l2_error = u.ComputeL2Error(sol);
+
+      if (Mpi::Root())
+      {
+         string analytic_solution,expected_mesh;
+         switch (dim)
+         {
+            case 1:
+               analytic_solution = "sin(π x)";
+               expected_mesh = "inline_segment.mesh";
+               break;
+            case 2:
+               analytic_solution = "sin(π x) sin(π y)";
+               expected_mesh = "inline_quad.mesh";
+               break;
+            default:
+               analytic_solution = "sin(π x) sin(π y) sin(π z)";
+               expected_mesh = "inline_hex.mesh";
+               break;
+         }
+
+         mfem::out << "\n" << string(80,'=')
+                   << "\n\nSolution Verification in "<< dim << "D \n\n"
+                   << "Analytic solution : " << analytic_solution << "\n"
+                   << "Expected mesh     : " << expected_mesh <<"\n"
+                   << "Your mesh         : " << mesh_file << "\n"
+                   << "L2 error          : " << l2_error << "\n\n"
+                   << string(80,'=') << endl;
       }
    }
 

fem/bilinearform.cpp

@@ -100,6 +100,7 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
 
    // Copy the pointers to the integrators
    domain_integs = bf->domain_integs;
+   domain_integs_marker = bf->domain_integs_marker;
 
    boundary_integs = bf->boundary_integs;
    boundary_integs_marker = bf->boundary_integs_marker;
@@ -112,6 +113,138 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    AllocMat();
 }
 
+BilinearForm::BilinearForm(BilinearForm &&other)
+   : Matrix(other.fes->GetVSize()), mat(other.mat), mat_e(other.mat_e),
+     fes(other.fes), assembly(other.assembly), batch(other.batch),
+     ext(other.ext), sequence(other.sequence), extern_bfs(other.extern_bfs),
+     element_matrices(other.element_matrices), static_cond(other.static_cond),
+     hybridization(other.hybridization), diag_policy(other.diag_policy),
+     precompute_sparsity(other.precompute_sparsity)
+{
+   // We swap stored integrators and markers with the moved nonlinear form
+   mfem::Swap(domain_integs, other.domain_integs);
+   mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+   mfem::Swap(boundary_integs, other.boundary_integs);
+   mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+   mfem::Swap(interior_face_integs, other.interior_face_integs);
+   mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+   mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+   /// Leave the moved nonlinear form in a state as if it was just constructed
+   /// with fes
+   other.sequence = fes->GetSequence();
+   other.mat = nullptr;
+   other.mat_e = nullptr;
+   other.extern_bfs = 0;
+   other.element_matrices = nullptr;
+   other.static_cond = nullptr;
+   other.hybridization = nullptr;
+   other.precompute_sparsity = 0;
+   other.diag_policy = DIAG_KEEP;
+
+   other.assembly = AssemblyLevel::LEGACY;
+   other.batch = 1;
+   other.ext = nullptr;
+}
+
+BilinearForm& BilinearForm::operator=(BilinearForm &&other)
+{
+   if (this != &other)
+   {
+      /// Cleanup current bilinear form first
+      delete mat_e;
+      delete mat;
+      delete element_matrices;
+      delete static_cond;
+      delete hybridization;
+      for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
+      for (int i = 0; i < boundary_integs.Size(); i++) { delete boundary_integs[i]; }
+      for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
+      for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
+      delete ext;
+
+      /// Null out all our integs and set size of their arrays to zero
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         domain_integs[k] = nullptr;
+      }
+      domain_integs.SetSize(0);
+      for (int k = 0; k < boundary_integs.Size(); k++)
+      {
+         boundary_integs[k] = nullptr;
+      }
+      boundary_integs.SetSize(0);
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k] = nullptr;
+      }
+      interior_face_integs.SetSize(0);
+      for (int k = 0; k < boundary_face_integs.Size(); ++k)
+      {
+         boundary_face_integs[k] = nullptr;
+      }
+      boundary_face_integs.SetSize(0);
+
+      /// Null out all our markers and set size of their arrays to zero
+      for (int k = 0; k < domain_integs_marker.Size(); ++k)
+      {
+         domain_integs_marker[k] = nullptr;
+      }
+      domain_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_integs_marker.Size(); ++k)
+      {
+         boundary_integs_marker[k] = nullptr;
+      }
+      boundary_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
+      {
+         boundary_face_integs_marker[k] = nullptr;
+      }
+      boundary_face_integs_marker.SetSize(0);
+
+      /// Now steal data from other bilinear form leaving it in a state as if
+      /// it was just constructed with fes
+      Matrix::operator=(std::move(other));
+
+      mat = other.mat;
+      other.mat = nullptr;
+      mat_e = other.mat_e;
+      other.mat_e = nullptr;
+      fes = other.fes;
+      assembly = other.assembly;
+      other.assembly = AssemblyLevel::LEGACY;
+      batch = other.batch;
+      other.batch = 1;
+      ext = other.ext;
+      other.ext = nullptr;
+      sequence = other.sequence;
+      other.sequence = fes->GetSequence();
+      extern_bfs = other.extern_bfs;
+      other.extern_bfs = 0;
+
+      // Swap our empty integ and marker arrays with the moved bilinear form
+      mfem::Swap(domain_integs, other.domain_integs);
+      mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+      mfem::Swap(boundary_integs, other.boundary_integs);
+      mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+      mfem::Swap(interior_face_integs, other.interior_face_integs);
+      mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+      mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+      element_matrices = other.element_matrices;
+      other.element_matrices = nullptr;
+      static_cond = other.static_cond;
+      other.static_cond = nullptr;
+      hybridization = other.hybridization;
+      other.hybridization = nullptr;
+      diag_policy = other.diag_policy;
+      other.diag_policy = DIAG_KEEP;
+      precompute_sparsity = other.precompute_sparsity;
+      other.precompute_sparsity = 0;
+   }
+   return *this;
+}
+
 void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
 {
    if (ext)

fem/bilinearform.hpp

@@ -90,23 +90,23 @@ protected:
    int extern_bfs;
 
    /// Set of Domain Integrators to be applied.
-   Array<BilinearFormIntegrator*> domain_integs;
+   Array<BilinearFormIntegrator*> domain_integs; // owned
    /// Element attribute marker (should be of length mesh->attributes.Max() or
    /// 0 if mesh->attributes is empty)
    /// Includes all by default.
    /// 0 - ignore attribute
    /// 1 - include attribute
-   Array<Array<int>*>             domain_integs_marker;
+   Array<Array<int>*>             domain_integs_marker; // not owned
 
    /// Set of Boundary Integrators to be applied.
-   Array<BilinearFormIntegrator*> boundary_integs;
+   Array<BilinearFormIntegrator*> boundary_integs; // owned
    Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
 
    /// Set of interior face Integrators to be applied.
-   Array<BilinearFormIntegrator*> interior_face_integs;
+   Array<BilinearFormIntegrator*> interior_face_integs; // owned
 
    /// Set of boundary face Integrators to be applied.
-   Array<BilinearFormIntegrator*> boundary_face_integs;
+   Array<BilinearFormIntegrator*> boundary_face_integs; // owned
    Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
 
    DenseMatrix elemmat;
@@ -141,13 +141,6 @@ protected:
       ext = NULL;
    }
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   BilinearForm(const BilinearForm &);
-
-   /// Copy assignment is not supported; body is undefined.
-   BilinearForm &operator=(const BilinearForm &);
-
 public:
    /// Creates bilinear form associated with FE space @a *f.
    /** The pointer @a f is not owned by the newly constructed object. */
@@ -165,6 +158,20 @@ public:
        #precompute_sparsity, see UsePrecomputedSparsity() for details. */
    BilinearForm(FiniteElementSpace *f, BilinearForm *bf, int ps = 0);
 
+   /// Explicitly prohibit copy construction/assignment of BilinearForms
+   BilinearForm(const BilinearForm &) = delete;
+   BilinearForm &operator=(const BilinearForm &) = delete;
+
+   /// Move constructor for BilinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       BilinearForm. */
+   BilinearForm(BilinearForm &&other);
+
+   /// Move assignment operator for BilinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other BilinearForm. */
+   BilinearForm& operator=(BilinearForm &&other);
+
    /// Get the size of the BilinearForm as a square matrix.
    int Size() const { return height; }
 
@@ -305,38 +312,66 @@ public:
    /// Finalizes the matrix initialization.
    virtual void Finalize(int skip_zeros = 1);
 
-   /// Returns a const reference to the sparse matrix.
+   /** @brief Returns a const reference to the sparse matrix:  \f$ M \f$
+
+       This will fail if HasSpMat() is false. */
    const SparseMatrix &SpMat() const
    {
       MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
       return *mat;
    }
 
-   /// Returns a reference to the sparse matrix:  \f$ M \f$
+   /** @brief Returns a reference to the sparse matrix:  \f$ M \f$
+
+       This will fail if HasSpMat() is false. */
    SparseMatrix &SpMat()
    {
       MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
       return *mat;
    }
 
+   /** @brief Returns true if the sparse matrix is not null, false otherwise.
+
+       @sa SpMat(). */
+   bool HasSpMat()
+   {
+      return mat != nullptr;
+   }
+
+
    /**  @brief Nullifies the internal matrix \f$ M \f$ and returns a pointer
         to it.  Used for transfering ownership. */
    SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; }
 
-   /// Returns a const reference to the sparse matrix of eliminated b.c.: \f$ M_e \f$
+   /** @brief Returns a const reference to the sparse matrix of eliminated b.c.:
+       \f$ M_e \f$
+
+       This will fail if HasSpMatElim() is false. */
    const SparseMatrix &SpMatElim() const
    {
       MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
       return *mat_e;
    }
 
-   /// Returns a reference to the sparse matrix of eliminated b.c.: \f$ M_e \f$
+   /** @brief Returns a reference to the sparse matrix of eliminated b.c.:
+       \f$ M_e \f$
+
+       This will fail if HasSpMatElim() is false. */
    SparseMatrix &SpMatElim()
    {
       MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
       return *mat_e;
    }
 
+   /**  @brief Returns true if the sparse matrix of eliminated b.c.s is not null,
+        false otherwise.
+
+        @sa SpMatElim(). */
+   bool HasSpMatElim()
+   {
+      return mat_e != nullptr;
+   }
+
    /// Adds new Domain Integrator. Assumes ownership of @a bfi.
    void AddDomainIntegrator(BilinearFormIntegrator *bfi);
    /// Adds new Domain Integrator restricted to certain elements specified by

fem/bilininteg_convection_pa.cpp

@@ -1497,6 +1497,7 @@ static void PAConvectionApply(const int dim,
    {
       switch ((D1D << 4 ) | Q1D)
       {
+         case 0x22: return SmemPAConvectionApply3D<2,2>(NE,B,G,Bt,Gt,op,x,y);
          case 0x23: return SmemPAConvectionApply3D<2,3>(NE,B,G,Bt,Gt,op,x,y);
          case 0x24: return SmemPAConvectionApply3D<2,4>(NE,B,G,Bt,Gt,op,x,y);
          case 0x26: return SmemPAConvectionApply3D<2,6>(NE,B,G,Bt,Gt,op,x,y);
@@ -1548,6 +1549,7 @@ static void PAConvectionApplyT(const int dim,
    {
       switch ((D1D << 4 ) | Q1D)
       {
+         case 0x22: return SmemPAConvectionApplyT3D<2,2>(NE,B,G,Bt,Gt,op,x,y);
          case 0x23: return SmemPAConvectionApplyT3D<2,3>(NE,B,G,Bt,Gt,op,x,y);
          case 0x24: return SmemPAConvectionApplyT3D<2,4>(NE,B,G,Bt,Gt,op,x,y);
          case 0x26: return SmemPAConvectionApplyT3D<2,6>(NE,B,G,Bt,Gt,op,x,y);

fem/bilininteg_dgtrace_pa.cpp

@@ -136,6 +136,9 @@ static void PADGTraceSetup(const int dim,
 
 void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
 {
+   const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
+                         Device::GetDeviceMemoryType() : pa_mt;
+
    nf = fes.GetNFbyType(type);
    if (nf==0) { return; }
    // Assumes tensor-product elements
@@ -153,7 +156,7 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    geom = mesh->GetFaceGeometricFactors(
              *ir,
              FaceGeometricFactors::DETERMINANTS |
-             FaceGeometricFactors::NORMALS, type);
+             FaceGeometricFactors::NORMALS, type, mt);
    maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
@@ -695,6 +698,7 @@ static void PADGTraceApply(const int dim,
    {
       switch ((D1D << 4 ) | Q1D)
       {
+         case 0x22: return SmemPADGTraceApply3D<2,2,1>(NF,B,Bt,op,x,y);
          case 0x23: return SmemPADGTraceApply3D<2,3,1>(NF,B,Bt,op,x,y);
          case 0x34: return SmemPADGTraceApply3D<3,4,2>(NF,B,Bt,op,x,y);
          case 0x45: return SmemPADGTraceApply3D<4,5,2>(NF,B,Bt,op,x,y);
@@ -1124,6 +1128,7 @@ static void PADGTraceApplyTranspose(const int dim,
    {
       switch ((D1D << 4 ) | Q1D)
       {
+         case 0x22: return SmemPADGTraceApplyTranspose3D<2,2>(NF,B,Bt,op,x,y);
          case 0x23: return SmemPADGTraceApplyTranspose3D<2,3>(NF,B,Bt,op,x,y);
          case 0x34: return SmemPADGTraceApplyTranspose3D<3,4>(NF,B,Bt,op,x,y);
          case 0x45: return SmemPADGTraceApplyTranspose3D<4,5>(NF,B,Bt,op,x,y);

fem/linearform.cpp

@@ -26,15 +26,134 @@ LinearForm::LinearForm(FiniteElementSpace *f, LinearForm *lf)
    ext = nullptr;
    extern_lfs = 1;
 
-   // Copy the pointers to the integrators
+   // Copy the pointers to the integrators and markers
    domain_integs = lf->domain_integs;
+   domain_integs_marker = lf->domain_integs_marker;
 
    domain_delta_integs = lf->domain_delta_integs;
 
    boundary_integs = lf->boundary_integs;
+   boundary_integs_marker = lf->boundary_integs_marker;
 
    boundary_face_integs = lf->boundary_face_integs;
    boundary_face_integs_marker = lf->boundary_face_integs_marker;
+
+   interior_face_integs = lf->interior_face_integs;
+}
+
+LinearForm::LinearForm(LinearForm &&other)
+   : Vector(std::move(other)), fes(other.fes), ext(other.ext),
+     extern_lfs(other.extern_lfs)
+{
+   // Linear forms are stored on the device
+   UseDevice(true);
+
+   // We swap stored integrators and markers with the moved LinearForm
+   mfem::Swap(domain_integs, other.domain_integs);
+   mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+   mfem::Swap(domain_delta_integs, other.domain_delta_integs);
+   mfem::Swap(boundary_integs, other.boundary_integs);
+   mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+   mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+   mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+   mfem::Swap(interior_face_integs, other.interior_face_integs);
+
+   other.fes = nullptr;
+
+   other.ext = nullptr;
+
+   // moved LinearForm now set to owns its integrators, though none are stored
+   // since we swapped our newly constructed empty Arrays with its Arrays
+   other.extern_lfs = 0;
+}
+
+LinearForm& LinearForm::operator=(LinearForm &&other)
+{
+   if (this != &other)
+   {
+      Vector::operator=(std::move(other));
+
+      /// If we own our integrators, delete them
+      if (extern_lfs == 0)
+      {
+         for (int k = 0; k < domain_delta_integs.Size(); k++)
+         { delete domain_delta_integs[k]; }
+         for (int k = 0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
+         for (int k = 0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
+         for (int k = 0; k < boundary_face_integs.Size(); k++)
+         { delete boundary_face_integs[k]; }
+         for (int k = 0; k < interior_face_integs.Size(); k++)
+         { delete interior_face_integs[k]; }
+      }
+
+      /// Null out all integs
+      for (int k = 0; k < domain_delta_integs.Size(); k++)
+      {
+         domain_delta_integs[k] = nullptr;
+      }
+      domain_delta_integs.SetSize(0);
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         domain_integs[k] = nullptr;
+      }
+      domain_integs.SetSize(0);
+      for (int k = 0; k < boundary_integs.Size(); k++)
+      {
+         boundary_integs[k] = nullptr;
+      }
+      boundary_integs.SetSize(0);
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
+      {
+         boundary_face_integs[k] = nullptr;
+      }
+      boundary_face_integs.SetSize(0);
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k] = nullptr;
+      }
+      interior_face_integs.SetSize(0);
+
+      /// Null out all markers
+      for (int k = 0; k < domain_integs_marker.Size(); ++k)
+      {
+         domain_integs_marker[k] = nullptr;
+      }
+      domain_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_integs_marker.Size(); ++k)
+      {
+         boundary_integs_marker[k] = nullptr;
+      }
+      boundary_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
+      {
+         boundary_face_integs_marker[k] = nullptr;
+      }
+      boundary_face_integs_marker.SetSize(0);
+
+      Vector::operator=(std::move(other));
+
+      // swap stored integrators and markers with the moved LinearForm
+      mfem::Swap(domain_integs, other.domain_integs);
+      mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+      mfem::Swap(domain_delta_integs, other.domain_delta_integs);
+      mfem::Swap(boundary_integs, other.boundary_integs);
+      mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+      mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+      mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+      mfem::Swap(interior_face_integs, other.interior_face_integs);
+
+      fes = other.fes;
+      other.fes = nullptr;
+
+      ext = other.ext;
+      other.ext = nullptr;
+
+      extern_lfs = other.extern_lfs;
+      // moved LinearForm now set to owns its integrators, though none are stored
+      // since we swapped our empty Arrays with its Arrays
+      other.extern_lfs = 0;
+   }
+   return *this;
 }
 
 void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)

fem/linearform.hpp

@@ -74,10 +74,6 @@ protected:
    /// Force (re)computation of delta locations.
    void ResetDeltaLocations() { domain_delta_integs_elem_id.SetSize(0); }
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   LinearForm(const LinearForm &);
-
 public:
    /// Creates linear form associated with FE space @a *f.
    /** The pointer @a f is not owned by the newly constructed object. */
@@ -107,6 +103,9 @@ public:
    LinearForm(FiniteElementSpace *f, double *data) : Vector(data, f->GetVSize())
    { fes = f; ext = nullptr; extern_lfs = 0; }
 
+   /// Explicitly prohibit copy construction of LinearForm
+   LinearForm(const LinearForm &other) = delete;
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
        have the same size.
@@ -116,6 +115,16 @@ public:
    LinearForm &operator=(const LinearForm &rhs)
    { return operator=((const Vector &)rhs); }
 
+   /// Move constructor for LinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       LinearForm. */
+   LinearForm(LinearForm &&other);
+
+   /// Move assignment operator for LinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other LinearForm. */
+   LinearForm& operator=(LinearForm &&other);
+
    /// (DEPRECATED) Return the FE space associated with the LinearForm.
    /** @deprecated Use FESpace() instead. */
    MFEM_DEPRECATED FiniteElementSpace *GetFES() { return fes; }

fem/nonlinearform.cpp

@@ -15,6 +15,92 @@
 namespace mfem
 {
 
+NonlinearForm::NonlinearForm(NonlinearForm &&other)
+   : Operator(other.fes->GetTrueVSize()), assembly(other.assembly),
+     ext(other.ext), fes(other.fes), Grad(other.Grad), cGrad(other.cGrad),
+     sequence(other.fes->GetSequence()), P(other.fes->GetProlongationMatrix()),
+     cP(dynamic_cast<const SparseMatrix*>(P))
+{
+   // We swap stored integrators and markers with the moved nonlinear form
+   mfem::Swap(domain_integs, other.domain_integs);
+   mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+   mfem::Swap(interior_face_integs, other.interior_face_integs);
+   mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+   mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+   /// Leave the moved nonlinear form in a state as if it was just constructed
+   /// with fes
+   other.ext = nullptr;
+   other.cGrad = nullptr;
+   other.Grad = nullptr;
+   other.assembly = AssemblyLevel::LEGACY;
+}
+
+NonlinearForm& NonlinearForm::operator=(NonlinearForm &&other)
+{
+   if (this != &other)
+   {
+      /// Cleanup current nonlinear form first
+      delete cGrad;
+      delete Grad;
+      for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
+      for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
+      for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
+      delete ext;
+
+      /// Null out all our integs and set size of their arrays to zero
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         domain_integs[k] = nullptr;
+      }
+      domain_integs.SetSize(0);
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
+      {
+         boundary_face_integs[k] = nullptr;
+      }
+      boundary_face_integs.SetSize(0);
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k] = nullptr;
+      }
+      interior_face_integs.SetSize(0);
+
+      /// Null out all our markers and set size of their arrays to zero
+      for (int k = 0; k < domain_integs_marker.Size(); ++k)
+      {
+         domain_integs_marker[k] = nullptr;
+      }
+      domain_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
+      {
+         boundary_face_integs_marker[k] = nullptr;
+      }
+      boundary_face_integs_marker.SetSize(0);
+
+      /// Now steal data from other nonlinear form leaving it in a state as if
+      /// it was just constructed with fes
+      Operator::operator=(std::move(other));
+
+      assembly = other.assembly;
+      other.assembly = AssemblyLevel::LEGACY;
+      Grad = other.Grad;
+      other.Grad = nullptr;
+      cGrad = other.cGrad;
+      other.cGrad = nullptr;
+
+      // Swap our empty integ and marker arrays with the moved nonlinear form
+      mfem::Swap(domain_integs, other.domain_integs);
+      mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+      mfem::Swap(interior_face_integs, other.interior_face_integs);
+      mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+      mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+      ext = other.ext;
+      other.ext = nullptr;
+   }
+   return *this;
+}
+
 void NonlinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
 {
    if (ext)
@@ -38,6 +124,38 @@ void NonlinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
    }
 }
 
+void NonlinearForm::AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
+{
+   domain_integs.Append(nlfi);
+   domain_integs_marker.Append(nullptr); // null marker means apply everywhere
+}
+
+void NonlinearForm::AddDomainIntegrator(NonlinearFormIntegrator *nlfi,
+                                        Array<int> &elem_marker)
+{
+   domain_integs.Append(nlfi);
+   domain_integs_marker.Append(&elem_marker);
+}
+
+void NonlinearForm::AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi)
+{
+   interior_face_integs.Append(nlfi);
+}
+
+void NonlinearForm::AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi)
+{
+   boundary_face_integs.Append(nlfi);
+   // null marker means apply everywhere
+   boundary_face_integs_marker.Append(nullptr);
+}
+
+void NonlinearForm::AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi,
+                                         Array<int> &bdr_marker)
+{
+   boundary_face_integs.Append(nlfi);
+   boundary_face_integs_marker.Append(&bdr_marker);
+}
+
 void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
                                    Vector *rhs)
 {
@@ -87,8 +205,10 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
 {
    if (ext)
    {
-      MFEM_VERIFY(!fnfi.Size(), "Interior faces terms not yet implemented!");
-      MFEM_VERIFY(!bfnfi.Size(), "Boundary face terms not yet implemented!");
+      MFEM_VERIFY(!interior_face_integs.Size(),
+                  "Interior faces terms not yet implemented!");
+      MFEM_VERIFY(!boundary_face_integs.Size(),
+                  "Boundary face terms not yet implemented!");
       return ext->GetGridFunctionEnergy(x);
    }
 
@@ -97,30 +217,47 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
    const FiniteElement *fe;
    ElementTransformation *T;
    DofTransformation *doftrans;
+   Mesh *mesh = fes->GetMesh();
    double energy = 0.0;
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != nullptr)
+         {
+            MFEM_VERIFY(mesh->attributes.Size() ==
+                        domain_integs_marker[k]->Size(),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         int elem_attr = mesh->GetAttribute(i);
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
          x.GetSubVector(vdofs, el_x);
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
-         for (int k = 0; k < dnfi.Size(); k++)
+         for (int k = 0; k < domain_integs.Size(); k++)
          {
-            energy += dnfi[k]->GetElementEnergy(*fe, *T, el_x);
+            if (domain_integs_marker[k] == nullptr ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               energy += domain_integs[k]->GetElementEnergy(*fe, *T, el_x);
+            }
          }
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from interior face terms");
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from boundary face terms");
    }
@@ -173,25 +310,41 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
 
    py = 0.0;
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != nullptr)
+         {
+            MFEM_VERIFY(mesh->attributes.Size() ==
+                        domain_integs_marker[k]->Size(),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         int elem_attr = mesh->GetAttribute(i);
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
          px.GetSubVector(vdofs, el_x);
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
-         for (int k = 0; k < dnfi.Size(); k++)
+         for (int k = 0; k < domain_integs.Size(); k++)
          {
-            dnfi[k]->AssembleElementVector(*fe, *T, el_x, el_y);
-            if (doftrans) {doftrans->TransformDual(el_y); }
-            py.AddElementVector(vdofs, el_y);
+            if (domain_integs_marker[k] == nullptr ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               domain_integs[k]->AssembleElementVector(*fe, *T, el_x, el_y);
+               if (doftrans) {doftrans->TransformDual(el_y); }
+               py.AddElementVector(vdofs, el_y);
+            }
          }
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -211,16 +364,16 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
             fe1 = fes->GetFE(tr->Elem1No);
             fe2 = fes->GetFE(tr->Elem2No);
 
-            for (int k = 0; k < fnfi.Size(); k++)
+            for (int k = 0; k < interior_face_integs.Size(); k++)
             {
-               fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+               interior_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
                py.AddElementVector(vdofs, el_y);
             }
          }
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -229,14 +382,14 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); k++)
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -262,12 +415,12 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
             // but we can't dereference a NULL pointer, and we don't want to
             // actually make a fake element.
             fe2 = fe1;
-            for (int k = 0; k < bfnfi.Size(); k++)
+            for (int k = 0; k < boundary_face_integs.Size(); k++)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
 
-               bfnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+               boundary_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
                py.AddElementVector(vdofs, el_y);
             }
          }
@@ -320,26 +473,42 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
       *Grad = 0.0;
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != nullptr)
+         {
+            MFEM_VERIFY(mesh->attributes.Size() ==
+                        domain_integs_marker[k]->Size(),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         int elem_attr = fes->GetMesh()->GetAttribute(i);
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
          px.GetSubVector(vdofs, el_x);
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
-         for (int k = 0; k < dnfi.Size(); k++)
+         for (int k = 0; k < domain_integs.Size(); k++)
          {
-            dnfi[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
-            if (doftrans) { doftrans->TransformDual(elmat); }
-            Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
-            // Grad->AddSubMatrix(vdofs, vdofs, elmat, 1);
+            if (domain_integs_marker[k] == nullptr ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               domain_integs[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
+               if (doftrans) { doftrans->TransformDual(elmat); }
+               Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
+               // Grad->AddSubMatrix(vdofs, vdofs, elmat, 1);
+            }
          }
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -359,16 +528,16 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
             fe1 = fes->GetFE(tr->Elem1No);
             fe2 = fes->GetFE(tr->Elem2No);
 
-            for (int k = 0; k < fnfi.Size(); k++)
+            for (int k = 0; k < interior_face_integs.Size(); k++)
             {
-               fnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
+               interior_face_integs[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
                Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
             }
          }
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -377,14 +546,14 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); k++)
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -410,12 +579,12 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
             // but we can't dereference a NULL pointer, and we don't want to
             // actually make a fake element.
             fe2 = fe1;
-            for (int k = 0; k < bfnfi.Size(); k++)
+            for (int k = 0; k < boundary_face_integs.Size(); k++)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
 
-               bfnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
+               boundary_face_integs[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
                Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
             }
          }
@@ -471,9 +640,9 @@ NonlinearForm::~NonlinearForm()
 {
    delete cGrad;
    delete Grad;
-   for (int i = 0; i <  dnfi.Size(); i++) { delete  dnfi[i]; }
-   for (int i = 0; i <  fnfi.Size(); i++) { delete  fnfi[i]; }
-   for (int i = 0; i < bfnfi.Size(); i++) { delete bfnfi[i]; }
+   for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
+   for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
+   for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
    delete ext;
 }
 
@@ -564,8 +733,8 @@ BlockNonlinearForm::BlockNonlinearForm(Array<FiniteElementSpace *> &f) :
 void BlockNonlinearForm::AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nfi,
                                               Array<int> &bdr_attr_marker)
 {
-   bfnfi.Append(nfi);
-   bfnfi_marker.Append(&bdr_attr_marker);
+   boundary_face_integs.Append(nfi);
+   boundary_face_integs_marker.Append(&bdr_attr_marker);
 }
 
 void BlockNonlinearForm::SetEssentialBC(
@@ -600,7 +769,7 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
       vdofs[i] = new Array<int>;
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
          T = fes[0]->GetElementTransformation(i);
@@ -612,9 +781,9 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
             if (doftrans) {doftrans->InvTransformPrimal(*el_x[s]); }
          }
 
-         for (int k = 0; k < dnfi.Size(); ++k)
+         for (int k = 0; k < domain_integs.Size(); ++k)
          {
-            energy += dnfi[k]->GetElementEnergy(fe, *T, el_x_const);
+            energy += domain_integs[k]->GetElementEnergy(fe, *T, el_x_const);
          }
       }
 
@@ -625,12 +794,12 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
       delete vdofs[i];
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from interior face terms");
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from boundary face terms");
    }
@@ -668,7 +837,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       vdofs2[s] = new Array<int>;
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
@@ -681,10 +850,10 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
             if (doftrans[s]) {doftrans[s]->InvTransformPrimal(*el_x[s]); }
          }
 
-         for (int k = 0; k < dnfi.Size(); ++k)
+         for (int k = 0; k < domain_integs.Size(); ++k)
          {
-            dnfi[k]->AssembleElementVector(fe, *T,
-                                           el_x_const, el_y);
+            domain_integs[k]->AssembleElementVector(fe, *T,
+                                                    el_x_const, el_y);
 
             for (int s=0; s<fes.Size(); ++s)
             {
@@ -696,7 +865,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       Mesh *mesh = fes[0]->GetMesh();
       FaceElementTransformations *tr;
@@ -719,10 +888,10 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
                bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
             }
 
-            for (int k = 0; k < fnfi.Size(); ++k)
+            for (int k = 0; k < interior_face_integs.Size(); ++k)
             {
 
-               fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
+               interior_face_integs[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
 
                for (int s=0; s<fes.Size(); ++s)
                {
@@ -734,7 +903,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       Mesh *mesh = fes[0]->GetMesh();
       FaceElementTransformations *tr;
@@ -742,14 +911,14 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); ++k)
+      for (int k = 0; k < boundary_face_integs.Size(); ++k)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -776,12 +945,12 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
                bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
             }
 
-            for (int k = 0; k < bfnfi.Size(); ++k)
+            for (int k = 0; k < boundary_face_integs.Size(); ++k)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
 
-               bfnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
+               boundary_face_integs[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
 
                for (int s=0; s<fes.Size(); ++s)
                {
@@ -886,7 +1055,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
@@ -899,9 +1068,9 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
             if (doftrans[s]) {doftrans[s]->InvTransformPrimal(*el_x[s]); }
          }
 
-         for (int k = 0; k < dnfi.Size(); ++k)
+         for (int k = 0; k < domain_integs.Size(); ++k)
          {
-            dnfi[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
+            domain_integs[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
 
             for (int j=0; j<fes.Size(); ++j)
             {
@@ -920,7 +1089,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       FaceElementTransformations *tr;
       Mesh *mesh = fes[0]->GetMesh();
@@ -941,9 +1110,9 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
             bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
          }
 
-         for (int k = 0; k < fnfi.Size(); ++k)
+         for (int k = 0; k < interior_face_integs.Size(); ++k)
          {
-            fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
+            interior_face_integs[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
             for (int j=0; j<fes.Size(); ++j)
             {
                for (int l=0; l<fes.Size(); ++l)
@@ -957,7 +1126,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       FaceElementTransformations *tr;
       Mesh *mesh = fes[0]->GetMesh();
@@ -966,14 +1135,14 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); ++k)
+      for (int k = 0; k < boundary_face_integs.Size(); ++k)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -1000,11 +1169,11 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
                bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
             }
 
-            for (int k = 0; k < bfnfi.Size(); ++k)
+            for (int k = 0; k < boundary_face_integs.Size(); ++k)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
-               bfnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
+               boundary_face_integs[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
                for (int l=0; l<fes.Size(); ++l)
                {
                   for (int j=0; j<fes.Size(); ++j)
@@ -1109,19 +1278,19 @@ BlockNonlinearForm::~BlockNonlinearForm()
       delete ess_tdofs[i];
    }
 
-   for (int i = 0; i < dnfi.Size(); ++i)
+   for (int i = 0; i < domain_integs.Size(); ++i)
    {
-      delete dnfi[i];
+      delete domain_integs[i];
    }
 
-   for (int i = 0; i < fnfi.Size(); ++i)
+   for (int i = 0; i < interior_face_integs.Size(); ++i)
    {
-      delete fnfi[i];
+      delete interior_face_integs[i];
    }
 
-   for (int i = 0; i < bfnfi.Size(); ++i)
+   for (int i = 0; i < boundary_face_integs.Size(); ++i)
    {
-      delete bfnfi[i];
+      delete boundary_face_integs[i];
    }
 
 }

fem/nonlinearform.hpp

@@ -36,14 +36,19 @@ protected:
    FiniteElementSpace *fes; // not owned
 
    /// Set of Domain Integrators to be assembled (added).
-   Array<NonlinearFormIntegrator*> dnfi; // owned
+   Array<NonlinearFormIntegrator*> domain_integs; // owned
+   /// Element attribute marker (should be of length mesh->attributes)
+   /// Includes all by default.
+   /// 0 - ignore attribute
+   /// 1 - include attribute
+   Array<Array<int>*>             domain_integs_marker; // not owned
 
    /// Set of interior face Integrators to be assembled (added).
-   Array<NonlinearFormIntegrator*> fnfi; // owned
+   Array<NonlinearFormIntegrator*> interior_face_integs; // owned
 
    /// Set of boundary face Integrators to be assembled (added).
-   Array<NonlinearFormIntegrator*> bfnfi; // owned
-   Array<Array<int>*>              bfnfi_marker; // not owned
+   Array<NonlinearFormIntegrator*> boundary_face_integs; // owned
+   Array<Array<int>*>              boundary_face_integs_marker; // not owned
 
    mutable SparseMatrix *Grad, *cGrad; // owned
    /// Gradient Operator when not assembled as a matrix.
@@ -77,6 +82,20 @@ public:
         cP(dynamic_cast<const SparseMatrix*>(P))
    { }
 
+   /// Explicitly prohibit copy construction/assignment of NonlinearForms
+   NonlinearForm(const NonlinearForm &other) = delete;
+   NonlinearForm& operator=(const NonlinearForm &other) = delete;
+
+   /// Move constructor for NonlinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       NonlinearForm. */
+   NonlinearForm(NonlinearForm &&other);
+
+   /// Move assignment operator for NonlinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other NonlinearForm. */
+   NonlinearForm& operator=(NonlinearForm &&other);
+
    /// Set the desired assembly level. The default is AssemblyLevel::LEGACY.
    /** For nonlinear operators, the "matrix" assembly levels usually do not make
        sense, so only LEGACY, NONE (matrix-free) and PARTIAL are supported.
@@ -106,37 +125,44 @@ public:
    FiniteElementSpace *FESpace() { return fes; }
    const FiniteElementSpace *FESpace() const { return fes; }
 
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
+   /** @brief Adds new domain integrator.
+       @note Assumes ownership of @a nlfi. */
+   void AddDomainIntegrator(NonlinearFormIntegrator *nlfi);
+
+   /** @brief Adds new domain integrator restricted to certain elements
+       specified by the @a elem_marker.
+       @note Assumes ownership of @a nlfi. The array @a elem_marker is stored
+       internally as a pointer to the given Array<int> object. */
+   void AddDomainIntegrator(NonlinearFormIntegrator *nlfi,
+                            Array<int> &elem_marker);
 
    /// Access all integrators added with AddDomainIntegrator().
-   Array<NonlinearFormIntegrator*> *GetDNFI() { return &dnfi; }
-   const Array<NonlinearFormIntegrator*> *GetDNFI() const { return &dnfi; }
+   Array<NonlinearFormIntegrator*> *GetDNFI() { return &domain_integs; }
+   const Array<NonlinearFormIntegrator*> *GetDNFI() const { return &domain_integs; }
 
    /// Adds new Interior Face Integrator.
-   void AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
+   void AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi);
 
    /** @brief Access all interior face integrators added with
        AddInteriorFaceIntegrator(). */
    const Array<NonlinearFormIntegrator*> &GetInteriorFaceIntegrators() const
-   { return fnfi; }
+   { return interior_face_integs; }
 
-   /// Adds new Boundary Face Integrator.
-   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
+   /** @brief Adds new boundary face integrator.
+       @note Assumes ownership of @a nlfi. */
+   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi);
 
-   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
-       attributes. */
-   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nfi,
-                             Array<int> &bdr_marker)
-   { bfnfi.Append(nfi); bfnfi_marker.Append(&bdr_marker); }
+   /** @brief Adds new boundary face integrator restricted to specific boundary
+       attributes by the @a bdr_marker.
+       @note Assumes ownership of @a nlfi. The array @a bdr_marker is stored
+       internally as a pointer to the given Array<int> object. */
+   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi,
+                             Array<int> &bdr_marker);
 
    /** @brief Access all boundary face integrators added with
        AddBdrFaceIntegrator(). */
    const Array<NonlinearFormIntegrator*> &GetBdrFaceIntegrators() const
-   { return bfnfi; }
+   { return boundary_face_integs; }
 
    /// Specify essential boundary conditions.
    /** This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
@@ -226,14 +252,14 @@ protected:
    Array<FiniteElementSpace*> fes;
 
    /// Set of Domain Integrators to be assembled (added).
-   Array<BlockNonlinearFormIntegrator*> dnfi;
+   Array<BlockNonlinearFormIntegrator*> domain_integs;
 
    /// Set of interior face Integrators to be assembled (added).
-   Array<BlockNonlinearFormIntegrator*> fnfi;
+   Array<BlockNonlinearFormIntegrator*> interior_face_integs;
 
    /// Set of Boundary Face Integrators to be assembled (added).
-   Array<BlockNonlinearFormIntegrator*> bfnfi;
-   Array<Array<int>*>           bfnfi_marker;
+   Array<BlockNonlinearFormIntegrator*> boundary_face_integs;
+   Array<Array<int>*>           boundary_face_integs_marker;
 
    /** Auxiliary block-vectors for wrapping input and output vectors or holding
        GridFunction-like block-vector data (e.g. in parallel). */
@@ -298,15 +324,15 @@ public:
 
    /// Adds new Domain Integrator.
    void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
+   { domain_integs.Append(nlfi); }
 
    /// Adds new Interior Face Integrator.
    void AddInteriorFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
+   { interior_face_integs.Append(nlfi); }
 
    /// Adds new Boundary Face Integrator.
    void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
+   { boundary_face_integs.Append(nlfi); boundary_face_integs_marker.Append(NULL); }
 
    /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
        attributes. */

fem/pbilinearform.cpp

@@ -19,6 +19,55 @@
 namespace mfem
 {
 
+ParBilinearForm::ParBilinearForm(ParBilinearForm &&other)
+   : BilinearForm(std::move(other)), pfes(other.pfes),
+     Xaux(other.pfes, other.Xaux.GetData()), Yaux(other.pfes, other.Yaux.GetData()),
+     Ytmp(std::move(other.Ytmp)), p_mat(other.p_mat), p_mat_e(other.p_mat_e),
+     keep_nbr_block(other.keep_nbr_block)
+{
+   other.Xaux.MakeRef(other.pfes, nullptr);
+   other.Yaux.MakeRef(other.pfes, nullptr);
+
+   p_mat.SetOperatorOwner();
+   other.p_mat.SetOperatorOwner(false);
+   other.p_mat.SetType(Operator::Hypre_ParCSR);
+
+   p_mat_e.SetOperatorOwner();
+   other.p_mat_e.SetOperatorOwner(false);
+   other.p_mat_e.SetType(Operator::Hypre_ParCSR);
+
+   other.keep_nbr_block = false;
+}
+
+ParBilinearForm& ParBilinearForm::operator=(ParBilinearForm &&other)
+{
+   if (this != &other)
+   {
+      BilinearForm::operator=(std::move(other));
+      pfes = other.pfes;
+      Xaux.MakeRef(other.pfes, other.Xaux.GetData());
+      other.Xaux.MakeRef(other.pfes, nullptr);
+      Yaux.MakeRef(other.pfes, other.Yaux.GetData());
+      other.Yaux.MakeRef(other.pfes, nullptr);
+
+      Ytmp = std::move(other.Ytmp);
+
+      p_mat = other.p_mat;
+      p_mat.SetOperatorOwner();
+      other.p_mat.SetOperatorOwner(false);
+      other.p_mat.SetType(Operator::Hypre_ParCSR);
+
+      p_mat_e = other.p_mat_e;
+      p_mat_e.SetOperatorOwner();
+      other.p_mat_e.SetOperatorOwner(false);
+      other.p_mat_e.SetType(Operator::Hypre_ParCSR);
+
+      keep_nbr_block = other.keep_nbr_block;
+      other.keep_nbr_block = false;
+   }
+   return *this;
+}
+
 void ParBilinearForm::pAllocMat()
 {
    int nbr_size = pfes->GetFaceNbrVSize();

fem/pbilinearform.hpp

@@ -44,13 +44,6 @@ protected:
 
    void AssembleSharedFaces(int skip_zeros = 1);
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   ParBilinearForm(const ParBilinearForm &);
-
-   /// Copy assignment is not supported; body is undefined.
-   ParBilinearForm &operator=(const ParBilinearForm &);
-
 public:
    /// Creates parallel bilinear form associated with the FE space @a *pf.
    /** The pointer @a pf is not owned by the newly constructed object. */
@@ -71,6 +64,13 @@ public:
         p_mat(Operator::Hypre_ParCSR), p_mat_e(Operator::Hypre_ParCSR)
    { keep_nbr_block = false; }
 
+   /// Explicitly prohibit copy construction/assignment of ParBilinearForm
+   ParBilinearForm(const ParBilinearForm &) = delete;
+   ParBilinearForm &operator=(const ParBilinearForm &) = delete;
+
+   ParBilinearForm(ParBilinearForm &&other);
+   ParBilinearForm& operator=(ParBilinearForm &&other);
+
    /** When set to true and the ParBilinearForm has interior face integrators,
        the local SparseMatrix will include the rows (in addition to the columns)
        corresponding to face-neighbor dofs. The default behavior is to disregard

fem/plinearform.cpp

@@ -18,6 +18,17 @@
 namespace mfem
 {
 
+ParLinearForm& ParLinearForm::operator=(ParLinearForm &&other)
+{
+   if (this != &other)
+   {
+      LinearForm::operator=(std::move(other));
+      pfes = other.pfes;
+      other.pfes = nullptr;
+   }
+   return *this;
+}
+
 void ParLinearForm::Update(ParFiniteElementSpace *pf)
 {
    if (pf) { pfes = pf; }

fem/plinearform.hpp

@@ -28,10 +28,6 @@ class ParLinearForm : public LinearForm
 protected:
    ParFiniteElementSpace *pfes; ///< Points to the same object as #fes
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   ParLinearForm(const ParLinearForm &);
-
 public:
    /** @brief Create an empty ParLinearForm without an associated
        ParFiniteElementSpace.
@@ -64,6 +60,9 @@ public:
    ParLinearForm(ParFiniteElementSpace *pf, ParLinearForm * plf)
       : LinearForm(pf, plf) { pfes = pf; }
 
+   /// Explicitly prohibit copy construction of ParLinearForm
+   ParLinearForm(const ParLinearForm &other) = delete;
+
    /// 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.
@@ -73,6 +72,18 @@ public:
    ParLinearForm &operator=(const ParLinearForm &rhs)
    { return operator=((const Vector &)rhs); }
 
+   /// Move constructor for ParLinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       ParLinearForm. */
+   ParLinearForm(ParLinearForm &&other)
+      : LinearForm(std::move(other)), pfes(other.pfes)
+   { other.pfes = nullptr; }
+
+   /// Move assignment operator for ParLinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other ParLinearForm. */
+   ParLinearForm& operator=(ParLinearForm &&other);
+
    ParFiniteElementSpace *ParFESpace() const { return pfes; }
 
    /// Update the object according to the given new FE space @a *pf.

fem/pnonlinearform.cpp

@@ -27,13 +27,45 @@ ParNonlinearForm::ParNonlinearForm(ParFiniteElementSpace *pf)
    MFEM_VERIFY(!Serial(), "internal MFEM error");
 }
 
+ParNonlinearForm::ParNonlinearForm(ParNonlinearForm &&other)
+   : NonlinearForm(std::move(other)),
+     X((ParFiniteElementSpace *)other.fes, other.X.GetData()),
+     Y((ParFiniteElementSpace *)other.fes, other.Y.GetData()), pGrad(other.pGrad)
+{
+   other.X.MakeRef(other.fes, nullptr);
+   other.Y.MakeRef(other.fes, nullptr);
+
+   pGrad.SetOperatorOwner();
+   other.pGrad.SetOperatorOwner(false);
+   other.pGrad.SetType(Operator::Hypre_ParCSR);
+}
+
+ParNonlinearForm& ParNonlinearForm::operator=(ParNonlinearForm &&other)
+{
+   if (this != &other)
+   {
+      NonlinearForm::operator=(std::move(other));
+      X.MakeRef(other.fes, other.X.GetData());
+      other.X.MakeRef(other.fes, nullptr);
+      Y.MakeRef(other.fes, other.Y.GetData());
+      other.Y.MakeRef(other.fes, nullptr);
+
+      pGrad = other.pGrad;
+      pGrad.SetOperatorOwner();
+      other.pGrad.SetOperatorOwner(false);
+      other.pGrad.SetType(Operator::Hypre_ParCSR);
+
+   }
+   return *this;
+}
+
 double ParNonlinearForm::GetParGridFunctionEnergy(const Vector &x) const
 {
    double loc_energy, glob_energy;
 
    loc_energy = GetGridFunctionEnergy(x);
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from shared faces");
    }
@@ -48,7 +80,7 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
 {
    NonlinearForm::Mult(x, y); // x --(P)--> aux1 --(A_local)--> aux2
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_VERIFY(!NonlinearForm::ext, "Not implemented (extensions + faces");
       // Terms over shared interior faces in parallel.
@@ -78,9 +110,9 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
          X.GetSubVector(vdofs1, el_x.GetData());
          X.FaceNbrData().GetSubVector(vdofs2, el_x.GetData() + vdofs1.Size());
 
-         for (int k = 0; k < fnfi.Size(); k++)
+         for (int k = 0; k < interior_face_integs.Size(); k++)
          {
-            fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+            interior_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
             aux2.AddElementVector(vdofs1, el_y.GetData());
          }
       }
@@ -116,7 +148,7 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 
    OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
 
-   if (fnfi.Size() == 0)
+   if (interior_face_integs.Size() == 0)
    {
       dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
                              pfes->GetDofOffsets(), Grad);
@@ -252,7 +284,7 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
 
    BlockNonlinearForm::MultBlocked(xs, ys);
 
-   if (fnfi.Size() > 0)
+   if (interior_face_integs.Size() > 0)
    {
       MFEM_ABORT("TODO: assemble contributions from shared face terms");
    }
@@ -332,7 +364,7 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
 
    GetLocalGradient(x); // gradients are stored in 'Grads'
 
-   if (fnfi.Size() > 0)
+   if (interior_face_integs.Size() > 0)
    {
       MFEM_ABORT("TODO: assemble contributions from shared face terms");
    }

fem/pnonlinearform.hpp

@@ -32,6 +32,13 @@ protected:
 public:
    ParNonlinearForm(ParFiniteElementSpace *pf);
 
+   /// Explicitly prohibit copy construction/assignment of ParNonlinearForm
+   ParNonlinearForm(const ParNonlinearForm &other) = delete;
+   ParNonlinearForm& operator=(const ParNonlinearForm &other) = delete;
+
+   ParNonlinearForm(ParNonlinearForm &&other);
+   ParNonlinearForm& operator=(ParNonlinearForm &&other);
+
    ParFiniteElementSpace *ParFESpace() const
    { return (ParFiniteElementSpace *)fes; }
 

fem/quadinterpolator_face.cpp

@@ -10,6 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #include "quadinterpolator_face.hpp"
+#include "../general/annotation.hpp"
 #include "../general/forall.hpp"
 #include "../linalg/dtensor.hpp"
 #include "../linalg/kernels.hpp"
@@ -217,9 +218,6 @@ void FaceQuadratureInterpolator::Eval3D(
                "Derivatives on the faces are not yet supported.");
    MFEM_FORALL(f, NF,
    {
-      const int ND1D = T_ND1D ? T_ND1D : nd1d;
-      const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
-      const int VDIM = T_VDIM ? T_VDIM : vdim;
       constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
       constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
       constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
@@ -355,6 +353,184 @@ void FaceQuadratureInterpolator::Eval3D(
    });
 }
 
+template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
+void FaceQuadratureInterpolator::SmemEval3D(
+   const int NF,
+   const int vdim,
+   const DofToQuad &maps,
+   const Array<bool> &signs,
+   const Vector &e_vec,
+   Vector &q_val,
+   Vector &q_der,
+   Vector &q_det,
+   Vector &q_nor,
+   const int eval_flags)
+{
+   MFEM_PERF_SCOPE("FaceQuadInterpolator::SmemEval3D");
+   const int nd1d = maps.ndof;
+   const int nq1d = maps.nqpt;
+   const int ND1D = T_ND1D ? T_ND1D : nd1d;
+   const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
+   const int VDIM = T_VDIM ? T_VDIM : vdim;
+   MFEM_VERIFY(ND1D <= MAX_ND1D, "");
+   MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
+   MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
+   auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
+   auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
+   auto F = Reshape(e_vec.Read(), ND1D, ND1D, VDIM, NF);
+   auto sign = signs.Read();
+   auto val = Reshape(q_val.Write(), NQ1D, NQ1D, VDIM, NF);
+   // auto der = Reshape(q_der.Write(), NQ1D, VDIM, 3, NF);
+   auto det = Reshape(q_det.Write(), NQ1D, NQ1D, NF);
+   auto nor = Reshape(q_nor.Write(), NQ1D, NQ1D, 3, NF);
+   MFEM_VERIFY(eval_flags | DERIVATIVES,
+               "Derivatives on the faces are not yet supported.");
+
+   MFEM_FORALL_3D(f, NF, NQ1D, NQ1D, VDIM,
+   {
+      constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
+      constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
+      constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
+
+      MFEM_SHARED double sm1[max_NQ1D*max_NQ1D*max_VDIM];
+      MFEM_SHARED double sm2[max_NQ1D*max_ND1D*max_VDIM];
+
+      auto s_F = (double(*)[max_ND1D][max_VDIM])sm1;
+      MFEM_FOREACH_THREAD(d1,x,ND1D)
+      {
+         MFEM_FOREACH_THREAD(d2,y,ND1D)
+         {
+            MFEM_FOREACH_THREAD(c,z,VDIM)
+            {
+               s_F[d1][d2][c] = F(d1,d2,c,f);
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      if (eval_flags & VALUES)
+      {
+         auto Bu = (double (*)[max_ND1D][max_VDIM])sm2;
+         MFEM_FOREACH_THREAD(d2,x,ND1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double thrdBu = 0.0;
+                  for (int d1 = 0; d1 < ND1D; ++d1)
+                  {
+                     thrdBu += B(q1,d1)*s_F[d1][d2][c];
+                  }
+                  Bu[q1][d2][c] = thrdBu;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         MFEM_FOREACH_THREAD(q2,x,NQ1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double v = 0.0;
+                  for (int d2 = 0; d2 < ND1D; ++d2)
+                  {
+                     v += B(q2,d2)*Bu[q1][d2][c];
+                  }
+                  val(q1,q2,c,f) = v;
+               }
+            }
+         }
+      }
+
+      if ((eval_flags & DERIVATIVES)
+          || (eval_flags & DETERMINANTS)
+          || (eval_flags & NORMALS))
+      {
+         // We only compute the tangential derivatives
+         auto Gu = (double (*)[max_ND1D][max_VDIM])sm2;
+         MFEM_SHARED double Bu[max_NQ1D][max_ND1D][max_VDIM];
+         MFEM_FOREACH_THREAD(d2,x,ND1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double thrdGu = 0;
+                  double thrdBu = 0;
+                  for (int d1 = 0; d1 < ND1D; ++d1)
+                  {
+                     const double u = s_F[d1][d2][c];
+                     thrdBu += B(q1,d1)*u;
+                     thrdGu += G(q1,d1)*u;
+                  }
+                  Gu[q1][d2][c] = thrdGu;
+                  Bu[q1][d2][c] = thrdBu;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         auto BGu = (double (*)[max_NQ1D][max_VDIM])sm1;
+         MFEM_SHARED double GBu[max_NQ1D][max_NQ1D][max_VDIM];
+         MFEM_FOREACH_THREAD(q2,x,NQ1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double thrdBGu = 0.0;
+                  double thrdGBu = 0.0;
+                  for (int d2 = 0; d2 < ND1D; ++d2)
+                  {
+                     thrdBGu += B(q2,d2)*Gu[q1][d2][c];
+                     thrdGBu += G(q2,d2)*Bu[q1][d2][c];
+                  }
+                  BGu[q2][q1][c] = thrdBGu;
+                  GBu[q2][q1][c] = thrdGBu;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         if (VDIM == 3 && ((eval_flags & NORMALS) ||
+                           (eval_flags & DETERMINANTS)))
+         {
+            double n[3];
+            MFEM_FOREACH_THREAD(q2,x,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(q1,y,NQ1D)
+               {
+                  if (MFEM_THREAD_ID(z) == 0)
+                  {
+                     const double s = sign[f] ? -1.0 : 1.0;
+                     n[0] = s*( BGu[q2][q1][1]*GBu[q2][q1][2]-GBu[q2][q1][1]*
+                                BGu[q2][q1][2] );
+                     n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
+                               BGu[q2][q1][2] );
+                     n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
+                                BGu[q2][q1][1] );
+
+                     const double norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+
+                     if (eval_flags & DETERMINANTS) { det(q1,q2,f) = norm; }
+
+                     if (eval_flags & NORMALS)
+                     {
+                        nor(q1,q2,0,f) = n[0]/norm;
+                        nor(q1,q2,1,f) = n[1]/norm;
+                        nor(q1,q2,2,f) = n[2]/norm;
+                     }
+                  }
+               }
+            }
+         }
+      }
+   });
+}
+
 void FaceQuadratureInterpolator::Mult(
    const Vector &e_vec, unsigned eval_flags,
    Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
@@ -414,21 +590,22 @@ void FaceQuadratureInterpolator::Mult(
          switch (10*nd1d + nq1d)
          {
             // Q0
-            case 11: eval_func = &Eval3D<1,1,1>; break;
-            case 12: eval_func = &Eval3D<1,1,2>; break;
+            case 11: eval_func = &SmemEval3D<1,1,1>; break;
+            case 12: eval_func = &SmemEval3D<1,1,2>; break;
             // Q1
-            case 22: eval_func = &Eval3D<1,2,2>; break;
-            case 23: eval_func = &Eval3D<1,2,3>; break;
+            case 22: eval_func = &SmemEval3D<1,2,2>; break;
+            case 23: eval_func = &SmemEval3D<1,2,3>; break;
+            case 24: eval_func = &SmemEval3D<1,2,4>; break;
             // Q2
-            case 33: eval_func = &Eval3D<1,3,3>; break;
-            case 34: eval_func = &Eval3D<1,3,4>; break;
+            case 33: eval_func = &SmemEval3D<1,3,3>; break;
+            case 34: eval_func = &SmemEval3D<1,3,4>; break;
             // Q3
-            case 44: eval_func = &Eval3D<1,4,4>; break;
-            case 45: eval_func = &Eval3D<1,4,5>; break;
-            case 46: eval_func = &Eval3D<1,4,6>; break;
+            case 44: eval_func = &SmemEval3D<1,4,4>; break;
+            case 45: eval_func = &SmemEval3D<1,4,5>; break;
+            case 46: eval_func = &SmemEval3D<1,4,6>; break;
             // Q4
-            case 55: eval_func = &Eval3D<1,5,5>; break;
-            case 56: eval_func = &Eval3D<1,5,6>; break;
+            case 55: eval_func = &SmemEval3D<1,5,5>; break;
+            case 56: eval_func = &SmemEval3D<1,5,6>; break;
          }
          if (nq1d >= 10 || !eval_func)
          {
@@ -468,18 +645,19 @@ void FaceQuadratureInterpolator::Mult(
          switch (10*nd1d + nq1d)
          {
             // Q1
-            case 22: eval_func = &Eval3D<3,2,2>; break;
-            case 23: eval_func = &Eval3D<3,2,3>; break;
+            case 22: eval_func = &SmemEval3D<3,2,2>; break;
+            case 23: eval_func = &SmemEval3D<3,2,3>; break;
+            case 24: eval_func = &SmemEval3D<3,2,4>; break;
             // Q2
-            case 33: eval_func = &Eval3D<3,3,3>; break;
-            case 34: eval_func = &Eval3D<3,3,4>; break;
+            case 33: eval_func = &SmemEval3D<3,3,3>; break;
+            case 34: eval_func = &SmemEval3D<3,3,4>; break;
             // Q3
-            case 44: eval_func = &Eval3D<3,4,4>; break;
-            case 45: eval_func = &Eval3D<3,4,5>; break;
-            case 46: eval_func = &Eval3D<3,4,6>; break;
+            case 44: eval_func = &SmemEval3D<3,4,4>; break;
+            case 45: eval_func = &SmemEval3D<3,4,5>; break;
+            case 46: eval_func = &SmemEval3D<3,4,6>; break;
             // Q4
-            case 55: eval_func = &Eval3D<3,5,5>; break;
-            case 56: eval_func = &Eval3D<3,5,6>; break;
+            case 55: eval_func = &SmemEval3D<3,5,5>; break;
+            case 56: eval_func = &SmemEval3D<3,5,6>; break;
          }
          if (nq1d >= 10 || !eval_func)
          {

fem/quadinterpolator_face.hpp

@@ -112,6 +112,18 @@ public:
                       Vector &q_det,
                       Vector &q_nor,
                       const int eval_flags);
+
+   template<const int T_VDIM = 0, const int T_ND = 0, const int T_NQ = 0>
+   static void SmemEval3D(const int NF,
+                          const int vdim,
+                          const DofToQuad &maps,
+                          const Array<bool> &signs,
+                          const Vector &e_vec,
+                          Vector &q_val,
+                          Vector &q_der,
+                          Vector &q_det,
+                          Vector &q_nor,
+                          const int eval_flags);
 };
 
 } // mfem namespace

linalg/solvers.cpp

@@ -1903,7 +1903,7 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
       mfem::out << "Newton: Number of iterations: " << final_iter << '\n'
                 << "   ||r|| = " << final_norm << '\n';
    }
-   if (print_options.summary || (!converged && print_options.warnings))
+   if (!converged && (print_options.summary || print_options.warnings))
    {
       mfem::out << "Newton: No convergence!\n";
    }

linalg/vector.hpp

@@ -307,6 +307,12 @@ public:
 
    Vector &operator+=(const Vector &v);
 
+   /// operator- is not supported. Use @ref subtract or @ref Add.
+   Vector &operator-(const Vector &v) = delete;
+
+   /// operator+ is not supported. Use @ref Add.
+   Vector &operator+(const Vector &v) = delete;
+
    /// (*this) += a * Va
    Vector &Add(const double a, const Vector &Va);
 

mesh/mesh.cpp

@@ -859,7 +859,7 @@ const GeometricFactors* Mesh::GetGeometricFactors(const IntegrationRule& ir,
 
 const FaceGeometricFactors* Mesh::GetFaceGeometricFactors(
    const IntegrationRule& ir,
-   const int flags, FaceType type)
+   const int flags, FaceType type, MemoryType d_mt)
 {
    for (int i = 0; i < face_geom_factors.Size(); i++)
    {
@@ -873,7 +873,8 @@ const FaceGeometricFactors* Mesh::GetFaceGeometricFactors(
 
    this->EnsureNodes();
 
-   FaceGeometricFactors *gf = new FaceGeometricFactors(this, ir, flags, type);
+   FaceGeometricFactors *gf = new FaceGeometricFactors(this, ir, flags, type,
+                                                       d_mt);
    face_geom_factors.Append(gf);
    return gf;
 }
@@ -12033,7 +12034,8 @@ void GeometricFactors::Compute(const GridFunction &nodes,
 
 FaceGeometricFactors::FaceGeometricFactors(const Mesh *mesh,
                                            const IntegrationRule &ir,
-                                           int flags, FaceType type)
+                                           int flags, FaceType type,
+                                           MemoryType d_mt)
    : type(type)
 {
    this->mesh = mesh;
@@ -12050,28 +12052,34 @@ FaceGeometricFactors::FaceGeometricFactors(const Mesh *mesh,
                                           ElementDofOrdering::LEXICOGRAPHIC,
                                           type,
                                           L2FaceValues::SingleValued );
-   Vector Fnodes(face_restr->Height());
+
+
+   MemoryType my_d_mt = (d_mt != MemoryType::DEFAULT) ? d_mt :
+                        Device::GetDeviceMemoryType();
+
+   Vector Fnodes(face_restr->Height(), my_d_mt);
    face_restr->Mult(*nodes, Fnodes);
 
    unsigned eval_flags = 0;
+
    if (flags & FaceGeometricFactors::COORDINATES)
    {
-      X.SetSize(vdim*NQ*NF);
+      X.SetSize(vdim*NQ*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::VALUES;
    }
    if (flags & FaceGeometricFactors::JACOBIANS)
    {
-      J.SetSize(vdim*vdim*NQ*NF);
+      J.SetSize(vdim*vdim*NQ*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::DERIVATIVES;
    }
    if (flags & FaceGeometricFactors::DETERMINANTS)
    {
-      detJ.SetSize(NQ*NF);
+      detJ.SetSize(NQ*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::DETERMINANTS;
    }
    if (flags & FaceGeometricFactors::NORMALS)
    {
-      normal.SetSize(vdim*NQ*NF);
+      normal.SetSize(vdim*NQ*NF, my_d_mt);
       eval_flags |= FaceQuadratureInterpolator::NORMALS;
    }
 

mesh/mesh.hpp

@@ -982,7 +982,8 @@ public:
        destructor). */
    const FaceGeometricFactors* GetFaceGeometricFactors(const IntegrationRule& ir,
                                                        const int flags,
-                                                       FaceType type);
+                                                       FaceType type,
+                                                       MemoryType d_mt = MemoryType::DEFAULT);
 
    /// Destroy all GeometricFactors stored by the Mesh.
    /** This method can be used to force recomputation of the GeometricFactors,
@@ -1882,7 +1883,7 @@ public:
    };
 
    FaceGeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags,
-                        FaceType type);
+                        FaceType type, MemoryType d_mt = MemoryType::DEFAULT);
 
    /// Mapped (physical) coordinates of all quadrature points.
    /** This array uses a column-major layout with dimensions (NQ x SDIM x NF)

tests/unit/README.md

@@ -0,0 +1,125 @@
+# MFEM Unit Tests
+
+This directory contains MFEM's suite of unit tests, using the
+[Catch2](https://github.com/catchorg/Catch2) unit testing framework.
+
+## Test executables
+
+MFEM's unit test suite includes a number of executables:
+
+* `unit_tests`
+* `cunit_tests` if MFEM is compiled with CUDA support
+* `sedov_tests_cpu`, `sedov_tests_debug` (and `sedov_tests_cuda` and
+  `sedov_tests_cuda_uvm` if CUDA is enabled), testing a Sedov hydrodynamics case
+* `tmop_pa_tests_cpu`, `tmop_pa_tests_debug` (and `tmop_pa_tests_cuda` if CUDA
+  is enabled), testing TMOP with partial assembly
+
+There are also parallel versions of these executables (prefixed with `p`), which
+are built if MFEM is compiled with MPI enabled.
+
+## Basics of using Catch
+
+To run the unit tests, any of the executables listed above can be run from the
+command line, for example
+```
+./unit_tests
+```
+which will run _all_ serial unit tests. If you want to run only a specific test
+case, it is possible to specify the name of the test case as an argument to the
+test executable (the names of the test cases are given as the first argument to
+the `TEST_CASE` macro in the source code). For example
+```
+./unit_tests "NCMesh PA diagonal"
+```
+will run only test case testing partial assembly of the diagonal on
+non-conforming meshes. Test cases are optionally given one or more _tags_, which
+can be used to group test cases together. Tags can also be specified on the
+command line, for example
+```
+./unit_tests "[NCMesh]"
+```
+will run all the tests that relate to non-conforming meshes.
+
+## Listing tests and tags
+
+It is possible to list all test cases with the `-l` flag, for example
+```
+./unit_tests -l
+```
+will list all serial unit tests, and
+```
+./unit_tests -l "[NCMesh]"
+```
+will list all test cases with the given tag.
+```
+./unit_tests -t
+```
+will list all available tags, along with the number of test cases that are
+assigned to each tag.
+
+## Special tags
+
+For the most part, tags are just used to group similar tests together according
+to their subject matter. However, several specific tags have special meanings,
+and those are:
+
+* `[Parallel]`, which indicates that a test will **not** be tested with the
+  serial test executables, and will only be tested with the parallel executable
+  (e.g. `punit_tests`). `punit_tests` will only run tests marked with
+  `[Parallel]`.
+* `[CUDA]`, which indicates that a test will be tested with the CUDA executables
+  (e.g. `cunit_tests`). These tests will still be run by the standard (CPU)
+  executables. `cunit_tests` will only run tests marked with `[CUDA]`, and its
+  parallel version `pcunit_tests` will only run tests marked with _both_
+  `[CUDA]` and `[Parallel]`.
+* `[MFEMData]`, which indicates that a test requires access to a clone of the
+  MFEM data repository (see the `--data` flag below), in order to run tests on
+  some larger mesh files. By default, tests tagged with this tag are skipped,
+  unless the `--data` flag is provided.
+
+## Special command line arguments
+
+In addition to the standard Catch command line arguments (which can be viewed
+with the `-h` or `--help` flag), MFEM's unit tests support two additional
+command line arguments:
+
+* `--all`, which enables some more thorough tests, at the expense of longer
+  runtimes. This sets the global variable `launch_all_non_regression_tests` to
+  true.
+* `--data`, which specifies a path to a clone of the MFEM [data
+  repository](https://github.com/mfem/data), which contains some larger mesh
+  files. If this argument is provided, tests tagged with `[MFEMData]` will be
+  run, and they will have access to the files in the data repo through the
+  `mfem_data_dir` global variable.
+
+## Test output and debug messages
+
+By default, MFEM's unit tests display relatively little output (a couple of info
+lines at the beginning, and a summary at the end with the number of test cases
+and assertions that were run). If a test fails, some additional information
+about the failing test will be printed. Output to `mfem::out` and `mfem::err` is
+suppressed by default.
+
+To enable more verbose test output, run the unit tests with the `-s` or
+`--success` flag, which will print a message for every successful test
+assertion, including some additional informational messages. With this option,
+output to `mfem::out` and `mfem::err` is enabled.
+
+## Writing unit tests
+
+The following are some guidelines for developers writing unit tests:
+
+* Give your test case a concise yet descriptive name. Whitespace is allowed.
+* Tag your test with the relevant tags. Look at similar tags to see what
+  relevant tags are in use. Class names are often used as tags. Another common
+  tag is `[PartialAssembly]`. See also the section on [special
+  tags](#special-tags).
+* Do not use `std::cout` or `std::cerr` in your test. Prefer the Catch macros
+  [`INFO`](https://github.com/catchorg/Catch2/blob/v2.x/docs/logging.md#top),
+  [`CAPTURE`](https://github.com/catchorg/Catch2/blob/v2.x/docs/logging.md#quickly-capture-value-of-variables-or-expressions),
+  and similar. If you need more control over the output, prefer `mfem::out` and
+  `mfem::err`.
+* Use the
+  [`GENERATE`](https://github.com/catchorg/Catch2/blob/v2.x/docs/generators.md#top)
+  macro instead of nested for-loops when testing many combinations of
+  parameters.

tests/unit/cunit_test_main.cpp

@@ -13,9 +13,6 @@
 #include "mfem.hpp"
 #include "run_unit_tests.hpp"
 
-bool launch_all_non_regression_tests = false;
-std::string mfem_data_dir;
-
 int main(int argc, char *argv[])
 {
    mfem::Device device("cuda");

tests/unit/fem/test_assemblediagonalpa.cpp

@@ -97,8 +97,8 @@ TEST_CASE("Mass Diagonal PA", "[PartialAssembly][AssembleDiagonal]")
    {
       for (int ne = 1; ne < 3; ++ne)
       {
-         std::cout << "Testing " << dimension << "D partial assembly mass diagonal: "
-                   << std::pow(ne, dimension) << " elements." << std::endl;
+         const int n_elements = pow(ne, dimension);
+         CAPTURE(dimension, n_elements);
          for (int order = 1; order < 5; ++order)
          {
             Mesh mesh;
@@ -131,7 +131,7 @@ TEST_CASE("Mass Diagonal PA", "[PartialAssembly][AssembleDiagonal]")
 
             assembly_diag -= pa_diag;
             double error = assembly_diag.Norml2();
-            std::cout << "    order: " << order << ", error norm: " << error << std::endl;
+            CAPTURE(order, error);
             REQUIRE(assembly_diag.Norml2() < 1.e-12);
 
             delete h1_fec;
@@ -146,9 +146,9 @@ TEST_CASE("Diffusion Diagonal PA", "[PartialAssembly][AssembleDiagonal]")
    {
       for (int ne = 1; ne < 3; ++ne)
       {
-         std::cout << "Testing " << dimension <<
-                   "D partial assembly diffusion diagonal: "
-                   << std::pow(ne, dimension) << " elements." << std::endl;
+         const int n_elements = pow(ne, dimension);
+         CAPTURE(dimension, n_elements);
+
          for (int order = 1; order < 5; ++order)
          {
             Mesh mesh;
@@ -224,8 +224,7 @@ TEST_CASE("Diffusion Diagonal PA", "[PartialAssembly][AssembleDiagonal]")
 
                assembly_diag -= pa_diag;
                double error = assembly_diag.Norml2();
-               std::cout << "    order: " << order << ", coefficient type "
-                         << coeffType << ", error norm: " << error << std::endl;
+               CAPTURE(order, coeffType, error);
                REQUIRE(assembly_diag.Norml2() < 1.e-12);
 
                delete coeff;
@@ -364,16 +363,8 @@ TEST_CASE("Hcurl/Hdiv diagonal PA",
             {
                for (int ne = 1; ne < 3; ++ne)
                {
-                  if (spaceType == Hcurl)
-                     std::cout << "Testing " << dimension <<
-                               "D partial assembly H(curl) diagonal for integrator " << integrator
-                               << " and coeffType " << coeffType << ": "
-                               << std::pow(ne, dimension) << " elements." << std::endl;
-                  else
-                     std::cout << "Testing " << dimension <<
-                               "D partial assembly H(div) diagonal for integrator " << integrator
-                               << " and coeffType " << coeffType << ": "
-                               << std::pow(ne, dimension) << " elements." << std::endl;
+                  const int n_elements = std::pow(ne, dimension);
+                  CAPTURE(dimension, spaceType, integrator, coeffType, n_elements);
 
                   int max_order = (dimension == 3) ? 2 : 3;
 
@@ -458,7 +449,7 @@ TEST_CASE("Hcurl/Hdiv diagonal PA",
 
                      assembly_diag -= pa_diag;
                      double error = assembly_diag.Norml2();
-                     std::cout << "    order: " << order << ", error norm: " << error << std::endl;
+                     CAPTURE(order, error);
                      REQUIRE(assembly_diag.Norml2() < 1.e-11);
 
                      delete fec;

tests/unit/fem/test_blocknonlinearform.cpp

@@ -150,7 +150,7 @@ TEST_CASE("ParBlockNonlinearForm",
       // Compute the energy: integral over 1/8 sphere = Pi*1*1*1/6
       double A4 = nf->GetEnergy(x);
 
-      std::cout << "Rank " << my_rank
+      mfem::out << "Rank " << my_rank
                 << ": ParBlockNonlinearForm::GetEnergy = " << A4
                 << ", expected = " << M_PI / 6.0
                 << ", diff = " << (A4 - M_PI / 6.0)

tests/unit/fem/test_calccurlshape.cpp

@@ -105,28 +105,26 @@ void TestCalcCurlShape(FiniteElement* fe, ElementTransformation * T, int res)
    }
 }
 
-TEST_CASE("CalcCurlShape for several ND FiniteElement instances",
+TEST_CASE("CalcCurlShape ND",
           "[ND_TriangleElement]"
           "[ND_QuadrilateralElement]"
           "[ND_TetrahedronElement]"
           "[ND_WedgeElement]"
           "[ND_HexahedronElement]")
 {
-   int maxOrder = 5;
-   int resolution = 10;
+   const int maxOrder = 5;
+   const int resolution = 10;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+
+   CAPTURE(order);
 
    SECTION("ND_TriangleElement")
    {
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TRIANGLE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_TriangleElement::CalcCurlShape() "
-                   << "for order " << order << std::endl;
-         ND_TriangleElement fe(order);
-         TestCalcCurlShape(&fe, &T, resolution);
-      }
+      ND_TriangleElement fe(order);
+      TestCalcCurlShape(&fe, &T, resolution);
    }
 
    SECTION("ND_QuadrilateralElement")
@@ -134,13 +132,8 @@ TEST_CASE("CalcCurlShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::QUADRILATERAL, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_QuadrilateralElement::CalcCurlShape() "
-                   << "for order " << order << std::endl;
-         ND_QuadrilateralElement fe(order);
-         TestCalcCurlShape(&fe, &T, resolution);
-      }
+      ND_QuadrilateralElement fe(order);
+      TestCalcCurlShape(&fe, &T, resolution);
    }
 
    SECTION("ND_TetrahedronElement")
@@ -148,13 +141,8 @@ TEST_CASE("CalcCurlShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TETRAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_TetrahedronElement::CalcCurlShape() "
-                   << "for order " << order << std::endl;
-         ND_TetrahedronElement fe(order);
-         TestCalcCurlShape(&fe, &T, resolution);
-      }
+      ND_TetrahedronElement fe(order);
+      TestCalcCurlShape(&fe, &T, resolution);
    }
 
    SECTION("ND_WedgeElement")
@@ -162,13 +150,8 @@ TEST_CASE("CalcCurlShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::WEDGE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_WedgeElement::CalcCurlShape() "
-                   << "for order " << order << std::endl;
-         ND_WedgeElement fe(order);
-         TestCalcCurlShape(&fe, &T, resolution);
-      }
+      ND_WedgeElement fe(order);
+      TestCalcCurlShape(&fe, &T, resolution);
    }
 
    SECTION("ND_HexahedronElement")
@@ -176,12 +159,7 @@ TEST_CASE("CalcCurlShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::HEXAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_HexahedronElement::CalcCurlShape() "
-                   << "for order " << order << std::endl;
-         ND_HexahedronElement fe(order);
-         TestCalcCurlShape(&fe, &T, resolution);
-      }
+      ND_HexahedronElement fe(order);
+      TestCalcCurlShape(&fe, &T, resolution);
    }
 }

tests/unit/fem/test_calcdivshape.cpp

@@ -95,28 +95,26 @@ void TestCalcDivShape(FiniteElement* fe, ElementTransformation * T, int res)
    }
 }
 
-TEST_CASE("CalcDivShape for several RT FiniteElement instances",
+TEST_CASE("CalcDivShape RT",
           "[RT_TriangleElement]"
           "[RT_QuadrilateralElement]"
           "[RT_TetrahedronElement]"
           "[RT_WedgeElement]"
           "[RT_HexahedronElement]")
 {
-   int maxOrder = 5;
-   int resolution = 10;
+   const int maxOrder = 5;
+   const int resolution = 10;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+
+   CAPTURE(order);
 
    SECTION("RT_TriangleElement")
    {
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TRIANGLE, T);
 
-      for (int order = 1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_TriangleElement::CalcDivShape() "
-                   << "for order " << order << std::endl;
-         RT_TriangleElement fe(order - 1);
-         TestCalcDivShape(&fe, &T, resolution);
-      }
+      RT_TriangleElement fe(order - 1);
+      TestCalcDivShape(&fe, &T, resolution);
    }
 
    SECTION("RT_QuadrilateralElement")
@@ -124,13 +122,8 @@ TEST_CASE("CalcDivShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::QUADRILATERAL, T);
 
-      for (int order = 1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_QuadrilateralElement::CalcDivShape() "
-                   << "for order " << order << std::endl;
-         RT_QuadrilateralElement fe(order - 1);
-         TestCalcDivShape(&fe, &T, resolution);
-      }
+      RT_QuadrilateralElement fe(order - 1);
+      TestCalcDivShape(&fe, &T, resolution);
    }
 
    SECTION("RT_TetrahedronElement")
@@ -138,13 +131,8 @@ TEST_CASE("CalcDivShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TETRAHEDRON, T);
 
-      for (int order = 1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_TetrahedronElement::CalcDivShape() "
-                   << "for order " << order << std::endl;
-         RT_TetrahedronElement fe(order - 1);
-         TestCalcDivShape(&fe, &T, resolution);
-      }
+      RT_TetrahedronElement fe(order - 1);
+      TestCalcDivShape(&fe, &T, resolution);
    }
 
    SECTION("RT_WedgeElement")
@@ -152,13 +140,8 @@ TEST_CASE("CalcDivShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::WEDGE, T);
 
-      for (int order = 1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_WedgeElement::CalcDivShape() "
-                   << "for order " << order << std::endl;
-         RT_WedgeElement fe(order - 1);
-         TestCalcDivShape(&fe, &T, resolution);
-      }
+      RT_WedgeElement fe(order - 1);
+      TestCalcDivShape(&fe, &T, resolution);
    }
 
    SECTION("RT_HexahedronElement")
@@ -166,12 +149,7 @@ TEST_CASE("CalcDivShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::HEXAHEDRON, T);
 
-      for (int order = 1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_HexahedronElement::CalcDivShape() "
-                   << "for order " << order << std::endl;
-         RT_HexahedronElement fe(order - 1);
-         TestCalcDivShape(&fe, &T, resolution);
-      }
+      RT_HexahedronElement fe(order - 1);
+      TestCalcDivShape(&fe, &T, resolution);
    }
 }

tests/unit/fem/test_calcdshape.cpp

@@ -109,7 +109,7 @@ void TestCalcDShape(FiniteElement* fe, ElementTransformation * T, int res)
    }
 }
 
-TEST_CASE("CalcDShape for several H1 FiniteElement instances",
+TEST_CASE("CalcDShape H1",
           "[H1_SegmentElement]"
           "[H1_TriangleElement]"
           "[H1_QuadrilateralElement]"
@@ -117,21 +117,19 @@ TEST_CASE("CalcDShape for several H1 FiniteElement instances",
           "[H1_WedgeElement]"
           "[H1_HexahedronElement]")
 {
-   int maxOrder = 5;
-   int resolution = 10;
+   const int maxOrder = 5;
+   const int resolution = 10;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+
+   CAPTURE(order);
 
    SECTION("H1_SegmentElement")
    {
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::SEGMENT, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_SegmentElement::CalcDShape() "
-                   << "for order " << order << std::endl;
-         H1_SegmentElement fe(order);
-         TestCalcDShape(&fe, &T, resolution);
-      }
+      H1_SegmentElement fe(order);
+      TestCalcDShape(&fe, &T, resolution);
    }
 
    SECTION("H1_TriangleElement")
@@ -139,13 +137,8 @@ TEST_CASE("CalcDShape for several H1 FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TRIANGLE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_TriangleElement::CalcDShape() "
-                   << "for order " << order << std::endl;
-         H1_TriangleElement fe(order);
-         TestCalcDShape(&fe, &T, resolution);
-      }
+      H1_TriangleElement fe(order);
+      TestCalcDShape(&fe, &T, resolution);
    }
 
    SECTION("H1_QuadrilateralElement")
@@ -153,13 +146,8 @@ TEST_CASE("CalcDShape for several H1 FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::QUADRILATERAL, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_QuadrilateralElement::CalcDShape() "
-                   << "for order " << order << std::endl;
-         H1_QuadrilateralElement fe(order);
-         TestCalcDShape(&fe, &T, resolution);
-      }
+      H1_QuadrilateralElement fe(order);
+      TestCalcDShape(&fe, &T, resolution);
    }
 
    SECTION("H1_TetrahedronElement")
@@ -167,13 +155,8 @@ TEST_CASE("CalcDShape for several H1 FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TETRAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_TetrahedronElement::CalcDShape() "
-                   << "for order " << order << std::endl;
-         H1_TetrahedronElement fe(order);
-         TestCalcDShape(&fe, &T, resolution);
-      }
+      H1_TetrahedronElement fe(order);
+      TestCalcDShape(&fe, &T, resolution);
    }
 
    SECTION("H1_WedgeElement")
@@ -181,13 +164,8 @@ TEST_CASE("CalcDShape for several H1 FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::WEDGE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_WedgeElement::CalcDShape() "
-                   << "for order " << order << std::endl;
-         H1_WedgeElement fe(order);
-         TestCalcDShape(&fe, &T, resolution);
-      }
+      H1_WedgeElement fe(order);
+      TestCalcDShape(&fe, &T, resolution);
    }
 
    SECTION("H1_HexahedronElement")
@@ -195,12 +173,7 @@ TEST_CASE("CalcDShape for several H1 FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::HEXAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_HexahedronElement::CalcDShape() "
-                   << "for order " << order << std::endl;
-         H1_HexahedronElement fe(order);
-         TestCalcDShape(&fe, &T, resolution);
-      }
+      H1_HexahedronElement fe(order);
+      TestCalcDShape(&fe, &T, resolution);
    }
 }

tests/unit/fem/test_calcshape.cpp

@@ -108,36 +108,32 @@ void TestCalcShape(FiniteElement* fe, int res, double tol=1e-12)
 }
 
 
-TEST_CASE("CalcShape for several Lagrange FiniteElement instances",
+TEST_CASE("CalcShape Lagrange",
           "[Lagrange1DFiniteElement]"
           "[BiLinear2DFiniteElement]"
           "[BiQuad2DFiniteElement]"
           "[LagrangeHexFiniteElement]")
 {
-   int maxOrder = 5;
-   int resolution = 10;
+   const int maxOrder = 5;
+   const int resolution = 10;
+
 
    SECTION("Lagrange1DFiniteElement")
    {
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing Lagrange1DFiniteElement::CalcShape() "
-                   << "for order " << order << std::endl;
-         Lagrange1DFiniteElement fe(order);
-         TestCalcShape(&fe, resolution);
-      }
+      auto order = GENERATE_COPY(range(1, maxOrder + 1));
+      CAPTURE(order);
+      Lagrange1DFiniteElement fe(order);
+      TestCalcShape(&fe, resolution);
    }
 
    SECTION("BiLinear2DFiniteElement")
    {
-      std::cout << "Testing BiLinear2DFiniteElement::CalcShape()" << std::endl;
       BiLinear2DFiniteElement fe;
       TestCalcShape(&fe, resolution);
    }
 
    SECTION("BiQuad2DFiniteElement")
    {
-      std::cout << "Testing BiQuad2DFiniteElement::CalcShape()" << std::endl;
       BiQuad2DFiniteElement fe;
       TestCalcShape(&fe, resolution);
    }
@@ -145,9 +141,6 @@ TEST_CASE("CalcShape for several Lagrange FiniteElement instances",
 
    SECTION("LagrangeHexFiniteElement")
    {
-      std::cout << "Testing LagrangeHexFiniteElement::CalcShape() "
-                << "for order 2" << std::endl;
-
       // Comments for LagrangeHexFiniteElement state
       // that only degree 2 is functional for this class
       LagrangeHexFiniteElement fe(2);
@@ -155,7 +148,7 @@ TEST_CASE("CalcShape for several Lagrange FiniteElement instances",
    }
 }
 
-TEST_CASE("CalcShape for several H1 FiniteElement instances",
+TEST_CASE("CalcShape H1",
           "[H1_SegmentElement]"
           "[H1_TriangleElement]"
           "[H1_QuadrilateralElement]"
@@ -163,73 +156,46 @@ TEST_CASE("CalcShape for several H1 FiniteElement instances",
           "[H1_HexahedronElement]"
           "[H1_WedgeElement]")
 {
-   int maxOrder = 5;
-   int resolution = 10;
+   const int maxOrder = 5;
+   const int resolution = 10;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+
+   CAPTURE(order);
 
    SECTION("H1_SegmentElement")
    {
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_SegmentElement::CalcShape() "
-                   << "for order " << order << std::endl;
-         H1_SegmentElement fe(order);
-         TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
-      }
+      H1_SegmentElement fe(order);
+      TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
    }
 
    SECTION("H1_TriangleElement")
    {
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_TriangleElement::CalcShape() "
-                   << "for order " << order << std::endl;
-         H1_TriangleElement fe(order);
-         TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
-      }
+      H1_TriangleElement fe(order);
+      TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
    }
 
    SECTION("H1_QuadrilateralElement")
    {
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_QuadrilateralElement::CalcShape() "
-                   << "for order " << order << std::endl;
-         H1_QuadrilateralElement fe(order);
-         TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
-      }
+      H1_QuadrilateralElement fe(order);
+      TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
    }
 
    SECTION("H1_TetrahedronElement")
    {
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_TetrahedronElement::CalcShape() "
-                   << "for order " << order << std::endl;
-         H1_TetrahedronElement fe(order);
-         TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
-      }
+      H1_TetrahedronElement fe(order);
+      TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
    }
 
    SECTION("H1_HexahedronElement")
    {
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_HexahedronElement::CalcShape() "
-                   << "for order " << order << std::endl;
-         H1_HexahedronElement fe(order);
-         TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
-      }
+      H1_HexahedronElement fe(order);
+      TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
    }
 
    SECTION("H1_WedgeElement")
    {
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing H1_WedgeElement::CalcShape() "
-                   << "for order " << order << std::endl;
-         H1_WedgeElement fe(order);
-         TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
-      }
+      H1_WedgeElement fe(order);
+      TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
    }
 
 }

tests/unit/fem/test_calcvshape.cpp

@@ -211,7 +211,7 @@ void TestCalcVShape(FiniteElement* fe, ElementTransformation * T, int res)
    }
 }
 
-TEST_CASE("CalcVShape for several ND FiniteElement instances",
+TEST_CASE("CalcVShape ND",
           "[ND_SegmentElement]"
           "[ND_TriangleElement]"
           "[ND_QuadrilateralElement]"
@@ -219,21 +219,19 @@ TEST_CASE("CalcVShape for several ND FiniteElement instances",
           "[ND_WedgeElement]"
           "[ND_HexahedronElement]")
 {
-   int maxOrder = 5;
-   int resolution = 10;
+   const int maxOrder = 5;
+   const int resolution = 10;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+
+   CAPTURE(order);
 
    SECTION("ND_SegmentElement")
    {
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::SEGMENT, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_SegmentElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         ND_SegmentElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      ND_SegmentElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("ND_TriangleElement")
@@ -241,13 +239,8 @@ TEST_CASE("CalcVShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TRIANGLE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_TriangleElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         ND_TriangleElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      ND_TriangleElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("ND_QuadrilateralElement")
@@ -255,13 +248,8 @@ TEST_CASE("CalcVShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::QUADRILATERAL, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_QuadrilateralElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         ND_QuadrilateralElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      ND_QuadrilateralElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("ND_TetrahedronElement")
@@ -269,13 +257,8 @@ TEST_CASE("CalcVShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TETRAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_TetrahedronElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         ND_TetrahedronElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      ND_TetrahedronElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("ND_WedgeElement")
@@ -283,13 +266,8 @@ TEST_CASE("CalcVShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::WEDGE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_WedgeElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         ND_WedgeElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      ND_WedgeElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("ND_HexahedronElement")
@@ -297,38 +275,31 @@ TEST_CASE("CalcVShape for several ND FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::HEXAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing ND_HexahedronElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         ND_HexahedronElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      ND_HexahedronElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 }
 
-TEST_CASE("CalcVShape for several RT FiniteElement instances",
+TEST_CASE("CalcVShape RT",
           "[RT_TriangleElement]"
           "[RT_QuadrilateralElement]"
           "[RT_TetrahedronElement]"
           "[RT_WedgeElement]"
           "[RT_HexahedronElement]")
 {
-   int maxOrder = 5;
-   int resolution = 10;
+   const int maxOrder = 5;
+   const int resolution = 10;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+
+   CAPTURE(order);
 
    SECTION("RT_TriangleElement")
    {
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TRIANGLE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_TriangleElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         RT_TriangleElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      RT_TriangleElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("RT_QuadrilateralElement")
@@ -336,13 +307,8 @@ TEST_CASE("CalcVShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::QUADRILATERAL, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_QuadrilateralElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         RT_QuadrilateralElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      RT_QuadrilateralElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("RT_TetrahedronElement")
@@ -350,13 +316,8 @@ TEST_CASE("CalcVShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::TETRAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_TetrahedronElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         RT_TetrahedronElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      RT_TetrahedronElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("RT_WedgeElement")
@@ -364,13 +325,8 @@ TEST_CASE("CalcVShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::WEDGE, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_WedgeElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         RT_WedgeElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      RT_WedgeElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 
    SECTION("RT_HexahedronElement")
@@ -378,12 +334,7 @@ TEST_CASE("CalcVShape for several RT FiniteElement instances",
       IsoparametricTransformation T;
       GetReferenceTransformation(Element::HEXAHEDRON, T);
 
-      for (int order =1; order <= maxOrder; ++order)
-      {
-         std::cout << "Testing RT_HexahedronElement::CalcVShape() "
-                   << "for order " << order << std::endl;
-         RT_HexahedronElement fe(order);
-         TestCalcVShape(&fe, &T, resolution);
-      }
+      RT_HexahedronElement fe(order);
+      TestCalcVShape(&fe, &T, resolution);
    }
 }

tests/unit/fem/test_datacollection.cpp

@@ -27,7 +27,6 @@ TEST_CASE("Save and load from collections", "[DataCollection]")
 {
    SECTION("VisIt data files")
    {
-      std::cout<<"Testing VisIt data files"<<std::endl;
       // Set up a small mesh and a couple of grid function on that mesh
       Mesh mesh = Mesh::MakeCartesian2D(2, 3, Element::QUADRILATERAL, 0, 2.0, 3.0);
       FiniteElementCollection *fec = new LinearFECollection;
@@ -59,8 +58,6 @@ TEST_CASE("Save and load from collections", "[DataCollection]")
 
       SECTION("Uncompressed MFEM format")
       {
-         std::cout<<"Testing uncompressed MFEM format"<<std::endl;
-
          // Collect the mesh and grid functions into a DataCollection and test that they got in there
          VisItDataCollection dc("base", &mesh);
          dc.RegisterField("u", u);
@@ -147,8 +144,6 @@ TEST_CASE("Save and load from collections", "[DataCollection]")
 #ifdef MFEM_USE_ZLIB
       SECTION("Compressed MFEM format")
       {
-         std::cout<<"Testing compressed MFEM format"<<std::endl;
-
          // Collect the mesh and grid functions into a DataCollection and test that they got in there
          VisItDataCollection dc("base", &mesh);
          dc.RegisterField("u", u);

tests/unit/fem/test_eigs.cpp

@@ -78,7 +78,7 @@ TEST_CASE("Laplacian Eigenvalues",
       H1_FECollection fec(order, dim);
       FiniteElementSpace fespace(mesh, &fec);
       int size = fespace.GetTrueVSize();
-      std::cout << mt << " Eigenvalue system size: " << size << std::endl;
+      CAPTURE(mt, size);
 
       Array<int> ess_bdr;
       if (mesh->bdr_attributes.Size())
@@ -149,8 +149,7 @@ TEST_CASE("Laplacian Eigenvalues",
          max_err = std::max(max_err, err);
          REQUIRE(err < 5.0);
       }
-      std::cout << mt << " Maximum relative error: " << max_err << "%"
-                << std::endl;
+      CAPTURE(mt, max_err);
 
       delete mesh;
    }
@@ -199,10 +198,7 @@ TEST_CASE("Laplacian Eigenvalues in Parallel",
       H1_FECollection fec(order, dim);
       ParFiniteElementSpace fespace(&pmesh, &fec);
       HYPRE_Int size = fespace.GlobalTrueVSize();
-      if (my_rank == 0)
-      {
-         std::cout << mt << " Eigenvalue system size: " << size << std::endl;
-      }
+      CAPTURE(mt, size);
 
       Array<int> ess_bdr;
       if (pmesh.bdr_attributes.Size())
@@ -260,11 +256,7 @@ TEST_CASE("Laplacian Eigenvalues in Parallel",
          max_err = std::max(max_err, err);
          REQUIRE(err < 5.0);
       }
-      if (my_rank == 0)
-      {
-         std::cout << mt << " Maximum relative error: " << max_err << "%"
-                   << std::endl;
-      }
+      CAPTURE(mt, max_err);
 
       delete A;
       delete M;

tests/unit/fem/test_face_elem_trans.cpp

@@ -20,7 +20,6 @@ namespace face_elem_trans
 TEST_CASE("3D FaceElementTransformations",
           "[FaceElementTransformations]")
 {
-   int log = 0;
    int n = 1;
    int order = 1;
 
@@ -32,9 +31,9 @@ TEST_CASE("3D FaceElementTransformations",
       int npts = 0;
       for (int f=0; f<mesh.GetNFaces(); f++)
       {
-         if (log > 0)
+         if (verbose_tests)
          {
-            std::cout << "Getting trans for face " << f << std::endl;
+            mfem::out << "Getting trans for face " << f << std::endl;
          }
          FaceElementTransformations *T =
             mesh.GetInteriorFaceTransformations(f);
@@ -43,9 +42,9 @@ TEST_CASE("3D FaceElementTransformations",
          {
             const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
                                                      2*order + 2);
-            if (log > 0)
+            if (verbose_tests)
             {
-               std::cout << f << " " << T->Elem1No
+               mfem::out << f << " " << T->Elem1No
                          << " " << T->Elem2No << std::endl;
             }
 
@@ -83,9 +82,9 @@ TEST_CASE("3D FaceElementTransformations",
                }
             }
          }
-         if (log > 0)
+         if (verbose_tests)
          {
-            std::cout << "Checked " << npts << " points within face "
+            mfem::out << "Checked " << npts << " points within face "
                       << f << std::endl;
          }
       }

tests/unit/fem/test_face_permutation.cpp

@@ -219,7 +219,6 @@ TEST_CASE("2D Face Permutation", "[Face Permutation]")
          delete mesh;
       }
    }
-   std::cout << "2D Face Permutation: max_err = " << max_err << '\n';
    REQUIRE(max_err < 1e-15);
 }
 
@@ -237,6 +236,5 @@ TEST_CASE("3D Face Permutation", "[Face Permutation]")
          delete mesh;
       }
    }
-   std::cout << "3D Face Permutation: max_err = " << max_err << '\n';
    REQUIRE(max_err < 1e-15);
 }

tests/unit/fem/test_inversetransform.cpp

@@ -113,8 +113,7 @@ TEST_CASE("InverseElementTransformation",
             max_err = std::max(max_err, std::abs(ipRev.y - ip.y));
          }
       }
-      std::cout << "Points found: " << pts_found << '/' << npts << '\n'
-                << "Maximum error: " << max_err << '\n';
+      CAPTURE(pts_found, npts, max_err);
       REQUIRE( pts_found == npts );
       REQUIRE( max_err <= tol );
    }
@@ -163,8 +162,7 @@ TEST_CASE("InverseElementTransformation",
             max_err = std::max(max_err, std::abs(ipRev.y - ip.y));
          }
       }
-      std::cout << "Points found: " << pts_found << '/' << npts << '\n'
-                << "Maximum error: " << max_err << '\n';
+      CAPTURE(pts_found, npts, max_err);
       REQUIRE( pts_found >= min_found_pts );
       REQUIRE( max_err <= tol );
    }

tests/unit/fem/test_operatorjacobismoother.cpp

@@ -14,19 +14,16 @@
 
 using namespace mfem;
 
-namespace operatorjacobismoother
-{
-
-TEST_CASE("operatorjacobismoother")
+TEST_CASE("OperatorJacobiSmoother", "[OperatorJacobiSmoother]")
 {
    for (int dimension = 2; dimension < 4; ++dimension)
    {
       for (int ne = 1; ne < 3; ++ne)
       {
-         std::cout << "Testing " << dimension << "D partial assembly smoother: "
-                   << std::pow(ne, dimension) << " elements." << std::endl;
+         const int n_elements = std::pow(ne, dimension);
          for (int order = 1; order < 5; ++order)
          {
+            CAPTURE(dimension, n_elements, order);
             Mesh mesh;
             if (dimension == 2)
             {
@@ -78,8 +75,6 @@ TEST_CASE("operatorjacobismoother")
             pa_smoother.Mult(xin, y_pa);
 
             y_fa -= y_pa;
-            double error = y_fa.Norml2();
-            std::cout << "    order: " << order << ", error norm: " << error << std::endl;
             REQUIRE(y_fa.Norml2() < 1.e-12);
 
             delete h1_fec;
@@ -88,15 +83,14 @@ TEST_CASE("operatorjacobismoother")
    }
 }
 
-TEST_CASE("operatorjacobifichera")
+TEST_CASE("OperatorJacobiSmoother Fichera", "[OperatorJacobiSmoother]")
 {
    const int dimension = 3;
    for (int refine = 1; refine < 4; ++refine)
    {
-      std::cout << "Testing " << 3 << "D partial assembly smoother: "
-                << "fichera mesh, refine level " << refine << std::endl;
       for (int order = 1; order < 5; ++order)
       {
+         CAPTURE(refine, order);
          Mesh mesh("../../data/fichera.mesh", 1, refine, true);
          FiniteElementCollection *h1_fec = new H1_FECollection(order, dimension);
          FiniteElementSpace h1_fespace(&mesh, h1_fec);
@@ -138,13 +132,9 @@ TEST_CASE("operatorjacobifichera")
          pa_smoother.Mult(xin, y_pa);
 
          y_fa -= y_pa;
-         double error = y_fa.Norml2();
-         std::cout << "    order: " << order << ", error norm: " << error << std::endl;
          REQUIRE(y_fa.Norml2() < 1.e-12);
 
          delete h1_fec;
       }
    }
 }
-
-} // namespace operatorjacobismoother

tests/unit/fem/test_pa_coeff.cpp

@@ -135,11 +135,9 @@ TEST_CASE("H1 PA Coefficient", "[PartialAssembly][Coefficient]")
          for (int integrator = 0; integrator < 2; ++integrator)
          {
             const int ne = 2;
-            std::cout << "Testing " << dimension << "D partial assembly with "
-                      << "coeffType " << coeffType << " and "
-                      << "integrator " << integrator << std::endl;
             for (int order = 1; order < 4; ++order)
             {
+               CAPTURE(dimension, coeffType, integrator, order);
                Mesh mesh = MakeCartesianNonaligned(dimension, ne);
 
                FiniteElementCollection* h1_fec =
@@ -251,15 +249,10 @@ TEST_CASE("H1 PA Coefficient", "[PartialAssembly][Coefficient]")
 
                y_pa -= y_mat;
                double pa_error = y_pa.Norml2();
-               std::cout << "  order: " << order
-                         << ", pa error norm: " << pa_error << std::endl;
                REQUIRE(pa_error < 1.e-12);
 
                y_assembly -= y_mat;
                double assembly_error = y_assembly.Norml2();
-               std::cout << "  order: " << order
-                         << ", assembly error norm: " << assembly_error
-                         << std::endl;
                REQUIRE(assembly_error < 1.e-12);
 
                delete coeff;
@@ -329,29 +322,10 @@ TEST_CASE("Hcurl/Hdiv PA Coefficient",
 
             for (int integrator = 0; integrator < numIntegrators; ++integrator)
             {
-               if (spaceType == Hcurl)
-                  std::cout << "Testing " << dimension
-                            << "D ND partial assembly with coeffType "
-                            << coeffType << " and integrator "
-                            << integrator << std::endl;
-               else if (spaceType == Hdiv)
-                  std::cout << "Testing " << dimension
-                            << "D RT partial assembly with coeffType "
-                            << coeffType << " and integrator "
-                            << integrator << std::endl;
-               else if (spaceType == HcurlHdiv)
-                  std::cout << "Testing " << dimension
-                            << "D ND x RT partial assembly with coeffType "
-                            << coeffType << " and integrator "
-                            << integrator << std::endl;
-               else  // HdivHcurl
-                  std::cout << "Testing " << dimension
-                            << "D RT x ND partial assembly with coeffType "
-                            << coeffType << " and integrator "
-                            << integrator << std::endl;
-
                for (int order = 1; order < 4; ++order)
                {
+                  CAPTURE(spaceType, dimension, coeffType, integrator, order);
+
                   FiniteElementCollection* fec = nullptr;
                   if (spaceType == Hcurl || spaceType == HcurlHdiv)
                   {
@@ -518,15 +492,10 @@ TEST_CASE("Hcurl/Hdiv PA Coefficient",
 
                         v_pa -= v_mat;
                         double pa_error = v_pa.Norml2();
-                        std::cout << "  order: " << order
-                                  << ", pa transpose error norm: " << pa_error << std::endl;
                         REQUIRE(pa_error < 1.e-12);
 
                         v_assembly -= v_mat;
                         double assembly_error = v_assembly.Norml2();
-                        std::cout << "  order: " << order
-                                  << ", assembly transpose error norm: " << assembly_error
-                                  << std::endl;
                         REQUIRE(assembly_error < 1.e-12);
                      }
 
@@ -615,15 +584,10 @@ TEST_CASE("Hcurl/Hdiv PA Coefficient",
 
                   y_pa -= y_mat;
                   double pa_error = y_pa.Norml2();
-                  std::cout << "  order: " << order
-                            << ", pa error norm: " << pa_error << std::endl;
                   REQUIRE(pa_error < 1.e-10);
 
                   y_assembly -= y_mat;
                   double assembly_error = y_assembly.Norml2();
-                  std::cout << "  order: " << order
-                            << ", assembly error norm: " << assembly_error
-                            << std::endl;
                   REQUIRE(assembly_error < 1.e-12);
 
                   delete fec;
@@ -687,21 +651,9 @@ TEST_CASE("Hcurl/Hdiv Mixed PA Coefficient",
             const int numIntegrators = (spaceType == HcurlL2 && dimension == 3) ? 2 : 1;
             for (int integrator = 0; integrator < numIntegrators; ++integrator)
             {
-               if (spaceType == HcurlH1)
-                  std::cout << "Testing " << dimension << "D ND H1 mixed partial assembly with "
-                            << "coeffType " << coeffType << " and "
-                            << "integrator " << integrator << std::endl;
-               else if (spaceType == HcurlL2)
-                  std::cout << "Testing " << dimension << "D ND L2 mixed partial assembly with "
-                            << "coeffType " << coeffType << " and "
-                            << "integrator " << integrator << std::endl;
-               else
-                  std::cout << "Testing " << dimension << "D RT L2 mixed partial assembly with "
-                            << "coeffType " << coeffType << " and "
-                            << "integrator " << integrator << std::endl;
-
                for (int order = 1; order < 4; ++order)
                {
+                  CAPTURE(spaceType, dimension, coeffType, integrator, order);
                   FiniteElementCollection* vec_fec = nullptr;
                   if (spaceType == HcurlH1 || spaceType == HcurlL2 || spaceType == HcurlH1_2D)
                   {
@@ -818,15 +770,10 @@ TEST_CASE("Hcurl/Hdiv Mixed PA Coefficient",
 
                   y_pa -= y_mat;
                   double pa_error = y_pa.Norml2();
-                  std::cout << "  order: " << order
-                            << ", pa error norm: " << pa_error << std::endl;
                   REQUIRE(pa_error < 1.e-12);
 
                   y_assembly -= y_mat;
                   double assembly_error = y_assembly.Norml2();
-                  std::cout << "  order: " << order
-                            << ", assembly error norm: " << assembly_error
-                            << std::endl;
                   REQUIRE(assembly_error < 1.e-12);
 
                   delete xin;
@@ -852,15 +799,10 @@ TEST_CASE("Hcurl/Hdiv Mixed PA Coefficient",
 
                      y_pa -= y_mat;
                      pa_error = y_pa.Norml2();
-                     std::cout << "  order: " << order
-                               << ", pa transpose error norm: " << pa_error << std::endl;
                      REQUIRE(pa_error < 1.e-12);
 
                      y_assembly -= y_mat;
                      assembly_error = y_assembly.Norml2();
-                     std::cout << "  order: " << order
-                               << ", assembly transpose error norm: " << assembly_error
-                               << std::endl;
                      REQUIRE(assembly_error < 1.e-12);
                   }
 

tests/unit/fem/test_project_bdr.cpp

@@ -29,7 +29,6 @@ TEST_CASE("3D ProjectBdrCoefficientTangent",
           "[GridFunction]"
           "[VectorGridFunctionCoefficient]")
 {
-   int log = 1;
    int n = 1;
    int dim = 3;
    int order = 1;
@@ -104,9 +103,9 @@ TEST_CASE("3D ProjectBdrCoefficientTangent",
 
                   nd_err += nd_dist;
 
-                  if (log > 0 && nd_dist > tol)
+                  if (verbose_tests && nd_dist > tol)
                   {
-                     std::cout << be << ":" << j << " nd ("
+                     mfem::out << be << ":" << j << " nd ("
                                << f_val[0] << "," << f_val[1] << "," << f_val[2]
                                << ") vs. ("
                                << nd_val[0] << "," << nd_val[1] << ","

tests/unit/fem/test_quadf_coef.cpp

@@ -77,9 +77,7 @@ TEST_CASE("Quadrature Function Coefficients",
 
    SECTION("Operators on VecQuadFuncCoeff")
    {
-      std::cout << "Testing VecQuadFuncCoeff: " << std::endl;
 #ifdef MFEM_USE_EXCEPTIONS
-      std::cout << " Setting Component" << std::endl;
       REQUIRE_THROWS(qfvc.SetComponent(3, 1));
       REQUIRE_THROWS(qfvc.SetComponent(-1, 1));
       REQUIRE_NOTHROW(qfvc.SetComponent(1, 2));
@@ -93,7 +91,6 @@ TEST_CASE("Quadrature Function Coefficients",
 
    SECTION("Operators on VectorQuadratureLFIntegrator")
    {
-      std::cout << "Testing VectorQuadratureLFIntegrator: " << std::endl;
       H1_FECollection    fec_h1(order_h1, dim);
       FiniteElementSpace fespace_h1(&mesh, &fec_h1, dim);
 
@@ -124,7 +121,6 @@ TEST_CASE("Quadrature Function Coefficients",
 
    SECTION("Operators on QuadratureLFIntegrator")
    {
-      std::cout << "Testing QuadratureLFIntegrator: " << std::endl;
       H1_FECollection    fec_h1(order_h1, dim);
       FiniteElementSpace fespace_h1(&mesh, &fec_h1, 1);
       FiniteElementSpace fespace_h3(&mesh, &fec_h1, 3);

tests/unit/fem/test_quadinterpolator.cpp

@@ -21,8 +21,7 @@ static bool testQuadratureInterpolator(const int dim,
                                        const int nx, const int ny, const int nz)
 {
    // Keep for debugging purposes:
-   const bool verbose = false;
-   if (verbose)
+   if (verbose_tests)
    {
       std::cout << "testQuadratureInterpolator(dim=" << dim
                 << ",p=" << p
@@ -145,21 +144,21 @@ static bool testQuadratureInterpolator(const int dim,
       norm = sq_val_f.Normlinf();
       sq_val_f -= sq_val_t;
       rel_error = sq_val_f.Normlinf()/norm;
-      if (verbose)
+      if (verbose_tests)
       { std::cout << "sq_val rel. error = " << rel_error << std::endl; }
       REQUIRE(rel_error <= rel_tol);
 
       norm = sq_der_f.Normlinf();
       sq_der_f -= sq_der_t;
       rel_error = sq_der_f.Normlinf()/norm;
-      if (verbose)
+      if (verbose_tests)
       { std::cout << "sq_der rel. error = " << rel_error << std::endl; }
       REQUIRE(rel_error <= rel_tol);
 
       norm = sq_pdr_f.Normlinf();
       sq_pdr_f -= sq_pdr_t;
       rel_error = sq_pdr_f.Normlinf()/norm;
-      if (verbose)
+      if (verbose_tests)
       { std::cout << "sq_pdr rel. error = " << rel_error << std::endl; }
       REQUIRE(rel_error <= rel_tol);
    }
@@ -209,28 +208,28 @@ static bool testQuadratureInterpolator(const int dim,
       norm = vq_val_f.Normlinf();
       vq_val_f -= vq_val_t;
       rel_error = vq_val_f.Normlinf()/norm;
-      if (verbose)
+      if (verbose_tests)
       { std::cout << "vq_val rel. error = " << rel_error << std::endl; }
       REQUIRE(rel_error <= rel_tol);
 
       norm = vq_der_f.Normlinf();
       vq_der_f -= vq_der_t;
       rel_error = vq_der_f.Normlinf()/norm;
-      if (verbose)
+      if (verbose_tests)
       { std::cout << "vq_der rel. error = " << rel_error << std::endl; }
       REQUIRE(rel_error <= rel_tol);
 
       norm = vq_det_f.Normlinf();
       vq_det_f -= vq_det_t;
       rel_error = vq_det_f.Normlinf()/norm;
-      if (verbose)
+      if (verbose_tests)
       { std::cout << "vq_det rel. error = " << rel_error << std::endl; }
       REQUIRE(rel_error <= rel_tol);
 
       norm = vq_pdr_f.Normlinf();
       vq_pdr_f -= vq_pdr_t;
       rel_error = vq_pdr_f.Normlinf()/norm;
-      if (verbose)
+      if (verbose_tests)
       { std::cout << "vq_pdr rel. error = " << rel_error << std::endl; }
       REQUIRE(rel_error <= rel_tol);
    }

tests/unit/fem/test_quadraturefunc.cpp

@@ -1315,8 +1315,6 @@ TEST_CASE("1D Quadrature Functions")
    mfem::QuadratureFunctions1D quad_func;
 
    // The tests will be reported in these sections.
-   // Each REQUIRE counts as an assertion.
-   // true = pass, false = fail
    SECTION("Gauss-Legendre")
    {
       const int np = 21;
@@ -1337,8 +1335,7 @@ TEST_CASE("1D Quadrature Functions")
             double err_w = std::fabs( ir.IntPoint(i).weight - w_tbl[i] );
             if ( (err_x > tol) || (err_w > tol) )
             {
-               std::cout << "Gauss Legendre with " << n << " points wrong\n";
-               REQUIRE(false);
+               FAIL("Gauss Legendre with " << n << " points wrong");
                break;
             }
          }
@@ -1365,8 +1362,7 @@ TEST_CASE("1D Quadrature Functions")
             double err_w = std::fabs( ir.IntPoint(i).weight - w_tbl[i] );
             if ( (err_x > tol) || (err_w > tol) )
             {
-               std::cout << "Gauss Lobatto with " << n << " points wrong\n";
-               REQUIRE(false);
+               FAIL("Gauss Lobatto with " << n << " points wrong");
                break;
             }
          }
@@ -1395,8 +1391,7 @@ TEST_CASE("1D Quadrature Functions")
             double err_w = std::fabs( ir.IntPoint(i).weight - w_tbl[i] );
             if ( (err_x > tol) || (err_w > tol) )
             {
-               std::cout << "Closed Newton-Cotes with " << n << " points wrong\n";
-               REQUIRE(false);
+               FAIL("Closed Newton-Cotes with " << n << " points wrong");
                break;
             }
          }
@@ -1425,8 +1420,7 @@ TEST_CASE("1D Quadrature Functions")
             double err_w = std::fabs( ir.IntPoint(i).weight - w_tbl[i] );
             if ( (err_x > tol) || (err_w > tol) )
             {
-               std::cout << "Open Newton-Cotes with " << n << " points wrong\n";
-               REQUIRE(false);
+               FAIL("Open Newton-Cotes with " << n << " points wrong");
                break;
             }
          }

tests/unit/linalg/test_chebyshev.cpp

@@ -57,7 +57,7 @@ TEST_CASE("OperatorChebyshevSmoother", "[Chebyshev symmetry]")
       double transpose_val = right * smooth;
 
       double error = fabs(forward_val - transpose_val) / fabs(forward_val);
-      std::cout << "Order " << order << " symmetry error: " << error << std::endl;
+      CAPTURE(order, error);
       REQUIRE(error < 1.e-13);
 
       delete smoother;

tests/unit/linalg/test_constrainedsolver.cpp

@@ -177,7 +177,7 @@ TEST_CASE("SerialConstrainedSolver", "[ConstrainedSolver]")
 TEST_CASE("ConstrainedSolver", "[Parallel], [ConstrainedSolver]")
 {
 #ifdef HYPRE_USING_GPU
-   std::cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
+   mfem::out << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
              << "is NOT supported with the GPU version of hypre.\n\n";
    return;
 #endif
@@ -382,7 +382,7 @@ void ParallelTestProblem::Penalty(double pen, Vector& serr, Vector& lerr)
 TEST_CASE("ParallelConstrainedSolver", "[Parallel], [ConstrainedSolver]")
 {
 #ifdef HYPRE_USING_GPU
-   std::cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
+   mfem::out << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
              << "is NOT supported with the GPU version of hypre.\n\n";
    return;
 #endif
@@ -672,7 +672,7 @@ void ParallelTestProblemTwo::Penalty(double pen, Vector& serr, Vector& lerr)
 TEST_CASE("ParallelConstrainedSolverTwo", "[Parallel], [ConstrainedSolver]")
 {
 #ifdef HYPRE_USING_GPU
-   std::cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
+   mfem::out << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
              << "is NOT supported with the GPU version of hypre.\n\n";
    return;
 #endif
@@ -812,7 +812,7 @@ void ZerosTestProblem::Elimination(Vector& serr, Vector& lerr, bool twoblocks)
 TEST_CASE("ZerosTestCase", "[Parallel], [ConstrainedSolver]")
 {
 #ifdef HYPRE_USING_GPU
-   std::cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
+   mfem::out << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
              << "is NOT supported with the GPU version of hypre.\n\n";
    return;
 #endif

tests/unit/linalg/test_hypre_ilu.cpp

@@ -23,7 +23,7 @@ TEST_CASE("HypreILU and HypreFGMRES wrappers",
           "[Parallel], [HypreILU], [HypreFGMRES]")
 {
 #ifdef HYPRE_USING_GPU
-   std::cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
+   mfem::out << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this unit test\n"
              << "is NOT supported with the GPU version of hypre.\n\n";
    return;
 #endif

tests/unit/linalg/test_matrix_block.cpp

@@ -76,8 +76,9 @@ TEST_CASE("BlockMatrix", "[BlockMatrix]")
       {
          if ( A->RowSize(i) != Amono->RowSize(i) )
          {
-            std::cout<< "BlockMatrix::RowSize failure: " << i <<"\t"<< A->RowSize(
-                        i) <<"\t" <<  Amono->RowSize(i) << "\n";
+            UNSCOPED_INFO("BlockMatrix::RowSize failure: " << i << '\t'
+                          << A->RowSize( i) << '\t'
+                          << Amono->RowSize(i));
             ++nfails;
          }
       }

tests/unit/linalg/test_matrix_hypre.cpp

@@ -88,7 +88,7 @@ TEST_CASE("HypreParMatrixAbsMult",  "[Parallel], [HypreParMatrixAbsMult]")
       Y1 -=Y0;
       double error = Y1.Norml2();
 
-      std::cout << "Testing AbsMult:   order: " << order
+      mfem::out << "Testing AbsMult:   order: " << order
                 << ", error norm on rank "
                 << rank << ": " << error << std::endl;
 
@@ -104,7 +104,7 @@ TEST_CASE("HypreParMatrixAbsMult",  "[Parallel], [HypreParMatrixAbsMult]")
       X1 -=X0;
 
       error = X1.Norml1();
-      std::cout << "Testing AbsMultT:  order: " << order
+      mfem::out << "Testing AbsMultT:  order: " << order
                 << ", error norm on rank "
                 << rank << ": " << error << std::endl;
 

tests/unit/linalg/test_matrix_rectangular.cpp

@@ -257,7 +257,7 @@ TEST_CASE("HypreParMatrixBlocksRectangular",
 
       yH -= yB;
       double error = yH.Norml2();
-      std::cout << "  order: " << order
+      mfem::out << "  order: " << order
                 << ", block matrix error norm on rank " << rank << ": " << error << std::endl;
       REQUIRE(error < 1.e-12);
 

tests/unit/linalg/test_matrix_sparse.cpp

@@ -21,6 +21,8 @@ TEST_CASE("SparseMatrixAbsMult", "[SparseMatrixAbsMult]")
    int ne = 4;
    for (int order = 1; order <= 3; ++order)
    {
+      CAPTURE(order);
+
       Mesh mesh = Mesh::MakeCartesian2D(
                      ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
       FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
@@ -56,10 +58,6 @@ TEST_CASE("SparseMatrixAbsMult", "[SparseMatrixAbsMult]")
       Y1 -=Y0;
       double error = Y1.Norml2();
 
-      std::cout << "Testing AbsMult:   order: " << order
-                << ", error norm: "
-                << error << std::endl;
-
       REQUIRE(error == MFEM_Approx(0.0));
 
       Y0.Randomize();
@@ -71,10 +69,6 @@ TEST_CASE("SparseMatrixAbsMult", "[SparseMatrixAbsMult]")
 
       error = X1.Norml2();
 
-      std::cout << "Testing AbsMultT:  order: " << order
-                << ", error norm: "
-                << error << std::endl;
-
       REQUIRE(error == MFEM_Approx(0.0));
 
       delete Aabs;

tests/unit/linalg/test_matrix_square.cpp

@@ -23,8 +23,8 @@ TEST_CASE("FormLinearSystem", "[FormLinearSystem]")
    {
       for (int ne = 1; ne <= 4; ++ne)
       {
-         std::cout << "Testing " << dim << "D partial assembly: "
-                   << std::pow(ne, dim) << " elements." << std::endl;
+         const int n_elements = std::pow(ne, dim);
+         CAPTURE(dim, n_elements);
          for (int order = 1; order <= 3; ++order)
          {
             Mesh mesh;
@@ -74,7 +74,7 @@ TEST_CASE("FormLinearSystem", "[FormLinearSystem]")
 
             x0 -= x1;
             double error = x0.Norml2();
-            std::cout << "    order: " << order << ", error norm: " << error << std::endl;
+            CAPTURE(error, order);
             REQUIRE(x0.Norml2() == MFEM_Approx(0.0, 1e2*EPS));
 
             delete fec;
@@ -91,8 +91,8 @@ TEST_CASE("ParallelFormLinearSystem", "[Parallel], [ParallelFormLinearSystem]")
    {
       for (int ne = 4; ne <= 5; ++ne)
       {
-         std::cout << "Testing " << dim << "D partial assembly: "
-                   << std::pow(ne, dim) << " elements." << std::endl;
+         const int n_elements = std::pow(ne, dim);
+         CAPTURE(dim, n_elements);
          for (int order = 1; order <= 3; ++order)
          {
             Mesh mesh;
@@ -159,7 +159,7 @@ TEST_CASE("ParallelFormLinearSystem", "[Parallel], [ParallelFormLinearSystem]")
 
             x0 -= x1;
             double error = x0.Norml2();
-            std::cout << "    order: " << order << ", error norm: " << error << std::endl;
+            CAPTURE(order, error);
             REQUIRE(x0.Norml2() == MFEM_Approx(0.0, 2e2*EPS));
 
             delete pmesh;
@@ -260,7 +260,7 @@ TEST_CASE("HypreParMatrixBlocksSquare",
 
       yH -= yB;
       double error = yH.Norml2();
-      std::cout << "  order: " << order
+      mfem::out << "  order: " << order
                 << ", block matrix error norm on rank " << rank << ": " << error << std::endl;
       REQUIRE(error < EPS);
 

tests/unit/linalg/test_ode.cpp

@@ -15,8 +15,7 @@
 
 using namespace mfem;
 
-TEST_CASE("First order ODE methods",
-          "[ODE1]")
+TEST_CASE("First order ODE methods", "[ODE]")
 {
    double tol = 0.1;
 
@@ -120,10 +119,10 @@ TEST_CASE("First order ODE methods",
          u +=u0;
          error[0] = u.Norml2();
 
-         std::cout<<std::setw(12)<<"Error"
+         mfem::out<<std::setw(12)<<"Error"
                   <<std::setw(12)<<"Ratio"
                   <<std::setw(12)<<"Order"<<std::endl;
-         std::cout<<std::setw(12)<<error[0]<<std::endl;
+         mfem::out<<std::setw(12)<<error[0]<<std::endl;
 
          std::vector<Vector> uh(ode_solver->GetMaxStateSize());
          for (int l = 1; l < levels; l++)
@@ -168,7 +167,7 @@ TEST_CASE("First order ODE methods",
 
             u += u0;
             error[l] = u.Norml2();
-            std::cout<<std::setw(12)<<error[l]
+            mfem::out<<std::setw(12)<<error[l]
                      <<std::setw(12)<<error[l-1]/error[l]
                      <<std::setw(12)<<log(error[l-1]/error[l])/log(2)
                      <<std::endl;
@@ -184,117 +183,100 @@ TEST_CASE("First order ODE methods",
    // Implicit L-stable methods
    SECTION("BackwardEuler")
    {
-      std::cout <<"\nTesting BackwardEuler" << std::endl;
       double conv_rate = check.order(new BackwardEulerSolver);
       REQUIRE(conv_rate + tol > 1.0);
    }
 
    SECTION("SDIRK23Solver(2)")
    {
-      std::cout <<"\nTesting SDIRK23Solver(2)" << std::endl;
       double conv_rate = check.order(new SDIRK23Solver(2));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("SDIRK33Solver")
    {
-      std::cout <<"\nTesting SDIRK33Solver" << std::endl;
       double conv_rate = check.order(new SDIRK33Solver);
       REQUIRE(conv_rate + tol > 3.0);
    }
 
    SECTION("ForwardEulerSolver")
    {
-      std::cout <<"\nTesting ForwardEulerSolver" << std::endl;
       double conv_rate = check.order(new ForwardEulerSolver);
       REQUIRE(conv_rate + tol > 1.0);
    }
 
    SECTION("RK2Solver(0.5)")
    {
-      std::cout <<"\nTesting RK2Solver(0.5)" << std::endl;
       double conv_rate = check.order(new RK2Solver(0.5));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("RK3SSPSolver")
    {
-      std::cout <<"\nTesting RK3SSPSolver" << std::endl;
       double conv_rate = check.order(new RK3SSPSolver);
       REQUIRE(conv_rate + tol > 3.0);
    }
 
    SECTION("RK4Solver")
    {
-      std::cout <<"\nTesting RK4Solver" << std::endl;
       double conv_rate = check.order(new RK4Solver);
       REQUIRE(conv_rate + tol > 4.0);
    }
 
    SECTION("ImplicitMidpointSolver")
    {
-      std::cout <<"\nTesting ImplicitMidpointSolver" << std::endl;
       double conv_rate = check.order(new ImplicitMidpointSolver);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("SDIRK23Solver")
    {
-      std::cout <<"\nTesting SDIRK23Solver" << std::endl;
       double conv_rate = check.order(new SDIRK23Solver);
       REQUIRE(conv_rate + tol > 3.0);
    }
 
    SECTION("SDIRK34Solver")
    {
-      std::cout <<"\nTesting SDIRK34Solver" << std::endl;
       double conv_rate = check.order(new SDIRK34Solver);
       REQUIRE(conv_rate + tol > 4.0);
    }
 
    SECTION("TrapezoidalRuleSolver")
    {
-      std::cout <<"\nTesting TrapezoidalRuleSolver" << std::endl;
       REQUIRE(check.order(new TrapezoidalRuleSolver) + tol > 2.0 );
    }
 
    SECTION("ESDIRK32Solver")
    {
-      std::cout <<"\nTesting ESDIRK32Solver" << std::endl;
       REQUIRE(check.order(new ESDIRK32Solver) + tol > 2.0 );
    }
 
    SECTION("ESDIRK33Solver")
    {
-      std::cout <<"\nTesting ESDIRK33Solver" << std::endl;
       REQUIRE(check.order(new ESDIRK33Solver) + tol > 3.0 );
    }
 
    // Generalized-alpha
    SECTION("GeneralizedAlphaSolver(1.0)")
    {
-      std::cout <<"\nTesting GeneralizedAlphaSolver(1.0)" << std::endl;
       double conv_rate = check.order(new GeneralizedAlphaSolver(1.0));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("GeneralizedAlphaSolver(0.5)")
    {
-      std::cout <<"\nTesting GeneralizedAlphaSolver(0.5)" << std::endl;
       double conv_rate = check.order(new GeneralizedAlphaSolver(0.5));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("GeneralizedAlphaSolver(0.5) - restart")
    {
-      std::cout <<"\nTesting GeneralizedAlphaSolver(0.5) - restart" << std::endl;
       double conv_rate = check.order(new GeneralizedAlphaSolver(0.5), true);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("GeneralizedAlphaSolver(0.0)")
    {
-      std::cout <<"\nTesting GeneralizedAlphaSolver(0.0)" << std::endl;
       double conv_rate = check.order(new GeneralizedAlphaSolver(0.0));
       REQUIRE(conv_rate + tol > 2.0);
    }
@@ -302,56 +284,48 @@ TEST_CASE("First order ODE methods",
    // Adams-Bashforth
    SECTION("AB1Solver()")
    {
-      std::cout <<"\nTesting AB1Solver()" << std::endl;
       double conv_rate = check.order(new AB1Solver());
       REQUIRE(conv_rate + tol > 1.0);
    }
 
    SECTION("AB1Solver() - restart")
    {
-      std::cout <<"\nTesting AB1Solver() - restart" << std::endl;
       double conv_rate = check.order(new AB1Solver(), true);
       REQUIRE(conv_rate + tol > 1.0);
    }
 
    SECTION("AB2Solver()")
    {
-      std::cout <<"\nTesting AB2Solver()" << std::endl;
       double conv_rate = check.order(new AB2Solver());
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("AB2Solver() - restart")
    {
-      std::cout <<"\nTesting AB2Solver() - restart" << std::endl;
       double conv_rate = check.order(new AB2Solver(), true);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("AB3Solver()")
    {
-      std::cout <<"\nTesting AB3Solver()" << std::endl;
       double conv_rate = check.order(new AB3Solver());
       REQUIRE(conv_rate + tol > 3.0);
    }
 
    SECTION("AB4Solver()")
    {
-      std::cout <<"\nTesting AB4Solver()" << std::endl;
       double conv_rate = check.order(new AB4Solver());
       REQUIRE(conv_rate + tol > 4.0);
    }
 
    SECTION("AB5Solver()")
    {
-      std::cout <<"\nTesting AB5Solver()" << std::endl;
       double conv_rate = check.order(new AB5Solver());
       REQUIRE(conv_rate + tol > 5.0);
    }
 
    SECTION("AB5Solver() - restart")
    {
-      std::cout <<"\nTesting AB5Solver() - restart" << std::endl;
       double conv_rate = check.order(new AB5Solver(), true);
       REQUIRE(conv_rate + tol > 5.0);
    }
@@ -359,56 +333,48 @@ TEST_CASE("First order ODE methods",
    // Adams-Moulton
    SECTION("AM0Solver()")
    {
-      std::cout <<"\nTesting AM0Solver()" << std::endl;
       double conv_rate = check.order(new AM0Solver());
       REQUIRE(conv_rate + tol > 1.0);
    }
 
    SECTION("AM1Solver()")
    {
-      std::cout <<"\nTesting AM1Solver()" << std::endl;
       double conv_rate = check.order(new AM1Solver());
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("AM1Solver() - restart")
    {
-      std::cout <<"\nTesting AM1Solver() - restart" << std::endl;
       double conv_rate = check.order(new AM1Solver(), true);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("AM2Solver()")
    {
-      std::cout <<"\nTesting AM2Solver()" << std::endl;
       double conv_rate = check.order(new AM2Solver());
       REQUIRE(conv_rate + tol > 3.0);
    }
 
    SECTION("AM2Solver() - restart")
    {
-      std::cout <<"\nTesting AM2Solver() - restart" << std::endl;
       double conv_rate = check.order(new AM2Solver(), true);
       REQUIRE(conv_rate + tol > 1.0);
    }
 
    SECTION("AM3Solver()")
    {
-      std::cout <<"\nTesting AM3Solver()" << std::endl;
       double conv_rate = check.order(new AM3Solver());
       REQUIRE(conv_rate + tol > 4.0);
    }
 
    SECTION("AM4Solver()")
    {
-      std::cout <<"\nTesting AM4Solver()" << std::endl;
       double conv_rate = check.order(new AM4Solver());
       REQUIRE(conv_rate + tol > 5.0);
    }
 
    SECTION("AM4Solver() - restart")
    {
-      std::cout <<"\nTesting AM4Solver() - restart" << std::endl;
       double conv_rate = check.order(new AM4Solver(),true);
       REQUIRE(conv_rate + tol > 5.0);
    }

tests/unit/linalg/test_ode2.cpp

@@ -15,8 +15,7 @@
 
 using namespace mfem;
 
-TEST_CASE("Second order ODE methods",
-          "[ODE2]")
+TEST_CASE("Second order ODE methods", "[ODE]")
 {
    double tol = 0.1;
 
@@ -116,13 +115,13 @@ TEST_CASE("Second order ODE methods",
          err_u[0] = u.Norml2();
          err_du[0] = du.Norml2();
 
-         std::cout<<std::setw(12)<<"Error u"
+         mfem::out<<std::setw(12)<<"Error u"
                   <<std::setw(12)<<"Error du"
                   <<std::setw(12)<<"Ratio u"
                   <<std::setw(12)<<"Ratio du"
                   <<std::setw(12)<<"Order u"
                   <<std::setw(12)<<"Order du"<<std::endl;
-         std::cout<<std::setw(12)<<err_u[0]
+         mfem::out<<std::setw(12)<<err_u[0]
                   <<std::setw(12)<<err_du[0]<<std::endl;
 
          std::vector<Vector> uh(ode_solver->GetMaxStateSize());
@@ -172,7 +171,7 @@ TEST_CASE("Second order ODE methods",
             du -= dudt0;
             err_u[l] = u.Norml2();
             err_du[l] = du.Norml2();
-            std::cout<<std::setw(12)<<err_u[l]
+            mfem::out<<std::setw(12)<<err_u[l]
                      <<std::setw(12)<<err_du[l]
                      <<std::setw(12)<<err_u[l-1]/err_u[l]
                      <<std::setw(12)<<err_du[l-1]/err_du[l]
@@ -190,28 +189,24 @@ TEST_CASE("Second order ODE methods",
    // Newmark-based solvers
    SECTION("Newmark")
    {
-      std::cout <<"\nTesting NewmarkSolver" << std::endl;
       double conv_rate = check.order(new NewmarkSolver);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("LinearAcceleration")
    {
-      std::cout <<"\nLinearAccelerationSolver" << std::endl;
       double conv_rate = check.order(new LinearAccelerationSolver);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("CentralDifference")
    {
-      std::cout <<"\nTesting CentralDifference" << std::endl;
       double conv_rate = check.order(new CentralDifferenceSolver);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("FoxGoodwin")
    {
-      std::cout <<"\nTesting FoxGoodwin" << std::endl;
       double conv_rate = check.order(new FoxGoodwinSolver);
       REQUIRE(conv_rate + tol > 4.0);
    }
@@ -219,28 +214,24 @@ TEST_CASE("Second order ODE methods",
    // Generalized-alpha based solvers
    SECTION("GeneralizedAlpha(0.0)")
    {
-      std::cout <<"\nTesting GeneralizedAlpha(0.0)" << std::endl;
       double conv_rate = check.order(new GeneralizedAlpha2Solver(0.0));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("GeneralizedAlpha(0.5)")
    {
-      std::cout <<"\nTesting GeneralizedAlpha(0.5)" << std::endl;
       double conv_rate = check.order(new GeneralizedAlpha2Solver(0.5));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("GeneralizedAlpha(0.5) - restart")
    {
-      std::cout <<"\nTesting GeneralizedAlpha(0.5) - restart" << std::endl;
       double conv_rate = check.order(new GeneralizedAlpha2Solver(0.5),true);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("GeneralizedAlpha(1.0)")
    {
-      std::cout <<"\nTesting GeneralizedAlpha(1.0)" << std::endl;
       double conv_rate = check.order(new GeneralizedAlpha2Solver(1.0));
       REQUIRE(conv_rate + tol > 2.0);
    }
@@ -248,49 +239,42 @@ TEST_CASE("Second order ODE methods",
 
    SECTION("AverageAcceleration")
    {
-      std::cout <<"\nTesting AverageAcceleration" << std::endl;
       double conv_rate = check.order(new AverageAccelerationSolver);
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("HHTAlpha(2/3)")
    {
-      std::cout <<"\nTesting HHTAlpha(2/3)" << std::endl;
       double conv_rate = check.order(new HHTAlphaSolver(2.0/3.0));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("HHTAlpha(0.75)")
    {
-      std::cout <<"\nTesting HHTAlpha(0.75)" << std::endl;
       double conv_rate = check.order(new HHTAlphaSolver(0.75));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("HHTAlpha(1.0)")
    {
-      std::cout <<"\nTesting HHTAlpha(1.0)" << std::endl;
       double conv_rate = check.order(new HHTAlphaSolver(1.0));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("WBZAlpha(0.0)")
    {
-      std::cout <<"\nTesting WBZAlpha(0.0)" << std::endl;
       double conv_rate = check.order(new WBZAlphaSolver(0.0));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("WBZAlpha(0.5)")
    {
-      std::cout <<"\nTesting WBZAlpha(0.5)" << std::endl;
       double conv_rate = check.order(new WBZAlphaSolver(0.5));
       REQUIRE(conv_rate + tol > 2.0);
    }
 
    SECTION("WBZAlpha(1.0)")
    {
-      std::cout <<"\nTesting WBZAlpha(1.0)" << std::endl;
       double conv_rate = check.order(new WBZAlphaSolver(1.0));
       REQUIRE(conv_rate + tol > 2.0);
    }

tests/unit/mesh/test_ncmesh.cpp

@@ -62,8 +62,7 @@ TEST_CASE("NCMesh PA diagonal", "[NCMesh]")
          nc_a.AssembleDiagonal(nc_diag);
 
          double error = fabs(diag.Norml2() - nc_diag.Norml2());
-         std::cout << "Testing quad NCMesh PA diag:    "
-                   "order: " << order << ", error: " << error << std::endl;
+         CAPTURE(order, error);
          REQUIRE(error == MFEM_Approx(0.0, EPS));
       }
    }
@@ -107,8 +106,7 @@ TEST_CASE("NCMesh PA diagonal", "[NCMesh]")
          nc_a.AssembleDiagonal(nc_diag);
 
          double error = fabs(diag.Sum() - nc_diag.Sum());
-         std::cout << "Testing hexa NCMesh PA diag:    "
-                   "order: " << order << ", error: " << error << std::endl;
+         CAPTURE(order, error);
          REQUIRE(error == MFEM_Approx(0.0, EPS));
       }
    }
@@ -174,11 +172,7 @@ TEST_CASE("pNCMesh PA diagonal",  "[Parallel], [NCMesh]")
          MPI_Allreduce(&nc_diag_lsum, &nc_diag_gsum, 1, MPI_DOUBLE, MPI_SUM,
                        MPI_COMM_WORLD);
          double error = fabs(diag_gsum - nc_diag_gsum);
-         if (rank==0)
-         {
-            std::cout << "Testing quad pNCMesh PA diag:    "
-                      "order: " << order << ", error: " << error << std::endl;
-         }
+         CAPTURE(order, error);
          REQUIRE(error == MFEM_Approx(0.0, EPS));
          MPI_Barrier(MPI_COMM_WORLD);
       }
@@ -232,11 +226,7 @@ TEST_CASE("pNCMesh PA diagonal",  "[Parallel], [NCMesh]")
          MPI_Allreduce(&nc_diag_lsum, &nc_diag_gsum, 1, MPI_DOUBLE, MPI_SUM,
                        MPI_COMM_WORLD);
          double error = fabs(diag_gsum - nc_diag_gsum);
-         if (rank==0)
-         {
-            std::cout << "Testing hexa pNCMesh PA diag:    "
-                      "order: " << order << ", error: " << error << std::endl;
-         }
+         CAPTURE(order, error);
          REQUIRE(error == MFEM_Approx(0.0, EPS));
          MPI_Barrier(MPI_COMM_WORLD);
       }

tests/unit/pcunit_test_main.cpp

@@ -14,9 +14,6 @@
 #include "mfem.hpp"
 #include "run_unit_tests.hpp"
 
-bool launch_all_non_regression_tests = false;
-std::string mfem_data_dir;
-
 #ifndef MFEM_USE_MPI
 #error "This test should be disabled without MFEM_USE_MPI!"
 #endif

tests/unit/punit_test_main.cpp

@@ -14,9 +14,6 @@
 #include "mfem.hpp"
 #include "run_unit_tests.hpp"
 
-bool launch_all_non_regression_tests = false;
-std::string mfem_data_dir;
-
 #ifndef MFEM_USE_MPI
 #error "This test should be disabled without MFEM_USE_MPI!"
 #endif

tests/unit/run_unit_tests.hpp

@@ -14,6 +14,10 @@
 
 #include "unit_tests.hpp"
 
+bool launch_all_non_regression_tests = false;
+bool verbose_tests = false;
+std::string mfem_data_dir;
+
 static int RunCatchSession(int argc, char *argv[],
                            const std::vector<std::string> &testsOrTags,
                            bool root=true)
@@ -23,9 +27,12 @@ static int RunCatchSession(int argc, char *argv[],
 
    // Build a new command line parser on top of Catch's
    using namespace Catch::clara;
+   bool enable_mfem_out = false;
    auto cli = session.cli()
               | Opt(launch_all_non_regression_tests) ["--all"] ("all tests")
-              | Opt(mfem_data_dir, "") ["--data"] ("mfem/data repository");
+              | Opt(mfem_data_dir, "") ["--data"] ("mfem/data repository")
+              | Opt(enable_mfem_out) ["--enable-output"] ("enable mfem::out")
+              | Opt(verbose_tests) ["--verbose"] ("verbose output");
    session.cli(cli);
 
    // For floating point comparisons, print 8 digits for single precision
@@ -43,6 +50,11 @@ static int RunCatchSession(int argc, char *argv[],
    if (mfem_data_dir == "") { cfg.testsOrTags.push_back("~[MFEMData]"); }
    session.useConfigData(cfg);
 
+   if (!verbose_tests && !enable_mfem_out)
+   {
+      mfem::out.Disable();
+   }
+
    if (root)
    {
       std::cout << "INFO: Test filter: ";

tests/unit/unit_test_main.cpp

@@ -13,9 +13,6 @@
 #include "mfem.hpp"
 #include "run_unit_tests.hpp"
 
-bool launch_all_non_regression_tests = false;
-std::string mfem_data_dir;
-
 int main(int argc, char *argv[])
 {
    // Exclude parallel tests.

tests/unit/unit_tests.hpp

@@ -22,6 +22,10 @@ extern bool launch_all_non_regression_tests;
     and tests tagged with [MFEMData] will be skipped. */
 extern std::string mfem_data_dir;
 
+/// Command line '--verbose' argument, false by default.
+/** Enables extra log messages and the mfem::out stream. */
+extern bool verbose_tests;
+
 /** @brief MFEM_Approx can be used to compare floating point values within an
     absolute tolerance of @a abs_tol (default value 1e-12) and relative
     tolerance of @a rel_tol (default value 1e-12). */