formatting and include desired qfunction

Added CMake dependencies on individual data files

CMakeLists.txt

@@ -723,6 +723,7 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})
 
 # Declaring the library
 mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
+target_compile_features(mfem PUBLIC cxx_std_${CMAKE_CXX_STANDARD})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
 target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
 if (TPL_TARGETS)
@@ -869,11 +870,12 @@ add_dependencies(exec
 #  - https://cmake.org/Bug/view.php?id=8438
 
 # Add a target to copy the mfem data directory to the build directory
-add_custom_command(OUTPUT data_is_copied
-  COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/data data
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying the data directory ...")
-add_custom_target(copy_data DEPENDS data_is_copied)
+# Implementable as a single copy_directory_if_different command w/ CMake >= 3.26
+file(GLOB DATA_FILES CONFIGURE_DEPENDS ${PROJECT_SOURCE_DIR}/data/*)
+add_custom_target(copy_data
+  COMMAND ${CMAKE_COMMAND} -E make_directory data
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${DATA_FILES} data
+  COMMENT "Syncing the data directory ...")
 # Add 'copy_data' as a prerequisite for all executables, if the source and the
 # build directories are not the same.
 if (NOT ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}"))
@@ -1005,6 +1007,10 @@ install(FILES
 install(EXPORT ${PROJECT_NAME_UC}Targets
     DESTINATION ${INSTALL_CMAKE_DIR})
 
+# Install the data directory if present, i.e. if the copy_data target is built
+install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
+  DESTINATION ${MFEM_INSTALL_DIR} OPTIONAL)
+
 #-------------------------------------------------------------------------------
 # Create 'config.mk' from 'config.mk.in' for the build and install locations and
 # define install rules for 'config.mk' and 'test.mk'

INSTALL

@@ -725,7 +725,9 @@ The specific libraries and their options are:
   URL: https://ginkgo-project.github.io
   Options: GINKGO_OPT, GINKGO_LIB, GINKGO_DIR, GINKGO_BUILD_TYPE (Release or
            Debug).
-  Versions: Ginkgo >= 1.9.0.
+  Versions: Ginkgo >= 1.9.0. When building Ginkgo with distributed support, a
+            recent version of the "develop" branch is required (1.11 as defined
+            in include/ginkgo/config.hpp).
 
 - AmgX (optional), used when MFEM_USE_AMGX = YES.
   URL: https://github.com/NVIDIA/AMGX

config/defaults.mk

@@ -407,7 +407,7 @@ AMGX_LIB = -L$(AMGX_DIR)/lib -lamgx -lcusparse -lcusolver -lcublas -lnvToolsExt
 # MAGMA library configuration
 MAGMA_DIR = @MFEM_DIR@/../magma
 MAGMA_OPT = -I$(MAGMA_DIR)/include
-MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a -lcublas -lcusparse $(LAPACK_LIB)
+MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a $(LAPACK_LIB)
 
 # GnuTLS library configuration
 GNUTLS_OPT =

doc/CodeDocumentation.dox

@@ -117,6 +117,8 @@ namespace mfem {
  * - <a class="el" href="ex39p_8cpp_source.html">Example 39p</a>: parallel named mesh attributes
  * - <a class="el" href="ex40_8cpp_source.html">Example 40</a>: eikonal equation
  * - <a class="el" href="ex40p_8cpp_source.html">Example 40p</a>: parallel eikonal equation
+ * - <a class="el" href="ex41_8cpp_source.html">Example 41</a>: DG/CG IMEX time-dependent advection-diffusion
+ * - <a class="el" href="ex41p_8cpp_source.html">Example 41p</a>: parallel DG/CG IMEX time-dependent advection-diffusion
  *
  * <H4>AmgX Examples</H4>
  * - Variants of Examples
@@ -234,7 +236,8 @@ namespace mfem {
  * - <a class="el" href="miniapps_2performance_2ex1_8cpp_source.html">HPC Example 1</a>: high-performance nodal H1 FEM for the Poisson problem
  * - <a class="el" href="miniapps_2performance_2ex1p_8cpp_source.html">HPC Example 1p</a>: high-performance parallel nodal H1 FEM for the Poisson problem
  * - <a class="el" href="generate__random__field_8cpp_source.html">SPDE Solvers</a>: SPDE solver random field generation
- * - <a class="el" href="contact-patch-test_8cpp_source.html">Contact</a>: mortar contact patch test for elasticity
+ * - <a class="el" href="contact-patch-test_8cpp_source.html">Tribol</a>: mortar contact patch test for elasticity
+ * - <a class="el" href="contact_8cpp_source.html">Contact</a>: Frictionless contact examples using <a class="el" href="classmfem_1_1IPSolver.html#details">IP optimization</a> and the <a class="el" href="classmfem_1_1AMGFSolver.html#details">AMGF solver</a>
  * - <a class="el" href="multidomain_8cpp_source.html">Multidomain miniapp</a>: Multidomain and Submesh demonstration miniapp
  * - <a class="el" href="pdiffusion_8cpp_source.html">DPG Diffusion example</a>: DPG formulation for the diffusion problem
  * - <a class="el" href="pmaxwell_8cpp_source.html">DPG Maxwell example</a>: DPG formulation for the indefinite Maxwell problem

examples/ex16.cpp

@@ -105,6 +105,7 @@ int main(int argc, char *argv[])
    bool visualization = true;
    bool visit = false;
    int vis_steps = 5;
+   bool solve_implicit_state = false;
 
    int precision = 8;
    cout.precision(precision);
@@ -126,6 +127,9 @@ int main(int argc, char *argv[])
                   "Alpha coefficient.");
    args.AddOption(&kappa, "-k", "--kappa",
                   "Kappa coefficient offset.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -179,6 +183,11 @@ int main(int argc, char *argv[])
 
    // 7. Initialize the conduction operator and the visualization.
    ConductionOperator oper(fespace, alpha, kappa, u);
+   using ImplicitVariableType = ConductionOperator::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   oper.SetImplicitVariableType(imp_var);
 
    u_gf.SetFromTrueDofs(u);
    {
@@ -316,11 +325,14 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 }
 
 void ConductionOperator::ImplicitSolve(const real_t dt,
-                                       const Vector &u, Vector &du_dt)
+                                       const Vector &u, Vector &k)
 {
    // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt, where K is linearized by using u from the previous timestep
+   //    M*k = -K(u + dt*k) for k = du/dt, if solving for stage-slope
+   // or
+   //    M*k = -dt*K(k) + M*u for k = u_s, if solving for stage-state
+   // where K is linearized by using u from the previous timestep, and
+   // the stage-state and slope relation: du/dt = (u_s - u)/dt.
    if (!T)
    {
       T = Add(1.0, Mmat, dt, Kmat);
@@ -328,9 +340,20 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
       T_solver.SetOperator(*T);
    }
    MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u_s
+      Mmat.Mult(u, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      Kmat.Mult(u, z);
+      z.Neg();
+   }
+   T_solver.Mult(z, k);
 }
 
 void ConductionOperator::SetParameters(const Vector &u)

examples/ex16p.cpp

@@ -115,6 +115,7 @@ int main(int argc, char *argv[])
    bool visit = false;
    int vis_steps = 5;
    bool adios2 = false;
+   bool solve_implicit_state = false;
 
    int precision = 8;
    cout.precision(precision);
@@ -138,6 +139,9 @@ int main(int argc, char *argv[])
                   "Alpha coefficient.");
    args.AddOption(&kappa, "-k", "--kappa",
                   "Kappa coefficient offset.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -212,6 +216,11 @@ int main(int argc, char *argv[])
 
    // 9. Initialize the conduction operator and the VisIt visualization.
    ConductionOperator oper(fespace, alpha, kappa, u);
+   using ImplicitVariableType = ConductionOperator::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   oper.SetImplicitVariableType(imp_var);
 
    u_gf.SetFromTrueDofs(u);
    {
@@ -407,11 +416,14 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 }
 
 void ConductionOperator::ImplicitSolve(const real_t dt,
-                                       const Vector &u, Vector &du_dt)
+                                       const Vector &u, Vector &k)
 {
    // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt, where K is linearized by using u from the previous timestep
+   //    M*k = -K(u + dt*k) for k = du/dt, if solving for stage-slope
+   // or
+   //    M*k = -dt*K(k) + M*u for k = u_s, if solving for stage-state
+   // where K is linearized by using u from the previous timestep, and
+   // the stage-state and slope relation: du/dt = (u_s - u)/dt.
    if (!T)
    {
       T = Add(1.0, Mmat, dt, Kmat);
@@ -419,9 +431,20 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
       T_solver.SetOperator(*T);
    }
    MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u
+      Mmat.Mult(u, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      Kmat.Mult(u, z);
+      z.Neg();
+   }
+   T_solver.Mult(z, k);
 }
 
 void ConductionOperator::SetParameters(const Vector &u)

examples/ex9.cpp

@@ -160,6 +160,7 @@ int main(int argc, char *argv[])
    bool paraview = false;
    bool binary = false;
    int vis_steps = 5;
+   bool solve_implicit_state = false;
 
    int precision = 8;
    cout.precision(precision);
@@ -187,6 +188,9 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -366,6 +370,11 @@ int main(int argc, char *argv[])
    //    right-hand side, and perform time-integration (looping over the time
    //    iterations, ti, with a time-step dt).
    FE_Evolution adv(m, k, b);
+   using ImplicitVariableType = FE_Evolution::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   adv.SetImplicitVariableType(imp_var);
 
    real_t t = 0.0;
    adv.SetTime(t);
@@ -459,7 +468,17 @@ void FE_Evolution::ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
 {
    MFEM_VERIFY(dg_solver != NULL,
                "Implicit time integration is not supported with partial assembly");
-   K.Mult(x, z);
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u
+      M.Mult(x, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      K.Mult(x, z);
+   }
    z += b;
    dg_solver->SetTimeStep(dt);
    dg_solver->Mult(z, k);

examples/ex9p.cpp

@@ -257,6 +257,7 @@ int main(int argc, char *argv[])
    bool adios2 = false;
    bool binary = false;
    int vis_steps = 5;
+   bool solve_implicit_state = false;
 #if MFEM_HYPRE_VERSION >= 21800
    PrecType prec_type = PrecType::AIR;
 #else
@@ -290,6 +291,9 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
    args.AddOption((int *)&prec_type, "-pt", "--prec-type", "Preconditioner for "
                   "implicit solves. 0 for ILU, 1 for pAIR-AMG.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -536,6 +540,11 @@ int main(int argc, char *argv[])
    //     right-hand side, and perform time-integration (looping over the time
    //     iterations, ti, with a time-step dt).
    FE_Evolution adv(*m, *k, *B, prec_type);
+   using ImplicitVariableType = FE_Evolution::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   adv.SetImplicitVariableType(imp_var);
 
    real_t t = 0.0;
    adv.SetTime(t);
@@ -676,7 +685,17 @@ FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
 //    (M - dt*K) d = K*u + b
 void FE_Evolution::ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
 {
-   K->Mult(x, z);
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u
+      M->Mult(x, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      K->Mult(x, z);
+   }
    z += b;
    dg_solver->SetTimeStep(dt);
    dg_solver->Mult(z, k);

examples/ginkgo/CMakeLists.txt

@@ -14,6 +14,12 @@ list(APPEND GINKGO_EXAMPLES_SRCS
   ex1.cpp
 )
 
+if (MFEM_USE_MPI AND GINKGO_BUILD_MPI)
+    list(APPEND GINKGO_EXAMPLES_SRCS
+    ex1p.cpp
+)
+endif()
+
 # Include the source directory where mfem.hpp and mfem-performance.hpp are.
 include_directories(BEFORE ${PROJECT_BINARY_DIR})
 

examples/ginkgo/ex1.cpp

@@ -207,7 +207,7 @@ int main(int argc, char *argv[])
             Ginkgo::IcPreconditioner ginkgo_precond(exec, "paric", 30);
             Ginkgo::CGSolver ginkgo_solver(exec, ginkgo_precond);
             ginkgo_solver.SetPrintLevel(print_lvl);
-            ginkgo_solver.SetRelTol(1e-12);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
             ginkgo_solver.SetAbsTol(0.0);
             ginkgo_solver.SetMaxIter(400);
             ginkgo_solver.SetOperator(*(A.Ptr()));
@@ -225,7 +225,7 @@ int main(int argc, char *argv[])
             Ginkgo::MFEMPreconditioner gko_M(exec, M);
             Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
             ginkgo_solver.SetPrintLevel(print_lvl);
-            ginkgo_solver.SetRelTol(1e-12);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
             ginkgo_solver.SetAbsTol(0.0);
             ginkgo_solver.SetMaxIter(400);
             ginkgo_solver.SetOperator(*(A.Ptr()));
@@ -283,7 +283,7 @@ int main(int argc, char *argv[])
                Ginkgo::MFEMPreconditioner gko_M(exec, M);
                Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
                ginkgo_solver.SetPrintLevel(print_lvl);
-               ginkgo_solver.SetRelTol(1e-12);
+               ginkgo_solver.SetRelTol(sqrt(1e-12));
                ginkgo_solver.SetAbsTol(0.0);
                ginkgo_solver.SetMaxIter(400);
                ginkgo_solver.SetOperator(*(A.Ptr()));

examples/ginkgo/ex1p.cpp

@@ -0,0 +1,436 @@
+//                       MFEM Example 1 - Parallel Version
+//                             GINKGO Modification
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex1p -m ../data/star.mesh
+//               mpirun -np 4 ex1p -m ../data/star-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/escher.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
+//               mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
+//               mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
+//               mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
+//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
+//
+// Device sample runs:
+//               mpirun -np 4 ex1p -pa -d cuda
+//               mpirun -np 4 ex1p -fa -d cuda
+//               mpirun -np 4 ex1p -pa -d occa-cuda
+//               mpirun -np 4 ex1p -pa -d raja-omp
+//               mpirun -np 4 ex1p -pa -d ceed-cpu
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
+//             * mpirun -np 4 ex1p -pa -d ceed-cuda
+//             * mpirun -np 4 ex1p -pa -d ceed-hip
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#ifndef MFEM_USE_GINKGO
+#error This example requires that MFEM is built with MFEM_USE_GINKGO=YES
+#endif
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int solver_config = 0;
+   int print_lvl = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&solver_config, "-s", "--solver-config",
+                  "Solver and preconditioner combination: \n\t"
+                  "   0 - Ginkgo solver and Ginkgo preconditioner, \n\t"
+                  "   1 - Ginkgo solver and MFEM preconditioner, \n\t"
+                  "   2 - MFEM solver and Ginkgo preconditioner, \n\t"
+                  "   3 - MFEM solver and MFEM preconditioner.");
+   args.AddOption(&print_lvl, "-pl", "--print-level",
+                  "Print level for iterative solver (1 prints every iteration).");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.SetGPUAwareMPI(true);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   bool delete_fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
+   }
+   else if (pmesh.GetNodes())
+   {
+      fec = pmesh.GetNodes()->OwnFEC();
+      delete_fec = false;
+      if (myid == 0)
+      {
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
+   }
+   ParFiniteElementSpace fespace(&pmesh, fec);
+   HYPRE_BigInt size = fespace.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid
+   //     function corresponding to fespace. Initialize x with initial guess of
+   //     zero, which satisfies the boundary conditions.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the
+   //     Diffusion domain integrator.
+   ParBilinearForm a(&fespace);
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   if (fa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::FULL);
+      // Sort the matrix column indices when running on GPU or with OpenMP (i.e.
+      // when Device::IsEnabled() returns true). This makes the results
+      // bit-for-bit deterministic at the cost of somewhat longer run time.
+      a.EnableSparseMatrixSorting(Device::IsEnabled());
+   }
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a.EnableStaticCondensation(); }
+   a.Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // 13. Solve the linear system A X = B.
+   if (!pa)
+   {
+      switch (solver_config)
+      {
+         // Solve the linear system with CG + Schwarz (with IC) from Ginkgo
+         case 0:
+         {
+            if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
+            Ginkgo::GinkgoExecutor exec(device);
+            Ginkgo::IcPreconditioner local_solver(exec, "exact");
+            Ginkgo::SchwarzPreconditioner gko_M(exec, MPI_COMM_WORLD, local_solver);
+            Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
+            ginkgo_solver.SetPrintLevel(print_lvl);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
+            ginkgo_solver.SetAbsTol(0.0);
+            ginkgo_solver.SetMaxIter(400);
+            ginkgo_solver.SetOperator(*(A.Ptr()));
+            ginkgo_solver.Mult(B, X);
+            break;
+         }
+
+         // Solve the linear system with CG from Ginkgo + MFEM preconditioner
+         case 1:
+         {
+            if (myid == 0) { cout << "Using Ginkgo solver + MFEM preconditioner...\n"; }
+            Ginkgo::GinkgoExecutor exec(device);
+            //Create MFEM preconditioner and wrap it for Ginkgo's use.
+            HypreBoomerAMG M((HypreParMatrix&)(*A));
+            Ginkgo::MFEMPreconditioner gko_M(exec, M, MPI_COMM_WORLD);
+            Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
+            ginkgo_solver.SetPrintLevel(print_lvl);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
+            ginkgo_solver.SetAbsTol(0.0);
+            ginkgo_solver.SetMaxIter(400);
+            ginkgo_solver.SetOperator(*(A.Ptr()));
+            ginkgo_solver.Mult(B, X);
+            break;
+         }
+
+         // Ginkgo Schwarz preconditioner (local ParIC) +  MFEM CG solver
+         case 2:
+         {
+            if (myid == 0) { cout << "Using MFEM solver + Ginkgo preconditioner...\n"; }
+            Ginkgo::GinkgoExecutor exec(device);
+            Ginkgo::IcPreconditioner local_M(exec, "exact");
+            Ginkgo::SchwarzPreconditioner M(exec, MPI_COMM_WORLD, local_M);
+            M.SetOperator(*(A.Ptr()));  // Generate the preconditioner for the matrix A.
+            CGSolver cg(MPI_COMM_WORLD);
+            cg.SetRelTol(sqrt(1e-12));
+            cg.SetMaxIter(400);
+            cg.SetPrintLevel(1);
+            cg.SetPreconditioner(M);
+            cg.SetOperator(*A);
+            cg.Mult(B, X);
+            break;
+         }
+
+         // MFEM solver + MFEM preconditioner
+         case 3:
+         {
+            if (myid == 0) { cout << "Using MFEM solver + MFEM preconditioner...\n"; }
+            HypreBoomerAMG M((HypreParMatrix&)(*A));
+            CGSolver cg(MPI_COMM_WORLD);
+            cg.SetRelTol(sqrt(1e-12));
+            cg.SetMaxIter(400);
+            cg.SetPrintLevel(1);
+            cg.SetPreconditioner(M);
+            cg.SetOperator(*A);
+            cg.Mult(B, X);
+            break;
+         }
+      } // End switch on solver_config
+   }
+   // Partial assembly mode.  Cannot use Ginkgo preconditioners, but can use Ginkgo
+   // solvers.
+   else
+   {
+      if (UsesTensorBasis(fespace))
+      {
+         // Use Jacobi preconditioning in partial assembly mode.
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         switch (solver_config)
+         {
+            case 0:
+            {
+               if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
+               MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
+                          "            Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
+               break;
+            }
+
+            // Use Ginkgo solver with MFEM preconditioner
+            case 1:
+            {
+               if (myid == 0) { cout << "Using Ginkgo solver + MFEM preconditioner...\n"; }
+               Ginkgo::GinkgoExecutor exec(device);
+               // Wrap MFEM preconditioner for Ginkgo's use.
+               Ginkgo::MFEMPreconditioner gko_M(exec, M, MPI_COMM_WORLD);
+               Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
+               ginkgo_solver.SetPrintLevel(print_lvl);
+               ginkgo_solver.SetRelTol(sqrt(1e-12));
+               ginkgo_solver.SetAbsTol(0.0);
+               ginkgo_solver.SetMaxIter(400);
+               ginkgo_solver.SetOperator(*(A.Ptr()));
+               ginkgo_solver.Mult(B, X);
+               break;
+            }
+
+            // No Ginkgo preconditioners work with matrix-free; error
+            case 2:
+            {
+               if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
+               MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
+                          "            Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
+               break;
+            }
+
+            // Use MFEM solver and preconditioner
+            case 3:
+            {
+               if (myid == 0) { cout << "Using MFEM solver + MFEM preconditioner...\n"; }
+               CGSolver cg(MPI_COMM_WORLD);
+               cg.SetRelTol(sqrt(1e-12));
+               cg.SetMaxIter(400);
+               cg.SetPrintLevel(1);
+               cg.SetPreconditioner(M);
+               cg.SetOperator(*A);
+               cg.Mult(B, X);
+               break;
+            }
+         } // End switch on solver_config
+      }
+      else // CG with no preconditioning
+      {
+         if (myid == 0) { cout << "Using MFEM solver + no preconditioner...\n"; }
+         CGSolver cg(MPI_COMM_WORLD);
+         cg.SetRelTol(sqrt(1e-12));
+         cg.SetMaxIter(400);
+         cg.SetPrintLevel(1);
+         cg.SetOperator(*A);
+         cg.Mult(B, X);
+      }
+   }
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a.RecoverFEMSolution(X, b, x);
+
+   // 15. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh.Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   if (delete_fec)
+   {
+      delete fec;
+   }
+
+   return 0;
+}

examples/ginkgo/makefile

@@ -20,9 +20,8 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-# Currently there are only serial Ginkgo examples
 SEQ_EXAMPLES = ex1
-PAR_EXAMPLES =
+PAR_EXAMPLES = ex1p
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
 else

fem/dfem/doperator.hpp

@@ -968,12 +968,17 @@ void DifferentiableOperator::AddIntegrator(
          // Quadrature point local derivative cache for each element, with data
          // layout:
          // [test_vdim, test_op_dim, trial_vdim, trial_op_dim, qp, num_entities].
-         derivative_qp_caches[derivative_id] = Vector(test_vdim * test_op_dim *
-                                                      trial_vdim *
-                                                      total_trial_op_dim * num_qp * num_entities);
+         size_t derivative_qp_size = test_vdim * test_op_dim *
+                                     trial_vdim *
+                                     total_trial_op_dim * num_qp * num_entities;
+         if (derivative_qp_caches.count(derivative_qp_size) == 0)
+         {
+            derivative_qp_caches[derivative_qp_size] = Vector(derivative_qp_size);
+         }
+
          // Create local references for MSVC lambda capture compatibility
          auto& fields_ref = this->fields;
-         auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_id];
+         auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_qp_size];
 
          // In each of the callbacks we're saving the derivatives in the quadrature point
          // caches. This trades memory with computational effort but also minimizes

fem/dfem/tuple.hpp

@@ -10,7 +10,7 @@
 // CONTRIBUTING.md for details.
 #pragma once
 
-// This is serac's tuple implementation
+// This is smith's tuple implementation
 
 #include <ostream>
 #include "../../config/config.hpp"

fem/lor/lor_dg.hpp

@@ -23,6 +23,8 @@ class BatchedLOR_DG : BatchedLORKernel
 {
    IntegrationRule ir_face; ///< Collocated Gauss-Lobatto face quadrature rule.
    real_t kappa; ///< DG penalty parameter.
+   bool has_bdr_integ; ///< Is there a boundary integrator?
+   const Array<int> *bdr_markers; ///< Boundary integrator markers.
 public:
    template <int ORDER, int SDIM> void Assemble2D();
    template <int ORDER> void Assemble3D();
@@ -38,8 +40,7 @@ public:
       ProjectLORCoefficient<MassIntegrator>(a, c1);
       ProjectLORCoefficient<DiffusionIntegrator>(a, c2);
 
-      auto *integ = GetInteriorFaceIntegrator<DGDiffusionIntegrator>(a);
-      if (integ)
+      if (auto *integ = GetInteriorFaceIntegrator<DGDiffusionIntegrator>(a))
       {
          kappa = integ->GetPenaltyParameter();
       }
@@ -47,6 +48,19 @@ public:
       {
          kappa = 0.0;
       }
+
+      has_bdr_integ = false;
+      auto *bdr_face_integs = a.GetBFBFI();
+      for (int i = 0; i < bdr_face_integs->Size(); ++i)
+      {
+         if (auto *integ = dynamic_cast<DGDiffusionIntegrator*>((*bdr_face_integs)[i]))
+         {
+            kappa = integ->GetPenaltyParameter();
+            bdr_markers = (*a.GetBFBFI_Marker())[i];
+            has_bdr_integ = true;
+            break;
+         }
+      }
    }
 
    /// @brief Compute and return the face info array.

fem/lor/lor_dg_impl.hpp

@@ -22,9 +22,13 @@ namespace mfem
 Array<int> BatchedLOR_DG::GetFaceInfo() const
 {
    Mesh &mesh = *fes_ho.GetMesh();
+   const Array<int> &bdr_face_attrs = mesh.GetBdrFaceAttributes();
    const int nf = mesh.GetNumFaces();
    Array<int> face_info(nf * 6); // (e0, f0, o0, e1, f1, o1)
    auto h_face_info = Reshape(face_info.HostWrite(), 6, nf);
+
+   int bdr_face_counter = 0;
+
    for (int f = 0; f < nf; ++f)
    {
       auto finfo = mesh.GetFaceInformation(f);
@@ -43,6 +47,19 @@ Array<int> BatchedLOR_DG::GetFaceInfo() const
          h_face_info(4, f) = -1;
          h_face_info(5, f) = -1;
       }
+
+      if (finfo.IsBoundary())
+      {
+         // Check if Neumann boundary; skip these when adding boundary penalties
+         const int bdr_attr = bdr_face_attrs[bdr_face_counter];
+         if (!has_bdr_integ || (bdr_markers && !(*bdr_markers)[bdr_attr - 1]))
+         {
+            h_face_info(0, f) = -1;
+            h_face_info(1, f) = -1;
+            h_face_info(2, f) = -1;
+         }
+         bdr_face_counter += 1;
+      }
    }
    return face_info;
 }
@@ -144,6 +161,7 @@ void BatchedLOR_DG::AssembleFaceTerms()
    {
       const int f_0 = d_face_info(1, f);
       const int f_1 = d_face_info(4, f);
+      if (f_0 < 0) { return; } // Skip Neumann boundary faces
       const int nsides = (f_1 >= 0) ? 2 : 1;
       for (int el_i = 0; el_i < nsides; ++el_i)
       {

general/device.cpp

@@ -22,6 +22,8 @@
 
 #include <unordered_map>
 #include <map>
+#include <sstream>
+#include <iomanip>
 
 namespace mfem
 {
@@ -716,6 +718,29 @@ void Device::DeviceMem(size_t *free, size_t *total)
 #endif
 }
 
+std::string Device::GetUUID(const int device_id)
+{
+   std::stringstream res;
+#if defined(MFEM_USE_CUDA)
+   cudaDeviceProp prop;
+   MFEM_GPU_CHECK(cudaGetDeviceProperties(&prop, device_id));
+   for (int i = 0; i < 16; ++i)
+   {
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(prop.uuid.bytes[i]);
+   }
+#elif defined(MFEM_USE_HIP)
+   hipUUID uuid;
+   MFEM_GPU_CHECK(hipDeviceGetUuid(&uuid, device_id));
+   for (int i = 0; i < 16; ++i)
+   {
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(uuid.bytes[i]);
+   }
+#endif
+   return res.str();
+}
+
 int Device::NumMultiprocessors(int dev)
 {
 #if defined(MFEM_USE_CUDA)

general/device.hpp

@@ -255,6 +255,10 @@ public:
    /// Get the number of available devices (may be called before configuration).
    static int GetDeviceCount();
 
+   /// Gets a string representation of the GPU UUID.
+   /// 0 <= @a device_id < GetDeviceCount()
+   static std::string GetUUID(const int device_id = 0);
+
    /** @brief Return true if any of the backends in the backend mask, @a b_mask,
        are allowed. */
    /** This method can be used with any of the Backend::Id constants, the

general/table.hpp

@@ -146,7 +146,7 @@ public:
    int *ReadWriteJ(bool on_dev = true) { return J.ReadWrite(on_dev); }
    const int *HostReadJ() const { return J.HostRead(); }
    int *HostWriteJ() { return J.HostWrite(); }
-   int *ReadWriteJ() { return J.HostReadWrite(); }
+   int *HostReadWriteJ() { return J.HostReadWrite(); }
 
    /// Sort the column (TYPE II) indices in each row.
    void SortRows();

linalg/hypre.cpp

@@ -1681,6 +1681,13 @@ void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_BigInt* &cmap) const
    cmap = A->col_map_offd;
 }
 
+void HypreParMatrix::GetOffdColMap(HYPRE_BigInt* &cmap,
+                                   HYPRE_Int &num_cols) const
+{
+   cmap = A->col_map_offd;
+   num_cols = hypre_CSRMatrixNumCols(A->offd);
+}
+
 void HypreParMatrix::MergeDiagAndOffd(SparseMatrix &merged)
 {
    HostRead();

linalg/hypre.hpp

@@ -665,6 +665,8 @@ public:
    void GetDiag(SparseMatrix &diag) const;
    /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
    void GetOffd(SparseMatrix &offd, HYPRE_BigInt* &cmap) const;
+   /// Get the global column mapping for the local off-diagonal block.
+   void GetOffdColMap(HYPRE_BigInt* &cmap, HYPRE_Int &num_cols) const;
    /** @brief Get a single SparseMatrix containing all rows from this processor,
        merged from the diagonal and off-diagonal blocks stored by the
        HypreParMatrix. */
@@ -959,6 +961,14 @@ public:
    const Memory<HYPRE_Int> &GetDiagMemoryJ() const { return mem_diag.J; }
    const Memory<real_t> &GetDiagMemoryData() const { return mem_diag.data; }
 
+   Memory<HYPRE_Int> &GetOffdMemoryI() { return mem_offd.I; }
+   Memory<HYPRE_Int> &GetOffdMemoryJ() { return mem_offd.J; }
+   Memory<real_t> &GetOffdMemoryData() { return mem_offd.data; }
+
+   const Memory<HYPRE_Int> &GetOffdMemoryI() const { return mem_offd.I; }
+   const Memory<HYPRE_Int> &GetOffdMemoryJ() const { return mem_offd.J; }
+   const Memory<real_t> &GetOffdMemoryData() const { return mem_offd.data; }
+
    /// @brief Prints the locally owned rows in parallel. The resulting files can
    /// be read with Read_IJMatrix().
    void Print(const std::string &fname, HYPRE_Int offi = 0,

linalg/ode.cpp

@@ -10,6 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #include "../general/communication.hpp"
+#include "../general/forall.hpp"
 #include "operator.hpp"
 #include "ode.hpp"
 
@@ -184,6 +185,23 @@ void ODESolver::Init(TimeDependentOperator &f_)
    mem_type = GetMemoryType(f_.GetMemoryClass());
 }
 
+void ODESolver::ComputeSlopeFromState(const real_t dt, const Vector &u,
+                                      Vector &k)
+{
+   // k currently holds state u(t+dt),
+   // convert to slope k = du/dt ~= (u(t+dt)-u(t))/dt
+   const int usz = u.Size();
+   real_t fac = 1.0/dt;
+   auto d_u = u.Read();
+   auto d_k = k.ReadWrite();
+
+   mfem::forall(usz, [=] MFEM_HOST_DEVICE (int i)
+   {
+      d_k[i] -= d_u[i];
+      d_k[i] *= fac;
+   });
+}
+
 void ForwardEulerSolver::Init(TimeDependentOperator &f_)
 {
    ODESolver::Init(f_);
@@ -629,6 +647,10 @@ void AdamsMoultonSolver::Step(Vector &x, real_t &t, real_t &dt)
       }
       state.ShiftStages();
       f->ImplicitSolve(a[0]*dt, x, state[0]);
+      if (f->ImplicitVarTypeIsState())
+      {
+         ComputeSlopeFromState(a[0]*dt, x, state[0]);
+      }
       x.Add(a[0]*dt, state[0]);
       t += dt;
    }
@@ -661,7 +683,15 @@ void BackwardEulerSolver::Step(Vector &x, real_t &t, real_t &dt)
 {
    f->SetTime(t + dt);
    f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
-   x.Add(dt, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      x = k; // x = u_{i+1}
+   }
+   else
+   {
+      x.Add(dt, k);
+   }
+
    t += dt;
 }
 
@@ -676,7 +706,16 @@ void ImplicitMidpointSolver::Step(Vector &x, real_t &t, real_t &dt)
 {
    f->SetTime(t + dt/2);
    f->ImplicitSolve(dt/2, x, k);
-   x.Add(dt, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      x.Neg();
+      x.Add(2.0, k);
+   }
+   else
+   {
+      x.Add(dt, k);
+   }
+
    t += dt;
 }
 
@@ -718,11 +757,19 @@ void SDIRK23Solver::Step(Vector &x, real_t &t, real_t &dt)
    // note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
    f->SetTime(t + gamma*dt);
    f->ImplicitSolve(gamma*dt, x, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(gamma*dt, x, k);
+   }
    add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
    x.Add(dt/2, k);
 
    f->SetTime(t + (1.-gamma)*dt);
    f->ImplicitSolve(gamma*dt, y, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(gamma*dt, y, k);
+   }
    x.Add(dt/2, k);
    t += dt;
 }
@@ -749,17 +796,29 @@ void SDIRK34Solver::Step(Vector &x, real_t &t, real_t &dt)
 
    f->SetTime(t + a*dt);
    f->ImplicitSolve(a*dt, x, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, x, k);
+   }
    add(x, (0.5-a)*dt, k, y);
    add(x,  (2.*a)*dt, k, z);
    x.Add(b*dt, k);
 
    f->SetTime(t + dt/2);
    f->ImplicitSolve(a*dt, y, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, y, k);
+   }
    z.Add((1.-4.*a)*dt, k);
    x.Add((1.-2.*b)*dt, k);
 
    f->SetTime(t + (1.-a)*dt);
    f->ImplicitSolve(a*dt, z, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, z, k);
+   }
    x.Add(b*dt, k);
    t += dt;
 }
@@ -785,15 +844,27 @@ void SDIRK33Solver::Step(Vector &x, real_t &t, real_t &dt)
 
    f->SetTime(t + a*dt);
    f->ImplicitSolve(a*dt, x, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, x, k);
+   }
    add(x, (c-a)*dt, k, y);
    x.Add(b*dt, k);
 
    f->SetTime(t + c*dt);
    f->ImplicitSolve(a*dt, y, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, y, k);
+   }
    x.Add((1.0-a-b)*dt, k);
 
    f->SetTime(t + dt);
    f->ImplicitSolve(a*dt, x, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, x, k);
+   }
    x.Add(a*dt, k);
    t += dt;
 }
@@ -818,6 +889,10 @@ void TrapezoidalRuleSolver::Step(Vector &x, real_t &t, real_t &dt)
 
    f->SetTime(t + dt);
    f->ImplicitSolve(dt/2.0, y, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(0.5*dt, y, k);
+   }
    x.Add(dt/2.0, k);
    t += dt;
 }
@@ -848,11 +923,19 @@ void ESDIRK32Solver::Step(Vector &x, real_t &t, real_t &dt)
 
    f->SetTime(t + (2.0*a)*dt);
    f->ImplicitSolve(a*dt, y, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, y, k);
+   }
    z.Add(b*dt, k);
    x.Add(b*dt, k);
 
    f->SetTime(t + dt);
    f->ImplicitSolve(a*dt, z, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, z, k);
+   }
    x.Add(a*dt, k);
    t += dt;
 }
@@ -885,11 +968,19 @@ void ESDIRK33Solver::Step(Vector &x, real_t &t, real_t &dt)
 
    f->SetTime(t + (2.0*a)*dt);
    f->ImplicitSolve(a*dt, y, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, y, k);
+   }
    z.Add(b*dt, k);
    x.Add(b_2*dt, k);
 
    f->SetTime(t + dt);
    f->ImplicitSolve(a*dt, z, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(a*dt, z, k);
+   }
    x.Add(b_3*dt, k);
    t += dt;
 }
@@ -955,6 +1046,10 @@ void GeneralizedAlphaSolver::Step(Vector &x, real_t &t, real_t &dt)
    real_t dt_eff = (gamma*alpha_f/alpha_m)*dt;
    f->SetTime(t + alpha_f*dt);
    f->ImplicitSolve(dt_eff, y, k);
+   if (f->ImplicitVarTypeIsState())
+   {
+      ComputeSlopeFromState(dt_eff, y, k);
+   }
 
    // Update x and xdot
    x.Add((1.0 - (gamma/alpha_m))*dt, state[0]);
@@ -1116,8 +1211,8 @@ void SecondOrderODESolver::EulerStep(Vector &x, Vector &dxdt, real_t &t,
    f->SetTime(t + dt);
    f->ImplicitSolve(0.5*dt*dt, dt, x, dxdt, state[0]);
 
-   x   .Add(0.5*dt*dt, state[0]);
-   dxdt.Add(dt,    state[0]);
+   x.Add(0.5*dt*dt, state[0]);
+   dxdt.Add(dt, state[0]);
    t += dt;
 }
 
@@ -1203,8 +1298,8 @@ void NewmarkSolver::Step(Vector &x, Vector &dxdt, real_t &t, real_t &dt)
    f->SetTime(t + dt);
    f->ImplicitSolve(fac3*dt*dt, fac4*dt, x, dxdt, state[0]);
 
-   x   .Add(fac3*dt*dt, state[0]);
-   dxdt.Add(fac4*dt,    state[0]);
+   x.Add(fac3*dt*dt, state[0]);
+   dxdt.Add(fac4*dt, state[0]);
    t += dt;
 }
 

linalg/ode.hpp

@@ -120,6 +120,7 @@ public:
 class ODESolver
 {
 protected:
+   using ImplicitVariableType = TimeDependentOperator::ImplicitVariableType;
    /// Pointer to the associated TimeDependentOperator.
    TimeDependentOperator *f;  // f(.,t) : R^n --> R^n
    MemoryType mem_type;
@@ -192,6 +193,22 @@ public:
    /// Returns how many State vectors the ODE requires
    virtual int GetStateSize() { return 0; };
 
+   ///@brief Returns @a true if the ODESolver supports the given
+   /// #ImplicitVariableType, @a var, and returns @a false otherwise.
+   ///@note Should be overriden in ODESolver that calls TimeDependentOperator::ImplicitSolve().
+   virtual bool SupportsImplicitVariableType(ImplicitVariableType var) const
+   { return false; };
+
+   /** @brief Compute the finite-difference slope, @a $\frac{du}{dt} \approx \frac{u(t+dt)-u(t)}{dt}$,
+    * and store it in @a k.
+    * @param [in] dt Finite difference step size.
+    * @param [in] u  state vector, @a u(t).
+    * @param [in,out] k   On input, @a k contains the state vector, @a u( @a t+ @a dt).
+    * On output, @a k contains the computed slope, @a du/dt.
+    * */
+   virtual void ComputeSlopeFromState(const real_t dt, const Vector &u,
+                                      Vector &k);
+
    // Help info for ODESolver options
    static MFEM_EXPORT std::string ExplicitTypes;
    static MFEM_EXPORT std::string ImplicitTypes;
@@ -361,6 +378,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -374,6 +397,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -395,6 +424,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -409,6 +444,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -423,6 +464,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -437,6 +484,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -451,6 +504,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -465,6 +524,12 @@ public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -490,6 +555,12 @@ public:
 
    ODEStateData& GetState() override { return state; }
    const ODEStateData& GetState() const override { return state; }
+
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 
@@ -606,6 +677,11 @@ public:
 
    ODEStateData& GetState() override { return state; }
    const ODEStateData& GetState() const override { return state; }
+   bool SupportsImplicitVariableType(ImplicitVariableType var) const override
+   {
+      return (var == ImplicitVariableType::STATE ||
+              var == ImplicitVariableType::SLOPE);
+   }
 };
 
 /** A 1-stage, 2nd order AM method. */

linalg/operator.hpp

@@ -381,11 +381,24 @@ public:
       ADDITIVE_TERM_2
    };
 
+   /** Used to specify the variable being returned by ImplicitSolve(). This can
+    * be queried by ODESolver to identify the variable being solved for.
+    * @warning Not all ODESolver may support all options. See ODESolver::SupportsImplicitVariableType() */
+   enum ImplicitVariableType
+   {
+      SLOPE, ///< stage slope, $k = \frac{du}{dt}$.
+      STATE  ///< stage state, $k = u$.
+   };
+
 protected:
    real_t t;  ///< Current time.
    Type type; /**< @brief Describes the form of the TimeDependentOperator, see
                    the documentation of #Type. */
    EvalMode eval_mode; ///< Current evaluation mode.
+   ImplicitVariableType implicit_variable_type =
+      ImplicitVariableType::SLOPE; /**< @brief
+                                                        Return variable for
+                                                        ImplicitSolve()*/
 
 public:
    /** @brief Construct a "square" TimeDependentOperator (u,t) -> k(u,t), where
@@ -429,6 +442,24 @@ public:
    virtual void SetEvalMode(const EvalMode new_eval_mode)
    { eval_mode = new_eval_mode; }
 
+   /** @brief Sets the #ImplicitVariableType for ImplicitSolve()*/
+   virtual void SetImplicitVariableType(const ImplicitVariableType variable_type)
+   { implicit_variable_type = variable_type; }
+
+   /** @brief Returns the #ImplicitVariableType for ImplicitSolve(). */
+   virtual ImplicitVariableType GetImplicitVariableType() const
+   { return implicit_variable_type; }
+
+   /** @brief Returns @a true if implicit variable is #STATE and @a false otherwise.
+    * Used by ODESolver to identify the stage variable returned by ImplicitSolve() */
+   virtual bool ImplicitVarTypeIsState() const
+   { return (implicit_variable_type == ImplicitVariableType::STATE); }
+
+   /** @brief Returns @a true if implicit variable is #SLOPE and @a false otherwise.
+    * Used by ODESolver to identify the stage variable returned by ImplicitSolve() */
+   virtual bool ImplicitVarTypeIsSlope() const
+   { return (implicit_variable_type == ImplicitVariableType::SLOPE); }
+
    /** @brief Perform the action of the explicit part of the operator, G:
        @a v = G(@a u, t) where t is the current time.
 
@@ -462,7 +493,8 @@ public:
 
    /** @brief Solve for the unknown @a k, at the current time t, the following
        equation:
-       F(@a u + @a gamma @a k, @a k, t) = G(@a u + @a gamma @a k, t).
+        1. $F( u + \gamma k, k, t) = G( u +  \gamma k, t)$, if solving for stage-slope (default)
+        2. $F( u , \frac{k-u}{\gamma}, t) = G(k, t)$, if solving for stage-state
 
        For solving an ordinary differential equation of the form
        $ M \frac{dy}{dt} = g(y,t) $, recall that F and G can be defined in
@@ -472,8 +504,9 @@ public:
          2. F(u,k,t) = M k and G(u,t) = g(u,t)
          3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0
 
-       Regardless of the choice of F and G, this function should solve for @a k
-       in M @a k = g(@a u + @a gamma @a k, t).
+        Regardless of the choice of F and G, this function should solve for @a k:
+        - $~Mk = g( u + \gamma k, t)~$, if solving for stage-slope.
+        - $~Mk  = \gamma g(k, t) + Mu~$, if solving for stage-state
 
        To see how @a k can be useful, consider the backward Euler method defined
        by $ y(t + \Delta t) = y(t) + \Delta t k_0 $ where
@@ -491,6 +524,7 @@ public:
        $ y(t) + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j $ and @a gamma set to
        $ a_{ii} \Delta t $, for $ k_i $. For example, see class SDIRK33Solver.
 
+       See SetImplicitVariableType() to switch between different variable modes.
        If not re-implemented, this method simply generates an error. */
    virtual void ImplicitSolve(const real_t gamma, const Vector &u, Vector &k);
 

linalg/sparsesmoothers.cpp

@@ -23,15 +23,31 @@ namespace mfem
 void SparseSmoother::SetOperator(const Operator &a)
 {
    oper = dynamic_cast<const SparseMatrix*>(&a);
-   if (oper == NULL)
-   {
-      mfem_error("SparseSmoother::SetOperator : not a SparseMatrix!");
-   }
+   MFEM_VERIFY(oper != nullptr, "Operator must be a SparseMatrix");
    height = oper->Height();
    width = oper->Width();
+
+   At.reset();
+   oper_T = nullptr;
+}
+
+void SparseSmoother::EnsureTranspose() const
+{
+   if (oper_T) { return; }
+
+   const real_t tol = 1e-14;
+   if (oper->IsSymmetric() > tol * oper->MaxNorm())
+   {
+      At.reset(Transpose(*oper));
+      oper_T = At.get();
+   }
+   else
+   {
+      At.reset();
+      oper_T = oper;
+   }
 }
 
-/// Matrix vector multiplication with GS Smoother.
 void GSSmoother::Mult(const Vector &x, Vector &y) const
 {
    if (!iterative_mode)
@@ -51,21 +67,33 @@ void GSSmoother::Mult(const Vector &x, Vector &y) const
    }
 }
 
-/// Create the Jacobi smoother.
-DSmoother::DSmoother(const SparseMatrix &a, int t, real_t s, int it)
-   : SparseSmoother(a)
+void GSSmoother::MultTranspose(const Vector &x, Vector &y) const
 {
-   type = t;
-   scale = s;
-   iterations = it;
+   EnsureTranspose();
+
+   if (!iterative_mode)
+   {
+      y = 0.0;
+   }
+
+   for (int i = 0; i < iterations; i++)
+   {
+      if (type != 1)
+      {
+         oper_T->Gauss_Seidel_forw(x, y);
+      }
+      if (type != 2)
+      {
+         oper_T->Gauss_Seidel_back(x, y);
+      }
+   }
 }
 
-/// Matrix vector multiplication with Jacobi smoother.
-void DSmoother::Mult(const Vector &x, Vector &y) const
+void DSmoother::Mult_(const SparseMatrix &A, const Vector &x, Vector &y) const
 {
    if (!iterative_mode && type == 0 && iterations == 1)
    {
-      oper->DiagScale(x, y, scale, use_abs_diag);
+      A.DiagScale(x, y, scale, use_abs_diag);
       return;
    }
 
@@ -90,22 +118,41 @@ void DSmoother::Mult(const Vector &x, Vector &y) const
    {
       if (type == 0)
       {
-         oper->Jacobi(x, *p, *r, scale, use_abs_diag);
+         A.Jacobi(x, *p, *r, scale, use_abs_diag);
       }
       else if (type == 1)
       {
-         oper->Jacobi2(x, *p, *r, scale);
+         A.Jacobi2(x, *p, *r, scale);
       }
       else if (type == 2)
       {
-         oper->Jacobi3(x, *p, *r, scale);
+         A.Jacobi3(x, *p, *r, scale);
       }
       else
       {
-         mfem_error("DSmoother::Mult wrong type");
+         MFEM_ABORT("Invalid type.");
       }
       Swap<Vector*>(r, p);
    }
 }
 
+void DSmoother::Mult(const Vector &x, Vector &y) const
+{
+   Mult_(*oper, x, y);
+}
+
+void DSmoother::MultTranspose(const Vector &x, Vector &y) const
+{
+   if (iterations == 1 && !iterative_mode)
+   {
+      Mult_(*oper, x, y);
+      return;
+   }
+
+   EnsureTranspose();
+   MFEM_VERIFY(type == 0 || !At, "l1 or lumped Jacobi transpose not implemented"
+               " for non-symmetric matrices");
+   Mult_(*oper_T, x, y);
+}
+
 }

linalg/sparsesmoothers.hpp

@@ -15,67 +15,159 @@
 #include "../config/config.hpp"
 #include "sparsemat.hpp"
 
+#include <memory>
+
 namespace mfem
 {
 
+/// Abstract base class for smoothers created from a SparseMatrix.
 class SparseSmoother : public MatrixInverse
 {
 protected:
-   const SparseMatrix *oper;
+   const SparseMatrix *oper = nullptr; ///< The underlying matrix.
+
+   /// Pointer to the transpose of the underlying matrix. If the matrix is
+   /// symmetric, this will be the same as @a oper. If the matrix is not
+   /// symmetric, the transpose will be formed and stored in @a At. The
+   /// transpose will only be formed if MultTranspose() is called.
+   mutable const SparseMatrix *oper_T = nullptr;
+
+   mutable std::unique_ptr<SparseMatrix> At; ///< Transpose of A, if needed.
+
+   void EnsureTranspose() const; ///< Ensure that the transpose is set.
 
 public:
-   SparseSmoother() { oper = NULL; }
+   SparseSmoother() = default;
 
-   SparseSmoother(const SparseMatrix &a)
-      : MatrixInverse(a) { oper = &a; }
+   SparseSmoother(const SparseMatrix &a) { SetOperator(a); }
 
+   /// Sets the underlying matrix. @a a must be a SparseMatrix.
    void SetOperator(const Operator &a) override;
 };
 
-/// Data type for Gauss-Seidel smoother of sparse matrix
+/// Gauss-Seidel smoother of a sparse matrix.
 class GSSmoother : public SparseSmoother
 {
+public:
+   enum GSType
+   {
+      SYMMETRIC, ///< Forward Gauss-Seidel, then backward.
+      FORWARD, ///< Forward Gauss-Seidel ($L^{-1}$).
+      BACKWARD ///< Backward Gauss-Seidel ($U^{-1}$).
+   };
 protected:
-   int type; // 0, 1, 2 - symmetric, forward, backward
-   int iterations;
+   GSType type; ///< Type of Gauss-Seidel, see GSSmoother::GSType.
+   int iterations; ///< Number of stationary iterations.
 
 public:
-   /// Create GSSmoother.
-   GSSmoother(int t = 0, int it = 1) { type = t; iterations = it; }
+   /// @brief Create a Gauss-Seidel smoother. SetOperator() will need to be
+   /// called with a SparseMatrix before first use.
+   ///
+   /// @param[in]  t        Type of GS smoother (see GSSmoother::GSType)
+   /// @param[in]  it       Number of stationary iterations to perform
+   GSSmoother(GSType t = SYMMETRIC, int it = 1) { type = t; iterations = it; }
 
-   /// Create GSSmoother.
-   GSSmoother(const SparseMatrix &a, int t = 0, int it = 1)
-      : SparseSmoother(a) { type = t; iterations = it; }
+   /// @brief Create a Jacobi smoother using the SparseMatrix @a a.
+   ///
+   /// @param[in]  a        The underlying SparseMatrix
+   /// @param[in]  t        Type of GS smoother (see GSSmoother::GSType)
+   /// @param[in]  it       Number of stationary iterations to perform
+   GSSmoother(const SparseMatrix &a, GSType t = SYMMETRIC, int it = 1)
+      : GSSmoother(t, it) { SetOperator(a); }
 
-   /// Matrix vector multiplication with GS Smoother.
+   /// Same as GSSmoother(GSType,int), for backwards compatibility.
+   GSSmoother(int t, int it = 1) : GSSmoother(GSType(t), it) { }
+
+   /// @brief Same as GSSmoother(const SparseMatrix&,GSType,int), for
+   /// backwards compatibility.
+   GSSmoother(const SparseMatrix &a, int t, int it = 1)
+      : GSSmoother(a, GSType(t), it) { }
+
+   /// @brief Application of the Gauss-Seidel smoother.
+   ///
+   /// Applies a stationary Gauss-Seidel iteration. If Solver::iterative_mode is
+   /// true, then @a y is used as the initial guess, and Gauss-Seidel is applied
+   /// to the residual $x - Ay$.
    void Mult(const Vector &x, Vector &y) const override;
+
+   /// Application of the transpose of the Gauss-Seidel smoother.
+   void MultTranspose(const Vector &x, Vector &y) const override;
 };
 
-/// Data type for scaled Jacobi-type smoother of sparse matrix
+/// Jacobi-type diagonal smoother of a sparse matrix.
 class DSmoother : public SparseSmoother
 {
+public:
+   enum JacobiType
+   {
+      JACOBI,       ///< Scale by the diagonal of the matrix.
+      L1_JACOBI,    ///< Scale by the l1-norm of the rows.
+      LUMPED_JACOBI ///< Scale by the sum of the rows.
+   };
 protected:
-   int type; // 0, 1, 2 - scaled Jacobi, scaled l1-Jacobi, scaled lumped-Jacobi
-   real_t scale;
-   int iterations;
-   /// Uses abs values of the diagonal entries. Relevant only when type = 0.
+   JacobiType type; ///< Type of diagonal scaling, see DSmoother::JacobiType.
+   real_t scale; ///< Scaling (damping) factor.
+   int iterations; ///< Number of stationary iterations to perform.
+
+   /// @brief Uses abs values of the diagonal entries. Relevant only with type
+   /// JacobiType::JACOBI.
    bool use_abs_diag = false;
 
-   mutable Vector z;
+   mutable Vector z; ///< Temporary work vector.
+
+   /// Apply the Jacobi smoother (used internally by Mult() and MultTranspose())
+   void Mult_(const SparseMatrix &A, const Vector &x, Vector &y) const;
 
 public:
-   /// Create Jacobi smoother.
-   DSmoother(int t = 0, real_t s = 1., int it = 1)
+   /// @brief Create a Jacobi smoother. SetOperator() will need to be called
+   /// with a SparseMatrix before first use.
+   ///
+   /// @param[in]  t        Type of Jacobi smoother (see DSmoother::JacobiType)
+   /// @param[in]  s        Scaling factor
+   /// @param[in]  it       Number of stationary iterations to perform
+   DSmoother(JacobiType t = JACOBI, real_t s = 1., int it = 1)
    { type = t; scale = s; iterations = it; }
 
-   /// Create Jacobi smoother.
-   DSmoother(const SparseMatrix &a, int t = 0, real_t s = 1., int it = 1);
+   /// @brief Create a Jacobi smoother using the SparseMatrix @a a.
+   ///
+   /// @param[in]  a        The underlying SparseMatrix
+   /// @param[in]  t        Type of Jacobi smoother (see DSmoother::JacobiType)
+   /// @param[in]  s        Scaling factor
+   /// @param[in]  it       Number of stationary iterations to perform
+   DSmoother(const SparseMatrix &a, JacobiType t = JACOBI, real_t s = 1.,
+             int it = 1) : DSmoother(t, s, it) { SetOperator(a); }
 
-   /// Replace diag entries with their abs values. Relevant only when type = 0.
+   /// @brief Same as DSmoother(JacobiType,real_t,int), for backwards compatbility.
+   DSmoother(int t, real_t s = 1., int it = 1)
+      : DSmoother(JacobiType(t), s, it) { }
+
+   /// @brief Same as DSmoother(const SparseMatrix&,JacobiType,real_t,int), for
+   /// backwards compatbility.
+   DSmoother(const SparseMatrix &a, int t, real_t s = 1., int it = 1)
+      : DSmoother(a, JacobiType(t), s, it) { }
+
+   /// @brief Replace diagonal entries with their absolute values. Relevant only
+   /// with JacobiType::JACOBI.
    void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
 
-   /// Matrix vector multiplication with Jacobi smoother.
+   /// @brief Apply the Jacobi smoother.
+   ///
+   /// Applies a stationary iteration with diagonal scaling. If
+   /// Solver::iterative_mode is true, then @a y is used as the initial guess
+   /// (and the diagonal scaling is applied to the residual $x - Ay$, giving
+   /// $D^{-1}(x - Ay)$).
+   ///
+   /// By default, Solver::iterative_mode is false and only one iteration is
+   /// performed, corresponding to $y = D^{-1}x$.
    void Mult(const Vector &x, Vector &y) const override;
+
+   /// @brief Apply the transpose of the Jacobi smoother.
+   ///
+   /// If the underlying matrix is symmetric, or if only one iteration is
+   /// performed with zero initial guess (Solver::iterative_mode is false), then
+   /// this is the same as Mult(). For non-symmetric matrices with iteration
+   /// count greater than one, only JacobiType::JACOBI is supported.
+   void MultTranspose(const Vector &x, Vector &y) const override;
 };
 
 }

mesh/mesh.cpp

@@ -36,6 +36,7 @@
 #include <numeric>
 #include <unordered_map>
 #include <unordered_set>
+#include <list>
 
 // Include the METIS header, if using version 5. If using METIS 4, the needed
 // declarations are inlined below, i.e. no header is needed.
@@ -4772,12 +4773,12 @@ Mesh::Mesh(real_t *vertices_, int num_vertices,
    FinalizeTopology();
 }
 
-Mesh::Mesh( const NURBSExtension& ext )
+Mesh::Mesh(const NURBSExtension& ext)
    : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
 {
    SetEmpty();
    /// make an internal copy of the NURBSExtension
-   NURBSext = new NURBSExtension( ext );
+   NURBSext = new NURBSExtension(ext);
 
    Dim              = NURBSext->Dimension();
    NumOfVertices    = NURBSext->GetNV();
@@ -6537,6 +6538,8 @@ void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_ukv)
       Array<int> ukv_to_rpkv;
       GetEdgeToUniqueKnotvector(edge_to_ukv, ukv_to_rpkv);
    }
+
+   CorrectPatchTopoOrientations(edge_to_ukv);
 }
 
 void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
@@ -6547,9 +6550,9 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
    const int NPKV = NP * dim;    // number of patch knotvectors
    constexpr int notset = -9999999;
    // Sign convention
-   auto sign = [](int i) { return -1 - i; };
-   auto unsign = [](int i) { return (i < 0) ? -1 - i : i; };
-   // Edge index -> dimension convention
+   auto flipSign = [](int i) { return -1 - i; };
+   auto unSign = [](int i) { return (i < 0) ? -1 - i : i; };
+   // Local edge index -> dimension convention
    auto edge_to_dim = [](int i) { return (i < 8) ? ((i & 1) ? 1 : 0) : 2; };
 
    Array<int> v(2); // vertices of an edge
@@ -6564,7 +6567,7 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
       {
          GetElementVertices(i, v);
          // Sign is based on the edge's vertex indices
-         edge_to_ukv[i] = (v[1] > v[0]) ? i : sign(i);
+         edge_to_ukv[i] = (v[1] > v[0]) ? i : flipSign(i);
          ukv_to_rpkv[i] = i;
       }
       return;
@@ -6614,14 +6617,14 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
          // We've set this edge already - link this index to it
          if (edge_to_pkv[edge] != notset)
          {
-            const int pkv_other = unsign(edge_to_pkv[edge]);
+            const int pkv_other = unSign(edge_to_pkv[edge]);
             unite(pkv, pkv_other);
          }
          else
          {
             GetEdgeVertices(edge, v);
             // Sign is based on the edge's vertex indices
-            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : sign(pkv);
+            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : flipSign(pkv);
          }
       }
    }
@@ -6648,11 +6651,255 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
    edge_to_ukv.SetSize(NumOfEdges);
    for (int i = 0; i < NumOfEdges; i++)
    {
-      const int pkv = unsign(edge_to_pkv[i]);
+      const int pkv = unSign(edge_to_pkv[i]);
       const int rpkv = pkv_to_rpkv[pkv];
       const int ukv = rpkv_to_ukv[rpkv];
-      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? sign(ukv) : ukv;
+      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? flipSign(ukv) : ukv;
    }
+
+   CorrectPatchTopoOrientations(edge_to_ukv);
+}
+
+void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
+{
+   const int dim = Dimension(); // Topological (not physical) dimension
+   if (dim == 1) { return; }
+
+   // Sign convention
+   auto flipSign = [](int i) { return -1 - i; };
+
+   const Table *face2elem = GetFaceToElementTable();
+   Array<int> pfaces, orient;
+   Array<int> fe, feo;
+
+   // Finds elements sharing a face containing knotvector kv.
+   auto faceNeighbors = [&](int p, int kv, std::unordered_set<int> &nghb)
+   {
+      if (dim == 2) { GetElementEdges(p, pfaces, orient); }
+      else { GetElementFaces(p, pfaces, orient); }
+
+      for (auto face : pfaces)
+      {
+         // Check whether this face contains kv.
+         GetFaceEdges(face, fe, feo);
+         bool hasKV = false;
+         for (auto e : fe)
+         {
+            const int skv = edge_to_ukv[e];
+            if (skv == kv || flipSign(skv) == kv) { hasKV = true; }
+         }
+         if (hasKV)
+         {
+            Array<int> row;
+            face2elem->GetRow(face, row);
+            for (auto elem : row) { nghb.insert(elem); }
+         }
+      }
+   };
+
+   std::vector<std::vector<int>> dir_edges;
+   if (dim == 2)
+   {
+      dir_edges =
+      {
+         {0,2},
+         {1,3}
+      };
+   }
+   else
+   {
+      dir_edges =
+      {
+         {0,2,4,6},
+         {1,3,5,7},
+         {8,9,10,11}
+      };
+   }
+
+   Array<int> ukvs((dim==2) ? 4 : 12);
+   Array<int> pe, oe;
+   bool initKV = false;
+
+   auto setPatchDirections = [&](int p, int kv, Array<bool> &edgeSet,
+                                 std::unordered_set<int> &visited)
+   {
+      // Edges and orientations for this patch
+      GetElementEdges(p, pe, oe);
+
+      // Get the signed unique knot vector indices
+      for (int i = 0; i < pe.Size(); i++)
+      {
+         ukvs[i] = edge_to_ukv[pe[i]];
+         ukvs[i] = (oe[i] < 0) ? flipSign(ukvs[i]) : ukvs[i];
+      }
+
+      // Find the direction with this kv.
+      int thisDir = -1;
+      for (int d=0; d<dim; ++d) // Loop over directions.
+      {
+         const int skv = edge_to_ukv[pe[dir_edges[d][0]]];
+         if (skv == kv || flipSign(skv) == kv)
+         {
+            thisDir = d;
+         }
+      }
+      MFEM_VERIFY(thisDir >= 0, "");
+
+      // For this direction, find any edge already set. If no edge is set, we
+      // arbitrarily take the first.
+      int ref_edge0 = dir_edges[thisDir][0];
+      for (auto ref_edge : dir_edges[thisDir])
+      {
+         const int edge = pe[ref_edge];
+         if (edgeSet[edge])
+         {
+            ref_edge0 = ref_edge;
+         }
+      }
+
+      if (initKV && !edgeSet[pe[ref_edge0]])
+      {
+         visited.erase(p);
+         return false; // There is no set edge in this direction on this patch.
+      }
+
+      initKV = true;
+
+      // Use ref_edge0 to set other edges in this direction.
+      edgeSet[pe[ref_edge0]] = true;
+      for (auto i : dir_edges[thisDir])
+      {
+         if (i == ref_edge0)
+         {
+            continue;
+         }
+
+         const int edge = pe[i];
+         if ((dim == 2 && ukvs[i] != flipSign(ukvs[ref_edge0])) ||
+             (dim == 3 && ukvs[i] == flipSign(ukvs[ref_edge0])))
+         {
+            // Flip the sign of this edge
+            MFEM_VERIFY(!edgeSet[edge], "");
+            edge_to_ukv[edge] = flipSign(edge_to_ukv[edge]);
+         }
+
+         edgeSet[edge] = true;
+      }
+
+      return true;
+   };
+
+   Array<bool> edgeSet(NumOfEdges); // Whether edge has orientation set
+   edgeSet = false;
+
+   std::unordered_set<int> unset; // Patches with an unset edge
+   for (int i=0; i<NumOfElements; ++i) { unset.insert(i); }
+
+   const int max_iter = 3 * NumOfElements;
+   for (int iter=0; iter<max_iter; ++iter)
+   {
+      // Iteratively choose an unset patch (meaning not all edges have
+      // orientation set), choose a knotvector index for which the corresponding
+      // edges on this patch are not set, and sweep over all patches containing
+      // this knotvector. The patch sweep is ordered, by maintaining an ordered
+      // list `nextPatches` set by finding face-neighbor patches of visited
+      // patches, where the common face contains the knotvector. When each patch
+      // is visited, the edge orientations are set consistently. This iteration
+      // terminates when all edges have been set on all patches.
+
+      std::list<int> nextPatches; // Next patches to visit, ordered
+      std::unordered_set<int> nextSet; // nextPatches as a set
+      std::unordered_set<int> visited; // Visit each patch only once
+
+      if (unset.size() == 0)
+      {
+         break;
+      }
+
+      const int p0 = *unset.begin();
+      nextPatches.push_back(p0); // Start from arbitrary unset patch
+      nextSet.insert(p0);
+
+      // Choose an arbitrary unset direction for the first patch.
+      GetElementEdges(p0, pe, oe);
+      int unsetDim = -1;
+      for (int d=0; d<dim; ++d) // Loop over dimensions.
+      {
+         if (!edgeSet[pe[dir_edges[d][0]]])
+         {
+            unsetDim = d;
+         }
+      }
+
+      if (unsetDim == -1)
+      {
+         unset.erase(p0);
+         continue;
+      }
+
+      const int kv_signed = edge_to_ukv[pe[dir_edges[unsetDim][0]]];
+      const int kv = kv_signed < 0 ? flipSign(kv_signed) : kv_signed;
+      MFEM_VERIFY(!edgeSet[pe[dir_edges[unsetDim][0]]], "");
+
+      initKV = false;
+
+      while (nextPatches.size() > 0)
+      {
+         const int p = nextPatches.front();
+         nextPatches.pop_front();
+         nextSet.erase(p);
+         visited.insert(p);
+
+         const bool somethingSet = setPatchDirections(p, kv, edgeSet, visited);
+         if (!somethingSet)
+         {
+            continue;
+         }
+
+         // Find neighbors of patch p sharing a conforming face, via face2elem.
+         std::unordered_set<int> neighbors;
+         faceNeighbors(p, kv, neighbors);
+
+         bool allSet = true;
+         GetElementEdges(p, pe, oe);
+         for (auto edge : pe)
+         {
+            if (!edgeSet[edge])
+            {
+               allSet = false;
+            }
+         }
+         if (allSet)
+         {
+            unset.erase(p);
+         }
+
+         // Add neighbors not done to nextPatches.
+         for (auto n : neighbors)
+         {
+            if (n != p && visited.count(n) == 0 && unset.count(n) > 0)
+            {
+               if (nextSet.count(n) == 0)
+               {
+                  nextPatches.push_back(n);
+                  nextSet.insert(n);
+               }
+            }
+         }
+      }
+   }
+
+   bool allSet = true;
+   for (auto eset : edgeSet)
+   {
+      if (!eset)
+      {
+         allSet = false;
+      }
+   }
+   MFEM_VERIFY(allSet && unset.size() == 0, "Some edge is not set");
+
+   delete face2elem;
 }
 
 void Mesh::LoadNonconformingPatchTopo(std::istream &input,
@@ -9563,6 +9810,8 @@ void Mesh::GetVertices(Vector &vert_coord) const
 
 void Mesh::SetVertices(const Vector &vert_coord)
 {
+   MFEM_VERIFY(vert_coord.Size() == spaceDim * NumOfVertices, "");
+   vertices.SetSize(NumOfVertices);
    for (int i = 0, nv = vertices.Size(); i < nv; i++)
       for (int j = 0; j < spaceDim; j++)
       {

mesh/mesh.hpp

@@ -527,6 +527,9 @@ protected:
    void PrintTopoEdges(std::ostream &out, const Array<int> &e_to_k,
                        bool vmap = false) const;
 
+   /// Set signs to ensure knotvectors are pointed in the same direction.
+   void CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const;
+
    /// Used in GetFaceElementTransformations (...)
    void GetLocalPtToSegTransformation(IsoparametricTransformation &,
                                       int i) const;
@@ -984,8 +987,8 @@ public:
 
    ///@}
 
-   /// Construct a Mesh from a NURBSExtension
-   explicit Mesh( const NURBSExtension& ext );
+   /// Construct a Mesh from a NURBSExtension, which is deep-copied.
+   explicit Mesh(const NURBSExtension& ext);
 
    /** @anchor mfem_Mesh_construction
        @name Methods for piecewise Mesh construction.
@@ -2538,13 +2541,16 @@ public:
        changing the mesh file itself. Examples in miniapps/nurbs/meshes. */
    void RefineNURBSFromFile(std::string ref_file);
 
-   /// For NURBS meshes, insert the new knots in @a kv, for each direction.
+   /// For NURBS meshes, insert the new knots in @a kv, for each KnotVector.
+   /// The size of @a kv should be the number of KnotVectors in NURBSExtension.
    void KnotInsert(Array<KnotVector*> &kv);
 
-   /// For NURBS meshes, insert the knots in @a kv, for each direction.
+   /// For NURBS meshes, insert the knots in @a kv, for each KnotVector.
+   /// The size of @a kv should be the number of KnotVectors in NURBSExtension.
    void KnotInsert(Array<Vector*> &kv);
 
-   /// For NURBS meshes, remove the knots in @a kv, for each direction.
+   /// For NURBS meshes, remove the knots in @a kv, for each KnotVector.
+   /// The size of @a kv should be the number of KnotVectors in NURBSExtension.
    void KnotRemove(Array<Vector*> &kv);
 
    /* For each knot vector:

mesh/ncnurbs.cpp

@@ -2827,6 +2827,8 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
       }
    }
 
+   delete face2elem;
+
    // For any unset entries of kvf, set to default refinement factor rf_default.
    for (size_t i=0; i<kvf.size(); ++i)
    {

mesh/nurbs.cpp

@@ -2278,8 +2278,9 @@ void NURBSExtension::Load(std::istream &input, bool spacing)
 {
    own_topo = true;
 
-   CheckPatches();
-   // CheckBdrPatches();
+   MFEM_VERIFY(CheckPatches(),
+               "NURBSExtension::CheckPatch"
+               "\n  Inconsistent edge-to-knotvector mapping!");
 
    skip_comment_lines(input, '#');
 
@@ -2676,7 +2677,9 @@ NURBSExtension::NURBSExtension(const Mesh *patch_topology,
    patchTopo->GetEdgeToUniqueKnotvector(edge_to_ukv, ukv_to_rpkv);
    own_topo = true;
 
-   CheckPatches(); // This is checking the edge_to_ukv mapping
+   MFEM_VERIFY(CheckPatches(),
+               "NURBSExtension::CheckPatch"
+               "\n  Inconsistent edge-to-knotvector mapping!");
 
    // Set number of unique (not comprehensive) knot vectors
    NumOfKnotVectors = ukv_to_rpkv.Size();
@@ -3246,9 +3249,9 @@ void NURBSExtension::MergeGridFunctions(
    }
 }
 
-void NURBSExtension::CheckPatches()
+bool NURBSExtension::CheckPatches()
 {
-   if (Dimension() == 1 ) { return; }
+   if (Dimension() == 1 ) { return true; }
 
    Array<int> edges, oedge;
 
@@ -3256,6 +3259,7 @@ void NURBSExtension::CheckPatches()
    {
       patchTopo->GetElementEdges(p, edges, oedge);
 
+      // Convert to ukv and apply sign-flip
       for (int i = 0; i < edges.Size(); i++)
       {
          edges[i] = edge_to_ukv[edges[i]];
@@ -3265,6 +3269,16 @@ void NURBSExtension::CheckPatches()
          }
       }
 
+      // In 2d - opposite edges must be same knotvector with opposite sign.
+      // In 3d - opposite edges must be same knotvector with same sign.
+      // This logic is the result of Mesh::GetElementEdges setting orientation
+      // for edges based on ascending vertex indices, using reference vertex
+      // ordering
+      // {0, 1}, {1, 2}, {2, 3}, {3, 0} for Geometry::SQUARE in 2D
+      // and
+      // {0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6},
+      // {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7} for Geometry::CUBE in 3D
+      // See fem/geom.cpp for these definitions.
       if ((Dimension() == 2 &&
            (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
 
@@ -3275,11 +3289,10 @@ void NURBSExtension::CheckPatches()
             edges[8] != edges[9] || edges[8] != edges[10] ||
             edges[8] != edges[11])))
       {
-         mfem::err << "NURBSExtension::CheckPatch (patch = " << p
-                   << ")\n  Inconsistent edge-to-knotvector mapping!";
-         mfem_error();
+         return false;
       }
    }
+   return true;
 }
 
 void NURBSExtension::CheckBdrPatches()
@@ -3757,8 +3770,8 @@ void NURBSExtension::GetPatchOffsets(int &meshCounter, int &spaceCounter)
 
       if (dim == 1)
       {
-         meshCounter  += KnotVec(0)->GetNE() - 1;
-         spaceCounter += KnotVec(0)->GetNCP() - 2;
+         meshCounter  += KnotVec(p)->GetNE() - 1;
+         spaceCounter += KnotVec(p)->GetNCP() - 2;
       }
       else if (dim == 2)
       {

mesh/nurbs.hpp

@@ -596,11 +596,8 @@ protected:
        if the KnotVector index associated with edge @a edge is negative. */
    inline const KnotVector *KnotVec(int edge, int oedge, int *okv) const;
 
-   /// Throw an error if any patch has an inconsistent edge_to_ukv mapping.
-   void CheckPatches();
-
    /// Throw an error if any boundary patch has invalid KnotVector orientation.
-   void CheckBdrPatches();
+   MFEM_DEPRECATED void CheckBdrPatches();
 
    /** @brief Return the directions in @a kvdir of the KnotVectors in patch @a p
        based on the patch edge orientations. Each entry of @a kvdir is -1 if the
@@ -794,6 +791,9 @@ public:
    void MergeGridFunctions(GridFunction *gf_array[], int num_pieces,
                            GridFunction &merged);
 
+   /// Returns false if any patch has an inconsistent edge_to_ukv mapping.
+   bool CheckPatches();
+
    /// Destroy a NURBSExtension.
    virtual ~NURBSExtension();
 

miniapps/fluids/navier/navier_bifurcation.cpp

@@ -127,7 +127,7 @@ int main(int argc, char *argv[])
    }
 
    // Load mesh + complete any serial refinements
-   Mesh mesh("../../data/channel-bifurcation-2d.mesh");
+   Mesh mesh("../../../data/channel-bifurcation-2d.mesh");
    for (int lev = 0; lev < ctx.rs_levels; lev++)
    {
       mesh.UniformRefinement();

miniapps/gslib/CMakeLists.txt

@@ -15,16 +15,14 @@ set(MESH_GF_FILES
     triple-pt-1.gf
     triple-pt-2.gf
 )
-
+# add target which keeps required mesh files in sync
 set(SRC_MESH_GF_FILES)
 foreach(MESH_GF_FILE ${MESH_GF_FILES})
   list(APPEND SRC_MESH_GF_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_GF_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
+add_custom_target(copy_miniapps_gslib_data
   COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_GF_FILES} .
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying gslib miniapps data files ...")
-add_custom_target(copy_miniapps_gslib_data DEPENDS data_is_copied)
+  COMMENT "Syncing gslib miniapps data files ...")
 
 if (MFEM_USE_GSLIB)
   add_mfem_miniapp(schwarz_ex1

miniapps/meshing/CMakeLists.txt

@@ -27,11 +27,9 @@ set(SRC_MESH_FILES)
 foreach(MESH_FILE ${MESH_FILES})
   list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
+add_custom_target(copy_miniapps_meshing_data
   COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} .
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying meshing miniapps data files ...")
-add_custom_target(copy_miniapps_meshing_data DEPENDS data_is_copied)
+  COMMENT "Syncing meshing miniapps data files ...")
 
 add_mfem_miniapp(klein-bottle
   MAIN klein-bottle.cpp

miniapps/meshing/reflector.cpp

@@ -33,8 +33,7 @@
 using namespace std;
 using namespace mfem;
 
-
-void ReflectPoint(Vector & p, Vector const& origin, Vector const& normal)
+void ReflectPoint(Vector &p, const Vector &origin, const Vector &normal)
 {
    Vector diff(3);
    Vector proj(3);
@@ -57,12 +56,12 @@ private:
    // Map from reflected to original mesh elements
    std::vector<int> *r2o;
 
-   std::vector<std::vector<int>> *perm;
+   std::vector<std::array<int, 8>> *perm;
 
 public:
    ReflectedCoefficient(VectorCoefficient &A, Vector const& origin_,
                         Vector const& normal_, std::vector<int> *r2o_,
-                        Mesh *mesh, std::vector<std::vector<int>> *refPerm) :
+                        Mesh *mesh, std::vector<std::array<int, 8>> *refPerm) :
       VectorCoefficient(3), a(&A), origin(origin_), normal(normal_),
       meshOrig(mesh), r2o(r2o_), perm(refPerm)
    { }
@@ -109,7 +108,7 @@ void ReflectedCoefficient::Eval(Vector &V, ElementTransformation &T,
       // give the columns of A.
 
       // Permutation p is such that hex_reflected[i] = hex_init[p[i]]
-      const std::vector<int>& p = (*perm)[elem];
+      const std::array<int, 8>& p = (*perm)[elem];
 
       // ip is on reflected hex. We map from the reflected hex to the initial
       // hex, in reference space. Thus we use y = Ax + b, where x is in the
@@ -164,7 +163,7 @@ void ReflectedCoefficient::Eval(Vector &V, ElementTransformation &T,
 
 // Find perm such that h1[i] = h2[perm[i]]
 void GetHexPermutation(Array<int> const& h1, Array<int> const& h2,
-                       std::vector<int> & perm)
+                       std::array<int, 8> &perm)
 {
    std::map<int, int> h2inv;
    const int n = perm.size();
@@ -236,7 +235,7 @@ public:
    int AddElement(Array<int> const& vertices, const bool reorder);
 
    Mesh *mesh;
-   std::vector<std::vector<int>> refPerm;
+   std::vector<std::array<int, 8>> refPerm;
 
 private:
    std::vector<std::vector<int>> faces;
@@ -279,13 +278,13 @@ int HexMeshBuilder::AddElement(Array<int> const& vertices, const bool reorder)
       }
       while (reordered);
 
-      std::vector<int> perm_e(8);
+      std::array<int, 8> perm_e;
       GetHexPermutation(rvert, vertices, perm_e);
       refPerm.push_back(perm_e);
    }
    else
    {
-      refPerm.push_back(std::vector<int> {0, 1, 2, 3, 4, 5, 6, 7});
+      refPerm.push_back(std::array<int, 8> {0, 1, 2, 3, 4, 5, 6, 7});
    }
 
    SaveHexFaces(mesh->GetNE(), rvert);
@@ -762,7 +761,10 @@ bool GetMeshElementOrder(Mesh const& mesh, Vector const& origin,
    return true;
 }
 
-Mesh* ReflectHighOrderMesh(Mesh & mesh, Vector origin, Vector normal)
+Mesh* ReflectHighOrderMesh(Mesh &mesh,
+                           const Vector &origin, const Vector &normal,
+                           std::vector<std::array<int, 8>> &hexPerm,
+                           std::vector<int> &elOrder)
 {
    MFEM_VERIFY(mesh.Dimension() == 3, "Only 3D meshes can be reflected");
 
@@ -837,7 +839,6 @@ Mesh* ReflectHighOrderMesh(Mesh & mesh, Vector origin, Vector normal)
       }
    }
 
-   std::vector<int> elOrder;
    const bool onPlane = GetMeshElementOrder(mesh, origin, normal, elOrder);
 
    for (int eidx=0; eidx<mesh.GetNE(); eidx++)
@@ -1006,6 +1007,133 @@ Mesh* ReflectHighOrderMesh(Mesh & mesh, Vector origin, Vector normal)
       *reflected_nodes = newReflectedNodes;
    }
 
+   hexPerm = builder.refPerm;
+
+   return reflected;
+}
+
+void ReorderHexArray(const std::array<int, 3> &dim,
+                     const array<int, 8> &hexperm,
+                     std::array<int, 3> &dir, std::array<int, 3> &dims,
+                     Array3D<int> &permArray);
+
+NURBSPatch* ReflectPatch(NURBSPatch *patch, int nx, int ny, int nz,
+                         const Vector &origin, const Vector &normal,
+                         const std::array<int, 8> &hexPerm)
+{
+   // The hexahedral element for this patch in the reflected patch topology mesh
+   // is the reflection of an original patch topology mesh element, with
+   // reference vertices permuted according to hexPerm. The original grid of
+   // (nx + 1) x (ny + 1) x (nz + 1)
+   // control points has a new size and ordering, depending on hexPerm. Now,
+   // ReorderHexArray finds the new dimensions of this grid in `dims`, maps the
+   // directions in `dir`, and sets the permutation of grid indices as triples
+   // in `permArray`.
+   std::array<int, 3> dims, dir;
+   Array3D<int> permArray;
+   ReorderHexArray({nx+1, ny+1, nz+1}, hexPerm, dir, dims, permArray);
+
+   const KnotVector *kv0 = patch->GetKV(dir[0]);
+   const KnotVector *kv1 = patch->GetKV(dir[1]);
+   const KnotVector *kv2 = patch->GetKV(dir[2]);
+
+   NURBSPatch *rpatch = new NURBSPatch(kv0, kv1, kv2, 4);
+
+   // Reflect the control points in this reflected patch `rpatch`.
+   Vector vr(3);
+   for (int i=0; i<dims[0]; ++i)
+   {
+      for (int j=0; j<dims[1]; ++j)
+      {
+         for (int k=0; k<dims[2]; ++k)
+         {
+            const int old = permArray(i,j,k);
+            const int i0 = old / ((ny + 1) * (nz + 1));
+            const int j0 = (old - (i0 * (ny + 1) * (nz + 1))) / (nz + 1);
+            const int k0 = old - (i0 * (ny + 1) * (nz + 1)) - (j0 * (nz+1));
+
+            const real_t w = (*patch)(i0,j0,k0,3); // Weight
+            for (int l=0; l<3; ++l) { vr[l] = (*patch)(i0,j0,k0,l) / w; }
+
+            ReflectPoint(vr, origin, normal);
+
+            for (int l=0; l<3; ++l) { (*rpatch)(i,j,k,l) = vr[l] * w; }
+            (*rpatch)(i,j,k,3) = w;
+         }
+      }
+   }
+
+   return rpatch;
+}
+
+Mesh* ReflectNURBSMesh(Mesh &mesh, const Vector &origin, const Vector &normal)
+{
+   MFEM_VERIFY(mesh.NURBSext && mesh.Dimension() == 3,
+               "Only 3D NURBS meshes can be reflected");
+
+   Mesh patchTopo = mesh.NURBSext->GetPatchTopology(); // Deep copy
+
+   Array<NURBSPatch*> patchesOriginal, patches;
+   mesh.GetNURBSPatches(patchesOriginal); // Deep copy
+
+   NURBSPatchMap p2g(mesh.NURBSext);
+   const KnotVector *kv[3];
+
+   const int pnv = patchTopo.GetNV();
+   Vector vert_coord(3 * patchTopo.GetNV());
+   for (int p=0; p<patchesOriginal.Size(); ++p)
+   {
+      p2g.SetPatchDofMap(p, kv);
+      const int nx = p2g.nx();
+      const int ny = p2g.ny();
+      const int nz = p2g.nz();
+
+      Array<int> vert;
+      patchTopo.GetElementVertices(p, vert);
+
+      for (int l=0; l<3; ++l)
+      {
+         const int os = l * pnv;
+         vert_coord[vert[0] + os] = (*patchesOriginal[p])(0,0,0,l);
+         vert_coord[vert[1] + os] = (*patchesOriginal[p])(nx,0,0,l);
+         vert_coord[vert[2] + os] = (*patchesOriginal[p])(nx,ny,0,l);
+         vert_coord[vert[3] + os] = (*patchesOriginal[p])(0,ny,0,l);
+         vert_coord[vert[4] + os] = (*patchesOriginal[p])(0,0,nz,l);
+         vert_coord[vert[5] + os] = (*patchesOriginal[p])(nx,0,nz,l);
+         vert_coord[vert[6] + os] = (*patchesOriginal[p])(nx,ny,nz,l);
+         vert_coord[vert[7] + os] = (*patchesOriginal[p])(0,ny,nz,l);
+      }
+   }
+
+   patchTopo.SetVertices(vert_coord);
+
+   std::vector<std::array<int, 8>> hexPerm;
+   std::vector<int> elOrder;
+   Mesh *reflectedPatchTopo = ReflectHighOrderMesh(patchTopo, origin, normal,
+                                                   hexPerm, elOrder);
+
+   // Construct reflected patches. Note that reflectedPatchTopo has patch
+   // ordering depending on patchTopo.
+   for (int p=0; p<patchesOriginal.Size(); ++p)
+   {
+      const int p_orig = elOrder[p]; // TODO: use r2o instead?
+      p2g.SetPatchDofMap(p_orig, kv);
+      const int nx = p2g.nx();
+      const int ny = p2g.ny();
+      const int nz = p2g.nz();
+
+      patches.Append(patchesOriginal[p_orig]);
+      patches.Append(ReflectPatch(patchesOriginal[p_orig], nx, ny, nz,
+                                  origin, normal, hexPerm[(2 * p) + 1]));
+   }
+
+   NURBSExtension *ne = new NURBSExtension(reflectedPatchTopo, patches);
+   delete reflectedPatchTopo;
+
+   for (auto patch : patches) { delete patch; }
+
+   Mesh *reflected = new Mesh(*ne);
+   delete ne;
    return reflected;
 }
 
@@ -1045,7 +1173,19 @@ int main(int argc, char *argv[])
 
    Mesh mesh(mesh_file, 0, 0);
 
-   Mesh *reflected = ReflectHighOrderMesh(mesh, origin, normal);
+   Mesh *reflected{nullptr};
+
+   //if (mesh.IsNURBS()) // TODO: available in PR 4936
+   if (mesh.NURBSext)
+   {
+      reflected = ReflectNURBSMesh(mesh, origin, normal);
+   }
+   else
+   {
+      std::vector<std::array<int, 8>> hexPerm;
+      std::vector<int> elOrder;
+      reflected = ReflectHighOrderMesh(mesh, origin, normal, hexPerm, elOrder);
+   }
 
    // Save the final mesh
    ofstream mesh_ofs("reflected.mesh");
@@ -1065,3 +1205,85 @@ int main(int argc, char *argv[])
 
    return 0;
 }
+
+void HexVertexIJK(const int idx, std::array<int, 3>& ijk)
+{
+   ijk[2] = idx / 4;
+   const int id2d = idx - (4 * ijk[2]);
+   ijk[1] = id2d / 2;
+   ijk[0] = (ijk[1] == 0) ? id2d : 3 - id2d;
+}
+
+void ReorderHexArray(const std::array<int, 3> &dim,
+                     const array<int, 8> &hexperm,
+                     std::array<int, 3> &dir, std::array<int, 3> &dims,
+                     Array3D<int> &permArray)
+{
+   int prinV[4] = {0, 1, 3, 4}; // Vertices in principal directions (after 0)
+   int newPrinV[4];
+
+   // newVertices[i] = oldVertices[hexperm[i]]
+   // Hence newPrinV[0] = hexperm[0] is the index
+   // of new vertex 0 in the old hex.
+
+   std::array<int, 3> newIJK[4];
+   for (int i = 0; i < 4; ++i)
+   {
+      newPrinV[i] = hexperm[prinV[i]];
+      HexVertexIJK(newPrinV[i], newIJK[i]);
+   }
+
+   // For direction i in the new hex, dir[i] is the direction in the old hex.
+   Array<bool> rev(3);
+   for (int i = 0; i < 3; ++i)
+   {
+      bool iset = false;
+      for (int j = 0; j < 3; ++j)
+      {
+         const int d = newIJK[i + 1][j] - newIJK[0][j];
+         if (d != 0)
+         {
+            MFEM_VERIFY(!iset, "");
+            MFEM_VERIFY(d == 1 || d == -1, "");
+            dir[i] = j;
+            rev[i] = (d == -1);
+            iset = true;
+         }
+      }
+
+      MFEM_VERIFY(iset, "");
+
+      dims[i] = dim[dir[i]];
+   }
+
+   MFEM_VERIFY(dir[0] + dir[1] + dir[2] == 3, "");
+
+   permArray.SetSize(dims[0], dims[1], dims[2]);
+
+   Array<int> old_ijk(3);
+   Array<int> new_ijk(3);
+   for (int i = 0; i < dims[0]; ++i)
+      for (int j = 0; j < dims[1]; ++j)
+         for (int k = 0; k < dims[2]; ++k)
+         {
+            new_ijk[0] = i;
+            new_ijk[1] = j;
+            new_ijk[2] = k;
+
+            for (int m = 0; m < 3; ++m)
+            {
+               const int d = dir[m]; // Old hex direction
+               if (rev[m])
+               {
+                  old_ijk[d] = dim[d] - 1 - new_ijk[m];
+               }
+               else
+               {
+                  old_ijk[d] = new_ijk[m];
+               }
+            }
+
+            permArray(i, j, k) =
+               old_ijk[2] + (old_ijk[1] * dim[2]) + (old_ijk[0] * dim[1] * dim[2]);
+         }
+}

miniapps/multidomain/CMakeLists.txt

@@ -19,11 +19,9 @@ set(SRC_MESH_FILES)
 foreach(MESH_FILE ${MESH_FILES})
   list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
+add_custom_target(copy_miniapps_multidomain_data
   COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} .
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying multidomain miniapps data files ...")
-add_custom_target(copy_miniapps_multidomain_data DEPENDS data_is_copied)
+  COMMENT "Syncing multidomain miniapps data files ...")
 
 # Parallel apps.
 if (MFEM_USE_MPI)

miniapps/nurbs/CMakeLists.txt

@@ -9,31 +9,12 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-set(MESH_FILES
-  cube-nurbs.mesh
-  ijk-hex-nurbs.mesh
-  plus-nurbs-permuted.mesh
-  plus-nurbs.mesh
-  square-nurbs.mesh
-  two-cubes-nurbs-autoedge.mesh
-  two-cubes-nurbs-rot.mesh
-  two-cubes-nurbs.mesh
-  two-squares-nurbs-autoedge.mesh
-  two-squares-nurbs-rot.mesh
-  two-squares-nurbs.mesh
-)
-# Add a target to copy the mesh files from the source directory; used by sample
-# runs.
-set(SRC_MESH_FILES)
-foreach(MESH_FILE ${MESH_FILES})
-  list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/meshes/${MESH_FILE})
-endforeach()
-add_custom_command(OUTPUT data_is_copied
+# add target which keeps required mesh files in sync
+file(GLOB SRC_MESH_FILES CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/meshes/*)
+add_custom_target(copy_miniapps_nurbs_data
   COMMAND ${CMAKE_COMMAND} -E make_directory meshes
-  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} meshes/
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying nurbs miniapps data files ...")
-add_custom_target(copy_miniapps_nurbs_data DEPENDS data_is_copied)
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} meshes
+  COMMENT "Syncing nurbs miniapps data directory ...")
 
 add_mfem_miniapp(nurbs_ex1
   MAIN nurbs_ex1.cpp

miniapps/nurbs/meshes/3patch-nurbs-flipedge.mesh

@@ -0,0 +1,76 @@
+MFEM NURBS mesh v1.0
+
+# Same as 3patch-nurbs.mesh but with some flipped edges
+# This will fail to load without CorrectPatchTopoOrientations
+
+dimension
+2
+
+elements
+3
+1 3 0 1 4 3
+1 3 1 2 5 4
+1 3 5 6 3 4
+
+boundary
+0
+
+edges
+9
+2 1 0
+1 1 4
+2 4 3
+1 0 3
+0 1 2
+1 5 2
+0 5 4
+2 5 6
+0 3 6
+
+
+vertices
+7
+
+patches
+
+knotvectors
+2
+1 2 0 0 1 1
+1 2 0 0 1 1
+
+dimension
+2
+
+controlpoints
+0 0 1
+1 0 1
+0 1 1
+1 1 1
+
+knotvectors
+2
+1 2 0 0 1 1
+1 2 0 0 1 1
+
+dimension
+2
+
+controlpoints
+1 0 1
+2 0 1
+1 1 1
+2 2 1
+
+knotvectors
+2
+1 2 0 0 1 1
+1 2 0 0 1 1
+
+dimension
+2
+
+controlpoints
+2 2 1
+1 2 1
+1 1 1
+0 1 1

miniapps/nurbs/meshes/3patch-nurbs.mesh

@@ -0,0 +1,72 @@
+MFEM NURBS mesh v1.0
+
+dimension
+2
+
+elements
+3
+1 3 0 1 4 3
+1 3 1 2 5 4
+1 3 5 6 3 4
+
+boundary
+0
+
+edges
+9
+2 0 1
+1 1 4
+2 3 4
+1 0 3
+0 1 2
+1 2 5
+0 4 5
+2 6 5
+0 3 6
+
+vertices
+7
+
+patches
+
+knotvectors
+2
+1 2 0 0 1 1
+1 2 0 0 1 1
+
+dimension
+2
+
+controlpoints
+0 0 1
+1 0 1
+0 1 1
+1 1 1
+
+knotvectors
+2
+1 2 0 0 1 1
+1 2 0 0 1 1
+
+dimension
+2
+
+controlpoints
+1 0 1
+2 0 1
+1 1 1
+2 2 1
+
+knotvectors
+2
+1 2 0 0 1 1
+1 2 0 0 1 1
+
+dimension
+2
+
+controlpoints
+2 2 1
+1 2 1
+1 1 1
+0 1 1

miniapps/tools/CMakeLists.txt

@@ -20,11 +20,9 @@ set(SRC_MESH_GF_FILES)
 foreach(MESH_FILE ${MESH_GF_FILES})
   list(APPEND SRC_MESH_GF_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
 endforeach()
-add_custom_command(OUTPUT data_is_copied
-  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_GF_FILES} ../gslib/.
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying tools miniapps data files ...")
-add_custom_target(copy_miniapps_tools_data DEPENDS data_is_copied)
+add_custom_target(copy_miniapps_tools_data
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_GF_FILES} .
+  COMMENT "Syncing tools miniapps data files ...")
 
 add_mfem_miniapp(display-basis
   MAIN display-basis.cpp

miniapps/tribol/README.md

@@ -27,8 +27,8 @@ accelerations, so the relationship between forces/contact pressures and
 deformations/contact gaps is linear and, therefore, the problem can be solved
 exactly with a single linear solve. The mortar implementation is based on [Puso
 and Laursen (2004)](https://doi.org/10.1016/j.cma.2003.10.010). A description of
-the Tribol implementation is available in [Serac
-documentation](https://serac.readthedocs.io/en/latest/sphinx/theory_reference/solid.html#contact-mechanics).
+the Tribol implementation is available in [smith
+documentation](https://llnlsmith.readthedocs.io/en/latest/sphinx/theory_reference/solid.html#contact-mechanics).
 Lagrange multipliers are used to solve for the pressure required to prevent
 violation of the contact constraints.
 

tests/unit/CMakeLists.txt

@@ -11,6 +11,24 @@
 
 project(mfem-unit-tests NONE)
 
+# Define a target that all examples and miniapps will depend on.
+set(MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME test_exec_prerequisites)
+add_custom_target(${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
+
+# Add a target to copy the mfem data directory to the build directory
+# Implementable as a single copy_directory_if_different command w/ CMake >= 3.26
+file(GLOB DATA_FILES CONFIGURE_DEPENDS ${PROJECT_SOURCE_DIR}/data/*)
+add_custom_target(copy_test_data
+  COMMAND ${CMAKE_COMMAND} -E make_directory data
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${DATA_FILES} data
+  COMMENT "Syncing the test data directory ...")
+# Add 'copy_test_data' as a prerequisite for test executables, if the source and the
+# build directories are not the same.
+if (NOT ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}"))
+  add_dependencies(${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME} copy_test_data)
+  add_dependencies(${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME} copy_data)
+endif()
+
 # Include the source directory for the unit tests - catch.hpp is there.
 include_directories(BEFORE ${CMAKE_CURRENT_SOURCE_DIR})
 
@@ -53,6 +71,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_ode2.cpp
   linalg/test_operator.cpp
   linalg/test_particlevector.cpp
+  linalg/test_sparsesmoothers.cpp
   linalg/test_vector.cpp
   mesh/mesh_test_utils.cpp
   mesh/test_exodus_reader.cpp
@@ -68,8 +87,6 @@ set(UNIT_TESTS_SRCS
   mesh/test_psubmesh.cpp
   mesh/test_submesh.cpp
   mesh/test_vtu.cpp
-  mesh/test_nurbs.cpp
-  mesh/test_exodus_writer.cpp
   fem/make_permuted_mesh.cpp
   fem/test_1d_bilininteg.cpp
   fem/test_2d_bilininteg.cpp
@@ -168,25 +185,19 @@ endif()
 
 add_library(unit_tests_srcs OBJECT ${UNIT_TESTS_SRCS})
 target_link_libraries(unit_tests_srcs PUBLIC mfem)
+add_dependencies(unit_tests_srcs
+      ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
 
 # All serial non-device unit tests are built into a single executable,
 # 'unit_tests'.
 mfem_add_executable(unit_tests unit_test_main.cpp)
 target_link_libraries(unit_tests unit_tests_srcs)
 add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} unit_tests)
-# Unit tests need the ../../data directory.
-add_dependencies(unit_tests copy_data)
 # ParSubMesh tests need meshes in ../../miniapps/multidomain
 add_dependencies(unit_tests copy_miniapps_multidomain_data)
 # NURBS tests need meshes in ../../miniapps/nurbs
 add_dependencies(unit_tests copy_miniapps_nurbs_data)
 
-# Copy data to the build directory.
-add_custom_command(TARGET unit_tests POST_BUILD
-COMMAND ${CMAKE_COMMAND} -E copy_directory
-        ${CMAKE_CURRENT_SOURCE_DIR}/data data
-        COMMENT "Copying the unit tests data directory ...")
-
 # Create a test called 'unit_tests' that runs the 'unit_tests' executable.
 # The unit tests can be built and run separately from the rest of the tests:
 #   make unit_tests
@@ -205,7 +216,6 @@ if (MFEM_USE_CUDA)
     set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
     mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS})
     target_link_libraries(gpu_unit_tests unit_tests_srcs)
-    add_dependencies(gpu_unit_tests copy_data)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
     if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
         add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
@@ -221,7 +231,6 @@ if (MFEM_USE_HIP)
     set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
     mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS})
     target_link_libraries(gpu_unit_tests unit_tests_srcs)
-    add_dependencies(gpu_unit_tests copy_data)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
     if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
         add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
@@ -248,6 +257,8 @@ function(add_serial_miniapp_test name test_uvm)
     endif(MFEM_USE_HIP)
 
     mfem_add_executable(${name}_tests_cpu ${${NAME}_TESTS_SRCS})
+    add_dependencies(${name}_tests_cpu copy_miniapps_meshing_data)
+    add_dependencies(${name}_tests_cpu ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
     target_compile_definitions(${name}_tests_cpu PUBLIC MFEM_${NAME}_DEVICE="cpu")
     target_link_libraries(${name}_tests_cpu mfem)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_cpu)
@@ -256,6 +267,8 @@ function(add_serial_miniapp_test name test_uvm)
     endif()
 
     mfem_add_executable(${name}_tests_debug ${${NAME}_TESTS_SRCS})
+    add_dependencies(${name}_tests_debug copy_miniapps_meshing_data)
+    add_dependencies(${name}_tests_debug ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
     target_compile_definitions(${name}_tests_debug PUBLIC MFEM_${NAME}_DEVICE="debug")
     target_link_libraries(${name}_tests_debug mfem)
     add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_debug)
@@ -265,6 +278,8 @@ function(add_serial_miniapp_test name test_uvm)
 
     if (MFEM_USE_CUDA OR MFEM_USE_HIP)
         mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
+        add_dependencies(${name}_tests_gpu copy_miniapps_meshing_data)
+        add_dependencies(${name}_tests_gpu ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
         target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
         target_link_libraries(${name}_tests_gpu mfem)
         add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
@@ -274,6 +289,8 @@ function(add_serial_miniapp_test name test_uvm)
 
         if (test_uvm)
             mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
+            add_dependencies(${name}_tests_gpu_uvm copy_miniapps_meshing_data)
+            add_dependencies(${name}_tests_gpu_uvm ${MFEM_TEST_EXEC_PREREQUISITES_TARGET_NAME})
             target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
                MFEM_${NAME}_DEVICE="gpu:uvm")
             target_link_libraries(${name}_tests_gpu_uvm mfem)
@@ -289,7 +306,6 @@ endfunction(add_serial_miniapp_test)
 add_serial_miniapp_test(sedov ON)  # UVM ON
 add_serial_miniapp_test(tmop_pa OFF)  # UVM OFF
 # TMOP tests need meshes in ../../miniapps/meshing
-add_dependencies(tmop_pa_tests_cpu copy_miniapps_meshing_data)
 
 #-----------------------------------------------------------
 # SERIAL CEED TESTS:
@@ -343,7 +359,6 @@ if (MFEM_USE_MPI)
         set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
         set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
         mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS})
-        add_dependencies(pgpu_unit_tests copy_data)
         target_link_libraries(pgpu_unit_tests unit_tests_srcs)
         add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
         foreach(np 1 ${MFEM_MPI_NP})
@@ -359,7 +374,6 @@ if (MFEM_USE_MPI)
         # pgpu_unit_tests
         set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
         mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS})
-        add_dependencies(pgpu_unit_tests copy_data)
         target_link_libraries(pgpu_unit_tests unit_tests_srcs)
         add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
         foreach(np 1 ${MFEM_MPI_NP})

tests/unit/dfem/test_diffusion.cpp

@@ -69,9 +69,9 @@ template <int DIM> struct Diffusion
 };
 
 template <int DIM>
-void diffusion(const char *filename, int p)
+void diffusion(const char *filename, int p, bool include_mass)
 {
-   CAPTURE(filename, DIM, p);
+   CAPTURE(filename, DIM, p, include_mass);
 
    Mesh smesh(filename);
    ParMesh pmesh(MPI_COMM_WORLD, smesh);
@@ -116,6 +116,10 @@ void diffusion(const char *filename, int p)
 
    ParBilinearForm blf_fa(&pfes);
    blf_fa.AddDomainIntegrator(new DiffusionIntegrator(rho_coeff, ir));
+   if(include_mass)
+   {
+      blf_fa.AddDomainIntegrator(new MassIntegrator(ir));
+   }
    blf_fa.SetAssemblyLevel(AssemblyLevel::FULL);
    blf_fa.Assemble();
    blf_fa.Finalize();
@@ -130,6 +134,16 @@ void diffusion(const char *filename, int p)
    static constexpr int U = 0, Coords = 1, Rho = 3;
    const auto sol = std::vector{ FieldDescriptor{ U, &pfes } };
 
+   // Define mass qfunction in case it's used.
+      const auto mf_mass_qf =
+         [] MFEM_HOST_DEVICE(const dscalar_t u, const real_t &rho /*Needed for some reason*/,
+                             const tensor<real_t, DIM, DIM> &J, const real_t &w)
+      { return tuple{u * w * det(J)}; };
+      const auto desired_mass_qf =
+         [] MFEM_HOST_DEVICE(const dscalar_t u,
+                             const tensor<real_t, DIM, DIM> &J, const real_t &w)
+      { return tuple{u * w * det(J)}; };
+
    SECTION("action")
    {
       DOperator dop_mf(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
@@ -139,6 +153,21 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr);
+      if(include_mass)
+      {
+         // This version works, but had to include unused Rho parameter in qfunction.
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{},
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr);
+         // Below version does not work.
+         // dop_mf.AddDomainIntegrator(desired_mass_qf,
+         //                            tuple{ Value<U>{},
+         //                                   Gradient<Coords>{}, Weight{} },
+         //                            tuple{ Value<U>{} }, *ir,
+         //                            all_domain_attr);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
 
       pfes.GetRestrictionMatrix()->Mult(x, X);
@@ -159,6 +188,9 @@ void diffusion(const char *filename, int p)
 
    SECTION("action partial assembly")
    {
+      // is there a better way to skip a test than this?
+      if(!include_mass)
+      {
       static constexpr int QData = 2;
       UniformParameterSpace qd_ps(pmesh, *ir, DIM * DIM);
       ParameterFunction qdata(qd_ps);
@@ -196,6 +228,7 @@ void diffusion(const char *filename, int p)
 
       REQUIRE(norm_global == MFEM_Approx(0.0));
       MPI_Barrier(MPI_COMM_WORLD);
+      }
    }
 
    SECTION("action linearized")
@@ -208,6 +241,14 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr, derivatives);
+      if(include_mass)
+      {
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{} /*need this along with dummy argument to q function for some reason*/,
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr, derivatives);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
       auto dRdU = dop_mf.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
 
@@ -258,6 +299,10 @@ void diffusion(const char *filename, int p)
          ConstantCoefficient one(1.0);
          ParBilinearForm vblf_fa(&vpfes);
          vblf_fa.AddDomainIntegrator(new VectorDiffusionIntegrator(one, ir));
+         // if(include_mass)
+         // {
+         //    vblf_fa.AddDomainIntegrator(new VectorMassIntegrator(one, ir));
+         // }
          vblf_fa.SetAssemblyLevel(AssemblyLevel::LEGACYFULL);
          vblf_fa.Assemble();
          vblf_fa.Finalize();
@@ -290,6 +335,14 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr, derivatives);
+      if(include_mass)
+      {
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{},
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr, derivatives);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
       auto dRdU = dop_mf.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
 
@@ -306,6 +359,8 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
 
    const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
 
+   const bool include_mass = GENERATE(false, true);
+
    SECTION("2d")
    {
       const auto filename2d =
@@ -316,7 +371,7 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
             "../../data/inline-quad.mesh",
             "../../data/periodic-square.mesh"
          );
-      diffusion<2>(filename2d, p);
+      diffusion<2>(filename2d, p, include_mass);
    }
 
    SECTION("3d")
@@ -329,7 +384,7 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
             "../../data/toroid-hex.mesh",
             "../../data/periodic-cube.mesh"
          );
-      diffusion<3>(filename3d, p);
+      diffusion<3>(filename3d, p, include_mass);
    }
 }
 

tests/unit/linalg/test_sparsesmoothers.cpp

@@ -0,0 +1,76 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+static void TestTranspose(const Operator &A)
+{
+   DenseMatrix A_dense(A.Height(), A.Width());
+
+   Vector e(A.Width());
+   e = 0.0;
+   for (int i = 0; i < A.Width(); ++i)
+   {
+      e[i] = 1.0;
+
+      Vector Ae(A.Height());
+      A.Mult(e, Ae);
+      A_dense.SetCol(i, Ae);
+
+      e[i] = 0.0;
+   }
+
+   Vector v(A.Height());
+   v.Randomize();
+
+   Vector w1(A.Width()), w2(A.Width());
+
+   A.MultTranspose(v, w1);
+   A_dense.MultTranspose(v, w2);
+
+   w1 -= w2;
+
+   REQUIRE(w1.Normlinf() == MFEM_Approx(0.0));
+}
+
+TEST_CASE("Sparse Smoothers Transposed", "[DSmoother][GSSmoother]")
+{
+   const bool sym = GENERATE(true, false);
+
+   constexpr int n = 10;
+   SparseMatrix A(n, n);
+
+   for (int i = 0; i < n; ++i)
+   {
+      for (int j = 0; j < n; ++j)
+      {
+         const real_t val = rand_real();
+         A.Set(i, j, val);
+         if (sym) { A.Set(j, i, val); }
+      }
+      A.Add(i, i, 10.0);
+   }
+   A.Finalize();
+
+   constexpr int nit = 2; // Number of smoother iterations
+   TestTranspose(DSmoother(A, 0, 1.0, nit)); // scaled
+   if (sym)
+   {
+      TestTranspose(DSmoother(A, 1, 1.0, nit)); // l1-Jacobi
+      TestTranspose(DSmoother(A, 2, 1.0, nit)); // lumped Jacobi
+   }
+   TestTranspose(GSSmoother(A, 0, nit)); // symmetric
+   TestTranspose(GSSmoother(A, 1, nit)); // forward
+   TestTranspose(GSSmoother(A, 2, nit)); // backward
+}

tests/unit/mesh/test_nurbs.cpp

@@ -132,6 +132,14 @@ TEST_CASE("NURBS mesh reconstruction", "[NURBS]")
    for (auto *p : patches) { delete p; }
 }
 
+TEST_CASE("NURBS knotvector orientation", "[NURBS]")
+{
+   // This will fail to load without CorrectPatchTopoOrientations
+   auto mesh_fname = "../../miniapps/nurbs/meshes/3patch-nurbs-flipedge.mesh";
+   Mesh mesh(mesh_fname, 1, 1);
+   REQUIRE(mesh.NURBSext->CheckPatches());
+}
+
 TEST_CASE("NURBS NC-patch mesh loading", "[NURBS]")
 {
    auto mesh_fname = GENERATE("../../data/nc3-nurbs.mesh",