Fixing copyright year

Removing partial assembly code from example
Updating other DiffusionIntegrator methods
2021-04-09 15:49:41 -07:00 · 2021-04-09 15:49:26 -07:00 · 2021-04-09 15:48:59 -07:00 · 2021-03-28 12:32:35 -04:00 · 2021-03-28 12:31:36 -04:00 · 2021-03-26 17:03:46 -07:00
45 changed files with 723 additions and 2495 deletions
@@ -8,24 +8,6 @@
                               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 17)
+set(CMAKE_CXX_STANDARD 11)
 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.1)
+set(${PROJECT_NAME}_VERSION 4.2.0)

 # 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 MUMPS PETSC SLEPC STRUMPACK PUMI)
+  set(PKGS_NEED_MPI SUPERLU PETSC SLEPC STRUMPACK PUMI)
  foreach(PKG IN LISTS PKGS_NEED_MPI)
    if (MFEM_USE_${PKG})
      message(STATUS "Disabling package ${PKG} - requires MPI")
@@ -264,15 +264,6 @@ 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)
@@ -391,7 +382,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 MUMPS STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
+    SLEPC MESQUITE SuperLUDist 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,12 +348,6 @@ 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
@@ -584,12 +578,6 @@ 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
@@ -849,7 +837,6 @@ 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,7 +32,6 @@ 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,9 +86,6 @@
 // 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

@@ -1,23 +0,0 @@
-# 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,9 +42,6 @@
 #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,10 +95,6 @@
 // #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,7 +33,6 @@ 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,7 +34,6 @@ 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)
@@ -119,15 +118,6 @@ 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,7 +123,6 @@ 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
@@ -172,7 +171,7 @@ HYPRE_OPT = -I$(HYPRE_DIR)/include
 HYPRE_LIB = -L$(HYPRE_DIR)/lib -lHYPRE

 # METIS library configuration
-ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
+ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK),NONO)
   ifeq ($(MFEM_USE_METIS_5),NO)
     METIS_DIR = @MFEM_DIR@/../metis-4.0
     METIS_OPT =
@@ -253,12 +252,12 @@ SCOTCH_OPT = -I$(SCOTCH_DIR)/include
 SCOTCH_LIB = -L$(SCOTCH_DIR)/lib -lptscotch -lptscotcherr -lscotch -lscotcherr\
 -lpthread

-# SCALAPACK library configuration (required by STRUMPACK and MUMPS)
+# SCALAPACK library configuration (required by STRUMPACK)
 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 or MUMPS
+# MPI Fortran library, needed e.g. by STRUMPACK
 # MPICH:
 MPI_FORTRAN_LIB = -lmpifort
 # OpenMPI:
@@ -266,12 +265,6 @@ 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)
@@ -1,33 +0,0 @@
-// 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,10 +42,6 @@ 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 ...)
@@ -54,13 +50,6 @@ $(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)\
@@ -73,18 +62,10 @@ check-smx:
 	$(call mfem-info, MFEM_USE_SIMMETRIX = $(MFEM_USE_SIMMETRIX))
 	$(eval export MFEM_USE_SIMMETRIX)

-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,)
+header: $(if $(MPI),get-hypre-version,) $(if $(SMX),check-smx)
 	$(call mfem-info, Writing $(CONFIG_HPP) ...)
 	@set -- && \
-	for def in $${MFEM_DEFINES} $(if $(MPI),MFEM_HYPRE_VERSION) $(SMX) \
-	           $(if $(MUMPS),MFEM_MUMPS_VERSION); do \
+	for def in $${MFEM_DEFINES} $(if $(MPI),MFEM_HYPRE_VERSION) $(SMX); 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.1
+PROJECT_NUMBER         = v4.2.0

 # 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,8 +164,6 @@ 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);
@@ -187,26 +185,10 @@ 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;
@@ -448,7 +430,6 @@ 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();
@@ -461,19 +442,7 @@ int main(int argc, char *argv[])
      superlu.Mult(B, X);
      delete A;
   }
-#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
+#else
   // 16a. Set up the parallel Bilinear form a(.,.) for the preconditioner
   //
   //    In Comp
@@ -481,7 +450,6 @@ 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);
@@ -531,6 +499,7 @@ 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.
@@ -0,0 +1,322 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+Mesh * GetMesh(int type);
+
+void trans(const Vector &x, Vector &r);
+
+void sigmaFunc(const Vector &x, DenseMatrix &s);
+
+double uExact(const Vector &x)
+{
+   return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int order = 3;
+   int mesh_type = 4; // Default to Quadrilateral mesh
+   int ref_levels = 0;
+   bool static_cond = false;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_type, "-mt", "--mesh-type",
+                  "Mesh type: 3 - Triangular, 4 - Quadrilateral.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = GetMesh(mesh_type);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   mesh->SetCurvature(3);
+   mesh->Transform(trans);
+
+   // 4. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   bool delete_fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+      delete_fec = false;
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
+   }
+   FiniteElementSpace fespace(mesh, fec);
+   cout << "Number of finite element unknowns: "
+        << fespace.GetTrueVSize() << endl;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking all
+   //    the boundary attributes from the mesh as essential (Dirichlet) and
+   //    converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   LinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 7. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //    domain integrator.
+   BilinearForm a(&fespace);
+   MatrixFunctionCoefficient sigma(3, sigmaFunc);
+   a.AddDomainIntegrator(new DiffusionIntegrator(sigma));
+
+   // 9. Assemble the bilinear form and the corresponding linear system,
+   //    applying any necessary transformations such as: eliminating boundary
+   //    conditions, applying conforming constraints for non-conforming AMR,
+   //    static condensation, etc.
+   if (static_cond) { a.EnableStaticCondensation(); }
+   a.Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   cout << "Size of linear system: " << A->Height() << endl;
+
+   // 10. Solve the linear system A X = B.
+   if (!pa)
+   {
+      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+   }
+   else // Jacobi preconditioning in partial assembly mode
+   {
+      if (UsesTensorBasis(fespace))
+      {
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
+      }
+      else
+      {
+         CG(*A, B, X, 1, 400, 1e-12, 0.0);
+      }
+   }
+
+   // 11. Recover the solution as a finite element grid function.
+   a.RecoverFEMSolution(X, b, x);
+
+   FunctionCoefficient uCoef(uExact);
+   double err = x.ComputeL2Error(uCoef);
+
+   mfem::out << "|u - u_h|_2 = " << err << endl;
+
+   // 12. Save the refined mesh and the solution. This output can be viewed later
+   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
+   ofstream mesh_ofs("refined.mesh");
+   mesh_ofs.precision(8);
+   mesh->Print(mesh_ofs);
+   ofstream sol_ofs("sol.gf");
+   sol_ofs.precision(8);
+   x.Save(sol_ofs);
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 14. Free the used memory.
+   if (delete_fec)
+   {
+      delete fec;
+   }
+   delete mesh;
+
+   return 0;
+}
+
+Mesh * GetMesh(int type)
+{
+   Mesh * mesh = NULL;
+
+   if (type == 3)
+   {
+      mesh = new Mesh(2, 12, 16, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+      mesh->AddVertex( 0.0, -1.0, 0.5);
+      mesh->AddVertex( 1.0,  0.0, 0.5);
+      mesh->AddVertex( 0.0,  1.0, 0.5);
+      mesh->AddVertex(-1.0,  0.0, 0.5);
+
+      mesh->AddTriangle(0, 1, 8);
+      mesh->AddTriangle(1, 5, 8);
+      mesh->AddTriangle(5, 4, 8);
+      mesh->AddTriangle(4, 0, 8);
+      mesh->AddTriangle(1, 2, 9);
+      mesh->AddTriangle(2, 6, 9);
+      mesh->AddTriangle(6, 5, 9);
+      mesh->AddTriangle(5, 1, 9);
+      mesh->AddTriangle(2, 3, 10);
+      mesh->AddTriangle(3, 7, 10);
+      mesh->AddTriangle(7, 6, 10);
+      mesh->AddTriangle(6, 2, 10);
+      mesh->AddTriangle(3, 0, 11);
+      mesh->AddTriangle(0, 4, 11);
+      mesh->AddTriangle(4, 7, 11);
+      mesh->AddTriangle(7, 3, 11);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else if (type == 4)
+   {
+      mesh = new Mesh(2, 8, 4, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+
+      mesh->AddQuad(0, 1, 5, 4);
+      mesh->AddQuad(1, 2, 6, 5);
+      mesh->AddQuad(2, 3, 7, 6);
+      mesh->AddQuad(3, 0, 4, 7);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else
+   {
+      MFEM_ABORT("Unrecognized mesh type " << type << "!");
+   }
+   mesh->FinalizeTopology();
+
+   return mesh;
+}
+
+void trans(const Vector &x, Vector &r)
+{
+   r.SetSize(3);
+
+   double tol = 1e-6;
+   double theta = 0.0;
+   if (fabs(x[1] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (x[0] - 2.0);
+   }
+   else if (fabs(x[0] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * x[1];
+   }
+   else if (fabs(x[1] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (2.0 - x[0]);
+   }
+   else if (fabs(x[0] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (4.0 - x[1]);
+   }
+   else
+   {
+      cout << "side not recognized "
+           << x[0] << " " << x[1] << " " << x[2] << endl;
+   }
+
+   r[0] = cos(theta);
+   r[1] = sin(theta);
+   r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
+}
+
+void sigmaFunc(const Vector &x, DenseMatrix &s)
+{
+   s.SetSize(3);
+   double a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
+   s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
+   s(0,1) = x[0] * x[1] * (8.0 / a - 0.5);
+   s(0,2) = 0.0;
+   s(1,0) = s(0,1);
+   s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] / a;
+   s(1,2) = 0.0;
+   s(2,0) = 0.0;
+   s(2,1) = 0.0;
+   s(2,2) = a / 32.0;
+}
@@ -552,13 +552,28 @@ void DiffusionIntegrator::AssembleElementMatrix
   bool square = (dim == spaceDim);
   double w;

+   if (VQ)
+   {
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
+   }
+
 #ifdef MFEM_THREAD_SAFE
-   DenseMatrix dshape(nd,dim), dshapedxt(nd,spaceDim), invdfdx(dim,spaceDim);
+   DenseMatrix dshape(nd, dim), dshapedxt(nd, spaceDim);
+   DenseMatrix dshapedxt_m(nd, MQ ? spaceDim : 0);
   Vector D(VQ ? VQ->GetVDim() : 0);
 #else
-   dshape.SetSize(nd,dim);
-   dshapedxt.SetSize(nd,spaceDim);
-   invdfdx.SetSize(dim,spaceDim);
+   dshape.SetSize(nd, dim);
+   dshapedxt.SetSize(nd, spaceDim);
+   dshapedxt_m.SetSize(nd, MQ ? spaceDim : 0);
+   M.SetSize(MQ ? spaceDim : 0);
   D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
   elmat.SetSize(nd);
@@ -579,10 +594,10 @@ void DiffusionIntegrator::AssembleElementMatrix
      Mult(dshape, Trans.AdjugateJacobian(), dshapedxt);
      if (MQ)
      {
-         MQ->Eval(invdfdx, Trans, ip);
-         invdfdx *= w;
-         Mult(dshapedxt, invdfdx, dshape);
-         AddMultABt(dshape, dshapedxt, elmat);
+         MQ->Eval(M, Trans, ip);
+         M *= w;
+         Mult(dshapedxt, M, dshapedxt_m);
+         AddMultABt(dshapedxt_m, dshapedxt, elmat);
      }
      else if (VQ)
      {
@@ -612,10 +627,25 @@ void DiffusionIntegrator::AssembleElementMatrix2(
   bool square = (dim == spaceDim);
   double w;

+   if (VQ)
+   {
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
+   }
+
 #ifdef MFEM_THREAD_SAFE
   DenseMatrix dshape(tr_nd, dim), dshapedxt(tr_nd, spaceDim);
   DenseMatrix te_dshape(te_nd, dim), te_dshapedxt(te_nd, spaceDim);
   DenseMatrix invdfdx(dim, spaceDim);
+   DenseMatrix dshapedxt_m(te_nd, MQ ? spaceDim : 0);
+   DenseMatrix M(MQ ? spaceDim : 0);
   Vector D(VQ ? VQ->GetVDim() : 0);
 #else
   dshape.SetSize(tr_nd, dim);
@@ -623,6 +653,8 @@ void DiffusionIntegrator::AssembleElementMatrix2(
   te_dshape.SetSize(te_nd, dim);
   te_dshapedxt.SetSize(te_nd, spaceDim);
   invdfdx.SetSize(dim, spaceDim);
+   dshapedxt_m.SetSize(te_nd, MQ ? spaceDim : 0);
+   M.SetSize(MQ ? spaceDim : 0);
   D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
   elmat.SetSize(te_nd, tr_nd);
@@ -645,10 +677,10 @@ void DiffusionIntegrator::AssembleElementMatrix2(
      // invdfdx, dshape, and te_dshape no longer needed
      if (MQ)
      {
-         MQ->Eval(invdfdx, Trans, ip);
-         invdfdx *= w;
-         Mult(te_dshapedxt, invdfdx, te_dshape);
-         AddMultABt(te_dshape, dshapedxt, elmat);
+         MQ->Eval(M, Trans, ip);
+         M *= w;
+         Mult(te_dshapedxt, M, dshapedxt_m);
+         AddMultABt(dshapedxt_m, dshapedxt, elmat);
      }
      else if (VQ)
      {
@@ -674,24 +706,34 @@ void DiffusionIntegrator::AssembleElementVector(
 {
   int nd = el.GetDof();
   int dim = el.GetDim();
+   int spaceDim = Tr.GetSpaceDim();
   double w;

   if (VQ)
   {
-      MFEM_VERIFY(VQ->GetVDim() == dim, "Unexpected dimension for VectorCoefficient");
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
   }

 #ifdef MFEM_THREAD_SAFE
-   DenseMatrix dshape(nd,dim), invdfdx(dim), mq(dim);
+   DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim), M(MQ ? spaceDim : 0);
   Vector D(VQ ? VQ->GetVDim() : 0);
 #else
   dshape.SetSize(nd,dim);
-   invdfdx.SetSize(dim);
-   mq.SetSize(dim);
+   invdfdx.SetSize(dim, spaceDim);
+   M.SetSize(MQ ? spaceDim : 0);
   D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
   vec.SetSize(dim);
-   pointflux.SetSize(dim);
+   vecdxt.SetSize((VQ || MQ) ? spaceDim : 0);
+   pointflux.SetSize(spaceDim);

   elvect.SetSize(nd);

@@ -718,19 +760,19 @@ void DiffusionIntegrator::AssembleElementVector(
      }
      else
      {
-         dshape.MultTranspose(elfun, pointflux);
-         invdfdx.MultTranspose(pointflux, vec);
+         dshape.MultTranspose(elfun, vec);
+         invdfdx.MultTranspose(vec, vecdxt);
         if (MQ)
         {
-            MQ->Eval(mq, Tr, ip);
-            mq.Mult(vec, pointflux);
+            MQ->Eval(M, Tr, ip);
+            M.Mult(vecdxt, pointflux);
         }
         else
         {
            VQ->Eval(D, Tr, ip);
-            for (int j=0; j<dim; ++j)
+            for (int j=0; j<spaceDim; ++j)
            {
-               pointflux[j] *= D[j];
+               pointflux[j] = D[j] * vecdxt[j];
            }
         }
      }
@@ -750,14 +792,32 @@ void DiffusionIntegrator::ComputeElementFlux
   dim = el.GetDim();
   spaceDim = Trans.GetSpaceDim();

+   if (VQ)
+   {
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
+   }
+
 #ifdef MFEM_THREAD_SAFE
   DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim);
+   DenseMatrix M(MQ ? spaceDim : 0);
+   Vector D(VQ ? VQ->GetVDim() : 0);
 #else
   dshape.SetSize(nd,dim);
   invdfdx.SetSize(dim, spaceDim);
+   M.SetSize(MQ ? spaceDim : 0);
+   D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
   vec.SetSize(dim);
-   pointflux.SetSize(spaceDim);
+   vecdxt.SetSize(spaceDim);
+   pointflux.SetSize(MQ ? spaceDim : 0);

   const IntegrationRule &ir = fluxelem.GetNodes();
   fnd = ir.GetNPoints();
@@ -771,28 +831,38 @@ void DiffusionIntegrator::ComputeElementFlux

      Trans.SetIntPoint (&ip);
      CalcInverse(Trans.Jacobian(), invdfdx);
-      invdfdx.MultTranspose(vec, pointflux);
+      invdfdx.MultTranspose(vec, vecdxt);

-      if (!MQ)
+      if (!MQ && !VQ)
      {
         if (Q && with_coef)
         {
-            pointflux *= Q->Eval(Trans,ip);
+            vecdxt *= Q->Eval(Trans,ip);
         }
         for (j = 0; j < spaceDim; j++)
         {
-            flux(fnd*j+i) = pointflux(j);
+            flux(fnd*j+i) = vecdxt(j);
         }
      }
      else
      {
-         // assuming dim == spaceDim
-         MFEM_ASSERT(dim == spaceDim, "TODO");
-         MQ->Eval(invdfdx, Trans, ip);
-         invdfdx.Mult(pointflux, vec);
+         if (MQ)
+         {
+            MQ->Eval(M, Trans, ip);
+            M.Mult(vecdxt, pointflux);
+         }
+         else
+         {
+            VQ->Eval(D, Trans, ip);
+            for (int j=0; j<spaceDim; ++j)
+            {
+               pointflux[j] = D[j] * vecdxt[j];
+            }
+
+         }
         for (j = 0; j < dim; j++)
         {
-            flux(fnd*j+i) = vec(j);
+            flux(fnd*j+i) = pointflux(j);
         }
      }
   }
@@ -807,13 +877,13 @@ double DiffusionIntegrator::ComputeFluxEnergy
   int spaceDim = Trans.GetSpaceDim();

 #ifdef MFEM_THREAD_SAFE
-   DenseMatrix mq;
+   DenseMatrix M;
 #endif

   shape.SetSize(nd);
   pointflux.SetSize(spaceDim);
-   if (d_energy) { vec.SetSize(dim); }
-   if (MQ) { mq.SetSize(dim); }
+   if (d_energy) { vec.SetSize(spaceDim); }
+   if (MQ) { M.SetSize(spaceDim); }

   int order = 2 * fluxelem.GetOrder(); // <--
   const IntegrationRule *ir = &IntRules.Get(fluxelem.GetGeomType(), order);
@@ -846,8 +916,9 @@ double DiffusionIntegrator::ComputeFluxEnergy
      }
      else
      {
-         MQ->Eval(mq, Trans, ip);
-         energy += w * mq.InnerProduct(pointflux, pointflux);
+         MFEM_ASSERT(dim == spaceDim, "TODO");
+         MQ->Eval(M, Trans, ip);
+         energy += w * M.InnerProduct(pointflux, pointflux);
      }

      if (d_energy)
@@ -1902,9 +1902,9 @@ protected:
   MatrixCoefficient *MQ;

 private:
-   Vector vec, pointflux, shape;
+   Vector vec, vecdxt, pointflux, shape;
 #ifndef MFEM_THREAD_SAFE
-   DenseMatrix dshape, dshapedxt, invdfdx, mq;
+   DenseMatrix dshape, dshapedxt, invdfdx, M, dshapedxt_m;
   DenseMatrix te_dshape, te_dshapedxt;
   Vector D;
 #endif
@@ -656,6 +656,7 @@ double TMOP_Metric_315::EvalW(const DenseMatrix &Jpt) const
   ie.SetJacobian(Jpt.GetData());
   const double c1 = ie.Get_I3b() - 1.0;
   return c1*c1;
+   //return c1*c1*c1*c1*c1*c1;
 }

 void TMOP_Metric_315::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
@@ -664,6 +665,7 @@ void TMOP_Metric_315::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
   // P = 2*(I3b - 1)*dI3b
   ie.SetJacobian(Jpt.GetData());
   P.Set(2*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
+   //P.Set(6*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
 }

 void TMOP_Metric_315::AssembleH(const DenseMatrix &Jpt,
@@ -677,6 +679,9 @@ void TMOP_Metric_315::AssembleH(const DenseMatrix &Jpt,
   ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
   ie.Assemble_TProd(2*weight, ie.Get_dI3b(), A.GetData());
   ie.Assemble_ddI3b(2*weight*(ie.Get_I3b() - 1.0), A.GetData());
+   //ie.Assemble_TProd(6*weight, ie.Get_dI3b(), A.GetData());
+   //ie.Assemble_ddI3b(6*weight*(ie.Get_I3b() - 1.0), A.GetData());
+   
 }

 double TMOP_Metric_316::EvalW(const DenseMatrix &Jpt) const
@@ -197,10 +197,19 @@ void SerialAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
   di_dt = 0.0;
   CGSolver lin_solver;
   DSmoother prec;
+   /*
+   FGMRESSolver lin_solver;
+   GMRESSolver prec;
+   prec.SetMaxIter(50);
+   prec.SetRelTol(0.0);
+   prec.SetAbsTol(0.0);    
+   prec.SetOperator(M.SpMat());
+   */
+   
   lin_solver.SetPreconditioner(prec);
   lin_solver.SetOperator(M.SpMat());
   lin_solver.SetRelTol(1e-12); lin_solver.SetAbsTol(0.0);
-   lin_solver.SetMaxIter(100);
+   lin_solver.SetMaxIter(200);
   lin_solver.SetPrintLevel(0);
   lin_solver.Mult(rhs, di_dt);
 }
@@ -83,11 +83,6 @@ 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)
@@ -2497,13 +2497,13 @@ void MultADBt(const DenseMatrix &A, const Vector &D,

 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
 {
-#ifdef MFEM_DEBUG
+  #ifdef MFEM_DEBUG
   if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
       A.Width() != B.Width())
   {
      mfem_error("AddMultABt(...): dimension mismatch");
   }
-#endif
+  #endif

 #ifdef MFEM_USE_LAPACK
   static char transa = 'N', transb = 'T';
@@ -49,10 +49,6 @@
 #include "hypre_parcsr.hpp"
 #include "hypre.hpp"

-#ifdef MFEM_USE_MUMPS
-#include "mumps.hpp"
-#endif
-
 #ifdef MFEM_USE_PETSC
 #include "petsc.hpp"
 #endif
@@ -1,422 +0,0 @@
-// 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
@@ -1,156 +0,0 @@
-// 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.0-a-b)*dt, k);
+   x.Add((1.-a-b)*dt, k);

   f->SetTime(t + dt);
   f->ImplicitSolve(a*dt, x, k);
@@ -685,101 +685,6 @@ 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,48 +460,6 @@ 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 = 40201
+MFEM_VERSION = 40200
 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 MUMPS STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
+   PKGS_NEED_MPI = SUPERLU 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 MUMPS METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
+MFEM_REQ_LIB_DEPS = SUPERLU 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,8 +332,7 @@ 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_USE_MUMPS MFEM_SOURCE_DIR\
- MFEM_INSTALL_DIR
+ MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX 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\
@@ -367,7 +366,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 MUMPS_OPT
+   export VERBOSE HYPRE_OPT PUMI_DIR
 endif

 # If we have 'install' target, export variables used by config/makefile
@@ -635,7 +634,6 @@ 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,4 +26,3 @@ add_subdirectory(tools)
 add_subdirectory(toys)
 add_subdirectory(nurbs)
 add_subdirectory(gslib)
-add_subdirectory(variationalform)
@@ -351,7 +351,6 @@ 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"
@@ -559,8 +558,6 @@ 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);
@@ -592,23 +589,10 @@ 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)
         {
@@ -708,7 +692,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 == 'J' || mk == 'p')
+          mk == 'k' || mk == 'p')
      {
         Array<int> bdr_part;
         Array<int> part(mesh->GetNE());
@@ -781,7 +765,7 @@ int main (int argc, char *argv[])
            h_max = -h_min;
            for (int i = 0; i < mesh->GetNE(); i++)
            {
-               Geometry::Type geom = mesh->GetElementBaseGeometry(i);
+               int geom = mesh->GetElementBaseGeometry(i);
               ElementTransformation *T = mesh->GetElementTransformation(i);
               T->SetIntPoint(&Geometries.GetCenter(geom));
               Geometries.JacToPerfJac(geom, T->Jacobian(), J);
@@ -806,7 +790,7 @@ int main (int argc, char *argv[])
            DenseMatrix J(dim);
            for (int i = 0; i < mesh->GetNE(); i++)
            {
-               Geometry::Type geom = mesh->GetElementBaseGeometry(i);
+               int geom = mesh->GetElementBaseGeometry(i);
               ElementTransformation *T = mesh->GetElementTransformation(i);
               T->SetIntPoint(&Geometries.GetCenter(geom));
               Geometries.JacToPerfJac(geom, T->Jacobian(), J);
@@ -814,48 +798,6 @@ 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;
@@ -1008,7 +950,7 @@ int main (int argc, char *argv[])
            else
            {
               sol_sock << "fem3d_gf_data_keys\n";
-               if (mk == 'v' || mk == 'h' || mk == 'k' || mk == 'J')
+               if (mk == 'v' || mk == 'h' || mk == 'k')
               {
                  mesh->Print(sol_sock);
               }
@@ -311,7 +311,7 @@ int main(int argc, char *argv[])
   x0 = x;

   // 11. Form the integrator that uses the chosen metric and target.
-   double tauval = -0.1;
+   double tauval = -0.002;
   TMOP_QualityMetric *metric = NULL;
   switch (metric_id)
   {
@@ -290,21 +290,6 @@ 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();
 }

@@ -533,18 +518,6 @@ 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)
@@ -1102,6 +1075,4 @@ NavierSolver::~NavierSolver()
   delete pfec;
   delete vfes;
   delete pfes;
-   delete vfec_filter;
-   delete vfes_filter;
 }
@@ -241,19 +241,6 @@ 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();
@@ -425,14 +412,6 @@ 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
@@ -1,15 +0,0 @@
-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)
@@ -1,23 +0,0 @@
-#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
@@ -1,94 +0,0 @@
-#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
@@ -1,109 +0,0 @@
-#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
@@ -1,164 +0,0 @@
-#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;
-}
@@ -1,118 +0,0 @@
-#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;
-}
@@ -1,404 +0,0 @@
-#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
@@ -1,347 +0,0 @@
-// 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;
-}
@@ -1,69 +0,0 @@
-#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.2
- *  Generated: 2020-10-07 11:32:53.302017
+ *  Catch v2.13.0
+ *  Generated: 2020-07-12 20:07:49.015950
 *  ----------------------------------------------------------
 *  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 2
+#define CATCH_VERSION_PATCH 0

 #ifdef __clang__
 #    pragma clang system_header
@@ -132,9 +132,13 @@ 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) && !defined(__CUDACC__)
+#if defined(__GNUC__) && !defined(__clang__) && !defined(__ICC)
 #    define CATCH_INTERNAL_START_WARNINGS_SUPPRESSION _Pragma( "GCC diagnostic push" )
 #    define CATCH_INTERNAL_STOP_WARNINGS_SUPPRESSION  _Pragma( "GCC diagnostic pop" )

@@ -158,7 +162,7 @@ namespace Catch {
 // ```
 //
 // Therefore, `CATCH_INTERNAL_IGNORE_BUT_WARN` is not implemented.
-#  if !defined(__ibmxl__) && !defined(__CUDACC__)
+#  if !defined(__ibmxl__)
 #    define CATCH_INTERNAL_IGNORE_BUT_WARN(...) (void)__builtin_constant_p(__VA_ARGS__) /* NOLINT(cppcoreguidelines-pro-type-vararg, hicpp-vararg) */
 #  endif

@@ -240,6 +244,10 @@ 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)
@@ -322,10 +330,7 @@ namespace Catch {

  // Check if byte is available and usable
  #  if __has_include(<cstddef>) && defined(CATCH_CPP17_OR_GREATER)
-  #    include <cstddef>
-  #    if __cpp_lib_byte > 0
-  #      define CATCH_INTERNAL_CONFIG_CPP17_BYTE
-  #    endif
+  #    define CATCH_INTERNAL_CONFIG_CPP17_BYTE
  #  endif // __has_include(<cstddef>) && defined(CATCH_CPP17_OR_GREATER)

  // Check if variant is available and usable
@@ -368,6 +373,10 @@ 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
@@ -1096,7 +1105,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++)... };/* NOLINT */ \
+                    (void)expander{(reg_test(Types{}, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index]), Tags } ), index++, 0)... };/* NOLINT */ \
                }\
            };\
            static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1142,7 +1151,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++)... };/* 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++, 0)... };/* NOLINT */\
                }                                                     \
            };                                                        \
            static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){ \
@@ -1186,7 +1195,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++)... };/* 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++, 0)... };/* NOLINT */\
            }                                                     \
        };\
        static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){ \
@@ -1220,7 +1229,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++)... };/* NOLINT */ \
+                    (void)expander{(reg_test(Types{}, #ClassName, Catch::NameAndTags{ Name " - " + std::string(tmpl_types[index]), Tags } ), index++, 0)... };/* NOLINT */ \
                }\
            };\
            static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1269,7 +1278,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++)... };/* 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++, 0)... };/* NOLINT */ \
                }\
            };\
            static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1316,7 +1325,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++)... };/* 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++, 0)... };/* NOLINT */ \
                }\
            };\
            static int INTERNAL_CATCH_UNIQUE_NAME( globalRegistrar ) = [](){\
@@ -1820,8 +1829,8 @@ namespace Catch {
 #endif

    namespace Detail {
-        template<typename InputIterator, typename Sentinel = InputIterator>
-        std::string rangeToString(InputIterator first, Sentinel last) {
+        template<typename InputIterator>
+        std::string rangeToString(InputIterator first, InputIterator last) {
            ReusableStringStream rss;
            rss << "{ ";
            if (first != last) {
@@ -12033,7 +12042,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");
@@ -14082,10 +14091,10 @@ namespace Catch {

    namespace {
        struct TestHasher {
-            explicit TestHasher(Catch::SimplePcg32& rng_instance) {
-                basis = rng_instance();
+            explicit TestHasher(Catch::SimplePcg32& rng) {
+                basis = rng();
                basis <<= 32;
-                basis |= rng_instance();
+                basis |= rng();
            }

            uint64_t basis;
@@ -15183,41 +15192,6 @@ 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 {
@@ -15264,7 +15238,7 @@ namespace Catch {
    }

    Version const& libraryVersion() {
-        static Version version( 2, 13, 2, "", 0 );
+        static Version version( 2, 13, 0, "", 0 );
        return version;
    }

@@ -17822,4 +17796,3 @@ using Catch::Detail::Approx;
 // end catch_reenable_warnings.h
 // end catch.hpp
 #endif // TWOBLUECUBES_SINGLE_INCLUDE_CATCH_HPP_INCLUDED
-
@@ -0,0 +1,219 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "catch.hpp"
+
+using namespace mfem;
+
+namespace surf_blf
+{
+
+Mesh * GetCylMesh(int type);
+
+void trans_skew_cyl(const Vector &x, Vector &r);
+
+void sigmaFunc(const Vector &x, DenseMatrix &s);
+
+double uExact(const Vector &x)
+{
+   return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
+}
+
+TEST_CASE("Embedded Surface Diffusion",
+          "[DiffusionIntegrator]")
+{
+   for (int type = (int) Element::TRIANGLE;
+        type <= (int) Element::QUADRILATERAL;
+        type++)
+   {
+      int order = 3;
+
+      Mesh *mesh = GetCylMesh(type);
+      int dim = mesh->Dimension();
+
+      mesh->SetCurvature(3);
+      mesh->Transform(trans_skew_cyl);
+
+      H1_FECollection fec(order, dim);
+      FiniteElementSpace fespace(mesh, &fec);
+
+      Array<int> ess_tdof_list;
+      if (mesh->bdr_attributes.Size())
+      {
+         Array<int> ess_bdr(mesh->bdr_attributes.Max());
+         ess_bdr = 1;
+         fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+
+      LinearForm b(&fespace);
+      ConstantCoefficient one(1.0);
+      b.AddDomainIntegrator(new DomainLFIntegrator(one));
+      b.Assemble();
+
+      GridFunction x(&fespace);
+      x = 0.0;
+
+      BilinearForm a(&fespace);
+      MatrixFunctionCoefficient sigma(3, sigmaFunc);
+      a.AddDomainIntegrator(new DiffusionIntegrator(sigma));
+      a.Assemble();
+
+      OperatorPtr A;
+      Vector B, X;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 0, 200, 1e-12, 0.0);
+
+      a.RecoverFEMSolution(X, b, x);
+
+      FunctionCoefficient uCoef(uExact);
+      double err = x.ComputeL2Error(uCoef);
+
+      REQUIRE(err < 1.5e-3);
+
+      delete mesh;
+   }
+}
+
+Mesh * GetCylMesh(int type)
+{
+   Mesh * mesh = NULL;
+
+   if (type == (int) Element::TRIANGLE)
+   {
+      mesh = new Mesh(2, 12, 16, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+      mesh->AddVertex( 0.0, -1.0, 0.5);
+      mesh->AddVertex( 1.0,  0.0, 0.5);
+      mesh->AddVertex( 0.0,  1.0, 0.5);
+      mesh->AddVertex(-1.0,  0.0, 0.5);
+
+      mesh->AddTriangle(0, 1, 8);
+      mesh->AddTriangle(1, 5, 8);
+      mesh->AddTriangle(5, 4, 8);
+      mesh->AddTriangle(4, 0, 8);
+      mesh->AddTriangle(1, 2, 9);
+      mesh->AddTriangle(2, 6, 9);
+      mesh->AddTriangle(6, 5, 9);
+      mesh->AddTriangle(5, 1, 9);
+      mesh->AddTriangle(2, 3, 10);
+      mesh->AddTriangle(3, 7, 10);
+      mesh->AddTriangle(7, 6, 10);
+      mesh->AddTriangle(6, 2, 10);
+      mesh->AddTriangle(3, 0, 11);
+      mesh->AddTriangle(0, 4, 11);
+      mesh->AddTriangle(4, 7, 11);
+      mesh->AddTriangle(7, 3, 11);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else if (type == (int) Element::QUADRILATERAL)
+   {
+      mesh = new Mesh(2, 8, 4, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+
+      mesh->AddQuad(0, 1, 5, 4);
+      mesh->AddQuad(1, 2, 6, 5);
+      mesh->AddQuad(2, 3, 7, 6);
+      mesh->AddQuad(3, 0, 4, 7);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else
+   {
+      MFEM_ABORT("Unrecognized mesh type " << type << "!");
+   }
+   mesh->FinalizeTopology();
+
+   return mesh;
+}
+
+void trans_skew_cyl(const Vector &x, Vector &r)
+{
+   r.SetSize(3);
+
+   double tol = 1e-6;
+   double theta = 0.0;
+   if (fabs(x[1] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (x[0] - 2.0);
+   }
+   else if (fabs(x[0] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * x[1];
+   }
+   else if (fabs(x[1] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (2.0 - x[0]);
+   }
+   else if (fabs(x[0] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (4.0 - x[1]);
+   }
+   else
+   {
+      mfem::out << "side not recognized "
+                << x[0] << " " << x[1] << " " << x[2] << std::endl;
+   }
+
+   r[0] = cos(theta);
+   r[1] = sin(theta);
+   r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
+}
+
+void sigmaFunc(const Vector &x, DenseMatrix &s)
+{
+   s.SetSize(3);
+   double a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
+   s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
+   s(0,1) = x[0] * x[1] * (8.0 / a - 0.5);
+   s(0,2) = 0.0;
+   s(1,0) = s(0,1);
+   s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] / a;
+   s(1,2) = 0.0;
+   s(2,0) = 0.0;
+   s(2,1) = 0.0;
+   s(2,2) = a / 32.0;
+}
+
+} // namespace surf_blf
@@ -251,24 +251,6 @@ 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
Stowell, Mark L	6c62756b28	Fixing copyright year	2021-04-09 15:49:41 -07:00
Stowell, Mark L	96e9dc09a6	Removing partial assembly code from example	2021-04-09 15:49:26 -07:00
Stowell, Mark L	52ccc6d980	Updating other DiffusionIntegrator methods	2021-04-09 15:48:59 -07:00
Syun'ichi Shiraiwa	eef5786aa8	Merge branch 'fix_embeded_diffusion_matrix_coeff' of github.com:mfem/mfem into fix_embeded_diffusion_matrix_coeff	2021-03-28 12:32:35 -04:00
Syun'ichi Shiraiwa	017460d517	added two assertions in routines which Veselin pointed out	2021-03-28 12:31:36 -04:00
Stowell, Mark L	62668b0e81	Adding unit test	2021-03-26 17:03:46 -07:00
Syun'ichi Shiraiwa	134c3b6c7f	added a DenseMatrix with proper size	2021-03-23 18:01:15 -04:00