update removeinternalboundariesnotadjacentto

Quadrature function size fix

CHANGELOG

@@ -8,6 +8,24 @@
                                https://mfem.org
 
 
+Version 4.2.1 (development)
+===========================
+
+- Added interface to MUMPS direct solver. Its usage is demonstrated in ex25p.
+  See http://mumps.enseeiht.fr/ for more details. Supported versions >= 5.1.1.
+
+- Added three ESDIRK time integrators: implicit trapezoid rule, L-stable
+  ESDIRK-32, and A-stable ESDIRK-33.
+
+- Added a "scaled Jacobian" visualization option in the Mesh Explorer miniapp to
+  help identify elements with poor mesh quality.
+
+- Upgraded the Catch unit test framework from version 2.13.0 to version 2.13.2.
+
+- Implemented a filter method for the Navier miniapp to stabilize highly
+  turbulent flows in direct numerical simulation.
+
+
 Version 4.2, released on October 30, 2020
 =========================================
 

CMakeLists.txt

@@ -50,7 +50,7 @@ project(mfem NONE)
 # Current version of MFEM, see also `makefile`.
 #   mfem_VERSION = (string)
 #   MFEM_VERSION = (int)   [automatically derived from mfem_VERSION]
-set(${PROJECT_NAME}_VERSION 4.2.0)
+set(${PROJECT_NAME}_VERSION 4.2.1)
 
 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -185,7 +185,7 @@ if (MFEM_USE_MPI)
     endif()
   endif()
 else()
-  set(PKGS_NEED_MPI SUPERLU PETSC SLEPC STRUMPACK PUMI)
+  set(PKGS_NEED_MPI SUPERLU MUMPS PETSC SLEPC STRUMPACK PUMI)
   foreach(PKG IN LISTS PKGS_NEED_MPI)
     if (MFEM_USE_${PKG})
       message(STATUS "Disabling package ${PKG} - requires MPI")
@@ -264,6 +264,15 @@ if (MFEM_USE_SUPERLU)
   endif()
 endif()
 
+# MUMPS can only be enabled in parallel
+if (MFEM_USE_MUMPS)
+  if (MFEM_USE_MPI)
+    find_package(MUMPS REQUIRED mumps_common pord)
+  else()
+    message(FATAL_ERROR " *** MUMPS requires that MPI be enabled.")
+  endif()
+endif()
+
 # STRUMPACK can only be enabled in parallel
 if (MFEM_USE_STRUMPACK)
   if (MFEM_USE_MPI)
@@ -382,7 +391,7 @@ endif()
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
-    SLEPC MESQUITE SuperLUDist STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
+    SLEPC MESQUITE SuperLUDist MUMPS STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
     MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2
     CUSPARSE MKL_CPARDISO AMGX)
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.

INSTALL

@@ -348,6 +348,12 @@ MFEM_USE_SUPERLU5 = YES/NO
    If SuperLU functionality is enabled, use the older 5.1.0 version rather than
    the more recent 6+ versions.
 
+MFEM_USE_MUMPS = YES/NO
+   Enable MFEM functionality based on the MUMPS library. Currently, this
+   option adds the class MUMPSSolver (a parallel sparse direct solver). 
+   When enabled, this option uses the MUMPS_* library options, see
+   below.
+   
 MFEM_USE_STRUMPACK = YES/NO
    Enable MFEM functionality based on the STRUMPACK sparse direct solver and
    preconditioner through the STRUMPACKSolver and STRUMPACKRowLocMatrix
@@ -578,6 +584,12 @@ The specific libraries and their options are:
   Options: SUPERLU_OPT, SUPERLU_LIB.
   Versions: SuperLU_DIST >= 5.1.0.
 
+- MUMPS (optional), used when MFEM_USE_MUMPS = YES. Note that MUMPS
+  requires LAPACK, SCALAPACK and a reordering package such as PORD or METIS.
+  URL: http://mumps.enseeiht.fr
+  Options: MUMPS_OPT, MUMPS_LIB.
+  Versions: MUMPS >= 5.1.1
+
 - STRUMPACK (optional), used when MFEM_USE_STRUMPACK = YES. Note that STRUMPACK
   requires the PT-Scotch and Scalapack libraries as well as ParMETIS, which
   includes METIS 5 in its distribution. Starting with STRUMPACK v2.2.0, ParMETIS
@@ -837,6 +849,7 @@ MFEM_TIMER_TYPE - Set automatically, can be overwritten.
 MFEM_USE_MESQUITE
 MFEM_USE_SUITESPARSE
 MFEM_USE_SUPERLU
+MFEM_USE_MUMPS
 MFEM_USE_STRUMPACK
 MFEM_USE_GINKGO
 MFEM_USE_AMGX

config/cmake/MFEMConfig.cmake.in

@@ -32,6 +32,7 @@ set(MFEM_USE_SUNDIALS @MFEM_USE_SUNDIALS@)
 set(MFEM_USE_MESQUITE @MFEM_USE_MESQUITE@)
 set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)
+set(MFEM_USE_MUMPS @MFEM_USE_MUMPS@)
 set(MFEM_USE_STRUMPACK @MFEM_USE_STRUMPACK@)
 set(MFEM_USE_GINKGO @MFEM_USE_GINKGO@)
 set(MFEM_USE_AMGX @MFEM_USE_AMGX@)

config/cmake/config.hpp.in

@@ -86,6 +86,9 @@
 // Enable MFEM functionality based on the SuperLU_DIST library.
 #cmakedefine MFEM_USE_SUPERLU
 
+// Enable MFEM functionality based on the MUMPS library.
+#cmakedefine MFEM_USE_MUMPS
+
 // Enable MFEM functionality based on the STRUMPACK library.
 #cmakedefine MFEM_USE_STRUMPACK
 

config/cmake/modules/FindMUMPS.cmake

@@ -0,0 +1,23 @@
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Sets the following variables:
+#   - MUMPS_FOUND
+#   - MUMPS_INCLUDE_DIRS
+#   - MUMPS_LIBRARIES
+
+include(MfemCmakeUtilities)
+mfem_find_package(MUMPS MUMPS MUMPS_DIR
+  "include" dmumps_c.h "lib" dmumps
+  "Paths to headers required by MUMPS."
+  "Libraries required by MUMPS."
+  ADD_COMPONENT mumps_common "include" dmumps_c.h "lib" mumps_common
+  ADD_COMPONENT pord "include" dmumps_c.h "lib" pord)

config/config.hpp

@@ -34,6 +34,12 @@
 // Macro needed to get defines like M_PI from <cmath>. (Visual Studio C++ only?)
 #define _USE_MATH_DEFINES
 #endif
+// On Cygwin the option -std=c++11 prevents the definition of M_PI. Defining
+// the following macro allows us to get M_PI and some needed functions, e.g.
+// posix_memalign(), strdup(), strerror_r().
+#ifdef __CYGWIN__
+#define _XOPEN_SOURCE 600
+#endif
 
 // Check dependencies:
 
@@ -42,6 +48,9 @@
 #ifdef MFEM_USE_SUPERLU
 #error Building with SuperLU_DIST (MFEM_USE_SUPERLU=YES) requires MPI (MFEM_USE_MPI=YES)
 #endif
+#ifdef MFEM_USE_MUMPS
+#error Building with MUMPS (MFEM_USE_MUMPS=YES) requires MPI (MFEM_USE_MPI=YES)
+#endif
 #ifdef MFEM_USE_STRUMPACK
 #error Building with STRUMPACK (MFEM_USE_STRUMPACK=YES) requires MPI (MFEM_USE_MPI=YES)
 #endif

config/config.hpp.in

@@ -95,6 +95,10 @@
 // #define MFEM_USE_SUPERLU
 // #define MFEM_USE_SUPERLU5
 
+// Enable MFEM functionality based on the MUMPS library.
+// #define MFEM_USE_MUMPS
+// #define MFEM_MUMPS_VERSION @MFEM_MUMPS_VERSION@
+
 // Enable MFEM functionality based on the STRUMPACK library.
 // #define MFEM_USE_STRUMPACK
 

config/config.mk.in

@@ -33,6 +33,7 @@ MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
+MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@
 MFEM_USE_STRUMPACK     = @MFEM_USE_STRUMPACK@
 MFEM_USE_GINKGO        = @MFEM_USE_GINKGO@
 MFEM_USE_AMGX          = @MFEM_USE_AMGX@

config/defaults.cmake

@@ -34,6 +34,7 @@ option(MFEM_USE_MESQUITE "Enable MESQUITE usage" OFF)
 option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" OFF)
 option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
 option(MFEM_USE_SUPERLU5 "Use the old SuperLU_DIST 5.1 version" OFF)
+option(MFEM_USE_MUMPS "Enable MUMPS usage" OFF)
 option(MFEM_USE_STRUMPACK "Enable STRUMPACK usage" OFF)
 option(MFEM_USE_GINKGO "Enable Ginkgo usage" OFF)
 option(MFEM_USE_AMGX "Enable AmgX usage" OFF)
@@ -118,6 +119,15 @@ set(SuperLUDist_DIR "${MFEM_DIR}/../SuperLU_DIST_6.3.1" CACHE PATH
 set(SuperLUDist_REQUIRED_PACKAGES "MPI" "BLAS" "ParMETIS" CACHE STRING
     "Additional packages required by SuperLU_DIST.")
 
+set(MUMPS_DIR "${MFEM_DIR}/../MUMPS_5.2.0" CACHE PATH
+    "Path to the MUMPS library.")
+# Packages required by MUMPS, depending on how it was compiled.
+set(MUMPS_REQUIRED_PACKAGES "MPI" "BLAS" "METIS" "ScaLAPACK" CACHE STRING
+    "Additional packages required by MUMPS.")    
+# If the MPI package does not find all required Fortran libraries:
+# set(MUMPS_REQUIRED_LIBRARIES "gfortran" "mpi_mpifh" CACHE STRING
+#     "Additional libraries required by MUMPS.")    
+
 set(STRUMPACK_DIR "${MFEM_DIR}/../STRUMPACK-build" CACHE PATH
     "Path to the STRUMPACK library.")
 # STRUMPACK may also depend on "OpenMP", depending on how it was compiled.

config/defaults.mk

@@ -123,6 +123,7 @@ MFEM_USE_MESQUITE      = NO
 MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
+MFEM_USE_MUMPS         = NO
 MFEM_USE_STRUMPACK     = NO
 MFEM_USE_GINKGO        = NO
 MFEM_USE_AMGX          = NO
@@ -171,7 +172,7 @@ HYPRE_OPT = -I$(HYPRE_DIR)/include
 HYPRE_LIB = -L$(HYPRE_DIR)/lib -lHYPRE
 
 # METIS library configuration
-ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK),NONO)
+ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
    ifeq ($(MFEM_USE_METIS_5),NO)
      METIS_DIR = @MFEM_DIR@/../metis-4.0
      METIS_OPT =
@@ -252,12 +253,12 @@ SCOTCH_OPT = -I$(SCOTCH_DIR)/include
 SCOTCH_LIB = -L$(SCOTCH_DIR)/lib -lptscotch -lptscotcherr -lscotch -lscotcherr\
  -lpthread
 
-# SCALAPACK library configuration (required by STRUMPACK)
+# SCALAPACK library configuration (required by STRUMPACK and MUMPS)
 SCALAPACK_DIR = @MFEM_DIR@/../scalapack-2.0.2
 SCALAPACK_OPT = -I$(SCALAPACK_DIR)/SRC
 SCALAPACK_LIB = -L$(SCALAPACK_DIR)/lib -lscalapack $(LAPACK_LIB)
 
-# MPI Fortran library, needed e.g. by STRUMPACK
+# MPI Fortran library, needed e.g. by STRUMPACK or MUMPS
 # MPICH:
 MPI_FORTRAN_LIB = -lmpifort
 # OpenMPI:
@@ -265,6 +266,12 @@ MPI_FORTRAN_LIB = -lmpifort
 # Additional Fortan library:
 # MPI_FORTRAN_LIB += -lgfortran
 
+# MUMPS library configuration
+MUMPS_DIR = @MFEM_DIR@/../MUMPS_5.2.0
+MUMPS_OPT = -I$(MUMPS_DIR)/include
+MUMPS_LIB = -Wl,-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps\
+ -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
+
 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build
 STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)

config/get_mumps_version.cpp

@@ -0,0 +1,33 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "dmumps_c.h"
+#include <string>
+#include <iostream>
+#include <algorithm>
+
+// Macros to expand a macro as a string
+#define STR_EXPAND(s) #s
+#define STR(s) STR_EXPAND(s)
+
+int main()
+{
+#ifdef MUMPS_VERSION
+   const char *ptr = STR(MUMPS_VERSION);
+   std::string s(ptr);
+   s.erase(std::remove(s.begin(), s.end(), '"'), s.end());
+   s.erase(std::remove(s.begin(), s.end(), '.'), s.end());
+   std::cout << s << "\n";
+   return 0;
+#else
+   return -1;
+#endif
+}

config/makefile

@@ -42,6 +42,10 @@ GHV_FLAGS = $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(HYPRE_OPT))
 SMX = $(if $(MFEM_USE_PUMI:NO=),MFEM_USE_SIMMETRIX)
 SMX_PATH = $(PUMI_DIR)/include/gmi_sim.h
 SMX_FILE = $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(SMX_PATH))
+MUMPS = $(MFEM_USE_MUMPS:NO=)
+GMV_CXX ?= $(MFEM_CXX)
+GMV = get_mumps_version
+GMV_FLAGS = $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(MUMPS_OPT))
 
 $(GHV): $(SRC)$(GHV).cpp
 	$(call mfem-info, Determining HYPRE version ...)
@@ -50,6 +54,13 @@ $(GHV).out: $(GHV)
 	./$(GHV) > $(GHV).out
 .INTERMEDIATE: $(GHV) $(GHV).out
 
+$(GMV): $(SRC)$(GMV).cpp
+	$(call mfem-info, Determining MUMPS version ...)
+	$(GMV_CXX) ${GMV_FLAGS} $(SRC)$(GMV).cpp -o $(GMV)
+$(GMV).out: $(GMV)
+	./$(GMV) > $(GMV).out
+.INTERMEDIATE: $(GMV) $(GMV).out
+
 get-hypre-version: $(GHV).out
 	$(eval MFEM_HYPRE_VERSION:=$(shell cat $(GHV).out))
 	$(if $(MFEM_HYPRE_VERSION),$(eval export MFEM_HYPRE_VERSION)\
@@ -62,10 +73,18 @@ check-smx:
 	$(call mfem-info, MFEM_USE_SIMMETRIX = $(MFEM_USE_SIMMETRIX))
 	$(eval export MFEM_USE_SIMMETRIX)
 
-header: $(if $(MPI),get-hypre-version,) $(if $(SMX),check-smx)
+get-mumps-version: $(GMV).out
+	$(eval MFEM_MUMPS_VERSION:=$(shell cat $(GMV).out))
+	$(if $(MFEM_MUMPS_VERSION),$(eval export MFEM_MUMPS_VERSION)\
+	$(info MUMPS version: $(MFEM_MUMPS_VERSION)),\
+	$(error Unable to determine MUMPS version))
+
+header: $(if $(MPI),get-hypre-version,) $(if $(SMX),check-smx,) \
+ $(if $(MUMPS),get-mumps-version,)
 	$(call mfem-info, Writing $(CONFIG_HPP) ...)
 	@set -- && \
-	for def in $${MFEM_DEFINES} $(if $(MPI),MFEM_HYPRE_VERSION) $(SMX); do \
+	for def in $${MFEM_DEFINES} $(if $(MPI),MFEM_HYPRE_VERSION) $(SMX) \
+	           $(if $(MUMPS),MFEM_MUMPS_VERSION); do \
 	   eval var=\$$$$def && \
 	   if [ "NO" != "$${var}" ]; then \
 	      set -- "$$@" -e "s|// \(#define $${def} \)|\1|" && \

doc/CodeDocumentation.conf.in

@@ -38,7 +38,7 @@ PROJECT_NAME           = "MFEM"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = v4.2.0
+PROJECT_NUMBER         = v4.2.1
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

examples/ex25p.cpp

@@ -164,6 +164,8 @@ int main(int argc, char *argv[])
    int iprob = 4;
    double freq = 5.0;
    bool herm_conv = true;
+   bool slu_solver  = false;
+   bool mumps_solver = false;
    bool visualization = 1;
 
    OptionsParser args(argc, argv);
@@ -185,10 +187,26 @@ int main(int argc, char *argv[])
                   "Frequency (in Hz).");
    args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
                   "--no-hermitian", "Use convention for Hermitian operators.");
+#ifdef MFEM_USE_SUPERLU
+   args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
+                  "--no-superlu", "Use the SuperLU Solver.");
+#endif
+#ifdef MFEM_USE_MUMPS
+   args.AddOption(&mumps_solver, "-mumps", "--mumps-solver", "-no-mumps",
+                  "--no-mumps-solver", "Use the MUMPS Solver.");
+#endif
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
    args.Parse();
+   if (slu_solver && mumps_solver)
+   {
+      if (myid == 0)
+         cout << "WARNING: Both SuperLU and MUMPS have been selected,"
+              << " please choose either one." << endl
+              << "         Defaulting to SuperLU." << endl;
+      mumps_solver = false;
+   }
 
    if (iprob > 4) { iprob = 4; }
    prob = (prob_type)iprob;
@@ -430,6 +448,7 @@ int main(int argc, char *argv[])
 
    // 15. Solve using a direct or an iterative solver
 #ifdef MFEM_USE_SUPERLU
+   if (slu_solver)
    {
       // Transform to monolithic HypreParMatrix
       HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
@@ -442,7 +461,19 @@ int main(int argc, char *argv[])
       superlu.Mult(B, X);
       delete A;
    }
-#else
+#endif
+#ifdef MFEM_USE_MUMPS
+   if (mumps_solver)
+   {
+      HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
+      MUMPSSolver mumps;
+      mumps.SetPrintLevel(0);
+      mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
+      mumps.SetOperator(*A);
+      mumps.Mult(B,X);
+      delete A;
+   }
+#endif
    // 16a. Set up the parallel Bilinear form a(.,.) for the preconditioner
    //
    //    In Comp
@@ -450,6 +481,7 @@ int main(int argc, char *argv[])
    //
    //    In PML:   1/mu (abs(1/det(J) J^T J) Curl E, Curl F)
    //              + omega^2 * epsilon (abs(det(J) * (J^T J)^-1) * E, F)
+   if (!slu_solver && !mumps_solver)
    {
       ConstantCoefficient absomeg(pow(omega, 2) * epsilon);
       RestrictedCoefficient restr_absomeg(absomeg,attr);
@@ -499,7 +531,6 @@ int main(int argc, char *argv[])
       gmres.SetPreconditioner(BlockAMS);
       gmres.Mult(B, X);
    }
-#endif
 
    // 17. Recover the parallel grid function corresponding to X. This is the
    //     local finite element solution on each processor.

fem/bilininteg.hpp

@@ -2050,8 +2050,6 @@ public:
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe,
                                          ElementTransformation &Trans);
-
-   void SetupPA(const FiniteElementSpace &fes);
 };
 
 /** Mass integrator (u, v) restricted to the boundary of a domain */

fem/bilininteg_vecmass.cpp

@@ -106,7 +106,7 @@ void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
 static void PAVectorMassApply2D(const int NE,
-                                const Array<double> &_B,
+                                const Array<double> &B_,
                                 const Array<double> &_Bt,
                                 const Vector &_op,
                                 const Vector &_x,
@@ -119,7 +119,7 @@ static void PAVectorMassApply2D(const int NE,
    constexpr int VDIM = 2;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto Bt = Reshape(_Bt.Read(), D1D, Q1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, NE);
    auto x = Reshape(_x.Read(), D1D, D1D, VDIM, NE);
@@ -203,7 +203,7 @@ static void PAVectorMassApply2D(const int NE,
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
 static void PAVectorMassApply3D(const int NE,
-                                const Array<double> &_B,
+                                const Array<double> &B_,
                                 const Array<double> &_Bt,
                                 const Vector &_op,
                                 const Vector &_x,
@@ -216,7 +216,7 @@ static void PAVectorMassApply3D(const int NE,
    constexpr int VDIM = 3;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto Bt = Reshape(_Bt.Read(), D1D, Q1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
    auto x = Reshape(_x.Read(), D1D, D1D, D1D, VDIM, NE);
@@ -381,7 +381,7 @@ void VectorMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
 static void PAVectorMassAssembleDiagonal2D(const int NE,
-                                           const Array<double> &_B,
+                                           const Array<double> &B_,
                                            const Array<double> &_Bt,
                                            const Vector &_op,
                                            Vector &_diag,
@@ -393,7 +393,7 @@ static void PAVectorMassAssembleDiagonal2D(const int NE,
    constexpr int VDIM = 2;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, NE);
    auto y = Reshape(_diag.ReadWrite(), D1D, D1D, VDIM, NE);
    MFEM_FORALL(e, NE,
@@ -433,7 +433,7 @@ static void PAVectorMassAssembleDiagonal2D(const int NE,
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
 static void PAVectorMassAssembleDiagonal3D(const int NE,
-                                           const Array<double> &_B,
+                                           const Array<double> &B_,
                                            const Array<double> &_Bt,
                                            const Vector &_op,
                                            Vector &_diag,
@@ -445,7 +445,7 @@ static void PAVectorMassAssembleDiagonal3D(const int NE,
    constexpr int VDIM = 3;
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(_B.Read(), Q1D, D1D);
+   auto B = Reshape(B_.Read(), Q1D, D1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
    auto y = Reshape(_diag.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
    MFEM_FORALL(e, NE,

fem/coefficient.cpp

@@ -485,32 +485,32 @@ VectorSumCoefficient::VectorSumCoefficient(int dim)
 }
 
 VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &_A,
-                                           VectorCoefficient &_B,
+                                           VectorCoefficient &B_,
                                            double _alpha, double _beta)
    : VectorCoefficient(_A.GetVDim()),
-     ACoef(&_A), BCoef(&_B),
+     ACoef(&_A), BCoef(&B_),
      A(_A.GetVDim()), B(_A.GetVDim()),
      alphaCoef(NULL), betaCoef(NULL),
      alpha(_alpha), beta(_beta)
 {
-   MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
+   MFEM_ASSERT(_A.GetVDim() == B_.GetVDim(),
                "VectorSumCoefficient:  "
                "Arguments must have the same dimension.");
 }
 
 VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &_A,
-                                           VectorCoefficient &_B,
+                                           VectorCoefficient &B_,
                                            Coefficient &_alpha,
                                            Coefficient &_beta)
    : VectorCoefficient(_A.GetVDim()),
-     ACoef(&_A), BCoef(&_B),
+     ACoef(&_A), BCoef(&B_),
      A(_A.GetVDim()),
      B(_A.GetVDim()),
      alphaCoef(&_alpha),
      betaCoef(&_beta),
      alpha(0.0), beta(0.0)
 {
-   MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
+   MFEM_ASSERT(_A.GetVDim() == B_.GetVDim(),
                "VectorSumCoefficient:  "
                "Arguments must have the same dimension.");
 }
@@ -697,6 +697,11 @@ void OuterProductCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
    }
 }
 
+CrossCrossCoefficient::CrossCrossCoefficient(double A, VectorCoefficient &K)
+   : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(A), a(NULL), k(&K),
+     vk(K.GetVDim())
+{}
+
 CrossCrossCoefficient::CrossCrossCoefficient(Coefficient &A,
                                              VectorCoefficient &K)
    : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(0.0), a(&A), k(&K),

fem/coefficient.hpp

@@ -1171,8 +1171,8 @@ public:
                         double _alpha = 1.0, double _beta = 1.0);
 
    /** Constructor with scalar coefficients.
-       Result is _alpha * _A + _beta * _B */
-   VectorSumCoefficient(VectorCoefficient &_A, VectorCoefficient &_B,
+       Result is _alpha * _A + _beta * B_ */
+   VectorSumCoefficient(VectorCoefficient &_A, VectorCoefficient &B_,
                         Coefficient &_alpha, Coefficient &_beta);
 
    /// Reset the first vector coefficient
@@ -1201,7 +1201,7 @@ public:
    const Vector & GetA() const { return A; }
 
    /// Reset the second vector as a constant
-   void SetB(const Vector &_B) { B = _B; BCoef = NULL; }
+   void SetB(const Vector &B_) { B = B_; BCoef = NULL; }
    /// Return the second vector constant
    const Vector & GetB() const { return B; }
 

fem/conduitdatacollection.cpp

@@ -981,6 +981,11 @@ ConduitDataCollection::SaveRootFile(int num_domains,
    n_root["file_pattern"]     = MeshFilePattern(relay_protocol);
    n_root["tree_pattern"]     = "";
 
+   // Add the time, time step, and cycle
+   n_root["blueprint_index/mesh/state/time"] = time;
+   n_root["blueprint_index/mesh/state/time_step"] = time_step;
+   n_root["blueprint_index/mesh/state/cycle"] = cycle;
+
    relay::io::save(n_root, RootFileName(), root_proto);
 }
 

fem/convergence.cpp

@@ -110,7 +110,8 @@ void ConvergenceStudy::AddL2Error(GridFunction *gf,
 
 void ConvergenceStudy::AddGf(GridFunction *gf, Coefficient *scalar_u,
                              VectorCoefficient *grad,
-                             Coefficient *ell_coeff, double Nu)
+                             Coefficient *ell_coeff,
+                             JumpScaling jump_scaling)
 {
    cont_type = gf->FESpace()->FEColl()->GetContType();
 
@@ -140,7 +141,7 @@ void ConvergenceStudy::AddGf(GridFunction *gf, Coefficient *scalar_u,
 
    if (cont_type == mfem::FiniteElementCollection::DISCONTINUOUS && ell_coeff)
    {
-      double DGErr = gf->ComputeDGFaceJumpError(scalar_u,ell_coeff,Nu);
+      double DGErr = gf->ComputeDGFaceJumpError(scalar_u,ell_coeff,jump_scaling);
       DGFaceErrors.Append(DGErr);
       // Compute the rate of convergence by:
       // rate = log (||u - u_h|| / ||u - u_{h/2}||)/log(2)
@@ -270,26 +271,26 @@ void ConvergenceStudy::Print(bool relative, std::ostream &out)
             }
             out << "\n";
          }
-         if (cont_type == 3 && fcounter)
+      }
+      if (cont_type == 3 && fcounter)
+      {
+         out << " -------------------------------------------" << "\n";
+         out << "            DG Face Jump Error          " << "\n";
+         out << " -------------------------------------------"
+             << "\n";
+         out << std::right<< std::setw(11)<< "DOFs "<< std::setw(13);
+         out << "Error ";
+         out << std::setw(15) << "Rate " << "\n";
+         out << " -------------------------------------------"
+             << "\n";
+         out << std::setprecision(4);
+         for (int i =0; i<fcounter; i++)
          {
-            out << " -------------------------------------------" << "\n";
-            out << "            DG Face Jump Error          " << "\n";
-            out << " -------------------------------------------"
-                << "\n";
-            out << std::right<< std::setw(11)<< "DOFs "<< std::setw(13);
-            out << "Error ";
-            out << std::setw(15) << "Rate " << "\n";
-            out << " -------------------------------------------"
-                << "\n";
-            out << std::setprecision(4);
-            for (int i =0; i<fcounter; i++)
-            {
-               out << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
-                   << std::scientific << DGFaceErrors[i] << std::setw(13)
-                   << std::fixed << DGFaceRates[i] << "\n";
-            }
-            out << "\n";
+            out << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
+                << std::scientific << DGFaceErrors[i] << std::setw(13)
+                << std::fixed << DGFaceRates[i] << "\n";
          }
+         out << "\n";
       }
    }
 }

fem/convergence.hpp

@@ -59,7 +59,8 @@ private:
                    VectorCoefficient *vector_u);
    void AddGf(GridFunction *gf, Coefficient *scalar_u,
               VectorCoefficient *grad=nullptr,
-              Coefficient *ell_coeff=nullptr, double Nu=1.0);
+              Coefficient *ell_coeff=nullptr,
+              JumpScaling jump_scaling = {1.0, JumpScaling::ONE_OVER_H});
    void AddGf(GridFunction *gf, VectorCoefficient *vector_u,
               VectorCoefficient *curl, Coefficient *div);
    // returns the L2-norm of scalar_u or vector_u
@@ -75,9 +76,10 @@ public:
    /// DG face jumps parameters
    void AddL2GridFunction(GridFunction *gf, Coefficient *scalar_u,
                           VectorCoefficient *grad=nullptr,
-                          Coefficient *ell_coeff=nullptr, double Nu=1.0)
+                          Coefficient *ell_coeff=nullptr,
+                          JumpScaling jump_scaling = {1.0, JumpScaling::ONE_OVER_H})
    {
-      AddGf(gf, scalar_u, grad, ell_coeff, Nu);
+      AddGf(gf, scalar_u, grad, ell_coeff, jump_scaling);
    }
 
    /// Add H1 GridFunction, the exact solution and possibly its gradient

fem/gridfunc.cpp

@@ -2777,10 +2777,11 @@ double GridFunction::ComputeDivError(
 }
 
 double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
-                                            Coefficient *ell_coeff, double Nu,
+                                            Coefficient *ell_coeff,
+                                            class JumpScaling jump_scaling,
                                             const IntegrationRule *irs[])  const
 {
-   int fdof, dim, intorder, k;
+   int fdof, intorder, k;
    Mesh *mesh;
    const FiniteElement *fe;
    ElementTransformation *transf;
@@ -2791,20 +2792,24 @@ double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
    double error = 0.0;
 
    mesh = fes->GetMesh();
-   dim = mesh->Dimension();
 
    for (int i = 0; i < mesh->GetNumFaces(); i++)
    {
-      face_elem_transf = mesh->GetFaceElementTransformations(i, 5);
-      int i1 = face_elem_transf->Elem1No;
-      int i2 = face_elem_transf->Elem2No;
+      int i1, i2;
+      mesh->GetFaceElements(i, &i1, &i2);
+      double h = mesh->GetElementSize(i1);
       intorder = fes->GetFE(i1)->GetOrder();
       if (i2 >= 0)
+      {
          if ( (k = fes->GetFE(i2)->GetOrder()) > intorder )
          {
             intorder = k;
          }
+         h = std::min(h, mesh->GetElementSize(i2));
+      }
+      int p = intorder;
       intorder = 2 * intorder;  // <-------------
+      face_elem_transf = mesh->GetFaceElementTransformations(i, 5);
       const IntegrationRule *ir;
       if (irs)
       {
@@ -2875,8 +2880,9 @@ double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
       {
          const IntegrationPoint &ip = ir->IntPoint(j);
          transf->SetIntPoint(&ip);
-         error += (ip.weight * Nu * ell_coeff_val(j) *
-                   pow(transf->Weight(), 1.0-1.0/(dim-1)) *
+         double nu = jump_scaling.Eval(h, p);
+         error += (ip.weight * nu * ell_coeff_val(j) *
+                   transf->Weight() *
                    err_val(j) * err_val(j));
       }
    }
@@ -2884,6 +2890,15 @@ double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
    return (error < 0.0) ? -sqrt(-error) : sqrt(error);
 }
 
+double GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
+                                            Coefficient *ell_coeff,
+                                            double Nu,
+                                            const IntegrationRule *irs[])  const
+{
+   return ComputeDGFaceJumpError(
+             exsol, ell_coeff, {Nu, JumpScaling::ONE_OVER_H}, irs);
+}
+
 double GridFunction::ComputeH1Error(Coefficient *exsol,
                                     VectorCoefficient *exgrad,
                                     Coefficient *ell_coef, double Nu,
@@ -2892,7 +2907,11 @@ double GridFunction::ComputeH1Error(Coefficient *exsol,
    double error1 = 0.0;
    double error2 = 0.0;
    if (norm_type & 1) { error1 = GridFunction::ComputeGradError(exgrad); }
-   if (norm_type & 2) { error2 = GridFunction::ComputeDGFaceJumpError(exsol,ell_coef,Nu); }
+   if (norm_type & 2)
+   {
+      error2 = GridFunction::ComputeDGFaceJumpError(
+                  exsol, ell_coef, {Nu, JumpScaling::ONE_OVER_H});
+   }
 
    return sqrt(error1 * error1 + error2 * error2);
 }
@@ -3670,7 +3689,7 @@ QuadratureFunction & QuadratureFunction::operator=(double value)
 
 QuadratureFunction & QuadratureFunction::operator=(const Vector &v)
 {
-   MFEM_ASSERT(qspace && v.Size() == qspace->GetSize(), "");
+   MFEM_ASSERT(qspace && v.Size() == this->Size(), "");
    Vector::operator=(v);
    return *this;
 }

fem/gridfunc.hpp

@@ -451,7 +451,17 @@ public:
    virtual double ComputeDivError(Coefficient *exdiv,
                                   const IntegrationRule *irs[] = NULL) const;
 
-   /// Returns the Face Jumps error for L2 elements
+   /// Returns the Face Jumps error for L2 elements. The error can be weighted
+   /// by a constant nu, by nu/h, or nu*p^2/h, depending on the value of
+   /// @a jump_scaling.
+   virtual double ComputeDGFaceJumpError(Coefficient *exsol,
+                                         Coefficient *ell_coeff,
+                                         class JumpScaling jump_scaling,
+                                         const IntegrationRule *irs[] = NULL)
+   const;
+
+   /// Returns the Face Jumps error for L2 elements, with 1/h scaling.
+   MFEM_DEPRECATED
    virtual double ComputeDGFaceJumpError(Coefficient *exsol,
                                          Coefficient *ell_coeff,
                                          double Nu,
@@ -664,6 +674,32 @@ public:
     derived class ParGridFunction */
 std::ostream &operator<<(std::ostream &out, const GridFunction &sol);
 
+/// Class used to specify how the jump terms in
+/// GridFunction::ComputeDGFaceJumpError are scaled.
+class JumpScaling
+{
+public:
+   enum JumpScalingType
+   {
+      CONSTANT,
+      ONE_OVER_H,
+      P_SQUARED_OVER_H
+   };
+private:
+   double nu;
+   JumpScalingType type;
+public:
+   JumpScaling(double nu_=1.0, JumpScalingType type_=CONSTANT)
+      : nu(nu_), type(type_) { }
+   double Eval(double h, int p) const
+   {
+      double val = nu;
+      if (type != CONSTANT) { val /= h; }
+      if (type == P_SQUARED_OVER_H) { val *= p*p; }
+      return val;
+   }
+};
+
 
 /** @brief Class representing a function through its values (scalar or vector)
     at quadrature points. */
@@ -752,7 +788,8 @@ public:
 
    /// Copy the data from @a v.
    /** The size of @a v must be equal to the size of the associated
-       QuadratureSpace #qspace. */
+       QuadratureSpace #qspace times the QuadratureFunction dimension
+       i.e. QuadratureFunction::Size(). */
    QuadratureFunction &operator=(const Vector &v);
 
    /// Copy assignment. Only the data of the base class Vector is copied.

fem/lininteg.cpp

@@ -10,8 +10,8 @@
 // CONTRIBUTING.md for details.
 
 
-#include <cmath>
 #include "fem.hpp"
+#include <cmath>
 
 namespace mfem
 {

fem/pgridfunc.cpp

@@ -657,12 +657,12 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
 
 double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
                                                Coefficient *ell_coeff,
-                                               double Nu,
+                                               JumpScaling jump_scaling,
                                                const IntegrationRule *irs[]) const
 {
    const_cast<ParGridFunction *>(this)->ExchangeFaceNbrData();
 
-   int fdof, dim, intorder, k;
+   int fdof, intorder, k;
    ElementTransformation *transf;
    Vector shape, el_dofs, err_val, ell_coeff_val;
    Array<int> vdofs;
@@ -670,7 +670,6 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
    double error = 0.0;
 
    ParMesh *mesh = pfes->GetParMesh();
-   dim = mesh->Dimension();
 
    std::map<int,int> local_to_shared;
    for (int i = 0; i < mesh->GetNSharedFaces(); ++i)
@@ -687,6 +686,7 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
       mesh->GetFaceElements(i, &iel1, &iel2);
       mesh->GetFaceInfos(i, &info1, &info2);
 
+      double h = mesh->GetElementSize(iel1);
       intorder = fes->GetFE(iel1)->GetOrder();
 
       FaceElementTransformations *face_elem_transf;
@@ -703,11 +703,10 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
          }
          shared_face = true;
          shared_face_factor = 0.5;
+         h = std::min(h, mesh->GetFaceNbrElementSize(iel2));
       }
       else
       {
-         face_elem_transf = mesh->GetFaceElementTransformations(i);
-
          if (iel2 >= 0)
          {
             fe2 = pfes->GetFE(iel2);
@@ -715,12 +714,15 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
             {
                intorder = k;
             }
+            h = std::min(h, mesh->GetElementSize(iel2));
          }
          else
          {
             fe2 = NULL;
          }
+         face_elem_transf = mesh->GetFaceElementTransformations(i);
       }
+      int p = intorder;
 
       intorder = 2 * intorder;  // <-------------
       const IntegrationRule *ir;
@@ -806,8 +808,9 @@ double ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
       {
          const IntegrationPoint &ip = ir->IntPoint(j);
          transf->SetIntPoint(&ip);
-         error += shared_face_factor*(ip.weight * Nu * ell_coeff_val(j) *
-                                      pow(transf->Weight(), 1.0-1.0/(dim-1)) *
+         double nu = jump_scaling.Eval(h, p);
+         error += shared_face_factor*(ip.weight * nu * ell_coeff_val(j) *
+                                      transf->Weight() *
                                       err_val(j) * err_val(j));
       }
    }

fem/pgridfunc.hpp

@@ -310,7 +310,7 @@ public:
    /// Returns the Face Jumps error for L2 elements
    virtual double ComputeDGFaceJumpError(Coefficient *exsol,
                                          Coefficient *ell_coeff,
-                                         double Nu,
+                                         JumpScaling jump_scaling,
                                          const IntegrationRule *irs[]=NULL)
    const;
 

general/mem_manager.cpp

@@ -482,7 +482,8 @@ public:
       HostMemorySpace(),
       name(mm.GetUmpireAllocatorHostName()),
       rm(umpire::ResourceManager::getInstance()),
-      h_allocator(rm.isAllocator(name)? rm.getAllocator(name):
+      h_allocator((!std::strcmp(name, "HOST") || rm.isAllocator(name)) ?
+                  rm.getAllocator(name) :
                   rm.makeAllocator<umpire::strategy::DynamicPool>
                   (name, rm.getAllocator("HOST"))),
       strat(h_allocator.getAllocationStrategy()) { }
@@ -506,7 +507,8 @@ public:
       DeviceMemorySpace(),
       name(mm.GetUmpireAllocatorDeviceName()),
       rm(umpire::ResourceManager::getInstance()),
-      d_allocator(rm.isAllocator(name)? rm.getAllocator(name):
+      d_allocator((!std::strcmp(name, "DEVICE") || rm.isAllocator(name)) ?
+                  rm.getAllocator(name) :
                   rm.makeAllocator<umpire::strategy::DynamicPool>
                   (name, rm.getAllocator("DEVICE"))) { }
    void Alloc(Memory &base) { base.d_ptr = d_allocator.allocate(base.bytes); }

general/text.hpp

@@ -12,6 +12,7 @@
 #ifndef MFEM_TEXT
 #define MFEM_TEXT
 
+#include "../config/config.hpp"
 #include <istream>
 #include <iomanip>
 #include <sstream>
@@ -24,6 +25,8 @@ namespace mfem
 
 // Utilities for text parsing
 
+using std::to_string;
+
 /// Check if the stream starts with @a comment_char. If so skip it.
 inline void skip_comment_lines(std::istream &is, const char comment_char)
 {
@@ -47,18 +50,6 @@ inline void filter_dos(std::string &line)
    }
 }
 
-/// Convert an integer to an std::string.
-inline std::string to_string(int i)
-{
-   std::stringstream ss;
-   ss << i;
-
-   // trim leading spaces
-   std::string out_str = ss.str();
-   out_str = out_str.substr(out_str.find_first_not_of(" \t"));
-   return out_str;
-}
-
 /// Convert an integer to a 0-padded string with the given number of @a digits
 inline std::string to_padded_string(int i, int digits)
 {

linalg/CMakeLists.txt

@@ -83,6 +83,11 @@ if (MFEM_USE_SUPERLU)
   list(APPEND HDRS superlu.hpp)
 endif()
 
+if (MFEM_USE_MUMPS)
+  list(APPEND SRCS mumps.cpp)
+  list(APPEND HDRS mumps.hpp)
+endif()
+
 if (MFEM_USE_STRUMPACK)
   list(APPEND SRCS strumpack.cpp)
   list(APPEND HDRS strumpack.hpp)

linalg/kernels.hpp

@@ -12,11 +12,6 @@
 #ifndef MFEM_KERNELS_HPP
 #define MFEM_KERNELS_HPP
 
-#ifdef _WIN32
-#define _USE_MATH_DEFINES
-#include <cmath>
-#endif
-
 #include "../config/config.hpp"
 #include "../general/backends.hpp"
 #include "../general/globals.hpp"

linalg/linalg.hpp

@@ -49,6 +49,10 @@
 #include "hypre_parcsr.hpp"
 #include "hypre.hpp"
 
+#ifdef MFEM_USE_MUMPS
+#include "mumps.hpp"
+#endif
+
 #ifdef MFEM_USE_PETSC
 #include "petsc.hpp"
 #endif

linalg/mumps.cpp

@@ -0,0 +1,422 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MUMPS
+#ifdef MFEM_USE_MPI
+
+#include "mumps.hpp"
+
+#ifdef HYPRE_BIGINT
+#error "MUMPSSolver requires HYPRE_Int == int, for now."
+#endif
+
+// macro s.t. indices match MUMPS documentation
+#define MUMPS_ICNTL(I) icntl[(I) -1]
+#define MUMPS_INFO(I) info[(I) -1]
+
+namespace mfem
+{
+
+void MUMPSSolver::SetOperator(const Operator &op)
+{
+   auto APtr = dynamic_cast<const HypreParMatrix *>(&op);
+
+   MFEM_VERIFY(APtr, "Not compatible matrix type");
+
+   height = op.Height();
+   width = op.Width();
+
+   comm = APtr->GetComm();
+   MPI_Comm_size(comm, &numProcs);
+   MPI_Comm_rank(comm, &myid);
+
+   auto parcsr_op = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>(*APtr);
+
+   hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+#if MFEM_HYPRE_VERSION >= 21600
+   hypre_CSRMatrixBigJtoJ(csr_op);
+#endif
+
+   int *Iptr = csr_op->i;
+   int *Jptr = csr_op->j;
+   int n_loc = csr_op->num_rows;
+
+   row_start = parcsr_op->first_row_index;
+
+   MUMPS_INT8 nnz = 0;
+   if (mat_type)
+   {
+      // count nnz in case of symmetric mode
+      int k = 0;
+      for (int i = 0; i < n_loc; i++)
+      {
+         for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+         {
+            int ii = row_start + i + 1;
+            int jj = Jptr[k] + 1;
+            k++;
+            if (ii>=jj) { nnz++; }
+         }
+      }
+   }
+   else
+   {
+      nnz = csr_op->num_nonzeros;
+   }
+
+   int * I = new int[nnz];
+   int * J = new int[nnz];
+
+   // Fill in I and J arrays for
+   // COO format in 1-based indexing
+   int k = 0;
+   double * data;
+   if (mat_type)
+   {
+      int l = 0;
+      data = new double[nnz];
+      for (int i = 0; i < n_loc; i++)
+      {
+         for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+         {
+            int ii = row_start + i + 1;
+            int jj = Jptr[k] + 1;
+            if (ii >= jj)
+            {
+               I[l] = ii;
+               J[l] = jj;
+               data[l++] = csr_op->data[k];
+            }
+            k++;
+         }
+      }
+   }
+   else
+   {
+      for (int i = 0; i < n_loc; i++)
+      {
+         for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+         {
+            I[k] = row_start + i + 1;
+            J[k] = Jptr[k] + 1;
+            k++;
+         }
+      }
+      data = csr_op->data;
+   }
+
+   // new MUMPS object
+   if (id)
+   {
+      id->job = -2;
+      dmumps_c(id);
+      delete id;
+   }
+   id = new DMUMPS_STRUC_C;
+   // C to Fortran communicator
+   id->comm_fortran = (MUMPS_INT) MPI_Comm_c2f(comm);
+
+   // Host is involved in computation
+   id->par = 1;
+
+   id->sym = mat_type;
+
+   // MUMPS init
+   id->job = -1;
+   dmumps_c(id);
+
+   // Set MUMPS default parameters
+   SetParameters();
+
+   id->n = parcsr_op->global_num_rows;
+
+   id->nnz_loc = nnz;
+
+   id->irn_loc = I;
+
+   id->jcn_loc = J;
+
+   id->a_loc = data;
+
+   // MUMPS Analysis
+   id->job = 1;
+   dmumps_c(id);
+
+   // MUMPS Factorization
+   id->job = 2;
+   dmumps_c(id);
+
+   hypre_CSRMatrixDestroy(csr_op);
+   delete [] I;
+   delete [] J;
+   if (mat_type) { delete [] data; }
+
+#if MFEM_MUMPS_VERSION >= 530
+   delete [] irhs_loc;
+   irhs_loc = new int[n_loc];
+   for (int i = 0; i < n_loc; i++)
+   {
+      irhs_loc[i] = row_start + i + 1;
+   }
+   row_starts.SetSize(numProcs);
+   MPI_Allgather(&row_start, 1, MPI_INT, row_starts, 1, MPI_INT, comm);
+#else
+   if (myid == 0)
+   {
+      delete [] rhs_glob;
+      delete [] recv_counts;
+      rhs_glob = new double[parcsr_op->global_num_rows];
+      recv_counts = new int[numProcs];
+   }
+   MPI_Gather(&n_loc, 1, MPI_INT, recv_counts, 1, MPI_INT, 0, comm);
+   if (myid == 0)
+   {
+      delete [] displs;
+      displs = new int[numProcs];
+      displs[0] = 0;
+      int s = 0;
+      for (int k = 0; k < numProcs-1; k++)
+      {
+         s += recv_counts[k];
+         displs[k+1] = s;
+      }
+   }
+#endif
+}
+
+void MUMPSSolver::Mult(const Vector &x, Vector &y) const
+{
+#if MFEM_MUMPS_VERSION >= 530
+
+   id->nloc_rhs = x.Size();
+   id->lrhs_loc = x.Size();
+   id->rhs_loc = x.GetData();
+   id->irhs_loc = irhs_loc;
+
+   id->lsol_loc = id->MUMPS_INFO(23);
+   id->isol_loc = new int[id->MUMPS_INFO(23)];
+   id->sol_loc = new double[id->MUMPS_INFO(23)];
+
+   // MUMPS solve
+   id->job = 3;
+   dmumps_c(id);
+
+   RedistributeSol(id->isol_loc, id->sol_loc, y.GetData());
+
+   delete [] id->sol_loc;
+   delete [] id->isol_loc;
+#else
+   MPI_Gatherv(x.GetData(), x.Size(), MPI_DOUBLE,
+               rhs_glob, recv_counts,
+               displs, MPI_DOUBLE, 0, comm);
+
+   if (myid == 0) { id->rhs = rhs_glob; }
+
+   // MUMPS solve
+   id->job = 3;
+   dmumps_c(id);
+
+   MPI_Scatterv(rhs_glob, recv_counts, displs,
+                MPI_DOUBLE, y.GetData(), y.Size(),
+                MPI_DOUBLE, 0, comm);
+#endif
+}
+
+void MUMPSSolver::MultTranspose(const Vector &x, Vector &y) const
+{
+   // Set flag for Transpose Solve
+   id->MUMPS_ICNTL(9) = 0;
+   Mult(x,y);
+   // Reset the flag
+   id->MUMPS_ICNTL(9) = 1;
+
+}
+
+void MUMPSSolver::SetPrintLevel(int print_lvl)
+{
+   print_level = print_lvl;
+}
+
+void MUMPSSolver::SetMatrixSymType(MatType mtype)
+{
+   mat_type = mtype;
+}
+
+MUMPSSolver::~MUMPSSolver()
+{
+   if (id)
+   {
+#if MFEM_MUMPS_VERSION >= 530
+      delete [] irhs_loc;
+#else
+      delete [] recv_counts;
+      delete [] displs;
+      delete [] rhs_glob;
+#endif
+      id->job = -2;
+      dmumps_c(id);
+      delete id;
+   }
+}
+
+void MUMPSSolver::SetParameters()
+{
+   // output stream for error messages
+   id->MUMPS_ICNTL(1) = 6;
+   // output stream for diagnosting printing local to each proc
+   id->MUMPS_ICNTL(2) = 6;
+   // output stream for global info
+   id->MUMPS_ICNTL(3) = 6;
+   // Level of error printing
+   id->MUMPS_ICNTL(4) = print_level;
+   //input matrix format (assembled)
+   id->MUMPS_ICNTL(5) = 0;
+   // Use A or A^T
+   id->MUMPS_ICNTL(9) = 1;
+   // Iterative refinement (disabled)
+   id->MUMPS_ICNTL(10) = 0;
+   // Error analysis-statistics (disabled)
+   id->MUMPS_ICNTL(11) = 0;
+   // Use of ScaLAPACK (Parallel factorization on root)
+   id->MUMPS_ICNTL(13) = 0;
+   // Percentage increase of estimated workspace (default = 20%)
+   id->MUMPS_ICNTL(14) = 20;
+   // Number of OpenMP threads (default)
+   id->MUMPS_ICNTL(16) = 0;
+   // Matrix input format (distributed)
+   id->MUMPS_ICNTL(18) = 3;
+   // Schur complement (no Schur complement matrix returned)
+   id->MUMPS_ICNTL(19) = 0;
+
+#if MFEM_MUMPS_VERSION >= 530
+   // Distributed RHS
+   id->MUMPS_ICNTL(20) = 10;
+   // Distributed Sol
+   id->MUMPS_ICNTL(21) = 1;
+#else
+   // Centralized RHS
+   id->MUMPS_ICNTL(20) = 0;
+   // Centralized Sol
+   id->MUMPS_ICNTL(21) = 0;
+#endif
+   // Out of core factorization and solve (disabled)
+   id->MUMPS_ICNTL(22) = 0;
+   // Max size of working memory (default = based on estimates)
+   id->MUMPS_ICNTL(23) = 0;
+}
+
+#if MFEM_MUMPS_VERSION >= 530
+int MUMPSSolver::GetRowRank(int i, const Array<int> &row_starts_) const
+{
+   if (row_starts_.Size() == 1)
+   {
+      return 0;
+   }
+   auto up = std::upper_bound(row_starts_.begin(), row_starts_.end(), i);
+   return std::distance(row_starts_.begin(), up) - 1;
+}
+
+void MUMPSSolver::RedistributeSol(const int * row_map,
+                                  const double * x, double * y) const
+{
+   int size = id->MUMPS_INFO(23);
+   int * send_count = new int[numProcs]();
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1;
+      int row_rank = GetRowRank(j, row_starts);
+      if (myid == row_rank) { continue; }
+      send_count[row_rank]++;
+   }
+
+   int * recv_count = new int[numProcs];
+   MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
+
+   int * send_displ = new int [numProcs]; send_displ[0] = 0;
+   int * recv_displ = new int [numProcs]; recv_displ[0] = 0;
+   int sbuff_size = send_count[numProcs-1];
+   int rbuff_size = recv_count[numProcs-1];
+   for (int k = 0; k < numProcs - 1; k++)
+   {
+      send_displ[k + 1] = send_displ[k] + send_count[k];
+      recv_displ[k + 1] = recv_displ[k] + recv_count[k];
+      sbuff_size += send_count[k];
+      rbuff_size += recv_count[k];
+   }
+
+   int * sendbuf_index = new int[sbuff_size];
+   double * sendbuf_values = new double[sbuff_size];
+   int * soffs = new int[numProcs]();
+
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1;
+      int row_rank = GetRowRank(j, row_starts);
+      if (myid == row_rank)
+      {
+         int local_index = j - row_start;
+         y[local_index] = x[i];
+      }
+      else
+      {
+         int k = send_displ[row_rank] + soffs[row_rank];
+         sendbuf_index[k] = j;
+         sendbuf_values[k] = x[i];
+         soffs[row_rank]++;
+      }
+   }
+
+   int * recvbuf_index = new int[rbuff_size];
+   double * recvbuf_values = new double[rbuff_size];
+   MPI_Alltoallv(sendbuf_index,
+                 send_count,
+                 send_displ,
+                 MPI_INT,
+                 recvbuf_index,
+                 recv_count,
+                 recv_displ,
+                 MPI_INT,
+                 comm);
+   MPI_Alltoallv(sendbuf_values,
+                 send_count,
+                 send_displ,
+                 MPI_DOUBLE,
+                 recvbuf_values,
+                 recv_count,
+                 recv_displ,
+                 MPI_DOUBLE,
+                 comm);
+
+   // Unpack recv buffer
+   for (int i = 0; i < rbuff_size; i++)
+   {
+      int local_index = recvbuf_index[i] - row_start;
+      y[local_index] = recvbuf_values[i];
+   }
+
+   delete [] recvbuf_values;
+   delete [] recvbuf_index;
+   delete [] soffs;
+   delete [] sendbuf_values;
+   delete [] sendbuf_index;
+   delete [] recv_displ;
+   delete [] send_displ;
+   delete [] recv_count;
+   delete [] send_count;
+}
+#endif
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MUMPS

linalg/mumps.hpp

@@ -0,0 +1,156 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_MUMPS
+#define MFEM_MUMPS
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MUMPS
+#ifdef MFEM_USE_MPI
+#include "operator.hpp"
+#include "hypre.hpp"
+
+#include <mpi.h>
+#include "dmumps_c.h"
+#include <vector>
+
+namespace mfem
+{
+
+/**
+ * @brief MUMPS: A Parallel Sparse Direct Solver
+ *
+ * Interface for the distributed MUMPS solver
+ */
+class MUMPSSolver : public mfem::Solver
+{
+public:
+   enum MatType
+   {
+      UNSYMMETRIC = 0,
+      SYMMETRIC_INDEFINITE = 1,
+      SYMMETRIC_POSITIVE_DEFINITE = 2
+   };
+
+   /**
+    * @brief Default Constructor
+    */
+   MUMPSSolver() {}
+
+   /**
+    * @brief Set the Operator and perform factorization
+    *
+    * @a op needs to be of type HypreParMatrix.
+    *
+    * @param op Operator used in factorization and solve
+    */
+   void SetOperator(const Operator &op);
+
+   /**
+    * @brief Solve y = Op^{-1} x.
+    *
+    * @param x RHS vector
+    * @param y Solution vector
+    */
+   void Mult(const Vector &x, Vector &y) const;
+
+   /**
+    * @brief Transpose Solve y = Op^{-T} x.
+    *
+    * @param x RHS vector
+    * @param y Solution vector
+    */
+   void MultTranspose(const Vector &x, Vector &y) const;
+
+   /**
+    * @brief Set the error print level for MUMPS
+    *
+    * @param print_lvl Print level
+    *
+    * @note This method has to be called before SetOperator.
+    */
+   void SetPrintLevel(int print_lvl);
+
+   /**
+    * @brief Set the matrix type
+    *
+    * Supported matrix types: General, symmetric indefinite and
+    * symmetric positive definite
+    *
+    * @param mtype Matrix type
+    *
+    * @note This method has to be called before SetOperator.
+    */
+   void SetMatrixSymType(MatType mtype);
+
+   // Destructor
+   ~MUMPSSolver();
+
+private:
+
+   // MPI communicator
+   MPI_Comm comm;
+
+   // Number of procs
+   int numProcs;
+
+   // local mpi id
+   int myid;
+
+   // parameter controling the matrix type
+   MatType mat_type = MatType::UNSYMMETRIC;
+
+   // parameter controling the printing level
+   int print_level = 0;
+
+   // local row offsets
+   int row_start;
+
+   // MUMPS object
+   DMUMPS_STRUC_C *id=nullptr;
+
+   // Method for setting MUMPS interal parameters
+   void SetParameters();
+
+#if MFEM_MUMPS_VERSION >= 530
+
+   // row offests array on all procs
+   Array<int> row_starts;
+
+   // row map
+   int * irhs_loc = nullptr;
+
+   // These two methods are needed to distribute the local solution
+   // vectors returned by MUMPS to the original MFEM parallel partition
+   int GetRowRank(int i, const Array<int> &row_starts_) const;
+
+   void RedistributeSol(const int * row_map,
+                        const double * x,
+                        double * y) const;
+#else
+
+   // Arrays needed for MPI_Gather and MPI_Scatter
+   int * recv_counts = nullptr;
+
+   int * displs = nullptr;
+
+   double * rhs_glob = nullptr;
+
+#endif
+
+}; // mfem::MUMPSSolver class
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MUMPS
+#endif // MFEM_MUMPS

linalg/ode.cpp

@@ -677,7 +677,7 @@ void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
 
    f->SetTime(t + c*dt);
    f->ImplicitSolve(a*dt, y, k);
-   x.Add((1.-a-b)*dt, k);
+   x.Add((1.0-a-b)*dt, k);
 
    f->SetTime(t + dt);
    f->ImplicitSolve(a*dt, x, k);
@@ -685,6 +685,101 @@ void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
    t += dt;
 }
 
+void TrapezoidalRuleSolver::Init(TimeDependentOperator &_f)
+{
+   ODESolver::Init(_f);
+   k.SetSize(f->Width(), mem_type);
+   y.SetSize(f->Width(), mem_type);
+}
+
+void TrapezoidalRuleSolver::Step(Vector &x, double &t, double &dt)
+{
+   //   0   |   0    0
+   //   1   |  1/2  1/2
+   // ------+-----------
+   //       |  1/2  1/2
+   f->SetTime(t);
+   f->Mult(x,k);
+   add(x, dt/2.0, k, y);
+   x.Add(dt/2.0, k);
+
+   f->SetTime(t + dt);
+   f->ImplicitSolve(dt/2.0, y, k);
+   x.Add(dt/2.0, k);
+   t += dt;
+}
+
+void ESDIRK32Solver::Init(TimeDependentOperator &_f)
+{
+   ODESolver::Init(_f);
+   k.SetSize(f->Width(), mem_type);
+   y.SetSize(f->Width(), mem_type);
+   z.SetSize(f->Width(), mem_type);
+}
+
+void ESDIRK32Solver::Step(Vector &x, double &t, double &dt)
+{
+   //   0   |    0      0    0
+   //   2a  |    a      a    0
+   //   1   |  1-b-a    b    a
+   // ------+--------------------
+   //       |  1-b-a    b    a
+   const double a = (2.0 - sqrt(2.0)) / 2.0;
+   const double b = (1.0 - 2.0*a) / (4.0*a);
+
+   f->SetTime(t);
+   f->Mult(x,k);
+   add(x, a*dt, k, y);
+   add(x, (1.0-b-a)*dt, k, z);
+   x.Add((1.0-b-a)*dt, k);
+
+   f->SetTime(t + (2.0*a)*dt);
+   f->ImplicitSolve(a*dt, y, k);
+   z.Add(b*dt, k);
+   x.Add(b*dt, k);
+
+   f->SetTime(t + dt);
+   f->ImplicitSolve(a*dt, z, k);
+   x.Add(a*dt, k);
+   t += dt;
+}
+
+void ESDIRK33Solver::Init(TimeDependentOperator &_f)
+{
+   ODESolver::Init(_f);
+   k.SetSize(f->Width(), mem_type);
+   y.SetSize(f->Width(), mem_type);
+   z.SetSize(f->Width(), mem_type);
+}
+
+void ESDIRK33Solver::Step(Vector &x, double &t, double &dt)
+{
+   //   0   |      0          0        0
+   //   2a  |      a          a        0
+   //   1   |    1-b-a        b        a
+   // ------+----------------------------
+   //       |  1-b_2-b_3     b_2      b_3
+   const double a   = (3.0 + sqrt(3.0)) / 6.0;
+   const double b   = (1.0 - 2.0*a) / (4.0*a);
+   const double b_2 = 1.0 / ( 12.0*a*(1.0 - 2.0*a) );
+   const double b_3 = (1.0 - 3.0*a) / ( 3.0*(1.0 - 2.0*a) );
+
+   f->SetTime(t);
+   f->Mult(x,k);
+   add(x, a*dt, k, y);
+   add(x, (1.0-b-a)*dt, k, z);
+   x.Add((1.0-b_2-b_3)*dt, k);
+
+   f->SetTime(t + (2.0*a)*dt);
+   f->ImplicitSolve(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);
+   x.Add(b_3*dt, k);
+   t += dt;
+}
 
 void GeneralizedAlphaSolver::Init(TimeDependentOperator &_f)
 {

linalg/ode.hpp

@@ -460,6 +460,48 @@ public:
 };
 
 
+/** Two stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
+    of order 2. A-stable. */
+class TrapezoidalRuleSolver : public ODESolver
+{
+protected:
+   Vector k, y;
+
+public:
+   virtual void Init(TimeDependentOperator &_f);
+
+   virtual void Step(Vector &x, double &t, double &dt);
+};
+
+
+/** Three stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
+    of order 2. L-stable. */
+class ESDIRK32Solver : public ODESolver
+{
+protected:
+   Vector k, y, z;
+
+public:
+   virtual void Init(TimeDependentOperator &_f);
+
+   virtual void Step(Vector &x, double &t, double &dt);
+};
+
+
+/** Three stage, explicit singly diagonal implicit Runge-Kutta (ESDIRK) method
+    of order 3. A-stable. */
+class ESDIRK33Solver : public ODESolver
+{
+protected:
+   Vector k, y, z;
+
+public:
+   virtual void Init(TimeDependentOperator &_f);
+
+   virtual void Step(Vector &x, double &t, double &dt);
+};
+
+
 /// Generalized-alpha ODE solver from "A generalized-α method for integrating
 /// the filtered Navier-Stokes equations with a stabilized finite element
 /// method" by K.E. Jansen, C.H. Whiting and G.M. Hulbert.

linalg/vector.cpp

@@ -22,6 +22,10 @@
 #endif
 #endif
 
+#ifdef MFEM_USE_OPENMP
+#include <omp.h>
+#endif
+
 #include <iostream>
 #include <iomanip>
 #include <cmath>
@@ -1076,6 +1080,30 @@ double Vector::operator*(const Vector &v) const
 #ifdef MFEM_USE_OPENMP
    if (Device::Allows(Backend::OMP_MASK))
    {
+#define MFEM_USE_OPENMP_DETERMINISTIC_DOT
+#ifdef MFEM_USE_OPENMP_DETERMINISTIC_DOT
+      // By default, use a deterministic way of computing the dot product
+      static Vector th_dot;
+      #pragma omp parallel
+      {
+         const int nt = omp_get_num_threads();
+         #pragma omp master
+         th_dot.SetSize(nt);
+         const int tid    = omp_get_thread_num();
+         const int stride = (size + nt - 1)/nt;
+         const int start  = tid*stride;
+         const int stop   = std::min(start + stride, size);
+         double my_dot = 0.0;
+         for (int i = start; i < stop; i++)
+         {
+            my_dot += m_data[i] * v_data[i];
+         }
+         #pragma omp barrier
+         th_dot(tid) = my_dot;
+      }
+      return th_dot.Sum();
+#else
+      // The standard way of computing the dot product is non-deterministic
       double prod = 0.0;
       #pragma omp parallel for reduction(+:prod)
       for (int i = 0; i < size; i++)
@@ -1083,8 +1111,9 @@ double Vector::operator*(const Vector &v) const
          prod += m_data[i] * v_data[i];
       }
       return prod;
+#endif // MFEM_USE_OPENMP_DETERMINISTIC_DOT
    }
-#endif
+#endif // MFEM_USE_OPENMP
    if (Device::Allows(Backend::DEBUG_DEVICE))
    {
       const int N = size;

linalg/vector.hpp

@@ -42,10 +42,19 @@ namespace mfem
 inline int CheckFinite(const double *v, const int n);
 
 /// Define a shortcut for std::numeric_limits<double>::infinity()
+#ifndef __CYGWIN__
 inline double infinity()
 {
    return std::numeric_limits<double>::infinity();
 }
+#else
+// On Cygwin math.h defines a function 'infinity()' which will conflict with the
+// above definition if we have 'using namespace mfem;' and try to use something
+// like 'double a = infinity();'. This 'infinity()' function is non-standard and
+// is defined by the Newlib C standard library implementation used by Cygwin,
+// see https://en.wikipedia.org/wiki/Newlib, http://www.sourceware.org/newlib.
+using ::infinity;
+#endif
 
 /// Vector data type.
 class Vector

makefile

@@ -10,7 +10,7 @@
 # CONTRIBUTING.md for details.
 
 # The current MFEM version as an integer, see also `CMakeLists.txt`.
-MFEM_VERSION = 40200
+MFEM_VERSION = 40201
 MFEM_VERSION_STRING = $(shell printf "%06d" $(MFEM_VERSION) | \
   sed -e 's/^0*\(.*.\)\(..\)\(..\)$$/\1.\2.\3/' -e 's/\.0/./g' -e 's/\.0$$//')
 
@@ -204,7 +204,7 @@ CXXFLAGS ?= $(OPTIM_FLAGS)
 # MPI configuration
 ifneq ($(MFEM_USE_MPI),YES)
    MFEM_HOST_CXX = $(CXX)
-   PKGS_NEED_MPI = SUPERLU STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
+   PKGS_NEED_MPI = SUPERLU MUMPS STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
    $(foreach mpidep,$(PKGS_NEED_MPI),$(if $(MFEM_USE_$(mpidep):NO=),\
      $(warning *** [MPI is OFF] setting MFEM_USE_$(mpidep) = NO)\
      $(eval override MFEM_USE_$(mpidep)=NO),))
@@ -265,7 +265,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
  OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX
 
@@ -332,7 +332,8 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT MFEM_USE_PUMI\
  MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_OCCA\
  MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2\
- MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX MFEM_USE_MUMPS MFEM_SOURCE_DIR\
+ MFEM_INSTALL_DIR
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -366,7 +367,7 @@ MFEM_INSTALL_DIR = $(abspath $(MFEM_PREFIX))
 # If we have 'config' target, export variables used by config/makefile
 ifneq (,$(filter config,$(MAKECMDGOALS)))
    export $(MFEM_DEFINES) MFEM_DEFINES $(MFEM_CONFIG_VARS) MFEM_CONFIG_VARS
-   export VERBOSE HYPRE_OPT PUMI_DIR
+   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT
 endif
 
 # If we have 'install' target, export variables used by config/makefile
@@ -634,6 +635,7 @@ status info:
 	$(info MFEM_USE_MESQUITE      = $(MFEM_USE_MESQUITE))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
+	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))
 	$(info MFEM_USE_GINKGO        = $(MFEM_USE_GINKGO))
 	$(info MFEM_USE_AMGX          = $(MFEM_USE_AMGX))

mesh/mesh.cpp

@@ -29,6 +29,7 @@
 #include <cstring>
 #include <ctime>
 #include <functional>
+#include <unordered_set>
 
 // Include the METIS header, if using version 5. If using METIS 4, the needed
 // declarations are inlined below, i.e. no header is needed.
@@ -71,10 +72,14 @@ void Mesh::GetElementCenter(int i, Vector &center)
    eltransf->Transform(Geometries.GetCenter(geom), center);
 }
 
-double Mesh::GetElementSize(int i, int type)
+double Mesh::GetElementSize(ElementTransformation *T, int type)
 {
    DenseMatrix J(Dim);
-   GetElementJacobian(i, J);
+
+   Geometry::Type geom = T->GetGeometryType();
+   T->SetIntPoint(&Geometries.GetCenter(geom));
+   Geometries.JacToPerfJac(geom, T->Jacobian(), J);
+
    if (type == 0)
    {
       return pow(fabs(J.Det()), 1./Dim);
@@ -89,6 +94,11 @@ double Mesh::GetElementSize(int i, int type)
    }
 }
 
+double Mesh::GetElementSize(int i, int type)
+{
+   return GetElementSize(GetElementTransformation(i), type);
+}
+
 double Mesh::GetElementSize(int i, const Vector &dir)
 {
    DenseMatrix J(Dim);
@@ -1463,6 +1473,13 @@ void Mesh::AddBdrQuadAsTriangles(const int *vi, int attr)
    }
 }
 
+int Mesh::AddBdrPoint(int v, int attr)
+{
+   CheckEnlarge(boundary, NumOfBdrElements);
+   boundary[NumOfBdrElements] = new Point(&v, attr);
+   return NumOfBdrElements++;
+}
+
 void Mesh::GenerateBoundaryElements()
 {
    int i, j;
@@ -10533,6 +10550,227 @@ void Mesh::RemoveInternalBoundaries()
    attribs.Copy(bdr_attributes);
 }
 
+void Mesh::RemoveInternalBoundaries(const Array<int> &_keep)
+{
+   if (NURBSext || ncmesh) { return; }
+
+   std::unordered_set<int> keep(_keep.GetData(),
+                                _keep.GetData() + _keep.Size());
+
+   int num_bdr_elem = 0;
+   int new_bel_to_edge_nnz = 0;
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      if (FaceIsInterior(GetBdrElementEdgeIndex(i)) && 
+         (keep.count(boundary[i]->GetAttribute()) == 0))
+      {
+         FreeElement(boundary[i]);
+      }
+      else
+      {
+         num_bdr_elem++;
+         if (Dim == 3)
+         {
+            new_bel_to_edge_nnz += bel_to_edge->RowSize(i);
+         }
+      }
+   }
+
+   if (num_bdr_elem == GetNBE()) { return; }
+
+   Array<Element *> new_boundary(num_bdr_elem);
+   Array<int> new_be_to_edge, new_be_to_face;
+   Table *new_bel_to_edge = NULL;
+   new_boundary.SetSize(0);
+   if (Dim == 2)
+   {
+      new_be_to_edge.Reserve(num_bdr_elem);
+   }
+   else if (Dim == 3)
+   {
+      new_be_to_face.Reserve(num_bdr_elem);
+      new_bel_to_edge = new Table;
+      new_bel_to_edge->SetDims(num_bdr_elem, new_bel_to_edge_nnz);
+   }
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      /// if it's not interior or is in the keep boundary list
+      if (!FaceIsInterior(GetBdrElementEdgeIndex(i)) || 
+          keep.count(boundary[i]->GetAttribute()))
+      {
+         new_boundary.Append(boundary[i]);
+         if (Dim == 2)
+         {
+            new_be_to_edge.Append(be_to_edge[i]);
+         }
+         else if (Dim == 3)
+         {
+            int row = new_be_to_face.Size();
+            new_be_to_face.Append(be_to_face[i]);
+            int *e = bel_to_edge->GetRow(i);
+            int ne = bel_to_edge->RowSize(i);
+            int *new_e = new_bel_to_edge->GetRow(row);
+            for (int j = 0; j < ne; j++)
+            {
+               new_e[j] = e[j];
+            }
+            new_bel_to_edge->GetI()[row+1] = new_bel_to_edge->GetI()[row] + ne;
+         }
+      }
+   }
+
+   NumOfBdrElements = new_boundary.Size();
+   mfem::Swap(boundary, new_boundary);
+
+   if (Dim == 2)
+   {
+      mfem::Swap(be_to_edge, new_be_to_edge);
+   }
+   else if (Dim == 3)
+   {
+      mfem::Swap(be_to_face, new_be_to_face);
+      delete bel_to_edge;
+      bel_to_edge = new_bel_to_edge;
+   }
+
+   Array<int> attribs(num_bdr_elem);
+   for (int i = 0; i < attribs.Size(); i++)
+   {
+      attribs[i] = GetBdrAttribute(i);
+   }
+   attribs.Sort();
+   attribs.Unique();
+   bdr_attributes.DeleteAll();
+   attribs.Copy(bdr_attributes);
+}
+
+void Mesh::RemoveInternalBoundariesNotAdjacentTo(const Array<int> &_regions)
+{
+   if (NURBSext || ncmesh) { return; }
+
+   std::unordered_set<int> regions(_regions.GetData(),
+                                   _regions.GetData() + _regions.Size());
+
+   int num_bdr_elem = 0;
+   int new_bel_to_edge_nnz = 0;
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      auto faceNo = GetBdrElementEdgeIndex(i);
+      bool interior = FaceIsInterior(faceNo);
+
+      auto face_info = faces_info[faceNo];
+      auto adj = regions.count(elements[face_info.Elem1No]->GetAttribute());
+      if (face_info.Elem2No > 0)
+         adj += regions.count(elements[face_info.Elem2No]->GetAttribute());
+
+      bool remove = interior & !adj;
+
+      // if (remove)
+      // {
+      //    std::cout << "remove face on: " << boundary[i]->GetAttribute()
+      //              << " adj to ("
+      //              << elements[face_info.Elem1No]->GetAttribute();
+      //    if (face_info.Elem2No > 0)
+      //    {
+      //       std::cout << ", "
+      //                 << elements[face_info.Elem2No]->GetAttribute() << ")\n";
+      //    }
+      //    else
+      //    {
+      //       std::cout << ")\n";
+      //    }
+      // }
+
+      if (remove)
+      {
+         FreeElement(boundary[i]);
+      }
+      else
+      {
+         num_bdr_elem++;
+         if (Dim == 3)
+         {
+            new_bel_to_edge_nnz += bel_to_edge->RowSize(i);
+         }
+      }
+   }
+
+   if (num_bdr_elem == GetNBE()) { return; }
+
+   Array<Element *> new_boundary(num_bdr_elem);
+   Array<int> new_be_to_edge, new_be_to_face;
+   Table *new_bel_to_edge = NULL;
+   new_boundary.SetSize(0);
+   if (Dim == 2)
+   {
+      new_be_to_edge.Reserve(num_bdr_elem);
+   }
+   else if (Dim == 3)
+   {
+      new_be_to_face.Reserve(num_bdr_elem);
+      new_bel_to_edge = new Table;
+      new_bel_to_edge->SetDims(num_bdr_elem, new_bel_to_edge_nnz);
+   }
+   for (int i = 0; i < GetNBE(); i++)
+   {
+      auto faceNo = GetBdrElementEdgeIndex(i);
+      bool interior = FaceIsInterior(faceNo);
+
+      auto face_info = faces_info[faceNo];
+      auto adj = regions.count(elements[face_info.Elem1No]->GetAttribute());
+      if (face_info.Elem2No > 0)
+         adj += regions.count(elements[face_info.Elem2No]->GetAttribute());
+
+      bool keep = !interior || adj;
+
+      if (keep)
+      {
+         new_boundary.Append(boundary[i]);
+         if (Dim == 2)
+         {
+            new_be_to_edge.Append(be_to_edge[i]);
+         }
+         else if (Dim == 3)
+         {
+            int row = new_be_to_face.Size();
+            new_be_to_face.Append(be_to_face[i]);
+            int *e = bel_to_edge->GetRow(i);
+            int ne = bel_to_edge->RowSize(i);
+            int *new_e = new_bel_to_edge->GetRow(row);
+            for (int j = 0; j < ne; j++)
+            {
+               new_e[j] = e[j];
+            }
+            new_bel_to_edge->GetI()[row+1] = new_bel_to_edge->GetI()[row] + ne;
+         }
+      }
+   }
+
+   NumOfBdrElements = new_boundary.Size();
+   mfem::Swap(boundary, new_boundary);
+
+   if (Dim == 2)
+   {
+      mfem::Swap(be_to_edge, new_be_to_edge);
+   }
+   else if (Dim == 3)
+   {
+      mfem::Swap(be_to_face, new_be_to_face);
+      delete bel_to_edge;
+      bel_to_edge = new_bel_to_edge;
+   }
+
+   Array<int> attribs(num_bdr_elem);
+   for (int i = 0; i < attribs.Size(); i++)
+   {
+      attribs[i] = GetBdrAttribute(i);
+   }
+   attribs.Sort();
+   attribs.Unique();
+   bdr_attributes.DeleteAll();
+   attribs.Copy(bdr_attributes);
+}
+
 void Mesh::FreeElement(Element *E)
 {
 #ifdef MFEM_USE_MEMALLOC

mesh/mesh.hpp

@@ -471,6 +471,8 @@ protected:
    void GetElementData(const Array<Element*> &elem_array, int geom,
                        Array<int> &elem_vtx, Array<int> &attr) const;
 
+   double GetElementSize(ElementTransformation *T, int type = 0);
+
 public:
 
    Mesh() { SetEmpty(); }
@@ -555,6 +557,8 @@ public:
    int AddBdrQuad(const int *vi, int attr = 1);
    void AddBdrQuadAsTriangles(const int *vi, int attr = 1);
 
+   int AddBdrPoint(int v, int attr = 1);
+
    void GenerateBoundaryElements();
    /// Finalize the construction of a triangular Mesh.
    void FinalizeTriMesh(int generate_edges = 0, int refine = 0,
@@ -1302,6 +1306,12 @@ public:
        have two adjacent faces in 3D, or edges in 2D. */
    void RemoveInternalBoundaries();
 
+   /** Remove boundary elements that lie in the interior of the mesh, except 
+       for those with boundary attributes in `keep` */
+   void RemoveInternalBoundaries(const Array<int> &keep);
+
+   void RemoveInternalBoundariesNotAdjacentTo(const Array<int> &regions);
+
    /** @brief Get the size of the i-th element relative to the perfect
        reference element. */
    double GetElementSize(int i, int type = 0);

mesh/pmesh.cpp

@@ -1735,6 +1735,11 @@ void ParMesh::GetFaceNbrElementTransformation(
    }
 }
 
+double ParMesh::GetFaceNbrElementSize(int i, int type)
+{
+   return GetElementSize(GetFaceNbrElementTransformation(i), type);
+}
+
 void ParMesh::DeleteFaceNbrData()
 {
    if (!have_face_nbr_data)

mesh/pmesh.hpp

@@ -305,6 +305,10 @@ public:
       return &FaceNbrTransformation;
    }
 
+   /// Get the size of the i-th face neighbor element relative to the reference
+   /// element.
+   double GetFaceNbrElementSize(int i, int type=0);
+
    /// Return the number of shared faces (3D), edges (2D), vertices (1D)
    int GetNSharedFaces() const;
 

miniapps/meshing/mesh-explorer.cpp

@@ -351,6 +351,7 @@ int main (int argc, char *argv[])
            "e) View elements\n"
            "h) View element sizes, h\n"
            "k) View element ratios, kappa\n"
+           "J) View scaled Jacobian\n"
            "l) Plot a function\n"
            "x) Print sub-element stats\n"
            "f) Find physical point in reference space\n"
@@ -558,6 +559,8 @@ int main (int argc, char *argv[])
          cin >> sd;
          Array<int> bad_elems_by_geom(Geometry::NumGeom);
          bad_elems_by_geom = 0;
+         // Only print so many to keep output compact
+         const int max_to_print = 10;
          for (int i = 0; i < mesh->GetNE(); i++)
          {
             Geometry::Type geom = mesh->GetElementBaseGeometry(i);
@@ -589,10 +592,23 @@ int main (int argc, char *argv[])
             max_det_J = fmax(max_det_J, max_det_J_z);
             if (min_det_J_z <= 0.0)
             {
+               if (nz < max_to_print)
+               {
+                  Vector center;
+                  mesh->GetElementCenter(i, center);
+                  cout << "det(J) < 0 = " << min_det_J_z << " in element "
+                       << i << ", centered at: ";
+                  center.Print();
+               }
                nz++;
                bad_elems_by_geom[geom]++;
             }
          }
+         if (nz >= max_to_print)
+         {
+            cout << "det(J) < 0 for " << nz - max_to_print << " more elements "
+                 << "not printed.\n";
+         }
          cout << "\nbad elements = " << nz;
          if (nz)
          {
@@ -692,7 +708,7 @@ int main (int argc, char *argv[])
 
       // These are most of the cases that open a new GLVis window
       if (mk == 'm' || mk == 'b' || mk == 'e' || mk == 'v' || mk == 'h' ||
-          mk == 'k' || mk == 'p')
+          mk == 'k' || mk == 'J' || mk == 'p')
       {
          Array<int> bdr_part;
          Array<int> part(mesh->GetNE());
@@ -765,7 +781,7 @@ int main (int argc, char *argv[])
             h_max = -h_min;
             for (int i = 0; i < mesh->GetNE(); i++)
             {
-               int geom = mesh->GetElementBaseGeometry(i);
+               Geometry::Type geom = mesh->GetElementBaseGeometry(i);
                ElementTransformation *T = mesh->GetElementTransformation(i);
                T->SetIntPoint(&Geometries.GetCenter(geom));
                Geometries.JacToPerfJac(geom, T->Jacobian(), J);
@@ -790,7 +806,7 @@ int main (int argc, char *argv[])
             DenseMatrix J(dim);
             for (int i = 0; i < mesh->GetNE(); i++)
             {
-               int geom = mesh->GetElementBaseGeometry(i);
+               Geometry::Type geom = mesh->GetElementBaseGeometry(i);
                ElementTransformation *T = mesh->GetElementTransformation(i);
                T->SetIntPoint(&Geometries.GetCenter(geom));
                Geometries.JacToPerfJac(geom, T->Jacobian(), J);
@@ -798,6 +814,48 @@ int main (int argc, char *argv[])
             }
          }
 
+         if (mk == 'J')
+         {
+            // The "scaled Jacobian" is the determinant of the Jacobian scaled
+            // by the l2 norms of its columns. It can be used to identify badly
+            // skewed elements, since it takes values between 0 and 1, with 0
+            // corresponding to a flat element, and 1 to orthogonal columns.
+            DenseMatrix J(dim);
+            int sd;
+            cout << "subdivision factor ---> " << flush;
+            cin >> sd;
+            for (int i = 0; i < mesh->GetNE(); i++)
+            {
+               Geometry::Type geom = mesh->GetElementBaseGeometry(i);
+               ElementTransformation *T = mesh->GetElementTransformation(i);
+
+               RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
+               IntegrationRule &ir = RefG->RefPts;
+
+               // For each element, find the minimal scaled Jacobian in a
+               // lattice of points with the given subdivision factor.
+               attr(i) = infinity();
+               for (int j = 0; j < ir.GetNPoints(); j++)
+               {
+                  T->SetIntPoint(&ir.IntPoint(j));
+                  Geometries.JacToPerfJac(geom, T->Jacobian(), J);
+
+                  // Jacobian determinant
+                  double sJ = J.Det();
+
+                  for (int k = 0; k < J.Width(); k++)
+                  {
+                     Vector col;
+                     J.GetColumnReference(k,col);
+                     // Scale by column norms
+                     sJ /= col.Norml2();
+                  }
+
+                  attr(i) = fmin(sJ, attr(i));
+               }
+            }
+         }
+
          if (mk == 'p')
          {
             int *partitioning = NULL, np;
@@ -950,7 +1008,7 @@ int main (int argc, char *argv[])
             else
             {
                sol_sock << "fem3d_gf_data_keys\n";
-               if (mk == 'v' || mk == 'h' || mk == 'k')
+               if (mk == 'v' || mk == 'h' || mk == 'k' || mk == 'J')
                {
                   mesh->Print(sol_sock);
                }

miniapps/navier/navier_solver.cpp

@@ -290,6 +290,21 @@ void NavierSolver::Setup(double dt)
 
    un_gf.GetTrueDofs(un);
 
+   if (filter_alpha != 0.0)
+   {
+      vfec_filter = new H1_FECollection(order - filter_cutoff_modes,
+                                        pmesh->Dimension());
+      vfes_filter = new ParFiniteElementSpace(pmesh,
+                                              vfec_filter,
+                                              pmesh->Dimension());
+
+      un_NM1_gf.SetSpace(vfes_filter);
+      un_NM1_gf = 0.0;
+
+      un_filtered_gf.SetSpace(vfes);
+      un_filtered_gf = 0.0;
+   }
+
    sw_setup.Stop();
 }
 
@@ -518,6 +533,18 @@ void NavierSolver::Step(double &time, double dt, int cur_step)
 
    un_gf.GetTrueDofs(un);
 
+   if (filter_alpha != 0.0)
+   {
+      un_NM1_gf.ProjectGridFunction(un_gf);
+      un_filtered_gf.ProjectGridFunction(un_NM1_gf);
+      const auto d_un_filtered_gf = un_filtered_gf.Read();
+      auto d_un_gf = un_gf.ReadWrite();
+      MFEM_FORALL(i,
+                  un_gf.Size(),
+                  d_un_gf[i] = (1.0 - filter_alpha) * d_un_gf[i]
+                               + filter_alpha * d_un_filtered_gf[i];);
+   }
+
    sw_step.Stop();
 
    if (verbose && pmesh->GetMyRank() == 0)
@@ -1075,4 +1102,6 @@ NavierSolver::~NavierSolver()
    delete pfec;
    delete vfes;
    delete pfes;
+   delete vfec_filter;
+   delete vfes_filter;
 }

miniapps/navier/navier_solver.hpp

@@ -241,6 +241,19 @@ public:
    /// Compute CFL
    double ComputeCFL(ParGridFunction &u, double dt);
 
+   /// Set the number of modes to cut off in the interpolation filter
+   void SetCutoffModes(int c) { filter_cutoff_modes = c; }
+
+   /// Set the interpolation filter parameter @a a
+   /**
+    * If @a a is > 0, the filtering algorithm for the velocity field after every
+    * time step from [1] is used. The parameter should be 0 > @a >= 1.
+    *
+    * [1] Paul Fischer, Julia Mullen (2001) Filter-based stabilization of
+    * spectral element methods
+    */
+   void SetFilterAlpha(double a) { filter_alpha = a; }
+
 protected:
    /// Print informations about the Navier version.
    void PrintInfo();
@@ -412,6 +425,14 @@ protected:
    OperatorHandle Mv_lor;
    OperatorHandle Sp_lor;
    OperatorHandle H_lor;
+
+   // Filter-based stabilization
+   int filter_cutoff_modes = 1;
+   double filter_alpha = 0.0;
+   FiniteElementCollection *vfec_filter = nullptr;
+   ParFiniteElementSpace *vfes_filter = nullptr;
+   ParGridFunction un_NM1_gf;
+   ParGridFunction un_filtered_gf;
 };
 
 } // namespace navier

tests/convergence/prates.cpp

@@ -26,6 +26,7 @@
 //               mpirun -np 4 prates -m ../../data/inline-hex.mesh -sr 0 -pr 1 -prob 1 -o 3
 //               mpirun -np 4 prates -m ../../data/square-disc.mesh -sr 1 -pr 2 -prob 1 -o 2
 //               mpirun -np 4 prates -m ../../data/star.mesh -sr 1 -pr 2 -prob 3 -o 2
+//               mpirun -np 4 prates -m ../../data/star.mesh -sr 1 -pr 2 -prob 3 -o 2 -j 0
 //               mpirun -np 4 prates -m ../../data/inline-hex.mesh -sr 1 -pr 1 -prob 3 -o 2
 //
 // Description:  This example code demonstrates the use of MFEM to define and
@@ -80,6 +81,7 @@ int main(int argc, char *argv[])
    bool visualization = 1;
    int sr = 1;
    int pr = 1;
+   int jump_scaling_type = 1;
    double sigma = -1.0;
    double kappa = -1.0;
    OptionsParser args(argc, argv);
@@ -95,6 +97,9 @@ int main(int argc, char *argv[])
    args.AddOption(&kappa, "-k", "--kappa",
                   "One of the two DG penalty parameters, should be positive."
                   " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&jump_scaling_type, "-j", "--jump-scaling",
+                  "Scaling of the jump error for DG methods: "
+                  "0: no scaling, 1: 1/h, 2: p^2/h");
    args.AddOption(&sr, "-sr", "--serial_ref",
                   "Number of serial refinements.");
    args.AddOption(&pr, "-pr", "--parallel_ref",
@@ -285,12 +290,15 @@ int main(int argc, char *argv[])
       delete solver;
 
       x = *X;
+      JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
+                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                        : JumpScaling::CONSTANT);
       switch (prob)
       {
          case 0: rates.AddH1GridFunction(&x,scalar_u,gradu); break;
          case 1: rates.AddHcurlGridFunction(&x,vector_u,curlu); break;
          case 2: rates.AddHdivGridFunction(&x,vector_u,divu);  break;
-         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one); break;
+         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one,js); break;
       }
 
       delete X;

tests/convergence/rates.cpp

@@ -26,6 +26,7 @@
 //               rates -m ../../data/inline-hex.mesh -sr 1 -prob 1 -o 2
 //               rates -m ../../data/square-disc.mesh -sr 2 -prob 1 -o 1
 //               rates -m ../../data/star.mesh -sr 2 -prob 3 -o 2
+//               rates -m ../../data/star.mesh -sr 2 -prob 3 -o 2 -j 0
 //               rates -m ../../data/inline-hex.mesh -sr 1 -prob 3 -o 1
 //
 // Description:  This example code demonstrates the use of MFEM to define and
@@ -73,6 +74,7 @@ int main(int argc, char *argv[])
    int order = 1;
    bool visualization = 1;
    int sr = 1;
+   int jump_scaling_type = 1;
    double sigma = -1.0;
    double kappa = -1.0;
    OptionsParser args(argc, argv);
@@ -88,6 +90,9 @@ int main(int argc, char *argv[])
    args.AddOption(&kappa, "-k", "--kappa",
                   "One of the two DG penalty parameters, should be positive."
                   " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&jump_scaling_type, "-j", "--jump-scaling",
+                  "Scaling of the jump error for DG methods: "
+                  "0: no scaling, 1: 1/h, 2: p^2/h");
    args.AddOption(&sr, "-sr", "--serial_ref",
                   "Number of serial refinements.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -218,12 +223,16 @@ int main(int argc, char *argv[])
          PCG(A, M, b, x, 0, 500, 1e-12, 0.0);
       }
 
+      JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
+                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                        : JumpScaling::CONSTANT);
+
       switch (prob)
       {
          case 0: rates.AddH1GridFunction(&x,scalar_u,gradu); break;
          case 1: rates.AddHcurlGridFunction(&x,vector_u,curlu); break;
          case 2: rates.AddHdivGridFunction(&x,vector_u,divu);  break;
-         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one); break;
+         case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one,js); break;
       }
 
       if (l==sr) break;

tests/unit/catch.hpp

@@ -1,6 +1,6 @@
 /*
- *  Catch v2.13.0
- *  Generated: 2020-07-12 20:07:49.015950
+ *  Catch v2.13.2
+ *  Generated: 2020-10-07 11:32:53.302017
  *  ----------------------------------------------------------
  *  This file has been merged from multiple headers. Please don't edit it directly
  *  Copyright (c) 2020 Two Blue Cubes Ltd. All rights reserved.
@@ -15,7 +15,7 @@
 
 #define CATCH_VERSION_MAJOR 2
 #define CATCH_VERSION_MINOR 13
-#define CATCH_VERSION_PATCH 0
+#define CATCH_VERSION_PATCH 2
 
 #ifdef __clang__
 #    pragma clang system_header
@@ -132,13 +132,9 @@ namespace Catch {
 
 #endif
 
-#if defined(__cpp_lib_uncaught_exceptions)
-#  define CATCH_INTERNAL_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS
-#endif
-
 // We have to avoid both ICC and Clang, because they try to mask themselves
 // as gcc, and we want only GCC in this block
-#if defined(__GNUC__) && !defined(__clang__) && !defined(__ICC)
+#if defined(__GNUC__) && !defined(__clang__) && !defined(__ICC) && !defined(__CUDACC__)
 #    define CATCH_INTERNAL_START_WARNINGS_SUPPRESSION _Pragma( "GCC diagnostic push" )
 #    define CATCH_INTERNAL_STOP_WARNINGS_SUPPRESSION  _Pragma( "GCC diagnostic pop" )
 
@@ -162,7 +158,7 @@ namespace Catch {
 // ```
 //
 // Therefore, `CATCH_INTERNAL_IGNORE_BUT_WARN` is not implemented.
-#  if !defined(__ibmxl__)
+#  if !defined(__ibmxl__) && !defined(__CUDACC__)
 #    define CATCH_INTERNAL_IGNORE_BUT_WARN(...) (void)__builtin_constant_p(__VA_ARGS__) /* NOLINT(cppcoreguidelines-pro-type-vararg, hicpp-vararg) */
 #  endif
 
@@ -244,10 +240,6 @@ namespace Catch {
 #  define CATCH_INTERNAL_START_WARNINGS_SUPPRESSION __pragma( warning(push) )
 #  define CATCH_INTERNAL_STOP_WARNINGS_SUPPRESSION  __pragma( warning(pop) )
 
-#  if _MSC_VER >= 1900 // Visual Studio 2015 or newer
-#    define CATCH_INTERNAL_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS
-#  endif
-
 // Universal Windows platform does not support SEH
 // Or console colours (or console at all...)
 #  if defined(WINAPI_FAMILY) && (WINAPI_FAMILY == WINAPI_FAMILY_APP)
@@ -330,7 +322,10 @@ namespace Catch {
 
   // Check if byte is available and usable
   #  if __has_include(<cstddef>) && defined(CATCH_CPP17_OR_GREATER)
-  #    define CATCH_INTERNAL_CONFIG_CPP17_BYTE
+  #    include <cstddef>
+  #    if __cpp_lib_byte > 0
+  #      define CATCH_INTERNAL_CONFIG_CPP17_BYTE
+  #    endif
   #  endif // __has_include(<cstddef>) && defined(CATCH_CPP17_OR_GREATER)
 
   // Check if variant is available and usable
@@ -373,10 +368,6 @@ namespace Catch {
 #  define CATCH_CONFIG_CPP17_OPTIONAL
 #endif
 
-#if defined(CATCH_INTERNAL_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS) && !defined(CATCH_CONFIG_NO_CPP17_UNCAUGHT_EXCEPTIONS) && !defined(CATCH_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS)
-#  define CATCH_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS
-#endif
-
 #if defined(CATCH_INTERNAL_CONFIG_CPP17_STRING_VIEW) && !defined(CATCH_CONFIG_NO_CPP17_STRING_VIEW) && !defined(CATCH_CONFIG_CPP17_STRING_VIEW)
 #  define CATCH_CONFIG_CPP17_STRING_VIEW
 #endif
@@ -1105,7 +1096,7 @@ struct AutoReg : NonCopyable {
                     int index = 0;                                    \
                     constexpr char const* tmpl_types[] = {CATCH_REC_LIST(INTERNAL_CATCH_STRINGIZE_WITHOUT_PARENS, __VA_ARGS__)};\
                     using expander = int[];\
-                    (void)expander{(reg_test(Types{}, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index]), Tags } ), index++, 0)... };/* NOLINT */ \
+                    (void)expander{(reg_test(Types{}, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index]), Tags } ), index++)... };/* NOLINT */ \
                 }\
             };\
             static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1151,7 +1142,7 @@ struct AutoReg : NonCopyable {
                     constexpr char const* tmpl_types[] = {CATCH_REC_LIST(INTERNAL_CATCH_STRINGIZE_WITHOUT_PARENS, INTERNAL_CATCH_REMOVE_PARENS(TmplTypes))};\
                     constexpr char const* types_list[] = {CATCH_REC_LIST(INTERNAL_CATCH_STRINGIZE_WITHOUT_PARENS, INTERNAL_CATCH_REMOVE_PARENS(TypesList))};\
                     constexpr auto num_types = sizeof(types_list) / sizeof(types_list[0]);\
-                    (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestFuncName<Types> ), CATCH_INTERNAL_LINEINFO, Catch::StringRef(), Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index / num_types]) + "<" + std::string(types_list[index % num_types]) + ">", Tags } ), index++, 0)... };/* NOLINT */\
+                    (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestFuncName<Types> ), CATCH_INTERNAL_LINEINFO, Catch::StringRef(), Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index / num_types]) + "<" + std::string(types_list[index % num_types]) + ">", Tags } ), index++)... };/* NOLINT */\
                 }                                                     \
             };                                                        \
             static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){ \
@@ -1195,7 +1186,7 @@ struct AutoReg : NonCopyable {
             void reg_tests() {                                          \
                 int index = 0;                                    \
                 using expander = int[];                           \
-                (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestFunc<Types> ), CATCH_INTERNAL_LINEINFO, Catch::StringRef(), Catch::NameAndTags{ Name " - " + std::string(INTERNAL_CATCH_STRINGIZE(TmplList)) + " - " + std::to_string(index), Tags } ), index++, 0)... };/* NOLINT */\
+                (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestFunc<Types> ), CATCH_INTERNAL_LINEINFO, Catch::StringRef(), Catch::NameAndTags{ Name " - " + std::string(INTERNAL_CATCH_STRINGIZE(TmplList)) + " - " + std::to_string(index), Tags } ), index++)... };/* NOLINT */\
             }                                                     \
         };\
         static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){ \
@@ -1229,7 +1220,7 @@ struct AutoReg : NonCopyable {
                     int index = 0;                                    \
                     constexpr char const* tmpl_types[] = {CATCH_REC_LIST(INTERNAL_CATCH_STRINGIZE_WITHOUT_PARENS, __VA_ARGS__)};\
                     using expander = int[];\
-                    (void)expander{(reg_test(Types{}, #ClassName, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index]), Tags } ), index++, 0)... };/* NOLINT */ \
+                    (void)expander{(reg_test(Types{}, #ClassName, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index]), Tags } ), index++)... };/* NOLINT */ \
                 }\
             };\
             static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1278,7 +1269,7 @@ struct AutoReg : NonCopyable {
                     constexpr char const* tmpl_types[] = {CATCH_REC_LIST(INTERNAL_CATCH_STRINGIZE_WITHOUT_PARENS, INTERNAL_CATCH_REMOVE_PARENS(TmplTypes))};\
                     constexpr char const* types_list[] = {CATCH_REC_LIST(INTERNAL_CATCH_STRINGIZE_WITHOUT_PARENS, INTERNAL_CATCH_REMOVE_PARENS(TypesList))};\
                     constexpr auto num_types = sizeof(types_list) / sizeof(types_list[0]);\
-                    (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestName<Types>::test ), CATCH_INTERNAL_LINEINFO, #ClassName, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index / num_types]) + "<" + std::string(types_list[index % num_types]) + ">", Tags } ), index++, 0)... };/* NOLINT */ \
+                    (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestName<Types>::test ), CATCH_INTERNAL_LINEINFO, #ClassName, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index / num_types]) + "<" + std::string(types_list[index % num_types]) + ">", Tags } ), index++)... };/* NOLINT */ \
                 }\
             };\
             static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1325,7 +1316,7 @@ struct AutoReg : NonCopyable {
                 void reg_tests(){\
                     int index = 0;\
                     using expander = int[];\
-                    (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestName<Types>::test ), CATCH_INTERNAL_LINEINFO, #ClassName, Catch::NameAndTags{ Name " - " + std::string(INTERNAL_CATCH_STRINGIZE(TmplList)) + " - " + std::to_string(index), Tags } ), index++, 0)... };/* NOLINT */ \
+                    (void)expander{(Catch::AutoReg( Catch::makeTestInvoker( &TestName<Types>::test ), CATCH_INTERNAL_LINEINFO, #ClassName, Catch::NameAndTags{ Name " - " + std::string(INTERNAL_CATCH_STRINGIZE(TmplList)) + " - " + std::to_string(index), Tags } ), index++)... };/* NOLINT */ \
                 }\
             };\
             static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1829,8 +1820,8 @@ namespace Catch {
 #endif
 
     namespace Detail {
-        template<typename InputIterator>
-        std::string rangeToString(InputIterator first, InputIterator last) {
+        template<typename InputIterator, typename Sentinel = InputIterator>
+        std::string rangeToString(InputIterator first, Sentinel last) {
             ReusableStringStream rss;
             rss << "{ ";
             if (first != last) {
@@ -12042,7 +12033,7 @@ namespace Catch {
         if (tmpnam_s(m_buffer)) {
             CATCH_RUNTIME_ERROR("Could not get a temp filename");
         }
-        if (fopen_s(&m_file, m_buffer, "w")) {
+        if (fopen_s(&m_file, m_buffer, "w+")) {
             char buffer[100];
             if (strerror_s(buffer, errno)) {
                 CATCH_RUNTIME_ERROR("Could not translate errno to a string");
@@ -14091,10 +14082,10 @@ namespace Catch {
 
     namespace {
         struct TestHasher {
-            explicit TestHasher(Catch::SimplePcg32& rng) {
-                basis = rng();
+            explicit TestHasher(Catch::SimplePcg32& rng_instance) {
+                basis = rng_instance();
                 basis <<= 32;
-                basis |= rng();
+                basis |= rng_instance();
             }
 
             uint64_t basis;
@@ -15192,6 +15183,41 @@ namespace Catch {
 // end catch_totals.cpp
 // start catch_uncaught_exceptions.cpp
 
+// start catch_config_uncaught_exceptions.hpp
+
+//              Copyright Catch2 Authors
+// Distributed under the Boost Software License, Version 1.0.
+//   (See accompanying file LICENSE_1_0.txt or copy at
+//        https://www.boost.org/LICENSE_1_0.txt)
+
+// SPDX-License-Identifier: BSL-1.0
+
+#ifndef CATCH_CONFIG_UNCAUGHT_EXCEPTIONS_HPP
+#define CATCH_CONFIG_UNCAUGHT_EXCEPTIONS_HPP
+
+#if defined(_MSC_VER)
+#  if _MSC_VER >= 1900 // Visual Studio 2015 or newer
+#    define CATCH_INTERNAL_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS
+#  endif
+#endif
+
+#include <exception>
+
+#if defined(__cpp_lib_uncaught_exceptions) \
+    && !defined(CATCH_INTERNAL_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS)
+
+#  define CATCH_INTERNAL_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS
+#endif // __cpp_lib_uncaught_exceptions
+
+#if defined(CATCH_INTERNAL_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS) \
+    && !defined(CATCH_CONFIG_NO_CPP17_UNCAUGHT_EXCEPTIONS) \
+    && !defined(CATCH_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS)
+
+#  define CATCH_CONFIG_CPP17_UNCAUGHT_EXCEPTIONS
+#endif
+
+#endif // CATCH_CONFIG_UNCAUGHT_EXCEPTIONS_HPP
+// end catch_config_uncaught_exceptions.hpp
 #include <exception>
 
 namespace Catch {
@@ -15238,7 +15264,7 @@ namespace Catch {
     }
 
     Version const& libraryVersion() {
-        static Version version( 2, 13, 0, "", 0 );
+        static Version version( 2, 13, 2, "", 0 );
         return version;
     }
 
@@ -17796,3 +17822,4 @@ using Catch::Detail::Approx;
 // end catch_reenable_warnings.h
 // end catch.hpp
 #endif // TWOBLUECUBES_SINGLE_INCLUDE_CATCH_HPP_INCLUDED
+

tests/unit/linalg/test_matrix_block.cpp

@@ -10,9 +10,9 @@
 // CONTRIBUTING.md for details.
 
 #include <fstream>
-using namespace std;
-
 #include "mfem.hpp"
+
+using namespace std;
 using namespace mfem;
 
 #include "unit_tests.hpp"

tests/unit/linalg/test_ode.cpp

@@ -251,6 +251,24 @@ TEST_CASE("First order ODE methods",
       REQUIRE(conv_rate + tol > 4.0);
    }
 
+   SECTION("TrapezoidalRuleSolver")
+   {
+      std::cout <<"\nTesting TrapezoidalRuleSolver" << std::endl;
+      REQUIRE(check.order(new TrapezoidalRuleSolver) + tol > 2.0 );
+   }
+
+   SECTION("ESDIRK32Solver")
+   {
+      std::cout <<"\nTesting ESDIRK32Solver" << std::endl;
+      REQUIRE(check.order(new ESDIRK32Solver) + tol > 2.0 );
+   }
+
+   SECTION("ESDIRK33Solver")
+   {
+      std::cout <<"\nTesting ESDIRK33Solver" << std::endl;
+      REQUIRE(check.order(new ESDIRK33Solver) + tol > 3.0 );
+   }
+
    // Generalized-alpha
    SECTION("GeneralizedAlphaSolver(1.0)")
    {

tests/unit/miniapps/test_sedov.cpp

@@ -98,14 +98,14 @@ struct Tensors1D
 
 template<int DIM, int D1D, int Q1D, int L1D, int H1D, int NBZ =1> static
 void kSmemForceMult2D(const int NE,
-                      const Array<double> &_B,
+                      const Array<double> &B_,
                       const Array<double> &_Bt,
                       const Array<double> &_Gt,
                       const DenseTensor &_sJit,
                       const Vector &_e,
                       Vector &_v)
 {
-   auto b = Reshape(_B.Read(), Q1D, L1D);
+   auto b = Reshape(B_.Read(), Q1D, L1D);
    auto bt = Reshape(_Bt.Read(), H1D, Q1D);
    auto gt = Reshape(_Gt.Read(), H1D, Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*NE*2*2),
@@ -241,14 +241,14 @@ void kSmemForceMult2D(const int NE,
 
 template<int DIM, int D1D, int Q1D, int L1D, int H1D> static
 void kSmemForceMult3D(const int NE,
-                      const Array<double> &_B,
+                      const Array<double> &B_,
                       const Array<double> &_Bt,
                       const Array<double> &_Gt,
                       const DenseTensor &_sJit,
                       const Vector &_e,
                       Vector &_v)
 {
-   auto b = Reshape(_B.Read(), Q1D, L1D);
+   auto b = Reshape(B_.Read(), Q1D, L1D);
    auto bt = Reshape(_Bt.Read(), H1D, Q1D);
    auto gt = Reshape(_Gt.Read(), H1D, Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*Q1D*NE*3*3),
@@ -493,14 +493,14 @@ static void kForceMult(const int DIM,
 template<int DIM, int D1D, int Q1D, int L1D, int H1D, int NBZ =1> static
 void kSmemForceMultTranspose2D(const int NE,
                                const Array<double> &_Bt,
-                               const Array<double> &_B,
+                               const Array<double> &B_,
                                const Array<double> &_G,
                                const DenseTensor &_sJit,
                                const Vector &_v,
                                Vector &_e)
 {
    MFEM_VERIFY(D1D==H1D,"");
-   auto b = Reshape(_B.Read(), Q1D,H1D);
+   auto b = Reshape(B_.Read(), Q1D,H1D);
    auto g = Reshape(_G.Read(), Q1D,H1D);
    auto bt = Reshape(_Bt.Read(), L1D,Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*NE*2*2),
@@ -633,14 +633,14 @@ void kSmemForceMultTranspose2D(const int NE,
 template<int DIM, int D1D, int Q1D, int L1D, int H1D> static
 void kSmemForceMultTranspose3D(const int NE,
                                const Array<double> &_Bt,
-                               const Array<double> &_B,
+                               const Array<double> &B_,
                                const Array<double> &_G,
                                const DenseTensor &_sJit,
                                const Vector &_v,
                                Vector &_e)
 {
    MFEM_VERIFY(D1D==H1D,"");
-   auto b = Reshape(_B.Read(), Q1D,H1D);
+   auto b = Reshape(B_.Read(), Q1D,H1D);
    auto g = Reshape(_G.Read(), Q1D,H1D);
    auto bt = Reshape(_Bt.Read(), L1D,Q1D);
    auto sJit = Reshape(Read(_sJit.GetMemory(), Q1D*Q1D*Q1D*NE*3*3),