renaming and alternative interface

example test
support coordinates in qfunc kernel
2020-12-03 08:49:51 -08:00 · 2020-11-30 13:43:45 -08:00 · 2020-11-25 09:59:55 -08:00 · 2020-11-24 14:29:27 -08:00 · 2020-11-24 12:05:45 -08:00 · 2020-11-23 18:10:52 -08:00
37 changed files with 2451 additions and 53 deletions
@@ -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
 =========================================

@@ -14,7 +14,7 @@ set(USER_CONFIG "${CMAKE_CURRENT_SOURCE_DIR}/config/user.cmake" CACHE PATH
  "Path to optional user configuration file.")

 # Require C++11 and disable compiler-specific extensions
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)

@@ -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.
@@ -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
@@ -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@)
@@ -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

@@ -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)
@@ -42,6 +42,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
@@ -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

@@ -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@
@@ -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.
@@ -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)
@@ -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
+}
@@ -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|" && \
@@ -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
@@ -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.
@@ -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)
@@ -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
@@ -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
@@ -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
@@ -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)
 {
@@ -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.
@@ -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))
@@ -26,3 +26,4 @@ add_subdirectory(tools)
 add_subdirectory(toys)
 add_subdirectory(nurbs)
 add_subdirectory(gslib)
+add_subdirectory(variationalform)
@@ -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);
               }
@@ -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;
 }
@@ -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
@@ -0,0 +1,15 @@
+add_library(parvariationalform parvariationalform.cpp)
+target_include_directories(parvariationalform PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+target_link_libraries(parvariationalform mfem)
+
+add_executable(parvariationalform_example parvariationalform_example.cpp)
+target_include_directories(parvariationalform_example PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+target_link_libraries(parvariationalform_example parvariationalform mfem)
+
+add_executable(parvariationalform_ex1 parvariationalform_ex1.cpp)
+target_include_directories(parvariationalform_ex1 PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+target_link_libraries(parvariationalform_ex1 parvariationalform mfem)
+
+add_executable(parvariationalform_unit_tests test_tensor_ad.cpp)
+target_include_directories(parvariationalform_unit_tests PRIVATE ${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_SOURCE_DIR}/tests/unit)
+target_link_libraries(parvariationalform_unit_tests parvariationalform mfem)
@@ -0,0 +1,23 @@
+#include "mfem.hpp"
+
+#pragma once
+
+namespace mfem
+{
+class GenericIntegrator
+{
+public:
+   GenericIntegrator(const IntegrationRule *ir = nullptr) : IntRule(ir) {}
+
+   virtual void Setup(const FiniteElementSpace &) = 0;
+
+   virtual void Apply(const Vector &x, Vector &y) const = 0;
+
+   virtual void ApplyGradient(const Vector &x,
+                              const Vector &v,
+                              Vector &y) const = 0;
+
+protected:
+   const IntegrationRule *IntRule;
+};
+} // namespace mfem
@@ -0,0 +1,94 @@
+#include "parvariationalform.hpp"
+#include "qfuncintegrator.hpp"
+
+namespace mfem
+{
+ParVariationalForm::ParVariationalForm(ParFiniteElementSpace *f)
+   : Operator(f->GetTrueVSize()), fes(f), P(f->GetProlongationMatrix()),
+     grad(*this)
+{
+   G = fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+   MFEM_ASSERT(G, "Some GetElementRestriction error");
+   x_local.SetSize(G->Height(), Device::GetMemoryType());
+   v_local.SetSize(G->Height(), Device::GetMemoryType());
+   y_local.SetSize(G->Height(), Device::GetMemoryType());
+}
+
+void ParVariationalForm::Mult(const Vector &x, Vector &y) const
+{
+   px.SetSize(P->Height());
+   py.SetSize(P->Height());
+
+   P->Mult(x, px);
+   G->Mult(px, x_local);
+
+   y_local = 0.0;
+   for (int i = 0; i < domain_integrators.Size(); ++i)
+   {
+      // y += F(x)
+      domain_integrators[i]->Apply(x_local, y_local);
+   }
+
+   G->MultTranspose(y_local, py);
+   P->MultTranspose(py, y);
+
+   y.HostReadWrite();
+   for (int i = 0; i < ess_tdof_list.Size(); i++)
+   {
+      y(ess_tdof_list[i]) = 0.0;
+   }
+}
+
+void ParVariationalForm::GradientMult(const Vector &v, Vector &y) const
+{
+   px.SetSize(P->Height());
+   py.SetSize(P->Height());
+   pv.SetSize(P->Height());
+
+   P->Mult(v, pv);
+   G->Mult(pv, v_local);
+
+   P->Mult(x_lin, px);
+   G->Mult(px, x_local);
+
+   y_local = 0.0;
+   for (int i = 0; i < domain_integrators.Size(); ++i)
+   {
+      // y += dF(x)/dx * v
+      if (is_linear)
+      {
+         // take care of RHS
+         // domain_integrators[i]->Apply(v_local, y_local);
+      }
+      domain_integrators[i]->ApplyGradient(x_local, v_local, y_local);
+   }
+
+   G->MultTranspose(y_local, py);
+   P->MultTranspose(py, y);
+
+   y.HostReadWrite();
+   for (int i = 0; i < ess_tdof_list.Size(); i++)
+   {
+      y(ess_tdof_list[i]) = v(ess_tdof_list[i]);
+   }
+}
+
+Operator &ParVariationalForm::GetGradient(const Vector &x) const
+{
+   x_lin = x;
+   return grad;
+}
+
+HypreParMatrix *ParVariationalForm::GetGradientMatrix(const Vector &x)
+{
+   delete gradient_matrix;
+   gradient_matrix = new HypreParMatrix;
+   return gradient_matrix;
+}
+
+void ParVariationalForm::SetEssentialBC(const Array<int> &ess_attr)
+{
+   fes->GetEssentialTrueDofs(ess_attr, ess_tdof_list);
+}
+
+} // namespace mfem
@@ -0,0 +1,109 @@
+#include "mfem.hpp"
+#include "genericintegrator.hpp"
+#include "qfuncintegrator.hpp"
+
+#pragma once
+
+namespace mfem
+{
+class ParVariationalForm : public Operator
+{
+   class Gradient : public Operator
+   {
+   public:
+      Gradient(ParVariationalForm &f) : Operator(f.Height()), form(f){};
+
+      void Mult(const Vector &x, Vector &y) const override
+      {
+         form.GradientMult(x, y);
+      }
+
+   private:
+      ParVariationalForm &form;
+   };
+
+public:
+   ParVariationalForm(ParFiniteElementSpace *f);
+
+   void Mult(const Vector &x, Vector &y) const override;
+
+   void AddDomainIntegrator(GenericIntegrator *i)
+   {
+      domain_integrators.Append(i);
+      i->Setup(*fes);
+   }
+
+   template<typename integrator_type,
+            typename qfunc_type,
+            typename... qfunc_args_type>
+   void AddDomainIntegrator(qfunc_type f, qfunc_args_type const &... fargs)
+   {
+      if constexpr (std::is_same_v<integrator_type, DomainLFIntegrator>)
+      {
+         auto i = new QFunctionIntegrator(
+            [&](auto... args) {
+               auto du = std::get<1>(std::tuple{args...});
+               return std::tuple{f(args...), decltype(du){}};
+            },
+            0,
+            *fes->GetParMesh(),
+            fargs...);
+
+         domain_integrators.Append(i);
+         i->Setup(*fes);
+      }
+
+      if constexpr (std::is_same_v<integrator_type, DiffusionIntegrator>)
+      {
+         auto i = new QFunctionIntegrator(
+            [&](auto... args) {
+               auto du = std::get<1>(std::tuple{args...});
+               return std::tuple{0.0, f(args...) * du};
+            },
+            0,
+            *fes->GetParMesh(),
+            fargs...);
+
+         domain_integrators.Append(i);
+         i->Setup(*fes);
+      }
+   }
+
+   void AssumeLinear() { is_linear = true; }
+
+   // Return an Operator that provides a Mult(x, y) which is the MatVec of the
+   // gradient of the ParVariationalForm wrt x. Acts as a "passthrough" to
+   // ::GradientMult in order to satisfy mfem interfaces.
+   Operator &GetGradient(const Vector &x) const override;
+
+   // Return an assmbled parallel matrix which represents the gradient of the
+   // ParVariationalForm wrt to x.
+   HypreParMatrix *GetGradientMatrix(const Vector &x);
+
+   void SetEssentialBC(const Array<int> &ess_attr);
+
+protected:
+   // y = F'(x_lin) * v
+   void GradientMult(const Vector &v, Vector &y) const;
+
+   bool is_linear = false;
+   ParFiniteElementSpace *fes;
+   Array<GenericIntegrator *> domain_integrators;
+   Array<int> ess_tdof_list;
+
+   // T -> L
+   const Operator *P;
+
+   // L -> E
+   const Operator *G;
+
+   mutable Vector x_local, y_local, v_local, px, py, pv;
+
+   // State to build the Gradient on, single source of the true state.
+   mutable Vector x_lin;
+
+   mutable Gradient grad;
+
+   HypreParMatrix *gradient_matrix = nullptr;
+};
+} // namespace mfem
@@ -0,0 +1,164 @@
+#include "mfem.hpp"
+#include "parvariationalform.hpp"
+#include "tensor.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order,
+                  "-o",
+                  "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   {
+      int ref_levels = (int) floor(log(10. / mesh.GetNE()) / log(2.) / dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+
+   FiniteElementCollection *fec;
+   bool delete_fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
+   }
+   else if (pmesh.GetNodes())
+   {
+      fec = pmesh.GetNodes()->OwnFEC();
+      delete_fec = false;
+      if (myid == 0)
+      {
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
+   }
+   ParFiniteElementSpace fespace(&pmesh, fec);
+   HYPRE_Int size = fespace.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   FunctionCoefficient u_excoeff([&](const Vector &coords) {
+      double x = coords(0);
+      double y = coords(1);
+      return x * x + y * y;
+   });
+
+   ParGridFunction x(&fespace);
+
+   x.ProjectBdrCoefficient(u_excoeff, ess_bdr);
+
+   ParVariationalForm form(&fespace);
+
+   auto b_coeff = [&](auto u, auto du, auto x) {
+      return 4.0 * (1.0 + 2.0 * x[0] * x[0] + 2.0 * x[1] * x[1]);
+   };
+
+   form.AddDomainIntegrator<DomainLFIntegrator>(b_coeff);
+
+   auto a_coeff = [&](auto u, auto du, auto x) { return 1.0 + u; };
+
+   form.AddDomainIntegrator<DiffusionIntegrator>(a_coeff);
+
+   form.SetEssentialBC(ess_bdr);
+
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-6);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+
+   NewtonSolver newton(MPI_COMM_WORLD);
+   newton.SetOperator(form);
+   newton.SetSolver(cg);
+   newton.SetPrintLevel(1);
+   newton.SetRelTol(1e-8);
+   newton.SetMaxIter(100);
+
+   Vector zero;
+   Vector X;
+   x.GetTrueDofs(X);
+   newton.Mult(zero, X);
+
+   x.Distribute(X);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+   if (delete_fec)
+   {
+      delete fec;
+   }
+   MPI_Finalize();
+
+   return 0;
+}
@@ -0,0 +1,118 @@
+#include "mfem.hpp"
+#include "parvariationalform.hpp"
+#include "qfuncintegrator.hpp"
+#include "tensor.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int order = 1;
+   int refinements = 0;
+   double p = 5.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&refinements, "-r", "--ref", "");
+   args.AddOption(&order, "-o", "--order", "");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   Mesh mesh(mesh_file, 1, 1);
+
+   int dim = mesh.Dimension();
+   {
+      for (int l = 0; l < refinements; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
+   auto fec = H1_FECollection(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   FunctionCoefficient u_excoeff([&](const Vector &coords) {
+      double x = coords(0);
+      double y = coords(1);
+
+      // return 1.0 - pow(sqrt(x * x + y * y), p / (p - 1.0));
+      return x * x + y * y;
+   });
+
+   ParGridFunction x(&fespace);
+   x.Randomize();
+
+   x.ProjectBdrCoefficient(u_excoeff, ess_bdr);
+
+   ParVariationalForm form(&fespace);
+
+   auto plaplacian = new QFunctionIntegrator([&](auto u, auto du, auto x) {
+      // auto f0 = -1.0;
+      // auto f1 = pow(norm(du), p - 2.0) * du;
+      auto f0 = 4.0 * (1.0 + 2.0 * x[0] * x[0] + 2.0 * x[1] * x[1]);
+      auto f1 = (1.0 + u) * du;
+      return std::tuple{f0, f1};
+   }, 0, pmesh);
+
+   form.AddDomainIntegrator(plaplacian);
+
+   form.SetEssentialBC(ess_bdr);
+
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-6);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+
+   NewtonSolver newton(MPI_COMM_WORLD);
+   newton.SetOperator(form);
+   newton.SetSolver(cg);
+   newton.SetPrintLevel(1);
+   newton.SetRelTol(1e-8);
+   newton.SetMaxIter(100);
+
+   Vector zero;
+   Vector X;
+   x.GetTrueDofs(X);
+   newton.Mult(zero, X);
+
+   x.Distribute(X);
+
+   // x.ProjectCoefficient(u_excoeff);
+   
+   char vishost[] = "localhost";
+   int visport = 19916;
+   socketstream sol_sock(vishost, visport);
+   sol_sock << "parallel " << num_procs << " " << myid << "\n";
+   sol_sock.precision(8);
+   sol_sock << "solution\n" << pmesh << x << flush;
+
+   MPI_Finalize();
+
+   return 0;
+}
@@ -0,0 +1,404 @@
+#include "mfem.hpp"
+#include "../../general/forall.hpp"
+#include "genericintegrator.hpp"
+#include "tensor.hpp"
+
+#pragma once
+
+namespace mfem
+{
+template<typename T>
+struct supported_type
+{
+   static constexpr bool value = false;
+};
+
+template<>
+struct supported_type<ParMesh>
+{
+   static constexpr bool value = true;
+};
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+class QFunctionIntegrator : public GenericIntegrator
+{
+public:
+   QFunctionIntegrator(qfunc_type f,
+                       qfunc_grad_type f_grad,
+                       qfunc_args_type const &... fargs);
+
+   QFunctionIntegrator(qfunc_type f, qfunc_args_type const &... fargs);
+
+   void Setup(const FiniteElementSpace &fes) override;
+
+   void Apply(const Vector &, Vector &) const override;
+
+   // y += F'(x) * v
+   void ApplyGradient(const Vector &x,
+                      const Vector &v,
+                      Vector &y) const override;
+
+protected:
+   template<int D1D, int Q1D>
+   void Apply2D(const Vector &u_in_, Vector &y_) const;
+
+   template<int D1D, int Q1D>
+   void ApplyGradient2D(const Vector &u_in_,
+                        const Vector &v_in_,
+                        Vector &y_) const;
+
+   auto EvaluateFargValue(const Mesh &m, const double qx, const double qy) const;
+
+   const FiniteElementSpace *fespace;
+   const DofToQuad *maps;        ///< Not owned
+   const GeometricFactors *geom; ///< Not owned
+   int dim, ne, nq, dofs1D, quad1D;
+
+   // Geometric factors
+   Vector J_;
+   Vector W_;
+
+   qfunc_type qf;
+   qfunc_grad_type qf_grad;
+   std::tuple<qfunc_args_type...> qf_farg_values;
+};
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
+   QFunctionIntegrator(qfunc_type f,
+                       qfunc_grad_type df,
+                       qfunc_args_type const &... fargs)
+   : GenericIntegrator(nullptr), maps(nullptr), geom(nullptr), qf(f),
+     qf_grad(df), qf_farg_values(std::tuple{fargs...})
+{
+   static_assert((supported_type<qfunc_args_type>::value && ...),
+                 "Type not supported for parameter expansion. See "
+                 "documentation for supported types.");
+}
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
+   QFunctionIntegrator(qfunc_type f, qfunc_args_type const &... fargs)
+   : GenericIntegrator(nullptr), maps(nullptr), geom(nullptr), qf(f),
+     qf_grad(qfunc_grad_type{}), qf_farg_values(std::tuple{fargs...})
+{
+   static_assert((supported_type<qfunc_args_type>::value && ...),
+                 "Type not supported for parameter expansion. See "
+                 "documentation for supported types.");
+}
+
+template<typename qfunc_type, typename... qfunc_args_type>
+QFunctionIntegrator(qfunc_type, qfunc_args_type const &...)
+   -> QFunctionIntegrator<qfunc_type, int, qfunc_args_type const &...>;
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::Setup(
+   const FiniteElementSpace &fes)
+{
+   // Assuming the same element type
+   fespace = &fes;
+   Mesh *mesh = fes.GetMesh();
+   if (mesh->GetNE() == 0)
+   {
+      return;
+   }
+   const FiniteElement &el = *fes.GetFE(0);
+   ElementTransformation *T = mesh->GetElementTransformation(0);
+   const IntegrationRule *ir = nullptr;
+   if (!IntRule)
+   {
+      IntRule = &IntRules.Get(el.GetGeomType(), el.GetOrder() * 2);
+   }
+   ir = IntRule;
+
+   dim = mesh->Dimension();
+   ne = fes.GetMesh()->GetNE();
+   nq = ir->GetNPoints();
+   geom = mesh->GetGeometricFactors(*ir,
+                                    GeometricFactors::COORDINATES
+                                       | GeometricFactors::JACOBIANS);
+   maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
+   dofs1D = maps->ndof;
+   quad1D = maps->nqpt;
+   //    pa_data.SetSize(ne * nq, Device::GetDeviceMemoryType());
+
+   W_.SetSize(nq, Device::GetDeviceMemoryType());
+   W_.GetMemory().CopyFrom(ir->GetWeights().GetMemory(), nq);
+
+   // J.SetSize(ne * nq, Device::GetDeviceMemoryType());
+   J_ = geom->J;
+}
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+auto QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
+   EvaluateFargValue(const Mesh &m, const double qx, const double qy) const
+{
+   Vector trip(3);
+   trip = 0.0;
+   ElementTransformation *tr = const_cast<Mesh &>(m).GetElementTransformation(0);
+   tr->Transform(IntRule->IntPoint(qx + quad1D * qy), trip);
+   return tensor<double, 3>{
+      {trip(0), trip(1), m.SpaceDimension() == 2 ? 0.0 : trip(2)}};
+}
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::Apply(
+   const Vector &x, Vector &y) const
+{
+   if (dim == 2)
+   {
+      switch ((dofs1D << 4) | quad1D)
+      {
+      case 0x22:
+         return Apply2D<2, 2>(x, y);
+      default:
+         MFEM_ASSERT(false, "NOPE");
+      }
+   }
+}
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+template<int D1D, int Q1D>
+void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::Apply2D(
+   const Vector &u_in_, Vector &y_) const
+{
+   int NE = ne;
+
+   auto v1d = Reshape(maps->B.Read(), Q1D, D1D);
+   auto dv1d_dX = Reshape(maps->G.Read(), Q1D, D1D);
+   // (NQ x SDIM x DIM x NE)
+   auto J = Reshape(J_.Read(), Q1D, Q1D, 2, 2, NE);
+   auto W = Reshape(W_.Read(), Q1D, Q1D);
+   auto u = Reshape(u_in_.Read(), D1D, D1D, NE);
+   auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
+
+   // MFEM_FORALL(e, NE, {
+   for (int e = 0; e < NE; e++)
+   {
+      // loop over quadrature points
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            double u_q = 0.0;
+            double du_dX_q[2] = {0.0};
+            for (int ix = 0; ix < D1D; ix++)
+            {
+               for (int iy = 0; iy < D1D; iy++)
+               {
+                  u_q += u(ix, iy, e) * v1d(qx, ix) * v1d(qy, iy);
+                  du_dX_q[0] += u(ix, iy, e) * dv1d_dX(qx, ix) * v1d(qy, iy);
+                  du_dX_q[1] += u(ix, iy, e) * v1d(qx, ix) * dv1d_dX(qy, iy);
+               }
+            }
+
+            // du_dx_q = invJ^T * du_dX_q
+            //         = (adjJ^T * du_dX_q) / detJ
+            double J_q[2][2] = {{J(qx, qy, 0, 0, e),
+                                 J(qx, qy, 0, 1, e)}, // J_q[0][0], J_q[0][1]
+                                {J(qx, qy, 1, 0, e),
+                                 J(qx, qy, 1, 1, e)}}; // J_q[1][0], J_q[1][1]
+
+            double detJ_q = (J_q[0][0] * J_q[1][1]) - (J_q[0][1] * J_q[1][0]);
+
+            double adjJ[2][2] = {{J_q[1][1], -J_q[0][1]},
+                                 {-J_q[1][0], J_q[0][0]}};
+
+            tensor<double, 2> du_dx_q
+               = {(adjJ[0][0] * du_dX_q[0] + adjJ[1][0] * du_dX_q[1]) / detJ_q,
+                  (adjJ[0][1] * du_dX_q[0] + adjJ[1][1] * du_dX_q[1]) / detJ_q};
+
+            auto processed_qf_farg_values = std::apply(
+               [=](auto &... a) {
+                  return std::make_tuple(u_q,
+                                         du_dx_q,
+                                         EvaluateFargValue(a, qx, qy)...);
+               },
+               qf_farg_values);
+
+            auto [f0, f1] = std::apply(qf, processed_qf_farg_values);
+
+            double f0_X = f0 * detJ_q;
+
+            // f1_X = invJ * f1 * detJ
+            //      = adjJ * f1
+            double f1_X[2] = {
+               adjJ[0][0] * f1[0] + adjJ[0][1] * f1[1],
+               adjJ[1][0] * f1[0] + adjJ[1][1] * f1[1],
+            };
+
+            for (int ix = 0; ix < D1D; ix++)
+            {
+               for (int iy = 0; iy < D1D; iy++)
+               {
+                  // accumulate v * f0 + dot(dv_dx, f1)
+                  y(ix, iy, e) += (f0_X * v1d(qx, ix) * v1d(qy, iy)
+                                   + f1_X[0] * dv1d_dX(qx, ix) * v1d(qy, iy)
+                                   + f1_X[1] * dv1d_dX(qy, iy) * v1d(qx, ix))
+                                  * W(qx, qy);
+               }
+            }
+         }
+      }
+   }
+   // });
+}
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
+   ApplyGradient(const Vector &x, const Vector &v, Vector &y) const
+{
+   ApplyGradient2D<2, 2>(x, v, y);
+}
+
+template<typename qfunc_type, typename qfunc_grad_type, typename... qfunc_args_type>
+template<int D1D, int Q1D>
+void QFunctionIntegrator<qfunc_type, qfunc_grad_type, qfunc_args_type...>::
+   ApplyGradient2D(const Vector &u_in_, const Vector &v_in_, Vector &y_) const
+{
+   int NE = ne;
+
+   auto v1d = Reshape(maps->B.Read(), Q1D, D1D);
+   auto dv1d_dX = Reshape(maps->G.Read(), Q1D, D1D);
+   // (NQ x SDIM x DIM x NE)
+   auto J = Reshape(J_.Read(), Q1D, Q1D, 2, 2, NE);
+   auto W = Reshape(W_.Read(), Q1D, Q1D);
+   auto u = Reshape(u_in_.Read(), D1D, D1D, NE);
+   auto v = Reshape(v_in_.Read(), D1D, D1D, NE);
+   auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
+
+   for (int e = 0; e < NE; e++)
+   {
+      // loop over quadrature points
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            double u_q = 0.0;
+            double du_dX_q[2] = {0.0};
+            double v_q = 0.0;
+            double dv_dX_q[2] = {0.0};
+            for (int ix = 0; ix < D1D; ix++)
+            {
+               for (int iy = 0; iy < D1D; iy++)
+               {
+                  u_q += u(ix, iy, e) * v1d(qx, ix) * v1d(qy, iy);
+
+                  du_dX_q[0] += u(ix, iy, e) * dv1d_dX(qx, ix) * v1d(qy, iy);
+                  du_dX_q[1] += u(ix, iy, e) * v1d(qx, ix) * dv1d_dX(qy, iy);
+
+                  v_q += v(ix, iy, e) * v1d(qx, ix) * v1d(qy, iy);
+
+                  dv_dX_q[0] += v(ix, iy, e) * dv1d_dX(qx, ix) * v1d(qy, iy);
+                  dv_dX_q[1] += v(ix, iy, e) * v1d(qx, ix) * dv1d_dX(qy, iy);
+               }
+            }
+
+            // du_dx_q = invJ^T * du_dX_q
+            //         = (adjJ^T * du_dX_q) / detJ
+            double J_q[2][2] = {{J(qx, qy, 0, 0, e),
+                                 J(qx, qy, 0, 1, e)}, // J_q[0][0], J_q[0][1]
+                                {J(qx, qy, 1, 0, e),
+                                 J(qx, qy, 1, 1, e)}}; // J_q[1][0], J_q[1][1]
+
+            double detJ_q = (J_q[0][0] * J_q[1][1]) - (J_q[0][1] * J_q[1][0]);
+
+            double adjJ[2][2] = {{J_q[1][1], -J_q[0][1]},
+                                 {-J_q[1][0], J_q[0][0]}};
+
+            tensor<double, 2> du_dx_q
+               = {(adjJ[0][0] * du_dX_q[0] + adjJ[1][0] * du_dX_q[1]) / detJ_q,
+                  (adjJ[0][1] * du_dX_q[0] + adjJ[1][1] * du_dX_q[1]) / detJ_q};
+
+            double dv_dx_q[2]
+               = {(adjJ[0][0] * dv_dX_q[0] + adjJ[1][0] * dv_dX_q[1]) / detJ_q,
+                  (adjJ[0][1] * dv_dX_q[0] + adjJ[1][1] * dv_dX_q[1]) / detJ_q};
+
+            // compute dF(u, du)/du
+            auto processed_qf_farg_values_u = std::apply(
+               [=](auto &... a) {
+                  return std::make_tuple(derivative_wrt(u_q),
+                                         du_dx_q,
+                                         EvaluateFargValue(a, qx, qy)...);
+               },
+               qf_farg_values);
+
+            auto [f0u, f1u] = std::apply(qf, processed_qf_farg_values_u);
+
+            double f00 = 0.0;
+            if constexpr (std::is_same_v<decltype(f0u), double>)
+            {
+               f00 = 0.0;
+            }
+            else
+            {
+               f00 = f0u.gradient;
+            }
+
+            tensor<double, 2> f10;
+            if constexpr (std::is_same_v<decltype(f1u), tensor<double, 2>>)
+            {
+               f10 = {0.0, 0.0};
+            }
+            else
+            {
+               f10 = {f1u[0].gradient, f1u[1].gradient};
+            }
+
+            // compute dF(u, du)/ddu
+            auto processed_qf_farg_values_du = std::apply(
+               [=](auto &... a) {
+                  return std::make_tuple(u_q,
+                                         derivative_wrt(du_dx_q),
+                                         EvaluateFargValue(a, qx, qy)...);
+               },
+               qf_farg_values);
+
+            auto [f0du, f1du] = std::apply(qf, processed_qf_farg_values_du);
+
+            tensor<double, 2> f01;
+            if constexpr (std::is_same_v<decltype(f0du), double>)
+            {
+               f01 = {0.0, 0.0};
+            }
+            else
+            {
+               f01 = f0du.gradient;
+            }
+            tensor<double, 2, 2> f11 = {
+               {{f1du[0].gradient[0], f1du[0].gradient[1]},
+                {f1du[1].gradient[0], f1du[1].gradient[1]}}};
+
+            double W0 = f00 * v_q + f01[0] * dv_dx_q[0] + f01[1] * dv_dx_q[1];
+
+            double W1[2] = {f10[0] * v_q + f11[0][0] * dv_dx_q[0]
+                               + f11[0][1] * dv_dx_q[1],
+                            f10[1] * v_q + f11[1][0] * dv_dx_q[0]
+                               + +f11[1][1] * dv_dx_q[1]};
+
+            double W0_X = W0 * detJ_q;
+
+            // W1_X = invJ * W1 * detJ
+            //      = adjJ * W1
+            double W1_X[2] = {
+               adjJ[0][0] * W1[0] + adjJ[0][1] * W1[1],
+               adjJ[1][0] * W1[0] + adjJ[1][1] * W1[1],
+            };
+
+            for (int ix = 0; ix < D1D; ix++)
+            {
+               for (int iy = 0; iy < D1D; iy++)
+               {
+                  // @TODO: proper comment
+                  y(ix, iy, e) += (W0_X * v1d(qx, ix) * v1d(qy, iy)
+                                   + W1_X[0] * dv1d_dX(qx, ix) * v1d(qy, iy)
+                                   + W1_X[1] * dv1d_dX(qy, iy) * v1d(qx, ix))
+                                  * W(qx, qy);
+               }
+            }
+         }
+      }
+   }
+}
+
+} // namespace mfem
@@ -0,0 +1,347 @@
+// tensor
+#include <iostream>
+
+#pragma once
+
+template<typename T, int... n>
+struct tensor;
+
+template<typename T>
+struct tensor<T, 1>
+{
+   static constexpr int shape[1] = {1};
+   operator T() { return value; }
+   T value;
+};
+
+template<typename T, int n>
+struct tensor<T, n>
+{
+   static constexpr int shape[1] = {n};
+   constexpr auto &operator[](int i) { return value[i]; };
+   constexpr auto operator[](int i) const { return value[i]; };
+   T value[n];
+};
+
+template<typename T, int first, int... rest>
+struct tensor<T, first, rest...>
+{
+   static constexpr int shape[1 + sizeof...(rest)] = {first, rest...};
+   constexpr auto &operator[](int i) { return value[i]; };
+   constexpr auto operator[](int i) const { return value[i]; };
+   tensor<T, rest...> value[first];
+};
+
+template<int n>
+constexpr int product(int (&values)[n])
+{
+   int p = 1;
+   for (int i = 0; i < n; i++)
+   {
+      p *= values[i];
+   }
+   return p;
+}
+
+template<typename S, typename T, int... n>
+auto operator+(tensor<S, n...> A, tensor<T, n...> B)
+{
+   tensor<decltype(S{} + T{}), n...> C{};
+   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
+   {
+      C[i] = A[i] + B[i];
+   }
+   return C;
+}
+
+template<typename S, typename T, int... n>
+auto operator-(tensor<S, n...> A, tensor<T, n...> B)
+{
+   tensor<decltype(S{} - T{}), n...> C{};
+   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
+   {
+      C[i] = A[i] - B[i];
+   }
+   return C;
+}
+
+template<typename S, typename T, int... n>
+auto operator*(S scale, tensor<T, n...> A)
+{
+   tensor<decltype(S{} * T{}), n...> C{};
+   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
+   {
+      C[i] = scale * A[i];
+   }
+   return C;
+}
+
+template<typename S, typename T, int... n>
+auto operator*(tensor<T, n...> A, S scale)
+{
+   tensor<decltype(T{} * S{}), n...> C{};
+   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
+   {
+      C[i] = A[i] * scale;
+   }
+   return C;
+}
+
+template<typename S, typename T, int... n>
+auto operator/(S scale, tensor<T, n...> A)
+{
+   tensor<decltype(S{} / T{}), n...> C{};
+   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
+   {
+      C[i] = scale / A[i];
+   }
+   return C;
+}
+
+template<typename S, typename T, int... n>
+auto operator/(tensor<T, n...> A, S scale)
+{
+   tensor<decltype(T{} / S{}), n...> C{};
+   for (int i = 0; i < tensor<T, n...>::shape[0]; i++)
+   {
+      C[i] = A[i] / scale;
+   }
+   return C;
+}
+
+template<typename S, typename T, int m, int n, int p>
+auto dot(tensor<S, m, n> A, tensor<T, n, p> B)
+{
+   tensor<decltype(S{} * T{}), m, p> AB{};
+   for (int i = 0; i < m; i++)
+   {
+      for (int j = 0; j < p; j++)
+      {
+         for (int k = 0; k < n; k++)
+         {
+            AB[i][j] = AB[i][j] + A[i][k] * B[k][j];
+         }
+      }
+   }
+   return AB;
+}
+
+template<typename T, int m, int n>
+auto inner(tensor<T, m, n> A, tensor<T, m, n> B)
+{
+   double value = 0.0;
+   for (int i = 0; i < m; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         value += A[i][j] * B[i][j];
+      }
+   }
+   return value;
+}
+
+auto inner(double a, double b)
+{
+   return a * b;
+}
+
+template<typename T, int n>
+auto tr(tensor<T, n, n> A)
+{
+   T trA{};
+   for (int i = 0; i < n; i++)
+   {
+      trA = trA + A[i][i];
+   }
+   return trA;
+}
+
+template<int dim>
+constexpr tensor<double, dim, dim> Identity()
+{
+   tensor<double, dim, dim> I{};
+   for (int i = 0; i < dim; i++)
+   {
+      for (int j = 0; j < dim; j++)
+      {
+         I[i][j] = (i == j);
+      }
+   }
+   return I;
+}
+
+template<typename T, int m, int n>
+auto transpose(const tensor<T, m, n> &A)
+{
+   tensor<T, n, m> AT{};
+   for (int i = 0; i < n; i++)
+   {
+      for (int j = 0; j < m; j++)
+      {
+         AT[i][j] = A[j][i];
+      }
+   }
+   return AT;
+}
+
+template<typename T, int n>
+auto norm(const tensor<T, n> &A)
+{
+   T r = {};
+   for (int i = 0; i < n; i++)
+   {
+      r = r + A[i] * A[i];
+   }
+   return pow(r, 0.5);
+}
+
+template<typename gradient_type>
+struct dual
+{
+   double value;
+   gradient_type gradient;
+};
+
+template<typename gradient_type>
+auto operator+(dual<gradient_type> a, double b)
+{
+   return dual<gradient_type>{a.value + b, a.gradient};
+}
+
+template<typename gradient_type>
+auto operator+(double a, dual<gradient_type> b)
+{
+   return dual<gradient_type>{a + b.value, b.gradient};
+}
+
+template<typename gradient_type>
+auto operator+(dual<gradient_type> a, dual<gradient_type> b)
+{
+   return dual<gradient_type>{a.value + b.value, a.gradient + b.gradient};
+}
+
+template<typename gradient_type>
+auto operator*(dual<gradient_type> a, double b)
+{
+   return dual<gradient_type>{a.value * b, a.gradient * b};
+}
+
+template<typename gradient_type>
+auto operator*(double a, dual<gradient_type> b)
+{
+   return dual<gradient_type>{a * b.value, a * b.gradient};
+}
+
+template<typename gradient_type>
+auto operator*(dual<gradient_type> a, dual<gradient_type> b)
+{
+   return dual<gradient_type>{a.value * b.value,
+                              b.value * a.gradient + a.value * b.gradient};
+}
+
+template<typename gradient_type>
+auto cos(dual<gradient_type> a)
+{
+   return dual<gradient_type>{cos(a.value), -a.gradient * sin(a.value)};
+}
+
+template<typename gradient_type>
+auto exp(dual<gradient_type> a)
+{
+   return dual<gradient_type>{exp(a.value), exp(a.value)};
+}
+
+template<typename gradient_type>
+auto log(dual<gradient_type> a)
+{
+   return dual<gradient_type>{log(a.value), a.gradient / a.value};
+}
+
+template<typename gradient_type>
+auto pow(dual<gradient_type> a, dual<gradient_type> b)
+{
+   double value = pow(a.value, b.value);
+   return dual<gradient_type>{value,
+                              value
+                                 * (a.gradient * (b.value / a.value)
+                                    + b.gradient * log(a.value))};
+}
+
+template<typename gradient_type>
+auto pow(double a, dual<gradient_type> b)
+{
+   double value = pow(a, b.value);
+   return dual<gradient_type>{value, value * b.gradient * log(a)};
+}
+
+template<typename gradient_type>
+auto pow(dual<gradient_type> a, double b)
+{
+   double value = pow(a.value, b);
+   return dual<gradient_type>{value, value * a.gradient * b / a.value};
+}
+
+template<typename T, int... n>
+auto &operator<<(std::ostream &out, dual<T> A)
+{
+   out << '(' << A.value << ' ' << A.gradient << ')';
+   return out;
+}
+
+template<typename T, int... n>
+auto &operator<<(std::ostream &out, tensor<T, n...> A)
+{
+   out << '{' << A[0];
+   for (int i = 1; i < tensor<T, n...>::shape[0]; i++)
+   {
+      out << ", " << A[i];
+   }
+   out << '}';
+   return out;
+}
+
+auto derivative_wrt(double a)
+{
+   return dual<double>{a, 1};
+}
+
+template<typename T, int m>
+auto derivative_wrt(tensor<T, m> A)
+{
+   tensor<dual<tensor<double, m>>, m> A_dual{};
+   for (int i = 0; i < m; i++)
+   {
+      A_dual[i].value = A[i];
+      A_dual[i].gradient[i] = 1.0;
+   }
+   return A_dual;
+}
+
+template<typename T, int m, int n>
+auto derivative_wrt(tensor<T, m, n> A)
+{
+   tensor<dual<tensor<double, m, n>>, m, n> A_dual{};
+   for (int i = 0; i < m; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         A_dual[i][j].value = A[i][j];
+         A_dual[i][j].gradient[i][j] = 1.0;
+      }
+   }
+   return A_dual;
+}
+
+template<typename grad_type, int nrows, int ncols>
+auto directional_derivative(tensor<dual<grad_type>, nrows, ncols> A, grad_type n)
+{
+   tensor<double, nrows, ncols> dA_dn{};
+   for (int i = 0; i < nrows; i++)
+   {
+      for (int j = 0; j < ncols; j++)
+      {
+         dA_dn[i][j] = inner(A[i][j].gradient, n);
+      }
+   }
+   return dA_dn;
+}
@@ -0,0 +1,69 @@
+#include "mfem.hpp"
+#include "tensor.hpp"
+
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
+using namespace mfem;
+
+TEST_CASE("Tensor tests", "[tensor]")
+{
+   SECTION("norm")
+   {
+      tensor<double, 5> a = {{1.0, 2.0, 3.0, 4.0, 5.0}};
+      REQUIRE(norm(a) - sqrt(55) == Approx(0.0));
+   }
+}
+
+TEST_CASE("Dual number tensor tests", "[DualNumber]")
+{
+   auto eps = std::numeric_limits<double>::epsilon();
+   double x = 0.5;
+
+   SECTION("cos")
+   {
+      auto xd = cos(derivative_wrt(x));
+      REQUIRE(abs(-sin(x) - xd.gradient) == Approx(0.0));
+   }
+
+   SECTION("exp")
+   {
+      auto xd = exp(derivative_wrt(x));
+      REQUIRE(abs(exp(x) - xd.gradient) == Approx(0.0));
+   }
+
+   SECTION("log")
+   {
+      auto xd = log(derivative_wrt(x));
+      REQUIRE(abs(1.0 / x - xd.gradient) == Approx(0.0));
+   }
+
+   SECTION("pow")
+   {
+      // f(x) = x^3/2
+      auto xd = pow(derivative_wrt(x), 1.5);
+      REQUIRE(abs(1.5 * pow(x, 0.5) - xd.gradient) == Approx(0.0));
+   }
+
+   SECTION("mixed operations")
+   {
+      auto xd = derivative_wrt(x);
+      auto r = cos(xd) * cos(xd);
+      REQUIRE(abs(-2.0 * sin(x) * cos(x) - r.gradient) == Approx(0.0));
+
+      r = exp(xd) * cos(xd);
+      REQUIRE(abs(exp(x) * (cos(x) - sin(x)) - r.gradient) < eps);
+
+      r = log(xd) * cos(xd);
+      REQUIRE(abs((cos(x) / x - log(x) * sin(x)) - r.gradient) < eps);
+
+      r = exp(xd) * pow(xd, 1.5);
+      REQUIRE(abs((exp(x) * (pow(x, 1.5) + 1.5 * pow(x, 0.5))) - r.gradient)
+              < eps);
+
+      tensor<double, 2> vx = {{0.5, 0.25}};
+      tensor<double, 2> vre = {{0.894427190999916, 0.4472135954999579}};
+      auto vr = norm(derivative_wrt(vx));
+      REQUIRE(norm(vr.gradient - vre) < eps);
+   }
+}
@@ -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
+
@@ -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)")
   {
Author	SHA1	Message	Date
Julian Andrej	7167757d39	renaming and alternative interface	2020-12-03 08:49:51 -08:00
Julian Andrej	60ef17a321	example test	2020-11-30 13:43:45 -08:00
Julian Andrej	97272d60e8	support coordinates in qfunc kernel	2020-11-25 09:59:55 -08:00
Julian Andrej	bb9eb4872b	fixed parallel	2020-11-24 14:29:27 -08:00
Julian Andrej	41798f7028	more tests and example updated	2020-11-24 12:05:45 -08:00
Julian Andrej	633ee29601	more unit tests	2020-11-23 18:10:52 -08:00
Julian Andrej	bb194bbcd5	plaplacian ready	2020-11-23 16:51:43 -08:00
Julian Andrej	9d5977165f	more ad functions	2020-11-23 16:03:57 -08:00
Julian Andrej	68fb8e80a9	add unit tests and work towards p-Laplacian	2020-11-23 14:51:46 -08:00
Julian Andrej	274416caf5	add temporary safeguards for independent variables during AD	2020-11-20 17:53:40 -08:00
Julian Andrej	5853221745	ad working	2020-11-20 16:50:04 -08:00
Sam Mish	c79cccdc08	changes to get variadic args working	2020-11-19 08:53:35 -08:00
Julian Andrej	91f30e88c9	refactor	2020-11-18 15:44:09 -08:00
Julian Andrej	1432ac14a6	working example	2020-11-18 08:26:51 -08:00
Julian Andrej	02cfe4ff0a	parameter pack expansion	2020-11-17 09:00:58 -08:00
Julian Andrej	447214ef57	switch to tuple return type	2020-11-16 10:50:46 -08:00
Julian Andrej	fba6b4d832	gradient matvec	2020-11-13 16:32:48 -08:00
Julian Andrej	12ec6db87c	working operator application	2020-11-11 18:18:38 -08:00
Julian Andrej	f117c0a4e7	first try	2020-11-11 16:04:23 -08:00
Tzanio Kolev	7cd6686e42	Merge pull request #1805 from mfem/navier-filter Implementation of a filter method for Navier	2020-11-03 12:39:54 -08:00
Tzanio	4101b62099	Merge branch 'master' into navier-filter	2020-11-03 12:37:38 -08:00
Julian Andrej	eb0df95218	changelog entry	2020-11-03 12:35:56 -08:00
Veselin Dobrev	a828dbb5fe	Merge pull request #1788 from mfem/mumps-solver Interface with the distributed MUMPS solver	2020-11-03 10:38:51 -08:00
Veselin Dobrev	538ac3a6ab	Merge branch 'master' into mumps-solver	2020-11-03 08:19:30 -08:00
Tzanio Kolev	3c0d346d0d	Merge pull request #964 from mfem/RK_solvers Added three ESDIRK schemes [RK_solvers]	2020-11-03 07:39:13 -08:00
Tzanio	4914e93d0f	Merge branch 'master' into navier-filter	2020-11-03 07:38:47 -08:00
Tzanio	6e63b81c94	Mentioned #964 in CHANGELOG	2020-11-03 07:37:55 -08:00
Tzanio	edbaa5b3fd	Merge branch 'master' into RK_solvers	2020-11-03 07:36:07 -08:00
Tzanio Kolev	eb1cf52ebd	Merge pull request #1835 from mfem/scaled-jacobian-plot-dev Scaled Jacobian plot and locations of bad elements [scaled-jacobian-plot-dev]	2020-11-03 07:31:57 -08:00
Tzanio Kolev	03aa166858	Merge branch 'master' into scaled-jacobian-plot-dev	2020-11-03 07:31:20 -08:00
Tzanio Kolev	37efde9e55	Merge pull request #1855 from mfem/catch-v2.13.2 Upgrade from Catch v2.13.0 to v2.13.2	2020-11-01 13:37:18 -08:00
Tzanio	a17a1f1a2d	Updated CHANGELOG	2020-11-01 13:32:02 -08:00
Tzanio	0640727b77	Merge branch 'master' into mumps-solver	2020-11-01 13:31:21 -08:00
Tzanio	01c0f4af68	Updated CHANGELOG	2020-11-01 13:29:51 -08:00
Tzanio	8de1d37871	Merge branch 'master' into scaled-jacobian-plot-dev	2020-11-01 13:27:11 -08:00
Tzanio	0359c410a9	Merge branch 'master' into RK_solvers	2020-11-01 13:26:02 -08:00
Tzanio	d3028ea057	Merge branch 'master' into navier-filter	2020-11-01 13:24:39 -08:00
Tzanio	9f8e2f8bc7	Updated CHANGELOG	2020-11-01 13:22:53 -08:00
Tzanio	2436fa67ca	Merge branch 'master' into catch-v2.13.2	2020-11-01 13:20:54 -08:00
Tzanio Kolev	b4ba111505	Merge pull request #1857 from mfem/new-dev-version-4.2.1 Update version numbers to 4.2.1 -- a new development version [new-dev-version-4.2.1]	2020-11-01 13:20:07 -08:00
Veselin Dobrev	7bf942176e	Merge branch 'master' into mumps-solver	2020-11-01 03:10:08 -08:00
Veselin Dobrev	4149fbc0d6	Update version numbers to 4.2.1 -- a new development version.	2020-11-01 02:18:15 -08:00
Will Pazner	949d43f80e	Upgrade from Catch v2.13.0 to v2.13.2 Adds a workaround for internal NVCC compiler error: nvcc segfaults when encountering a call to __builtin_constant_p with an argument that involves a new-expression trying to invoke a user-provided constructor We would encounter this problem with complex REQUIRE statements (e.g. in test_ode.cpp).	2020-10-30 09:42:10 -07:00
Tzanio	6be78c57e1	Added comments and a refinement for the scaled Jacobian option in mesh explorer.	2020-10-24 14:46:22 -07:00
Walter I. Nissen	ca37dc2e67	Fixed brace style	2020-10-23 15:25:49 -07:00
Walter I. Nissen	eb7833533f	Slight comments	2020-10-23 15:04:39 -07:00
Walter I. Nissen	9d0a7ac938	Tweak to output to use mfem-style Print function and summary of elements not printed.	2020-10-23 15:04:39 -07:00
Walter I. Nissen	f34b48335a	Add scaled jacobians with the 'j' key	2020-10-23 15:04:39 -07:00
Walter I. Nissen	4e570ad11c	Thanks to Tzanio for telling me where to put this.	2020-10-23 15:04:39 -07:00
psocratis	4b09bc532f	moving macros to .cpp file	2020-10-22 12:06:11 -07:00
Veselin Dobrev	3ab4dd4194	Some formatting edits.	2020-10-22 11:54:40 -07:00
psocratis	ca3f53b37b	fixed version in changelog	2020-10-22 11:04:10 -07:00
psocratis	799069babc	renaming info and icntl to avoid conflicts with catch. INFO->MUMPS_INFO, ICNTL->MUMPS_ICNTL	2020-10-21 16:18:01 -07:00
psocratis	a14b018692	adding check for HYPRE_BIGINT	2020-10-21 13:49:42 -07:00
Julian Andrej	c10ede8fba	clarify docs	2020-10-19 16:33:42 -07:00
Julian Andrej	81dcca3e12	reference formatting	2020-10-19 11:28:45 -07:00
Julian Andrej	bdd237c7ea	spelling	2020-10-19 11:27:41 -07:00
Julian Andrej	cacc93d8dc	mention filter algorithm reference	2020-10-19 11:27:00 -07:00
psocratis	3b9c4b28ad	style	2020-10-18 15:43:19 -07:00
psocratis	3df5658695	fixing comments for the set* methods	2020-10-18 15:42:40 -07:00
psocratis	8cc5feb057	resolving conflicts with master. Minor changes addressing reviewers' comments	2020-10-18 15:33:05 -07:00
psocratis	84d4013a50	fix conflicts	2020-10-12 18:09:01 -07:00
Julian Andrej	c507d3dc3d	make style	2020-10-12 13:43:33 -07:00
Ben Southworth	62d78501a3	Change to three space=indent	2020-10-12 13:25:37 -07:00
Ben Southworth	5fbe7bd995	Added unit tests	2020-10-12 13:01:00 -07:00
bensworth	e4a8344e94	Merge branch 'master' into RK_solvers	2020-10-09 09:41:31 -06:00
Socratis	cb5f1830de	checked cmake both on linux and macos	2020-10-08 16:54:57 -07:00
psocratis	d91e31458c	minor fix in cmake	2020-10-08 11:31:27 -07:00
Julian Andrej	7893c819a9	implementation of a filter method	2020-10-07 16:03:33 -07:00
Julian Andrej	561daa66b6	add dependencies to mumps	2020-10-02 16:52:07 -07:00
Julian Andrej	cd1cacd1a6	correct findmumps cmake	2020-10-02 16:48:33 -07:00
psocratis	2b94b8bd43	configuring cmake	2020-10-02 14:56:38 -07:00
Socratis	b6e4244806	minor tweaks	2020-10-02 11:32:40 -07:00
Socratis	59eaee6f13	fixing travis doc errors	2020-10-01 19:31:03 -07:00
Socratis	9fd56d22a2	fixed newly introduced valgrind mem leak	2020-10-01 19:11:19 -07:00
psocratis	fc3babf85d	minor, initialization of int arrays	2020-10-01 19:03:35 -07:00
Socratis	1d0b509827	mumps_solver valgrind issue free. Added mumps_solver to ex25p	2020-10-01 18:18:58 -07:00
psocratis	50b525907f	minor	2020-09-30 19:11:28 -07:00
psocratis	5137ef5604	clean up	2020-09-29 19:42:44 -07:00
psocratis	8f1d08dc93	Adding symmetric and SPD matrix case	2020-09-29 19:04:17 -07:00
Julian Andrej	e63f9887f4	add MUMPS version check to makefile build system	2020-09-29 14:13:14 -07:00
psocratis	2819c0e361	more cleaning up	2020-09-29 13:23:43 -07:00
psocratis	30b0ae9581	cleanup	2020-09-28 17:47:08 -07:00
psocratis	59a32942f2	make style	2020-09-28 15:06:23 -07:00
psocratis	90af7527d0	mumps interface works with any combination of distributed RHS and Solution or not... still needs some cleanup and a valgrind run	2020-09-28 15:05:18 -07:00
psocratis	ac6ff275ba	Adding MUMPSSolver interface. Finished the case where RHS and Solution are gathered on host	2020-09-27 18:09:25 -07:00
Julian Andrej	b78402803f	bugs everywhere	2020-09-26 18:21:36 -07:00
Julian Andrej	6e5bb7e5e1	more distributed solution vector fiddling	2020-09-26 16:25:48 -07:00
Julian Andrej	99a3803b4f	distributed rhs and solution	2020-09-26 15:48:07 -07:00
Julian Andrej	c0495229fc	add MUMPS to the makefile buildsystem	2020-09-26 14:08:56 -07:00
psocratis	adce9a5755	Adding interface for distributed MUMPS still needs implementation for distributed RHS	2020-09-25 20:22:36 -07:00
Ben Southworth	bb28ac7435	Added three ESDIRK schemes	2019-06-18 07:49:15 -07:00