fixed comment

Fix GitHub actions failure on MacOS

.github/workflows/build-container.yaml

@@ -62,7 +62,7 @@ jobs:
 
       - name: GHCR Login
         if: (github.event_name != 'pull_request')
-        uses: docker/login-action@v1
+        uses: docker/login-action@v2
         with:
           registry: ghcr.io
           username: ${{ github.actor }}

.github/workflows/builds-and-tests.yml

@@ -94,7 +94,7 @@ jobs:
     # This external action allows to interrupt a workflow already running on
     # the same branch to save resource
     - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.9.0
+      uses: styfle/cancel-workflow-action@0.11.0
       with:
         access_token: ${{ github.token }}
 
@@ -102,7 +102,7 @@ jobs:
     # /home/runner/work/mfem/mfem/mfem
     # Note: Done now to access "install-hypre" and "install-metis" actions.
     - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         path: ${{ env.MFEM_TOP_DIR }}
         # Fetch the complete history for codecov to access commits ID
@@ -115,25 +115,25 @@ jobs:
       if: matrix.mpi == 'par' && matrix.os == 'ubuntu-latest'
       run: |
         sudo apt-get install mpich libmpich-dev
-        export MAKE_CXX_FLAG="MPICXX=mpic++"
 
     - name: get lcov (Linux)
       if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-latest'
       run: |
         sudo apt-get install lcov
 
-    - name: Set up Homebrew
-      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
-      uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3
+    # Keep the following section in case we need it again in the future,
+    # see: https://github.com/mfem/mfem/pull/3385#discussion_r1058013032
+    # - name: Set up Homebrew
+    #   if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
+    #   uses: Homebrew/actions/setup-homebrew@master
 
     - name: get MPI (MacOS)
       if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install openmpi
-        export MAKE_CXX_FLAG="MPICXX=mpic++"
 
-    - name: get MPI (MacOS)
+    - name: get lcov (MacOS)
       if: matrix.codecov == 'YES' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
@@ -141,14 +141,14 @@ jobs:
 
     - name: get MPI (Windows)
       if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
-      uses: mpi4py/setup-mpi@v1.0.3
+      uses: mpi4py/setup-mpi@v1.1.2
 
     # Get Hypre through cache, or build it.
     # Install will only run on cache miss.
     - name: cache hypre
       id: hypre-cache
       if: matrix.mpi == 'par'
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       with:
         path: ${{ env.HYPRE_TOP_DIR }}
         key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2
@@ -176,7 +176,7 @@ jobs:
     - name: cache metis
       id: metis-cache
       if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       with:
         path: ${{ env.METIS_TOP_DIR }}
         key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
@@ -190,12 +190,13 @@ jobs:
 
     - name: cache vcpkg (Windows)
       id: vcpkg-cache
+      if: matrix.os == 'windows-latest'
       uses: actions/cache@v3
       with:
         path: vcpkg_cache
         key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
 
-    - name: prepare binary cache location
+    - name: prepare vcpkg binary cache location (Windows)
       if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
       run: |
         mkdir -p vcpkg_cache
@@ -230,11 +231,7 @@ jobs:
       run: |
         cd ${{ env.MFEM_TOP_DIR }}  && make check
 
-    - name: unit tests
-      if: matrix.build-system == 'make' && matrix.target == 'opt'
-      run: |
-        cd ${{ env.MFEM_TOP_DIR }} && make unittest
-
+    # Note: 'tests' include the unit tests
     - name: tests
       if: matrix.build-system == 'make' && matrix.target == 'opt'
       run: |
@@ -247,7 +244,7 @@ jobs:
         cd ${{ env.MFEM_TOP_DIR }} && cmake --build build --target check --config ${CTEST_CONFIG}
       shell: bash
 
-    - name: cmake unit tests (Ubuntu 20.04)
+    - name: cmake unit tests (Ubuntu)
       if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
       run: |
         CTEST_CONFIG="Release"
@@ -265,7 +262,7 @@ jobs:
     # Code coverage (process and upload reports)
     - name: codecov
       if: matrix.codecov == 'YES'
-      uses: mfem/github-actions/upload-coverage@v2.0
+      uses: mfem/github-actions/upload-coverage@v2.2
       with:
         name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
         project_dir: ${{ env.MFEM_TOP_DIR }}

.github/workflows/mfem-analysis.yml

@@ -35,23 +35,22 @@ jobs:
 
     steps:
     - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.9.0
+      uses: styfle/cancel-workflow-action@0.11.0
       with:
         access_token: ${{ github.token }}
 
     - name: checkout MFEM
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         path: mfem
 
     - name: Get MPI (Linux)
       run: |
         sudo apt-get install mpich libmpich-dev
-        export MAKE_CXX_FLAG="MPICXX=mpic++"
 
     - name: Cache Hypre Install
       id: hypre-cache
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       with:
         path: ${{ env.HYPRE_TOP_DIR }}
         key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2.2
@@ -66,7 +65,7 @@ jobs:
 
     - name: Cache Metis Install
       id: metis-cache
-      uses: actions/cache@v2
+      uses: actions/cache@v3
       with:
         path: ${{ env.METIS_TOP_DIR }}
         key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2

.github/workflows/repo-check.yml

@@ -34,12 +34,12 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
     - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.9.0
+      uses: styfle/cancel-workflow-action@0.11.0
       with:
         access_token: ${{ github.token }}
 
     - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: copyright check
       id: copyright
@@ -84,7 +84,7 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
     - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: get astyle
       run: |
@@ -101,11 +101,14 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
     - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: get doxygen and graphviz
       run: |
         sudo apt-get install doxygen graphviz
+
+    - name: update doxygen config file
+      run: |
         cd doc
         doxygen -u CodeDocumentation.conf.in
 
@@ -123,7 +126,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
     - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         fetch-depth: 0
 

CHANGELOG

@@ -8,6 +8,51 @@
                                https://mfem.org
 
 
+
+
+Version 4.5.1 (development)
+===========================
+
+Meshing improvements
+--------------------
+- Added support for pyramids in non-conforming meshes. Currently only isotropic
+  refinement is supported in this case.
+
+Discretization improvements
+---------------------------
+- TBD
+
+Linear and nonlinear solvers
+----------------------------
+- Added a fast normalization-based distance solver, see the Distance miniapp
+  in the miniapps/shifted/ directory.
+
+New and updated examples and miniapps
+-------------------------------------
+- TBD
+
+Integrations, testing and documentation
+---------------------------------------
+- Removed the support for the Mesquite toolkit. We recommend using MFEM's TMOP
+  functionality instead for mesh optimization. See the mesh-optimizer miniapp.
+
+Miscellaneous
+-------------
+- VisItDataCollection now correctly handles collection names containing
+  underscores.
+
+- VisItDataCollection::SetPadDigits() no longer alters the number of digits
+  used to represent the MPI rank because VisIt seems to require 6 digits.
+  This parameter can still be explicitly overridden with
+  VisItDataCollection::SetPadDigitsRank().
+
+API changes
+-----------
+- The implicit cast methods of class Vector to 'double *' and 'const double *'
+  have been deprecated and generate deprecation warnings if used. They will be
+  removed in a future release.
+
+
 Version 4.5, released on October 22, 2022
 =========================================
 

CMakeLists.txt

@@ -10,7 +10,9 @@
 # CONTRIBUTING.md for details.
 
 # The variable CMAKE_CXX_STANDARD and related were introduced in CMake v3.1
-cmake_minimum_required(VERSION 3.1)
+# Version 3.8 fixes the handling of CMAKE_CXX_STANDARD for try_compile.
+# Version 3.8 or newer is required for direct CUDA support.
+cmake_minimum_required(VERSION 3.8)
 set(USER_CONFIG "${CMAKE_CURRENT_SOURCE_DIR}/config/user.cmake" CACHE PATH
   "Path to optional user configuration file.")
 
@@ -51,7 +53,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.5.0)
+set(${PROJECT_NAME}_VERSION 4.5.1)
 
 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -66,8 +68,7 @@ set(XSDK_ENABLE_C OFF)
 set(XSDK_ENABLE_Fortran OFF)
 
 # Check if we need to enable C or Fortran.
-if (CMAKE_VERSION VERSION_LESS 3.2 OR
-    MFEM_USE_CONDUIT OR
+if (MFEM_USE_CONDUIT OR
     MFEM_USE_SIDRE OR
     MFEM_USE_PETSC)
    # This seems to be needed by:
@@ -81,11 +82,9 @@ if (MFEM_USE_STRUMPACK)
   # Just needed to find the MPI_Fortran libraries to link with
   set(XSDK_ENABLE_Fortran ON)
 endif()
-# SUNDIALS >= 6.4.0 requires C++14:
-if (MFEM_USE_SUNDIALS AND ("${CMAKE_CXX_STANDARD}" LESS "14"))
-  set(CMAKE_CXX_STANDARD 14)
-endif()
-if (MFEM_USE_GINKGO AND ("${CMAKE_CXX_STANDARD}" LESS "14"))
+# SUNDIALS, STRUMPACK, and Ginkgo require C++14:
+if ((MFEM_USE_SUNDIALS OR MFEM_USE_STRUMPACK OR MFEM_USE_GINKGO) AND
+    ("${CMAKE_CXX_STANDARD}" LESS "14"))
   set(CMAKE_CXX_STANDARD 14)
 endif()
 
@@ -103,8 +102,6 @@ if (MFEM_USE_CUDA)
    if (MFEM_USE_HIP)
       message(FATAL_ERROR " *** MFEM_USE_HIP cannot be combined with MFEM_USE_CUDA.")
    endif()
-   # MFEM_USE_CUDA requires CMake 3.8 or newer (for direct CUDA support)
-   cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
    # Use ${CMAKE_CXX_COMPILER} as the cuda host compiler.
    if (NOT CMAKE_CUDA_HOST_COMPILER)
       set(CMAKE_CUDA_HOST_COMPILER ${CMAKE_CXX_COMPILER})
@@ -199,6 +196,26 @@ if (MFEM_USE_HIP)
   find_package(HIPSPARSE REQUIRED)
 endif()
 
+# OpenMP
+if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
+  if (NOT MFEM_THREAD_SAFE AND MFEM_USE_LEGACY_OPENMP)
+    message(FATAL_ERROR " *** MFEM_USE_LEGACY_OPENMP requires MFEM_THREAD_SAFE=ON.")
+  endif()
+  find_package(OpenMP REQUIRED)
+  set(OPENMP_LIBRARIES ${OpenMP_CXX_LIBRARIES})
+  if(APPLE)
+    # On macOS, the compiler needs additional help to find the <omp.h> header.
+    # See issue #2642 for more information.
+    set(OPENMP_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS})
+  endif(APPLE)
+  if (OPENMP_FOUND)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+    if (MFEM_USE_CUDA)
+      set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -Xcompiler=${OpenMP_CXX_FLAGS}")
+    endif()
+  endif()
+endif()
+
 # MPI -> hypre; PETSc (optional)
 if (MFEM_USE_MPI)
   find_package(MPI REQUIRED)
@@ -264,20 +281,6 @@ if (MFEM_USE_LAPACK)
   find_package(LAPACK REQUIRED)
 endif()
 
-# OpenMP
-if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
-  if (NOT MFEM_THREAD_SAFE AND MFEM_USE_LEGACY_OPENMP)
-    message(FATAL_ERROR " *** MFEM_USE_LEGACY_OPENMP requires MFEM_THREAD_SAFE=ON.")
-  endif()
-  find_package(OpenMP REQUIRED)
-  set(OPENMP_LIBRARIES ${OpenMP_CXX_LIBRARIES})
-  if(APPLE)
-    # On macOS, the compiler needs additional help to find the <omp.h> header.
-    # See issue #2642 for more information.
-    set(OPENMP_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS})
-  endif(APPLE)
-endif()
-
 # SuiteSparse (before SUNDIALS which may depend on KLU)
 if (MFEM_USE_SUITESPARSE)
   find_package(SuiteSparse REQUIRED
@@ -296,11 +299,6 @@ if (MFEM_USE_SUNDIALS)
   find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS})
 endif()
 
-# Mesquite
-if (MFEM_USE_MESQUITE)
-  find_package(Mesquite REQUIRED)
-endif()
-
 # SuperLU_DIST can only be enabled in parallel
 if (MFEM_USE_SUPERLU)
   if (MFEM_USE_MPI)
@@ -488,6 +486,10 @@ endif()
 #   an ALIAS target is missing?
 # Call Stack (most recent call first):
 #   CMakeLists.txt:474 (mfem_add_library)
+#
+# NOTE: We need to figure out which TPL library adds the dependency on
+#       "Threads::Threads" and call the next line only when that TPL library is
+#       enabled.  -V. Dobrev
 find_package(Threads REQUIRED)
 
 # List all possible libraries in order of dependencies.
@@ -495,8 +497,8 @@ find_package(Threads REQUIRED)
 #    With newer versions of SuiteSparse which include METIS header using 64-bit
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
-set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist METIS SuiteSparse SUNDIALS
-    PETSC SLEPC MESQUITE MUMPS STRUMPACK AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
+set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
+    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
     NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
     ADIOS2 CUBLAS CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
     MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)
@@ -511,18 +513,12 @@ foreach(TPL IN LISTS MFEM_TPLS)
     list(APPEND TPL_INCLUDE_DIRS ${${TPL}_INCLUDE_DIRS})
   endif()
 endforeach(TPL)
+list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
+list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_INCLUDE_DIRS)
 # message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")
 
-if (OPENMP_FOUND)
-  message(STATUS "MFEM: using package OpenMP")
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-  if (MFEM_USE_CUDA)
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -Xcompiler=${OpenMP_CXX_FLAGS}")
-  endif()
-endif()
-
 message(STATUS "MFEM build type: CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 message(STATUS "MFEM version: v${MFEM_VERSION_STRING}")
 message(STATUS "MFEM git string: ${MFEM_GIT_STRING}")

INSTALL

@@ -337,10 +337,6 @@ MFEM_USE_SUNDIALS = YES/NO
    library. When enabled, this option uses the SUNDIALS_* library options,
    see below.
 
-MFEM_USE_MESQUITE = YES/NO
-   Enable MFEM functionality based on the Mesquite library. When enabled, this
-   option uses the MESQUITE_* library options, see below.
-
 MFEM_USE_SUITESPARSE = YES/NO
    Enable MFEM functionality based on the SuiteSparse library. Currently, this
    option adds the classes UMFPackSolver and KLUSolver (both sparse serial
@@ -636,11 +632,6 @@ The specific libraries and their options are:
   Options: SUNDIALS_OPT, SUNDIALS_LIB.
   Versions: SUNDIALS >= 5.0.0, SUNDIALS >= 5.4.0 for CUDA support.
 
-- Mesquite (optional), used when MFEM_USE_MESQUITE = YES.
-  URL: http://trilinos.org/oldsite/packages/mesquite
-  Options: MESQUITE_OPT, MESQUITE_LIB.
-  The Mesquite support is deprecated and will be removed in the future.
-
 - SuiteSparse (optional), used when MFEM_USE_SUITESPARSE = YES.
   URL: http://faculty.cse.tamu.edu/davis/suitesparse.html
   Options: SUITESPARSE_OPT, SUITESPARSE_LIB.
@@ -969,7 +960,6 @@ MFEM_USE_LEGACY_OPENMP
 MFEM_USE_OPENMP
 MFEM_USE_MEMALLOC
 MFEM_TIMER_TYPE - Set automatically, can be overwritten.
-MFEM_USE_MESQUITE
 MFEM_USE_SUITESPARSE
 MFEM_USE_SUPERLU
 MFEM_USE_MUMPS
@@ -1034,7 +1024,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - HYPRE
  - METIS - The option MFEM_USE_METIS_5 is auto-detected.
  - ParMETIS
- - MESQUITE
  - SuiteSparse
  - SuperLUDist, STRUMPACK
  - Ginkgo

config/XSDKDefaults.cmake

@@ -212,10 +212,6 @@ IF (DEFINED TPL_ENABLE_SUNDIALS)
   SET(MFEM_USE_SUNDIALS ${TPL_ENABLE_SUNDIALS} CACHE BOOL "Enable SUNDIALS usage" FORCE)
 ENDIF()
 
-IF (DEFINED TPL_ENABLE_MESQUITE)
-  SET(MFEM_USE_MESQUITE ${TPL_ENABLE_MESQUITE} CACHE BOOL "Enable MESQUITE usage" FORCE)
-ENDIF()
-
 IF (DEFINED TPL_ENABLE_SUITESPARSE)
   SET(MFEM_USE_SUITESPARSE ${TPL_ENABLE_SUITESPARSE} CACHE BOOL "Enable SuiteSparse usage" FORCE)
 ENDIF()

config/cmake/MFEMConfig.cmake.in

@@ -29,7 +29,6 @@ set(MFEM_USE_LEGACY_OPENMP @MFEM_USE_LEGACY_OPENMP@)
 set(MFEM_USE_MEMALLOC @MFEM_USE_MEMALLOC@)
 set(MFEM_TIMER_TYPE @MFEM_TIMER_TYPE@)
 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@)

config/cmake/config.hpp.in

@@ -77,9 +77,6 @@
 // Internal MFEM option: enable group/batch allocation for some small objects.
 #cmakedefine MFEM_USE_MEMALLOC
 
-// Enable MFEM functionality based on the Mesquite library.
-#cmakedefine MFEM_USE_MESQUITE
-
 // Enable MFEM functionality based on the SuiteSparse library.
 #cmakedefine MFEM_USE_SUITESPARSE
 

config/cmake/modules/FindCaliper.cmake

@@ -16,7 +16,22 @@
 
 include(MfemCmakeUtilities)
 mfem_find_package(Caliper CALIPER CALIPER_DIR
-       	"include" "caliper/cali.h"
-       	"lib" "caliper"
-  	"Paths to headers required by Caliper."
-       	"Libraries required by Caliper.")
+      "include" "caliper/cali.h"
+      "lib" "caliper"
+      "Paths to headers required by Caliper."
+      "Libraries required by Caliper.")
+
+# Append adiak path/lib if the user provided ADIAK_DIR
+if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
+    find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})
+    list(APPEND CALIPER_INCLUDE_DIRS ${adiak_INCLUDE_DIRS})
+    list(APPEND CALIPER_LIBRARIES ${adiak_LIBRARIES})
+endif()
+
+# Append gotcha path/lib if the user provided GOTCHA_DIR
+if(GOTCHA_DIR AND EXISTS ${GOTCHA_DIR})
+    find_package(gotcha NO_DEFAULT_PATH REQUIRED PATHS ${GOTCHA_DIR}/lib/cmake/gotcha ${GOTCHA_DIR})
+    list(APPEND CALIPER_INCLUDE_DIRS ${gotcha_INCLUDE_DIRS})
+    list(APPEND CALIPER_LIBRARIES ${gotcha_LIBRARIES})
+endif()
+

config/cmake/modules/FindHDF5.cmake

@@ -14,9 +14,21 @@
 #   - HDF5_LIBRARIES - The HDF5 libraries
 #   - HDF5_INCLUDE_DIRS - The HDF5 include directories
 
+# NOTE: Using this FindHDF5.cmake instead of the CMake provided version may lead
+#       to issues with some TPL libraries that depend (or may depend) on HDF5.
+#       For this reason, we should consider removing this file, or at least
+#       making it use the CMake provided version by default and apply the logic
+#       below only when specifically requested by a user.  -V. Dobrev
+
 # First Check for HDF5_DIR
 if(NOT HDF5_DIR)
-    MESSAGE(FATAL_ERROR "Could not find HDF5. HDF5 support needs explicit HDF5_DIR")
+  message(FATAL_ERROR
+    "Could not find HDF5. HDF5 support needs explicit HDF5_DIR")
+endif()
+
+if (NOT HDF5_FIND_QUIETLY)
+  message(STATUS "Looking for HDF5 ...")
+  message(STATUS "   in HDF5_DIR = ${HDF5_DIR}")
 endif()
 
 # Find includes
@@ -50,8 +62,9 @@ include(FindPackageHandleStandardArgs)
 
 # Handle the QUIETLY and REQUIRED arguments and set HDF5_FOUND to TRUE if all
 # listed variables are TRUE
-find_package_handle_standard_args(HDF5  DEFAULT_MSG
+find_package_handle_standard_args(HDF5
+                                  " *** HDF5 not found. Please set HDF5_DIR."
+                                  HDF5_LIBRARIES
                                   HDF5_INCLUDE_DIRS
                                   __HDF5_LIBRARY
-                                  __HDF5_HL_LIBRARY
-                                  HDF5_LIBRARIES )
+                                  __HDF5_HL_LIBRARY)

config/cmake/modules/FindMesquite.cmake

@@ -1,20 +0,0 @@
-# Copyright (c) 2010-2022, 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.
-
-# Defines the following variables:
-#   - MESQUITE_FOUND
-#   - MESQUITE_LIBRARIES
-#   - MESQUITE_INCLUDE_DIRS
-
-include(MfemCmakeUtilities)
-mfem_find_package(Mesquite MESQUITE MESQUITE_DIR
-  "include" "Mesquite_all_headers.hpp" "lib" "mesquite"
-  "Paths to headers required by Mesquite." "Libraries required by Mesquite.")

config/cmake/modules/FindNetCDF.cmake

@@ -17,18 +17,24 @@
 include(MfemCmakeUtilities)
 
 # FindHDF5.cmake uses HDF5_ROOT, so we "translate" from the MFEM convention
-set(HDF5_ROOT ${HDF5_DIR} CACHE PATH "")
+# (MFEM's FindHDF5.cmake does not need HDF5_ROOT)
+# set(HDF5_ROOT ${HDF5_DIR} CACHE PATH "")
+
 # We need to guard against the case where HDF5 was already found but without
 # the HL extensions (in which case mfem_find_package will treat the package
 # as already having been found), so we reset the variable to force FindHDF5.cmake
 # to be called for a second time
 set(HDF5_FOUND OFF)
 enable_language(C) # FindHDF5.cmake uses the C compiler
+
 mfem_find_package(NetCDF NETCDF NETCDF_DIR "include" netcdf.h "lib" netcdf
   "Paths to headers required by NetCDF." "Libraries required by NetCDF.")
+
 # The HL extension libraries are in a separate variable and must precede
 # the "regular" hdf5 library, as hdf5_hl depends on hdf5
 # The netcdf library will always be the first element of NETCDF_LIBRARIES
 # and we need to insert after that library but before the hdf5 library, so
 # position 1 is used
-list(INSERT NETCDF_LIBRARIES 1 ${HDF5_C_LIBRARY_hdf5_hl})
+# (MFEM's FindHDF5.cmake does not set HDF5_C_LIBRARY_hdf5_hl and the HL library
+# is already added to NETCDF_LIBRARIES)
+# list(INSERT NETCDF_LIBRARIES 1 ${HDF5_C_LIBRARY_hdf5_hl})

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -46,6 +46,10 @@ endfunction()
 # Wrapper for add_executable
 macro(mfem_add_executable NAME)
   add_executable(${NAME} ${ARGN})
+  if (MFEM_USE_CUDA)
+    set_target_properties(${NAME} PROPERTIES
+      CUDA_RESOLVE_DEVICE_SYMBOLS ON)
+  endif()
 endmacro()
 
 # Wrapper for add_library
@@ -158,27 +162,12 @@ macro(add_mfem_miniapp MFEM_EXE_NAME)
 
   # Append the additional libraries and options
   if (LIBRARIES_LIST)
-    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
-      target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${LIBRARIES_LIST})
-    else()
-      target_link_libraries(${MFEM_EXE_NAME} ${LIBRARIES_LIST})
-    endif()
+    target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${LIBRARIES_LIST})
   endif()
   if (EXTRA_OPTIONS_LIST)
     string(REPLACE ";" " " EXTRA_OPTIONS_STRING "${EXTRA_OPTIONS_LIST}")
     message(STATUS "${MFEM_EXE_NAME}: add flags \"${EXTRA_OPTIONS_STRING}\"")
-    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
-      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${EXTRA_OPTIONS_LIST})
-    else()
-      get_target_property(THIS_COMPILE_FLAGS ${MFEM_EXE_NAME} COMPILE_FLAGS)
-      if (THIS_COMPILE_FLAGS)
-        set(THIS_COMPILE_FLAGS "${THIS_COMPILE_FLAGS} ${EXTRA_OPTIONS_STRING}")
-      else()
-        set(THIS_COMPILE_FLAGS "${EXTRA_OPTIONS_STRING}")
-      endif()
-      set_target_properties(${MFEM_EXE_NAME}
-        PROPERTIES COMPILE_FLAGS ${THIS_COMPILE_FLAGS})
-    endif()
+    target_compile_options(${MFEM_EXE_NAME} PRIVATE ${EXTRA_OPTIONS_LIST})
   endif()
   if (EXTRA_DEFINES_LIST)
     target_compile_definitions(${MFEM_EXE_NAME} PRIVATE ${EXTRA_DEFINES_LIST})
@@ -187,17 +176,15 @@ macro(add_mfem_miniapp MFEM_EXE_NAME)
   # Handle the MPI separately
   if (MFEM_USE_MPI)
     # Add MPI_CXX_LIBRARIES, in case this target does not link with mfem.
-    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
-      target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_LIBRARIES})
-    else()
-      target_link_libraries(${MFEM_EXE_NAME} ${MPI_CXX_LIBRARIES})
-    endif()
+    target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_LIBRARIES})
 
     if (MPI_CXX_INCLUDE_PATH)
       target_include_directories(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_INCLUDE_PATH})
     endif()
     if (MPI_CXX_COMPILE_FLAGS)
-      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_COMPILE_FLAGS})
+      separate_arguments(MPI_CXX_COMPILE_ARGS UNIX_COMMAND
+        "${MPI_CXX_COMPILE_FLAGS}")
+      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_COMPILE_ARGS})
     endif()
 
     if (MPI_CXX_LINK_FLAGS)
@@ -878,7 +865,7 @@ function(mfem_export_mk_files)
   set(CONFIG_MK_BOOL_VARS MFEM_USE_MPI MFEM_USE_METIS MFEM_USE_METIS_5
       MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB MFEM_USE_LIBUNWIND
       MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_LEGACY_OPENMP MFEM_USE_OPENMP
-      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
+      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_SUITESPARSE
       MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS MFEM_USE_STRUMPACK
       MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_GNUTLS MFEM_USE_NETCDF
       MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS

config/config.hpp.in

@@ -85,9 +85,6 @@
 // Enable MFEM functionality based on the SUNDIALS libraries.
 // #define MFEM_USE_SUNDIALS
 
-// Enable MFEM functionality based on the Mesquite library.
-// #define MFEM_USE_MESQUITE
-
 // Enable MFEM functionality based on the SuiteSparse library.
 // #define MFEM_USE_SUITESPARSE
 

config/config.mk.in

@@ -29,7 +29,6 @@ MFEM_USE_OPENMP        = @MFEM_USE_OPENMP@
 MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
-MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@

config/defaults.cmake

@@ -30,7 +30,6 @@ option(MFEM_USE_OPENMP "Enable the OpenMP backend" OFF)
 option(MFEM_USE_LEGACY_OPENMP "Enable legacy OpenMP usage" OFF)
 option(MFEM_USE_MEMALLOC "Enable the internal MEMALLOC option." ON)
 option(MFEM_USE_SUNDIALS "Enable SUNDIALS usage" OFF)
-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)
@@ -124,9 +123,6 @@ set(SUNDIALS_DIR "${MFEM_DIR}/../sundials-5.0.0/instdir" CACHE PATH
 # set(SUNDIALS_REQUIRED_PACKAGES "SuiteSparse/KLU/AMD/BTF/COLAMD/config"
 #     CACHE STRING "Additional packages required by SUNDIALS.")
 
-set(MESQUITE_DIR "${MFEM_DIR}/../mesquite-2.99" CACHE PATH
-    "Path to the Mesquite library.")
-
 set(SuiteSparse_DIR "${MFEM_DIR}/../SuiteSparse" CACHE PATH
     "Path to the SuiteSparse library.")
 set(SuiteSparse_REQUIRED_PACKAGES "BLAS" "METIS"
@@ -188,6 +184,7 @@ set(GNUTLS_DIR "" CACHE PATH "Path to the GnuTLS library.")
 
 set(GSLIB_DIR "" CACHE PATH "Path to the GSLIB library.")
 
+set(HDF5_DIR "/usr" CACHE PATH "Path to the HDF5 library.")
 set(NETCDF_DIR "" CACHE PATH "Path to the NetCDF library.")
 set(NetCDF_REQUIRED_PACKAGES "HDF5/C/HL" CACHE STRING
     "Additional packages required by NetCDF.")

config/defaults.mk

@@ -131,7 +131,6 @@ MFEM_USE_LEGACY_OPENMP = NO
 MFEM_USE_MEMALLOC      = YES
 MFEM_TIMER_TYPE        = $(if $(NOTMAC),2,4)
 MFEM_USE_SUNDIALS      = NO
-MFEM_USE_MESQUITE      = NO
 MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
@@ -270,11 +269,6 @@ endif
 # If SUNDIALS was built with KLU:
 # MFEM_USE_SUITESPARSE = YES
 
-# MESQUITE library configuration
-MESQUITE_DIR = @MFEM_DIR@/../mesquite-2.99
-MESQUITE_OPT = -I$(MESQUITE_DIR)/include
-MESQUITE_LIB = -L$(MESQUITE_DIR)/lib -lmesquite
-
 # SuiteSparse library configuration
 LIB_RT = $(if $(NOTMAC),-lrt,)
 SUITESPARSE_DIR = @MFEM_DIR@/../SuiteSparse
@@ -324,6 +318,9 @@ MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps\
 
 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build
+ifeq ($(MFEM_USE_STRUMPACK),YES)
+   BASE_FLAGS = -std=c++14
+endif
 STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)
 # If STRUMPACK was build with OpenMP support, the following may be need:
 # STRUMPACK_OPT += $(OPENMP_OPT)
@@ -476,7 +473,17 @@ OCCA_LIB = $(XLINKER)-rpath,$(OCCA_DIR)/lib -L$(OCCA_DIR)/lib -locca
 # CALIPER library configuration
 CALIPER_DIR = @MFEM_DIR@/../caliper
 CALIPER_OPT = -I$(CALIPER_DIR)/include
-CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib64 -lcaliper
+CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 $(XLINKER)-rpath,$(CALIPER_DIR)/lib -L$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib -lcaliper
+
+ifdef ADIAK_DIR
+   CALIPER_OPT += -I$(ADIAK_DIR)/include
+   CALIPER_LIB += $(XLINKER)-rpath,$(ADIAK_DIR)/lib64 $(XLINKER)-rpath,$(ADIAK_DIR)/lib -L$(ADIAK_DIR)/lib64 -L$(ADIAK_DIR)/lib -ladiak
+endif
+ifdef GOTCHA_DIR
+   CALIPER_OPT += -I$(GOTCHA_DIR)/include
+   CALIPER_LIB += $(XLINKER)-rpath,$(GOTCHA_DIR)/lib64 $(XLINKER)-rpath,$(GOTCHA_DIR)/lib -L$(GOTCHA_DIR)/lib64 -L$(GOTCHA_DIR)/lib -lgotcha
+endif
+
 
 # BLITZ library configuration
 BLITZ_DIR = @MFEM_DIR@/../blitz

doc/CodeDocumentation.conf.in

@@ -38,7 +38,7 @@ PROJECT_NAME           = "MFEM"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = v4.5.0
+PROJECT_NUMBER         = v4.5.1
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
@@ -763,36 +763,54 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                          @MFEM_SOURCE_DIR@/config \
                          @MFEM_SOURCE_DIR@/general \
                          @MFEM_SOURCE_DIR@/linalg \
+                         @MFEM_SOURCE_DIR@/linalg/simd \
                          @MFEM_SOURCE_DIR@/mesh \
+                         @MFEM_SOURCE_DIR@/mesh/submesh \
                          @MFEM_SOURCE_DIR@/fem \
-                         @MFEM_SOURCE_DIR@/fem/moonolith \
+                         @MFEM_SOURCE_DIR@/fem/ceed \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/convection \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/diffusion \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/mass \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/nlconvection \
+                         @MFEM_SOURCE_DIR@/fem/ceed/interface \
+                         @MFEM_SOURCE_DIR@/fem/ceed/solvers \
                          @MFEM_SOURCE_DIR@/fem/fe \
                          @MFEM_SOURCE_DIR@/fem/lor \
+                         @MFEM_SOURCE_DIR@/fem/moonolith \
+                         @MFEM_SOURCE_DIR@/fem/qinterp \
+                         @MFEM_SOURCE_DIR@/fem/tmop \
                          @MFEM_SOURCE_DIR@/examples \
-                         @MFEM_SOURCE_DIR@/examples/caliper \
                          @MFEM_SOURCE_DIR@/examples/amgx \
+                         @MFEM_SOURCE_DIR@/examples/caliper \
                          @MFEM_SOURCE_DIR@/examples/ginkgo \
-                         @MFEM_SOURCE_DIR@/examples/moonolith \
                          @MFEM_SOURCE_DIR@/examples/hiop \
+                         @MFEM_SOURCE_DIR@/examples/moonolith \
                          @MFEM_SOURCE_DIR@/examples/petsc \
                          @MFEM_SOURCE_DIR@/examples/pumi \
                          @MFEM_SOURCE_DIR@/examples/sundials \
                          @MFEM_SOURCE_DIR@/examples/superlu \
                          @MFEM_SOURCE_DIR@/miniapps/adjoint \
+                         @MFEM_SOURCE_DIR@/miniapps/autodiff \
                          @MFEM_SOURCE_DIR@/miniapps/common \
                          @MFEM_SOURCE_DIR@/miniapps/electromagnetics \
                          @MFEM_SOURCE_DIR@/miniapps/gslib \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/kernels \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/materials \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/operators \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/preconditioners \
                          @MFEM_SOURCE_DIR@/miniapps/meshing \
                          @MFEM_SOURCE_DIR@/miniapps/mtop \
-                         @MFEM_SOURCE_DIR@/miniapps/autodiff \
+                         @MFEM_SOURCE_DIR@/miniapps/multidomain \
                          @MFEM_SOURCE_DIR@/miniapps/navier \
                          @MFEM_SOURCE_DIR@/miniapps/nurbs \
+                         @MFEM_SOURCE_DIR@/miniapps/parelag \
                          @MFEM_SOURCE_DIR@/miniapps/performance \
                          @MFEM_SOURCE_DIR@/miniapps/shifted \
                          @MFEM_SOURCE_DIR@/miniapps/solvers \
                          @MFEM_SOURCE_DIR@/miniapps/tools \
-                         @MFEM_SOURCE_DIR@/miniapps/toys \
-                         @MFEM_SOURCE_DIR@/miniapps/parelag
+                         @MFEM_SOURCE_DIR@/miniapps/toys
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses

examples/caliper/CMakeLists.txt

@@ -33,13 +33,13 @@ add_mfem_examples(CALIPER_EXE_SRCS ${PREFIX})
 if (MFEM_ENABLE_TESTING)
   foreach(SRC_FILE ${CALIPER_EXE_SRCS})
     get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
-    string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
+    string(REPLACE ".cpp" "" TEST_NAME ${PREFIX}${SRC_FILENAME})
 
     set(THIS_TEST_OPTIONS "-no-vis")
 
     if (NOT (${TEST_NAME} MATCHES ".*p$"))
       add_test(NAME ${TEST_NAME}_ser
-        COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
+        COMMAND $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS})
     else()
       add_test(NAME ${TEST_NAME}_np=${MFEM_MPI_NP}
         COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}

examples/ex1.cpp

@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               ex1 -pa -d cuda
-//             * ex1 -fa -d cuda
+//               ex1 -fa -d cuda
 //               ex1 -pa -d raja-cuda
 //             * ex1 -pa -d raja-hip
 //               ex1 -pa -d occa-cuda
@@ -192,7 +192,14 @@ int main(int argc, char *argv[])
    //    domain integrator.
    BilinearForm a(&fespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   if (fa) { a.SetAssemblyLevel(AssemblyLevel::FULL); }
+   if (fa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::FULL);
+      // Sort the matrix column indices when running on GPU or with OpenMP (i.e.
+      // when Device::IsEnabled() returns true). This makes the results
+      // bit-for-bit deterministic at the cost of somewhat longer run time.
+      a.EnableSparseMatrixSorting(Device::IsEnabled());
+   }
    a.AddDomainIntegrator(new DiffusionIntegrator(one));
 
    // 10. Assemble the bilinear form and the corresponding linear system,

examples/ex1p.cpp

@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               mpirun -np 4 ex1p -pa -d cuda
-//             * mpirun -np 4 ex1p -fa -d cuda
+//               mpirun -np 4 ex1p -fa -d cuda
 //               mpirun -np 4 ex1p -pa -d occa-cuda
 //               mpirun -np 4 ex1p -pa -d raja-omp
 //               mpirun -np 4 ex1p -pa -d ceed-cpu
@@ -219,7 +219,14 @@ int main(int argc, char *argv[])
    //     Diffusion domain integrator.
    ParBilinearForm a(&fespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   if (fa) { a.SetAssemblyLevel(AssemblyLevel::FULL); }
+   if (fa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::FULL);
+      // Sort the matrix column indices when running on GPU or with OpenMP (i.e.
+      // when Device::IsEnabled() returns true). This makes the results
+      // bit-for-bit deterministic at the cost of somewhat longer run time.
+      a.EnableSparseMatrixSorting(Device::IsEnabled());
+   }
    a.AddDomainIntegrator(new DiffusionIntegrator(one));
 
    // 12. Assemble the parallel bilinear form and the corresponding linear

examples/petsc/rc_ex1p

@@ -14,4 +14,5 @@
 -mg_levels_esteig_ksp_type cg
 -mg_levels_esteig_ksp_max_it 10
 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05
+-pc_gamg_use_sa_esteig 0
 -mg_levels_pc_type sor

fem/bilinearform.cpp

@@ -124,6 +124,7 @@ void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
       case AssemblyLevel::LEGACY:
          break;
       case AssemblyLevel::FULL:
+         SetDiagonalPolicy( DIAG_ONE ); // Only diagonal policy supported on device
          ext = new FABilinearFormExtension(this);
          break;
       case AssemblyLevel::ELEMENT:

fem/bilinearform.hpp

@@ -80,6 +80,9 @@ protected:
    /** @brief Extension for supporting Full Assembly (FA), Element Assembly (EA),
        Partial Assembly (PA), or Matrix Free assembly (MF). */
    BilinearFormExtension *ext;
+   /** Indicates if the sparse matrix is sorted after assembly when using
+       Full Assembly (FA). */
+   bool sort_sparse_matrix = false;
 
    /** @brief Indicates the Mesh::sequence corresponding to the current state of
        the BilinearForm. */
@@ -181,6 +184,21 @@ public:
        If used, this method must be called before assembly. */
    void SetAssemblyLevel(AssemblyLevel assembly_level);
 
+   /** @brief Force the sparse matrix column indices to be sorted when using
+       AssemblyLevel::FULL.
+
+       When assembling on device the assembly algorithm uses atomic operations
+       to insert values in the sparse matrix, which can result in different
+       column index orderings across runs. Calling this method with @a enable_it
+       set to @a true forces a sorting algorithm to be called at the end of the
+       assembly procedure to ensure sorted column indices (and therefore
+       deterministic results).
+   */
+   void EnableSparseMatrixSorting(bool enable_it)
+   {
+      sort_sparse_matrix = enable_it;
+   }
+
    /// Returns the assembly level
    AssemblyLevel GetAssemblyLevel() const { return assembly; }
 

fem/bilinearform_ext.cpp

@@ -955,6 +955,10 @@ void FABilinearFormExtension::Assemble()
       }
       a->mat = mat;
    }
+   if ( a->sort_sparse_matrix )
+   {
+      a->mat->SortColumnIndices();
+   }
 }
 
 

fem/bilininteg.cpp

@@ -2338,7 +2338,7 @@ void MixedCurlIntegrator::AssembleElementMatrix2(
    if (spaceH1)
    {
       dshape.SetSize(trial_dof,dim);
-      curlshape.SetSize(dim*trial_dof,1);
+      curlshape.SetSize(trial_dof,dim);
       dimc = dim;
    }
    else
@@ -2367,7 +2367,7 @@ void MixedCurlIntegrator::AssembleElementMatrix2(
       if (spaceH1)
       {
          trial_fe.CalcPhysDShape(Trans, dshape);
-         dshape.GradToCurl(curlshape);
+         dshape.GradToVectorCurl2D(curlshape);
       }
       else
       {

fem/bilininteg.hpp

@@ -3059,8 +3059,8 @@ public:
 
 /** Integrator for the DG form:
 
-    - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} >
-    + kappa < {h^{-1} Q} [u], [v] >,
+        - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} >
+        + kappa < {h^{-1} Q} [u], [v] >
 
     where Q is a scalar or matrix diffusion coefficient and u, v are the trial
     and test spaces, respectively. The parameters sigma and kappa determine the

fem/ceed/interface/restriction.hpp

@@ -47,6 +47,7 @@ void InitRestrictionWithIndices(const FiniteElementSpace &fes,
 
 /** @brief Initialize a strided CeedElemRestriction
 
+    @param[in] fes Input finite element space.
     @param[in] nelem is the number of elements.
     @param[in] nqpts is the total number of quadrature points.
     @param[in] qdatasize is the number of data per quadrature point.

fem/ceed/interface/util.cpp

@@ -140,7 +140,11 @@ int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
    CeedOperator *subops;
    if (isComposite)
    {
+#if CEED_VERSION_GE(0, 10, 2)
+      ierr = CeedCompositeOperatorGetSubList(oper, &subops); CeedChk(ierr);
+#else
       ierr = CeedOperatorGetSubList(oper, &subops); CeedChk(ierr);
+#endif
       ierr = CeedOperatorGetQFunction(subops[0], &qf); CeedChk(ierr);
    }
    else

fem/ceed/solvers/algebraic.cpp

@@ -275,9 +275,14 @@ CeedOperator CoarsenCeedCompositeOperator(
                                       &op_coarse); PCeedChk(ierr);
 
    int nsub;
-   ierr = CeedOperatorGetNumSub(op, &nsub); PCeedChk(ierr);
    CeedOperator *subops;
+#if CEED_VERSION_GE(0, 10, 2)
+   ierr = CeedCompositeOperatorGetNumSub(op, &nsub); PCeedChk(ierr);
+   ierr = CeedCompositeOperatorGetSubList(op, &subops); PCeedChk(ierr);
+#else
+   ierr = CeedOperatorGetNumSub(op, &nsub); PCeedChk(ierr);
    ierr = CeedOperatorGetSubList(op, &subops); PCeedChk(ierr);
+#endif
    for (int isub=0; isub<nsub; ++isub)
    {
       CeedOperator subop = subops[isub];

fem/ceed/solvers/full-assembly.cpp

@@ -310,8 +310,13 @@ int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat)
    {
       CeedInt numsub;
       CeedOperator *subops;
+#if CEED_VERSION_GE(0, 10, 2)
+      CeedCompositeOperatorGetNumSub(op, &numsub);
+      ierr = CeedCompositeOperatorGetSubList(op, &subops); CeedChk(ierr);
+#else
       CeedOperatorGetNumSub(op, &numsub);
       ierr = CeedOperatorGetSubList(op, &subops); CeedChk(ierr);
+#endif
       for (int i = 0; i < numsub; ++i)
       {
          ierr = CeedSingleOperatorFullAssemble(subops[i], out); CeedChk(ierr);

fem/ceed/solvers/solvers-atpmg.hpp

@@ -66,6 +66,7 @@ int CeedBasisATPMGCoarsen(CeedBasis basisin, CeedBasis* basisout,
     @param[in] coarse_er         CeedElemRestriction for coarse operator
                                  (see CeedATPMGElemRestriction)
     @param[out] coarse_basis_out CeedBasis for coarser operator
+    @param[out] basis_ctof_out   CeedBasis describing interpolation from coarse to fine
     @param[out] out              coarsened CeedOperator
 */
 int CeedATPMGOperator(CeedOperator oper, int order_reduction,

fem/datacollection.cpp

@@ -724,7 +724,7 @@ void VisItDataCollection::ParseVisItRootString(const std::string& json)
 
    // Set the DataCollection::name using the mesh path
    std::string path = mesh.get("path").get<std::string>();
-   size_t right_sep = path.find('_');
+   size_t right_sep = path.rfind('_');
    if (right_sep == std::string::npos)
    {
       error = READ_ERROR;

fem/datacollection.hpp

@@ -338,11 +338,12 @@ public:
    /// Set the precision (number of digits) used for the text output of doubles
    void SetPrecision(int prec) { precision = prec; }
    /// Set the number of digits used for both the cycle and the MPI rank
-   void SetPadDigits(int digits) { pad_digits_cycle=pad_digits_rank = digits; }
+   virtual void SetPadDigits(int digits)
+   { pad_digits_cycle=pad_digits_rank = digits; }
    /// Set the number of digits used for the cycle
-   void SetPadDigitsCycle(int digits) { pad_digits_cycle = digits; }
+   virtual void SetPadDigitsCycle(int digits) { pad_digits_cycle = digits; }
    /// Set the number of digits used for the MPI rank in filenames
-   void SetPadDigitsRank(int digits) { pad_digits_rank = digits; }
+   virtual void SetPadDigitsRank(int digits) { pad_digits_rank = digits; }
    /// Set the desired output mesh and data format.
    /** See the enumeration #Format for valid options. Derived classes can define
        their own format enumerations and override this method to perform input
@@ -441,21 +442,29 @@ public:
 #endif
 
    /// Set/change the mesh associated with the collection
-   virtual void SetMesh(Mesh *new_mesh);
+   virtual void SetMesh(Mesh *new_mesh) override;
 
 #ifdef MFEM_USE_MPI
    /// Set/change the mesh associated with the collection.
-   virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh);
+   virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh) override;
 #endif
 
    /// Add a grid function to the collection and update the root file
-   virtual void RegisterField(const std::string& field_name, GridFunction *gf);
+   virtual void RegisterField(const std::string& field_name,
+                              GridFunction *gf) override;
 
    /// Add a quadrature function to the collection and update the root file.
    /** Visualization of quadrature function is not supported in VisIt(3.12).
        A patch has been sent to VisIt developers in June 2020. */
    virtual void RegisterQField(const std::string& q_field_name,
-                               QuadratureFunction *qf);
+                               QuadratureFunction *qf) override;
+
+   /// Set the number of digits used for both the cycle and the MPI rank
+   /// @note VisIt seems to require 6 pad digits for the MPI rank. Therefore,
+   /// this function uses this default value. This behavior can be overridden
+   /// by calling SetPadDigitsCycle() and SetPadDigitsRank() instead.
+   virtual void SetPadDigits(int digits) override
+   { pad_digits_cycle=digits; pad_digits_rank=6; }
 
    /// Set VisIt parameter: default levels of detail for the MultiresControl
    void SetLevelsOfDetail(int levels_of_detail);
@@ -468,13 +477,13 @@ public:
    void DeleteAll();
 
    /// Save the collection and a VisIt root file
-   virtual void Save();
+   virtual void Save() override;
 
    /// Save a VisIt root file for the collection
    void SaveRootFile();
 
    /// Load the collection based on its VisIt data (described in its root file)
-   virtual void Load(int cycle_ = 0);
+   virtual void Load(int cycle_ = 0) override;
 
    /// We will delete the mesh and fields if we own them
    virtual ~VisItDataCollection() {}

fem/doftrans.cpp

@@ -339,8 +339,11 @@ ND_TriDofTransformation::TransformDual(double *v) const
 void
 ND_TriDofTransformation::InvTransformDual(double *v) const
 {
-   int nedofs = order; // number of DoFs per edge
-   int nfdofs = order*(order-1); // number of DoFs per face
+   // Return immediately when no face DoFs are present
+   if (nfdofs < 2) { return; }
+
+   MFEM_VERIFY(Fo.Size() >= 1,
+               "Face orientations are unset in ND_TriDofTransformation");
 
    double data[2];
    Vector v2(data, 2);
@@ -432,8 +435,11 @@ ND_TetDofTransformation::TransformDual(double *v) const
 void
 ND_TetDofTransformation::InvTransformDual(double *v) const
 {
-   int nedofs = order; // number of DoFs per edge
-   int nfdofs = order*(order-1); // number of DoFs per face
+   // Return immediately when no face DoFs are present
+   if (nfdofs < 2) { return; }
+
+   MFEM_VERIFY(Fo.Size() >= 4,
+               "Face orientations are unset in ND_TetDofTransformation");
 
    double data[2];
    Vector v2(data, 2);

fem/fe/fe_base.cpp

@@ -618,7 +618,7 @@ void InvertLinearTrans(ElementTransformation &trans,
 
    double store[3];
    Vector v(store, x.Size());
-   pt.Get(v, x.Size());
+   pt.Get(store, x.Size());
    v -= x;
 
    trans.InverseJacobian().Mult(v, x);
@@ -954,8 +954,8 @@ void VectorFiniteElement::Project_RT(
    {
       Trans.SetIntPoint(&Nodes.IntPoint(k));
       // dof_k = nk^t adj(J) xk
-      Vector vk(vc.GetData()+k*sdim, sdim);
-      dofs(k) = Trans.AdjugateJacobian().InnerProduct(vk, nk + d2n[k]*dim);
+      dofs(k) = Trans.AdjugateJacobian().InnerProduct(
+                   &vc[k*sdim], nk + d2n[k]*dim);
       if (!square_J) { dofs(k) /= Trans.Weight(); }
    }
 }
@@ -1171,9 +1171,8 @@ void VectorFiniteElement::Project_ND(
    for (int k = 0; k < dof; k++)
    {
       Trans.SetIntPoint(&Nodes.IntPoint(k));
-      Vector vk(vc.GetData()+k*dim, dim);
       // dof_k = xk^t J tk
-      dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, vk);
+      dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, &vc[k*dim]);
    }
 }
 
@@ -1320,7 +1319,7 @@ void VectorFiniteElement::LocalL2Projection_RT(
       double w = ip.weight;
       this->CalcVShape(ip, fine_shape);
       Trans.Transform(ip, v);
-      tr_ip.Set(v, dim);
+      tr_ip.Set(v.GetData(), dim);
       cfe.CalcVShape(tr_ip, coarse_shape);
 
       AddMult_a_AAt(w, fine_shape, fine_mass);
@@ -1407,7 +1406,7 @@ void VectorFiniteElement::LocalL2Projection_ND(
       const IntegrationPoint &ip = ir.IntPoint(i);
       this->CalcVShape(ip, fine_shape);
       Trans.Transform(ip, v);
-      tr_ip.Set(v, dim);
+      tr_ip.Set(v.GetData(), dim);
       cfe.CalcVShape(tr_ip, coarse_shape);
 
       AddMult_a_AAt(ip.weight, fine_shape, fine_mass);

fem/fe/fe_base.hpp

@@ -256,7 +256,7 @@ protected:
 
 public:
    /// Enumeration for range_type and deriv_range_type
-   enum RangeType { SCALAR, VECTOR };
+   enum RangeType { UNKNOWN_RANGE_TYPE = -1, SCALAR, VECTOR };
 
    /** @brief Enumeration for MapType: defines how reference functions are
        mapped to physical space.
@@ -270,6 +270,8 @@ public:
    */
    enum MapType
    {
+      UNKNOWN_MAP_TYPE = -1, /**< Used to distinguish an unset MapType variable
+                                  from the known values below. */
       VALUE,     /**< For scalar fields; preserves point values
                           \f$ u(x) = \hat u(\hat x) \f$ */
       INTEGRAL,  /**< For scalar fields; preserves volume integrals
@@ -1091,6 +1093,11 @@ public:
    // { CalcLegendre(p, x, u); }
    { CalcChebyshev(p, x, u); }
 
+   /** @brief Evaluate the values of a hierarchical 1D basis at point x
+       hierarchical = k-th basis function is degree k polynomial */
+   static void CalcBasis(const int p, const double x, Vector &u)
+   { CalcBasis(p, x, u.GetData()); }
+
    /// Evaluate the values and derivatives of a hierarchical 1D basis at point @a x
    static void CalcBasis(const int p, const double x, double *u, double *d)
    // { CalcMono(p, x, u, d); }
@@ -1098,6 +1105,11 @@ public:
    // { CalcLegendre(p, x, u, d); }
    { CalcChebyshev(p, x, u, d); }
 
+   /** @brief Evaluate the values and derivatives of a hierarchical 1D basis at
+       point @a x. */
+   static void CalcBasis(const int p, const double x, Vector &u, Vector &d)
+   { CalcBasis(p, x, u.GetData(), d.GetData()); }
+
    /// Evaluate the values, derivatives and second derivatives of a hierarchical 1D basis at point x
    static void CalcBasis(const int p, const double x, double *u, double *d,
                          double *dd)
@@ -1106,6 +1118,12 @@ public:
    // { CalcLegendre(p, x, u, d); }
    { CalcChebyshev(p, x, u, d, dd); }
 
+   /** @brief Evaluate the values, derivatives and second derivatives of a
+       hierarchical 1D basis at point @a x. */
+   static void CalcBasis(const int p, const double x, Vector &u, Vector &d,
+                         Vector &dd)
+   { CalcBasis(p, x, u.GetData(), d.GetData(), dd.GetData()); }
+
    /// Evaluate a representation of a Delta function at point x
    static double CalcDelta(const int p, const double x)
    { return pow(x, (double) p); }
@@ -1135,12 +1153,24 @@ public:
    static void CalcBernstein(const int p, const double x, double *u)
    { CalcBinomTerms(p, x, 1. - x, u); }
 
+   /** @brief Compute the values of the Bernstein basis functions of order
+       @a p at coordinate @a x and store the results in the already allocated
+       @a u array. */
+   static void CalcBernstein(const int p, const double x, Vector &u)
+   { CalcBernstein(p, x, u.GetData()); }
+
    /** @brief Compute the values and derivatives of the Bernstein basis functions
        of order @a p at coordinate @a x and store the results in the already allocated
        @a u and @a d arrays. */
    static void CalcBernstein(const int p, const double x, double *u, double *d)
    { CalcBinomTerms(p, x, 1. - x, u, d); }
 
+   /** @brief Compute the values and derivatives of the Bernstein basis
+       functions of order @a p at coordinate @a x and store the results in the
+       already allocated @a u and @a d arrays. */
+   static void CalcBernstein(const int p, const double x, Vector &u, Vector &d)
+   { CalcBernstein(p, x, u.GetData(), d.GetData()); }
+
    static void CalcLegendre(const int p, const double x, double *u);
    static void CalcLegendre(const int p, const double x, double *u, double *d);
 

fem/fe/fe_nd.cpp

@@ -319,6 +319,9 @@ void ND_HexahedronElement::CalcVShape(const IntegrationPoint &ip,
 
    if (obasis1d.IsIntegratedType())
    {
+#ifdef MFEM_THREAD_SAFE
+      Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
+#endif
       cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
       cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
       cbasis1d.Eval(ip.z, shape_cz, dshape_cz);
@@ -661,6 +664,9 @@ void ND_QuadrilateralElement::CalcVShape(const IntegrationPoint &ip,
 
    if (obasis1d.IsIntegratedType())
    {
+#ifdef MFEM_THREAD_SAFE
+      Vector dshape_cx(p + 1), dshape_cy(p + 1);
+#endif
       cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
       cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
       obasis1d.ScaleIntegrated(false);

fem/fe/fe_pos.cpp

@@ -212,7 +212,7 @@ void BiQuadPos2DFiniteElement::GetLocalInterpolation(
 void BiQuadPos2DFiniteElement::Project(
    Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
 {
-   double *d = dofs;
+   double *d = dofs.GetData();
 
    for (int i = 0; i < 9; i++)
    {
@@ -382,8 +382,8 @@ void H1Pos_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
    Vector shape_x(p+1), shape_y(p+1);
 #endif
 
-   Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData() );
-   Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData() );
+   Poly_1D::CalcBernstein(p, ip.x, shape_x);
+   Poly_1D::CalcBernstein(p, ip.y, shape_y);
 
    // Reorder so that vertices are at the beginning of the list
    for (int o = 0, j = 0; j <= p; j++)
@@ -402,8 +402,8 @@ void H1Pos_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
    Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
 #endif
 
-   Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData(), dshape_x.GetData() );
-   Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData(), dshape_y.GetData() );
+   Poly_1D::CalcBernstein(p, ip.x, shape_x, dshape_x);
+   Poly_1D::CalcBernstein(p, ip.y, shape_y, dshape_y);
 
    // Reorder so that vertices are at the beginning of the list
    for (int o = 0, j = 0; j <= p; j++)

fem/fe/fe_rt.cpp

@@ -150,6 +150,9 @@ void RT_QuadrilateralElement::CalcVShape(const IntegrationPoint &ip,
 
    if (obasis1d.IsIntegratedType())
    {
+#ifdef MFEM_THREAD_SAFE
+      Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
+#endif
       cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
       cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
       obasis1d.ScaleIntegrated(false);
@@ -478,6 +481,9 @@ void RT_HexahedronElement::CalcVShape(const IntegrationPoint &ip,
 
    if (obasis1d.IsIntegratedType())
    {
+#ifdef MFEM_THREAD_SAFE
+      Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1), dshape_cz(pp1 + 1);
+#endif
       cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
       cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
       cbasis1d.Eval(ip.z, shape_cz, dshape_cz);

fem/fe_coll.cpp

@@ -22,6 +22,70 @@ namespace mfem
 
 using namespace std;
 
+const FiniteElement *
+FiniteElementCollection::FiniteElementForDim(int dim) const
+{
+   for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
+   {
+      const FiniteElement *fe = FiniteElementForGeometry((Geometry::Type)g);
+      if (fe != NULL)
+      {
+         return fe;
+      }
+   }
+   return NULL;
+}
+
+int FiniteElementCollection::GetRangeType(int dim) const
+{
+   const FiniteElement *fe = FiniteElementForDim(dim);
+   if (fe)
+   {
+      return fe->GetRangeType();
+   }
+   return FiniteElement::UNKNOWN_RANGE_TYPE;
+}
+
+int FiniteElementCollection::GetDerivRangeType(int dim) const
+{
+   const FiniteElement *fe = FiniteElementForDim(dim);
+   if (fe)
+   {
+      return fe->GetDerivRangeType();
+   }
+   return FiniteElement::UNKNOWN_RANGE_TYPE;
+}
+
+int FiniteElementCollection::GetMapType(int dim) const
+{
+   const FiniteElement *fe = FiniteElementForDim(dim);
+   if (fe)
+   {
+      return fe->GetMapType();
+   }
+   return FiniteElement::UNKNOWN_MAP_TYPE;
+}
+
+int FiniteElementCollection::GetDerivType(int dim) const
+{
+   const FiniteElement *fe = FiniteElementForDim(dim);
+   if (fe)
+   {
+      return fe->GetDerivType();
+   }
+   return FiniteElement::NONE;
+}
+
+int FiniteElementCollection::GetDerivMapType(int dim) const
+{
+   const FiniteElement *fe = FiniteElementForDim(dim);
+   if (fe)
+   {
+      return fe->GetDerivMapType();
+   }
+   return FiniteElement::UNKNOWN_MAP_TYPE;
+}
+
 int FiniteElementCollection::HasFaceDofs(Geometry::Type geom, int p) const
 {
    switch (geom)

fem/fe_coll.hpp

@@ -51,6 +51,14 @@ public:
    virtual const FiniteElement *
    FiniteElementForGeometry(Geometry::Type GeomType) const = 0;
 
+   /** @brief Returns the first non-NULL FiniteElement for the given dimension
+
+       @note Repeatedly calls FiniteElementForGeometry in the order defined in
+       the Geometry::Type enumeration.
+   */
+   virtual const FiniteElement *
+   FiniteElementForDim(int dim) const;
+
    virtual int DofForGeometry(Geometry::Type GeomType) const = 0;
 
    /** @brief Returns an array, say p, that maps a local permuted index i to a
@@ -66,6 +74,17 @@ public:
 
    virtual int GetContType() const = 0;
 
+   /** @note The following methods provide the same information as the
+       corresponding methods of the FiniteElement base class.
+       @{
+   */
+   virtual int GetRangeType(int dim) const;
+   virtual int GetDerivRangeType(int dim) const;
+   virtual int GetMapType(int dim) const;
+   virtual int GetDerivType(int dim) const;
+   virtual int GetDerivMapType(int dim) const;
+   /** @} */
+
    int HasFaceDofs(Geometry::Type geom, int p) const;
 
    virtual const FiniteElement *TraceFiniteElementForGeometry(

fem/fmsconvert.cpp

@@ -1920,7 +1920,7 @@ DataCollectionToFmsDataCollection(DataCollection *mfem_dc,
                                              FMS_NODAL_GAUSS_CLOSED, 1);
       err |= FmsDataCollectionAddField(*dc, "Coords", &fcoords);
       err |= FmsFieldSet(fcoords, fdcoords, mmesh->SpaceDimension(), FMS_BY_NODES,
-                         FMS_DOUBLE, mverts);
+                         FMS_DOUBLE, mverts.HostRead());
       err |= FmsComponentSetCoordinates(volume, fcoords);
    }
 

fem/gridfunc.cpp

@@ -416,13 +416,24 @@ void GridFunction::GetNodalValues(int i, Array<double> &nval, int vdim) const
 
    if (FElem->GetRangeType() == FiniteElement::SCALAR)
    {
-      MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
-                  "invalid FE map type");
       Vector shape(dof);
-      for (k = 0; k < n; k++)
+      if (FElem->GetMapType() == FiniteElement::VALUE)
       {
-         FElem->CalcShape(ElemVert->IntPoint(k), shape);
-         nval[k] = shape * ((const double *)loc_data + dof * vdim);
+         for (k = 0; k < n; k++)
+         {
+            FElem->CalcShape(ElemVert->IntPoint(k), shape);
+            nval[k] = shape * (&loc_data[dof * vdim]);
+         }
+      }
+      else
+      {
+         ElementTransformation *Tr = fes->GetElementTransformation(i);
+         for (k = 0; k < n; k++)
+         {
+            Tr->SetIntPoint(&ElemVert->IntPoint(k));
+            FElem->CalcPhysShape(*Tr, shape);
+            nval[k] = shape * (&loc_data[dof * vdim]);
+         }
       }
    }
    else
@@ -495,7 +506,7 @@ void GridFunction::GetVectorValue(int i, const IntegrationPoint &ip,
       val.SetSize(vdim);
       for (int k = 0; k < vdim; k++)
       {
-         val(k) = shape * ((const double *)loc_data + dof * k);
+         val(k) = shape * (&loc_data[dof * k]);
       }
    }
    else
@@ -1027,7 +1038,7 @@ void GridFunction::GetVectorValue(ElementTransformation &T,
       val.SetSize(vdim);
       for (int k = 0; k < vdim; k++)
       {
-         val(k) = shape * ((const double *)loc_data + dof * k);
+         val(k) = shape * (&loc_data[dof * k]);
       }
    }
    else
@@ -1078,7 +1089,7 @@ void GridFunction::GetVectorValues(ElementTransformation &T,
 
          for (int k = 0; k < vdim; k++)
          {
-            vals(k,j) = shape * ((const double *)loc_data + dof * k);
+            vals(k,j) = shape * (&loc_data[dof * k]);
          }
       }
    }
@@ -1181,8 +1192,7 @@ void GridFunction::GetValuesFrom(const GridFunction &orig_func)
          orig_fe->CalcShape(ip, shape);
          for (d = 0; d < vdim; d++)
          {
-            loc_values(d*dof+j) =
-               shape * ((const double *)orig_loc_values + d * odof) ;
+            loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
          }
       }
       if (doftrans)
@@ -1224,8 +1234,7 @@ void GridFunction::GetBdrValuesFrom(const GridFunction &orig_func)
          orig_fe->CalcShape(ip, shape);
          for (d = 0; d < vdim; d++)
          {
-            loc_values(d*dof+j) =
-               shape * ((const double *)orig_loc_values + d * odof);
+            loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
          }
       }
       SetSubVector(vdofs, loc_values);
@@ -4099,16 +4108,16 @@ void TensorProductLegendre(int dim,                // input
    // Map x to [0, 1] to use CalcLegendre since it uses shifted Legendre Polynomials.
    double x1 = (x(0) - xmin(0))/(xmax(0)-xmin(0)), x2, x3;
    Vector poly_x(order+1), poly_y(order+1), poly_z(order+1);
-   poly1d.CalcLegendre(order, x1, poly_x);
+   poly1d.CalcLegendre(order, x1, poly_x.GetData());
    if (dim > 1)
    {
       x2 = (x(1)-xmin(1))/(xmax(1)-xmin(1));
-      poly1d.CalcLegendre(order, x2, poly_y);
+      poly1d.CalcLegendre(order, x2, poly_y.GetData());
    }
    if (dim == 3)
    {
       x3 = (x(2)-xmin(2))/(xmax(2)-xmin(2));
-      poly1d.CalcLegendre(order, x3, poly_z);
+      poly1d.CalcLegendre(order, x3, poly_z.GetData());
    }
 
    int basis_dimension = static_cast<int>(pow(order+1,dim));

fem/gslib.cpp

@@ -179,24 +179,27 @@ void FindPointsGSLIB::FindPoints(const Vector &point_pos,
    gsl_ref.SetSize(points_cnt * dim);
    gsl_dist.SetSize(points_cnt);
 
-   const double *xv_base[dim];
-   unsigned xv_stride[dim];
-   for (int d = 0; d < dim; d++)
+   auto xvFill = [&](const double *xv_base[], unsigned xv_stride[], int dim)
    {
-      if (point_pos_ordering == Ordering::byNODES)
+      for (int d = 0; d < dim; d++)
       {
-         xv_base[d] = point_pos.GetData() + d*points_cnt;
-         xv_stride[d] = sizeof(double);
+         if (point_pos_ordering == Ordering::byNODES)
+         {
+            xv_base[d] = point_pos.GetData() + d*points_cnt;
+            xv_stride[d] = sizeof(double);
+         }
+         else
+         {
+            xv_base[d] = point_pos.GetData() + d;
+            xv_stride[d] = dim*sizeof(double);
+         }
       }
-      else
-      {
-         xv_base[d] = point_pos.GetData() + d;
-         xv_stride[d] = dim*sizeof(double);
-      }
-   }
-
+   };
    if (dim == 2)
    {
+      const double *xv_base[2];
+      unsigned xv_stride[2];
+      xvFill(xv_base, xv_stride, dim);
       findpts_2(gsl_code.GetData(), sizeof(unsigned int),
                 gsl_proc.GetData(), sizeof(unsigned int),
                 gsl_elem.GetData(), sizeof(unsigned int),
@@ -204,8 +207,11 @@ void FindPointsGSLIB::FindPoints(const Vector &point_pos,
                 gsl_dist.GetData(), sizeof(double),
                 xv_base, xv_stride, points_cnt, fdata2D);
    }
-   else
+   else  // dim == 3
    {
+      const double *xv_base[3];
+      unsigned xv_stride[3];
+      xvFill(xv_base, xv_stride, dim);
       findpts_3(gsl_code.GetData(), sizeof(unsigned int),
                 gsl_proc.GetData(), sizeof(unsigned int),
                 gsl_elem.GetData(), sizeof(unsigned int),
@@ -1080,7 +1086,7 @@ void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
 
             sarray_transfer(struct send_pt, sendpt, proc, 1, cr);
             sdpt = (struct send_pt *)sendpt->ptr;
-            for (int index = 0; index < sendpt->n; index++)
+            for (int index = 0; index < static_cast<int>(sendpt->n); index++)
             {
                int idx = field_in.FESpace()->GetOrdering() == Ordering::byNODES ?
                          sdpt->index + j*nptorig :
@@ -1194,24 +1200,27 @@ void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
    gsl_ref.SetSize(points_cnt * dim);
    gsl_dist.SetSize(points_cnt);
 
-   const double *xv_base[dim];
-   unsigned xv_stride[dim];
-   for (int d = 0; d < dim; d++)
+   auto xvFill = [&](const double *xv_base[], unsigned xv_stride[], int dim)
    {
-      if (point_pos_ordering == Ordering::byNODES)
+      for (int d = 0; d < dim; d++)
       {
-         xv_base[d] = point_pos.GetData() + d*points_cnt;
-         xv_stride[d] = sizeof(double);
+         if (point_pos_ordering == Ordering::byNODES)
+         {
+            xv_base[d] = point_pos.GetData() + d*points_cnt;
+            xv_stride[d] = sizeof(double);
+         }
+         else
+         {
+            xv_base[d] = point_pos.GetData() + d;
+            xv_stride[d] = dim*sizeof(double);
+         }
       }
-      else
-      {
-         xv_base[d] = point_pos.GetData() + d;
-         xv_stride[d] = dim*sizeof(double);
-      }
-   }
-
+   };
    if (dim == 2)
    {
+      const double *xv_base[2];
+      unsigned xv_stride[2];
+      xvFill(xv_base, xv_stride, dim);
       findptsms_2(gsl_code.GetData(), sizeof(unsigned int),
                   gsl_proc.GetData(), sizeof(unsigned int),
                   gsl_elem.GetData(), sizeof(unsigned int),
@@ -1221,8 +1230,11 @@ void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
                   point_id.GetData(), sizeof(unsigned int), &match,
                   points_cnt, fdata2D);
    }
-   else
+   else  // dim == 3
    {
+      const double *xv_base[3];
+      unsigned xv_stride[3];
+      xvFill(xv_base, xv_stride, dim);
       findptsms_3(gsl_code.GetData(), sizeof(unsigned int),
                   gsl_proc.GetData(), sizeof(unsigned int),
                   gsl_elem.GetData(), sizeof(unsigned int),

fem/lininteg.cpp

@@ -1032,12 +1032,12 @@ void WhiteGaussianNoiseDomainLFIntegrator::AssembleRHSElementVect
       massinteg.AssembleElementMatrix(el, Tr, *M);
       CholeskyFactors chol(M->Data());
       chol.Factor(M->Height());
-      chol.LMult(n,1,elvect);
+      chol.LMult(n,1,elvect.GetData());
    }
    else
    {
       CholeskyFactors chol(L[iel]->Data());
-      chol.LMult(n,1,elvect);
+      chol.LMult(n,1,elvect.GetData());
    }
 }
 

fem/pfespace.cpp

@@ -2606,7 +2606,7 @@ int ParFiniteElementSpace
       if (dump < 10)
       {
          char fname[100];
-         sprintf(fname, "dofs%02d.txt", MyRank);
+         snprintf(fname, 100, "dofs%02d.txt", MyRank);
          std::ofstream f(fname);
          DebugDumpDOFs(f, deps, dof_group, dof_owner, finalized);
          dump++;

fem/pgridfunc.cpp

@@ -151,14 +151,14 @@ void ParGridFunction::ParallelAverage(Vector &tv) const
 {
    MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
    pfes->GetProlongationMatrix()->MultTranspose(*this, tv);
-   pfes->DivideByGroupSize(tv);
+   pfes->DivideByGroupSize(tv.HostReadWrite());
 }
 
 void ParGridFunction::ParallelAverage(HypreParVector &tv) const
 {
    MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
    pfes->GetProlongationMatrix()->MultTranspose(*this, tv);
-   pfes->DivideByGroupSize(tv);
+   pfes->DivideByGroupSize(tv.HostReadWrite());
 }
 
 HypreParVector *ParGridFunction::ParallelAverage() const
@@ -353,7 +353,7 @@ void ParGridFunction::GetVectorValue(int i, const IntegrationPoint &ip,
          val.SetSize(vdim);
          for (int k = 0; k < vdim; k++)
          {
-            val(k) = shape * ((const double *)loc_data + dof * k);
+            val(k) = shape * (&loc_data[dof * k]);
          }
       }
       else
@@ -468,7 +468,7 @@ void ParGridFunction::GetVectorValue(ElementTransformation &T,
       val.SetSize(vdim);
       for (int k = 0; k < vdim; k++)
       {
-         val(k) = shape * ((const double *)loc_data + dof * k);
+         val(k) = shape * (&loc_data[dof * k]);
       }
    }
    else
@@ -644,9 +644,9 @@ void ParGridFunction::ProjectBdrCoefficient(
 
    // Count the values globally.
    GroupCommunicator &gcomm = pfes->GroupComm();
-   gcomm.Reduce<int>(values_counter, GroupCommunicator::Sum);
+   gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
    // Accumulate the values globally.
-   gcomm.Reduce<double>(values, GroupCommunicator::Sum);
+   gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
    // Only the values in the master are guaranteed to be correct!
    for (int i = 0; i < values.Size(); i++)
    {
@@ -682,9 +682,9 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
 
    // Count the values globally.
    GroupCommunicator &gcomm = pfes->GroupComm();
-   gcomm.Reduce<int>(values_counter, GroupCommunicator::Sum);
+   gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
    // Accumulate the values globally.
-   gcomm.Reduce<double>(values, GroupCommunicator::Sum);
+   gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
    // Only the values in the master are guaranteed to be correct!
    for (int i = 0; i < values.Size(); i++)
    {
@@ -1109,11 +1109,11 @@ void ParGridFunction::ComputeFlux(
    SumFluxAndCount(blfi, flux, count, wcoef, subdomain);
 
    // Accumulate flux and counts in parallel
-   ffes->GroupComm().Reduce<double>(flux, GroupCommunicator::Sum);
-   ffes->GroupComm().Bcast<double>(flux);
+   ffes->GroupComm().Reduce<double>(flux.HostReadWrite(), GroupCommunicator::Sum);
+   ffes->GroupComm().Bcast<double>(flux.HostReadWrite());
 
-   ffes->GroupComm().Reduce<int>(count, GroupCommunicator::Sum);
-   ffes->GroupComm().Bcast<int>(count);
+   ffes->GroupComm().Reduce<int>(count.HostReadWrite(), GroupCommunicator::Sum);
+   ffes->GroupComm().Bcast<int>(count.HostReadWrite());
 
    // complete averaging
    for (int i = 0; i < count.Size(); i++)

fem/qfunction.cpp

@@ -44,31 +44,6 @@ QuadratureFunction::QuadratureFunction(Mesh *mesh, std::istream &in)
    Load(in, vdim*qspace->GetSize());
 }
 
-void QuadratureFunction::SetSpace(QuadratureSpaceBase *qspace_, int vdim_)
-{
-   if (qspace_ != qspace)
-   {
-      if (own_qspace) { delete qspace; }
-      qspace = qspace_;
-      own_qspace = false;
-   }
-   vdim = (vdim_ < 0) ? vdim : vdim_;
-   SetSize(vdim*qspace->GetSize());
-}
-
-void QuadratureFunction::SetSpace(
-   QuadratureSpaceBase *qspace_, double *qf_data, int vdim_)
-{
-   if (qspace_ != qspace)
-   {
-      if (own_qspace) { delete qspace; }
-      qspace = qspace_;
-      own_qspace = false;
-   }
-   vdim = (vdim_ < 0) ? vdim : vdim_;
-   NewDataAndSize(qf_data, vdim*qspace->GetSize());
-}
-
 void QuadratureFunction::Save(std::ostream &os) const
 {
    GetSpace()->Save(os);

fem/qfunction.hpp

@@ -262,6 +262,34 @@ inline void QuadratureFunction::GetValues(
    }
 }
 
+
+inline void QuadratureFunction::SetSpace(QuadratureSpaceBase *qspace_,
+                                         int vdim_)
+{
+   if (qspace_ != qspace)
+   {
+      if (own_qspace) { delete qspace; }
+      qspace = qspace_;
+      own_qspace = false;
+   }
+   vdim = (vdim_ < 0) ? vdim : vdim_;
+   SetSize(vdim*qspace->GetSize());
+}
+
+inline void QuadratureFunction::SetSpace(
+   QuadratureSpaceBase *qspace_, double *qf_data, int vdim_)
+{
+   if (qspace_ != qspace)
+   {
+      if (own_qspace) { delete qspace; }
+      qspace = qspace_;
+      own_qspace = false;
+   }
+   vdim = (vdim_ < 0) ? vdim : vdim_;
+   NewDataAndSize(qf_data, vdim*qspace->GetSize());
+}
+
+
 } // namespace mfem
 
 #endif

fem/qspace.cpp

@@ -16,7 +16,7 @@ namespace mfem
 
 QuadratureSpaceBase::QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
                                          const IntegrationRule &ir)
-   : mesh(mesh_)
+   : mesh(mesh_), order(ir.GetOrder())
 {
    for (int g = 0; g < Geometry::NumGeom; g++)
    {

fem/staticcond.cpp

@@ -349,7 +349,7 @@ void StaticCondensation::ReduceRHS(const Vector &b, Vector &sc_b) const
 
       LUFactors lu(const_cast<double*>((const double*)A_data) + A_offsets[i],
                    const_cast<int*>((const int*)A_ipiv) + A_ipiv_offsets[i]);
-      lu.LSolve(npd, 1, b_p);
+      lu.LSolve(npd, 1, b_p.GetData());
 
       if (symm)
       {
@@ -527,8 +527,9 @@ void StaticCondensation::ComputeSolution(
 
       LUFactors lu(const_cast<double*>((const double*)A_data) + A_offsets[i],
                    const_cast<int*>((const int*)A_ipiv) + A_ipiv_offsets[i]);
-      lu.LSolve(npd, 1, b_p);
-      lu.BlockBackSolve(npd, ned, 1, lu.data + npd*npd, s_e, b_p);
+      lu.LSolve(npd, 1, b_p.GetData());
+      lu.BlockBackSolve(npd, ned, 1, lu.data + npd*npd, s_e.GetData(),
+                        b_p.GetData());
 
       for (int j = 0; j < npd; j++)
       {

fem/tmop.hpp

@@ -1105,6 +1105,59 @@ public:
    virtual ~TMOP_QuadraticLimiter() { }
 };
 
+/// Exponential limiter function in TMOP_Integrator.
+class TMOP_ExponentialLimiter : public TMOP_LimiterFunction
+{
+public:
+   virtual double Eval(const Vector &x, const Vector &x0, double dist) const
+   {
+      MFEM_ASSERT(x.Size() == x0.Size(), "Bad input.");
+
+      return  exp(10.0*((x.DistanceSquaredTo(x0) / (dist * dist))-1.0));
+   }
+
+   virtual void Eval_d1(const Vector &x, const Vector &x0, double dist,
+                        Vector &d1) const
+   {
+      MFEM_ASSERT(x.Size() == x0.Size(), "Bad input.");
+
+      d1.SetSize(x.Size());
+      double dist_squared = dist*dist;
+      subtract(20.0*exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared) - 1.0)) /
+               dist_squared, x, x0, d1);
+   }
+
+   virtual void Eval_d2(const Vector &x, const Vector &x0, double dist,
+                        DenseMatrix &d2) const
+   {
+      MFEM_ASSERT(x.Size() == x0.Size(), "Bad input.");
+      Vector tmp;
+      tmp.SetSize(x.Size());
+      double dist_squared = dist*dist;
+      double dist_squared_squared = dist_squared*dist_squared;
+      double f = exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared)-1.0));
+
+      subtract(x,x0,tmp);
+      d2.SetSize(x.Size());
+      d2(0,0) = ((400.0*tmp(0)*tmp(0)*f)/dist_squared_squared)+(20.0*f/dist_squared);
+      d2(1,1) = ((400.0*tmp(1)*tmp(1)*f)/dist_squared_squared)+(20.0*f/dist_squared);
+      d2(0,1) = (400.0*tmp(0)*tmp(1)*f)/dist_squared_squared;
+      d2(1,0) = d2(0,1);
+
+      if (x.Size() == 3)
+      {
+         d2(0,2) = (400.0*tmp(0)*tmp(2)*f)/dist_squared_squared;
+         d2(1,2) = (400.0*tmp(1)*tmp(2)*f)/dist_squared_squared;
+         d2(2,0) = d2(0,2);
+         d2(2,1) = d2(1,2);
+         d2(2,2) = ((400.0*tmp(2)*tmp(2)*f)/dist_squared_squared)+(20.0*f/dist_squared);
+      }
+
+   }
+
+   virtual ~TMOP_ExponentialLimiter() { }
+};
+
 class FiniteElementCollection;
 class FiniteElementSpace;
 class ParFiniteElementSpace;

fem/tmop/tmop_pa.cpp

@@ -58,8 +58,9 @@ void TMOP_Integrator::AssemblePA_Limiting()
 
    MFEM_VERIFY(PA.enabled, "AssemblePA_Limiting but PA is not enabled!");
    MFEM_VERIFY(lim_func, "No TMOP_LimiterFunction specification!")
-   MFEM_VERIFY(dynamic_cast<TMOP_QuadraticLimiter*>(lim_func),
-               "Only TMOP_QuadraticLimiter is supported");
+   MFEM_VERIFY(dynamic_cast<TMOP_QuadraticLimiter*>(lim_func) ||
+               dynamic_cast<TMOP_ExponentialLimiter*>(lim_func),
+               "Only TMOP_QuadraticLimiter and TMOP_ExponentialLimiter are supported");
 
    const FiniteElementSpace *fes = PA.fes;
    const int NE = PA.ne;

fem/tmop/tmop_pa_h2s_c0.cpp

@@ -25,8 +25,12 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
                            const int NE,
                            const DenseTensor &j_,
                            const Array<double> &w_,
+                           const Array<double> &b_,
                            const Array<double> &bld_,
+                           const Vector &x0_,
+                           const Vector &x1_,
                            Vector &h0_,
+                           const bool exp_lim,
                            const int d1d,
                            const int q1d)
 {
@@ -41,8 +45,11 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
                    Reshape(c0_.Read(), Q1D, Q1D, NE);
    const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
-   const auto W = Reshape(w_.Read(), Q1D, Q1D);
+   const auto b = Reshape(b_.Read(), Q1D, D1D);
    const auto bld = Reshape(bld_.Read(), Q1D, D1D);
+   const auto W = Reshape(w_.Read(), Q1D, Q1D);
+   const auto X0 = Reshape(x0_.Read(), D1D, D1D, DIM, NE);
+   const auto X1 = Reshape(x1_.Read(), D1D, D1D, DIM, NE);
 
    auto H0 = Reshape(h0_.Write(), DIM, DIM, Q1D, Q1D, NE);
 
@@ -54,19 +61,37 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
       constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
 
+      MFEM_SHARED double B[MQ1*MD1];
       MFEM_SHARED double BLD[MQ1*MD1];
 
       MFEM_SHARED double XY[NBZ][MD1*MD1];
       MFEM_SHARED double DQ[NBZ][MD1*MQ1];
       MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
 
-      kernels::internal::LoadX<MD1,NBZ>(e,D1D,LD,XY);
+      MFEM_SHARED double XY0[2][NBZ][MD1*MD1];
+      MFEM_SHARED double DQ0[2][NBZ][MD1*MQ1];
+      MFEM_SHARED double QQ0[2][NBZ][MQ1*MQ1];
 
+      MFEM_SHARED double XY1[2][NBZ][MD1*MD1];
+      MFEM_SHARED double DQ1[2][NBZ][MD1*MQ1];
+      MFEM_SHARED double QQ1[2][NBZ][MQ1*MQ1];
+
+      kernels::internal::LoadX<MD1,NBZ>(e,D1D,LD,XY);
+      kernels::internal::LoadX<MD1,NBZ>(e,D1D,X0,XY0);
+      kernels::internal::LoadX<MD1,NBZ>(e,D1D,X1,XY1);
+
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
       kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
 
       kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,BLD,XY,DQ);
       kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,BLD,DQ,QQ);
 
+      kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY0,DQ0);
+      kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ0,QQ0);
+
+      kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY1,DQ1);
+      kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ1,QQ1);
+
       MFEM_FOREACH_THREAD(qy,y,Q1D)
       {
          MFEM_FOREACH_THREAD(qx,x,Q1D)
@@ -77,15 +102,37 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
             const double coeff0 = const_c0 ? C0(0,0,0) : C0(qx,qy,e);
             const double weight_m = weight * lim_normal * coeff0;
 
-            double D;
+            double D, p0[2], p1[2];
             kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ,D);
+            kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ0,p0);
+            kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ1,p1);
+
             const double dist = D; // GetValues, default comp set to 0
 
             // lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
-            // d2.Diag(1.0 / (dist * dist), x.Size());
-            const double c = 1.0 / (dist * dist);
             double grad_grad[4];
-            kernels::Diag<2>(c, grad_grad);
+
+            if (!exp_lim)
+            {
+               // d2.Diag(1.0 / (dist * dist), x.Size());
+               const double c = 1.0 / (dist * dist);
+               kernels::Diag<2>(c, grad_grad);
+            }
+            else
+            {
+               double tmp[2];
+               kernels::Subtract<2>(1.0, p1, p0, tmp);
+               double dsq = kernels::DistanceSquared<2>(p1,p0);
+               double dist_squared = dist*dist;
+               double dist_squared_squared = dist_squared*dist_squared;
+               double f = exp(10.0*((dsq / dist_squared)-1.0));
+               grad_grad[0] = ((400.0*tmp[0]*tmp[0]*f)/dist_squared_squared)+
+                              (20.0*f/dist_squared);
+               grad_grad[1] = (400.0*tmp[0]*tmp[1]*f)/dist_squared_squared;
+               grad_grad[2] = grad_grad[1];
+               grad_grad[3] = ((400.0*tmp[1]*tmp[1]*f)/dist_squared_squared)+
+                              (20.0*f/dist_squared);
+            }
             ConstDeviceMatrix gg(grad_grad,DIM,DIM);
 
             for (int i = 0; i < DIM; i++)
@@ -111,11 +158,17 @@ void TMOP_Integrator::AssembleGradPA_C0_2D(const Vector &X) const
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
    const Array<double> &W   = PA.ir->GetWeights();
+   const Array<double> &B   = PA.maps->B;
    const Array<double> &BLD = PA.maps_lim->B;
    const Vector &C0 = PA.C0;
+   const Vector &X0 = PA.X0;
    Vector &H0 = PA.H0;
 
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,BLD,H0);
+   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
+   const bool exp_lim = (el) ? true : false;
+
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,H0,
+                           exp_lim);
 }
 
 } // namespace mfem

fem/tmop/tmop_pa_h3s_c0.cpp

@@ -25,8 +25,12 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
                            const int NE,
                            const DenseTensor &j_,
                            const Array<double> &w_,
+                           const Array<double> &b_,
                            const Array<double> &bld_,
+                           const Vector &x0_,
+                           const Vector &x1_,
                            Vector &h0_,
+                           const bool exp_lim,
                            const int d1d,
                            const int q1d)
 {
@@ -40,8 +44,11 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
                    Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
    const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
-   const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
+   const auto b = Reshape(b_.Read(), Q1D, D1D);
    const auto bld = Reshape(bld_.Read(), Q1D, D1D);
+   const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
+   const auto X0 = Reshape(x0_.Read(), D1D, D1D, D1D, DIM, NE);
+   const auto X1 = Reshape(x1_.Read(), D1D, D1D, D1D, DIM, NE);
 
    auto H0 = Reshape(h0_.Write(), DIM, DIM, Q1D, Q1D, Q1D, NE);
 
@@ -54,6 +61,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
       constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
 
+      MFEM_SHARED double B[MQ1*MD1];
       MFEM_SHARED double sBLD[MQ1*MD1];
       kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
       ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
@@ -65,12 +73,34 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
       DeviceCube DQQ(sm0, MD1,MQ1,MQ1);
       DeviceCube QQQ(sm1, MQ1,MQ1,MQ1);
 
+      MFEM_SHARED double DDD0[3][MD1*MD1*MD1];
+      MFEM_SHARED double DDQ0[3][MD1*MD1*MQ1];
+      MFEM_SHARED double DQQ0[3][MD1*MQ1*MQ1];
+      MFEM_SHARED double QQQ0[3][MQ1*MQ1*MQ1];
+
+      MFEM_SHARED double DDD1[3][MD1*MD1*MD1];
+      MFEM_SHARED double DDQ1[3][MD1*MD1*MQ1];
+      MFEM_SHARED double DQQ1[3][MD1*MQ1*MQ1];
+      MFEM_SHARED double QQQ1[3][MQ1*MQ1*MQ1];
+
       kernels::internal::LoadX(e,D1D,LD,DDD);
+      kernels::internal::LoadX<MD1>(e,D1D,X0,DDD0);
+      kernels::internal::LoadX<MD1>(e,D1D,X1,DDD1);
+
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
 
       kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
       kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);
       kernels::internal::EvalZ(D1D,Q1D,BLD,DQQ,QQQ);
 
+      kernels::internal::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD0,DDQ0);
+      kernels::internal::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ0,DQQ0);
+      kernels::internal::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ0,QQQ0);
+
+      kernels::internal::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD1,DDQ1);
+      kernels::internal::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ1,DQQ1);
+      kernels::internal::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ1,QQQ1);
+
       MFEM_FOREACH_THREAD(qz,z,Q1D)
       {
          MFEM_FOREACH_THREAD(qy,y,Q1D)
@@ -83,15 +113,44 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
                const double coeff0 = const_c0 ? C0(0,0,0,0) : C0(qx,qy,qz,e);
                const double weight_m = weight * lim_normal * coeff0;
 
-               double D;
+               double D, p0[3], p1[3];
                kernels::internal::PullEval(qx,qy,qz,QQQ,D);
+               kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ0,p0);
+               kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ1,p1);
+
                const double dist = D; // GetValues, default comp set to 0
 
                // lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
-               // d2.Diag(1.0 / (dist * dist), x.Size());
-               const double c = 1.0 / (dist * dist);
+
                double grad_grad[9];
-               kernels::Diag<3>(c, grad_grad);
+
+               if (!exp_lim)
+               {
+                  // d2.Diag(1.0 / (dist * dist), x.Size());
+                  const double c = 1.0 / (dist * dist);
+                  kernels::Diag<3>(c, grad_grad);
+               }
+               else
+               {
+                  double tmp[3];
+                  kernels::Subtract<3>(1.0, p1, p0, tmp);
+                  double dsq = kernels::DistanceSquared<3>(p1,p0);
+                  double dist_squared = dist*dist;
+                  double dist_squared_squared = dist_squared*dist_squared;
+                  double f = exp(10.0*((dsq / dist_squared)-1.0));
+                  grad_grad[0] = ((400.0*tmp[0]*tmp[0]*f)/dist_squared_squared)+
+                                 (20.0*f/dist_squared);
+                  grad_grad[1] = (400.0*tmp[0]*tmp[1]*f)/dist_squared_squared;
+                  grad_grad[2] = (400.0*tmp[0]*tmp[2]*f)/dist_squared_squared;
+                  grad_grad[3] = grad_grad[1];
+                  grad_grad[4] = ((400.0*tmp[1]*tmp[1]*f)/dist_squared_squared)+
+                                 (20.0*f/dist_squared);
+                  grad_grad[5] = (400.0*tmp[1]*tmp[2]*f)/dist_squared_squared;
+                  grad_grad[6] = grad_grad[2];
+                  grad_grad[7] = grad_grad[5];
+                  grad_grad[8] = ((400.0*tmp[2]*tmp[2]*f)/dist_squared_squared)+
+                                 (20.0*f/dist_squared);
+               }
                ConstDeviceMatrix gg(grad_grad,DIM,DIM);
 
                for (int i = 0; i < DIM; i++)
@@ -117,11 +176,17 @@ void TMOP_Integrator::AssembleGradPA_C0_3D(const Vector &X) const
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
    const Array<double> &W = PA.ir->GetWeights();
+   const Array<double> &B   = PA.maps->B;
    const Array<double> &BLD = PA.maps_lim->B;
    const Vector &C0 = PA.C0;
+   const Vector &X0 = PA.X0;
    Vector &H0 = PA.H0;
 
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,BLD,H0);
+   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
+   const bool exp_lim = (el) ? true : false;
+
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,
+                           H0,exp_lim);
 }
 
 } // namespace mfem

fem/tmop/tmop_pa_p2_c0.cpp

@@ -30,6 +30,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
                            const Vector &x0_,
                            const Vector &x1_,
                            Vector &y_,
+                           const bool exp_lim,
                            const int d1d,
                            const int q1d)
 {
@@ -110,14 +111,35 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
             const double dist = ld; // GetValues, default comp set to 0
 
             double d1[2];
-            // Eval_d1
+            // Eval_d1 (Quadratic Limiter)
             // subtract(1.0 / (dist * dist), x, x0, d1);
             // z = a * (x - y)
             // grad = a * (x - x0)
-            const double a = 1.0 / (dist * dist);
+
+            // Eval_d1 (Exponential Limiter)
+            // double dist_squared = dist*dist;
+            // subtract(20.0*exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared) - 1.0)) /
+            // dist_squared, x, x0, d1);
+            // z = a * (x - y)
+            // grad = a * (x - x0)
+
+            double a = 0.0;
             const double w = weight * lim_normal * coeff0;
+            const double dist_squared = dist * dist;
+
+            if (!exp_lim)
+            {
+               a =  1.0 / dist_squared;
+            }
+            else
+            {
+               double dsq = kernels::DistanceSquared<2>(p1,p0) / dist_squared;
+               a = 20.0*exp(10.0*(dsq - 1.0))/dist_squared;
+            }
             kernels::Subtract<2>(w*a, p1, p0, d1);
             kernels::internal::PushEval<MQ1,NBZ>(Q1D,qx,qy,d1,QQ0);
+
+
          }
       }
       MFEM_SYNC_THREAD;
@@ -143,8 +165,11 @@ void TMOP_Integrator::AddMultPA_C0_2D(const Vector &X, Vector &Y) const
    MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
    const Vector &X0 = PA.X0;
    const Vector &C0 = PA.C0;
+   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
+   const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y);
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y,
+                           exp_lim);
 }
 
 } // namespace mfem

fem/tmop/tmop_pa_p3_c0.cpp

@@ -30,6 +30,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
                            const Vector &x0_,
                            const Vector &x1_,
                            Vector &y_,
+                           const bool exp_lim,
                            const int d1d,
                            const int q1d)
 {
@@ -117,13 +118,32 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
                kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ1,p1);
 
                double d1[3];
-               // Eval_d1
+               // Eval_d1 (Quadratic Limiter)
                // subtract(1.0 / (dist * dist), x, x0, d1);
                // z = a * (x - y)
                // grad = a * (x - x0)
+
+               // Eval_d1 (Exponential Limiter)
+               // double dist_squared = dist*dist;
+               // subtract(20.0*exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared) - 1.0)) /
+               // dist_squared, x, x0, d1);
+               // z = a * (x - y)
+               // grad = a * (x - x0)
                const double dist = D; // GetValues, default comp set to 0
-               const double a = 1.0 / (dist * dist);
+               double a = 0.0;
                const double w = weight * lim_normal * coeff0;
+               const double dist_squared = dist * dist;
+
+               if (!exp_lim)
+               {
+                  a =  1.0 / dist_squared;
+               }
+               else
+               {
+                  double dsq = kernels::DistanceSquared<3>(p1,p0) / dist_squared;
+                  a = 20.0*exp(10.0*(dsq - 1.0))/dist_squared;
+               }
+
                kernels::Subtract<3>(w*a, p1, p0, d1);
                kernels::internal::PushEval<MQ1>(Q1D,qx,qy,qz,d1,QQQ0);
             }
@@ -153,8 +173,11 @@ void TMOP_Integrator::AddMultPA_C0_3D(const Vector &X, Vector &Y) const
    MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
    const Vector &X0 = PA.X0;
    const Vector &C0 = PA.C0;
+   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
+   const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y);
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y,
+                           exp_lim);
 }
 
 } // namespace mfem

fem/tmop/tmop_pa_w2_c0.cpp

@@ -31,6 +31,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
                            const Vector &x1_,
                            const Vector &ones,
                            Vector &energy,
+                           const bool exp_lim,
                            const int d1d,
                            const int q1d)
 {
@@ -107,9 +108,20 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
             kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ0,p0);
             kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ1,p1);
             const double dist = ld; // GetValues, default comp set to 0
-            const double id2 = 0.5 / (dist*dist);
-            const double dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
-            E(qx,qy,e) = weight * lim_normal * dsq * coeff0;
+            double id2 = 0.0;
+            double dsq = 0.0;
+            if (!exp_lim)
+            {
+               id2 = 0.5 / (dist*dist);
+               dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
+               E(qx,qy,e) = weight * lim_normal * dsq * coeff0;
+            }
+            else
+            {
+               id2 = 1.0 / (dist*dist);
+               dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
+               E(qx,qy,e) = weight * lim_normal * exp(10.0*(dsq-1.0)) * coeff0;
+            }
          }
       }
    });
@@ -135,7 +147,11 @@ double TMOP_Integrator::GetLocalStateEnergyPA_C0_2D(const Vector &X) const
    const Vector &O = PA.O;
    Vector &E = PA.E;
 
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E);
+   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
+   const bool exp_lim = (el) ? true : false;
+
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E,
+                           exp_lim);
 }
 
 } // namespace mfem

fem/tmop/tmop_pa_w3_c0.cpp

@@ -31,6 +31,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
                            const Vector &x1_,
                            const Vector &ones,
                            Vector &energy,
+                           const bool exp_lim,
                            const int d1d,
                            const int q1d)
 {
@@ -118,10 +119,20 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
                kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ1,p1);
 
                const double dist = D; // GetValues, default comp set to 0
-               const double id2 = 0.5 / (dist*dist);
-
-               const double dsq = kernels::DistanceSquared<3>(p1,p0) * id2;
-               E(qx,qy,qz,e) = weight * lim_normal * dsq * coeff0;
+               double id2 = 0.0;
+               double dsq = 0.0;
+               if (!exp_lim)
+               {
+                  id2 = 0.5 / (dist*dist);
+                  dsq = kernels::DistanceSquared<3>(p1,p0) * id2;
+                  E(qx,qy,qz,e) = weight * lim_normal * dsq * coeff0;
+               }
+               else
+               {
+                  id2 = 1.0 / (dist*dist);
+                  dsq = kernels::DistanceSquared<3>(p1,p0) * id2;
+                  E(qx,qy,qz,e) = weight * lim_normal * exp(10.0*(dsq-1.0)) * coeff0;
+               }
             }
          }
       }
@@ -148,7 +159,11 @@ double TMOP_Integrator::GetLocalStateEnergyPA_C0_3D(const Vector &X) const
    const Vector &O = PA.O;
    Vector &E = PA.E;
 
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E);
+   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
+   const bool exp_lim = (el) ? true : false;
+
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E,
+                           exp_lim);
 }
 
 } // namespace mfem

general/error.hpp

@@ -61,6 +61,10 @@ void mfem_error(const char *msg = NULL);
 /// Function called by the macro MFEM_WARNING.
 void mfem_warning(const char *msg = NULL);
 
+#ifdef MFEM_USE_ENZYME
+static void* __enzyme_inactive_global_err = (void*)mfem_error;
+static void* __enzyme_inactive_global_warn = (void*)mfem_warning;
+#endif
 }
 
 #ifndef _MFEM_FUNC_NAME

general/mem_manager.cpp

@@ -310,10 +310,11 @@ inline uintptr_t MmuLengthP(const void *ptr, const size_t bytes)
 /// The protected access error, used for the host
 static void MmuError(int, siginfo_t *si, void*)
 {
+   constexpr size_t buf_size = 64;
    fflush(0);
-   char str[64];
+   char str[buf_size];
    const void *ptr = si->si_addr;
-   sprintf(str, "Error while accessing address %p!", ptr);
+   snprintf(str, buf_size, "Error while accessing address %p!", ptr);
    mfem::out << std::endl << "An illegal memory access was made!";
    MFEM_ABORT(str);
 }
@@ -907,7 +908,7 @@ void MemoryManager::SetDeviceMemoryType_(void *h_ptr, unsigned flags,
    }
 }
 
-MemoryType MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
+void MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
 {
    const bool alias = flags & Mem::ALIAS;
    const bool registered = flags & Mem::REGISTERED;
@@ -924,7 +925,7 @@ MemoryType MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
    MFEM_ASSERT(registered || !(owns_host || owns_device || owns_internal) ||
                (!(owns_device || owns_internal) && h_ptr == nullptr),
                "invalid Memory state");
-   if (!mm.exists || !registered) { return h_mt; }
+   if (!mm.exists || !registered) { return; }
    if (alias)
    {
       if (owns_internal)
@@ -945,7 +946,6 @@ MemoryType MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
          mm.Erase(h_ptr, owns_device);
       }
    }
-   return h_mt;
 }
 
 void MemoryManager::DeleteDevice_(void *h_ptr, unsigned & flags)

general/mem_manager.hpp

@@ -644,9 +644,8 @@ private: // Static methods used by the Memory<T> class
    static void SetDeviceMemoryType_(void *h_ptr, unsigned flags,
                                     MemoryType d_mt);
 
-   /// Un-register and free memory identified by its host pointer. Returns the
-   /// memory type of the host pointer.
-   static MemoryType Delete_(void *h_ptr, MemoryType mt, unsigned flags);
+   /// Un-register and free memory identified by its host pointer.
+   static void Delete_(void *h_ptr, MemoryType mt, unsigned flags);
 
    /// Free device memory identified by its host pointer
    static void DeleteDevice_(void *h_ptr, unsigned & flags);
@@ -828,6 +827,23 @@ public:
 
    static MemoryType GetHostMemoryType() { return host_mem_type; }
    static MemoryType GetDeviceMemoryType() { return device_mem_type; }
+
+#ifdef MFEM_USE_ENZYME
+   static void myfree(void* mem, MemoryType MT, unsigned &flags)
+   {
+      MemoryManager::Delete_(mem, MT, flags);
+   }
+   __attribute__((used))
+   inline static void* __enzyme_allocation_like1[4] = {(void*)static_cast<void*(*)(void*, size_t, MemoryType, unsigned&)>(MemoryManager::New_),
+                                                       (void*)1, (void*)"-1,2,3", (void*)myfree
+                                                      };
+   __attribute__((used))
+   inline static void* __enzyme_allocation_like2[4] = {(void*)static_cast<void*(*)(void*, size_t, MemoryType, MemoryType, unsigned, unsigned&)>(MemoryManager::New_),
+                                                       (void*)1, (void*)"-1,2,4", (void*)MemoryManager::Delete_
+                                                      };
+   __attribute__((used))
+   inline static void* __enzyme_function_like[2] = {(void*)MemoryManager::Delete_, (void*)"free"};
+#endif
 };
 
 
@@ -1007,8 +1023,12 @@ inline void Memory<T>::Delete()
    const bool mt_host = h_mt == MemoryType::HOST;
    const bool std_delete = !registered && mt_host;
 
-   if (std_delete ||
-       MemoryManager::Delete_((void*)h_ptr, h_mt, flags) == MemoryType::HOST)
+   if (!std_delete)
+   {
+      MemoryManager::Delete_((void*)h_ptr, h_mt, flags);
+   }
+
+   if (mt_host)
    {
       if (flags & OWNS_HOST) { delete [] h_ptr; }
    }

general/version.cpp

@@ -112,9 +112,6 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_MEMALLOC
       "MFEM_USE_MEMALLOC\n"
 #endif
-#ifdef MFEM_USE_MESQUITE
-      "MFEM_USE_MESQUITE\n"
-#endif
 #ifdef MFEM_USE_METIS
       "MFEM_USE_METIS\n"
 #endif

general/zstr.hpp

@@ -60,9 +60,9 @@ namespace strict_fstream
 {
 
 // Overloaded error checks to handle POSIX and GNU strerror_r
-inline char* check_strerror_r(int r, char* buff, int err)
+inline char* check_strerror_r(int r, char* buff, size_t buff_size)
 {
-   if (r) { sprintf(buff, "unknown error: %d", err); }
+   if (r) { snprintf(buff, buff_size, "unknown error: %d", r); }
    return buff;
 }
 

linalg/constraints.cpp

@@ -51,14 +51,14 @@ Eliminator::Eliminator(const SparseMatrix& B, const Array<int>& lagrange_tdofs_,
 void Eliminator::Eliminate(const Vector& vin, Vector& vout) const
 {
    Bp.Mult(vin, vout);
-   Bsinverse.Solve(Bs.Height(), 1, vout);
+   Bsinverse.Solve(Bs.Height(), 1, vout.GetData());
    vout *= -1.0;
 }
 
 void Eliminator::EliminateTranspose(const Vector& vin, Vector& vout) const
 {
    Vector work(vin);
-   BsTinverse.Solve(Bs.Height(), 1, work);
+   BsTinverse.Solve(Bs.Height(), 1, work.GetData());
    Bp.MultTranspose(work, vout);
    vout *= -1.0;
 }
@@ -66,14 +66,14 @@ void Eliminator::EliminateTranspose(const Vector& vin, Vector& vout) const
 void Eliminator::LagrangeSecondary(const Vector& vin, Vector& vout) const
 {
    vout = vin;
-   Bsinverse.Solve(Bs.Height(), 1, vout);
+   Bsinverse.Solve(Bs.Height(), 1, vout.GetData());
 }
 
 void Eliminator::LagrangeSecondaryTranspose(const Vector& vin,
                                             Vector& vout) const
 {
    vout = vin;
-   BsTinverse.Solve(Bs.Height(), 1, vout);
+   BsTinverse.Solve(Bs.Height(), 1, vout.GetData());
 }
 
 void Eliminator::ExplicitAssembly(DenseMatrix& mat) const

linalg/cpardiso.cpp

@@ -55,7 +55,9 @@ void CPardisoSolver::SetOperator(const Operator &op)
    auto parcsr_op = static_cast<hypre_ParCSRMatrix *>(
                        const_cast<HypreParMatrix &>(hypreParMat));
 
+   hypreParMat.HostRead();
    hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+   hypreParMat.HypreRead();
 #if MFEM_HYPRE_VERSION >= 21600
    hypre_CSRMatrixBigJtoJ(csr_op);
 #endif

linalg/densemat.cpp

@@ -175,12 +175,26 @@ void DenseMatrix::Mult(const double *x, double *y) const
    kernels::Mult(height, width, Data(), x, y);
 }
 
+void DenseMatrix::Mult(const double *x, Vector &y) const
+{
+   MFEM_ASSERT(height == y.Size(), "incompatible dimensions");
+
+   Mult(x, y.GetData());
+}
+
+void DenseMatrix::Mult(const Vector &x, double *y) const
+{
+   MFEM_ASSERT(width == x.Size(), "incompatible dimensions");
+
+   Mult(x.GetData(), y);
+}
+
 void DenseMatrix::Mult(const Vector &x, Vector &y) const
 {
    MFEM_ASSERT(height == y.Size() && width == x.Size(),
                "incompatible dimensions");
 
-   Mult((const double *)x, (double *)y);
+   Mult(x.GetData(), y.GetData());
 }
 
 double DenseMatrix::operator *(const DenseMatrix &m) const
@@ -213,12 +227,26 @@ void DenseMatrix::MultTranspose(const double *x, double *y) const
    }
 }
 
+void DenseMatrix::MultTranspose(const double *x, Vector &y) const
+{
+   MFEM_ASSERT(width == y.Size(), "incompatible dimensions");
+
+   MultTranspose(x, y.GetData());
+}
+
+void DenseMatrix::MultTranspose(const Vector &x, double *y) const
+{
+   MFEM_ASSERT(height == x.Size(), "incompatible dimensions");
+
+   MultTranspose(x.GetData(), y);
+}
+
 void DenseMatrix::MultTranspose(const Vector &x, Vector &y) const
 {
    MFEM_ASSERT(height == x.Size() && width == y.Size(),
                "incompatible dimensions");
 
-   MultTranspose((const double *)x, (double *)y);
+   MultTranspose(x.GetData(), y.GetData());
 }
 
 void DenseMatrix::AddMult(const Vector &x, Vector &y) const
@@ -226,8 +254,8 @@ void DenseMatrix::AddMult(const Vector &x, Vector &y) const
    MFEM_ASSERT(height == y.Size() && width == x.Size(),
                "incompatible dimensions");
 
-   const double *xp = x, *d_col = data;
-   double *yp = y;
+   const double *xp = x.GetData(), *d_col = data;
+   double *yp = y.GetData();
    for (int col = 0; col < width; col++)
    {
       double x_col = xp[col];
@@ -262,8 +290,8 @@ void DenseMatrix::AddMult_a(double a, const Vector &x, Vector &y) const
    MFEM_ASSERT(height == y.Size() && width == x.Size(),
                "incompatible dimensions");
 
-   const double *xp = x, *d_col = data;
-   double *yp = y;
+   const double *xp = x.GetData(), *d_col = data;
+   double *yp = y.GetData();
    for (int col = 0; col < width; col++)
    {
       const double x_col = a*xp[col];
@@ -1181,7 +1209,7 @@ void DenseMatrix::SingularValues(Vector &sv) const
    int         n            = Width();
    double      *a           = copy_of_this.data;
    sv.SetSize(min(m, n));
-   double      *s           = sv;
+   double      *s           = sv.GetData();
    double      *u           = NULL;
    double      *vt          = NULL;
    double      *work        = NULL;
@@ -1459,10 +1487,10 @@ void DenseMatrix::GradToCurl(DenseMatrix &curl)
          int j = i+n;
 
          // curl of (Ui,0)
-         curl(i,0) = y;
+         curl(i,0) = -y;
 
          // curl of (0,Ui)
-         curl(j,0) = -x;
+         curl(j,0) = x;
       }
    }
    else
@@ -1495,6 +1523,20 @@ void DenseMatrix::GradToCurl(DenseMatrix &curl)
    }
 }
 
+void DenseMatrix::GradToVectorCurl2D(DenseMatrix &curl)
+{
+   MFEM_VERIFY(Width() == 2,
+               "DenseMatrix::GradToVectorCurl2D(...): dimension must be 2")
+
+   int n = Height();
+   // rotate gradient
+   for (int i = 0; i < n; i++)
+   {
+      curl(i,0) = (*this)(i,1);
+      curl(i,1) = -(*this)(i,0);
+   }
+}
+
 void DenseMatrix::GradToDiv(Vector &div)
 {
    MFEM_ASSERT(Width()*Height() == div.Size(), "incompatible Vector 'div'!");
@@ -4176,10 +4218,10 @@ const
 {
    int n = SizeI(), ne = SizeK();
    const int *I = elem_dof.GetI(), *J = elem_dof.GetJ(), *dofs;
-   const double *d_col = tdata;
+   const double *d_col = mfem::HostRead(tdata, n*SizeJ()*ne);
    double *yp = y.HostReadWrite();
    double x_col;
-   const double *xp = x;
+   const double *xp = x.HostRead();
    // the '4' here can be tuned for given platform and compiler
    if (n <= 4)
    {

linalg/densemat.hpp

@@ -141,12 +141,24 @@ public:
    /// Matrix vector multiplication.
    void Mult(const double *x, double *y) const;
 
+   /// Matrix vector multiplication.
+   void Mult(const double *x, Vector &y) const;
+
+   /// Matrix vector multiplication.
+   void Mult(const Vector &x, double *y) const;
+
    /// Matrix vector multiplication.
    virtual void Mult(const Vector &x, Vector &y) const;
 
    /// Multiply a vector with the transpose matrix.
    void MultTranspose(const double *x, double *y) const;
 
+   /// Multiply a vector with the transpose matrix.
+   void MultTranspose(const double *x, Vector &y) const;
+
+   /// Multiply a vector with the transpose matrix.
+   void MultTranspose(const Vector &x, double *y) const;
+
    /// Multiply a vector with the transpose matrix.
    virtual void MultTranspose(const Vector &x, Vector &y) const;
 
@@ -180,7 +192,7 @@ public:
 
    /// Compute y^t A x
    double InnerProduct(const Vector &x, const Vector &y) const
-   { return InnerProduct((const double *)x, (const double *)y); }
+   { return InnerProduct(x.GetData(), y.GetData()); }
 
    /// Returns a pointer to the inverse matrix.
    virtual MatrixInverse *Inverse() const;
@@ -236,6 +248,13 @@ public:
    /// Take the 2-norm of the columns of A and store in v
    void Norm2(double *v) const;
 
+   /// Take the 2-norm of the columns of A and store in v
+   void Norm2(Vector &v) const
+   {
+      MFEM_ASSERT(v.Size() == Width(), "incompatible Vector size!");
+      Norm2(v.GetData());
+   }
+
    /// Compute the norm ||A|| = max_{ij} |A_{ij}|
    double MaxMaxNorm() const;
 
@@ -325,8 +344,13 @@ public:
    /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
        CurlShape matrix. If *this is a N by D matrix, then curl is a D*N by
        D*(D-1)/2 matrix. The size of curl must be set outside. The dimension D
-       can be either 2 or 3. */
+       can be either 2 or 3. In 2D this computes the scalar-valued curl of a
+       2D vector field */
    void GradToCurl(DenseMatrix &curl);
+   /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
+       CurlShape matrix. This computes the vector-valued curl of a scalar field.
+       *this is N by 2 matrix and curl is N by 2 matrix as well. */
+   void GradToVectorCurl2D(DenseMatrix &curl);
    /** Given a DShape matrix (from a scalar FE), stored in *this,
        returns the DivShape vector. If *this is a N by dim matrix,
        then div is a dim*N vector. The size of div must be set

linalg/hiop.hpp

@@ -92,6 +92,8 @@ public:
     *  internally in HiOp. */
    virtual bool get_starting_point(const hiop::size_type &n, double *x0);
 
+   using hiop::hiopInterfaceBase::get_starting_point;
+
    virtual bool get_vars_info(const hiop::size_type &n, double *xlow, double* xupp,
                               NonlinearityType* type);
 
@@ -136,6 +138,8 @@ public:
                           const hiop::index_type *idx_cons,
                           const double *x, bool new_x, double *cons);
 
+   using hiop::hiopInterfaceBase::eval_cons;
+
    /** Evaluates the Jacobian of the subset of constraints indicated by
     *  idx_cons. The idx_cons is assumed to be of size num_cons.
     *  Example: if cons[c] = C(x)[idx_cons[c]] where c = 0 .. num_cons-1, then
@@ -152,6 +156,8 @@ public:
                               const hiop::index_type *idx_cons,
                               const double *x, bool new_x, double *Jac);
 
+   using hiop::hiopInterfaceDenseConstraints::eval_Jac_cons;
+
    /** Specifies column partitioning for distributed memory vectors.
     *  Process p owns vector entries with indices cols[p] to cols[p+1]-1,
     *  where p = 0 .. nranks-1. The cols array is of size nranks + 1.

linalg/mumps.cpp

@@ -42,7 +42,9 @@ void MUMPSSolver::SetOperator(const Operator &op)
 
    auto parcsr_op = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>(*APtr);
 
+   APtr->HostRead();
    hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+   APtr->HypreRead();
 #if MFEM_HYPRE_VERSION >= 21600
    hypre_CSRMatrixBigJtoJ(csr_op);
 #endif

linalg/petsc.cpp

@@ -1335,7 +1335,12 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
       PetscBool ismatis;
 #endif
 
+#if PETSC_VERSION_LT(3,18,0)
       ierr = PetscObjectTypeCompare((PetscObject)(pA->A),MATTRANSPOSEMAT,&istrans);
+#else
+      ierr = PetscObjectTypeCompare((PetscObject)(pA->A),MATTRANSPOSEVIRTUAL,
+                                    &istrans);
+#endif
       CCHKERRQ(pA->GetComm(),ierr);
       if (!istrans)
       {
@@ -2850,37 +2855,43 @@ void PetscBCHandler::ZeroBC(const Vector &x, Vector &y)
 // PetscLinearSolver methods
 
 PetscLinearSolver::PetscLinearSolver(MPI_Comm comm, const std::string &prefix,
-                                     bool wrapin)
-   : PetscSolver(), Solver(), wrap(wrapin)
+                                     bool wrapin, bool iter_mode)
+   : PetscSolver(), Solver(0,iter_mode), wrap(wrapin)
 {
    KSP ksp;
    ierr = KSPCreate(comm,&ksp); CCHKERRQ(comm,ierr);
    obj  = (PetscObject)ksp;
    ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
    ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
+   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
+   PCHKERRQ(ksp, ierr);
 }
 
 PetscLinearSolver::PetscLinearSolver(const PetscParMatrix &A,
-                                     const std::string &prefix)
-   : PetscSolver(), Solver(), wrap(false)
+                                     const std::string &prefix, bool iter_mode)
+   : PetscSolver(), Solver(0,iter_mode), wrap(false)
 {
    KSP ksp;
    ierr = KSPCreate(A.GetComm(),&ksp); CCHKERRQ(A.GetComm(),ierr);
    obj  = (PetscObject)ksp;
    ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
    ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
+   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
+   PCHKERRQ(ksp, ierr);
    SetOperator(A);
 }
 
 PetscLinearSolver::PetscLinearSolver(const HypreParMatrix &A, bool wrapin,
-                                     const std::string &prefix)
-   : PetscSolver(), Solver(), wrap(wrapin)
+                                     const std::string &prefix, bool iter_mode)
+   : PetscSolver(), Solver(0,iter_mode), wrap(wrapin)
 {
    KSP ksp;
    ierr = KSPCreate(A.GetComm(),&ksp); CCHKERRQ(A.GetComm(),ierr);
    obj  = (PetscObject)ksp;
    ierr = PetscObjectGetClassId(obj, &cid); PCHKERRQ(obj, ierr);
    ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
+   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
+   PCHKERRQ(ksp, ierr);
    SetOperator(A);
 }
 
@@ -3089,12 +3100,12 @@ void PetscLinearSolver::MultKernel(const Vector &b, Vector &x, bool trans) const
       }
    }
    B->PlaceMemory(b.GetMemory());
-   X->PlaceMemory(x.GetMemory(),iterative_mode);
 
    Customize();
 
-   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
-   PCHKERRQ(ksp, ierr);
+   PetscBool flg;
+   ierr = KSPGetInitialGuessNonzero(ksp, &flg);
+   X->PlaceMemory(x.GetMemory(),flg);
 
    // Solve the system.
    if (trans)
@@ -3129,8 +3140,9 @@ PetscLinearSolver::~PetscLinearSolver()
 
 // PetscPCGSolver methods
 
-PetscPCGSolver::PetscPCGSolver(MPI_Comm comm, const std::string &prefix)
-   : PetscLinearSolver(comm,prefix)
+PetscPCGSolver::PetscPCGSolver(MPI_Comm comm, const std::string &prefix,
+                               bool iter_mode)
+   : PetscLinearSolver(comm,prefix,iter_mode)
 {
    KSP ksp = (KSP)obj;
    ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
@@ -3138,8 +3150,9 @@ PetscPCGSolver::PetscPCGSolver(MPI_Comm comm, const std::string &prefix)
    ierr = KSPSetNormType(ksp,KSP_NORM_NATURAL); PCHKERRQ(ksp,ierr);
 }
 
-PetscPCGSolver::PetscPCGSolver(PetscParMatrix& A, const std::string &prefix)
-   : PetscLinearSolver(A,prefix)
+PetscPCGSolver::PetscPCGSolver(PetscParMatrix& A, const std::string &prefix,
+                               bool iter_mode)
+   : PetscLinearSolver(A,prefix,iter_mode)
 {
    KSP ksp = (KSP)obj;
    ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
@@ -3148,8 +3161,8 @@ PetscPCGSolver::PetscPCGSolver(PetscParMatrix& A, const std::string &prefix)
 }
 
 PetscPCGSolver::PetscPCGSolver(HypreParMatrix& A, bool wrap,
-                               const std::string &prefix)
-   : PetscLinearSolver(A,wrap,prefix)
+                               const std::string &prefix, bool iter_mode)
+   : PetscLinearSolver(A,wrap,prefix,iter_mode)
 {
    KSP ksp = (KSP)obj;
    ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);

linalg/petsc.hpp

@@ -744,16 +744,16 @@ private:
 
 public:
    PetscLinearSolver(MPI_Comm comm, const std::string &prefix = std::string(),
-                     bool wrap = true);
+                     bool wrap = true, bool iter_mode = false);
    PetscLinearSolver(const PetscParMatrix &A,
-                     const std::string &prefix = std::string());
+                     const std::string &prefix = std::string(), bool iter_mode = false);
    /// Constructs a solver using a HypreParMatrix.
    /** If @a wrap is true, then the MatMult ops of HypreParMatrix are wrapped.
        No preconditioner can be automatically constructed from PETSc. If
        @a wrap is false, the HypreParMatrix is converted into a the AIJ
        PETSc format, which is suitable for most preconditioning methods. */
    PetscLinearSolver(const HypreParMatrix &A, bool wrap = true,
-                     const std::string &prefix = std::string());
+                     const std::string &prefix = std::string(), bool iter_mode = false);
    virtual ~PetscLinearSolver();
 
    /// Sets the operator to be used for mat-vec operations and
@@ -780,10 +780,12 @@ public:
 class PetscPCGSolver : public PetscLinearSolver
 {
 public:
-   PetscPCGSolver(MPI_Comm comm, const std::string &prefix = std::string());
-   PetscPCGSolver(PetscParMatrix &A, const std::string &prefix = std::string());
-   PetscPCGSolver(HypreParMatrix &A,bool wrap=true,
-                  const std::string &prefix = std::string());
+   PetscPCGSolver(MPI_Comm comm, const std::string &prefix = std::string(),
+                  bool iter_mode = false);
+   PetscPCGSolver(PetscParMatrix &A, const std::string &prefix = std::string(),
+                  bool iter_mode = false);
+   PetscPCGSolver(HypreParMatrix &A, bool wrap = true,
+                  const std::string &prefix = std::string(), bool iter_mode = false);
 };
 
 

linalg/solvers.cpp

@@ -3034,7 +3034,7 @@ void BlockILU::Mult(const Vector &b, Vector &x) const
       }
       LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
       // x_i = D_ii^{-1} x_i
-      A_ii_inv.Solve(block_size, 1, xi);
+      A_ii_inv.Solve(block_size, 1, xi.GetData());
    }
 }
 
@@ -3197,7 +3197,8 @@ void UMFPackSolver::Mult(const Vector &b, Vector &x) const
    {
       int status =
          umfpack_di_solve(UMFPACK_At, mat->HostReadI(), mat->HostReadJ(),
-                          mat->HostReadData(), x, b, Numeric, Control, Info);
+                          mat->HostReadData(), x.HostWrite(), b.HostRead(),
+                          Numeric, Control, Info);
       umfpack_di_report_info(Control, Info);
       if (status < 0)
       {
@@ -3208,8 +3209,9 @@ void UMFPackSolver::Mult(const Vector &b, Vector &x) const
    else
    {
       SuiteSparse_long status =
-         umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->HostReadData(), x, b,
-                          Numeric, Control, Info);
+         umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->HostReadData(),
+                          x.HostWrite(), b.HostRead(), Numeric, Control,
+                          Info);
       umfpack_dl_report_info(Control, Info);
       if (status < 0)
       {
@@ -3230,7 +3232,8 @@ void UMFPackSolver::MultTranspose(const Vector &b, Vector &x) const
    {
       int status =
          umfpack_di_solve(UMFPACK_A, mat->HostReadI(), mat->HostReadJ(),
-                          mat->HostReadData(), x, b, Numeric, Control, Info);
+                          mat->HostReadData(), x.HostWrite(), b.HostRead(),
+                          Numeric, Control, Info);
       umfpack_di_report_info(Control, Info);
       if (status < 0)
       {
@@ -3242,8 +3245,9 @@ void UMFPackSolver::MultTranspose(const Vector &b, Vector &x) const
    else
    {
       SuiteSparse_long status =
-         umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->HostReadData(), x, b,
-                          Numeric, Control, Info);
+         umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->HostReadData(),
+                          x.HostWrite(), b.HostRead(), Numeric, Control,
+                          Info);
       umfpack_dl_report_info(Control, Info);
       if (status < 0)
       {

linalg/solvers.hpp

@@ -899,6 +899,346 @@ public:
    virtual void Mult(const Vector &xt, Vector &x) const;
 };
 
+/** Defines operators and constraints for the following nonlinear contact
+ *  problem:
+ *
+ *    Find x in R^n and lambda in R^m such that
+ *    K(x) - \nabla g(x)^T lambda = f
+ *    0 <= g(x) \perp lambda >= 0
+ *    
+ *    Above, the perpendicularity ("\perp") equation is understood pointwise,
+ *    i.e., g_i(x) lambda_i = 0, for all i=1,2,...,m.
+ *
+ *  This is class is to be derived by concrete contact problems and will be
+ * called by the contact solver to evaluate the contact problem.
+ */
+class ContactProblem
+{
+protected:
+   // Problem sizes
+   const int m, n;
+   // Right-hand side of the contact problem
+   mutable Vector f;
+  
+public:
+   /// Constructor taking the problem sizes as parameters.
+   ContactProblem(int m, int n);
+
+   inline int GetNumConstraints() const { return m; }
+   inline int GetNumDOFs () { return n; }
+
+   // [TODO:] an alternative design would be via a pure virtual method `Eval()`
+   // a members+getters for the problem's rhs, K, gradK, g, and gradG. The
+   // downside is that for nonlinear IPMs, sometimes g and K are evaluated without
+   // needing to evaluate their derivatives, gradG and gradK.
+  
+   /// Returns the rhs of the (first set of equations of the ) contact problem.
+   virtual void GetRhs(Vector &f_out) { f_out = f; }
+  
+   /** This method evaluates the nonlinear function K from the first set of contact
+    * equations. The contact solver calls this function. Concrete contact 
+    * problems need to implement this callback. The method should return false
+    * when the evaluation failed, otherwise should return true. */
+   virtual bool EvalK(const Vector &x, Vector &Katx) = 0;
+
+   /** This method evaluates the gradient of the function K at x. The contact 
+    * solver calls this function. Concrete contact problems need to implement this
+    * callback. The method should return false when the evaluation  failed, otherwise 
+    * should return true. */
+   virtual bool EvalGradK(const Vector &x, Operator &gradKatx) = 0;
+  
+   /** This method evaluates the contact gap function g at x. 
+    * The contact solver calls this function and concrete contact problems will
+    * implement this callback. The method should return false when the 
+    * evaluation failed, otherwise should return true. */
+   virtual bool Evalg(const Vector &x, Vector &g) = 0;
+
+   /** This method evaluates the gradient (nonlinear operator) of the gap 
+    * function g at x. The contact solver calls this function. Concrete contact 
+    * problems need to implement this callback. The method should return false
+    * when the evaluation failed, otherwise should return true. */
+   virtual bool EvalGradg(const Vector &x, Operator &gradG) = 0;
+
+   //TODO: Hessian of the Lagrangian
+};
+
+//forward declaration of the "contact" linear solver
+class ContactLinearSolver;
+
+/// Abstract solver for ContactProblems.
+class ContactSolver : public IterativeSolver
+{
+protected:
+   ContactProblem *problem;
+   ContactLinearSolver* linear_solver;
+public:
+   ContactSolver()
+     : IterativeSolver(), problem(nullptr), linear_solver(nullptr) { }
+#ifdef MFEM_USE_MPI
+   ContactSolver(MPI_Comm comm_)
+     : IterativeSolver(comm_), problem(nullptr), linear_solver(nullptr) { }
+#endif
+   virtual ~ContactSolver() { }
+  
+   /** This method is virtual as solvers might need to perform some initial
+    *  actions (e.g., validation) with the ContactProblem. */
+   virtual void SetContactProblem(ContactProblem &prob)
+   { problem = &prob; }
+
+   virtual void SetLinearSolver(ContactLinearSolver &lsolver)
+   { linear_solver = &lsolver; }
+  
+   /** This method performs the numerical solve of the complementarity 
+    * problem. Potential implementations of this purely virtual method 
+    * will be available via derived classes, e.g. implementations of 
+    * Uzawa (Augmented Lagrangian), interior-point method, etc. 
+    * 
+    * [TODO: work around the two input parameters to pass the initial 
+    * point and return the solution. Are there multivector vectors in 
+    * MFEM?] */
+   ///[TODO: clarify the const `Mult`;
+   virtual void Mult(const Vector &xt, Vector &x) const = 0;
+
+   virtual void SetPreconditioner(Solver &pr)
+   { MFEM_ABORT("Not meaningful for this solver."); }
+   virtual void SetOperator(const Operator &op)
+   { MFEM_ABORT("Not meaningful for this solver."); }
+};
+
+#if 0
+// illustration of the contact nonlinear (outer) loop and interplay with
+// linear solver and contact problem
+void ContactSolverIPM::Mult(const Vector &xin, Vector &xout) const
+{
+  x_iter.SetToZero();
+  s_iter.SetToZero();
+  l_iter.SetToZero();
+  R.SetToZero();
+  mu = 1.0;
+  int num_iter = 0;
+  
+  //solve loop
+  while(true) {
+    //evaluate contact problem at x_iter
+    problem->GetRhs(f);
+    problem->EvalK(x_iter, K);
+    problem->EvalGradK(x_iter, GradK);
+    problem->Evalg(x_iter, g);
+    problem->EvalGradG(x_iter, gradG);
+
+    //evaluate residuals of the contact problem (illustration for IPM)
+    // rx = f + gradG^T l_iter - K
+    // rs = s - g
+    // rc = mu*ones(m,1) - l_iter .* s_iter
+    this->ComputeResiduals(rx, rs, rc);
+
+    //check residual norms, max # of iterations, etc.
+    if(this->StoppingCriteriaMet(num_iter, rx, rs, rc)) {
+      break;
+    }
+    num_iter++;
+    
+    
+    //
+    //set up linear system and solve to obtain search
+    //direction/step  dx, ds, dl
+    //
+
+    //right hand side for linear system for the 
+    // predictor step rc_pred = - l_iter .* s_iter
+    this->ComputePredictorResidual(rc_pred)
+    
+    linear_solver->SetHessian(K);
+    linear_solver->SetGradG(gradG)
+    linear_solver->SetRhs(rx, rs, rc_pred);
+    linear_solver->SetLambda(l_iter);
+    linear_solver->SetS(s_iter);
+    linear_solver->SetPrimalReg(R);
+
+    linear_solver->Mult(dx, ds, dl);
+    //probably more like
+    //linear_solver->Mult(rx, rs, rc_pred, dx, ds, dl);
+    
+    //was the linear solve successful?
+    while(!lin_sol_ok) {
+      this->ComputeRegularization(R);
+      linear_solver->SetPrimalReg(R);
+      linear_solver->Mult(dx, ds, dl);
+      //check linear solve ok
+    }
+
+    //uses x_iter, s_iter, l_iter, dx, ds, dl
+    double alpha = this->LineSearchBlocking();
+    //uses x_iter, s_iter, l_iter, dx, ds, dl, alpha;
+    double sigma = this->ComputeCenteringParameter(alpha); 
+    
+    //right hand side for linear system for the 
+    // corrector step rc_pred = mu*sigma*ones(m,1) + ds.*dl - l_iter .* s_iter
+    //also uses s_iter, l_iter, ds, dl,
+    this->ComputeCorrectorResidual(mu, sigma, rc_corr);
+
+    linear_solver->SetRhs(rx, rs, rc_corr);
+    linear_solver->Mult(dx, ds, dl);
+    
+    //was the linear solve successful?
+    while(!lin_sol_ok) {
+      this->ComputeRegularization(R);
+      linear_solver->SetPrimalReg(R);
+      linear_solver->Mult(dx, ds, dl);
+      //check linear solve ok
+    }
+
+    //uses x_iter, s_iter, l_iter, dx, ds, dl
+    double alpha_corr = this->LineSearchCorrector();
+
+    //x = x+alpha*dx  s=s+alpha*ds  l=l+alpha*dl
+    this->UpdateIter(alpha_corr);
+  } //end of contact loop
+
+  //return x_iter, s_iter, l_iter as the solution
+  //[TODO:]
+}
+
+#endif
+
+#if 0
+//illustration of the "main" driver
+int main()
+{
+  //dream big
+  //
+  mfem::ContactProblem* problem = new mfem::apps::ContactProblemMortar(/*mesh_files*/);
+  //this is curently doable
+  //mfem::ContactProblem* problem = new mfem::apps::ContactLinearProblemFromFiles(/*matrix_files*/);
+
+  //linear solver based on options
+  mfem::ContactLinearSolver* lin_solver = nullptr;
+  if(that) {
+    lin_solver = new mfem::ContactLinearSolverCondensed();
+  } else {
+    //testing
+    lin_solver = new mfem::ContactLinearSolverCondensendDirect();
+  }
+  mfem::ContactSolver* solver = new mfem::ContactSolverIPM();
+
+  solver->SetContactProblem(*problem);
+  solver->SetLinearSolver(*lin_solver);
+
+  solver->Mult();
+  //error control = ?
+
+  //[TODO: what does the interface(s) need to accomdate solving a sequence of (related) problems
+  
+  //clean up
+  delete lin_solver;
+  delete problem;
+}
+#endif
+
+/** Abstract linear solver for internal use in ContactSolver. Implementations
+ * of this class will solve linearizations of the contact problem in the form
+ * shown below. The unknowns are vectors dx, ds, and dlambda
+ * 
+ *    (H+R)dx - gradG^T dlambda = r1
+ *    gradG dx - ds = r2
+ *    Lambda ds + S dlambda = r3
+ * 
+ * Above, H is an nxn symm. psd matrix, for example can be
+ *    H = gradK(x)-\sum \lambda_i \nabla^2 g_i(x), 
+ * or
+ *    H = gradK(x),
+ * depending on the contact solver. gradG is an mxn matrix and Lambda and
+ * S are diagonal matrices with nonnegative entries. The regularization R 
+ * is also a diagonal matrix with nonnegative entries.
+ * 
+ */
+class ContactLinearSolver : public Solver
+{
+public:
+  // [TODO: better names?]
+  void SetHessian(const Operator &gradK);
+  void SetRhs(const Vector &r1, const Vector &r2, const Vector &r3);
+  void SetGradG(const Operator& gradG);
+  void SetLambda(const Vector& Lambda);
+  void SetS(const Vector& S);
+  void SetPrimalReg(const Vector& R);
+  //to be continued...
+   
+   virtual void Mult(const Vector &x_in, Vector &x_out) const = 0;
+protected:
+};
+
+/** Implementation of the contact linear solve in the form of
+ * 
+ *    (H + R + gradG^T S^{-1} Lambda gradG) dx = 
+ *                r1 + gradG^T S^{-1} (Lambda r2 + r3)
+ *    ds = gradG dx - r2
+ *    dlambda = -S^{-1} Lambda ds + S^{-1}r3
+ * 
+ * Essentially, the solve is a Gauss elimination of the linearization
+ * described in ContactLinearSolver. The matrix from the first set
+ * of equations is generally symm. positive definite. The solve
+ * done by this class can/is allowed to fail, in which case the
+ * subsequent linear solve will be done with a "more" p.d. R.
+ * [TODO: ask about error control]
+ *
+ * This applies to both interior-point and Uzawa linearizations. To
+ * use this class with Uzawa, one needs to set S to the identity 
+ * matrix and Lambda to a diagonal matrix having entry (i,i) equal to
+ * one when the corresponding gap is active (i.e., g_i(x)=0) or to 
+ * zero  otherwise, for all i=1,2,...,m. [TODO: what about r3?]
+ * 
+ * The solver uses AMG [TODO: update me].
+ */
+class ContactLinearSolverCondensed : public ContactLinearSolver
+{
+public:
+   virtual void Mult(const Vector &dx_in, Vector &dx_out) const
+   {
+     MFEM_ASSERT(false, "to be implemented");
+   }
+protected:
+};
+
+/** Implementation of the contact linear solve in the form of the
+ * general symmetric indefinite linear system shown above in the
+ * ContactLinearSolver class using direct linear solvers. This is
+ * mostly for testing and debugging purposes.
+ * 
+ * [TODO: update me] Using a symmetric direct linear solver such
+ * as MA57, Pardiso, WSMP, etc. would do it. LU solvers OK.
+ */
+class ContactLinearSolverDirect : public ContactLinearSolver
+{
+public:
+   virtual void Mult(const Vector &x_in, Vector &x_out) const
+   {
+     MFEM_ASSERT(false, "to be implemented");
+   }
+protected:
+};
+
+/** Implementation of the contact linear solve in the form of the
+ * symmetric p.d. linear system shown above in the
+ * ContactLinearSolverCondensed class using direct linear solvers. 
+ * This is mostly for testing and debugging purposes.
+ * 
+ * [TODO: update me] Using a Cholesky direct linear solver such
+ * as Cholmod would do it. LU or symmetric indefinite direct 
+ * solvers OK.
+ * 
+ */
+class ContactLinearSolverCondensendDirect : public ContactLinearSolverCondensed
+{
+public:
+   virtual void Mult(const Vector &x_in, Vector &x_out) const
+   {
+     MFEM_ASSERT(false, "to be implemented");
+   }
+protected:
+};
+
+
 /** Block ILU solver:
  *  Performs a block ILU(k) approximate factorization with specified block
  *  size. Currently only k=0 is supported. This is useful as a preconditioner

linalg/sparsemat.cpp

@@ -15,6 +15,7 @@
 #include "../general/forall.hpp"
 #include "../general/table.hpp"
 #include "../general/sort_pairs.hpp"
+#include "../general/backends.hpp"
 
 #include <iostream>
 #include <iomanip>
@@ -462,25 +463,108 @@ void SparseMatrix::SortColumnIndices()
       return;
    }
 
-   const int * Ip=HostReadI();
-   HostReadWriteJ();
-   HostReadWriteData();
-
-   Array<Pair<int,double> > row;
-   for (int j = 0, i = 0; i < height; i++)
+#ifdef MFEM_USE_CUDA_OR_HIP
+   if ( Device::Allows( Backend::CUDA_MASK ))
    {
-      int end = Ip[i+1];
-      row.SetSize(end - j);
-      for (int k = 0; k < row.Size(); k++)
+#if defined(MFEM_USE_CUDA)
+      size_t pBufferSizeInBytes = 0;
+      void *pBuffer = NULL;
+
+      const int n = Height();
+      const int m = Width();
+      const int nnzA = J.Capacity();
+      double * d_a_sorted = ReadWriteData();
+      const int * d_ia = ReadI();
+      int * d_ja_sorted = ReadWriteJ();
+      csru2csrInfo_t sortInfoA;
+
+      cusparseMatDescr_t matA_descr;
+      cusparseCreateMatDescr( &matA_descr );
+      cusparseSetMatIndexBase( matA_descr, CUSPARSE_INDEX_BASE_ZERO );
+      cusparseSetMatType( matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL );
+
+      cusparseCreateCsru2csrInfo( &sortInfoA );
+
+      cusparseDcsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
+                                       d_ja_sorted, sortInfoA,
+                                       &pBufferSizeInBytes);
+
+      CuMemAlloc( &pBuffer, pBufferSizeInBytes );
+
+      cusparseDcsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
+                         d_ja_sorted, sortInfoA, pBuffer);
+
+      cusparseDestroyCsru2csrInfo( sortInfoA );
+      cusparseDestroyMatDescr( matA_descr );
+
+      CuMemFree( pBuffer );
+#endif
+   }
+   else if ( Device::Allows( Backend::HIP_MASK ))
+   {
+#if defined(MFEM_USE_HIP)
+      size_t pBufferSizeInBytes = 0;
+      void *pBuffer = NULL;
+      int *P = NULL;
+
+      const int n = Height();
+      const int m = Width();
+      const int nnzA = J.Capacity();
+      double * d_a_sorted = ReadWriteData();
+      const int * d_ia = ReadI();
+      int * d_ja_sorted = ReadWriteJ();
+
+      hipsparseMatDescr_t descrA;
+      hipsparseCreateMatDescr( &descrA );
+      // FIXME: There is not in-place version of csr sort in hipSPARSE currently, so we make
+      //        a temporary copy of the data for gthr, sort that, and then copy the sorted values
+      //        back to the array being returned. Where there is an in-place version available,
+      //        we should use it.
+      Array< double > a_tmp( nnzA );
+      double *d_a_tmp = a_tmp.Write();
+
+      hipsparseXcsrsort_bufferSizeExt(handle, n, m, nnzA, d_ia, d_ja_sorted,
+                                      &pBufferSizeInBytes);
+
+      HipMemAlloc( &pBuffer, pBufferSizeInBytes );
+      HipMemAlloc( (void**)&P, nnzA * sizeof(int) );
+
+      hipsparseCreateIdentityPermutation(handle, nnzA, P);
+      hipsparseXcsrsort(handle, n, m, nnzA, descrA, d_ia, d_ja_sorted, P, pBuffer);
+
+      hipsparseDgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
+                     HIPSPARSE_INDEX_BASE_ZERO);
+
+      A.CopyFrom( a_tmp.GetMemory(), nnzA );
+      hipsparseDestroyMatDescr( descrA );
+
+      HipMemFree( pBuffer );
+      HipMemFree( P );
+#endif
+   }
+   else
+#endif // MFEM_USE_CUDA_OR_HIP
+   {
+      const int * Ip=HostReadI();
+      HostReadWriteJ();
+      HostReadWriteData();
+
+      Array<Pair<int,double> > row;
+      for (int j = 0, i = 0; i < height; i++)
       {
-         row[k].one = J[j+k];
-         row[k].two = A[j+k];
-      }
-      row.Sort();
-      for (int k = 0; k < row.Size(); k++, j++)
-      {
-         J[j] = row[k].one;
-         A[j] = row[k].two;
+         int end = Ip[i+1];
+         row.SetSize(end - j);
+         for (int k = 0; k < row.Size(); k++)
+         {
+            row[k].one = J[j+k];
+            row[k].two = A[j+k];
+         }
+         row.Sort();
+         for (int k = 0; k < row.Size(); k++, j++)
+         {
+            J[j] = row[k].one;
+            A[j] = row[k].two;
+         }
       }
    }
    isSorted = true;

linalg/sparsemat.hpp

@@ -23,7 +23,11 @@
 #include "densemat.hpp"
 
 #if defined(MFEM_USE_HIP)
+#if (HIP_VERSION_MAJOR * 100 + HIP_VERSION_MINOR) < 502
 #include <hipsparse.h>
+#else
+#include <hipsparse/hipsparse.h>
+#endif
 #endif
 
 

linalg/strumpack.cpp

@@ -56,7 +56,9 @@ STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(const HypreParMatrix & hypParMat)
 
    // Create the CSRMatrixMPI A_ by borrowing the internal data from a
    // hypre_CSRMatrix.
+   hypParMat.HostRead();
    hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+   hypParMat.HypreRead();
    hypre_CSRMatrixSetDataOwner(csr_op,0);
 #if MFEM_HYPRE_VERSION >= 21600
    // For now, this method assumes that HYPRE_Int is int. Also, csr_op->num_cols
@@ -201,8 +203,8 @@ void STRUMPACKSolver::Mult( const Vector & x, Vector & y ) const
    MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
                << ", expected size = " << Height());
 
-   double*  yPtr = (double*)y;
-   double*  xPtr = (double*)(const_cast<Vector&>(x));
+   double*  yPtr = y.HostWrite();
+   const double*  xPtr = x.HostRead();
 
    solver_->options().set_verbose( factor_verbose_ );
    ReturnCode ret = solver_->factor();
@@ -219,6 +221,10 @@ void STRUMPACKSolver::Mult( const Vector & x, Vector & y ) const
          MFEM_ABORT("STRUMPACK:  Matrix reordering failed!");
       }
       break;
+      default:
+      {
+         MFEM_ABORT("STRUMPACK: 'factor()' error code = " << ret);
+      }
    }
    solver_->options().set_verbose( solve_verbose_ );
    solver_->solve(xPtr, yPtr);

linalg/superlu.cpp

@@ -149,7 +149,9 @@ SuperLURowLocMatrix::SuperLURowLocMatrix( const HypreParMatrix & hypParMat )
 
    // Create the SuperMatrix A by borrowing the internal data from a
    // hypre_CSRMatrix.
+   hypParMat.HostRead();
    hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+   hypParMat.HypreRead();
    hypre_CSRMatrixSetDataOwner(csr_op,0);
 #if MFEM_HYPRE_VERSION >= 21600
    // For now, this method assumes that HYPRE_BigInt is int. Also, csr_op->num_cols

linalg/vector.cpp

@@ -786,16 +786,13 @@ void Vector::Randomize(int seed)
 {
    const double max = (double)(RAND_MAX) + 1.;
 
-   if (!global_seed_set)
+   if (seed == 0)
    {
-      if (seed == 0)
-      {
-         seed = (int)time(0);
-      }
-
-      srand((unsigned)seed);
+      seed = (int)time(0);
    }
 
+   srand((unsigned)seed);
+
    HostWrite();
    for (int i = 0; i < size; i++)
    {
@@ -803,12 +800,6 @@ void Vector::Randomize(int seed)
    }
 }
 
-void Vector::SetGlobalSeed(int gseed)
-{
-   srand((unsigned)gseed);
-   global_seed_set = true;
-}
-
 double Vector::Norml2() const
 {
    // Scale entries of Vector on the fly, using algorithms from

linalg/vector.hpp

@@ -63,7 +63,6 @@ protected:
 
    Memory<double> data;
    int size;
-   bool global_seed_set = false;
 
 public:
 
@@ -209,15 +208,11 @@ public:
    inline double *GetData() const
    { return const_cast<double*>((const double*)data); }
 
-   /// Conversion to `double *`.
-   /** @note This conversion function makes it possible to use [] for indexing
-       in addition to the overloaded operator()(int). */
-   inline operator double *() { return data; }
+   /// Conversion to `double *`. Deprecated.
+   MFEM_DEPRECATED inline operator double *() { return data; }
 
-   /// Conversion to `const double *`.
-   /** @note This conversion function makes it possible to use [] for indexing
-       in addition to the overloaded operator()(int). */
-   inline operator const double *() const { return data; }
+   /// Conversion to `const double *`. Deprecated.
+   MFEM_DEPRECATED inline operator const double *() const { return data; }
 
    /// STL-like begin.
    inline double *begin() { return data; }
@@ -269,6 +264,14 @@ public:
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
    inline const double &operator()(int i) const;
 
+   /// Access Vector entries using [] for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline double &operator[](int i) { return (*this)(i); }
+
+   /// Read only access to Vector entries using [] for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline const double &operator[](int i) const { return (*this)(i); }
+
    /// Dot product with a `double *` array.
    double operator*(const double *) const;
 
@@ -307,12 +310,6 @@ public:
 
    Vector &operator+=(const Vector &v);
 
-   /// operator- is not supported. Use @ref subtract or @ref Add.
-   Vector &operator-(const Vector &v) = delete;
-
-   /// operator+ is not supported. Use @ref Add.
-   Vector &operator+(const Vector &v) = delete;
-
    /// (*this) += a * Va
    Vector &Add(const double a, const Vector &Va);
 
@@ -417,8 +414,6 @@ public:
 
    /// Set random values in the vector.
    void Randomize(int seed = 0);
-   /// Set global seed for random values in sequential calls to Randomize().
-   void SetGlobalSeed(int gseed);
    /// Returns the l2 norm of the vector.
    double Norml2() const;
    /// Returns the l_infinity norm of the vector.
@@ -435,8 +430,12 @@ public:
    double Sum() const;
    /// Compute the square of the Euclidean distance to another vector.
    inline double DistanceSquaredTo(const double *p) const;
+   /// Compute the square of the Euclidean distance to another vector.
+   inline double DistanceSquaredTo(const Vector &p) const;
    /// Compute the Euclidean distance to another vector.
    inline double DistanceTo(const double *p) const;
+   /// Compute the Euclidean distance to another vector.
+   inline double DistanceTo(const Vector &p) const;
 
    /** @brief Count the number of entries in the Vector for which isfinite
        is false, i.e. the entry is a NaN or +/-Inf. */
@@ -646,16 +645,33 @@ inline double Distance(const double *x, const double *y, const int n)
    return std::sqrt(DistanceSquared(x, y, n));
 }
 
+inline double Distance(const Vector &x, const Vector &y)
+{
+   return x.DistanceTo(y);
+}
+
 inline double Vector::DistanceSquaredTo(const double *p) const
 {
    return DistanceSquared(data, p, size);
 }
 
+inline double Vector::DistanceSquaredTo(const Vector &p) const
+{
+   MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
+   return DistanceSquared(data, p.data, size);
+}
+
 inline double Vector::DistanceTo(const double *p) const
 {
    return Distance(data, p, size);
 }
 
+inline double Vector::DistanceTo(const Vector &p) const
+{
+   MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
+   return Distance(data, p.data, size);
+}
+
 /// Returns the inner product of x and y
 /** In parallel this computes the inner product of the local vectors,
     producing different results on each MPI rank.

makefile

@@ -10,7 +10,7 @@
 # CONTRIBUTING.md for details.
 
 # The current MFEM version as an integer, see also `CMakeLists.txt`.
-MFEM_VERSION = 40500
+MFEM_VERSION = 40501
 MFEM_VERSION_STRING = $(shell printf "%06d" $(MFEM_VERSION) | \
   sed -e 's/^0*\(.*.\)\(..\)\(..\)$$/\1.\2.\3/' -e 's/\.0/./g' -e 's/\.0$$//')
 
@@ -275,7 +275,7 @@ endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = ENZYME SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
- MESQUITE SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP\
+ SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER PARELAG BENCHMARK\
  MOONOLITH ALGOIM
 
@@ -338,7 +338,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_METIS MFEM_USE_METIS_5 MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB\
  MFEM_USE_LIBUNWIND MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_OPENMP\
  MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS\
- MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU\
+ MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
@@ -673,9 +673,9 @@ status info:
 	$(info MFEM_USE_MEMALLOC      = $(MFEM_USE_MEMALLOC))
 	$(info MFEM_TIMER_TYPE        = $(MFEM_TIMER_TYPE))
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
-	$(info MFEM_USE_MESQUITE      = $(MFEM_USE_MESQUITE))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
+	$(info MFEM_USE_SUPERLU5      = $(MFEM_USE_SUPERLU5))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))
 	$(info MFEM_USE_GINKGO        = $(MFEM_USE_GINKGO))

mesh/CMakeLists.txt

@@ -72,13 +72,6 @@ if (MFEM_USE_MPI)
     submesh/psubmesh.hpp)
 endif()
 
-if (MFEM_USE_MESQUITE)
-  list(APPEND SRCS mesquite.cpp)
-  # If this list (HDRS -> HEADERS) is used for install, we probably want the
-  # header added all the time.
-  list(APPEND HDRS mesquite.hpp)
-endif()
-
 if (MFEM_USE_PUMI)
   list(APPEND SRCS
     pumi.cpp)

mesh/mesh.cpp

@@ -1627,6 +1627,13 @@ int Mesh::AddVertex(const double *coords)
    return NumOfVertices++;
 }
 
+int Mesh::AddVertex(const Vector &coords)
+{
+   MFEM_ASSERT(coords.Size() >= spaceDim,
+               "invalid 'coords' size: " << coords.Size());
+   return AddVertex(coords.GetData());
+}
+
 void Mesh::AddVertexParents(int i, int p1, int p2)
 {
    tmp_vertex_parents.Append(Triple<int, int, int>(i, p1, p2));
@@ -7025,6 +7032,11 @@ int *Mesh::CartesianPartitioning(int nxyz[])
    return partitioning;
 }
 
+void FindPartitioningComponents(Table &elem_elem,
+                                const Array<int> &partitioning,
+                                Array<int> &component,
+                                Array<int> &num_comp);
+
 int *Mesh::GeneratePartitioning(int nparts, int part_method)
 {
 #ifdef MFEM_USE_METIS
@@ -7106,6 +7118,15 @@ int *Mesh::GeneratePartitioning(int nparts, int part_method)
 #else
       METIS_SetDefaultOptions(options);
       options[METIS_OPTION_CONTIG] = 1; // set METIS_OPTION_CONTIG
+      // If the mesh is disconnected, disable METIS_OPTION_CONTIG.
+      {
+         Array<int> part(partitioning, NumOfElements);
+         part = 0; // single part for the whole mesh
+         Array<int> component; // size will be set to num. elem.
+         Array<int> num_comp;  // size will be set to num. parts (1)
+         FindPartitioningComponents(*el_to_el, part, component, num_comp);
+         if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
+      }
 #endif
 
       // Sort the neighbor lists
@@ -10504,8 +10525,8 @@ void Mesh::PrintVTU(std::ostream &os, int ref, VTKFormat format,
    };
 
    int ne = bdr_elements ? GetNBE() : GetNE();
-   // count the points, cells, size
-   int np = 0, nc_ref = 0, size = 0;
+   // count the number of points and cells
+   int np = 0, nc_ref = 0;
    for (int i = 0; i < ne; i++)
    {
       Geometry::Type geom = get_geom(i);
@@ -10513,7 +10534,6 @@ void Mesh::PrintVTU(std::ostream &os, int ref, VTKFormat format,
       RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
       np += RefG->RefPts.GetNPoints();
       nc_ref += RefG->RefGeoms.Size() / nv;
-      size += (RefG->RefGeoms.Size() / nv) * (nv + 1);
    }
 
    os << "<Piece NumberOfPoints=\"" << np << "\" NumberOfCells=\""
@@ -11636,6 +11656,7 @@ void Mesh::Transform(void (*f)(const Vector&, Vector&))
       xnew.ProjectCoefficient(f_pert);
       *Nodes = xnew;
    }
+   NodesUpdated();
 }
 
 void Mesh::Transform(VectorCoefficient &deformation)
@@ -11660,6 +11681,7 @@ void Mesh::Transform(VectorCoefficient &deformation)
       xnew.ProjectCoefficient(deformation);
       *Nodes = xnew;
    }
+   NodesUpdated();
 }
 
 void Mesh::RemoveUnusedVertices()

mesh/mesh.hpp

@@ -272,8 +272,8 @@ public:
    NURBSExtension *NURBSext; ///< Optional NURBS mesh extension.
    NCMesh *ncmesh;           ///< Optional nonconforming mesh extension.
    Array<GeometricFactors*> geom_factors; ///< Optional geometric factors.
-   Array<FaceGeometricFactors*>
-   face_geom_factors; ///< Optional face geometric factors.
+   Array<FaceGeometricFactors*> face_geom_factors; /**< Optional face geometric
+                                                        factors. */
 
    // Global parameter that can be used to control the removal of unused
    // vertices performed when reading a mesh in MFEM format. The default value
@@ -380,7 +380,14 @@ protected:
    void InitRefinementTransforms();
    int FindCoarseElement(int i);
 
-   /// Update the nodes of a curved mesh after refinement
+   /** @brief Update the nodes of a curved mesh after the topological part of a
+       Mesh::Operation, such as refinement, has been performed. */
+   /** If Nodes GridFunction is defined, i.e. not NULL, this method calls
+       NodesUpdated().
+
+       @note Unlike the similarly named public method NodesUpdated() this
+       method modifies the mesh nodes (if they exist) and calls NodesUpdated().
+   */
    void UpdateNodes();
 
    /// Helper to set vertex coordinates given a high-order curvature function.
@@ -705,6 +712,7 @@ public:
 
    int AddVertex(double x, double y = 0.0, double z = 0.0);
    int AddVertex(const double *coords);
+   int AddVertex(const Vector &coords);
    /// Mark vertex @a i as nonconforming, with parent vertices @a p1 and @a p2.
    void AddVertexParents(int i, int p1, int p2);
 
@@ -964,10 +972,19 @@ public:
 
        The IntegrationRule used with GetGeometricFactors needs to remain valid
        until the internally stored GeometricFactors objects are destroyed (by
-       either calling Mesh::DeleteGeometricFactors or the Mesh destructor). If
-       the device MemoryType parameter @a d_mt is specified, then the returned
-       object will use that type unless it was previously allocated with a
-       different type. */
+       calling Mesh::DeleteGeometricFactors(), Mesh::NodesUpdated(), or the Mesh
+       destructor).
+
+       If the device MemoryType parameter @a d_mt is specified, then the
+       returned object will use that type unless it was previously allocated
+       with a different type.
+
+       The returned pointer points to an internal object that may be invalidated
+       by mesh operations such as refinement, vertex/node movement, etc. Since
+       not all such modifications can be tracked by the Mesh class (e.g. when
+       using the pointer returned by GetNodes() to change the nodes) one needs
+       to account for such changes by calling the method NodesUpdated() which,
+       in particular, will call DeleteGeometricFactors(). */
    const GeometricFactors* GetGeometricFactors(
       const IntegrationRule& ir,
       const int flags,
@@ -978,22 +995,41 @@ public:
 
        The IntegrationRule used with GetFaceGeometricFactors needs to remain
        valid until the internally stored FaceGeometricFactors objects are
-       destroyed (by either calling Mesh::DeleteGeometricFactors or the Mesh
-       destructor). */
-   const FaceGeometricFactors* GetFaceGeometricFactors(const IntegrationRule& ir,
-                                                       const int flags,
-                                                       FaceType type,
-                                                       MemoryType d_mt = MemoryType::DEFAULT);
+       destroyed (by either calling Mesh::DeleteGeometricFactors(),
+       Mesh::NodesUpdated(), or the Mesh destructor).
+
+       If the device MemoryType parameter @a d_mt is specified, then the
+       returned object will use that type unless it was previously allocated
+       with a different type.
+
+       The returned pointer points to an internal object that may be invalidated
+       by mesh operations such as refinement, vertex/node movement, etc. Since
+       not all such modifications can be tracked by the Mesh class (e.g. when
+       using the pointer returned by GetNodes() to change the nodes) one needs
+       to account for such changes by calling the method NodesUpdated() which,
+       in particular, will call DeleteGeometricFactors(). */
+   const FaceGeometricFactors* GetFaceGeometricFactors(
+      const IntegrationRule& ir,
+      const int flags,
+      FaceType type,
+      MemoryType d_mt = MemoryType::DEFAULT);
 
    /// Destroy all GeometricFactors stored by the Mesh.
    /** This method can be used to force recomputation of the GeometricFactors,
-       for example, after the mesh nodes are modified externally. */
+       for example, after the mesh nodes are modified externally.
+
+       @note In general, the preferred method for resetting the GeometricFactors
+       should be to call NodesUpdated(). */
    void DeleteGeometricFactors();
 
-   /// @brief This function should be called after the mesh node coordinates
-   /// have changed, e.g. after the mesh has moved.
-   /** It updates internal quantities derived from the node coordinates, such
-       as the GeometricFactors. */
+   /** @brief This function should be called after the mesh node coordinates
+       have been updated externally, e.g. by modifying the internal nodal
+       GridFunction returned by GetNodes(). */
+   /** It deletes internal quantities derived from the node coordinates,
+       such as the (Face)GeometricFactors.
+
+       @note Unlike the similarly named protected method UpdateNodes() this
+       method does not modify the nodes. */
    void NodesUpdated() { DeleteGeometricFactors(); }
 
    /// Equals 1 + num_holes - num_loops
@@ -1535,6 +1571,7 @@ public:
    // mesh is not curved (i.e. Nodes == NULL).
    void MoveNodes(const Vector &displacements);
    void GetNodes(Vector &node_coord) const;
+   /// Updates the vertex/node locations. Invokes NodesUpdated().
    void SetNodes(const Vector &node_coord);
 
    /// Return a pointer to the internal node GridFunction (may be NULL).
@@ -1545,9 +1582,11 @@ public:
    /// Set the mesh nodes ownership flag.
    void SetNodesOwner(bool nodes_owner) { own_nodes = nodes_owner; }
    /// Replace the internal node GridFunction with the given GridFunction.
+   /** Invokes NodesUpdated(). */
    void NewNodes(GridFunction &nodes, bool make_owner = false);
-   /** Swap the internal node GridFunction pointer and ownership flag members
-       with the given ones. */
+   /** @brief Swap the internal node GridFunction pointer and ownership flag
+       members with the given ones. */
+   /** Invokes NodesUpdated(). */
    void SwapNodes(GridFunction *&nodes, int &own_nodes_);
 
    /// Return the mesh nodes/vertices projected on the given GridFunction.
@@ -1783,8 +1822,6 @@ public:
       PrintCharacteristics(NULL, NULL, os);
    }
 
-   void MesquiteSmooth(const int mesquite_option = 0);
-
    /** @brief Find the ids of the elements that contain the given points, and
        their corresponding reference coordinates.
 

mesh/mesh_headers.hpp

@@ -33,10 +33,6 @@
 #include "wedge.hpp"
 #include "pyramid.hpp"
 
-#ifdef MFEM_USE_MESQUITE
-#include "mesquite.hpp"
-#endif
-
 #ifdef MFEM_USE_MPI
 #include "pncmesh.hpp"
 #include "pmesh.hpp"

mesh/mesh_readers.cpp

@@ -2851,9 +2851,10 @@ void Mesh::ReadCubit(const char *filename, int &curved, int &read_gf)
    int retval;
 
    // dummy string
-   char str_dummy[256];
+   constexpr size_t buf_size = 256;
+   char str_dummy[buf_size];
 
-   char temp_str[256];
+   char temp_str[buf_size];
    int temp_id;
 
    // open the file.
@@ -2897,14 +2898,14 @@ void Mesh::ReadCubit(const char *filename, int &curved, int &read_gf)
    int previous_num_node_per_el = 0;
    for (int i = 0; i < (int) num_el_blk; i++)
    {
-      sprintf(temp_str, "num_el_in_blk%d", i+1);
+      snprintf(temp_str, buf_size, "num_el_in_blk%d", i+1);
       if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
           (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_el_in_blk[i])))
       {
          MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
       }
 
-      sprintf(temp_str, "num_nod_per_el%d", i+1);
+      snprintf(temp_str, buf_size, "num_nod_per_el%d", i+1);
       if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
           (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_node_per_el)))
       {
@@ -3051,7 +3052,7 @@ void Mesh::ReadCubit(const char *filename, int &curved, int &read_gf)
    size_t *num_side_in_ss  = new size_t[num_side_sets];
    for (int i = 0; i < (int) num_side_sets; i++)
    {
-      sprintf(temp_str, "num_side_ss%d", i+1);
+      snprintf(temp_str, buf_size, "num_side_ss%d", i+1);
       if ((retval = nc_inq_dimid(ncid, temp_str, &temp_id)) ||
           (retval = nc_inq_dim(ncid, temp_id, str_dummy, &num_side_in_ss[i])))
       {
@@ -3086,7 +3087,7 @@ void Mesh::ReadCubit(const char *filename, int &curved, int &read_gf)
    for (int i = 0; i < (int) num_el_blk; i++)
    {
       elem_blk[i] = new int[num_el_in_blk[i] * num_node_per_el];
-      sprintf(temp_str, "connect%d", i+1);
+      snprintf(temp_str, buf_size, "connect%d", i+1);
       if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
           (retval = nc_get_var_int(ncid, temp_id, elem_blk[i])))
       {
@@ -3110,14 +3111,14 @@ void Mesh::ReadCubit(const char *filename, int &curved, int &read_gf)
       elem_ss[i] = new int[num_side_in_ss[i]];
       side_ss[i] = new int[num_side_in_ss[i]];
 
-      sprintf(temp_str, "elem_ss%d", i+1);
+      snprintf(temp_str, buf_size, "elem_ss%d", i+1);
       if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
           (retval = nc_get_var_int(ncid, temp_id, elem_ss[i])))
       {
          MFEM_ABORT("Fatal NetCDF error: " << nc_strerror(retval));
       }
 
-      sprintf(temp_str,"side_ss%d",i+1);
+      snprintf(temp_str, buf_size,"side_ss%d",i+1);
       if ((retval = nc_inq_varid(ncid, temp_str, &temp_id)) ||
           (retval = nc_get_var_int(ncid, temp_id, side_ss[i])))
       {

mesh/mesquite.hpp

@@ -1,358 +0,0 @@
-// Copyright (c) 2010-2022, 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_MESQUITE
-#define MFEM_MESQUITE
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_MESQUITE
-
-#include "mesh.hpp"
-
-#include "Mesquite_all_headers.hpp"
-
-namespace mfem
-{
-
-using namespace Mesquite;
-
-class MesquiteMesh : public Mesquite::Mesh
-{
-   // tagging interface definition
-private:
-
-   struct MfemTagDescription
-   {
-
-      std::string name;       //!< Tag name
-      Mesh::TagType type;     //!< Tag data type
-      size_t size;            //!< Size of tag data (sizeof(type)*array_length)
-
-      inline MfemTagDescription( std::string n,
-                                 Mesh::TagType t,
-                                 size_t s)
-         : name(n), type(t), size(s) {}
-
-      inline MfemTagDescription( )
-         : type(Mesh::BYTE), size(0) {}
-
-      inline bool operator==(const MfemTagDescription& o) const
-      { return name == o.name && type == o.type && size == o.size; }
-      inline bool operator!=(const MfemTagDescription& o) const
-      { return name != o.name || type != o.type || size != o.size; }
-   };
-
-   /**\class MeshTags
-    *
-    * Store tags and tag data for MFEM's mesh representation.
-    * Stores for each tag: properties, element data, and vertex data.
-    * The tag element and vertex data sets are maps between some element
-    * or vertex index and a tag value.
-    */
-   class MeshTags
-   {
-   public:
-
-      ~MeshTags() { clear(); }
-
-      /** \class TagData
-       * Store data for a single tag
-       */
-      struct TagData
-      {
-
-         //! tag meta data
-         const MfemTagDescription desc;
-
-         //! per-element data, or NULL if none has been set.
-         void* elementData;
-
-         //! number of entries in elementData
-         size_t elementCount;
-
-         //! per-vertex data, or NULL if none has been set.
-         void* vertexData;
-
-         //! number of entries in vertexData
-         size_t vertexCount;
-
-         //! Default value for tag
-         void* defaultValue;
-
-         /** \brief Construct tag
-          *\param name Tag name
-          *\param type Tag data type
-          *\param length Tag array length (1 for scalar/non-array)
-          *\param default_val Default value for tag
-          */
-         inline TagData( const std::string& name,
-                         Mesh::TagType type, unsigned length,
-                         void* default_val = 0)
-            : desc(name, type, length*size_from_tag_type(type)),
-              elementData(0), elementCount(0),
-              vertexData(0), vertexCount(0),
-              defaultValue(default_val) {}
-
-         /** \brief Construct tag
-          *\param descr Tag description object
-          */
-         inline TagData( const MfemTagDescription& descr )
-            : desc(descr), elementData(0), elementCount(0),
-              vertexData(0), vertexCount(0),
-              defaultValue(0) {}
-
-         ~TagData();
-      };
-
-      /** \brief Get the size of the passed data type */
-      static size_t size_from_tag_type( Mesh::TagType type );
-
-      /** \brief Clear all data */
-      void clear();
-
-      /** \brief Get tag index from name */
-      size_t handle( const std::string& name, MsqError& err ) const;
-
-      /** \brief Get tag properties */
-      const MfemTagDescription& properties( size_t tag_handle, MsqError& err ) const;
-
-      /** \brief Create a new tag
-       *
-       * Create a new tag with the passed properties
-       *\param name Tag name (must be unique)
-       *\param type Tag data type
-       *\param length Number of values in tag (array length, 1 for scalar)
-       *\param defval Optional default value for tag
-       */
-      size_t create( const std::string& name,
-                     Mesh::TagType type,
-                     unsigned length,
-                     const void* defval,
-                     MsqError& err );
-
-      /** \brief Create a new tag
-       *
-       * Create a new tag with the passed properties
-       */
-      size_t create( const MfemTagDescription& desc,
-                     const void* defval,
-                     MsqError& err );
-
-      /**\brief Remove a tag */
-      void destroy( size_t tag_index, MsqError& err );
-
-      /**\brief Set tag data on elements */
-      void set_element_data( size_t tag_handle,
-                             size_t num_indices,
-                             const size_t* elem_indices,
-                             const void* tag_data,
-                             MsqError& err );
-
-      /**\brief Set tag data on vertices */
-      void set_vertex_data( size_t tag_handle,
-                            size_t num_indices,
-                            const size_t* elem_indices,
-                            const void* tag_data,
-                            MsqError& err );
-
-      /**\brief Get tag data on elements */
-      void get_element_data( size_t tag_handle,
-                             size_t num_indices,
-                             const size_t* elem_indices,
-                             void* tag_data,
-                             MsqError& err ) const;
-
-      /**\brief Get tag data on vertices */
-      void get_vertex_data( size_t tag_handle,
-                            size_t num_indices,
-                            const size_t* elem_indices,
-                            void* tag_data,
-                            MsqError& err ) const;
-
-      /**\class TagIterator
-       *
-       * Iterate over list of valid tag handles
-       */
-      class TagIterator
-      {
-      public:
-         TagIterator() : tags(0), index(0) {}
-         TagIterator( MeshTags* d, size_t i ) : tags(d), index(i) {}
-         size_t operator*() const  { return index+1; }
-         TagIterator operator++();
-         TagIterator operator--();
-         TagIterator operator++(int);
-         TagIterator operator--(int);
-         bool operator==(TagIterator other) const { return index == other.index; }
-         bool operator!=(TagIterator other) const { return index != other.index; }
-      private:
-         MeshTags* tags;
-         size_t index;
-      };
-      TagIterator tag_begin();
-      TagIterator tag_end()   { return TagIterator(this,tagList.size()); }
-
-      /**\brief Check if any vertices have tag */
-      bool tag_has_vertex_data( size_t index, MsqError& err ) ;
-      /**\brief Check if any elements have tag */
-      bool tag_has_element_data( size_t index, MsqError& err ) ;
-
-   private:
-
-      friend class MeshTags::TagIterator;
-
-      std::vector<TagData*> tagList;
-   }; // class MeshTags
-
-   // data members
-private:
-   int ndofs;                      // number of nodes (or vertices) in mesh
-   int nelems;                     // number of elements in mesh
-   mfem::Mesh *mesh;               // pointer to mfem mesh object
-   mfem::Element *elem;            // pointer to mfem element object
-   mfem::GridFunction *nodes;      // pointer to mfem grid function object
-   // for nodes
-   mfem::FiniteElementSpace *fes;  // pointer to mfem finite element
-   // space object
-   mfem::Table *dof_elem;          // dof to element table
-   std::vector<char> mByte;        // length = ndofs
-   std::vector<bool> mFixed;       // length = ndofs
-
-   MeshTags* myTags;
-
-public:
-
-   // The constructor
-   MesquiteMesh(mfem::Mesh *mfem_mesh);
-
-   // The mesh dimension
-   int get_geometric_dimension(MsqError &err);
-
-   // The handles are just pointers to the indexes of the nodes/elements
-   void get_all_elements(std::vector<ElementHandle>& elements,
-                         MsqError& err );
-   void get_all_vertices(std::vector<VertexHandle>& vertices,
-                         MsqError& err );
-
-   // Get/set vertex coordinates
-   void vertices_get_coordinates(const VertexHandle vert_array[],
-                                 MsqVertex* coordinates,
-                                 size_t num_vtx,
-                                 MsqError &err );
-   void vertex_set_coordinates(VertexHandle vertex,
-                               const Vector3D &coordinates,
-                               MsqError &err );
-
-   // These are internal markers for Mesquite that we should allocate for them
-   void vertex_set_byte(VertexHandle vertex,
-                        unsigned char byte,
-                        MsqError &err);
-   void vertices_set_byte(const VertexHandle *vert_array,
-                          const unsigned char *byte_array,
-                          size_t array_size,
-                          MsqError &err );
-   void vertices_get_fixed_flag(const VertexHandle vert_array[],
-                                std::vector<bool>& fixed_flag_array,
-                                size_t num_vtx,
-                                MsqError &err );
-   void vertices_set_fixed_flag(const VertexHandle vert_array[],
-                                const std::vector< bool > &fixed_flag_array,
-                                size_t num_vtx,
-                                MsqError &err );
-   void vertex_get_byte( const VertexHandle vertex,
-                         unsigned char *byte,
-                         MsqError &err );
-   void vertices_get_byte( const VertexHandle *vertex,
-                           unsigned char *byte_array,
-                           size_t array_size,
-                           MsqError &err );
-   // The dof_elem table
-   void vertices_get_attached_elements(const VertexHandle* vertex_array,
-                                       size_t num_vertex,
-                                       std::vector<ElementHandle>& elements,
-                                       std::vector<size_t>& offsets,
-                                       MsqError& err );
-   // The elem_dof table
-   void elements_get_attached_vertices(const ElementHandle *elem_handles,
-                                       size_t num_elems,
-                                       std::vector<VertexHandle>& vert_handles,
-                                       std::vector<size_t>& offsets,
-                                       MsqError &err);
-
-   // The topology of the elements: tri/tet/quad/hex...
-   void elements_get_topologies(const ElementHandle *element_handle_array,
-                                EntityTopology *element_topologies,
-                                size_t num_elements,
-                                MsqError &err);
-   // The destructor
-   ~MesquiteMesh();
-
-   // tag the attributes on the elements from the underlying mfem mesh
-   void tag_attributes();
-
-   // Clean these up .....
-   void vertices_get_slaved_flag( const VertexHandle vert_array[],
-                                  std::vector<bool>& slaved_flag_array,
-                                  size_t num_vtx,
-                                  MsqError &err ) {};
-
-   TagHandle tag_create( const std::string& tag_name,
-                         TagType type, unsigned length,
-                         const void* default_value,
-                         MsqError &err);
-
-   void tag_delete( TagHandle handle, MsqError& err );
-
-   TagHandle tag_get( const std::string& name,
-                      MsqError& err );
-
-   void tag_properties( TagHandle handle,
-                        std::string& name_out,
-                        TagType& type_out,
-                        unsigned& length_out,
-                        MsqError& err );
-
-   void tag_set_element_data( TagHandle handle,
-                              size_t num_elems,
-                              const ElementHandle* elem_array,
-                              const void* tag_data,
-                              MsqError& err );
-
-   void tag_set_vertex_data ( TagHandle handle,
-                              size_t num_elems,
-                              const VertexHandle* node_array,
-                              const void* tag_data,
-                              MsqError& err );
-
-   void tag_get_element_data( TagHandle handle,
-                              size_t num_elems,
-                              const ElementHandle* elem_array,
-                              void* tag_data,
-                              MsqError& err );
-
-   void tag_get_vertex_data ( TagHandle handle,
-                              size_t num_elems,
-                              const VertexHandle* node_array,
-                              void* tag_data,
-                              MsqError& err );
-
-   void release_entity_handles(const EntityHandle *handle_array,
-                               size_t num_handles,
-                               MsqError &err) {};
-   void release() {};
-};
-
-}
-
-#endif
-
-#endif

mesh/ncmesh.cpp

@@ -33,10 +33,11 @@ void NCMesh::GeomInfo::InitGeom(Geometry::Type geom)
    {
       case Geometry::CUBE: elem = new Hexahedron; break;
       case Geometry::PRISM: elem = new Wedge; break;
-      case Geometry::SQUARE: elem = new Quadrilateral; break;
-      case Geometry::SEGMENT: elem = new Segment; break;
-      case Geometry::TRIANGLE: elem = new Triangle; break;
       case Geometry::TETRAHEDRON: elem = new Tetrahedron; break;
+      case Geometry::PYRAMID: elem = new Pyramid; break;
+      case Geometry::SQUARE: elem = new Quadrilateral; break;
+      case Geometry::TRIANGLE: elem = new Triangle; break;
+      case Geometry::SEGMENT: elem = new Segment; break;
       default: MFEM_ABORT("unsupported geometry " << geom);
    }
 
@@ -62,7 +63,7 @@ void NCMesh::GeomInfo::InitGeom(Geometry::Type geom)
       }
    }
 
-   // in 1D & 2D we pretend to have faces too, so we can use NCMesh::Face::elem[2]
+   // in 1D/2D we pretend to have faces too, so we can use NCMesh::Face::elem[2]
    if (!nf)
    {
       if (ne)
@@ -97,7 +98,7 @@ static void CheckSupportedGeom(Geometry::Type geom)
    MFEM_VERIFY(geom == Geometry::SEGMENT ||
                geom == Geometry::TRIANGLE || geom == Geometry::SQUARE ||
                geom == Geometry::CUBE || geom == Geometry::PRISM ||
-               geom == Geometry::TETRAHEDRON,
+               geom == Geometry::PYRAMID || geom == Geometry::TETRAHEDRON,
                "Element type " << geom << " is not supported by NCMesh.");
 }
 
@@ -119,6 +120,12 @@ NCMesh::NCMesh(const Mesh *mesh)
       CheckSupportedGeom(geom);
       GI[geom].InitGeom(geom);
 
+      // if we have pyramids we will need tets after refinement
+      if (geom == Geometry::PYRAMID)
+      {
+         GI[Geometry::TETRAHEDRON].InitGeom(Geometry::TETRAHEDRON);
+      }
+
       // create NCMesh::Element for this mfem::Element
       int root_id = AddElement(Element(geom, elem->GetAttribute()));
       MFEM_ASSERT(root_id == i, "");
@@ -454,9 +461,10 @@ NCMesh::Element::Element(Geometry::Type geom, int attr)
    : geom(geom), ref_type(0), tet_type(0), flag(0), index(-1)
    , rank(0), attribute(attr), parent(-1)
 {
-   for (int i = 0; i < 8; i++) { node[i] = -1; }
+   for (int i = 0; i < MaxElemNodes; i++) { node[i] = -1; }
+   for (int i = 0; i < MaxElemChildren; i++) { child[i] = -1; }
 
-   // NOTE: in 2D the 8-element node/child arrays are not optimal, however,
+   // NOTE: in 2D the 8/10-element node/child arrays are not optimal, however,
    // testing shows we would only save 17% of the total NCMesh memory if
    // 4-element arrays were used (e.g. through templates); we thus prefer to
    // keep the code as simple as possible.
@@ -476,7 +484,7 @@ int NCMesh::NewHexahedron(int n0, int n1, int n2, int n3,
    el.node[4] = n4, el.node[5] = n5, el.node[6] = n6, el.node[7] = n7;
 
    // get faces and assign face attributes
-   Face* f[6];
+   Face* f[MaxElemFaces];
    const GeomInfo &gi_hex = GI[Geometry::CUBE];
    for (int i = 0; i < gi_hex.nf; i++)
    {
@@ -549,6 +557,35 @@ int NCMesh::NewTetrahedron(int n0, int n1, int n2, int n3, int attr,
 
    return new_id;
 }
+int NCMesh::NewPyramid(int n0, int n1, int n2, int n3, int n4, int attr,
+                       int fattr0, int fattr1, int fattr2, int fattr3,
+                       int fattr4)
+{
+   // create new element, initialize nodes
+   int new_id = AddElement(Element(Geometry::PYRAMID, attr));
+   Element &el = elements[new_id];
+
+   el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
+   el.node[4] = n4;
+
+   // get faces and assign face attributes
+   Face* f[5];
+   const GeomInfo &gi_pyr = GI[Geometry::PYRAMID];
+   for (int i = 0; i < gi_pyr.nf; i++)
+   {
+      const int* fv = gi_pyr.faces[i];
+      f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
+                       el.node[fv[2]], el.node[fv[3]]);
+   }
+
+   f[0]->attribute = fattr0;
+   f[1]->attribute = fattr1;
+   f[2]->attribute = fattr2;
+   f[3]->attribute = fattr3;
+   f[4]->attribute = fattr4;
+
+   return new_id;
+}
 
 int NCMesh::NewQuadrilateral(int n0, int n1, int n2, int n3,
                              int attr,
@@ -895,7 +932,7 @@ void NCMesh::RefineElement(int elem, char ref_type)
       char remaining = ref_type & ~el.ref_type;
 
       // do the remaining splits on the children
-      for (int i = 0; i < 8; i++)
+      for (int i = 0; i < MaxElemChildren; i++)
       {
          if (el.child[i] >= 0) { RefineElement(el.child[i], remaining); }
       }
@@ -912,11 +949,11 @@ void NCMesh::RefineElement(int elem, char ref_type)
    int* no = el.node;
    int attr = el.attribute;
 
-   int child[8];
-   for (int i = 0; i < 8; i++) { child[i] = -1; }
+   int child[MaxElemChildren];
+   for (int i = 0; i < MaxElemChildren; i++) { child[i] = -1; }
 
    // get parent's face attributes
-   int fa[6];
+   int fa[MaxElemFaces];
    GeomInfo& gi = GI[el.Geom()];
    for (int i = 0; i < gi.nf; i++)
    {
@@ -1404,6 +1441,65 @@ void NCMesh::RefineElement(int elem, char ref_type)
                                    -1, -1, fa[3], -1);
       }
    }
+   else if (el.Geom() == Geometry::PYRAMID)
+   {
+      // Pyramid vertex numbering:
+      //
+      //    4
+      //     + \_                       Faces: 0 bottom (3,2,1,0)
+      //     |\\_ \_                           1 front (0, 1, 4)
+      //     ||  \_ \__                        2 right (1, 2, 4)
+      //     | \   \_   \__                    3 back (2, 3, 4)
+      //     |  +____\_ ____\                  4 left (3, 0, 4)
+      //     | /3      \_   2   Z  Y
+      //     |/          \ /    | /
+      //     +------------+     *--X
+      //    0              1
+
+      ref_type = Refinement::XYZ; // for consistence
+
+      int mid01 = GetMidEdgeNode(no[0], no[1]);
+      int mid12 = GetMidEdgeNode(no[1], no[2]);
+      int mid23 = GetMidEdgeNode(no[2], no[3]);
+      int mid03 = GetMidEdgeNode(no[0], no[3]);
+      int mid04 = GetMidEdgeNode(no[0], no[4]);
+      int mid14 = GetMidEdgeNode(no[1], no[4]);
+      int mid24 = GetMidEdgeNode(no[2], no[4]);
+      int mid34 = GetMidEdgeNode(no[3], no[4]);
+      int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid03);
+
+      child[0] = NewPyramid(no[0], mid01, midf0, mid03, mid04,
+                            attr, fa[0], fa[1], -1, -1, fa[4]);
+
+      child[1] = NewPyramid(mid01, no[1], mid12, midf0, mid14,
+                            attr, fa[0], fa[1], fa[2], -1, -1);
+
+      child[2] = NewPyramid(midf0, mid12, no[2], mid23, mid24,
+                            attr, fa[0], -1, fa[2], fa[3], -1);
+
+      child[3] = NewPyramid(mid03, midf0, mid23, no[3], mid34,
+                            attr, fa[0], -1, -1, fa[3], fa[4]);
+
+      child[4] = NewPyramid(mid24, mid14, mid04, mid34, midf0,
+                            attr, -1, -1, -1, -1, -1);
+
+      child[5] = NewPyramid(mid04, mid14, mid24, mid34, no[4],
+                            attr, -1, fa[1], fa[2], fa[3], fa[4]);
+
+      child[6] = NewTetrahedron(mid01, midf0, mid04, mid14,
+                                attr, -1, -1, -1, fa[1]);
+
+      child[7] = NewTetrahedron(midf0, mid14, mid12, mid24,
+                                attr, -1, -1, fa[2], -1);
+
+      child[8] = NewTetrahedron(midf0, mid23, mid34, mid24,
+                                attr, -1, -1, fa[3], -1);
+
+      child[9] = NewTetrahedron(mid03, mid04, midf0, mid34,
+                                attr, -1, fa[4], -1, -1);
+
+      CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid03, midf0);
+   }
    else if (el.Geom() == Geometry::SQUARE)
    {
       ref_type &= 0x3; // ignore Z bit
@@ -1486,20 +1582,20 @@ void NCMesh::RefineElement(int elem, char ref_type)
    }
 
    // start using the nodes of the children, create edges & faces
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       ReferenceElement(child[i]);
    }
 
-   int buf[6];
-   Array<int> parentFaces(buf, 6);
+   int buf[MaxElemFaces];
+   Array<int> parentFaces(buf, MaxElemFaces);
    parentFaces.SetSize(0);
 
    // sign off of all nodes of the parent, clean up unused nodes, but keep faces
    UnreferenceElement(elem, parentFaces);
 
    // register the children in their faces
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       RegisterFaces(child[i]);
    }
@@ -1508,7 +1604,7 @@ void NCMesh::RefineElement(int elem, char ref_type)
    DeleteUnusedFaces(parentFaces);
 
    // make the children inherit our rank; set the parent element
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       Element &ch = elements[child[i]];
       ch.rank = el.rank;
@@ -1556,6 +1652,8 @@ void NCMesh::Refine(const Array<Refinement>& refinements)
 #if defined(MFEM_DEBUG) && !defined(MFEM_USE_MPI)
    mfem::out << "Refined " << refinements.Size() << " + " << nforced
              << " elements" << std::endl;
+#else
+   MFEM_CONTRACT_VAR(nforced);
 #endif
 
    ref_stack.DeleteAll();
@@ -1586,6 +1684,12 @@ int NCMesh::RetrieveNode(const Element &el, int index)
          ch = el.child[ch];
          break;
 
+      case Geometry::PYRAMID:
+         ch = pyramid_deref_table[el.ref_type - 1][index];
+         MFEM_ASSERT(ch != -1, "");
+         ch = el.child[ch];
+         break;
+
       case Geometry::SQUARE:
          ch = el.child[quad_deref_table[el.ref_type - 1][index]];
          break;
@@ -1608,11 +1712,11 @@ void NCMesh::DerefineElement(int elem)
    Element &el = elements[elem];
    if (!el.ref_type) { return; }
 
-   int child[8];
+   int child[MaxElemChildren];
    std::memcpy(child, el.child, sizeof(child));
 
    // first make sure that all children are leaves, derefine them if not
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       if (elements[child[i]].ref_type)
       {
@@ -1620,10 +1724,10 @@ void NCMesh::DerefineElement(int elem)
       }
    }
 
-   int faces_attribute[6];
+   int faces_attribute[MaxElemFaces];
    int ref_type_key = el.ref_type - 1;
 
-   for (int i = 0; i < 8; i++) { el.node[i] = -1; }
+   for (int i = 0; i < MaxElemNodes; i++) { el.node[i] = -1; }
 
    // retrieve original corner nodes and face attributes from the children
    if (el.Geom() == Geometry::CUBE)
@@ -1678,6 +1782,34 @@ void NCMesh::DerefineElement(int elem)
                               ->attribute;
       }
    }
+   else if (el.Geom() == Geometry::PYRAMID)
+   {
+      MFEM_ASSERT(pyramid_deref_table[ref_type_key][0] != -1,
+                  "invalid pyramid refinement");
+      constexpr int nb_pyramid_childs = 5;
+      for (int i = 0; i < nb_pyramid_childs; i++)
+      {
+         const int child_local_index = pyramid_deref_table[ref_type_key][i];
+         const int child_global_index = child[child_local_index];
+         Element &ch = elements[child_global_index];
+         el.node[i] = ch.node[i];
+      }
+      el.node[5] = el.node[6] = el.node[7] = -1;
+
+
+      constexpr int nb_pyramid_faces = 5;
+      for (int i = 0; i < nb_pyramid_faces; i++)
+      {
+         const int child_local_index = pyramid_deref_table[ref_type_key]
+                                       [i + nb_pyramid_childs];
+         const int child_global_index = child[child_local_index];
+         Element &ch = elements[child_global_index];
+         const int* fv = GI[el.Geom()].faces[i];
+         faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
+                                         ch.node[fv[2]], ch.node[fv[3]])
+                              ->attribute;
+      }
+   }
    else if (el.Geom() == Geometry::TETRAHEDRON)
    {
       for (int i = 0; i < 4; i++)
@@ -1745,13 +1877,13 @@ void NCMesh::DerefineElement(int elem)
    // sign in to all nodes
    ReferenceElement(elem);
 
-   int buf[8*6];
-   Array<int> childFaces(buf, 8*6);
+   int buf[MaxElemChildren*MaxElemFaces];
+   Array<int> childFaces(buf, MaxElemChildren*MaxElemFaces);
    childFaces.SetSize(0);
 
    // delete children, determine rank
    el.rank = std::numeric_limits<int>::max();
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       el.rank = std::min(el.rank, elements[child[i]].rank);
       UnreferenceElement(child[i], childFaces);
@@ -1775,7 +1907,7 @@ void NCMesh::CollectDerefinements(int elem, Array<Connection> &list)
    if (!el.ref_type) { return; }
 
    int total = 0, ref = 0, ghost = 0;
-   for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
    {
       total++;
       Element &ch = elements[el.child[i]];
@@ -1787,7 +1919,7 @@ void NCMesh::CollectDerefinements(int elem, Array<Connection> &list)
    {
       // can be derefined, add to list
       int next_row = list.Size() ? (list.Last().from + 1) : 0;
-      for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+      for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
       {
          Element &ch = elements[el.child[i]];
          list.Append(Connection(next_row, ch.index));
@@ -1795,7 +1927,7 @@ void NCMesh::CollectDerefinements(int elem, Array<Connection> &list)
    }
    else
    {
-      for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+      for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
       {
          CollectDerefinements(el.child[i], list);
       }
@@ -1905,7 +2037,7 @@ void NCMesh::SetDerefMatrixCodes(int parent, Array<int> &fine_coarse)
 {
    // encode the ref_type and child number for GetDerefinementTransforms()
    Element &prn = elements[parent];
-   for (int i = 0; i < 8 && prn.child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && prn.child[i] >= 0; i++)
    {
       Element &ch = elements[prn.child[i]];
       if (ch.index >= 0)
@@ -1977,7 +2109,7 @@ void NCMesh::CollectLeafElements(int elem, int state, Array<int> &ghosts,
       }
       else // no space filling curve tables yet for remaining cases
       {
-         for (int i = 0; i < 8; i++)
+         for (int i = 0; i < MaxElemChildren; i++)
          {
             if (el.child[i] >= 0)
             {
@@ -2231,7 +2363,8 @@ void NCMesh::InitRootState(int root_count)
       if (v_in < 0) { v_in = 0; }
 
       // determine which nodes are shared with the next element
-      bool shared[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
+      bool shared[MaxElemNodes];
+      for (int ni = 0; ni < MaxElemNodes; ++ni) { shared[ni] = 0; }
       if (i+1 < root_count)
       {
          Element &next = elements[i+1];
@@ -2265,6 +2398,7 @@ mfem::Element* NCMesh::NewMeshElement(int geom) const
    {
       case Geometry::CUBE: return new mfem::Hexahedron;
       case Geometry::PRISM: return new mfem::Wedge;
+      case Geometry::PYRAMID: return new mfem::Pyramid;
       case Geometry::TETRAHEDRON: return new mfem::Tetrahedron;
       case Geometry::SQUARE: return new mfem::Quadrilateral;
       case Geometry::TRIANGLE: return new mfem::Triangle;
@@ -2347,7 +2481,8 @@ void NCMesh::GetMeshComponents(Mesh &mesh) const
          if (face->Boundary())
          {
             if ((nc_elem.geom == Geometry::CUBE) ||
-                (nc_elem.geom == Geometry::PRISM && nfv == 4))
+                ((nc_elem.geom == Geometry::PRISM ||
+                  nc_elem.geom == Geometry::PYRAMID) && nfv == 4))
             {
                auto* quad = (Quadrilateral*) mesh.NewElement(Geometry::SQUARE);
                quad->SetAttribute(face->attribute);
@@ -2358,6 +2493,7 @@ void NCMesh::GetMeshComponents(Mesh &mesh) const
                mesh.boundary.Append(quad);
             }
             else if (nc_elem.geom == Geometry::PRISM ||
+                     nc_elem.geom == Geometry::PYRAMID ||
                      nc_elem.geom == Geometry::TETRAHEDRON)
             {
                MFEM_ASSERT(nfv == 3, "");
@@ -2603,7 +2739,7 @@ bool NCMesh::TriFaceSplit(int v1, int v2, int v3, int mid[3]) const
 
 int NCMesh::find_node(const Element &el, int node)
 {
-   for (int i = 0; i < 8; i++)
+   for (int i = 0; i < MaxElemNodes; i++)
    {
       if (el.node[i] == node) { return i; }
    }
@@ -3625,7 +3761,7 @@ void NCMesh::FindNeighbors(int elem, Array<int> &neighbors,
       }
       else
       {
-         for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+         for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
          {
             stack.Append(el.child[i]);
          }
@@ -3737,8 +3873,8 @@ int NCMesh::GetVertexRootCoord(int elem, RefCoord coord[3]) const
       MFEM_ASSERT(pa.ref_type, "internal error");
 
       int ch = 0;
-      while (ch < 8 && pa.child[ch] != elem) { ch++; }
-      MFEM_ASSERT(ch < 8, "internal error");
+      while (ch < MaxElemChildren && pa.child[ch] != elem) { ch++; }
+      MFEM_ASSERT(ch < MaxElemChildren, "internal error");
 
       MFEM_ASSERT(geom_parent[el.Geom()], "unsupported geometry");
       const RefTrf &tr = geom_parent[el.Geom()][(int) pa.ref_type][ch];
@@ -3768,6 +3904,11 @@ static bool RefPointInside(Geometry::Type geom, const RefCoord pt[3])
          return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE) &&
                 (pt[2] >= 0) && (pt[2] <= T_ONE);
 
+      case Geometry::PYRAMID:
+         return (pt[0] >= 0) && (pt[1] >= 0) && (pt[2] >= 0.0) &&
+                (pt[0] + pt[2] <= T_ONE) && (pt[1] + pt[2] <= T_ONE) &&
+                (pt[2] <= T_ONE);
+
       default:
          MFEM_ABORT("unsupported geometry");
          return false;
@@ -3785,7 +3926,7 @@ void NCMesh::CollectIncidentElements(int elem, const RefCoord coord[3],
    }
 
    RefCoord tcoord[3];
-   for (int ch = 0; ch < 8 && el.child[ch] >= 0; ch++)
+   for (int ch = 0; ch < MaxElemChildren && el.child[ch] >= 0; ch++)
    {
       const RefTrf &tr = geom_child[el.Geom()][(int) el.ref_type][ch];
       tr.Apply(coord, tcoord);
@@ -3856,6 +3997,10 @@ NCMesh::PointMatrix NCMesh::pm_prism_identity(
    Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0),
    Point(0, 0, 1), Point(1, 0, 1), Point(0, 1, 1)
 );
+NCMesh::PointMatrix NCMesh::pm_pyramid_identity(
+   Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0),
+   Point(0, 1, 0), Point(0, 0, 1)
+);
 NCMesh::PointMatrix NCMesh::pm_hex_identity(
    Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
    Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
@@ -3870,6 +4015,7 @@ const NCMesh::PointMatrix& NCMesh::GetGeomIdentity(Geometry::Type geom)
       case Geometry::SQUARE:      return pm_quad_identity;
       case Geometry::TETRAHEDRON: return pm_tet_identity;
       case Geometry::PRISM:       return pm_prism_identity;
+      case Geometry::PYRAMID:     return pm_pyramid_identity;
       case Geometry::CUBE:        return pm_hex_identity;
       default:
          MFEM_ABORT("unsupported geometry " << geom);
@@ -4174,6 +4320,55 @@ void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
             }
          }
       }
+      else if (geom == Geometry::PYRAMID)
+      {
+         Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
+         Point mid03(pm(0), pm(3)), mid12(pm(1), pm(2));
+         Point mid04(pm(0), pm(4)), mid14(pm(1), pm(4));
+         Point mid24(pm(2), pm(4)), mid34(pm(3), pm(4));
+         Point midf0(mid23, mid12, mid01, mid03);
+
+         if (child == 0)        // Pyramid
+         {
+            pm = PointMatrix(pm(0), mid01, midf0, mid03, mid04);
+         }
+         else if (child == 1)   // Pyramid
+         {
+            pm = PointMatrix(mid01, pm(1), mid12, midf0, mid14);
+         }
+         else if (child == 2)   // Pyramid
+         {
+            pm = PointMatrix(midf0, mid12, pm(2), mid23, mid24);
+         }
+         else if (child == 3)   // Pyramid
+         {
+            pm = PointMatrix(mid03, midf0, mid23, pm(3), mid34);
+         }
+         else if (child == 4)   // Pyramid
+         {
+            pm = PointMatrix(mid24, mid14, mid04, mid34, midf0);
+         }
+         else if (child == 5)   // Pyramid
+         {
+            pm = PointMatrix(mid04, mid14, mid24, mid34, pm(4));
+         }
+         else if (child == 6)   // Tet
+         {
+            pm = PointMatrix(mid01, midf0, mid04, mid14);
+         }
+         else if (child == 7)   // Tet
+         {
+            pm = PointMatrix(midf0, mid14, mid12, mid24);
+         }
+         else if (child == 8)   // Tet
+         {
+            pm = PointMatrix(midf0, mid23, mid34, mid24);
+         }
+         else if (child == 9)   // Tet
+         {
+            pm = PointMatrix(mid03, mid04, midf0, mid34);
+         }
+      }
       else if (geom == Geometry::TETRAHEDRON)
       {
          Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid02(pm(2), pm(0));
@@ -4333,7 +4528,7 @@ void NCMesh::TraverseRefinements(int elem, int coarse_index,
       ref_path.push_back(el.ref_type);
       ref_path.push_back(0);
 
-      for (int i = 0; i < 8; i++)
+      for (int i = 0; i < MaxElemChildren; i++)
       {
          if (el.child[i] >= 0)
          {
@@ -5091,14 +5286,14 @@ void NCMesh::CountSplits(int elem, int splits[3]) const
    const int* node = el.node;
    GeomInfo& gi = GI[el.Geom()];
 
-   int elevel[12];
+   int elevel[MaxElemEdges];
    for (int i = 0; i < gi.ne; i++)
    {
       const int* ev = gi.edges[i];
       elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
    }
 
-   int flevel[6][2];
+   int flevel[MaxElemFaces][2];
    if (Dim >= 3)
    {
       for (int i = 0; i < gi.nf; i++)
@@ -5142,6 +5337,18 @@ void NCMesh::CountSplits(int elem, int splits[3]) const
       splits[2] = max6(flevel[2][1], flevel[3][1], flevel[4][1],
                        elevel[6], elevel[7], elevel[8]);
    }
+   else if (el.Geom() == Geometry::PYRAMID)
+   {
+      splits[0] = std::max(
+                     max6(flevel[0][0], flevel[1][0], 0,
+                          flevel[2][0], flevel[3][0], flevel[4][0]),
+                     max8(elevel[0], elevel[1], elevel[2],
+                          elevel[3], elevel[4], elevel[5],
+                          elevel[6], elevel[7]));
+
+      splits[1] = splits[0];
+      splits[2] = splits[0];
+   }
    else if (el.Geom() == Geometry::TETRAHEDRON)
    {
       splits[0] = std::max(
@@ -5404,7 +5611,8 @@ void NCMesh::Print(std::ostream &os) const
       "# SQUARE      = 3\n"
       "# TETRAHEDRON = 4\n"
       "# CUBE        = 5\n"
-      "# PRISM       = 6\n";
+      "# PRISM       = 6\n"
+      "# PYRAMID     = 7\n";
 
    os << "\ndimension\n" << Dim << "\n";
 
@@ -5425,7 +5633,7 @@ void NCMesh::Print(std::ostream &os) const
       if (el.parent == -2) { os << "-1\n"; continue; } // unused element
 
       os << int(el.geom) << " " << int(el.ref_type);
-      for (int j = 0; j < 8 && el.node[j] >= 0; j++)
+      for (int j = 0; j < MaxElemNodes && el.node[j] >= 0; j++)
       {
          os << " " << el.node[j];
       }
@@ -5482,7 +5690,7 @@ void NCMesh::InitRootElements()
       Element &el = elements[i];
       if (el.ref_type)
       {
-         for (int j = 0; j < 8 && el.child[j] >= 0; j++)
+         for (int j = 0; j < MaxElemChildren && el.child[j] >= 0; j++)
          {
             int child = el.child[j];
             MFEM_VERIFY(child < elements.Size(), "invalid mesh file: "
@@ -5704,7 +5912,7 @@ void NCMesh::CopyElements(int elem,
    Element &el = elements[elem];
    if (el.ref_type)
    {
-      for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+      for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
       {
          int old_id = el.child[i];
          // here we know 'free_element_ids' is empty
@@ -5984,9 +6192,9 @@ void NCMesh::Trim()
    // maybe also of 'nodes' and 'faces'.
 }
 
-std::size_t NCMesh::NCList::MemoryUsage() const
+long NCMesh::NCList::MemoryUsage() const
 {
-   std::size_t pm_size = 0;
+   int pm_size = 0;
    for (int i = 0; i < Geometry::NumGeom; i++)
    {
       for (int j = 0; j < point_matrices[i].Size(); i++)
@@ -6002,9 +6210,9 @@ std::size_t NCMesh::NCList::MemoryUsage() const
           pm_size;
 }
 
-std::size_t CoarseFineTransformations::MemoryUsage() const
+long CoarseFineTransformations::MemoryUsage() const
 {
-   std::size_t mem = embeddings.MemoryUsage();
+   long mem = embeddings.MemoryUsage();
    for (int i = 0; i < Geometry::NumGeom; i++)
    {
       mem += point_matrices[i].MemoryUsage();
@@ -6012,7 +6220,7 @@ std::size_t CoarseFineTransformations::MemoryUsage() const
    return mem;
 }
 
-std::size_t NCMesh::MemoryUsage() const
+long NCMesh::MemoryUsage() const
 {
    return nodes.MemoryUsage() +
           faces.MemoryUsage() +
@@ -6073,7 +6281,7 @@ void NCMesh::DebugLeafOrder(std::ostream &os) const
       {
          double sum = 0.0;
          int count = 0;
-         for (int k = 0; k < 8; k++)
+         for (int k = 0; k < MaxElemNodes; k++)
          {
             if (elem->node[k] >= 0)
             {

mesh/ncmesh.hpp

@@ -84,7 +84,7 @@ struct CoarseFineTransformations
 
    void Clear();
    bool IsInitialized() const;
-   std::size_t MemoryUsage() const;
+   long MemoryUsage() const;
 
    MFEM_DEPRECATED
    void GetCoarseToFineMap(const Mesh &fine_mesh, Table &coarse_to_fine) const
@@ -240,7 +240,7 @@ public:
       void Clear();
       bool Empty() const { return !conforming.Size() && !masters.Size(); }
       long TotalSize() const;
-      std::size_t MemoryUsage() const;
+      long MemoryUsage() const;
 
       const MeshId& LookUp(int index, int *type = NULL) const;
 
@@ -390,7 +390,7 @@ public:
    virtual void Trim();
 
    /// Return total number of bytes allocated.
-   std::size_t MemoryUsage() const;
+   long MemoryUsage() const;
 
    int PrintMemoryDetail() const;
 
@@ -421,6 +421,15 @@ protected: // implementation
    int Geoms; ///< bit mask of element geometries present, see InitGeomFlags()
    bool Legacy; ///< true if the mesh was loaded from the legacy v1.1 format
 
+   static const int MaxElemNodes =
+      8;       ///< Number of nodes of an element can have
+   static const int MaxElemEdges =
+      12;      ///< Number of edges of an element can have
+   static const int MaxElemFaces =
+      6;       ///< Number of faces of an element can have
+   static const int MaxElemChildren =
+      10;      ///< Number of children of an element can have
+
    /** A Node can hold a vertex, an edge, or both. Elements directly point to
        their corner nodes, but edge nodes also exist and can be accessed using
        a hash-table given their two end-point node IDs. All nodes can be
@@ -482,8 +491,8 @@ protected: // implementation
       int attribute;
       union
       {
-         int node[8];  ///< element corners (if ref_type == 0)
-         int child[8]; ///< 2-8 children (if ref_type != 0)
+         int node[MaxElemNodes];  ///< element corners (if ref_type == 0)
+         int child[MaxElemChildren]; ///< 2-10 children (if ref_type != 0)
       };
       int parent; ///< parent element, -1 if this is a root element, -2 if free'd
 
@@ -586,6 +595,9 @@ protected: // implementation
    /// Return true if the mesh contains prism elements.
    bool HavePrisms() const { return Geoms & (1 << Geometry::PRISM); }
 
+   /// Return true if the mesh contains pyramid elements.
+   bool HavePyramids() const { return Geoms & (1 << Geometry::PYRAMID); }
+
    /// Return true if the mesh contains tetrahedral elements.
    bool HaveTets() const   { return Geoms & (1 << Geometry::TETRAHEDRON); }
 
@@ -642,6 +654,10 @@ protected: // implementation
    int NewTetrahedron(int n0, int n1, int n2, int n3, int attr,
                       int fattr0, int fattr1, int fattr2, int fattr3);
 
+   int NewPyramid(int n0, int n1, int n2, int n3, int n4, int attr,
+                  int fattr0, int fattr1, int fattr2, int fattr3,
+                  int fattr4);
+
    int NewQuadrilateral(int n0, int n1, int n2, int n3, int attr,
                         int eattr0, int eattr1, int eattr2, int eattr3);
 
@@ -845,7 +861,7 @@ protected: // implementation
    struct PointMatrix
    {
       int np;
-      Point points[8];
+      Point points[MaxElemNodes];
 
       PointMatrix() : np(0) {}
 
@@ -858,6 +874,13 @@ protected: // implementation
       PointMatrix(const Point& p0, const Point& p1, const Point& p2, const Point& p3)
       { np = 4; points[0] = p0; points[1] = p1; points[2] = p2; points[3] = p3; }
 
+      PointMatrix(const Point& p0, const Point& p1, const Point& p2,
+                  const Point& p3, const Point& p4)
+      {
+         np = 5;
+         points[0] = p0; points[1] = p1; points[2] = p2;
+         points[3] = p3; points[4] = p4;
+      }
       PointMatrix(const Point& p0, const Point& p1, const Point& p2,
                   const Point& p3, const Point& p4, const Point& p5)
       {
@@ -887,6 +910,7 @@ protected: // implementation
    static PointMatrix pm_quad_identity;
    static PointMatrix pm_tet_identity;
    static PointMatrix pm_prism_identity;
+   static PointMatrix pm_pyramid_identity;
    static PointMatrix pm_hex_identity;
 
    static const PointMatrix& GetGeomIdentity(Geometry::Type geom);
@@ -976,9 +1000,9 @@ protected: // implementation
    struct GeomInfo
    {
       int nv, ne, nf;   // number of: vertices, edges, faces
-      int edges[12][2]; // edge vertices (up to 12 edges)
-      int faces[6][4];  // face vertices (up to 6 faces)
-      int nfv[6];       // number of face vertices
+      int edges[MaxElemEdges][2]; // edge vertices (up to 12 edges)
+      int faces[MaxElemFaces][4];  // face vertices (up to 6 faces)
+      int nfv[MaxElemFaces];       // number of face vertices
 
       bool initialized;
       GeomInfo() : initialized(false) {}

mesh/ncmesh_tables.hpp

@@ -19,6 +19,10 @@ const int ref_type_num_children[8] = { 0, 2, 2, 4, 2, 4, 4, 8 };
 
 
 // derefinement tables
+// The first n numbers in each line are the refined elements that contain
+// the vertices of the parent element.  The next m numbers in each line
+// are the refined elements that contain the faces attributes of the parent
+// element.
 
 const int quad_deref_table[3][4 + 4] =
 {
@@ -49,6 +53,16 @@ const int prism_deref_table[7][6 + 5] =
    { 0, 1, 2, 4, 5, 6, /**/  0, 5, 0, 5, 0 }  // 7 - iso
 };
 
+const int pyramid_deref_table[7][5 + 5] =
+{
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 1
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 2
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 3
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 4
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 5
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 6
+   { 0, 1, 2, 3, 5, /**/  0, 5, 5, 5, 5 }  // 7 - iso
+};
 
 // child ordering tables
 
@@ -140,8 +154,17 @@ const RefCoord prism_corners[6][3] =
    {    0, T_ONE, T_ONE}
 };
 
+const RefCoord pyramid_corners[5][3] =
+{
+   {    0,     0,     0},
+   {T_ONE,     0,     0},
+   {T_ONE, T_ONE,     0},
+   {    0, T_ONE,     0},
+   {    0,     0, T_ONE}
+};
+
 typedef RefCoord RefPoint[3];
-const RefPoint* geom_corners[7] =
+const RefPoint* geom_corners[8] =
 {
    NULL, // point
    NULL, // segment
@@ -149,7 +172,8 @@ const RefPoint* geom_corners[7] =
    quad_corners,
    NULL, // tetrahedron
    hex_corners,
-   prism_corners
+   prism_corners,
+   pyramid_corners
 };
 
 // reference domain transform: 3 scales, 3 translations

mesh/pmesh.cpp

@@ -6099,8 +6099,10 @@ void ParMesh::GetBoundingBox(Vector &gp_min, Vector &gp_max, int ref)
    gp_min.SetSize(sdim);
    gp_max.SetSize(sdim);
 
-   MPI_Allreduce(p_min.GetData(), gp_min, sdim, MPI_DOUBLE, MPI_MIN, MyComm);
-   MPI_Allreduce(p_max.GetData(), gp_max, sdim, MPI_DOUBLE, MPI_MAX, MyComm);
+   MPI_Allreduce(p_min.GetData(), gp_min.GetData(), sdim, MPI_DOUBLE,
+                 MPI_MIN, MyComm);
+   MPI_Allreduce(p_max.GetData(), gp_max.GetData(), sdim, MPI_DOUBLE,
+                 MPI_MAX, MyComm);
 }
 
 void ParMesh::GetCharacteristics(double &gh_min, double &gh_max,

mesh/pncmesh.cpp

@@ -2798,14 +2798,6 @@ static std::size_t arrays_memory_usage(const Array<Type> (&arrays)[Size])
 
 std::size_t ParNCMesh::MemoryUsage(bool with_base) const
 {
-   std::size_t total_groups_owners = 0;
-   for (int i = 0; i < 3; i++)
-   {
-      total_groups_owners += entity_owner[i].MemoryUsage() +
-                             entity_pmat_group[i].MemoryUsage() +
-                             entity_index_rank[i].MemoryUsage();
-   }
-
    return (with_base ? NCMesh::MemoryUsage() : 0) +
           GroupsMemoryUsage() +
           arrays_memory_usage(entity_owner) +

mesh/submesh/psubmesh.cpp

@@ -248,8 +248,11 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
       Transfer(*pn, *n);
    }
 
-   el_to_edge = new Table;
-   NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+   if (Dim > 1)
+   {
+      el_to_edge = new Table;
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+   }
 
    SetAttributes();
    Finalize();
@@ -811,4 +814,4 @@ ParTransferMap ParSubMesh::CreateTransferMap(const ParGridFunction &src,
 
 } // namespace mfem
 
-#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MPI

mesh/submesh/ptransfermap.cpp

@@ -183,7 +183,7 @@ void ParTransferMap::CommunicateSharedVdofs(Vector &f) const
    }
 
    // TODO: do the reduce only on dofs of interest
-   root_gc_->Reduce<double>(f, GroupCommunicator::Sum);
+   root_gc_->Reduce<double>(f.HostReadWrite(), GroupCommunicator::Sum);
 
    // Indices that were set from this rank or other ranks have been summed up
    // and therefore need to be "averaged". Note that this results in the exact
@@ -211,7 +211,7 @@ void ParTransferMap::CommunicateSharedVdofs(Vector &f) const
       }
    }
 
-   root_gc_->Bcast<double>(f);
+   root_gc_->Bcast<double>(f.HostReadWrite());
 }
 
 ParTransferMap::~ParTransferMap()

mesh/submesh/submesh_utils.cpp

@@ -95,6 +95,8 @@ void BuildVdofToVdofMap(const FiniteElementSpace& subfes,
    auto *m = subfes.GetMesh();
    vdof_to_vdof_map.SetSize(subfes.GetVSize());
 
+   const int vdim = parentfes.GetVDim();
+
    IntegrationPointTransformation Tr;
    DenseMatrix T;
    Array<int> z1;
@@ -137,14 +139,20 @@ void BuildVdofToVdofMap(const FiniteElementSpace& subfes,
 
             parentfes.GetElementVDofs(parent_volel_id, z1);
 
-            parent_vdofs.SetSize(parentfes.GetVDim() * T.Height());
+            parent_vdofs.SetSize(vdim * T.Height());
             for (int j = 0; j < T.Height(); j++)
             {
-               for (int k = 0; k < parentfes.GetVDim() * T.Width(); k++)
+               for (int k = 0; k < T.Width(); k++)
                {
                   if (T(j, k) != 0.0)
                   {
-                     parent_vdofs[j] = z1[static_cast<int>(k)];
+                     for (int vd=0; vd<vdim; vd++)
+                     {
+                        int sub_vdof = j + T.Height() * vd;
+                        int parent_vdof = k + T.Width() * vd;
+                        parent_vdofs[sub_vdof] =
+                           z1[static_cast<int>(parent_vdof)];
+                     }
                   }
                }
             }

miniapps/common/fem_extras.cpp

@@ -58,7 +58,7 @@ RT_FESpace::~RT_FESpace()
 
 void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
                    Mesh &mesh, const char *title,
-                   int x, int y, int w, int h, const char * keys, bool vec)
+                   int x, int y, int w, int h, const char * keys)
 {
    bool newly_opened = false;
    int connection_failed;
@@ -71,7 +71,7 @@ void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
          sock.precision(8);
          newly_opened = true;
       }
-      sock << "solution\n";
+      sock << "mesh\n";
 
       mesh.Print(sock);
 
@@ -81,8 +81,6 @@ void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
               << "window_geometry "
               << x << " " << y << " " << w << " " << h << "\n";
          if ( keys ) { sock << "keys " << keys << "\n"; }
-         else { sock << "keys maaAc\n"; }
-         if ( vec ) { sock << "vvv"; }
          sock << endl;
       }