minor

AI policy

.github/actions/sanitize/mfem/action.yml

@@ -25,7 +25,7 @@ runs:
   steps:
     - uses: ./.github/actions/sanitize/config
 
-    - uses: actions/cache@v4
+    - uses: actions/cache@v5
       if: ${{env.DEBUG == 'true'}}
       id: debug
       with:

.github/actions/sanitize/restore/action.yml

@@ -36,7 +36,7 @@ runs:
   steps:
     - uses: ./.github/actions/sanitize/config
 
-    - uses: actions/cache@v4
+    - uses: actions/cache@v5
       if: ${{env.DEBUG == 'true' && inputs.cache-skip != 'true'}}
       id: debug
       with:

.github/actions/sanitize/setup/action.yml

@@ -23,7 +23,7 @@ inputs:
 runs:
   using: 'composite'
   steps:
-    - uses: actions/cache/restore@v4  # Cache for LLVM libcxx
+    - uses: actions/cache/restore@v5  # Cache for LLVM libcxx
       with:
         path: ${{env.LLVM_DIR}}
         fail-on-cache-miss: true
@@ -32,14 +32,14 @@ runs:
     - uses: ./.github/actions/sanitize/mpi
       if: ${{inputs.par == 'true'}}
 
-    - uses: actions/cache/restore@v4  # Cache for Hypre
+    - uses: actions/cache/restore@v5  # Cache for Hypre
       if: ${{inputs.par == 'true'}}
       with:
         path: ${{env.HYPRE_DIR}}
         fail-on-cache-miss: true
         key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5
 
-    - uses: actions/cache/restore@v4  # Cache for Metis
+    - uses: actions/cache/restore@v5  # Cache for Metis
       if: ${{inputs.par == 'true'}}
       with:
         path: ${{env.METIS_DIR}}
@@ -51,13 +51,13 @@ runs:
       run: ln -s -f ${{env.HYPRE_DIR}} hypre && ln -s -f ${{env.METIS_DIR}} metis-4.0
       shell: bash
 
-    - uses: actions/cache/restore@v4  # Cache for LSAN suppression file
+    - uses: actions/cache/restore@v5  # Cache for LSAN suppression file
       with:
         path: ${{env.LSAN_DIR}}
         fail-on-cache-miss: true
         key: build-lsan-suppression-file
 
-    - uses: actions/checkout@v4  # Checkout the repository
+    - uses: actions/checkout@v6  # Checkout the repository
       with:
         path: mfem
         # ref: ${{env.BRANCH}}

.github/workflows/build-container.yml

@@ -43,7 +43,7 @@ jobs:
           remove-docker-images: 'true'
 
       - name: Checkout
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
 
       # It's easier to reference named variables than indexes of the matrix
       - name: Set Environment

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

@@ -153,7 +153,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@v4
+        uses: actions/checkout@v6
         with:
           path: ${{ env.MFEM_TOP_DIR }}
           # Fetch the complete history for codecov to access commits ID
@@ -225,7 +225,7 @@ jobs:
       - name: cache hypre
         id: hypre-cache
         if: matrix.mpi == 'par'
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{ env.HYPRE_TOP_DIR }}
           key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-${{ matrix.precision }}-v2.5
@@ -255,7 +255,7 @@ jobs:
       - name: cache metis
         id: metis-cache
         if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{ env.METIS_TOP_DIR }}
           key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5
@@ -270,7 +270,7 @@ jobs:
       - name: cache vcpkg (Windows)
         id: vcpkg-cache
         if: matrix.os == 'windows-latest'
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: vcpkg_cache
           key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
@@ -295,7 +295,8 @@ jobs:
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
           brew install enzyme
-          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
+          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required.*:.*\(llvm[^ ]*\).*/\1/p')
+          echo "ENZYME_LLVM=$ENZYME_LLVM"
           LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
           echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
           echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV

.github/workflows/codeql-analysis.yml

@@ -40,11 +40,11 @@ jobs:
 
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
 
       # Initializes the CodeQL tools for scanning.
       - name: Initialize CodeQL
-        uses: github/codeql-action/init@v2
+        uses: github/codeql-action/init@v4
         with:
           languages: ${{ matrix.language }}
           # If you wish to specify custom queries, you can do so here or in a config file.
@@ -57,7 +57,7 @@ jobs:
       # Autobuild attempts to build any compiled languages  (C/C++, C#, or Java).
       # If this step fails, then you should remove it and run the build manually (see below)
       - name: Autobuild
-        uses: github/codeql-action/autobuild@v2
+        uses: github/codeql-action/autobuild@v4
 
       # ℹ️ Command-line programs to run using the OS shell.
       # 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun
@@ -70,4 +70,4 @@ jobs:
       #   ./location_of_script_within_repo/buildscript.sh
 
       - name: Perform CodeQL Analysis
-        uses: github/codeql-action/analyze@v2
+        uses: github/codeql-action/analyze@v4

.github/workflows/mfem-analysis.yml

@@ -39,7 +39,7 @@ jobs:
 
     steps:
       - name: checkout MFEM
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
         with:
           path: mfem
 
@@ -50,7 +50,7 @@ jobs:
 
       - name: Cache Hypre Install
         id: hypre-cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{ env.HYPRE_TOP_DIR }}
           key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-v2.5
@@ -65,7 +65,7 @@ jobs:
 
       - name: Cache Metis Install
         id: metis-cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{ env.METIS_TOP_DIR }}
           key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5

.github/workflows/repo-check.yml

@@ -38,7 +38,7 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
 
       - name: copyright check
         id: copyright
@@ -93,7 +93,7 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
 
       - name: get astyle
         run: |
@@ -110,7 +110,7 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
 
       - name: get doxygen and graphviz
         run: |
@@ -135,7 +135,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
         with:
           fetch-depth: 0
 

.github/workflows/sanitize-build-hypre.yml

@@ -17,11 +17,11 @@ jobs:
     runs-on: ubuntu-latest
     name: 2.19.0
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/config
       - name: Cache
         id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{env.HYPRE_DIR}}
           key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5

.github/workflows/sanitize-build-libcxx.yml

@@ -27,13 +27,13 @@ jobs:
             llvm_use_sanitizer: "Undefined"
     name: ${{matrix.sanitizer}}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/config
         with:
           NO_FLAGS: true
       - name: Cache
         id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{env.LLVM_DIR}}
           key: build-libcxx-${{env.LLVM_VER}}-${{matrix.sanitizer}}

.github/workflows/sanitize-build-lsan.yml

@@ -17,11 +17,11 @@ jobs:
     runs-on: ubuntu-latest
     name: lsan.supp
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/config
       - name: Cache
         id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{env.LSAN_DIR}}
           key: build-lsan-suppression-file

.github/workflows/sanitize-build-metis.yml

@@ -17,11 +17,11 @@ jobs:
     runs-on: ubuntu-latest
     name: 4.0.3
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/config
       - name: Cache
         id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{env.METIS_DIR}}
           key: ${{runner.os}}-build-${{env.METIS_DIR}}-v2.5

.github/workflows/sanitize-tests.yml

@@ -28,7 +28,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/mfem
         with:
           par: ${{inputs.par}}
@@ -40,7 +40,7 @@ jobs:
     env:
       ex: ${{inputs.par && 'ex1p' || 'ex1'}}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -58,7 +58,7 @@ jobs:
     env:
       exclude: ${{inputs.par && '-E "_ser"' || ''}}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -82,7 +82,7 @@ jobs:
     env:
       exclude: ${{inputs.par && '-E "_ser"' || ''}}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -107,7 +107,7 @@ jobs:
       run: ${{inputs.par && '-R "_cpu_np"' || ''}}
       exclude: ${{inputs.par && '"unit_tests|debug"' || '"^unit_tests$|debug"'}}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -131,7 +131,7 @@ jobs:
     env:
       unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -165,7 +165,7 @@ jobs:
       unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
       np: ${{inputs.par && '_np=2' || ''}}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:

.gitignore

@@ -443,6 +443,10 @@ miniapps/diag-smoothers/mg-abs-l1-jacobi
 miniapps/contact/contact
 miniapps/contact/ParaView
 
+miniapps/plasma/pic/electrostatic-*
+!miniapps/plasma/pic/electrostatic-*.cpp
+miniapps/plasma/pic/*.csv
+
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests

.gitlab/jobs/dane.yml

@@ -85,3 +85,8 @@ opt_par_gcc_10_pumi:
   extends: .mfem_job_on_dane
   variables:
     SPEC: "%gcc@10.3.1 +pumi"
+
+opt_par_gcc_10_gslib:
+  extends: .mfem_job_on_dane
+  variables:
+    SPEC: "%gcc@10.3.1 +gslib"

.gitlab/jobs/matrix.yml

@@ -63,3 +63,8 @@ opt_mpi_cuda_hypre_cuda_gcc:
   extends: .mfem_job_on_matrix
   variables:
     SPEC: "%gcc@10.3.1 +mpi +cuda cuda_arch=90 ^hypre+cuda"
+
+opt_mpi_cuda_gcc_gslib:
+  extends: .mfem_job_on_matrix
+  variables:
+    SPEC: "%gcc@10.3.1 +mpi +cuda +gslib cuda_arch=90 ^hypre+cuda"

.gitlab/scripts/baseline

@@ -32,9 +32,9 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
 
 # run
 if [[ "${MACHINE_NAME}" == "dane" ]]; then
-  salloc --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  srun --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
-  salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  srun --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
   echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
   exit 1

CHANGELOG

@@ -11,8 +11,20 @@
 Version 4.9.1 (development)
 ===========================
 
+- Policy for AI-assisted contribution added to CONTRIBUTING.md
+
 Discretization improvements
 ---------------------------
+- Replaced legacy simplex quadrature rules with symmetric positive-weight
+  rules for triangles (orders 0-25) and tetrahedra (orders 0-20). These
+  rules guarantee all-positive weights and interior quadrature points,
+  improving numerical stability. Higher orders fall back to Grundmann-Moller.
+  Triangle rules: Witherden & Vincent, Comput. Math. Appl. 69(10):1232-1241,
+  2015.
+  Tet rules (d=1-13): Witherden & Vincent (ibid).
+  Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
+  2022.
+
 - Improved the gridfunction projection routines. Projections work for Scalar,
   Vector and VectorFE, also NURBS versions. Optionally different types of
   projections can be selected, default behaviour has not changed.

CMakeLists.txt

@@ -652,6 +652,8 @@ foreach(TPL IN LISTS MFEM_TPLS)
   endif()
 endforeach(TPL)
 
+# reverse to remove the first instance of entries in TPL_LIBRARIES
+# so later duplicates are kept (for dependency ordering)
 list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
 list(REVERSE TPL_LIBRARIES)
@@ -1015,5 +1017,7 @@ install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
 # Create 'config.mk' from 'config.mk.in' for the build and install locations and
 # define install rules for 'config.mk' and 'test.mk'
 #-------------------------------------------------------------------------------
-
+if (MFEM_USE_CUDA OR MFEM_USE_HIP)
+  option(MFEM_EXPORT_GPU_CONFIG "Export config.mk for GPU-enabled downstream packages" ON)
+endif()
 mfem_export_mk_files()

CONTRIBUTING.md

@@ -24,6 +24,14 @@ must be made under this license.
 Note also that MFEM has a [Code of Conduct](CODE_OF_CONDUCT.md). By participating
 in the MFEM community, you agree to abide by its rules.
 
+## AI Policy
+- Use of AI code generation in MFEM is allowed but must be disclosed, e.g. by
+  selecting the `AI-assisted` label on the PR.
+- By submitting a PR, the author acknowledges that they have reviewed and
+  understand the changes they are proposing.
+- PR authors are still responsible for correctness, licensing, and attribution
+  of all changes.
+
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull

config/cmake/MFEMConfig.cmake.in

@@ -109,6 +109,10 @@ if (MFEM_USE_RAJA)
    find_dependency(RAJA)
 endif()
 
+if (MFEM_USE_UMPIRE)
+   find_dependency(umpire)
+endif()
+
 if (NOT TARGET mfem)
   include(${CMAKE_CURRENT_LIST_DIR}/MFEMTargets.cmake)
 endif (NOT TARGET mfem)

config/cmake/modules/FindUMPIRE.cmake

@@ -14,12 +14,12 @@
 #   - UMPIRE_LIBRARIES
 #   - UMPIRE_INCLUDE_DIRS
 
-if (NOT umpire_DIR AND UMPIRE_DIR)
-  set(umpire_DIR ${UMPIRE_DIR}/lib/cmake/umpire)
+if (NOT umpire_ROOT AND UMPIRE_DIR)
+  set(umpire_ROOT ${UMPIRE_DIR})
 endif()
 message(STATUS "Looking for UMPIRE ...")
 message(STATUS "   in UMPIRE_DIR = ${UMPIRE_DIR}")
-message(STATUS "      umpire_DIR = ${umpire_DIR}")
+message(STATUS "      umpire_ROOT = ${umpire_ROOT}")
 find_package(umpire CONFIG)
 set(UMPIRE_FOUND ${umpire_FOUND})
 set(UMPIRE_LIBRARIES "umpire")

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -701,7 +701,6 @@ endfunction(mfem_find_library)
 # Extract compile and link options needed by the given target.
 #
 function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
-
   if (NOT TARGET ${Target})
     return()
   endif()
@@ -799,7 +798,12 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
           # message(STATUS "Lib = ${Lib}")
           # Filter-out generator expressions
           if (NOT ("${Lib}" MATCHES "^\\$"))
-            list(APPEND LinkOpts "${Lib}")
+            if(NOT ("${Lib}" STREQUAL "dl"))
+              list(APPEND LinkOpts "${Lib}")
+            else()
+              # for some reason libdl doesn't include the "-l"
+              list(APPEND LinkOpts "-ldl")
+            endif()
           endif()
         else()
           mfem_get_target_options(${Lib} COpts LOpts)
@@ -888,9 +892,18 @@ function(mfem_export_mk_files)
       set(${var} NO)
     endif()
   endforeach()
-  # TODO: Add support for MFEM_USE_CUDA=YES
-  set(MFEM_CXX ${CMAKE_CXX_COMPILER})
-  set(MFEM_HOST_CXX ${MFEM_CXX})
+  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
+    set(MFEM_CXX ${CMAKE_CUDA_COMPILER})
+    if(MFEM_CUDA_COMPILER_IS_NVCC)
+      set(MFEM_HOST_CXX ${CMAKE_CUDA_HOST_COMPILER})
+    else()
+      set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
+    endif()
+  else()
+    # mfem doesn't use enable_language(HIP)
+    set(MFEM_CXX ${CMAKE_CXX_COMPILER})
+    set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
+  endif()
   set(MFEM_CPPFLAGS "")
   get_target_property(cxx_std mfem CXX_STANDARD)
   # For now, we ignore the setting of the CXX_EXTENSIONS property. If this
@@ -900,6 +913,50 @@ function(mfem_export_mk_files)
   string(STRIP
          "${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
          MFEM_CXXFLAGS)
+  if(MFEM_EXPORT_GPU_CONFIG)
+    if (MFEM_USE_CUDA)
+      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} ${CMAKE_CUDA_FLAGS}")
+      if (MFEM_CUDA_COMPILER_IS_NVCC)
+        set(MFEM_CXXFLAGS "-x=cu ${MFEM_CXXFLAGS} -ccbin ${CMAKE_CXX_COMPILER} --forward-unknown-to-host-compiler")
+        # The following intentionally hides CUDA deprecation warnings
+        foreach(ENTRY IN LISTS CUDAToolkit_INCLUDE_DIRS)
+          set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -isystem ${ENTRY}")
+        endforeach()
+        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
+          # architecture flags not part of CMAKE_CUDA_FLAGS
+          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
+            set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -arch=${CMAKE_CUDA_ARCHITECTURES}")
+          else()
+            foreach (ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
+              set(MFEM_CXXFLAGS
+                "${MFEM_CXXFLAGS} -gencode arch=compute_${ENTRY},code=sm_${ENTRY}")
+            endforeach()
+          endif()
+        endif()
+      else()
+        set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xcuda --cuda-path=${CUDAToolkit_LIBRARY_ROOT}")
+        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
+          # architecture flags not part of CMAKE_CUDA_FLAGS
+          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
+            # TODO: not supported
+          else()
+            foreach(ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
+              set(MFEM_CXXFLAGS "-cuda-gpu-arch=sm_${ENTRY} ${MFEM_CXXFLAGS}")
+            endforeach()
+          endif()
+        endif()
+      endif()
+    elseif (MFEM_USE_HIP)
+      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xhip")
+      foreach(ENTRY IN LISTS CMAKE_HIP_ARCHITECTURES)
+        set(MFEM_CXXFLAGS "--offload-arch=${ENTRY} ${MFEM_CXXFLAGS}")
+      endforeach()
+    endif()
+  endif()
   set(MFEM_TPLFLAGS "")
   foreach(dir ${TPL_INCLUDE_DIRS})
     set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
@@ -930,6 +987,9 @@ function(mfem_export_mk_files)
     set(MFEM_SHARED NO)
     set(MFEM_STATIC YES)
   endif()
+  if (MFEM_USE_CUDA)
+    set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} -lcudart")
+  endif()
   set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
   set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
   # For the next 4 variables, these are the values for the build-tree version of
@@ -938,8 +998,15 @@ function(mfem_export_mk_files)
   set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
   set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
   set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
-  # TODO: CUDA/HIP support:
-  set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
+  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
+    if (MFEM_CUDA_COMPILER_IS_NVCC)
+      set(MFEM_XLINKER "-Xlinker=")
+    else()
+      set(MFEM_XLINKER "${CMAKE_CUDA_LINKER_WRAPPER_FLAG}")
+    endif()
+  else()
+    set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
+  endif()
   set(MFEM_MPIEXEC ${MPIEXEC})
   if (NOT MFEM_MPIEXEC)
     set(MFEM_MPIEXEC "mpirun")
@@ -987,16 +1054,21 @@ function(mfem_export_mk_files)
     # handle interfaces (e.g., SCOREC::apf)
     if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*" OR "${lib}" MATCHES "ParMoonolith::.*")
     elseif (TARGET "${lib}")
-      mfem_get_target_options(${lib} CompileOpts LinkOpts)
+      mfem_get_target_options(${lib} CompileOpts2 LinkOpts2)
+      # remove generator expressions
+      string(GENEX_STRIP "${CompileOpts2}" CompileOpts)
+      string(GENEX_STRIP "${LinkOpts2}" LinkOpts)
       # Removing duplicates may lead to issues:
       # list(REMOVE_DUPLICATES CompileOpts)
       # list(REMOVE_DUPLICATES LinkOpts)
-      string(REPLACE ";" " " COpts "${CompileOpts}")
-      string(REPLACE ";" " " LOpts "${LinkOpts}")
-      # message(STATUS "${lib}[COpts]: '${COpts}'")
-      # message(STATUS "${lib}[LOpts]: '${LOpts}'")
-      set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpts}")
-      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpts}")
+      # message(WARNING "${lib}[LinkOpts]: ${LinkOpts}")      
+      # message(WARNING "${lib}[CompileOpts]: ${CompileOpts}")
+      foreach(LOpt IN LISTS LinkOpts)
+        set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpt}")
+      endforeach()
+      foreach(COpt IN LISTS CompileOpts)
+        set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpt}")
+      endforeach()
       # message(FATAL_ERROR "***** interface lib found ... exiting *****")
       # handle static and shared libs
     elseif ("${suffix}" STREQUAL "${CMAKE_SHARED_LIBRARY_SUFFIX}")
@@ -1004,7 +1076,7 @@ function(mfem_export_mk_files)
       get_filename_component(fullLibName ${lib} NAME_WE)
       string(REGEX REPLACE "^lib" "" libname ${fullLibName})
       set(MFEM_EXT_LIBS
-          "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
+        "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
     else()
       set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${lib}")
     endif()
@@ -1013,7 +1085,7 @@ function(mfem_export_mk_files)
   # Create the build-tree version of 'config.mk'
   configure_file(
     "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config.mk")
+    "${PROJECT_BINARY_DIR}/config/config.mk" @ONLY)
   # Copy 'test.mk' from the source-tree to the build-tree
   configure_file(
     "${PROJECT_SOURCE_DIR}/config/test.mk"
@@ -1031,7 +1103,7 @@ function(mfem_export_mk_files)
   # Create the install-tree version of 'config.mk'
   configure_file(
     "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config-install.mk")
+    "${PROJECT_BINARY_DIR}/config/config-install.mk" @ONLY)
 
   # Install rules for 'config.mk' and 'test.mk'
   install(FILES ${PROJECT_SOURCE_DIR}/config/test.mk

examples/CMakeLists.txt

@@ -215,7 +215,7 @@ if (MFEM_ENABLE_TESTING)
     add_test(NAME ex1p_ceed_np=${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:ex1p> "-no-vis" "-d ceed-cpu" "-pa" "-a"
+      $<TARGET_FILE:ex1p> "-no-vis" "-d" "ceed-cpu" "-pa" "-a"
       ${MPIEXEC_POSTFLAGS})
   endif()
 endif()

examples/amgx/makefile

@@ -64,7 +64,7 @@ PARALLEL_NAME := Parallel AMGX example
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not build)
 
-clean: clean-build
+clean: clean-build clean-exec
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)

examples/caliper/makefile

@@ -64,12 +64,12 @@ ex1p-test-par: ex1p
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-clean: clean-build clean-exec $(SUBDIRS_CLEAN)
+clean: clean-build clean-exec
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
-	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
+	@rm -f refined.mesh mesh.*
+	@rm -f sol.*

examples/ex12p.cpp

@@ -5,9 +5,9 @@
 // Sample runs:
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh
-//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 462 -n 10 -o 2 -elast
+//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 464 -n 10 -o 2 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-hex.mesh -s 3878
-//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 81
+//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 82
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh -s 3877 -o 2 -sys
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh -s 4544 -n 6 -o 3 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh

examples/ex22p.cpp

@@ -302,15 +302,21 @@ int main(int argc, char *argv[])
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
       sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
       sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
       sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                  << "window_title 'Exact: Real Part'" << flush;
+      // Make sure all ranks have sent their real solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i.precision(8);
       sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                  << "window_title 'Exact: Imaginary Part'" << flush;
+      // Make sure all ranks have sent their imaginary solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
    }
 
    // 11. Set up the parallel sesquilinear form a(.,.) on the finite element
@@ -534,15 +540,21 @@ int main(int argc, char *argv[])
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
       sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
       sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
       sol_sock_r << "solution\n" << *pmesh << u.real()
                  << "window_title 'Solution: Real Part'" << flush;
+      // Make sure all ranks have sent their real solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i.precision(8);
       sol_sock_i << "solution\n" << *pmesh << u.imag()
                  << "window_title 'Solution: Imaginary Part'" << flush;
+      // Make sure all ranks have sent their imaginary solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
    }
    if (visualization && exact_sol)
    {
@@ -551,15 +563,21 @@ int main(int argc, char *argv[])
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
       sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
       sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
       sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                  << "window_title 'Error: Real Part'" << flush;
+      // Make sure all ranks have sent their real solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i.precision(8);
       sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                  << "window_title 'Error: Imaginary Part'" << flush;
+      // Make sure all ranks have sent their imaginary solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
    }
    if (visualization)
    {

examples/ex37.cpp

@@ -5,8 +5,8 @@
 // Sample runs:
 //     ex37 -alpha 10
 //     ex37 -alpha 10 -pv
-//     ex37 -lambda 0.1 -mu 0.1
-//     ex37 -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
+//     ex37 -lambda 0.1 -mu 0.1 -growth 1
+//     ex37 -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
 //     ex37 -r 6 -o 1 -alpha 25.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -55,53 +55,6 @@
 using namespace std;
 using namespace mfem;
 
-/**
- * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
- *        ∫_Ω ρ dx = θ vol(Ω) as follows:
- *
- *        1. Compute the root of the R → R function
- *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
- *        2. Set ψ ← ψ + c.
- *
- * @param psi a GridFunction to be updated
- * @param target_volume θ vol(Ω)
- * @param tol Newton iteration tolerance
- * @param max_its Newton maximum iteration number
- * @return real_t Final volume, ∫_Ω sigmoid(ψ)
- */
-real_t proj(GridFunction &psi, real_t target_volume, real_t tol=1e-12,
-            int max_its=10)
-{
-   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
-   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
-
-   LinearForm int_sigmoid_psi(psi.FESpace());
-   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
-   LinearForm int_der_sigmoid_psi(psi.FESpace());
-   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
-                                              der_sigmoid_psi));
-   bool done = false;
-   for (int k=0; k<max_its; k++) // Newton iteration
-   {
-      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
-      const real_t f = int_sigmoid_psi.Sum() - target_volume;
-
-      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
-      const real_t df = int_der_sigmoid_psi.Sum();
-
-      const real_t dc = -f/df;
-      psi += dc;
-      if (abs(dc) < tol) { done = true; break; }
-   }
-   if (!done)
-   {
-      mfem_warning("Projection reached maximum iteration without converging. "
-                   "Result may not be accurate.");
-   }
-   int_sigmoid_psi.Assemble();
-   return int_sigmoid_psi.Sum();
-}
-
 /*
  * ---------------------------------------------------------------
  *                      ALGORITHM PREAMBLE
@@ -180,10 +133,11 @@ int main(int argc, char *argv[])
    int ref_levels = 5;
    int order = 2;
    real_t alpha = 1.0;
+   real_t growth = 2;
    real_t epsilon = 0.01;
    real_t vol_fraction = 0.5;
    int max_it = 1e3;
-   real_t itol = 1e-1;
+   real_t itol = 1e-2;
    real_t ntol = 1e-4;
    real_t rho_min = 1e-6;
    real_t lambda = 1.0;
@@ -198,6 +152,8 @@ int main(int argc, char *argv[])
                   "Order (degree) of the finite elements.");
    args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                   "Step length for gradient descent.");
+   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
+                  "Growth rate of step length for gradient descent.");
    args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                   "Length scale for ρ.");
    args.AddOption(&max_it, "-mi", "--max-it",
@@ -332,6 +288,7 @@ int main(int argc, char *argv[])
    }
    FilterSolver->SetEssentialBoundary(ess_bdr_filter);
    FilterSolver->SetupFEM();
+   FilterSolver->AssembleDiffusionBilinear();
 
    BilinearForm mass(&control_fes);
    mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -385,7 +342,7 @@ int main(int argc, char *argv[])
    // 11. Iterate:
    for (int k = 1; k <= max_it; k++)
    {
-      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }
+      if (k > 1) { alpha = std::pow((real_t) k,growth); }
 
       mfem::out << "\nStep = " << k << std::endl;
 
@@ -422,7 +379,9 @@ int main(int argc, char *argv[])
 
       // Step 5 - Update design variable ψ ← proj(ψ - αG)
       psi.Add(-alpha, grad);
-      const real_t material_volume = proj(psi, target_volume);
+      GridFunction alpha_grad(grad);
+      alpha_grad *= alpha;
+      const real_t material_volume = proj(psi, alpha_grad, target_volume);
 
       // Compute ||ρ - ρ_old|| in control fes.
       real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);

examples/ex37.hpp

@@ -137,7 +137,7 @@ public:
         exponent(exponent_), rho_min(rho_min_)
    {
       MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
-      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be > 1");
+      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be < 1");
       MFEM_ASSERT(u, "displacement field is not set");
       MFEM_ASSERT(rho_filter, "density field is not set");
    }
@@ -231,9 +231,12 @@ private:
    FiniteElementCollection * fec = nullptr;
    FiniteElementSpace * fes = nullptr;
    Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
    Array<int> neumann_bdr;
    GridFunction * u = nullptr;
    LinearForm * b = nullptr;
+   BilinearForm * a = nullptr;
+   OperatorPtr A;
    bool parallel;
 #ifdef MFEM_USE_MPI
    ParMesh * pmesh = nullptr;
@@ -267,6 +270,8 @@ public:
    void ResetFEM();
    void SetupFEM();
 
+   void UpdateEssentialTDofs();
+   void AssembleDiffusionBilinear(bool update_ess_tdofs=true);
    void Solve();
    GridFunction * GetFEMSolution();
    LinearForm * GetLinearForm() {return b;}
@@ -371,6 +376,130 @@ public:
 
 };
 
+/**
+ * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
+ *        ∫_Ω ρ dx = θ vol(Ω) as follows:
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
+ *           using the Illinois method
+ *        2. Set ψ ← ψ + c.
+ *
+ * @param psi a GridFunction to be updated
+ * @param alpha_grad alpha multiplied by gradient
+ * @param target_volume θ vol(Ω)
+ * @param tol Illinois iteration tolerance
+ * @param max_its Illinois maximum iteration number
+ * @return real_t Final volume (∫_Ω sigmoid(ψ) dx)
+ */
+real_t proj(GridFunction &psi, GridFunction &alpha_grad, real_t target_volume,
+            real_t tol = 1e-12, int max_its = 100)
+{
+#ifdef MFEM_USE_MPI
+   FiniteElementSpace *fes = psi.FESpace();
+   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+#endif
+   ConstantCoefficient zero_cf(0.0);
+   real_t a = -alpha_grad.ComputeMaxError(zero_cf);
+   real_t b = -a;
+   real_t y = 0.0;
+
+   MappedGridFunctionCoefficient sigmoid_psi(
+   &psi, [&y](const real_t x) { return sigmoid(x + y); });
+   std::unique_ptr<LinearForm> int_sigmoid_psi;
+#ifdef MFEM_USE_MPI
+   ParGridFunction *par_psi = dynamic_cast<ParGridFunction *>(&psi);
+   if (par_psi)
+   {
+      int_sigmoid_psi.reset(new ParLinearForm(par_psi->ParFESpace()));
+   }
+   else
+   {
+      int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
+   }
+#else
+   int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
+#endif
+   int_sigmoid_psi->AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+
+   y = a;
+   int_sigmoid_psi->Assemble();
+   real_t f_a = int_sigmoid_psi->Sum();  // f_a := f(a) + θ vol(Ω)
+
+   y = b;
+   int_sigmoid_psi->Assemble();
+   real_t f_b = int_sigmoid_psi->Sum();  // f_b := f(b) + θ vol(Ω)
+#ifdef MFEM_USE_MPI
+   if (pfes)
+   {
+      MPI_Allreduce(MPI_IN_PLACE, &f_a, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+      MPI_Allreduce(MPI_IN_PLACE, &f_b, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+   }
+#endif
+   f_a -= target_volume; // f_a := f(a)
+   f_b -= target_volume; // f_b := f(b)
+   real_t c = 0.0;
+   real_t f_c = 0.0;
+   int side = 0;
+
+   bool done = false;
+   for (int k=0; k < max_its; k++)
+   {
+      c = (f_a * b - f_b * a) / (f_a - f_b);
+
+      if (abs(b - a) < tol * abs(b + a)) { done = true; break; }
+
+      y = c;
+      int_sigmoid_psi->Assemble();
+      f_c = int_sigmoid_psi->Sum(); // f_c := f(c) + θ vol(Ω)
+#ifdef MFEM_USE_MPI
+      if (pfes)
+      {
+         MPI_Allreduce(MPI_IN_PLACE, &f_c, 1, MPITypeMap<real_t>::mpi_type,
+                       MPI_SUM, MPI_COMM_WORLD);
+      }
+#endif
+      f_c -= target_volume; // f_c := f(c)
+
+      if (f_c * f_b > 0)
+      {
+         b = c;
+         f_b = f_c;
+         if (side == -1) { f_a /= 2.0; }
+         side = -1;
+      }
+      else if (f_c * f_a > 0)
+      {
+         a = c;
+         f_a = f_c;
+         if (side == 1) { f_b /= 2.0; }
+         side = 1;
+      }
+      else
+      {
+         done = true; break;
+      }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. "
+                   "Result may not be accurate.");
+   }
+   y = 0.0;
+   psi += c;
+   int_sigmoid_psi->Assemble();
+   real_t material_volume = int_sigmoid_psi->Sum();
+#ifdef MFEM_USE_MPI
+   if (pfes)
+   {
+      MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
+                    MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
+   }
+#endif
+   return material_volume;
+}
 
 // Poisson solver
 
@@ -422,12 +551,8 @@ void DiffusionSolver::SetupFEM()
    }
 }
 
-void DiffusionSolver::Solve()
+void DiffusionSolver::UpdateEssentialTDofs()
 {
-   OperatorPtr A;
-   Vector B, X;
-   Array<int> ess_tdof_list;
-
 #ifdef MFEM_USE_MPI
    if (parallel)
    {
@@ -440,7 +565,39 @@ void DiffusionSolver::Solve()
 #else
    fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
 #endif
-   *u=0.0;
+}
+
+void DiffusionSolver::AssembleDiffusionBilinear(bool update_ess_tdofs)
+{
+   if (update_ess_tdofs)
+   {
+      UpdateEssentialTDofs();
+   }
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      a = new ParBilinearForm(pfes);
+   }
+   else
+   {
+      a = new BilinearForm(fes);
+   }
+#else
+   a = new BilinearForm(fes);
+#endif
+   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
+   if (masscf)
+   {
+      a->AddDomainIntegrator(new MassIntegrator(*masscf));
+   }
+   a->Assemble();
+   a->FormSystemMatrix(ess_tdof_list, A);
+}
+
+void DiffusionSolver::Solve()
+{
+   Vector B, X;
+
    if (b)
    {
       delete b;
@@ -475,31 +632,33 @@ void DiffusionSolver::Solve()
 
    b->Assemble();
 
-   BilinearForm * a = nullptr;
-
-#ifdef MFEM_USE_MPI
-   if (parallel)
-   {
-      a = new ParBilinearForm(pfes);
-   }
-   else
-   {
-      a = new BilinearForm(fes);
-   }
-#else
-   a = new BilinearForm(fes);
-#endif
-   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
-   if (masscf)
-   {
-      a->AddDomainIntegrator(new MassIntegrator(*masscf));
-   }
-   a->Assemble();
+   *u=0.0;
    if (essbdr_cf)
    {
       u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
    }
-   a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B);
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      X.SetSize(pfes->TrueVSize());
+      B.SetSize(pfes->TrueVSize());
+      dynamic_cast<ParGridFunction*>(u)->ParallelAssemble(X);
+      dynamic_cast<ParLinearForm*>(b)->ParallelAssemble(B);
+      dynamic_cast<ParBilinearForm*>(a)->ParallelEliminateTDofsInRHS(
+         ess_tdof_list, X, B);
+   }
+   else
+   {
+      X.NewDataAndSize(u->GetData(), u->Size());
+      B.NewDataAndSize(b->GetData(), b->Size());
+      a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+   }
+#else
+   X.NewDataAndSize(u->GetData(), u->Size());
+   B.NewDataAndSize(b->GetData(), b->Size());
+   a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+#endif
 
    CGSolver * cg = nullptr;
    Solver * M = nullptr;
@@ -528,7 +687,6 @@ void DiffusionSolver::Solve()
    delete M;
    delete cg;
    a->RecoverFEMSolution(X, *b, *u);
-   delete a;
 }
 
 GridFunction * DiffusionSolver::GetFEMSolution()
@@ -560,6 +718,8 @@ DiffusionSolver::~DiffusionSolver()
 #endif
    delete fec; fec = nullptr;
    delete b;
+   A.Clear();
+   delete a;
 }
 
 

examples/ex37p.cpp

@@ -4,8 +4,8 @@
 //
 // Sample runs:
 //    mpirun -np 4 ex37p -alpha 10 -pv
-//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1
-//    mpirun -np 4 ex37p -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
+//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1 -growth 1
+//    mpirun -np 4 ex37p -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
 //    mpirun -np 4 ex37p -r 6 -o 2 -alpha 10.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -54,61 +54,6 @@
 using namespace std;
 using namespace mfem;
 
-/**
- * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
- *        ∫_Ω ρ dx = θ vol(Ω) as follows:
- *
- *        1. Compute the root of the R → R function
- *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
- *        2. Set ψ ← ψ + c.
- *
- * @param psi a GridFunction to be updated
- * @param target_volume θ vol(Ω)
- * @param tol Newton iteration tolerance
- * @param max_its Newton maximum iteration number
- * @return real_t Final volume, ∫_Ω sigmoid(ψ)
- */
-real_t proj(ParGridFunction &psi, real_t target_volume, real_t tol=1e-12,
-            int max_its=10)
-{
-   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
-   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
-
-   ParLinearForm int_sigmoid_psi(psi.ParFESpace());
-   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
-   ParLinearForm int_der_sigmoid_psi(psi.ParFESpace());
-   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
-                                              der_sigmoid_psi));
-   bool done = false;
-   for (int k=0; k<max_its; k++) // Newton iteration
-   {
-      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
-      real_t f = int_sigmoid_psi.Sum();
-      MPI_Allreduce(MPI_IN_PLACE, &f, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-      f -= target_volume;
-
-      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
-      real_t df = int_der_sigmoid_psi.Sum();
-      MPI_Allreduce(MPI_IN_PLACE, &df, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-
-      const real_t dc = -f/df;
-      psi += dc;
-      if (abs(dc) < tol) { done = true; break; }
-   }
-   if (!done)
-   {
-      mfem_warning("Projection reached maximum iteration without converging. "
-                   "Result may not be accurate.");
-   }
-   int_sigmoid_psi.Assemble();
-   real_t material_volume = int_sigmoid_psi.Sum();
-   MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
-                 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
-   return material_volume;
-}
-
 /*
  * ---------------------------------------------------------------
  *                      ALGORITHM PREAMBLE
@@ -193,10 +138,11 @@ int main(int argc, char *argv[])
    int ref_levels = 5;
    int order = 2;
    real_t alpha = 1.0;
+   real_t growth = 2;
    real_t epsilon = 0.01;
    real_t vol_fraction = 0.5;
    int max_it = 1e3;
-   real_t itol = 1e-1;
+   real_t itol = 1e-2;
    real_t ntol = 1e-4;
    real_t rho_min = 1e-6;
    real_t lambda = 1.0;
@@ -211,6 +157,8 @@ int main(int argc, char *argv[])
                   "Order (degree) of the finite elements.");
    args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                   "Step length for gradient descent.");
+   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
+                  "Growth rate of step length for gradient descent.");
    args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                   "Length scale for ρ.");
    args.AddOption(&max_it, "-mi", "--max-it",
@@ -359,6 +307,7 @@ int main(int argc, char *argv[])
    }
    FilterSolver->SetEssentialBoundary(ess_bdr_filter);
    FilterSolver->SetupFEM();
+   FilterSolver->AssembleDiffusionBilinear();
 
    ParBilinearForm mass(&control_fes);
    mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -412,7 +361,7 @@ int main(int argc, char *argv[])
    // 11. Iterate:
    for (int k = 1; k <= max_it; k++)
    {
-      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }
+      if (k > 1) { alpha = std::pow((real_t) k,growth); }
 
       if (myid == 0)
       {
@@ -452,7 +401,9 @@ int main(int argc, char *argv[])
 
       // Step 5 - Update design variable ψ ← proj(ψ - αG)
       psi.Add(-alpha, grad);
-      const real_t material_volume = proj(psi, target_volume);
+      ParGridFunction alpha_grad(grad);
+      alpha_grad *= alpha;
+      const real_t material_volume = proj(psi, alpha_grad, target_volume);
 
       // Compute ||ρ - ρ_old|| in control fes.
       real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);

examples/ginkgo/makefile

@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh sol.gf
+	@rm -f refined.mesh sol.gf mesh.* sol.*

examples/makefile

@@ -71,6 +71,7 @@ endif
 
 SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
 SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
+SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(SUBDIRS))
 SUBDIRS_CLEAN = $(addsuffix /clean,$(SUBDIRS))
 SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
@@ -87,8 +88,9 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
 all: $(EXAMPLES) $(SUBDIRS_ALL)
 
-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
+   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
 	$(MAKE) -C $(@D) $(@F)
 $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
@@ -107,6 +109,7 @@ endif
 MFEM_TESTS = EXAMPLES
 include $(MFEM_TEST_MK)
 test: $(SUBDIRS_TEST)
+test-noclean: $(SUBDIRS_TEST_NOCLEAN)
 test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Parallel vs. serial runs
@@ -157,6 +160,8 @@ ex37-test-seq: ex37
 	@$(call mfem-test,$<,, Serial example,-mi 3)
 ex37p-test-par: ex37p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-mi 3)
+ex39-test-seq: ex39
+	@$(call mfem-test,$<,, Serial example,-m ../data/compass.mesh)
 ex41-test-seq: ex41
 	@$(call mfem-test,$<,, Serial example,-tf 1.0)
 ex41p-test-par: ex41p

fem/bilinearform.cpp

@@ -729,7 +729,8 @@ void BilinearForm::Assemble(int skip_zeros)
          tr = mesh -> GetBdrFaceTransformations (i);
          if (tr != NULL)
          {
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
+            mfem::DofTransformation doftrans;
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
             fe1 = fes -> GetFE (tr -> Elem1No);
             // The fe2 object is really a dummy and not used on the boundaries,
             // but we can't dereference a NULL pointer, and we don't want to
@@ -743,6 +744,7 @@ void BilinearForm::Assemble(int skip_zeros)
 
                boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
                                                               elemmat);
+               doftrans.TransformDual(elemmat);
                mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
             }
          }
@@ -1723,6 +1725,7 @@ void MixedBilinearForm::Assemble(int skip_zeros)
          }
       }
 
+      DofTransformation dom_dof_trans, ran_dof_trans;
       for (int i = 0; i < trial_fes -> GetNBE(); i++)
       {
          const int bdr_attr = mesh->GetBdrAttribute(i);
@@ -1731,8 +1734,8 @@ void MixedBilinearForm::Assemble(int skip_zeros)
          ftr = mesh -> GetBdrFaceTransformations (i);
          if (ftr != NULL)
          {
-            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
-            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
+            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs, dom_dof_trans);
+            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs, ran_dof_trans);
             trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
             test_fe1 = test_fes->GetFE(ftr->Elem1No);
             // The test_fe2 object is really a dummy and not used on the
@@ -1748,6 +1751,7 @@ void MixedBilinearForm::Assemble(int skip_zeros)
                boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
                                                            *test_fe2,
                                                            *ftr, elemmat);
+               TransformDual(ran_dof_trans, dom_dof_trans, elemmat);
                mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
             }
          }

fem/bilininteg.hpp

@@ -2710,7 +2710,7 @@ public:
 
 
 /** Integrator for $(-Q u, \nabla v)$ for Nedelec ($u$) and $H^1$ ($v$) elements.
-    This is equivalent to a weak divergence of the $H(curl$ basis functions. */
+    This is equivalent to a weak divergence of the $H(curl)$ basis functions. */
 class VectorFEWeakDivergenceIntegrator: public BilinearFormIntegrator
 {
 protected:

fem/bounds.cpp

@@ -41,9 +41,14 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
    tol = tol_i;
    lbound.SetSize(ncp, nb);
    ubound.SetSize(ncp, nb);
+   lbound_t.SetSize(nb, ncp);
+   ubound_t.SetSize(nb, ncp);
    nodes.SetSize(nb);
    weights.SetSize(nb);
    control_points.SetSize(ncp);
+   xhat.SetSize(nb);
+   what.SetSize(nb);
+   cphat.SetSize(ncp);
 
    auto scalenodes = [](const Vector &in, const real_t a, const real_t b) -> Vector
    {
@@ -90,6 +95,10 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
       MFEM_ABORT("Unsupported interval points. Use [0,1].\n");
    }
    control_points = scalenodes(control_points, 0.0, 1.0); // rescale to [0,1]
+   for (int i = 0; i < ncp; i++)
+   {
+      cphat(i) = 2.0*control_points(i) - 1.0;
+   }
 
    Poly_1D::Basis &basis1d(poly1d.GetBasis(nb-1, b_type));
 
@@ -145,6 +154,8 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
                lbound(j,i) = std::max(lbound(j,i),0_r);
             }
          }
+         lbound_t(i,j) = lbound(j,i);
+         ubound_t(i,j) = ubound(j,i);
       }
    }
 
@@ -176,6 +187,11 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
          nodes(i) = irule.IntPoint(i).x;
       }
    }
+   for (int i = 0; i < nb; i++)
+   {
+      xhat(i) = 2.0*nodes(i) - 1.0;
+      what(i) = 2.0*weights(i);
+   }
 
    if (b_type == 2)
    {
@@ -755,4 +771,4 @@ void PLBound::Print(std::ostream &outp) const
    ubound.Print(outp);
 }
 
-}
+}

fem/bounds.hpp

@@ -13,6 +13,7 @@
 #define MFEM_BOUNDS
 
 #include "../config/config.hpp"
+#include "../general/forall.hpp"
 #include "fespace.hpp"
 
 namespace mfem
@@ -60,7 +61,9 @@ private:
    bool proj = true; // Use linear projection to compute bounds.
    real_t tol = 0.0; // offset bounds to avoid round-off errors
    Vector nodes, weights, control_points;
+   Vector xhat, what, cphat;
    DenseMatrix lbound, ubound; // ncp x nb matrices with bounds of all bases
+   DenseMatrix lbound_t, ubound_t; // nb x ncp transposes for device kernel
    // Some auxillary storage for computing the bounds with Bernstein
    DenseMatrix basisMatNodes; // Bernstein bases at equispaced nodes
    DenseMatrix basisMatInt;   // Bernstein bases at GLL nodes
@@ -113,7 +116,10 @@ public:
     *  @details This projection increases the computational cost but results in
     *  tighter bounds.
     */
-   void SetProjectionFlagForBounding(bool proj_) { proj = proj_; }
+   void SetProjectionFlagForBounding(bool proj_)
+   {
+      proj = proj_;
+   }
 
    /** @brief Compute piecewise linear bounds for the lexicographically-ordered
     *  nodal coefficients in @a coeff in 1D/2D/3D.
@@ -137,9 +143,23 @@ public:
    /// Get number of control points used to compute the bounds.
    int GetNControlPoints() const { return ncp; }
 
+   /// Get the underlying 1D basis type.
+   int GetBasisType() const { return b_type; }
+
    /// Get 1D control point locations (lexicographic order) in [0,1].
    const Vector &GetControlPoints() const { return control_points; }
 
+   /** @brief Compute element-wise bounds from a lexicographic E-vector.
+    *
+    *  @details The expected layout of @a e_vec is `ND x VDIM x NE`, where
+    *  `ND = nb^rdim`, `VDIM = fes_vdim`, and `NE` is the number of elements.
+    *  The output layout matches GridFunction::GetElementBounds:
+    *  `NE x active_vdim`, with the element index varying fastest.
+    */
+   void GetElementBoundsKernel(const int rdim, const int fes_vdim,
+                               const Vector &e_vec, Vector &lower,
+                               Vector &upper, const int vdim = 0) const;
+
    /** @brief Get lower and upper bounding matrix (ncp^dim x nb^dim)
     *
     *  @details The matrices can be used to compute the bounds at control points
@@ -183,6 +203,600 @@ private:
               const int cp_type_i, const real_t tol_i);
 };
 
+namespace internal
+{
+
+struct PLBoundDeviceData
+{
+   int nb;
+   int ncp;
+   const real_t *xhat;
+   const real_t *what;
+   const real_t *cphat;
+   const real_t *lbound;
+   const real_t *ubound;
+};
+
+template<int T_NB = 0, bool T_PROJ = true>
+inline void GetElementBoundsKernel1D(const PLBoundDeviceData &data,
+                                     const int fes_vdim,
+                                     const int ne,
+                                     const Vector &e_vec,
+                                     Vector &lower,
+                                     Vector &upper,
+                                     const int comp0,
+                                     const int ncomp)
+{
+   constexpr int GENERIC_MAX_ND = 32;
+   constexpr int MAX_ND = T_NB ? T_NB : GENERIC_MAX_ND;
+   constexpr int BLOCK_X = 2*MAX_ND;
+
+   const int nd = T_NB ? T_NB : data.nb;
+   MFEM_VERIFY(nd <= MAX_ND,
+               "Device element bounds kernel supports up to 32 "
+               "1D degrees of freedom.");
+
+   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
+   auto L = Reshape(lower.Write(), ne, ncomp);
+   auto U = Reshape(upper.Write(), ne, ncomp);
+
+   mfem::forall_2D<BLOCK_X>(ne*ncomp, BLOCK_X, 1,
+                            [=] MFEM_HOST_DEVICE (int ec)
+   {
+      const int e = ec % ne;
+      const int c = ec / ne;
+      const int vc = comp0 + c;
+      const real_t *coeff = &E(0, vc, e);
+      const int tid = MFEM_THREAD_ID(x);
+
+      MFEM_SHARED real_t sproj[MAX_ND];
+      MFEM_SHARED real_t ssum0[MAX_ND];
+      MFEM_SHARED real_t ssum1[MAX_ND];
+      MFEM_SHARED real_t smin[BLOCK_X];
+      MFEM_SHARED real_t smax[BLOCK_X];
+      MFEM_SHARED real_t sa0;
+      MFEM_SHARED real_t sa1;
+
+      MFEM_FOREACH_THREAD(i, x, nd)
+      {
+         if constexpr (T_PROJ)
+         {
+            const real_t x = data.xhat[i];
+            const real_t w = data.what[i];
+            ssum0[i] = 0.5*coeff[i]*w;
+            ssum1[i] = 1.5*coeff[i]*w*x;
+         }
+         else
+         {
+            ssum0[i] = 0.0;
+            ssum1[i] = 0.0;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(ii, x, 1)
+      {
+         sa0 = 0.0;
+         sa1 = 0.0;
+         for (int i = 0; i < nd; i++)
+         {
+            sa0 += ssum0[i];
+            sa1 += ssum1[i];
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(i, x, nd)
+      {
+         if constexpr (T_PROJ)
+         {
+            const real_t x = data.xhat[i];
+            sproj[i] = coeff[i] - sa0 - sa1*x;
+         }
+         else
+         {
+            sproj[i] = coeff[i];
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      real_t lower_local = HUGE_VAL;
+      real_t upper_local = -HUGE_VAL;
+      MFEM_FOREACH_THREAD(j, x, data.ncp)
+      {
+         real_t lo = 0.0;
+         real_t hi = 0.0;
+         if constexpr (T_PROJ)
+         {
+            const real_t xcp = data.cphat[j];
+            lo = sa0 + sa1*xcp;
+            hi = lo;
+         }
+
+         for (int i = 0; i < nd; i++)
+         {
+            const real_t val = sproj[i];
+            const real_t lv = data.lbound[j + i*data.ncp]*val;
+            const real_t uv = data.ubound[j + i*data.ncp]*val;
+            lo += lv < uv ? lv : uv;
+            hi += lv > uv ? lv : uv;
+         }
+         lower_local = lower_local < lo ? lower_local : lo;
+         upper_local = upper_local > hi ? upper_local : hi;
+      }
+
+      smin[tid] = lower_local;
+      smax[tid] = upper_local;
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(ii, x, 1)
+      {
+         real_t lower_ec = smin[0];
+         real_t upper_ec = smax[0];
+         const int nthreads = MFEM_THREAD_SIZE(x);
+         const int nactive = data.ncp < nthreads ? data.ncp : nthreads;
+         for (int t = 1; t < nactive; t++)
+         {
+            lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
+            upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
+         }
+         L(e, c) = lower_ec;
+         U(e, c) = upper_ec;
+      }
+   });
+}
+
+template<int T_NB = 0, int T_NCP = 0, bool T_PROJ = true>
+inline void GetElementBoundsKernel2D(const PLBoundDeviceData &data,
+                                     const int fes_vdim,
+                                     const int ne,
+                                     const Vector &e_vec,
+                                     Vector &lower,
+                                     Vector &upper,
+                                     const int comp0,
+                                     const int ncomp)
+{
+   constexpr int DEFAULT_MAX_NB = 8;
+   constexpr int DEFAULT_MAX_CP = 3*DEFAULT_MAX_NB;
+   constexpr int MAX_NB = T_NB ? T_NB : DEFAULT_MAX_NB;
+   constexpr int MAX_CP = T_NCP ? T_NCP : DEFAULT_MAX_CP;
+   constexpr int MAX_THREADS = MAX_CP*MAX_CP;
+
+   const int nb = data.nb;
+   const int ncp = data.ncp;
+   const int nd = nb*nb;
+   MFEM_VERIFY(nb <= MAX_NB,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "1D degree bound.");
+   MFEM_VERIFY(ncp <= MAX_CP,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "control-point bound.");
+   MFEM_VERIFY(ncp*ncp <= MAX_THREADS,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "thread-block bound.");
+
+   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
+   auto L = Reshape(lower.Write(), ne, ncomp);
+   auto U = Reshape(upper.Write(), ne, ncomp);
+
+   mfem::forall_2D<MAX_THREADS>(ne*ncomp, ncp, ncp,
+                                [=] MFEM_HOST_DEVICE (int ec)
+   {
+      const int e = ec % ne;
+      const int c = ec / ne;
+      const int vc = comp0 + c;
+      const real_t *coeff = &E(0, vc, e);
+      const int tx = MFEM_THREAD_ID(x);
+      const int ty = MFEM_THREAD_ID(y);
+
+      MFEM_SHARED real_t sproj[MAX_NB*MAX_NB];
+      MFEM_SHARED real_t srow_min[MAX_NB*MAX_CP];
+      MFEM_SHARED real_t srow_max[MAX_NB*MAX_CP];
+      MFEM_SHARED real_t srow_a0[MAX_NB];
+      MFEM_SHARED real_t srow_a1[MAX_NB];
+      MFEM_SHARED real_t sa0[MAX_CP];
+      MFEM_SHARED real_t sa1[MAX_CP];
+      MFEM_SHARED real_t smin[MAX_THREADS];
+      MFEM_SHARED real_t smax[MAX_THREADS];
+
+      // Stage 1a: for each nodal row, form the per-node contributions to the
+      // row-wise linear fit used by the first 1D bounding solve.
+      MFEM_FOREACH_THREAD(jrow, y, nb)
+      {
+         const real_t *row_coeff = coeff + jrow*nb;
+         const int row_ncp_off = jrow*MAX_CP;
+         MFEM_FOREACH_THREAD(i, x, nb)
+         {
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[i];
+               const real_t w = data.what[i];
+               srow_min[row_ncp_off + i] = 0.5*row_coeff[i]*w;
+               srow_max[row_ncp_off + i] = 1.5*row_coeff[i]*w*x;
+            }
+            else
+            {
+               srow_min[row_ncp_off + i] = 0.0;
+               srow_max[row_ncp_off + i] = 0.0;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 1b: reduce the row-wise projection coefficients a0/a1.
+      if constexpr (T_PROJ)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_ncp_off = jrow*MAX_CP;
+            real_t a0 = 0.0;
+            real_t a1 = 0.0;
+            MFEM_FOREACH_THREAD(ii, x, 1)
+            {
+               for (int i = 0; i < nb; i++)
+               {
+                  a0 += srow_min[row_ncp_off + i];
+                  a1 += srow_max[row_ncp_off + i];
+               }
+               srow_a0[jrow] = a0;
+               srow_a1[jrow] = a1;
+            }
+         }
+         MFEM_SYNC_THREAD;
+      }
+
+      // Stage 1c: subtract the row-wise linear fit once and cache the
+      // projected row coefficients for reuse across all x-control points.
+      MFEM_FOREACH_THREAD(jrow, y, nb)
+      {
+         const real_t *row_coeff = coeff + jrow*nb;
+         MFEM_FOREACH_THREAD(i, x, nb)
+         {
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[i];
+               sproj[jrow*MAX_NB + i] = row_coeff[i]
+                                        - srow_a0[jrow] - srow_a1[jrow]*x;
+            }
+            else
+            {
+               sproj[jrow*MAX_NB + i] = row_coeff[i];
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 1d: solve the first 1D bounding problem along each nodal row and
+      // store bounds at every x-direction control point.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_cp_off = jrow*ncp;
+            real_t lo = 0.0;
+            real_t hi = 0.0;
+            if constexpr (T_PROJ)
+            {
+               const real_t xcp = data.cphat[icp];
+               lo = srow_a0[jrow] + srow_a1[jrow]*xcp;
+               hi = lo;
+            }
+            for (int i = 0; i < nb; i++)
+            {
+               const real_t val = sproj[jrow*MAX_NB + i];
+               const real_t lv = data.lbound[icp + i*data.ncp]*val;
+               const real_t uv = data.ubound[icp + i*data.ncp]*val;
+               lo += lv < uv ? lv : uv;
+               hi += lv > uv ? lv : uv;
+            }
+            srow_min[row_cp_off + icp] = lo;
+            srow_max[row_cp_off + icp] = hi;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2a: from the row bounds, form the per-row contributions to the
+      // second 1D projection solve in the y-direction.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_cp_off = jrow*ncp;
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[jrow];
+               const real_t w = data.what[jrow];
+               const real_t t = 0.5*(srow_min[row_cp_off + icp] +
+                                     srow_max[row_cp_off + icp]);
+               smin[row_cp_off + icp] = 0.5*t*w;
+               smax[row_cp_off + icp] = 1.5*t*w*x;
+            }
+            else
+            {
+               smin[row_cp_off + icp] = 0.0;
+               smax[row_cp_off + icp] = 0.0;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2b: reduce the y-direction projection coefficients for each
+      // x-control-point column.
+      MFEM_FOREACH_THREAD(jj, y, 1)
+      {
+         MFEM_FOREACH_THREAD(icp, x, ncp)
+         {
+            real_t a0 = 0.0;
+            real_t a1 = 0.0;
+            for (int jrow = 0; jrow < nb; jrow++)
+            {
+               a0 += smin[jrow*ncp + icp];
+               a1 += smax[jrow*ncp + icp];
+            }
+            sa0[icp] = a0;
+            sa1[icp] = a1;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2c: subtract the y-direction linear fit from the intermediate
+      // row bounds so the final tensor-product bound uses the perturbation.
+      if constexpr (T_PROJ)
+      {
+         MFEM_FOREACH_THREAD(icp, x, ncp)
+         {
+            MFEM_FOREACH_THREAD(jrow, y, nb)
+            {
+               const int row_cp_off = jrow*ncp;
+               const real_t x = data.xhat[jrow];
+               const real_t t = sa0[icp] + sa1[icp]*x;
+               srow_min[row_cp_off + icp] -= t;
+               srow_max[row_cp_off + icp] -= t;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 3: each thread now owns one 2D control point (icp, kcp) and
+      // accumulates its final lower/upper bound from the row-bound data.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(kcp, y, ncp)
+         {
+            real_t lo = 0.0;
+            real_t hi = 0.0;
+            if constexpr (T_PROJ)
+            {
+               const real_t xcp = data.cphat[kcp];
+               lo = sa0[icp] + sa1[icp]*xcp;
+               hi = lo;
+            }
+            for (int jrow = 0; jrow < nb; jrow++)
+            {
+               const real_t w0 = srow_min[jrow*ncp + icp];
+               const real_t w1 = srow_max[jrow*ncp + icp];
+               const real_t lb = data.lbound[kcp + jrow*data.ncp];
+               const real_t ub = data.ubound[kcp + jrow*data.ncp];
+               const real_t v0 = lb*w0;
+               const real_t v1 = ub*w0;
+               const real_t v2 = lb*w1;
+               const real_t v3 = ub*w1;
+               real_t vlo = v0 < v1 ? v0 : v1;
+               real_t vhi = v0 > v1 ? v0 : v1;
+               vlo = vlo < v2 ? vlo : v2;
+               vlo = vlo < v3 ? vlo : v3;
+               vhi = vhi > v2 ? vhi : v2;
+               vhi = vhi > v3 ? vhi : v3;
+               lo += vlo;
+               hi += vhi;
+            }
+            const int slot = kcp*ncp + icp;
+            smin[slot] = lo;
+            smax[slot] = hi;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      const int lane = ty*ncp + tx;
+      const int nactive = ncp*ncp;
+      const int nthreads = MFEM_THREAD_SIZE(x)*MFEM_THREAD_SIZE(y);
+
+      // Reduce all 2D control-point bounds to one lower/upper pair per
+      // (element, component).
+      if (nthreads == 1)
+      {
+         if (tx == 0 && ty == 0)
+         {
+            real_t lower_ec = smin[0];
+            real_t upper_ec = smax[0];
+            for (int t = 1; t < nactive; t++)
+            {
+               lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
+               upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
+            }
+            L(e, c) = lower_ec;
+            U(e, c) = upper_ec;
+         }
+      }
+      else
+      {
+         for (int stride = (nactive + 1)/2; stride > 0;
+              stride = (stride + 1)/2)
+         {
+            if (lane < stride && lane + stride < nactive)
+            {
+               smin[lane] = smin[lane] < smin[lane + stride] ?
+                            smin[lane] : smin[lane + stride];
+               smax[lane] = smax[lane] > smax[lane + stride] ?
+                            smax[lane] : smax[lane + stride];
+            }
+            MFEM_SYNC_THREAD;
+            if (stride == 1) { break; }
+         }
+
+         if (lane == 0)
+         {
+            L(e, c) = smin[0];
+            U(e, c) = smax[0];
+         }
+      }
+   });
+}
+
+} // namespace internal
+
+inline void PLBound::GetElementBoundsKernel(const int rdim, const int fes_vdim,
+                                            const Vector &e_vec,
+                                            Vector &lower, Vector &upper,
+                                            const int vdim) const
+{
+   MFEM_VERIFY(b_type != BasisType::Positive,
+               "Bernstein device bounds are not implemented.");
+   if (rdim == 3)
+   {
+      MFEM_ABORT("Device element bounds kernel currently only supports 1D/2D.");
+   }
+   MFEM_VERIFY(rdim == 1 || rdim == 2, "Invalid element dimension.");
+   MFEM_VERIFY(vdim >= -1 && vdim <= fes_vdim, "Invalid vector component.");
+   const int nd = static_cast<int>(std::pow(nb, rdim));
+   const int ne = e_vec.Size()/(nd*fes_vdim);
+   const int ncomp = (vdim > 0) ? 1 : fes_vdim;
+
+   lower.SetSize(ne*ncomp, e_vec);
+   upper.SetSize(ne*ncomp, e_vec);
+   lower.UseDevice(true);
+   upper.UseDevice(true);
+
+   if (!proj)
+   {
+      MFEM_ABORT("Device element bounds kernel currently requires projection "
+                 "enabled.");
+   }
+
+   const real_t *dxhat = xhat.Read();
+   const real_t *dwhat = what.Read();
+   const real_t *dcphat = cphat.Read();
+   const real_t *dlbound = lbound.Read();
+   const real_t *dubound = ubound.Read();
+
+   internal::PLBoundDeviceData data
+   {
+      nb,
+      ncp,
+      dxhat,
+      dwhat,
+      dcphat,
+      dlbound,
+      dubound
+   };
+
+   const int comp0 = (vdim > 0) ? (vdim - 1) : 0;
+
+   if (rdim == 1)
+   {
+      switch (nb)
+      {
+         case 2: return internal::GetElementBoundsKernel1D<2, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 3: return internal::GetElementBoundsKernel1D<3, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 4: return internal::GetElementBoundsKernel1D<4, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 5: return internal::GetElementBoundsKernel1D<5, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 6: return internal::GetElementBoundsKernel1D<6, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 7: return internal::GetElementBoundsKernel1D<7, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 8: return internal::GetElementBoundsKernel1D<8, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 9: return internal::GetElementBoundsKernel1D<9, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 10: return internal::GetElementBoundsKernel1D<10, true>(data, fes_vdim, ne,
+                                                                       e_vec, lower, upper,
+                                                                       comp0, ncomp);
+         default: return internal::GetElementBoundsKernel1D<0, true>(data, fes_vdim, ne,
+                                                                      e_vec, lower, upper,
+                                                                      comp0, ncomp);
+      }
+   }
+#define MFEM_PLBOUND_2D_DISPATCH(NB, NCP) \
+   return internal::GetElementBoundsKernel2D<NB, NCP, true>(data, fes_vdim, ne, \
+                                                            e_vec, lower, upper, \
+                                                            comp0, ncomp)
+   switch (nb)
+   {
+      case 2:
+         switch (ncp)
+         {
+            case 4: MFEM_PLBOUND_2D_DISPATCH(2, 4);
+            case 6: MFEM_PLBOUND_2D_DISPATCH(2, 6);
+            case 8: MFEM_PLBOUND_2D_DISPATCH(2, 8);
+         }
+         break;
+      case 3:
+         switch (ncp)
+         {
+            case 6: MFEM_PLBOUND_2D_DISPATCH(3, 6);
+            case 9: MFEM_PLBOUND_2D_DISPATCH(3, 9);
+            case 12: MFEM_PLBOUND_2D_DISPATCH(3, 12);
+         }
+         break;
+      case 4:
+         switch (ncp)
+         {
+            case 8: MFEM_PLBOUND_2D_DISPATCH(4, 8);
+            case 12: MFEM_PLBOUND_2D_DISPATCH(4, 12);
+            case 16: MFEM_PLBOUND_2D_DISPATCH(4, 16);
+         }
+         break;
+      case 5:
+         switch (ncp)
+         {
+            case 10: MFEM_PLBOUND_2D_DISPATCH(5, 10);
+            case 15: MFEM_PLBOUND_2D_DISPATCH(5, 15);
+            case 20: MFEM_PLBOUND_2D_DISPATCH(5, 20);
+         }
+         break;
+      case 6:
+         switch (ncp)
+         {
+            case 12: MFEM_PLBOUND_2D_DISPATCH(6, 12);
+            case 18: MFEM_PLBOUND_2D_DISPATCH(6, 18);
+            case 24: MFEM_PLBOUND_2D_DISPATCH(6, 24);
+         }
+         break;
+      case 7:
+         switch (ncp)
+         {
+            case 14: MFEM_PLBOUND_2D_DISPATCH(7, 14);
+            case 21: MFEM_PLBOUND_2D_DISPATCH(7, 21);
+            case 28: MFEM_PLBOUND_2D_DISPATCH(7, 28);
+         }
+         break;
+      case 8:
+         switch (ncp)
+         {
+            case 16: MFEM_PLBOUND_2D_DISPATCH(8, 16);
+            case 24: MFEM_PLBOUND_2D_DISPATCH(8, 24);
+            case 32: MFEM_PLBOUND_2D_DISPATCH(8, 32);
+         }
+         break;
+   }
+#undef MFEM_PLBOUND_2D_DISPATCH
+   return internal::GetElementBoundsKernel2D<0, 0, true>(data, fes_vdim, ne,
+                                                         e_vec, lower, upper,
+                                                         comp0, ncomp);
+}
+
 } // namespace mfem
 
 #endif // MFEM_BOUNDS

fem/coefficient.hpp

@@ -52,6 +52,9 @@ public:
    /// Get the time for time dependent coefficients
    real_t GetTime() { return time; }
 
+   /// Returns dimension of the vector.
+   int GetVDim() { return 1; }
+
    /** @brief Evaluate the coefficient in the element described by @a T at the
        point @a ip. */
    /** @note When this method is called, the caller must make sure that the

fem/complex_fem.hpp

@@ -82,6 +82,25 @@ public:
    /// underlying #fes
    int VectorDim() const;
 
+   /// Copy assignment. Only the data of the base class Vector is copied.
+   /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
+       have the same size.
+
+       @note Defining this method overwrites the implicitly defined copy
+       assignment operator. */
+   ComplexGridFunction &operator=(const ComplexGridFunction &rhs)
+   { return operator=((const Vector &)rhs); }
+
+   /// Copy the data from @a v.
+   /** The size of @a v must be equal to double of the size of the associated
+       FiniteElementSpace #fes. */
+   ComplexGridFunction &operator=(const Vector &v)
+   {
+      MFEM_ASSERT(fes && v.Size() == 2*fes->GetVSize(), "");
+      Vector::operator=(v);
+      return *this;
+   }
+
    /// Assign constant values to the ComplexGridFunction data.
    ComplexGridFunction &operator=(const std::complex<real_t> & value)
    { *gfr = value.real(); *gfi = value.imag(); return *this; }

fem/datacollection.cpp

@@ -492,6 +492,8 @@ void VisItDataCollection::SaveRootFile()
                            to_padded_string(cycle, pad_digits_cycle) +
                            ".mfem_root";
    std::ofstream root_file(root_name);
+   MFEM_VERIFY(root_file.is_open(),
+               "Failed to open ofstream " << root_name);
    root_file << GetVisItRootString();
    if (!root_file)
    {
@@ -977,7 +979,10 @@ void ParaViewDataCollection::Save()
    // Save the local part of the mesh and grid functions fields to the local
    // VTU file. Also save coefficient fields.
    {
-      std::ofstream os(vtu_prefix + GenerateVTUFileName("proc", myid));
+      std::string os_str = vtu_prefix + GenerateVTUFileName("proc", myid);
+      std::ofstream os(os_str);
+      MFEM_VERIFY(os.is_open(),
+                  "Failed to open ofstream " << os_str);
       os.precision(precision);
       SaveDataVTU(os, levels_of_detail);
    }
@@ -989,7 +994,10 @@ void ParaViewDataCollection::Save()
                   "QuadratureFunction output is not supported for "
                   "ParaViewDataCollection on domain boundary!");
       const std::string &field_name = qfield.first;
-      std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
+      std::string os_str = vtu_prefix + GenerateVTUFileName(field_name, myid);
+      std::ofstream os(os_str);
+      MFEM_VERIFY(os.is_open(),
+                  "Failed to open ofstream " << os_str);
       qfield.second->SaveVTU(os, pv_data_format, GetCompressionLevel(), field_name);
    }
 
@@ -1000,7 +1008,10 @@ void ParaViewDataCollection::Save()
    {
       // Create the main PVTU file
       {
-         std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
+         std::string os_str = vtu_prefix + GeneratePVTUFileName("data");
+         std::ofstream pvtu_out(os_str);
+         MFEM_VERIFY(pvtu_out.is_open(),
+                     "Failed to open ofstream " << os_str);
          WritePVTUHeader(pvtu_out);
 
          // Grid function fields and coefficient fields
@@ -1055,8 +1066,10 @@ void ParaViewDataCollection::Save()
          const std::string &q_field_name = q_field.first;
          std::string q_fname = GeneratePVTUPath() + "/"
                                + GeneratePVTUFileName(q_field_name);
-
-         std::ofstream pvtu_out(col_path + "/" + q_fname);
+         std::string os_str = col_path + "/" + q_fname;
+         std::ofstream pvtu_out(os_str);
+         MFEM_VERIFY(pvtu_out.is_open(),
+                     "Failed to open ofstream " << os_str);
          WritePVTUHeader(pvtu_out);
          int vec_dim = q_field.second->GetVDim();
          pvtu_out << "<PPointData>\n";

fem/dfem/integrate.hpp

@@ -90,8 +90,8 @@ void map_quadrature_data_to_fields_impl(
    }
    else
    {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        " for this field descriptor");
    }
 }
 
@@ -169,8 +169,9 @@ void map_quadrature_data_to_fields_tensor_impl_1d(
    }
    else
    {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor with sum factorization on tensor product elements");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        "for this field descriptor with sum factorization on"
+                        " tensor product elements");
    }
 }
 
@@ -306,8 +307,9 @@ void map_quadrature_data_to_fields_tensor_impl_2d(
    }
    else
    {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor with sum factorization on tensor product elements");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        " for this field descriptor with sum factorization on"
+                        " tensor product elements");
    }
 }
 
@@ -492,8 +494,9 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
    }
    else
    {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor with sum factorization on tensor product elements");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        " for this field descriptor with sum factorization on"
+                        " tensor product elements");
    }
 }
 

fem/dgmassinv_kernels.hpp

@@ -57,7 +57,7 @@ void DGMassApply(const int e,
       }
       else if (DIM == 3)
       {
-         SmemPAMassApply3D_Element<TD1D,TQ1D,ACCUM>(e, NE, B, pa_data, x, y);
+         SmemPAMassApply3D_Element<TD1D,TQ1D,NBZ,ACCUM>(e, NE, B, pa_data, x, y);
       }
       else
       {

fem/estimators.hpp

@@ -320,8 +320,8 @@ public:
     error estimation procedure where the flux averaging is replaced by a global
     L2 projection (requiring a mass matrix solve).
 
-    The required BilinearFormIntegrator must implement the methods
-    ComputeElementFlux() and ComputeFluxEnergy().
+    The required BilinearFormIntegrator must implement the method
+    ComputeElementFlux().
 
     Implemented for the parallel case only.
  */
@@ -357,8 +357,8 @@ protected:
 
 public:
    /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the methods
-                       ComputeElementFlux() and ComputeFluxEnergy().
+       @param integ    This BilinearFormIntegrator must implement the method
+                       ComputeElementFlux().
        @param sol      The solution field whose error is to be estimated.
        @param flux_fes The L2ZienkiewiczZhuEstimator assumes ownership of this
                        FiniteElementSpace and will call its Update() method when
@@ -382,8 +382,8 @@ public:
    { }
 
    /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the methods
-                       ComputeElementFlux() and ComputeFluxEnergy().
+       @param integ    This BilinearFormIntegrator must implement the method
+                       ComputeElementFlux().
        @param sol      The solution field whose error is to be estimated.
        @param flux_fes The L2ZienkiewiczZhuEstimator does NOT assume ownership
                        of this FiniteElementSpace; will call its Update() method

fem/fe/fe_base.cpp

@@ -663,58 +663,59 @@ const
    #pragma omp critical (DofToQuad)
 #endif
    {
-      // If the new Dof2Quad is already present, e.g. added in a previous call
-      // or added by another omp thread, return.
+      // Do not run if the new Dof2Quad is already present, e.g. added in a
+      // previous call or added by another omp thread.
       if (DofToQuad::SearchArray(dof2quad_array, ir,
-                                 DofToQuad::LEXICOGRAPHIC_FULL))
-      { return; }
-
-      // Undo the native ordering which is what FiniteElement::GetDofToQuad
-      // returns.
-      auto *d2q_new = new DofToQuad(d2q);
-      d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
-      const int nqpt = ir.GetNPoints();
-
-      const int b_dim = (range_type == VECTOR) ? dim : 1;
-
-      for (int i = 0; i < nqpt; i++)
+                                 DofToQuad::LEXICOGRAPHIC_FULL) == nullptr)
       {
-         for (int d = 0; d < b_dim; d++)
+         // Undo the native ordering which is what FiniteElement::GetDofToQuad
+         // returns.
+         auto *d2q_new = new DofToQuad(d2q);
+         d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
+         const int nqpt = ir.GetNPoints();
+
+         const int b_dim = (range_type == VECTOR) ? dim : 1;
+
+         for (int i = 0; i < nqpt; i++)
          {
-            for (int j = 0; j < dof; j++)
+            for (int d = 0; d < b_dim; d++)
             {
-               const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
-               d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
-               d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
+                  d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
+                  d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               }
             }
          }
-      }
 
-      const int g_dim = [this]()
-      {
-         switch (deriv_type)
+         const int g_dim = [this]()
          {
-            case GRAD: return dim;
-            case DIV: return 1;
-            case CURL: return cdim;
-            default: return 0;
-         }
-      }();
-
-      for (int i = 0; i < nqpt; i++)
-      {
-         for (int d = 0; d < g_dim; d++)
-         {
-            for (int j = 0; j < dof; j++)
+            switch (deriv_type)
             {
-               const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
-               d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
-               d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               case GRAD: return dim;
+               case DIV: return 1;
+               case CURL: return cdim;
+               default: return 0;
+            }
+         }();
+
+         for (int i = 0; i < nqpt; i++)
+         {
+            for (int d = 0; d < g_dim; d++)
+            {
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
+                  d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
+                  d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               }
             }
          }
+
+         dof2quad_array.Append(d2q_new);
       }
 
-      dof2quad_array.Append(d2q_new);
    }
 }
 
@@ -1043,9 +1044,50 @@ void VectorFiniteElement::SetDerivMembers()
    switch (map_type)
    {
       case H_DIV:
-         deriv_type = DIV;
-         deriv_range_type = SCALAR;
-         deriv_map_type = INTEGRAL;
+         switch (dim)
+         {
+            case 3: // div: 3D H_DIV -> 3D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            case 2: // div: 2D H_DIV -> 2D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_DIV_R2D:
+         switch (dim)
+         {
+            case 2: // div: 2D H_DIV_R2D -> 2D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            case 1: // div: 1D H_DIV_R2D -> 1D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_DIV_R1D:
+         switch (dim)
+         {
+            case 1: // div: 1D H_DIV_R1D -> 1D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
          break;
       case H_CURL:
          switch (dim)
@@ -1063,13 +1105,49 @@ void VectorFiniteElement::SetDerivMembers()
                break;
             case 1:
                deriv_type = NONE;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
+               deriv_range_type = UNKNOWN_RANGE_TYPE;
+               deriv_map_type = UNKNOWN_MAP_TYPE;
                break;
             default:
                MFEM_ABORT("Invalid dimension, Dim = " << dim);
          }
          break;
+      case H_CURL_R2D:
+         switch (dim)
+         {
+            case 2:
+               // curl: 2D H_CURL_R2D -> H_DIV_R2D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R2D;
+               break;
+            case 1:
+               // curl: 1D H_CURL_R2D -> H_DIV_R2D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R2D;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_CURL_R1D:
+         switch (dim)
+         {
+            case 1:
+               // curl: 1D H_CURL_R1D -> H_DIV_R1D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R1D;
+               break;
+            case 0:
+               deriv_type = NONE;
+               deriv_range_type = UNKNOWN_RANGE_TYPE;
+               deriv_map_type = UNKNOWN_MAP_TYPE;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
       default:
          MFEM_ABORT("Invalid MapType = " << map_type);
    }

fem/fe/fe_base.hpp

@@ -295,10 +295,20 @@ public:
                           $ u(x) = (1/w) \hat u(\hat x) $ */
       H_DIV,     /**< For vector fields; preserves surface integrals of the
                           normal component $ u(x) = (J/w) \hat u(\hat x) $ */
-      H_CURL     /**< For vector fields; preserves line integrals of the
+      H_CURL,    /**< For vector fields; preserves line integrals of the
                           tangential component
                           $ u(x) = J^{-t} \hat u(\hat x) $ (square J),
                           $ u(x) = J(J^t J)^{-1} \hat u(\hat x) $ (general J) */
+      H_DIV_R2D, /**< For 3-component vector fields in 2D; equivalent to a
+                          direct sum of an H_DIV basis and an INTEGRAL basis */
+      H_CURL_R2D,/**< For 3-component vector fields in 2D; equivalent to a
+                          direct sum of an H_CURL basis and a VALUE basis */
+      H_DIV_R1D, /**< For 3-component vector fields in 1D; equivalent to a
+                          direct sum of a VALUE basis and a pair of INTEGRAL
+                          bases */
+      H_CURL_R1D /**< For 3-component vector fields in 1D; equivalent to a
+                          direct sum of an INTEGRAL basis and a pair of VALUE
+                          bases */
    };
 
    /** @brief Enumeration for DerivType: defines which derivative method
@@ -330,12 +340,28 @@ public:
    int GetDim() const { return dim; }
 
    /** @brief Returns the vector dimension for vector-valued finite elements,
-       which is also the dimension of the interpolation operation. */
+       which is also the dimension of the interpolation operation and the
+       width of the DenseMatrix argument in
+       CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape). */
    int GetRangeDim() const { return vdim; }
 
-   /// Returns the dimension of the curl for vector-valued finite elements.
+   /** @brief Returns the vector dimension, in physical space, for
+       vector-valued finite elements, which is also the width of the
+       DenseMatrix argument in
+       CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape). */
+   virtual int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
+
+   /** Returns the dimension of the curl for vector-valued finite elements,
+       which is also the width of the DenseMatrix argument in
+       CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape). */
    int GetCurlDim() const { return cdim; }
 
+   /** Returns the dimension, in physical space, of the curl for vector-valued
+       finite elements, which is also the width of the DenseMatrix argument in
+       CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape).
+   */
+   virtual int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
+
    /// Returns the Geometry::Type of the reference element.
    Geometry::Type GetGeomType() const { return geom_type; }
 
@@ -990,6 +1016,8 @@ protected:
 public:
    VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M,
                        int F = FunctionSpace::Pk);
+
+   int GetPhysRangeDim(int space_dim) const override { return space_dim; }
 };
 
 /// @brief Class for computing 1D special polynomials and their associated basis

fem/fe/fe_h1.cpp

@@ -589,7 +589,7 @@ void H1_TriangleElement::CalcHessian(const IntegrationPoint &ip,
    Vector   shape_x(p + 1),   shape_y(p + 1),   shape_l(p + 1);
    Vector  dshape_x(p + 1),  dshape_y(p + 1),  dshape_l(p + 1);
    Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
-   DenseMatrix ddu(dof, dim);
+   DenseMatrix ddu(dof, (dim*(dim+1))/2);
 #endif
 
    poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);

fem/fe/fe_nd.cpp

@@ -2531,7 +2531,7 @@ void ND_FuentesPyramidElement::calcCurlBasis(const int p,
 
 ND_R1D_PointElement::ND_R1D_PointElement(int p)
    : VectorFiniteElement(1, Geometry::POINT, 2, p,
-                         H_CURL, FunctionSpace::Pk)
+                         H_CURL_R1D, FunctionSpace::Pk)
 {
    // VectorFiniteElement::SetDerivMembers doesn't support 0D H_CURL elements
    // so we mimic a 1D element and then correct the dimension here.
@@ -2562,7 +2562,7 @@ ND_R1D_SegmentElement::ND_R1D_SegmentElement(const int p,
                                              const int cb_type,
                                              const int ob_type)
    : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 2, p,
-                         H_CURL, FunctionSpace::Pk),
+                         H_CURL_R1D, FunctionSpace::Pk),
      dof2tk(dof),
      cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
      obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -2839,7 +2839,7 @@ ND_R2D_SegmentElement::ND_R2D_SegmentElement(const int p,
                                              const int cb_type,
                                              const int ob_type)
    : VectorFiniteElement(1, Geometry::SEGMENT, 2 * p + 1, p,
-                         H_CURL, FunctionSpace::Pk),
+                         H_CURL_R2D, FunctionSpace::Pk),
      dof2tk(dof),
      cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
      obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -3023,7 +3023,7 @@ void ND_R2D_SegmentElement::Project(VectorCoefficient &vc,
 ND_R2D_FiniteElement::ND_R2D_FiniteElement(int p, Geometry::Type G, int Do,
                                            const real_t *tk_fe)
    : VectorFiniteElement(2, G, Do, p,
-                         H_CURL, FunctionSpace::Pk),
+                         H_CURL_R2D, FunctionSpace::Pk),
      tk(tk_fe),
      dof_map(dof),
      dof2tk(dof)

fem/fe/fe_nd.hpp

@@ -663,6 +663,9 @@ public:
                          const int cb_type = BasisType::GaussLobatto,
                          const int ob_type = BasisType::GaussLegendre);
 
+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { return 1; }
+
    void CalcVShape(const IntegrationPoint &ip,
                    DenseMatrix &shape) const override;
 
@@ -705,6 +708,9 @@ private:
                            DenseMatrix &I) const;
 
 public:
+   int GetPhysRangeDim(int space_dim) const override { return 3; }
+   int GetPhysCurlDim(int space_dim) const override { return 3; }
+
    using FiniteElement::CalcVShape;
    using FiniteElement::CalcPhysCurlShape;
 

fem/fe/fe_rt.cpp

@@ -2006,7 +2006,7 @@ RT_R1D_SegmentElement::RT_R1D_SegmentElement(const int p,
                                              const int cb_type,
                                              const int ob_type)
    : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 4, p + 1,
-                         H_DIV, FunctionSpace::Pk),
+                         H_DIV_R1D, FunctionSpace::Pk),
      dof2nk(dof),
      cbasis1d(poly1d.GetBasis(p + 1, VerifyClosed(cb_type))),
      obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
@@ -2281,7 +2281,7 @@ const real_t RT_R2D_SegmentElement::nk[2] = { 0.,1.};
 RT_R2D_SegmentElement::RT_R2D_SegmentElement(const int p,
                                              const int ob_type)
    : VectorFiniteElement(1, Geometry::SEGMENT, p + 1, p + 1,
-                         H_DIV, FunctionSpace::Pk),
+                         H_DIV_R2D, FunctionSpace::Pk),
      dof2nk(dof),
      obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
 {
@@ -2392,7 +2392,7 @@ void RT_R2D_SegmentElement::LocalInterpolation(const VectorFiniteElement &cfe,
 RT_R2D_FiniteElement::RT_R2D_FiniteElement(int p, Geometry::Type G, int Do,
                                            const real_t *nk_fe)
    : VectorFiniteElement(2, G, Do, p + 1,
-                         H_DIV, FunctionSpace::Pk),
+                         H_DIV_R2D, FunctionSpace::Pk),
      nk(nk_fe),
      dof_map(dof),
      dof2nk(dof)

fem/fe/fe_rt.hpp

@@ -510,6 +510,9 @@ public:
    RT_R2D_SegmentElement(const int p,
                          const int ob_type = BasisType::GaussLegendre);
 
+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { return 0; }
+
    void CalcVShape(const IntegrationPoint &ip,
                    DenseMatrix &shape) const override;
 
@@ -547,6 +550,9 @@ private:
                            DenseMatrix &I) const;
 
 public:
+   int GetPhysRangeDim(int space_dim) const override { return 3; }
+   int GetPhysCurlDim(int space_dim) const override { return 0; }
+
    using FiniteElement::CalcVShape;
 
    void CalcVShape(ElementTransformation &Trans,

fem/fespace.cpp

@@ -282,14 +282,7 @@ int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs_) const
 void FiniteElementSpace::AdjustVDofs(Array<int> &vdofs)
 {
    int n = vdofs.Size(), *vdof = vdofs;
-   for (int i = 0; i < n; i++)
-   {
-      int j;
-      if ((j = vdof[i]) < 0)
-      {
-         vdof[i] = -1-j;
-      }
-   }
+   for (int i = 0; i < n; i++) { vdof[i] = UnsignIndex(vdof[i]); }
 }
 
 void FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs,
@@ -483,13 +476,14 @@ void FiniteElementSpace::ReorderElementToDofTable()
    for (int k = 0, dof_counter = 0; k < nnz; k++)
    {
       const int sdof = J[k]; // signed dof
-      const int dof = (sdof < 0) ? -1-sdof : sdof;
+      const int dof = UnsignIndex(sdof);
       int new_dof = dof_marker[dof];
       if (new_dof < 0)
       {
          dof_marker[dof] = new_dof = dof_counter++;
       }
-      J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
+      // Preserve the sign of sdof
+      J[k] = (sdof < 0) ? FlipIndexSign(new_dof) : new_dof;
    }
 }
 
@@ -547,7 +541,7 @@ void MarkDofs(const Array<int> &dofs, Array<int> &mark_array)
 {
    for (auto d : dofs)
    {
-      mark_array[d >= 0 ? d : -1 - d] = -1;
+      mark_array[UnsignIndex(d)] = -1;
    }
 }
 
@@ -931,7 +925,7 @@ void FiniteElementSpace::AddDependencies(
             if (std::abs(coef) > 1e-12)
             {
                const int mdof = master_dofs[j];
-               if (mdof != sdof && mdof != (-1-sdof))
+               if (mdof != sdof && mdof != FlipIndexSign(sdof))
                {
                   deps.Add(sdof, mdof, coef);
                }
@@ -1024,7 +1018,7 @@ int FiniteElementSpace::GetDegenerateFaceDofs(int index, Array<int> &dofs,
    // FiniteElementSpace::AddDependencies.
 
    Array<int> edof;
-   int order = GetEdgeDofs(-1 - index, edof, variant);
+   int order = GetEdgeDofs(FlipIndexSign(index), edof, variant);
 
    int nv = fec->DofForGeometry(Geometry::POINT);
    int ne = fec->DofForGeometry(Geometry::SEGMENT);
@@ -1710,8 +1704,8 @@ SparseMatrix *FiniteElementSpace::RefinementMatrix_main(
 
          for (int i = 0; i < fine_ldof; i++)
          {
-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);
 
             if (!mark[m])
             {
@@ -1772,7 +1766,7 @@ SparseMatrix *FiniteElementSpace::VariableOrderRefinementMatrix(
          for (int i = 0; i < fine_ldof; i++)
          {
             const int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int m = UnsignIndex(r);
 
             if (!mark[m])
             {
@@ -2482,8 +2476,8 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
          {
             if (!std::isfinite(lR(i, 0))) { continue; }
 
-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);
 
             if (is_dg || !mark[m])
             {
@@ -3201,7 +3195,7 @@ void FiniteElementSpace::CalcEdgeFaceVarOrders(
             else
             {
                // degenerate face (i.e., edge-face constraint)
-               slave_orders |= edge_orders[-1 - slave.index];
+               slave_orders |= edge_orders[FlipIndexSign(slave.index)];
             }
          }
 
@@ -3940,6 +3934,16 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const
    return BE;
 }
 
+const FiniteElement *FiniteElementSpace::GetTypicalBE() const
+{
+   if (mesh->GetNBE() > 0) { return GetBE(0); }
+
+   Geometry::Type geom = mesh->GetTypicalFaceGeometry();
+   const FiniteElement *be = fec->FiniteElementForGeometry(geom);
+   MFEM_VERIFY(be != nullptr, "Could not determine a typical BE!");
+   return be;
+}
+
 const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
 {
    MFEM_VERIFY(!IsVariableOrder(), "not implemented");
@@ -3970,6 +3974,11 @@ const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
    return fe;
 }
 
+const FiniteElement *FiniteElementSpace::GetTypicalFaceElement() const
+{
+   return fec->FiniteElementForGeometry(mesh->GetTypicalFaceGeometry());
+}
+
 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
                                                         int variant) const
 {

fem/fespace.hpp

@@ -839,7 +839,7 @@ public:
        Note: For vector-valued elements, the results pads up the range dimension
        to the spatial dimension. E.g., consider a stack of 5 vector-valued
        elements each representing 2D vectors, living in a 3 dimensional space.
-       Then this fucntion would give 15, not 10.
+       Then this function would give 15, not 10.
        */
    int GetVectorDim() const;
 
@@ -1150,7 +1150,7 @@ public:
 
    /// Helper to return the DOF associated with a sign encoded DOF
    static inline int DecodeDof(int dof)
-   { return (dof >= 0) ? dof : (-1 - dof); }
+   { return UnsignIndex(dof); }
 
    /// Helper to determine the DOF and sign of a sign encoded DOF
    static inline int DecodeDof(int dof, real_t& sign)
@@ -1323,12 +1323,24 @@ public:
         associated with i'th boundary face in the mesh object. */
    const FiniteElement *GetBE(int i) const;
 
+   /// @brief Return a "typical" boundary element.
+   ///
+   /// This can be used in situations where the local mesh partition may be
+   /// empty.
+   const FiniteElement *GetTypicalBE() const;
+
    /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
         associated with i'th face in the mesh object.  Faces in this case refer
         to the MESHDIM-1 primitive so in 2D they are segments and in 1D they are
         points.*/
    const FiniteElement *GetFaceElement(int i) const;
 
+   /// @brief Return a "typical" face element.
+   ///
+   /// This can be used in situations where the local mesh partition may be
+   /// empty.
+   const FiniteElement *GetTypicalFaceElement() const;
+
    /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
         associated with i'th edge in the mesh object. */
    const FiniteElement *GetEdgeElement(int i, int variant = 0) const;

fem/gridfunc.cpp

@@ -345,27 +345,6 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi,
    }
 }
 
-int GridFunction::VectorDim() const
-{
-   const FiniteElement *fe = fes->GetTypicalFE();
-   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
-   {
-      return fes->GetVDim();
-   }
-   return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
-                                  fe->GetRangeDim());
-}
-
-int GridFunction::CurlDim() const
-{
-   const FiniteElement *fe = fes->GetTypicalFE();
-   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
-   {
-      return 2 * fes->GetMesh()->SpaceDimension() - 3;
-   }
-   return fes->GetVDim()*fe->GetCurlDim();
-}
-
 void GridFunction::GetTrueDofs(Vector &tv) const
 {
    const SparseMatrix *R = fes->GetRestrictionMatrix();
@@ -2050,6 +2029,18 @@ void GridFunction::AccumulateAndCountBdrValues(
    Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
    Array<int> &values_counter)
 {
+   if (vcoeff)
+   {
+      MFEM_VERIFY(fes->GetVDim() == vcoeff->GetVDim(),
+                  "vcoeff vdim != fes VDim");
+      MFEM_VERIFY(fes->GetTypicalBE()->GetMapType() == FiniteElement::VALUE &&
+                  fes->GetTypicalBE()->GetRangeType() ==
+                  FiniteElement::SCALAR,
+                  "Can only call ProjectBdrCoefficient on scalar value-type "
+                  "boundary elements. "
+                  "Did you intended to call ProjectBdrCoefficientNormal or "
+                  "ProjectBdrCoefficientTangent for vector finite elements?");
+   }
    Array<int> vdofs;
    Vector vc;
 
@@ -2202,6 +2193,9 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
    VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
    Array<int> &values_counter)
 {
+   MFEM_VERIFY(fes->GetTypicalBE()->GetPhysRangeDim(
+                  fes->GetMesh()->SpaceDimension()) == vcoeff.GetVDim(),
+               "vcoeff vdim != PhysRangeDim");
    const FiniteElement *fe;
    ElementTransformation *T;
    Array<int> dofs;
@@ -2355,6 +2349,9 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
 
 void GridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
+   MFEM_VERIFY(
+      VectorDim() == 1,
+      "Cannot project scalar Coefficient onto vector GridFunction");
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
    DofTransformation doftrans;
    Array<int> vdofs;
@@ -2630,6 +2627,7 @@ void GridFunction::ProjectCoefficient(
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff,
                                       ProjectType type)
 {
+   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    Array<int> vdofs;
    Vector vals;
    DofTransformation doftrans;
@@ -2945,6 +2943,7 @@ void GridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 void GridFunction::ProjectCoefficient(
    VectorCoefficient &vcoeff, Array<int> &dofs)
 {
+   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int el = -1;
    ElementTransformation *T = NULL;
    const FiniteElement *fe = NULL;
@@ -2974,6 +2973,7 @@ void GridFunction::ProjectCoefficient(
 
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff, int attribute)
 {
+   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int i;
    Array<int> vdofs;
    Vector vals;
@@ -3030,9 +3030,14 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
    }
 }
 
-void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
-                                          Array<int> &dof_attr)
+void GridFunction::ProjectDiscCoefficient(
+   std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr)
 {
+   std::visit([&](auto* c)
+   {
+      MFEM_VERIFY(VectorDim() == c->GetVDim(), "coeff vdim != VectorDim()");
+   }, coeff);
+
    Array<int> vdofs;
    Vector vals;
 
@@ -3046,7 +3051,10 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
    {
       fes->GetElementVDofs(i, vdofs);
       vals.SetSize(vdofs.Size());
-      fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
+      std::visit([&](auto* c)
+      {
+         fes->GetFE(i)->Project(*c, *fes->GetElementTransformation(i), vals);
+      }, coeff);
 
       // the values in shared dofs are determined from the element with maximal
       // attribute
@@ -3062,17 +3070,15 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
    }
 }
 
-void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
-{
-   Array<int> dof_attr;
-   ProjectDiscCoefficient(coeff, dof_attr);
-}
-
 void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
+   MFEM_VERIFY(
+      VectorDim() == 1,
+      "Cannot project a scalar coefficient onto a vector GridFunction");
+
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -3082,6 +3088,7 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
                                           AvgType type)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -3137,52 +3144,33 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
 }
 
 void GridFunction::ProjectBdrCoefficientNormal(
-   VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
+   Coefficient *coeff, VectorCoefficient *vcoeff, const Array<int> &bdr_attr)
 {
-#if 0
-   // implementation for the case when the face dofs are integrals of the
-   // normal component.
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Array<int> dofs;
-   int dim = vcoeff.GetVDim();
-   Vector vc(dim), nor(dim), lvec, shape;
-
-   for (int i = 0; i < fes->GetNBE(); i++)
+   MFEM_VERIFY(fes->GetVDim() == 1, "fespace VDim != 1");
+   MFEM_VERIFY(fes->GetTypicalBE()->GetRangeType() == FiniteElement::SCALAR &&
+               fes->GetTypicalBE()->GetMapType() == FiniteElement::INTEGRAL,
+               "Not an RT FE space!");
+   if (vcoeff)
    {
-      if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
-      {
-         continue;
-      }
-      fe = fes->GetBE(i);
-      T = fes->GetBdrElementTransformation(i);
-      int intorder = 2*fe->GetOrder(); // !!!
-      const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
-      int nd = fe->GetDof();
-      lvec.SetSize(nd);
-      shape.SetSize(nd);
-      lvec = 0.0;
-      for (int j = 0; j < ir.GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(j);
-         T->SetIntPoint(&ip);
-         vcoeff.Eval(vc, *T, ip);
-         CalcOrtho(T->Jacobian(), nor);
-         fe->CalcShape(ip, shape);
-         lvec.Add(ip.weight * (vc * nor), shape);
-      }
-      fes->GetBdrElementDofs(i, dofs);
-      SetSubVector(dofs, lvec);
+      MFEM_VERIFY(vcoeff->GetVDim() == fes->GetMesh()->SpaceDimension(),
+                  "vcoeff vdim (" << vcoeff->GetVDim()
+                  << ") != SpaceDimension ("
+                  << fes->GetMesh()->SpaceDimension() << ")");
    }
-#else
+
    // implementation for the case when the face dofs are scaled point
    // values of the normal component.
    const FiniteElement *fe;
    ElementTransformation *T;
    Array<int> dofs;
-   int dim = vcoeff.GetVDim();
-   Vector vc(dim), nor(dim), lvec;
+   Vector vc, nor, lvec;
    DofTransformation doftrans;
+   if (vcoeff)
+   {
+      const int dim = vcoeff->GetVDim();
+      vc.SetSize(dim);
+      nor.SetSize(dim);
+   }
 
    for (int i = 0; i < fes->GetNBE(); i++)
    {
@@ -3198,15 +3186,22 @@ void GridFunction::ProjectBdrCoefficientNormal(
       {
          const IntegrationPoint &ip = ir.IntPoint(j);
          T->SetIntPoint(&ip);
-         vcoeff.Eval(vc, *T, ip);
-         CalcOrtho(T->Jacobian(), nor);
-         lvec(j) = (vc * nor);
+         if (coeff)
+         {
+            const real_t c = coeff->Eval(*T, ip);
+            lvec(j) = c * T->Weight();
+         }
+         else if (vcoeff)
+         {
+            vcoeff->Eval(vc, *T, ip);
+            CalcOrtho(T->Jacobian(), nor);
+            lvec(j) = (vc * nor);
+         }
       }
       fes->GetBdrElementDofs(i, dofs, doftrans);
       doftrans.TransformPrimal(lvec);
       SetSubVector(dofs, lvec);
    }
-#endif
 }
 
 void GridFunction::ProjectBdrCoefficientTangent(
@@ -5007,6 +5002,14 @@ real_t ExtrudeCoefficient::Eval(ElementTransformation &T,
    return sol_in.Eval(*T_in, ip);
 }
 
+void VectorExtrudeCoefficient::Eval(Vector &v, ElementTransformation &T,
+                                    const IntegrationPoint &ip)
+{
+   ElementTransformation *T_in =
+      mesh_in->GetElementTransformation(T.ElementNo / n);
+   T_in->SetIntPoint(&ip);
+   sol_in.Eval(v, *T_in, ip);
+}
 
 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                     GridFunction *sol, const int ny)
@@ -5057,10 +5060,17 @@ GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
       return NULL;
    }
    FiniteElementSpace *solfes2d;
-   // assuming sol is scalar
-   solfes2d = new FiniteElementSpace(mesh2d, solfec2d);
+   const int vdim = sol->FESpace()->GetVDim();
+   solfes2d = new FiniteElementSpace(mesh2d, solfec2d, vdim);
    sol2d = new GridFunction(solfes2d);
    sol2d->MakeOwner(solfec2d);
+   if (vdim > 1)
+   {
+      VectorGridFunctionCoefficient vcsol(sol);
+      VectorExtrudeCoefficient vc2d(mesh, vcsol, ny);
+      sol2d->ProjectCoefficient(vc2d);
+   }
+   else
    {
       GridFunctionCoefficient csol(sol);
       ExtrudeCoefficient c2d(mesh, csol, ny);
@@ -5242,6 +5252,30 @@ void GridFunction::GetElementBounds(const PLBound &plb,
                                     Vector &lower, Vector &upper,
                                     const int vdim) const
 {
+   if (UseDevice() && Device::Allows(Backend::DEVICE_MASK) &&
+       plb.GetBasisType() != BasisType::Positive &&
+       UsesTensorBasis(*fes))
+   {
+      const FiniteElement &fe = *fes->GetTypicalFE();
+      const int rdim = fe.GetDim();
+      const int fes_dim = fes->GetVDim();
+      const int nel = fes->GetNE();
+      const int nd = fe.GetDof();
+
+      Vector e_vec(nd*fes_dim*nel, Device::GetDeviceMemoryType());
+      e_vec.UseDevice(true);
+      const ElementRestrictionOperator *elem_restr =
+         fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      MFEM_VERIFY(elem_restr != nullptr,
+                  "Element restriction is required for device bounds.");
+      elem_restr->Mult(*this, e_vec);
+
+      plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower, upper, vdim);
+      lower.HostRead();
+      upper.HostRead();
+      return;
+   }
+
    int nel  = fes->GetNE();
    int fes_dim  = fes->GetVDim();
    lower.SetSize(nel*(vdim > 0 ? 1 :fes_dim));
@@ -5758,4 +5792,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
    return std::make_pair(global_max_lower, global_max_upper);
 }
 
-}
+}

fem/gridfunc.hpp

@@ -23,6 +23,7 @@
 #include <limits>
 #include <ostream>
 #include <string>
+#include <variant>
 
 namespace mfem
 {
@@ -79,10 +80,18 @@ protected:
                         bool wcoef,
                         int subdomain);
 
-   /** Project a discontinuous vector coefficient in a continuous space and
-       return in dof_attr the maximal attribute of the elements containing each
-       degree of freedom. */
-   void ProjectDiscCoefficient(VectorCoefficient &coeff, Array<int> &dof_attr);
+   /** @brief Project a discontinuous (vector) coefficient as a grid function on
+       a continuous finite element space. Return in dof_attr the maximal
+       attribute of the elements containing each degree of freedom. */
+   virtual void ProjectDiscCoefficient(
+      std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr);
+
+   /** @brief Project a discontinuous (vector) coefficient as a grid function on
+       a continuous finite element space. The values in shared dofs are
+       determined from the element with maximal attribute. */
+   virtual void ProjectDiscCoefficient(
+      std::variant<Coefficient*, VectorCoefficient*> coeff)
+   { Array<int> dof_attr; ProjectDiscCoefficient(coeff, dof_attr); };
 
    /** Helper function for ProjectCoefficientElementL2 */
    void ProjectCoefficientElementL2_(Coefficient &coeff, Vector &sol, Vector &Va);
@@ -150,11 +159,13 @@ public:
 
    FiniteElementCollection *OwnFEC() { return fec_owned; }
 
-   /// Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying #fes
-   int VectorDim() const;
+   /** @brief Shortcut for calling FiniteElementSpace::GetVectorDim() on the
+       underlying #fes */
+   int VectorDim() const { return fes->GetVectorDim(); }
 
-   /// Shortcut for calling FiniteElementSpace::GetCurlDim() on the underlying #fes
-   int CurlDim() const;
+   /** @brief Shortcut for calling FiniteElementSpace::GetCurlDim() on the
+       underlying #fes */
+   int CurlDim() const { return fes->GetCurlDim(); }
 
    /// Read only access to the (optional) internal true-dof Vector.
    const Vector &GetTrueVector() const
@@ -513,10 +524,17 @@ public:
        but using an array of scalar coefficients for each component. */
    void ProjectCoefficient(Coefficient *coeff[]);
 
+   /** @brief Project a discontinuous coefficient as a grid function on
+       a continuous finite element space. The values in shared dofs are
+       determined from the element with maximal attribute. */
+   virtual void ProjectDiscCoefficient(Coefficient &coeff)
+   { ProjectDiscCoefficient(&coeff); }
+
    /** @brief Project a discontinuous vector coefficient as a grid function on
        a continuous finite element space. The values in shared dofs are
        determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff);
+   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
+   { ProjectDiscCoefficient(&coeff); }
 
    enum AvgType {ARITHMETIC, HARMONIC};
    /** @brief Projects a discontinuous coefficient so that the values in shared
@@ -532,6 +550,9 @@ public:
    std::unique_ptr<GridFunction> ProlongateToMaxOrder() const;
 
 protected:
+   void ProjectBdrCoefficientNormal(Coefficient *coeff, VectorCoefficient *vcoeff,
+                                    const Array<int> &attr);
+
    /** @brief Accumulates (depending on @a type) the values of @a coeff at all
        shared vdofs and counts in how many zones each vdof appears. */
    void AccumulateAndCountZones(Coefficient &coeff, AvgType type,
@@ -656,15 +677,26 @@ public:
    virtual void ProjectBdrCoefficient(Coefficient *coeff[],
                                       const Array<int> &attr);
 
-   /** Project the normal component of the given VectorCoefficient on
-       the boundary. Only boundary attributes that are marked in
-       'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */
+   /** @brief Project the normal component of the given VectorCoefficient on
+       the boundary. */
+   /** Only boundary attributes that are marked in @a bdr_attr are
+       projected. Assumes RT-type vector finite element GridFunction. */
    void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff,
-                                    const Array<int> &bdr_attr);
+                                    const Array<int> &bdr_attr)
+   { ProjectBdrCoefficientNormal(NULL, &vcoeff, bdr_attr); }
+
+   /** @brief Project the given Coefficient in the normal direction on the
+       boundary. */
+   /** Only boundary attributes that are marked in @a bdr_attr are projected.
+       Assumes RT-type vector finite element GridFunction. */
+   void ProjectBdrCoefficientNormal(Coefficient &coeff,
+                                    const Array<int> &bdr_attr)
+   { ProjectBdrCoefficientNormal(&coeff, NULL, bdr_attr); }
 
    /** @brief Project the tangential components of the given VectorCoefficient
-       on the boundary. Only boundary attributes that are marked in @a bdr_attr
-       are projected. Assumes ND-type VectorFE GridFunction. */
+       on the boundary. */
+   /** Only boundary attributes that are marked in @a bdr_attr
+       are projected. Assumes ND-type vector finite element GridFunction. */
    virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
                                              const Array<int> &bdr_attr);
 
@@ -1914,7 +1946,7 @@ real_t ComputeElementLpDistance(real_t p, int i,
                                 GridFunction& gf1, GridFunction& gf2);
 
 
-/// Class used for extruding scalar GridFunctions
+/// Class used for extruding a scalar coefficient
 class ExtrudeCoefficient : public Coefficient
 {
 private:
@@ -1922,13 +1954,53 @@ private:
    Mesh *mesh_in;
    Coefficient &sol_in;
 public:
+   /// Constructs an instance of VectorExtrudeCoefficient
+   /**
+    * @param m      1D mesh
+    * @param s      1D vector coefficient
+    * @param n_     number of transverse elements of the extruded mesh
+    */
    ExtrudeCoefficient(Mesh *m, Coefficient &s, int n_)
-      : n(n_), mesh_in(m), sol_in(s) { }
+      : n(n_), mesh_in(m), sol_in(s)
+   { MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
+
    real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
+
    virtual ~ExtrudeCoefficient() { }
 };
 
-/// Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
+/// Class used for extruding a vector coefficient
+class VectorExtrudeCoefficient : public VectorCoefficient
+{
+private:
+   int n;
+   Mesh *mesh_in;
+   VectorCoefficient &sol_in;
+public:
+   /// Constructs an instance of VectorExtrudeCoefficient
+   /**
+    * @param m      1D mesh
+    * @param s      1D vector coefficient
+    * @param n_     number of transverse elements of the extruded mesh
+    */
+   VectorExtrudeCoefficient(Mesh *m, VectorCoefficient &s, int n_)
+      : VectorCoefficient(s.GetVDim()), n(n_), mesh_in(m), sol_in(s)
+   { MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
+
+   void Eval(Vector &v, ElementTransformation &T,
+             const IntegrationPoint &ip) override;
+   using VectorCoefficient::Eval;
+
+   virtual ~VectorExtrudeCoefficient() { }
+};
+
+/// Extrude a 1D GridFunction, after extruding the mesh with Extrude1D()
+/**
+ * @param mesh      1D mesh
+ * @param mesh2d    extruded mesh
+ * @param sol       grid function
+ * @param ny        number of transverse elements of the extruded mesh
+ */
 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                     GridFunction *sol, const int ny);
 

fem/gslib.cpp

@@ -106,9 +106,7 @@ FindPointsGSLIB::FindPointsGSLIB()
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   int initialized = 0;
-   MPI_Initialized(&initialized);
-   if (!initialized) { MPI_Init(NULL, NULL); }
+   if (!Mpi::IsInitialized()) { Mpi::Init(); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else
@@ -490,7 +488,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
    }
    DEV.find_device = true;
 
-   const int id = gsl_comm->id, np = gsl_comm->np;
+   const unsigned int id = gsl_comm->id, np = gsl_comm->np;
 
    gsl_mfem_ref.SetSize(points_cnt * dim);
    gsl_mfem_elem.SetSize(points_cnt);
@@ -652,7 +650,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
          {
             const int pp = hash_offset[i];
             /* don't send back to where it just came from */
-            if (pp == p->proc)
+            if (static_cast<unsigned>(pp) == p->proc)
             {
                continue;
             }
@@ -1068,7 +1066,7 @@ void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
          sarray_transfer(struct evalOutPt_t, &outpt, proc, 1, cr);
 
          opt = (evalOutPt_t *)outpt.ptr;
-         for (int index = 0; index < outpt.n; index++)
+         for (size_t index = 0; index < outpt.n; index++)
          {
             int idx = ordering == Ordering::byNODES ?
                       opt->index + i*points_cnt :
@@ -1413,7 +1411,7 @@ void FindPointsGSLIB::SetupSplitMeshesAndIntegrationRules(const int order)
 {
    MFEM_VERIFY(mesh, "Setup FindPointsGSLIB with mesh first.");
    const int dof1D = order+1;
-   const int dim = mesh->Dimension();
+   dim = mesh->Dimension();
 
    SetupSplitMeshes();
    if (dim == 2)
@@ -2254,7 +2252,8 @@ void FindPointsGSLIB::DistributeInterpolatedValues(const Vector &int_vals,
       sarray_transfer(struct out_pt, outpt, proc, 1, cr);
 
       // Store received data
-      MFEM_VERIFY(outpt->n == points_cnt, "Incompatible size. Number of points "
+      MFEM_VERIFY(outpt->n == static_cast<size_t>(points_cnt),
+                  "Incompatible size. Number of points "
                   "received does not match the number of points originally "
                   "found using FindPoints.");
 
@@ -2623,9 +2622,7 @@ GSOPGSLIB::GSOPGSLIB(Array<long long> &ids)
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   int initialized;
-   MPI_Initialized(&initialized);
-   if (!initialized) { MPI_Init(NULL, NULL); }
+   if (!Mpi::IsInitialized()) { Mpi::Init(); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else

fem/gslib.hpp

@@ -202,13 +202,19 @@ protected:
                             const int dof1dsol, const int ordering);
 
 public:
+   /// Serial constructor
    FindPointsGSLIB();
+
+   /// Serial constructor + setup with given Mesh (see \ref Setup)
    FindPointsGSLIB(Mesh &mesh_in, const double bb_t = 0.1,
                    const double newt_tol = 1.0e-12,
                    const int npt_max = 256);
 
 #ifdef MFEM_USE_MPI
+   /// Constructor for ParMesh
    FindPointsGSLIB(MPI_Comm comm_);
+
+   /// Constructor + setup with given ParMesh (see \ref Setup)
    FindPointsGSLIB(ParMesh &mesh_in, const double bb_t = 0.1,
                    const double newt_tol = 1.0e-12,
                    const int npt_max = 256);

fem/gslib/findpts_local_2.cpp

@@ -254,7 +254,7 @@ get_edge(const double *elx[2], const double *wtend, int ei,
       edge.dxdn[d] = workspace + (2 + d) * pN; //dxdn and dydn at DOFs along edge
    }
 
-   if (side_init != (1u << ei))
+   if (static_cast<unsigned>(side_init) != (1u << ei))
    {
 #define ELX(d, j, k) elx[d][j + k * pN] // assumes lexicographic ordering
       for (int d = 0; d < 2; ++d)

fem/gslib/findpts_local_3.cpp

@@ -294,7 +294,7 @@ get_face(const double *elx[3], const double *wtend, int fi, double *workspace,
       face.dxdn[d] = workspace+(3+d)*p_Nfr;
    }
 
-   if (side_init != (1u << fi))
+   if (static_cast<unsigned>(side_init) != (1u << fi))
    {
       const int e_stride[3] = {1, pN, pN*pN};
 #define ELX(d, j, k, l) elx[d][j*e_stride[d1]+k*e_stride[d2]+l*e_stride[dn]]
@@ -342,7 +342,7 @@ get_edge(const double *elx[3], const double *wtend, int ei, double *workspace,
 
    if (jidx >= 3*pN) { return edge; }
 
-   if (side_init != (64u << ei))
+   if (static_cast<unsigned>(side_init) != (64u << ei))
    {
       const int e_stride[3] = {1, pN, pN*pN};
 #define ELX(d, j, k, l) elx[d][j*e_stride[de]+k*e_stride[dn1]+l*e_stride[dn2]]

fem/integ/bilininteg_hdiv_ea.cpp

@@ -197,15 +197,21 @@ static void EAHdivAssemble3D(const int NE,
       // Assemble (one row per thread)
       MFEM_FOREACH_THREAD(idx_i, x, NDOF)
       {
+         // NOTE: due to an llvm backend bug, usage of the modulus operator
+         // has been removed from this foreach section.
          const int ic = idx_i / NDOF_C;
-         const int idx_ii = idx_i % NDOF_C;
+         const int idx_ii = idx_i - ic * NDOF_C; // idx_i % NDOF_C
 
          const int nx_i = (ic == 0) ? D1D : D1D-1;
          const int ny_i = (ic == 1) ? D1D : D1D-1;
 
-         const int ix = idx_ii % nx_i;
-         const int iy = (idx_ii / nx_i) % ny_i;
-         const int iz = (idx_ii / nx_i) / ny_i;
+         const int qx_i = idx_ii / nx_i;
+         const int ix = idx_ii - qx_i * nx_i; // idx_ii % nx_i
+
+         const int qy_i = qx_i / ny_i;
+         const int iy = qx_i - qy_i * ny_i; // (idx_ii / nx_i) % ny_i
+
+         const int iz = qy_i; // (idx_ii / nx_i) / ny_i
 
          const real_t (&Bi1)[MQ1][MD1] = (ic == 0) ? r_Bc : r_Bo;
          const real_t (&Bi2)[MQ1][MD1] = (ic == 1) ? r_Bc : r_Bo;
@@ -214,14 +220,18 @@ static void EAHdivAssemble3D(const int NE,
          for (int idx_j = 0; idx_j < NDOF; ++idx_j)
          {
             const int jc = idx_j / NDOF_C;
-            const int idx_jj = idx_j % NDOF_C;
+            const int idx_jj = idx_j - jc * NDOF_C; // idx_j % NDOF_C
 
             const int nx_j = (jc == 0) ? D1D : D1D-1;
             const int ny_j = (jc == 1) ? D1D : D1D-1;
 
-            const int jx = idx_jj % nx_j;
-            const int jy = (idx_jj / nx_j) % ny_j;
-            const int jz = (idx_jj / nx_j) / ny_j;
+            const int qx_j = idx_jj / nx_j;
+            const int jx = idx_jj - qx_j * nx_j; // idx_jj % nx_j
+
+            const int qy_j = qx_j / ny_j;
+            const int jy = qx_j - qy_j * ny_j; // (idx_jj / nx_j) % ny_j
+
+            const int jz = qy_j; // (idx_jj / nx_j) / ny_j
 
             const real_t (&Bj1)[MQ1][MD1] = (jc == 0) ? r_Bc : r_Bo;
             const real_t (&Bj2)[MQ1][MD1] = (jc == 1) ? r_Bc : r_Bo;

fem/integ/bilininteg_mass_kernels.hpp

@@ -181,6 +181,12 @@ constexpr int NBZ(int D1D)
 {
    return ipow(2, D(D1D) >= 0 ? D(D1D) : 0);
 }
+constexpr int NBZ3D(int MDQ)
+{
+   return MDQ > 0 ? std::min<int>(
+             (128 + MDQ * MDQ * MDQ - 1) / (MDQ * MDQ * MDQ), 64)
+          : 1;
+}
 }
 
 // Shared memory PA Mass Diagonal 2D kernel
@@ -804,19 +810,23 @@ void PAMassApply3D_Element(const int e,
    }
 }
 
-template<int T_D1D, int T_Q1D, bool ACCUMULATE = true>
-MFEM_HOST_DEVICE inline
-void SmemPAMassApply3D_Element(const int e,
-                               const int NE,
-                               const real_t *b_,
-                               const real_t *d_,
-                               const real_t *x_,
-                               real_t *y_,
-                               const int d1d = 0,
-                               const int q1d = 0)
+template <int T_D1D, int T_Q1D, int TBATCH, bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline void
+SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
+                          const real_t *d_, const real_t *x_, real_t *y_,
+                          int d1d = 0, int q1d = 0)
 {
-   constexpr int D1D = T_D1D ? T_D1D : d1d;
-   constexpr int Q1D = T_Q1D ? T_Q1D : q1d;
+   static_assert(TBATCH > 0, "TBATCH must be positive");
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+   constexpr int tbatch = TBATCH;
+   const int tidz = MFEM_THREAD_ID(z);
+#else
+   // host always batch size 1
+   constexpr int tbatch = 1;
+   constexpr int tidz = 0;
+#endif
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
    constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
    constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
    constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
@@ -829,33 +839,37 @@ void SmemPAMassApply3D_Element(const int e,
    MFEM_SHARED real_t sDQ[MQ1*MD1];
    real_t (*B)[MD1] = (real_t (*)[MD1]) sDQ;
    real_t (*Bt)[MQ1] = (real_t (*)[MQ1]) sDQ;
-   MFEM_SHARED real_t sm0[MDQ*MDQ*MDQ];
-   MFEM_SHARED real_t sm1[MDQ*MDQ*MDQ];
-   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) sm0;
-   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) sm1;
-   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm0;
-   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm1;
-   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) sm0;
-   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) sm1;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_SHARED real_t sm0[tbatch][MDQ*MDQ*MDQ];
+   MFEM_SHARED real_t sm1[tbatch][MDQ*MDQ*MDQ];
+   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) (sm0+tidz);
+   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) (sm1+tidz);
+   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm0+tidz);
+   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm1+tidz);
+   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) (sm0+tidz);
+   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) (sm1+tidz);
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          MFEM_UNROLL(MD1)
          for (int dz = 0; dz < D1D; ++dz)
          {
-            X[dz][dy][dx] = x(dx,dy,dz,e);
+            X[dz][dy][dx] = x(dx, dy, dz, e);
          }
       }
-      MFEM_FOREACH_THREAD(dx,x,Q1D)
+      MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
+   }
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(dy, y, D1D)
       {
-         B[dx][dy] = b(dx,dy);
+         MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -880,9 +894,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -907,9 +921,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -929,22 +943,22 @@ void SmemPAMassApply3D_Element(const int e,
          MFEM_UNROLL(MQ1)
          for (int qz = 0; qz < Q1D; qz++)
          {
-            QQQ[qz][qy][qx] = u[qz] * d(qx,qy,qz,e);
+            QQQ[qz][qy][qx] = u[qz] * d(qx, qy, qz, e);
          }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(di,y,D1D)
+   if (tidz == 0)
    {
-      MFEM_FOREACH_THREAD(q,x,Q1D)
+      MFEM_FOREACH_THREAD(di, y, D1D)
       {
-         Bt[di][q] = b(q,di);
+         MFEM_FOREACH_THREAD(q, x, Q1D) { Bt[di][q] = b(q, di); }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -969,9 +983,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -996,9 +1010,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -1020,11 +1034,11 @@ void SmemPAMassApply3D_Element(const int e,
          {
             if (ACCUMULATE)
             {
-               y(dx,dy,dz,e) += u[dz];
+               y(dx, dy, dz, e) += u[dz];
             }
             else
             {
-               y(dx,dy,dz,e) = u[dz];
+               y(dx, dy, dz, e) = u[dz];
             }
          }
       }
@@ -1115,8 +1129,8 @@ inline void PAMassApply3D(const int NE,
    });
 }
 
-// Shared memory PA Mass Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
+// Shared memory PA Mass Apply 3D kernel
+template<int T_D1D = 0, int T_Q1D = 0, int TBATCH=1>
 inline void SmemPAMassApply3D(const int NE,
                               const Array<real_t> &b_,
                               const Array<real_t> &bt_,
@@ -1126,6 +1140,9 @@ inline void SmemPAMassApply3D(const int NE,
                               const int d1d = 0,
                               const int q1d = 0)
 {
+   static_assert(T_D1D > 0, "T_D1D must be positive");
+   static_assert(T_Q1D > 0, "T_Q1D must be positive");
+   static_assert(TBATCH > 0, "TBATCH must be positive");
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1137,9 +1154,11 @@ inline void SmemPAMassApply3D(const int NE,
    const auto d = d_.Read();
    const auto x = x_.Read();
    auto y = y_.ReadWrite();
-   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall_2D_batch<T_Q1D * T_Q1D * TBATCH>(NE, Q1D, Q1D, TBATCH,
+                                                 [=] MFEM_HOST_DEVICE(int e)
    {
-      internal::SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
+      internal::SmemPAMassApply3D_Element<T_D1D, T_Q1D, TBATCH>(e, NE, b, d, x,
+                                                                y, d1d, q1d);
    });
 }
 
@@ -1394,7 +1413,16 @@ ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
 {
    if constexpr (DIM == 1) { return internal::PAMassApply1D; }
    else if constexpr (DIM == 2) { return internal::SmemPAMassApply2D<T_D1D,T_Q1D>; }
-   else if constexpr (DIM == 3) { return internal::SmemPAMassApply3D<T_D1D, T_Q1D>; }
+   else if constexpr (DIM == 3)
+   {
+      constexpr int MDQ = T_D1D >= T_Q1D ? T_D1D : T_Q1D;
+      // max 64 threads in z limit in cuda and hip
+      if constexpr (MDQ > 0)
+      {
+         return internal::SmemPAMassApply3D<T_D1D, T_Q1D,
+                internal::mass::NBZ3D(MDQ)>;
+      }
+   }
    MFEM_ABORT("");
 }
 

fem/intrules.hpp

@@ -43,56 +43,52 @@ public:
       index = i;
    }
 
+   void Set3w(const real_t x1, const real_t x2, const real_t x3, const real_t w)
+   { x = x1; y = x2; z = x3; weight = w; }
+   void Set2w(const real_t x1, const real_t x2, const real_t w)
+   { x = x1; y = x2; weight = w; }
+   void Set1w(const real_t x1, const real_t w)
+   { x = x1; weight = w; }
+
+   void Set3w(const real_t *p) { Set3w(p[0], p[1], p[2], p[3]); }
+   void Set2w(const real_t *p) { Set2w(p[0], p[1], p[2]); }
+   void Set1w(const real_t *p) { Set1w(p[0], p[1]); }
+
+   void Set3(const real_t x1, const real_t x2, const real_t x3)
+   { x = x1; y = x2; z = x3; }
+   void Set2(const real_t x1, const real_t x2)
+   { x = x1; y = x2; }
+   void Set1(const real_t x1)
+   { x = x1; }
+
+   void Set3(const real_t *p) { Set3(p[0], p[1], p[2]); }
+   void Set2(const real_t *p) { Set2(p[0], p[1]); }
+   void Set1(const real_t *p) { Set1(p[0]); }
+
+   void Set(const real_t x1, const real_t x2, const real_t x3, const real_t w)
+   { Set3w(x1, x2, x3, w); }
+
    void Set(const real_t *p, const int dim)
    {
       MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
-      x = p[0];
-      if (dim > 1)
+      switch (dim)
       {
-         y = p[1];
-         if (dim > 2)
-         {
-            z = p[2];
-         }
+         case 3: Set3(p); break;
+         case 2: Set2(p); break;
+         case 1: Set1(p); break;
       }
    }
 
    void Get(real_t *p, const int dim) const
    {
       MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
-      p[0] = x;
-      if (dim > 1)
+      switch (dim)
       {
-         p[1] = y;
-         if (dim > 2)
-         {
-            p[2] = z;
-         }
+         case 3: p[2] = z;
+         case 2: p[1] = y;
+         case 1: p[0] = x;
       }
    }
-
-   void Set(const real_t x1, const real_t x2, const real_t x3, const real_t w)
-   { x = x1; y = x2; z = x3; weight = w; }
-
-   void Set3w(const real_t *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; }
-
-   void Set3(const real_t x1, const real_t x2, const real_t x3)
-   { x = x1; y = x2; z = x3; }
-
-   void Set3(const real_t *p) { x = p[0]; y = p[1]; z = p[2]; }
-
-   void Set2w(const real_t x1, const real_t x2, const real_t w)
-   { x = x1; y = x2; weight = w; }
-
-   void Set2w(const real_t *p) { x = p[0]; y = p[1]; weight = p[2]; }
-
-   void Set2(const real_t x1, const real_t x2) { x = x1; y = x2; }
-
-   void Set2(const real_t *p) { x = p[0]; y = p[1]; }
-
-   void Set1w(const real_t x1, const real_t w) { x = x1; weight = w; }
-
-   void Set1w(const real_t *p) { x = p[0]; weight = p[1]; }
 };
 
 /// Class for an integration rule - an Array of IntegrationPoint.
@@ -125,18 +121,6 @@ private:
    void AddTriPoints3b(const int off, const real_t b, const real_t weight)
    { AddTriPoints3(off, (1. - b)/2., b, weight); }
 
-   void AddTriPoints3R(const int off, const real_t a, const real_t b,
-                       const real_t c, const real_t weight)
-   {
-      IntPoint(off + 0).Set2w(a, b, weight);
-      IntPoint(off + 1).Set2w(c, a, weight);
-      IntPoint(off + 2).Set2w(b, c, weight);
-   }
-
-   void AddTriPoints3R(const int off, const real_t a, const real_t b,
-                       const real_t weight)
-   { AddTriPoints3R(off, a, b, 1. - a - b, weight); }
-
    void AddTriPoints6(const int off, const real_t a, const real_t b,
                       const real_t c, const real_t weight)
    {
@@ -183,14 +167,6 @@ private:
       AddTetPoints3(off + 1, a, 1. - 3.*a, weight);
    }
 
-   // given b, add the permutations of (a,a,a,b), where 3*a + b = 1
-   void AddTetPoints4b(const int off, const real_t b, const real_t weight)
-   {
-      const real_t a = (1. - b)/3.;
-      IntPoint(off).Set(a, a, a, weight);
-      AddTetPoints3(off + 1, a, b, weight);
-   }
-
    // add the permutations of (a,a,b,b), 2*(a + b) = 1
    void AddTetPoints6(const int off, const real_t a, const real_t weight)
    {
@@ -209,14 +185,37 @@ private:
       AddTetPoints6(off + 6, a, bc, cb, weight);
    }
 
-   // given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
-   void AddTetPoints12bc(const int off, const real_t b, const real_t c,
-                         const real_t weight)
+   // add all 24 permutations of (a,b,c,d) where a+b+c+d = 1, all distinct
+   void AddTetPoints24(const int off, const real_t a, const real_t b,
+                       const real_t c, const real_t weight)
    {
-      const real_t a = (1. - b - c)/2.;
-      AddTetPoints3(off,     a, b, weight);
-      AddTetPoints3(off + 3, a, c, weight);
-      AddTetPoints6(off + 6, a, b, c, weight);
+      const real_t d = 1. - a - b - c;
+      // all 24 permutations of 4 distinct barycentric coordinates
+      // permuting which coordinate goes to x, y, z (4th is 1-x-y-z)
+      IntPoint(off +  0).Set(a, b, c, weight);
+      IntPoint(off +  1).Set(a, b, d, weight);
+      IntPoint(off +  2).Set(a, c, b, weight);
+      IntPoint(off +  3).Set(a, c, d, weight);
+      IntPoint(off +  4).Set(a, d, b, weight);
+      IntPoint(off +  5).Set(a, d, c, weight);
+      IntPoint(off +  6).Set(b, a, c, weight);
+      IntPoint(off +  7).Set(b, a, d, weight);
+      IntPoint(off +  8).Set(b, c, a, weight);
+      IntPoint(off +  9).Set(b, c, d, weight);
+      IntPoint(off + 10).Set(b, d, a, weight);
+      IntPoint(off + 11).Set(b, d, c, weight);
+      IntPoint(off + 12).Set(c, a, b, weight);
+      IntPoint(off + 13).Set(c, a, d, weight);
+      IntPoint(off + 14).Set(c, b, a, weight);
+      IntPoint(off + 15).Set(c, b, d, weight);
+      IntPoint(off + 16).Set(c, d, a, weight);
+      IntPoint(off + 17).Set(c, d, b, weight);
+      IntPoint(off + 18).Set(d, a, b, weight);
+      IntPoint(off + 19).Set(d, a, c, weight);
+      IntPoint(off + 20).Set(d, b, a, weight);
+      IntPoint(off + 21).Set(d, b, c, weight);
+      IntPoint(off + 22).Set(d, c, a, weight);
+      IntPoint(off + 23).Set(d, c, b, weight);
    }
 
 public:

fem/linearform.cpp

@@ -297,7 +297,8 @@ void LinearForm::Assemble()
          tr = mesh->GetBdrFaceTransformations(i);
          if (tr != NULL)
          {
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
+            mfem::DofTransformation doftrans;
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
             for (int k = 0; k < boundary_face_integs.Size(); k++)
             {
                if (boundary_face_integs_marker[k] &&
@@ -307,6 +308,7 @@ void LinearForm::Assemble()
                boundary_face_integs[k]->
                AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
                                       *tr, elemvect);
+               doftrans.TransformDual(elemvect);
                AddElementVector (vdofs, elemvect);
             }
          }

fem/lininteg.hpp

@@ -164,8 +164,8 @@ private:
 
 public:
    /// Constructs the domain integrator $ (Q, \nabla v) $
-   DomainLFGradIntegrator(VectorCoefficient &QF)
-      : DeltaLFIntegrator(QF), Q(QF) { }
+   DomainLFGradIntegrator(VectorCoefficient &QF, const IntegrationRule *ir = NULL)
+      : DeltaLFIntegrator(QF, ir), Q(QF) { }
 
    bool SupportsDevice() const override { return true; }
 

fem/lor/lor.cpp

@@ -158,15 +158,16 @@ void LORBase::ConstructLocalDofPermutation(Array<int> &perm_) const
                int i;
                i = dofmap_lor[off_lor + i1 + i2*2];
                int s1 = i < 0 ? -1 : 1;
-               int idof_lor = vdof_lor[absdof(i)];
+               int idof_lor = vdof_lor[UnsignIndex(i)];
                i = dofmap_ho[off_ho + i1*n1 + i2*n2];
                int s2 = i < 0 ? -1 : 1;
-               int idof_ho = vdof_ho[absdof(i)];
+               int idof_ho = vdof_ho[UnsignIndex(i)];
                int s3 = idof_lor < 0 ? -1 : 1;
                int s4 = idof_ho < 0 ? -1 : 1;
                int s = s1*s2*s3*s4;
-               i = absdof(idof_ho);
-               perm_[absdof(idof_lor)] = s < 0 ? -1-absdof(i) : absdof(i);
+               i = UnsignIndex(idof_ho);
+               perm_[UnsignIndex(idof_lor)] = s < 0 ? -1-UnsignIndex(i) :
+                                              UnsignIndex(i);
             }
          }
       };
@@ -232,7 +233,7 @@ void LORBase::ConstructDofPermutation() const
          int j = l_perm[i];
          int s = j < 0 ? -1 : 1;
          int t_i = pfes_lor->GetLocalTDofNumber(i);
-         int t_j = pfes_ho->GetLocalTDofNumber(absdof(j));
+         int t_j = pfes_ho->GetLocalTDofNumber(UnsignIndex(j));
          // Either t_i and t_j both -1, or both non-negative
          if ((t_i < 0 && t_j >=0) || (t_j < 0 && t_i >= 0))
          {

fem/lor/lor.hpp

@@ -57,8 +57,6 @@ private:
    /// values (after temporarily changing them for LOR assembly).
    void ResetIntegrationRules(GetIntegratorsFn get_integrators);
 
-   static inline int absdof(int i) { return i < 0 ? -1-i : i; }
-
 protected:
    enum FESpaceType { H1, ND, RT, L2, INVALID };
 

fem/pfespace.cpp

@@ -424,7 +424,7 @@ void ParFiniteElementSpace::GetGroupComm(
             {
                if (ind[l] < 0)
                {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                   if (g_ldof_sign)
                   {
                      (*g_ldof_sign)[dofs[l]] = -1;
@@ -462,7 +462,7 @@ void ParFiniteElementSpace::GetGroupComm(
             {
                if (ind[l] < 0)
                {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                   if (g_ldof_sign)
                   {
                      (*g_ldof_sign)[dofs[l]] = -1;
@@ -500,7 +500,7 @@ void ParFiniteElementSpace::GetGroupComm(
             {
                if (ind[l] < 0)
                {
-                  dofs[l] = m + (-1-ind[l]);
+                  dofs[l] = m + FlipIndexSign(ind[l]);
                   if (g_ldof_sign)
                   {
                      (*g_ldof_sign)[dofs[l]] = -1;
@@ -538,16 +538,16 @@ void ParFiniteElementSpace::ApplyLDofSigns(Array<int> &dofs) const
    {
       if (dofs[i] < 0)
       {
-         if (ldof_sign[-1-dofs[i]] < 0)
+         if (ldof_sign[FlipIndexSign(dofs[i])] < 0)
          {
-            dofs[i] = -1-dofs[i];
+            dofs[i] = FlipIndexSign(dofs[i]);
          }
       }
       else
       {
          if (ldof_sign[dofs[i]] < 0)
          {
-            dofs[i] = -1-dofs[i];
+            dofs[i] = FlipIndexSign(dofs[i]);
          }
       }
    }
@@ -699,7 +699,8 @@ void ParFiniteElementSpace::GetSharedEdgeDofs(
       for (int i = 0; i < dofs.Size(); i++)
       {
          const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = di >= 0 ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
       }
    }
 }
@@ -723,7 +724,8 @@ void ParFiniteElementSpace::GetSharedTriangleDofs(
       for (int i = 0; i < dofs.Size(); i++)
       {
          const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = di >= 0 ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
       }
    }
 }
@@ -747,7 +749,8 @@ void ParFiniteElementSpace::GetSharedQuadrilateralDofs(
       for (int i = 0; i < dofs.Size(); i++)
       {
          const int di = dofs[i];
-         dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
+         dofs[i] = (di >= 0) ? rdofs[di] :
+                   FlipIndexSign(rdofs[FlipIndexSign(di)]);
       }
    }
 }
@@ -1487,7 +1490,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
          GetElementVDofs(my_elems[i], ldofs);
          for (int j = 0; j < ldofs.Size(); j++)
          {
-            int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
+            int ldof = UnsignIndex(ldofs[j]);
 
             if (ldof_marker[ldof] != fn)
             {
@@ -1548,7 +1551,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
          GetElementVDofs(my_elems[i], ldofs);
          for (int j = 0; j < ldofs.Size(); j++)
          {
-            int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
+            int ldof = UnsignIndex(ldofs[j]);
 
             if (ldof_marker[ldof] != fn)
             {
@@ -1573,14 +1576,15 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
 
       for (int i = 0; i < num_ldofs; i++)
       {
-         int ldof = (ldofs_fn[i] >= 0 ? ldofs_fn[i] : -1-ldofs_fn[i]);
+         int ldof = UnsignIndex(ldofs_fn[i]);
          ldof_marker[ldof] = i;
       }
 
       for ( ; j < j_end; j++)
       {
-         int ldof = (send_J[j] >= 0 ? send_J[j] : -1-send_J[j]);
-         send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] : -1-ldof_marker[ldof]);
+         const int ldof = UnsignIndex(send_J[j]);
+         send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] :
+                      FlipIndexSign(ldof_marker[ldof]));
       }
    }
 
@@ -1672,12 +1676,7 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
    {
       for (int j_end = face_nbr_ldof.GetI()[fn+1]; j < j_end; j++)
       {
-         int ldof = face_nbr_ldof.GetJ()[j];
-         if (ldof < 0)
-         {
-            ldof = -1-ldof;
-         }
-
+         const int ldof = UnsignIndex(face_nbr_ldof.GetJ()[j]);
          face_nbr_glob_dof_map[j] = dof_face_nbr_offsets[fn] + ldof;
       }
    }
@@ -1721,7 +1720,7 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
    MFEM_ASSERT(Nonconforming() && i >= pmesh->GetNumFaces(), "");
    int el1, el2, inf1, inf2;
    pmesh->GetFaceElements(i, &el1, &el2);
-   el2 = -1 - el2;
+   el2 = FlipIndexSign(el2);
    pmesh->GetFaceInfos(i, &inf1, &inf2);
    MFEM_ASSERT(0 <= el2 && el2 < face_nbr_element_dof.Size(), "");
    const int nd = face_nbr_element_dof.RowSize(el2);
@@ -1737,7 +1736,8 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
    for (int j = 0; j < vdofs.Size(); j++)
    {
       const int ldof = vdofs[j];
-      vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] : -1-vol_vdofs[-1-ldof];
+      vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] :
+                 FlipIndexSign(vol_vdofs[FlipIndexSign(ldof)]);
    }
 }
 
@@ -2061,8 +2061,8 @@ void ParFiniteElementSpace::GetGhostFaceDofs(const MeshId &face_id,
 
          for (int j = 0; j < ne; j++)
          {
-            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
-                             /*         */ : (-1 - (first + (-1 - ind[j])));
+            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j]) :
+                             FlipIndexSign(first + FlipIndexSign(ind[j]));
          }
       }
       else
@@ -2072,8 +2072,8 @@ void ParFiniteElementSpace::GetGhostFaceDofs(const MeshId &face_id,
          const int *ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[i]);
          for (int j = 0; j < ne; j++)
          {
-            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
-                             /*         */ : (-1 - (first + (-1 - ind[j])));
+            dofs[offset++] = (ind[j] >= 0) ? (first + ind[j]) :
+                             FlipIndexSign(first + FlipIndexSign(ind[j]));
          }
       }
    }
@@ -2866,7 +2866,7 @@ void NeighborRowMessage::Encode(int rank)
 
             if (ind && (edof = ind[edof]) < 0)
             {
-               edof = -1 - edof;
+               edof = FlipIndexSign(edof);
                s = -1;
             }
 
@@ -3067,10 +3067,10 @@ void NeighborRowMessage::Decode(int rank)
 
          // If edof arrived with a negative index, flip it, and the scaling.
          real_t s = (edof < 0) ? -1.0 : 1.0;
-         edof = (edof < 0) ? -1 - edof : edof;
+         edof = UnsignIndex(edof);
          if (ind && (edof = ind[edof]) < 0)
          {
-            edof = -1 - edof;
+            edof = FlipIndexSign(edof);
             s *= -1.0;
          }
 
@@ -3121,10 +3121,10 @@ void NeighborRowMessage::Decode(int rank)
 
             // If edof arrived with a negative index, flip it, and the scaling.
             s = (edof < 0) ? -1.0 : 1.0;
-            edof = (edof < 0) ? -1 - edof : edof;
+            edof = UnsignIndex(edof);
             if (ind && (edof = ind[edof]) < 0)
             {
-               edof = -1 - edof;
+               edof = FlipIndexSign(edof);
                s *= -1.0;
             }
 
@@ -4405,12 +4405,9 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
          {
             for (int j = 0; j < dofs.Size(); j++)
             {
-               int row = DofToVDof(dofs[j], vd);
-               if (row < 0) { row = -1 - row; }
-
-               int col = DofToVDof(old_dofs[j], vd, old_ndofs);
-               if (col < 0) { col = -1 - col; }
-
+               const int row = UnsignIndex(DofToVDof(dofs[j], vd));
+               const int col = UnsignIndex(DofToVDof(old_dofs[j], vd,
+                                                     old_ndofs));
                i_diag[row] = col;
             }
          }
@@ -4435,9 +4432,7 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
       {
          for (int j = 0; j < dofs.Size(); j++)
          {
-            int row = DofToVDof(dofs[j], vd);
-            if (row < 0) { row = -1 - row; }
-
+            const int row = UnsignIndex(DofToVDof(dofs[j], vd));
             if (i_diag[row] == i_diag[row+1]) // diag row empty?
             {
                i_offd[row] = old_dofs[j + vd * dofs.Size()];
@@ -4546,9 +4541,9 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
    {
       const Embedding &emb = dtrans.embeddings[k];
 
-      int fine_rank = old_ranks[k];
-      int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
-                        : old_pncmesh->ElementRank(emb.parent);
+      const int fine_rank = old_ranks[k];
+      const int coarse_rank = (emb.parent < 0) ? FlipIndexSign(emb.parent)
+                              : old_pncmesh->ElementRank(emb.parent);
 
       if (coarse_rank != MyRank && fine_rank == MyRank)
       {
@@ -4636,8 +4631,8 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
             {
                if (!std::isfinite(lR(i, 0))) { continue; }
 
-               int r = DofToVDof(dofs[i], vd);
-               int m = (r >= 0) ? r : (-1 - r);
+               const int r = DofToVDof(dofs[i], vd);
+               const int m = UnsignIndex(r);
 
                if (is_dg || !mark[m])
                {
@@ -4686,8 +4681,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
             {
                if (!std::isfinite(lR(i, 0))) { continue; }
 
-               int r = DofToVDof(dofs[i], vd);
-               int m = (r >= 0) ? r : (-1 - r);
+               const int m = UnsignIndex(DofToVDof(dofs[i], vd));
 
                if (is_dg || !mark[m])
                {

fem/pgridfunc.cpp

@@ -545,6 +545,8 @@ void ParGridFunction::GetElementDofValues(int el, Vector &dof_vals) const
 
 void ParGridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
+   MFEM_VERIFY(VectorDim() == 1,
+               "Cannot project scalar coefficient onto vector ParGridFunction");
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
 
    if (delta_c == NULL)
@@ -715,7 +717,8 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 }
 
 
-void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
+void ParGridFunction::ProjectDiscCoefficient(
+   std::variant<Coefficient*, VectorCoefficient*> coeff)
 {
    // local maximal element attribute for each dof
    Array<int> ldof_attr;
@@ -761,6 +764,9 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 
 void ParGridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
+   MFEM_VERIFY(
+      VectorDim() == 1,
+      "Cannot project scalar coefficient onto a vector ParGridFunction");
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
@@ -786,6 +792,8 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &vcoeff,
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
+   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
+
    // Number of zones that contain a given dof.
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(vcoeff, type, zones_per_vdof);
@@ -858,6 +866,12 @@ void ParGridFunction::ProjectBdrCoefficient(
 #endif
 }
 
+void ParGridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
+                                            const Array<int> &attr)
+{
+   ProjectBdrCoefficient(NULL, &vcoeff, attr);
+}
+
 void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
                                                    const Array<int> &bdr_attr)
 {

fem/pgridfunc.hpp

@@ -63,6 +63,12 @@ protected:
    void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
                               const Array<int> &attr);
 
+   /** @brief Project a discontinuous (vector) coefficient as a grid function on
+       a continuous finite element space. The values in shared dofs are
+       determined from the element with maximal attribute. */
+   virtual void ProjectDiscCoefficient(
+      std::variant<Coefficient*, VectorCoefficient*> coeff) override;
+
 public:
    ParGridFunction() { pfes = NULL; }
 
@@ -268,11 +274,6 @@ public:
                            ProjectType type = ProjectType::DEFAULT) override;
 
    using GridFunction::ProjectDiscCoefficient;
-   /** @brief Project a discontinuous vector coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   void ProjectDiscCoefficient(VectorCoefficient &coeff) override;
-
    void ProjectDiscCoefficient(Coefficient &coeff, AvgType type) override;
 
    void ProjectDiscCoefficient(VectorCoefficient &vcoeff, AvgType type) override;
@@ -280,8 +281,7 @@ public:
    using GridFunction::ProjectBdrCoefficient;
 
    void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
-                              const Array<int> &attr) override
-   { ProjectBdrCoefficient(NULL, &vcoeff, attr); }
+                              const Array<int> &attr) override;
 
    void ProjectBdrCoefficient(Coefficient *coeff[],
                               const Array<int> &attr) override

fem/prestriction.cpp

@@ -321,12 +321,17 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
    const int vd = vdim;
    const bool t = byvdim;
    const int threshold = ndofs;
-   const int nsdofs = pfes.GetFaceNbrVSize();
+   const int nsdofs = pfes.GetFaceNbrVSize() / vd;
    auto d_indices1 = scatter_indices1.Read();
    auto d_indices2 = scatter_indices2.Read();
    auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
-   auto d_x_shared = Reshape(face_nbr_data.Read(),
-                             t?vd:nsdofs, t?nsdofs:vd);
+   const int ne_shared = nsdofs / elem_dofs;
+   const int nedof = elem_dofs;
+   // Note: the shape of face_nbr_data, as determined by
+   // ParFiniteElementSpace::ExchangeFaceNbrData, is (elem_dofs, vdim,
+   // ne_shared), independent of the ordering (byNODES or byVDIM) of the finite
+   // element space.
+   auto d_x_shared = Reshape(face_nbr_data.Read(), elem_dofs, vd, ne_shared);
    auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
    mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
    {
@@ -346,8 +351,9 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
          }
          else if (idx2>=threshold) // shared boundary
          {
-            d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
-                                              t?(idx2-threshold):c);
+            const int e_shared = (idx2 - threshold) / nedof;
+            const int i_shared = (idx2 - threshold) % nedof;
+            d_y(dof, c, 1, face) = d_x_shared(i_shared,c,e_shared);
          }
          else // true boundary
          {

fem/qfunction.hpp

@@ -271,10 +271,7 @@ inline void QuadratureFunction::GetValues(
    const int s_offset = qspace->Offset(idx);
    const int sl_size = qspace->Offset(idx + 1) - s_offset;
    // Make the values matrix memory an alias of the quadrature function memory
-   Memory<real_t> &values_mem = values.GetMemory();
-   values_mem.Delete();
-   values_mem.MakeAlias(GetMemory(), vdim*s_offset, vdim*sl_size);
-   values.SetSize(vdim, sl_size);
+   values.MakeRef(GetMemory(), vdim*s_offset, vdim, sl_size);
 }
 
 inline void QuadratureFunction::GetValues(

fem/restriction.cpp

@@ -844,8 +844,6 @@ void ConformingFaceRestriction::ComputeGatherIndices(
    gather_offsets[0] = 0;
 }
 
-static inline int absdof(int i) { return i < 0 ? -1-i : i; }
-
 void ConformingFaceRestriction::SetFaceDofsScatterIndices(
    const Mesh::FaceInformation &face,
    const int face_index,
@@ -868,9 +866,9 @@ void ConformingFaceRestriction::SetFaceDofsScatterIndices(
    {
       const int lex_volume_dof = face_map[face_dof];
       const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof]; // signed
-      const int volume_dof = absdof(s_volume_dof);
+      const int volume_dof = UnsignIndex(s_volume_dof);
       const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
-      const int global_dof = absdof(s_global_dof);
+      const int global_dof = UnsignIndex(s_global_dof);
       const int restriction_dof = face_dofs*face_index + face_dof;
       scatter_indices[restriction_dof] = s_global_dof;
       ++gather_offsets[global_dof + 1];
@@ -897,10 +895,10 @@ void ConformingFaceRestriction::SetFaceDofsGatherIndices(
    {
       const int lex_volume_dof = face_map[face_dof];
       const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof];
-      const int volume_dof = absdof(s_volume_dof);
+      const int volume_dof = UnsignIndex(s_volume_dof);
       const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
       const int sgn = (s_global_dof >= 0) ? 1 : -1;
-      const int global_dof = absdof(s_global_dof);
+      const int global_dof = UnsignIndex(s_global_dof);
       const int restriction_dof = face_dofs*face_index + face_dof;
       const int s_restriction_dof = (sgn >= 0) ? restriction_dof : -1 -
                                     restriction_dof;
@@ -1400,20 +1398,17 @@ void L2FaceRestriction::PermuteAndSetSharedFaceDofsScatterIndices2(
    const int dim = fes.GetMesh()->Dimension();
    const int dof1d = fes.GetTypicalFE()->GetOrder()+1;
    fes.GetTypicalFE()->GetFaceMap(face_id2, face_map);
-   Array<int> face_nbr_dofs;
-   const ParFiniteElementSpace &pfes =
-      static_cast<const ParFiniteElementSpace&>(this->fes);
-   pfes.GetFaceNbrElementVDofs(elem_index, face_nbr_dofs);
 
    for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
    {
       const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
                                                orientation, dof1d, face_dof_elem1);
       const int volume_dof_elem2 = face_map[face_dof_elem2];
-      const int global_dof_elem2 = face_nbr_dofs[volume_dof_elem2];
+      // Encode the volume DOF index and element index
+      const int global_dof_elem2 = elem_index*elem_dofs + volume_dof_elem2;
       const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
       // Trick to differentiate dof location inter/shared
-      scatter_indices2[restriction_dof_elem2] = ndofs+global_dof_elem2;
+      scatter_indices2[restriction_dof_elem2] = ndofs + global_dof_elem2;
    }
 #endif
 }

fem/tmop.cpp

@@ -4102,8 +4102,11 @@ void TMOP_Integrator::GetSurfaceFittingErrors(const Vector &d_loc,
 #ifdef MFEM_USE_MPI
       // Don't count the overlapping DOFs in parallel.
       // The pfes might be ordered byVDIM, while the loop goes consecutively.
-      const int dof_i = pfes->DofToVDof(i, 0);
-      if (parallel && pfes->GetLocalTDofNumber(dof_i) < 0) { continue; }
+      if (parallel)
+      {
+         const int dof_i = pfes->DofToVDof(i, 0);
+         if (pfes->GetLocalTDofNumber(dof_i) < 0) { continue; }
+      }
 #endif
 
       dof_cnt++;

general/array.hpp

@@ -114,10 +114,22 @@ public:
    Array<T> &operator=(const Array<T> &src) { src.Copy(*this); return *this; }
 
    /// Move assignment operator
+   /** If *this is a non-owning view (e.g., from MakeRef()), the data is copied
+       so that the base is also modified. */
    Array<T> &operator=(Array<T> &&src)
    {
       if (this == &src) { return *this; }
-      Swap(src);  // Swap does not use move assignment!
+      // If *this is a non-owning view (alias), and its capacity is sufficient
+      // to contain src, then copy into *this so that the alias's base memory is
+      // modified.
+      if (!OwnsData() && Capacity() >= src.Size())
+      {
+         *this = src; // Copy assignment.
+      }
+      else
+      {
+         Swap(src); // Swap the pointers only.
+      }
       src.DeleteAll();
       return *this;
    }
@@ -251,6 +263,9 @@ public:
    /// Make this Array a reference to 'master'.
    inline void MakeRef(const Array &master);
 
+   /// Make this Array a reference to the given sub-Memory of @a base.
+   inline void MakeRef(Memory<T> &base, int offset, int size_);
+
    /// Reset the Array to use the given external Memory @a mem and size @a s.
    /** If @a own_mem is false, the Array will not own any of the pointers of
        @a mem.
@@ -1073,6 +1088,14 @@ inline void Array<T>::MakeRef(const Array &master)
    data.MakeAlias(master.GetMemory(), 0, size);
 }
 
+template <class T>
+inline void Array<T>::MakeRef(Memory<T> &base, int offset, int size_)
+{
+   data.Delete();
+   size = size_;
+   data.MakeAlias(base, offset, size_);
+}
+
 template <class T>
 inline void Array<T>::NewMemoryAndSize(
    const Memory<T> &mem, int s, bool own_mem)

general/arrays_by_name.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "array.hpp"
+#include "text.hpp"
 
 #include <iostream>
 #include <map>
@@ -247,7 +248,8 @@ inline void ArraysByName<T>::Print(std::ostream &os, int width) const
    os << data.size() << '\n';
    for (auto const &it : data)
    {
-      os << '"' << it.first << '"' << '\n' << it.second.Size() << '\n';
+      // Note: The method Load() can read any string formatted with std::quoted.
+      os << std::quoted(it.first) << '\n' << it.second.Size() << '\n';
       it.second.Print(os, width > 0 ? width : it.second.Size());
    }
 }
@@ -258,31 +260,36 @@ void ArraysByName<T>::Load(std::istream &in)
    int NumArrays;
    in >> NumArrays;
 
-   std::string ArrayLine, ArrayName;
-   for (int i=0; i < NumArrays; i++)
+   for (int i = 0; i < NumArrays; i++)
    {
       in >> std::ws;
-      getline(in, ArrayLine);
-
-      std::size_t q0 = ArrayLine.find('"');
-      std::size_t q1 = ArrayLine.rfind('"');
-
-      if (q0 != std::string::npos && q1 > q0)
+      // Read the name:
+      // - If the stream 'in' starts with " then parse it with the function
+      //   parse_quoted_string() from text.hpp. In this case, the name can be
+      //   empty. Note: this case allows for reading any string formatted using
+      //   std::quoted, e.g. as in the method Print().
+      // - If the name does not start with " then the name ends with the first
+      //   white space character (and the white space character is not included
+      //   in the name). Since white space characters are skipped before reading
+      //   the name, there will be at least one non-white-space character in the
+      //   name in this case.
+      std::string ArrayName;
+      if (in.peek() == '"')
       {
-         // Locate set name between first and last double quote
-         ArrayName = ArrayLine.substr(q0+1,q1-q0-1);
+         if (parse_quoted_string(ArrayName, in) != 0)
+         {
+            MFEM_ABORT("error parsing input!");
+         }
       }
       else
       {
-         // If no double quotes found locate set name using white space
-         q1 = ArrayLine.find(' ');
-         ArrayName = ArrayLine.substr(0,q1-1);
+         in >> ArrayName;
+         MFEM_VERIFY(in.good(), "error parsing input!");
       }
 
-      // Ignore the remainder of the line which may contain explanatory comments
-      data[ArrayName].Load(in, 0);
+      // Read the array
+      data[ArrayName].Load(in);
    }
-
 }
 
 }

general/device.cpp

@@ -726,16 +726,16 @@ std::string Device::GetUUID(const int device_id)
    MFEM_GPU_CHECK(cudaGetDeviceProperties(&prop, device_id));
    for (int i = 0; i < 16; ++i)
    {
-      res << std::setfill('0') << std::setw(2) << std::hex
-          << static_cast<unsigned>(prop.uuid.bytes[i]);
+      const unsigned b = static_cast<unsigned char>(prop.uuid.bytes[i]);
+      res << std::setfill('0') << std::setw(2) << std::hex << b;
    }
 #elif defined(MFEM_USE_HIP)
    hipUUID uuid;
    MFEM_GPU_CHECK(hipDeviceGetUuid(&uuid, device_id));
    for (int i = 0; i < 16; ++i)
    {
-      res << std::setfill('0') << std::setw(2) << std::hex
-          << static_cast<unsigned>(uuid.bytes[i]);
+      const unsigned b = static_cast<unsigned char>(uuid.bytes[i]);
+      res << std::setfill('0') << std::setw(2) << std::hex << b;
    }
 #endif
    return res.str();

general/forall.hpp

@@ -1090,6 +1090,12 @@ inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
    ForallWrap<2>(true, N, body, X, Y, BZ);
 }
 
+template<int MAX_THREADS_PER_BLOCK, typename lambda>
+inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
+{
+   ForallWrap<2, MAX_THREADS_PER_BLOCK>(true, N, body, X, Y, BZ);
+}
+
 template<typename lambda>
 inline void forall_3D(int N, int X, int Y, int Z, lambda &&body)
 {

general/globals.hpp

@@ -113,6 +113,10 @@ void SetGlobalMPI_Comm(MPI_Comm comm);
 /// to suppress the warning.
 const char* GetEnv(const char* name);
 
+/// Signed indices i -> -1 - i are used as a convention to encode orientation.
+inline MFEM_HOST_DEVICE int FlipIndexSign(int i) { return -1 - i; }
+inline MFEM_HOST_DEVICE int UnsignIndex(int i) { return i < 0 ? -1 - i : i; }
+
 } // namespace mfem
 
 #endif

general/mem_manager.cpp

@@ -759,7 +759,9 @@ private:
    {
       switch (mt)
       {
-         case MT::HOST_DEBUG: return new MmuHostMemorySpace();
+         case MT::HOST_DEBUG:
+            if (GetEnv("MFEM_MMU_STD")) { return new StdHostMemorySpace(); }
+            return new MmuHostMemorySpace();
 #ifdef MFEM_USE_UMPIRE
          case MT::HOST_UMPIRE:
             return new UmpireHostMemorySpace(
@@ -788,7 +790,9 @@ private:
          case MT::DEVICE_UMPIRE: return new NoDeviceMemorySpace();
          case MT::DEVICE_UMPIRE_2: return new NoDeviceMemorySpace();
 #endif
-         case MT::DEVICE_DEBUG: return new MmuDeviceMemorySpace();
+         case MT::DEVICE_DEBUG:
+            if (GetEnv("MFEM_MMU_STD")) { return new StdDeviceMemorySpace(); }
+            return new MmuDeviceMemorySpace();
          case MT::DEVICE:
          {
 #if defined(MFEM_USE_CUDA)

general/text.hpp

@@ -50,6 +50,48 @@ inline void filter_dos(std::string &line)
    }
 }
 
+/** @brief Read a string formatted using std::quoted. Return nonzero on error.
+
+    The stream @a in must begin with @a delim. After clearing @a result and
+    extracting the opening @a delim, characters are extracted from @a in and
+    processed as follows:
+    - if the character is @a delim, return 0;
+    - if the character is different from @a escape, it is appended to @a result;
+    - if the character is @a escape, the next character from @a in is extracted
+      and if it is one of @a delim or @a escape, it is appended to @a result;
+      otherwise, both @a escape and the character after it are appended to
+      @a result; note that the latter case is not possible if the input was
+      formatted with std::quoted with the same @a delim and @a escape
+      characters.
+
+    If the stream @a in does not begin with @a delim, error code 1 is returned.
+    If reading the stream fails, error code 2 is returned. On success, zero is
+    returned and the closing @a delim character is the last character extracted
+    from @a in. */
+inline int parse_quoted_string(std::string &result, std::istream &in,
+                               char delim = '"', char escape = '\\')
+{
+   using tt = std::string::traits_type;  // std::char_traits<char>
+   auto equal = [](tt::int_type c1, tt::char_type c2) -> bool
+   {
+      return tt::eq_int_type(c1, tt::to_int_type(c2));
+   };
+   result.clear();
+   if (!equal(in.peek(), delim)) { return 1; }
+   in.get();  // extract delim
+   for (auto c = in.get(); !equal(c, delim); c = in.get())
+   {
+      if (equal(c, escape))
+      {
+         c = in.get();
+         if (!equal(c, escape) && !equal(c, delim)) { result += escape; }
+      }
+      if (!in) { return 2; }
+      result += tt::to_char_type(c);
+   }
+   return 0;
+}
+
 /// Convert an integer to a 0-padded string with the given number of @a digits
 inline std::string to_padded_string(int i, int digits)
 {

linalg/densemat.hpp

@@ -82,6 +82,15 @@ public:
       }
    }
 
+   /// Make the DenseMatrix to reference the given sub-Memory of @a base.
+   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
+       not delete the @a base Memory. */
+   void MakeRef(Memory<real_t> &base, int offset, int h, int w)
+   {
+      data.MakeRef(base, offset, h*w);
+      height = h; width = w;
+   }
+
    /// Change the data array and the size of the DenseMatrix.
    /** The DenseMatrix does not assume ownership of the data array, i.e. it will
        not delete the data array @a d. */

linalg/hypre.cpp

@@ -317,6 +317,9 @@ void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
 
 Vector * HypreParVector::GlobalVector() const
 {
+   MFEM_VERIFY(size > 0,
+               "GlobalVector method can only be called on vectors wherein each "
+               "process owns one or more entries");
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
    Vector *v = new Vector(hv->data, internal::to_int(hv->size));
    v->MakeDataOwner();

linalg/petsc.cpp

@@ -38,6 +38,13 @@
 #if PETSC_VERSION_LT(3,19,0)
 #define PETSC_SUCCESS 0
 #endif
+#if PETSC_VERSION_LT(3,23,0)
+#define PetscContainerSetCtxDestroy(A,B) PetscContainerSetUserDestroy(A,B)
+typedef PetscErrorCode (PetscCtxDestroyFn)(void**);
+#endif
+#if PETSC_VERSION_LT(3,24,0)
+typedef PetscErrorCode KSPMonitorFn(KSP,PetscInt,PetscReal,void*);
+#endif
 
 #include <fstream>
 #include <iomanip>
@@ -77,13 +84,17 @@ static PetscErrorCode __mfem_mat_shell_apply_transpose(Mat,Vec,Vec);
 static PetscErrorCode __mfem_mat_shell_destroy(Mat);
 static PetscErrorCode __mfem_mat_shell_copy(Mat,Mat,MatStructure);
 #if PETSC_VERSION_LT(3,23,0)
-static PetscErrorCode __mfem_array_container_destroy(void*);
-static PetscErrorCode __mfem_matarray_container_destroy(void *);
-#else
-static PetscErrorCode __mfem_array_container_destroy(void**);
-static PetscErrorCode __mfem_matarray_container_destroy(void**);
+typedef void *PetscCtxRt;
+#elif PETSC_VERSION_LT(3,25,0)
+typedef void **PetscCtxRt;
 #endif
+static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt);
+static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt);
+#if PETSC_VERSION_LT(3,23,0)
 static PetscErrorCode __mfem_monitor_ctx_destroy(void**);
+#else
+static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt);
+#endif
 
 // auxiliary functions
 static PetscErrorCode Convert_Array_IS(MPI_Comm,bool,const mfem::Array<int>*,
@@ -1317,11 +1328,7 @@ BlockDiagonalConstructor(MPI_Comm comm,
 
          ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
          ierr = PetscContainerSetPointer(c,ptrs[i]); CCHKERRQ(comm,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-         ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
          ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
          CCHKERRQ(comm,ierr);
          ierr = PetscObjectCompose((PetscObject)A,names[i],(PetscObject)c);
          CCHKERRQ(comm,ierr);
@@ -1648,11 +1655,7 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
          PetscContainer c;
          ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
          ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-         ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
-#else
          ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
-#endif
          PCHKERRQ(c,ierr);
          ierr = PetscObjectCompose((PetscObject)(*A),"_MatIS_PtAP_l2l",(PetscObject)c);
          PCHKERRQ((*A),ierr);
@@ -1748,11 +1751,7 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
 
             ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); PCHKERRQ(B,ierr);
             ierr = PetscContainerSetPointer(c,ptrs[i]); PCHKERRQ(B,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-            ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
             ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
             PCHKERRQ(B,ierr);
             ierr = PetscObjectCompose((PetscObject)(B),names[i],(PetscObject)c);
             PCHKERRQ(B,ierr);
@@ -2198,11 +2197,7 @@ PetscParMatrix * RAP(PetscParMatrix *Rt, PetscParMatrix *A, PetscParMatrix *P)
          ierr = PetscContainerCreate(PetscObjectComm((PetscObject)B),&c);
          PCHKERRQ(B,ierr);
          ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
-#if PETSC_VERSION_LT(3,23,0)
-         ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
-#else
          ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
-#endif
          PCHKERRQ(c,ierr);
          ierr = PetscObjectCompose((PetscObject)B,"_MatIS_PtAP_l2l",(PetscObject)c);
          PCHKERRQ(B,ierr);
@@ -2485,7 +2480,6 @@ void PetscSolver::SetMaxIter(int max_iter)
 
 void PetscSolver::SetPrintLevel(int plev)
 {
-   typedef PetscErrorCode (*myPetscFunc)(void**);
    PetscViewerAndFormat *vf = NULL;
    PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm(obj));
 
@@ -2498,7 +2492,6 @@ void PetscSolver::SetPrintLevel(int plev)
    {
       // there are many other options, see the function KSPSetFromOptions() in
       // src/ksp/ksp/interface/itcl.c
-      typedef PetscErrorCode (*myMonitor)(KSP,PetscInt,PetscReal,void*);
       KSP ksp = (KSP)obj;
       if (plev >= 0)
       {
@@ -2507,29 +2500,29 @@ void PetscSolver::SetPrintLevel(int plev)
       if (plev == 1)
       {
 #if PETSC_VERSION_LT(3,15,0)
-         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorDefault,vf,
+         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorDefault,vf,
 #else
-         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorResidual,vf,
+         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorResidual,vf,
 #endif
-                              (myPetscFunc)PetscViewerAndFormatDestroy);
+                              (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
          PCHKERRQ(ksp,ierr);
       }
       else if (plev > 1)
       {
          ierr = KSPSetComputeSingularValues(ksp,PETSC_TRUE); PCHKERRQ(ksp,ierr);
-         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorSingularValue,vf,
-                              (myPetscFunc)PetscViewerAndFormatDestroy);
+         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorSingularValue,vf,
+                              (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
          PCHKERRQ(ksp,ierr);
          if (plev > 2)
          {
             ierr = PetscViewerAndFormatCreate(viewer,PETSC_VIEWER_DEFAULT,&vf);
             PCHKERRQ(viewer,ierr);
 #if PETSC_VERSION_LT(3,15,0)
-            ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidualNorm,vf,
+            ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidualNorm,vf,
 #else
-            ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidual,vf,
+            ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidual,vf,
 #endif
-                                 (myPetscFunc)PetscViewerAndFormatDestroy);
+                                 (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
             PCHKERRQ(ksp,ierr);
          }
       }
@@ -2545,7 +2538,7 @@ void PetscSolver::SetPrintLevel(int plev)
       if (plev > 0)
       {
          ierr = SNESMonitorSet(snes,(myMonitor)SNESMonitorDefault,vf,
-                               (myPetscFunc)PetscViewerAndFormatDestroy);
+                               (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
          PCHKERRQ(snes,ierr);
       }
    }
@@ -4163,20 +4156,31 @@ void PetscNonlinearSolver::SetUpdate(void (*update)(Operator *,int,
 void PetscNonlinearSolver::Mult(const Vector &b, Vector &x) const
 {
    SNES snes = (SNES)obj;
+   MPI_Comm comm = PetscObjectComm(obj);
 
-   bool b_nonempty = b.Size();
-   if (!B) { B = new PetscParVector(PetscObjectComm(obj), *this, true); }
-   if (!X) { X = new PetscParVector(PetscObjectComm(obj), *this, false, false); }
+   // Reduction needed: some processes may have null local size while others don't,
+   // and VecPlaceArray (used by PlaceMemory) is a logically collective operation.
+   PetscBool b_nonempty = b.Size() ? PETSC_TRUE : PETSC_FALSE;
+#if PETSC_VERSION_LT(3,24,0)
+   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPIU_BOOL,MPI_LOR,comm);
+#else
+   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPI_C_BOOL,MPI_LOR,comm);
+#endif
+   CCHKERRQ(comm,mpiierr);
+
+   // Always create B with allocate=false so that PlaceMemory can be called on
+   // it regardless of whether b was empty on a previous call.
+   if (!B) { B = new PetscParVector(comm, *this, true, false); }
+   if (!X) { X = new PetscParVector(comm, *this, false, false); }
    X->PlaceMemory(x.GetMemory(),iterative_mode);
    if (b_nonempty) { B->PlaceMemory(b.GetMemory()); }
-   else { *B = 0.0; }
 
    Customize();
 
    if (!iterative_mode) { *X = 0.; }
 
-   // Solve the system.
-   ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
+   // Solve the system. Pass nullptr for b when empty (PETSc treats it as zero RHS).
+   ierr = SNESSolve(snes, b_nonempty ? B->x : nullptr, X->x); PCHKERRQ(snes, ierr);
    X->ResetMemory();
    if (b_nonempty) { B->ResetMemory(); }
 }
@@ -5329,21 +5333,27 @@ static PetscErrorCode __mfem_pc_shell_destroy(PC pc)
    PetscFunctionReturn(PETSC_SUCCESS);
 }
 
+static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt ptr)
+{
+   PetscErrorCode ierr;
+
+   PetscFunctionBeginUser;
 #if PETSC_VERSION_LT(3,23,0)
-
-static PetscErrorCode __mfem_array_container_destroy(void *ptr)
-{
-   PetscErrorCode ierr;
-
-   PetscFunctionBeginUser;
    ierr = PetscFree(ptr); CHKERRQ(ierr);
+#else
+   ierr = PetscFree(*(void**)ptr); CHKERRQ(ierr);
+#endif
    PetscFunctionReturn(PETSC_SUCCESS);
 }
 
-static PetscErrorCode __mfem_matarray_container_destroy(void *ptr)
+static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt ptr)
 {
+#if PETSC_VERSION_LT(3,23,0)
    mfem::Array<Mat> *a = (mfem::Array<Mat>*)ptr;
-   PetscErrorCode   ierr;
+#else
+   mfem::Array<Mat> *a = *(mfem::Array<Mat>**)ptr;
+#endif
+   PetscErrorCode ierr;
 
    PetscFunctionBeginUser;
    for (int i=0; i<a->Size(); i++)
@@ -5356,41 +5366,16 @@ static PetscErrorCode __mfem_matarray_container_destroy(void *ptr)
    PetscFunctionReturn(PETSC_SUCCESS);
 }
 
+#if PETSC_VERSION_LT(3,23,0)
+static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
 #else
-
-static PetscErrorCode __mfem_array_container_destroy(void **ptr)
+static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt ctx)
+#endif
 {
    PetscErrorCode ierr;
 
    PetscFunctionBeginUser;
-   ierr = PetscFree(*ptr); CHKERRQ(ierr);
-   PetscFunctionReturn(PETSC_SUCCESS);
-}
-
-static PetscErrorCode __mfem_matarray_container_destroy(void **ptr)
-{
-   mfem::Array<Mat> *a = (mfem::Array<Mat>*)*ptr;
-   PetscErrorCode   ierr;
-
-   PetscFunctionBeginUser;
-   for (int i=0; i<a->Size(); i++)
-   {
-      Mat M = (*a)[i];
-      MPI_Comm comm = PetscObjectComm((PetscObject)M);
-      ierr = MatDestroy(&M); CCHKERRQ(comm,ierr);
-   }
-   delete a;
-   PetscFunctionReturn(PETSC_SUCCESS);
-}
-
-#endif
-
-static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
-{
-   PetscErrorCode  ierr;
-
-   PetscFunctionBeginUser;
-   ierr = PetscFree(*ctx); CHKERRQ(ierr);
+   ierr = PetscFree(*(void**)ctx); CHKERRQ(ierr);
    PetscFunctionReturn(PETSC_SUCCESS);
 }
 
@@ -5635,11 +5620,7 @@ static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix* hA,Mat* pA)
 
       ierr = PetscContainerCreate(comm,&c); CHKERRQ(ierr);
       ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
-#if PETSC_VERSION_LT(3,23,0)
-      ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
       ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
       CHKERRQ(ierr);
       ierr = PetscObjectCompose((PetscObject)(*pA),names[i],(PetscObject)c);
       CHKERRQ(ierr);
@@ -5733,11 +5714,7 @@ static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix* hA,Mat* pA)
 
       ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); CHKERRQ(ierr);
       ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
-#if PETSC_VERSION_LT(3,23,0)
-      ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
-#else
       ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-#endif
       CHKERRQ(ierr);
       ierr = PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
       CHKERRQ(ierr);

makefile

@@ -123,15 +123,20 @@ EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
-MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
+MINIAPP_ALL_SUBDIRS = common electromagnetics meshing performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
  hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
- fluids/navier fluids/schrodinger-flow
+ fluids/navier fluids/schrodinger-flow plasma plasma/pic
+MINIAPP_RECURSIVE_SUBDIRS = plasma/pic
+MINIAPP_SUBDIRS := $(filter-out \
+   $(MINIAPP_RECURSIVE_SUBDIRS),$(MINIAPP_ALL_SUBDIRS))
+MINIAPP_ALL_DIRS := $(addprefix miniapps/,$(MINIAPP_ALL_SUBDIRS))
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys shifted dpg diag-smoothers fluids/navier)
+ toys gslib shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
 
+EM_ALL_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_ALL_DIRS)
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
 
 TEST_SUBDIRS = unit
@@ -146,7 +151,7 @@ MFEM_BUILD_DIR ?= .
 BUILD_DIR := $(MFEM_BUILD_DIR)
 BUILD_REAL_DIR := $(abspath $(BUILD_DIR))
 ifneq ($(BUILD_REAL_DIR),$(MFEM_REAL_DIR))
-   BUILD_SUBDIRS = $(DIRS) config $(EM_DIRS) doc $(TEST_DIRS)
+   BUILD_SUBDIRS = $(DIRS) config $(EM_ALL_DIRS) doc $(TEST_DIRS)
    CONFIG_FILE_DEF = -DMFEM_CONFIG_FILE='"$(BUILD_REAL_DIR)/config/_config.hpp"'
    BLD := $(if $(BUILD_REAL_DIR:$(CURDIR)=),$(BUILD_DIR)/,)
    $(if $(word 2,$(BLD)),$(error Spaces in BLD = "$(BLD)" are not supported))
@@ -483,10 +488,10 @@ $(OBJECT_FILES): $(BLD)%.o: $(SRC)%.cpp $(CONFIG_MK)
 
 all: examples miniapps $(TEST_DIRS)
 
-.PHONY: miniapps $(EM_DIRS) $(TEST_DIRS)
+.PHONY: miniapps $(EM_ALL_DIRS) $(TEST_DIRS)
 miniapps: $(MINIAPP_DIRS)
 $(MINIAPP_USE_COMMON): miniapps/common
-$(EM_DIRS) $(TEST_DIRS): lib
+$(EM_ALL_DIRS) $(TEST_DIRS): lib
 	$(MAKE) -C $(BLD)$(@)
 
 .PHONY: doc
@@ -694,7 +699,7 @@ local-config:
 .PHONY: build-config
 build-config:
 	for d in $(BUILD_SUBDIRS); do mkdir -p $(BLD)$${d}; done
-	for dir in "" $(addsuffix /,config $(EM_DIRS) doc $(TEST_DIRS)); do \
+	for dir in "" $(addsuffix /,config $(EM_ALL_DIRS) doc $(TEST_DIRS)); do\
 	   printf "# Auto-generated file.\n%s\n%s\n" \
 	      "MFEM_DIR = $(MFEM_REAL_DIR)" \
 	      "include \$$(MFEM_DIR)/$${dir}makefile" \
@@ -796,13 +801,15 @@ status info:
 
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 3.1"
-FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
+FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_ALL_DIRS) config,$(dir)/*.?pp)
 TESTS_SUBDIRS = unit benchmarks convergence mem_manager par-mesh-format
-UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme
-MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners hooke/materials hooke/kernels
+UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme dfem
+MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners \
+   hooke/materials hooke/kernels
 FORMAT_FILES += $(foreach dir,$(TESTS_SUBDIRS),tests/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(UNIT_TESTS_SUBDIRS),tests/unit/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(MINIAPPS_SUBDIRS),miniapps/$(dir)/*.?pp)
+FORMAT_FILES += config/cmake/config.hpp.in config/config.hpp.in mfem*.hpp
 FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp
 FORMAT_LIST = $(filter-out $(FORMAT_EXCLUDE),$(wildcard $(FORMAT_FILES)))
 
@@ -833,14 +840,29 @@ mfem_check_command = \
 # Verify the C++ code styling in MFEM and check that std::cout and std::cerr are
 # not used in the library (use mfem::out and mfem::err instead).
 style:
-	@echo "Applying C++ code style..."
 	@astyle_version="$$($(ASTYLE_BIN) --version)";\
 	 if [ "$$astyle_version" != $(ASTYLE_VER) ]; then\
 	    printf "%s\n" "Invalid astyle version: '$$astyle_version'"\
 	           "Please use: '"$(ASTYLE_VER)"'";\
 	    exit 1;\
 	 fi
-	@err_code=0;\
+	@err_code=0; \
+	if command -v git 2>&1 > /dev/null && [ -d $(MFEM_DIR)/.git ]; then \
+	   echo "Checking if all git files are selected for formatting ..."; \
+	   ls -1 $(FORMAT_FILES) | sort > format-files-make.txt; \
+	   git -C $(MFEM_DIR) ls-files '*.[ch]pp*' | sort \
+	      > format-files-git.txt; \
+	   cat format-files-make.txt format-files-git.txt | sort | uniq \
+	      > format-files-make-plus-git.txt; \
+	   rm -f format-files-git.txt; \
+	   $(call mfem_check_command,\
+	      diff format-files-make.txt format-files-make-plus-git.txt | \
+	      grep "^> ",\
+	      "All git files are selected for formatting",\
+	      "The above git files are NOT selected for formatting"); \
+	   rm -f format-files-make.txt format-files-make-plus-git.txt; \
+	fi; \
+	echo "Applying C++ code style...";\
 	$(call mfem_check_command,\
 	    $(ASTYLE) $(FORMAT_LIST) | grep Formatted,\
 	    "No source files were changed",\

mesh/mesh.cpp

@@ -493,8 +493,7 @@ void Mesh::GetBdrElementTransformation(int i,
          {
             for (int j = 0; j < n; j++)
             {
-               int idx = vdofs[n*k+j];
-               pm(k,j) = nodes((idx<0)? -1-idx:idx);
+               pm(k,j) = nodes(UnsignIndex(vdofs[n*k+j]));
             }
          }
          ElTr->SetFE(bdr_el);
@@ -1356,7 +1355,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
                face.element[0].conformity = ElementConformity::Coincident;
                face.element[1].conformity = ElementConformity::Coincident;
                face.element[1].location = ElementLocation::FaceNbr;
-               face.element[1].index = -1 - e2;
+               face.element[1].index = FlipIndexSign(e2);
                face.element[1].orientation = inf2%64;
             }
          }
@@ -1379,7 +1378,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
                face.element[1].location = ElementLocation::FaceNbr;
                face.element[0].conformity = ElementConformity::Coincident;
                face.element[1].conformity = ElementConformity::Superset;
-               face.element[1].index = -1 - e2;
+               face.element[1].index = FlipIndexSign(e2);
                face.element[1].orientation = inf2%64;
             }
             face.point_matrix = nc_faces_info[ncface].PointMatrix;
@@ -1405,7 +1404,7 @@ Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
          face.element[1].location = ElementLocation::FaceNbr;
          face.element[0].conformity = ElementConformity::Superset;
          face.element[1].conformity = ElementConformity::Coincident;
-         face.element[1].index = -1 - e2;
+         face.element[1].index = FlipIndexSign(e2);
          face.element[1].orientation = inf2%64;
          face.point_matrix = nc_faces_info[ncface].PointMatrix;
       }
@@ -1438,7 +1437,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
          break;
       case FaceInfoTag::SharedConforming:
          res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
          res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
          res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
          break;
@@ -1448,7 +1447,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
          break;
       case FaceInfoTag::SharedSlaveNonconforming:
          res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
          res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
          res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
          break;
@@ -1456,7 +1455,7 @@ Mesh::FaceInformation::operator Mesh::FaceInfo() const
          break;
       case FaceInfoTag::GhostSlave:
          res.Elem1No = element[0].index;
-         res.Elem2No = -1 - element[1].index;
+         res.Elem2No = FlipIndexSign(element[1].index);
          res.Elem1Inf = element[0].orientation + element[0].local_face_id*64;
          res.Elem2Inf = element[1].orientation + element[1].local_face_id*64;
          break;
@@ -1617,7 +1616,9 @@ Element::Type Mesh::GetFaceElementType(int Face) const
 
 Array<int> Mesh::GetFaceToBdrElMap() const
 {
-   Array<int> face_to_be(Dim == 2 ? NumOfEdges : NumOfFaces);
+   Array<int> face_to_be(Dim == 1 ? NumOfVertices :
+                         Dim == 2 ? NumOfEdges :
+                         Dim == 3 ? NumOfFaces : 0);
    face_to_be = -1;
    for (int i = 0; i < NumOfBdrElements; i++)
    {
@@ -6514,7 +6515,7 @@ void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_ukv)
          input >> edge_to_ukv[j] >> v[0] >> v[1];
          if (v[0] > v[1])
          {
-            edge_to_ukv[j] = -1 - edge_to_ukv[j];
+            edge_to_ukv[j] = FlipIndexSign(edge_to_ukv[j]);
          }
       }
    }
@@ -6551,9 +6552,6 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
    const int NP = NumOfElements; // number of patches
    const int NPKV = NP * dim;    // number of patch knotvectors
    constexpr int notset = -9999999;
-   // Sign convention
-   auto flipSign = [](int i) { return -1 - i; };
-   auto unSign = [](int i) { return (i < 0) ? -1 - i : i; };
    // Local edge index -> dimension convention
    auto edge_to_dim = [](int i) { return (i < 8) ? ((i & 1) ? 1 : 0) : 2; };
 
@@ -6569,7 +6567,7 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
       {
          GetElementVertices(i, v);
          // Sign is based on the edge's vertex indices
-         edge_to_ukv[i] = (v[1] > v[0]) ? i : flipSign(i);
+         edge_to_ukv[i] = (v[1] > v[0]) ? i : FlipIndexSign(i);
          ukv_to_rpkv[i] = i;
       }
       return;
@@ -6619,14 +6617,14 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
          // We've set this edge already - link this index to it
          if (edge_to_pkv[edge] != notset)
          {
-            const int pkv_other = unSign(edge_to_pkv[edge]);
+            const int pkv_other = UnsignIndex(edge_to_pkv[edge]);
             unite(pkv, pkv_other);
          }
          else
          {
             GetEdgeVertices(edge, v);
             // Sign is based on the edge's vertex indices
-            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : flipSign(pkv);
+            edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : FlipIndexSign(pkv);
          }
       }
    }
@@ -6653,10 +6651,10 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
    edge_to_ukv.SetSize(NumOfEdges);
    for (int i = 0; i < NumOfEdges; i++)
    {
-      const int pkv = unSign(edge_to_pkv[i]);
+      const int pkv = UnsignIndex(edge_to_pkv[i]);
       const int rpkv = pkv_to_rpkv[pkv];
       const int ukv = rpkv_to_ukv[rpkv];
-      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? flipSign(ukv) : ukv;
+      edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? FlipIndexSign(ukv) : ukv;
    }
 
    CorrectPatchTopoOrientations(edge_to_ukv);
@@ -6667,9 +6665,6 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
    const int dim = Dimension(); // Topological (not physical) dimension
    if (dim == 1) { return; }
 
-   // Sign convention
-   auto flipSign = [](int i) { return -1 - i; };
-
    const Table *face2elem = GetFaceToElementTable();
    Array<int> pfaces, orient;
    Array<int> fe, feo;
@@ -6688,7 +6683,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          for (auto e : fe)
          {
             const int skv = edge_to_ukv[e];
-            if (skv == kv || flipSign(skv) == kv) { hasKV = true; }
+            if (skv == kv || FlipIndexSign(skv) == kv) { hasKV = true; }
          }
          if (hasKV)
          {
@@ -6718,7 +6713,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       };
    }
 
-   Array<int> ukvs((dim==2) ? 4 : 12);
+   Array<int> ukvs((dim == 2) ? 4 : 12);
    Array<int> pe, oe;
    bool initKV = false;
 
@@ -6732,7 +6727,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       for (int i = 0; i < pe.Size(); i++)
       {
          ukvs[i] = edge_to_ukv[pe[i]];
-         ukvs[i] = (oe[i] < 0) ? flipSign(ukvs[i]) : ukvs[i];
+         ukvs[i] = (oe[i] < 0) ? FlipIndexSign(ukvs[i]) : ukvs[i];
       }
 
       // Find the direction with this kv.
@@ -6740,12 +6735,19 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       for (int d=0; d<dim; ++d) // Loop over directions.
       {
          const int skv = edge_to_ukv[pe[dir_edges[d][0]]];
-         if (skv == kv || flipSign(skv) == kv)
+         if (skv == kv || FlipIndexSign(skv) == kv)
          {
-            thisDir = d;
+            for (auto e : dir_edges[d])
+               if (!edgeSet[pe[e]])
+               {
+                  thisDir = d;
+               }
          }
       }
-      MFEM_VERIFY(thisDir >= 0, "");
+      if (thisDir == -1)
+      {
+         return false;
+      }
 
       // For this direction, find any edge already set. If no edge is set, we
       // arbitrarily take the first.
@@ -6777,12 +6779,12 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          }
 
          const int edge = pe[i];
-         if ((dim == 2 && ukvs[i] != flipSign(ukvs[ref_edge0])) ||
-             (dim == 3 && ukvs[i] == flipSign(ukvs[ref_edge0])))
+         if ((dim == 2 && ukvs[i] != FlipIndexSign(ukvs[ref_edge0])) ||
+             (dim == 3 && ukvs[i] == FlipIndexSign(ukvs[ref_edge0])))
          {
             // Flip the sign of this edge
-            MFEM_VERIFY(!edgeSet[edge], "");
-            edge_to_ukv[edge] = flipSign(edge_to_ukv[edge]);
+            MFEM_ASSERT(!edgeSet[edge], "");
+            edge_to_ukv[edge] = FlipIndexSign(edge_to_ukv[edge]);
          }
 
          edgeSet[edge] = true;
@@ -6827,10 +6829,11 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       int unsetDim = -1;
       for (int d=0; d<dim; ++d) // Loop over dimensions.
       {
-         if (!edgeSet[pe[dir_edges[d][0]]])
-         {
-            unsetDim = d;
-         }
+         for (auto e : dir_edges[d])
+            if (!edgeSet[pe[e]])
+            {
+               unsetDim = d;
+            }
       }
 
       if (unsetDim == -1)
@@ -6839,9 +6842,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          continue;
       }
 
-      const int kv_signed = edge_to_ukv[pe[dir_edges[unsetDim][0]]];
-      const int kv = kv_signed < 0 ? flipSign(kv_signed) : kv_signed;
-      MFEM_VERIFY(!edgeSet[pe[dir_edges[unsetDim][0]]], "");
+      const int kv = UnsignIndex(edge_to_ukv[pe[dir_edges[unsetDim][0]]]);
 
       initKV = false;
 
@@ -6891,6 +6892,7 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
       }
    }
 
+#ifdef MFEM_DEBUG
    bool allSet = true;
    for (auto eset : edgeSet)
    {
@@ -6899,7 +6901,8 @@ void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
          allSet = false;
       }
    }
-   MFEM_VERIFY(allSet && unset.size() == 0, "Some edge is not set");
+   MFEM_ASSERT(allSet && unset.size() == 0, "Some edge is not set");
+#endif
 
    delete face2elem;
 }
@@ -6941,7 +6944,7 @@ void Mesh::LoadNonconformingPatchTopo(std::istream &input,
 
       if (v[0] > v[1])
       {
-         ukv = -1 - ukv;
+         ukv = FlipIndexSign(ukv);
       }
       edge_to_ukv[j] = ukv;
    }
@@ -12438,11 +12441,7 @@ void Mesh::PrintTopoEdges(std::ostream &os, const Array<int> &e_to_k,
    for (int i = 0; i < NumOfEdges; i++)
    {
       edge_vertex->GetRow(i, vert);
-      int ki = e_to_k[i];
-      if (ki < 0)
-      {
-         ki = -1 - ki;
-      }
+      const int ki = UnsignIndex(e_to_k[i]);
 
       if (vmap)
       {
@@ -15751,9 +15750,18 @@ Mesh PartitionMPI(int dim, int mpi_cnt, int elem_per_mpi, bool print,
 {
    MFEM_VERIFY(dim > 1, "Not implemented for 1D meshes.");
 
-   auto factor = [&](int N)
+   // Closest int divisor to the cubit root, going down.
+   auto factor3 = [](int N)
    {
-      for (int i = static_cast<int>(sqrt(N)); i > 0; i--)
+      for (int i = static_cast<int>(round(cbrt(N))); i > 0; i--)
+      {  if (N % i == 0) { return i; } }
+      return 1;
+   };
+
+   // Closest int divisor to the square root, going down.
+   auto factor2 = [](int N)
+   {
+      for (int i = static_cast<int>(round(sqrt(N))); i > 0; i--)
       { if (N % i == 0) { return i; } }
       return 1;
    };
@@ -15777,22 +15785,22 @@ Mesh PartitionMPI(int dim, int mpi_cnt, int elem_per_mpi, bool print,
    int el0_x, el0_y, el0_z;
    if (dim == 2)
    {
-      mpi_x = factor(mpi_cnt);
+      mpi_x = factor2(mpi_cnt);
       mpi_y = mpi_cnt / mpi_x;
 
       // Switch order for better balance.
-      el0_y = factor(el0);
+      el0_y = factor2(el0);
       el0_x = el0 / el0_y;
    }
    else
    {
-      mpi_x = factor(mpi_cnt);
-      mpi_y = factor(mpi_cnt / mpi_x);
+      mpi_x = factor3(mpi_cnt);
+      mpi_y = factor2(mpi_cnt / mpi_x);
       mpi_z = mpi_cnt / mpi_x / mpi_y;
 
       // Switch order for better balance.
-      el0_z = factor(el0);
-      el0_y = factor(el0 / el0_z);
+      el0_z = factor3(el0);
+      el0_y = factor2(el0 / el0_z);
       el0_x = el0 / el0_y / el0_z;
    }
 

mesh/mesh.hpp

@@ -3206,10 +3206,22 @@ public:
 
 
 /// Extrude a 1D mesh
+/**
+ * @param mesh      1D mesh
+ * @param ny        number of transverse elements of the extruded mesh
+ * @param sy        physical size in the direction of extrusion
+ * @param closed    if false, only the original boundaries are extruded,
+ *                  otherwise boundaries are generated all around the domain
+ */
 Mesh *Extrude1D(Mesh *mesh, const int ny, const real_t sy,
                 const bool closed = false);
 
 /// Extrude a 2D mesh
+/**
+ * @param mesh      2D mesh
+ * @param nz        number of transverse elements of the extruded mesh
+ * @param sz        physical size in the direction of extrusion
+ */
 Mesh *Extrude2D(Mesh *mesh, const int nz, const real_t sz);
 
 /** @brief Constructs the smallest possible [0,1]^dim serial mesh that can be

mesh/mesh_operators.cpp

@@ -63,7 +63,6 @@ ThresholdRefiner::ThresholdRefiner(ErrorEstimator &est)
 
    threshold = 0.0;
    num_marked_elements = 0LL;
-   current_sequence = -1;
 
    non_conforming = -1;
    nc_limit = 0;
@@ -87,7 +86,6 @@ int ThresholdRefiner::MarkWithoutRefining(Mesh & mesh,
    threshold = 0.0;
    num_marked_elements = 0LL;
    refinements.SetSize(0);
-   current_sequence = mesh.GetSequence();
 
    const long long num_elements = mesh.GetGlobalNE();
    if (num_elements >= max_elements) { return STOP; }
@@ -143,13 +141,12 @@ int ThresholdRefiner::ApplyImpl(Mesh &mesh)
    if (num_marked_elements == 0LL) { return STOP; }
 
    mesh.GeneralRefinement(marked_elements, non_conforming, nc_limit);
-   return CONTINUE + REFINED;
+   return static_cast<int>(CONTINUE) + static_cast<int>(REFINED);
 }
 
 void ThresholdRefiner::Reset()
 {
    estimator.Reset();
-   current_sequence = -1;
    num_marked_elements = 0LL;
    // marked_elements.SetSize(0); // not necessary
 }
@@ -162,7 +159,7 @@ int ThresholdDerefiner::ApplyImpl(Mesh &mesh)
    const Vector &local_err = estimator.GetLocalErrors();
    bool derefs = mesh.DerefineByError(local_err, threshold, nc_limit, op);
 
-   return derefs ? CONTINUE + DEREFINED : NONE;
+   return derefs ? static_cast<int>(CONTINUE) + static_cast<int>(DEREFINED) : NONE;
 }
 
 
@@ -290,7 +287,7 @@ int CoefficientRefiner::PreprocessMesh(Mesh &mesh, int max_it)
    }
    delete l2fes;
    delete gf;
-   return CONTINUE + REFINED;
+   return static_cast<int>(CONTINUE) + static_cast<int>(REFINED);
 
 }
 
@@ -310,7 +307,7 @@ int Rebalancer::ApplyImpl(Mesh &mesh)
    if (pmesh && pmesh->Nonconforming())
    {
       pmesh->Rebalance();
-      return CONTINUE + REBALANCED;
+      return static_cast<int>(CONTINUE) + static_cast<int>(REBALANCED);
    }
 #endif
    return NONE;

mesh/mesh_operators.hpp

@@ -188,7 +188,6 @@ protected:
    long long num_marked_elements;
 
    Array<Refinement> marked_elements;
-   long current_sequence;
 
    int non_conforming;
    int nc_limit;

mesh/mesh_readers.cpp

@@ -1516,12 +1516,15 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
 void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
 {
    string buff;
-   real_t version;
+   string version;
    int binary, dsize;
    input >> version >> binary >> dsize;
-   if (version < 2.2)
+   if (version != "2.2")
    {
-      MFEM_ABORT("Gmsh file version < 2.2");
+      MFEM_ABORT("Gmsh file version must be 2.2, found version "
+                 << version << ".\n"
+                 "To convert your mesh to the required format, use:\n"
+                 "  gmsh -format msh22 -save -o output.msh input.msh");
    }
    if (dsize != sizeof(double))
    {

mesh/ncmesh.cpp

@@ -3542,7 +3542,7 @@ void NCMesh::TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
 
             // create a slave face record with a degenerate point matrix
             face_list.slaves.Append(
-               Slave(-1 - enode.edge_index,
+               Slave(FlipIndexSign(enode.edge_index),
                      eid[0].element, eid[0].local, Geometry::SQUARE));
             Slave &sl = face_list.slaves.Last();
 
@@ -3589,7 +3589,7 @@ void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
          // non-slave edge is really a (face-)slave itself.
          const MeshId &eid = *eid_and_type.id;
          face_list.slaves.Append(
-            Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
+            Slave(FlipIndexSign(eid.index), eid.element, eid.local, Geometry::TRIANGLE));
 
          int v0index = nodes[vn0].vert_index;
          int v1index = nodes[vn1].vert_index;

mesh/ncnurbs.cpp

@@ -93,7 +93,7 @@ void NCNURBSExtension::GetMasterEdgeEntities(
       }
       else
       {
-         const int auxEdge = -1 - edge_i;
+         const int auxEdge = FlipIndexSign(edge_i);
          GetAuxEdgeVertices(auxEdge, sverts);
       }
 
@@ -159,7 +159,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
             {
                if (edge < 0)
                {
-                  sideAuxEdges[s].Append(-1 - edge);
+                  sideAuxEdges[s].Append(FlipIndexSign(edge));
                }
                else
                {
@@ -456,7 +456,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                                        == afverts[j], "");
                         }
 
-                        ori_f2 = -1 - ori_f2;
+                        ori_f2 = FlipIndexSign(ori_f2);
                      }
                      else
                      {
@@ -468,7 +468,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                      }
 
                      facePairs.emplace_back(FacePairInfo{fverts[vMinID], f,
-                     SlaveFaceInfo{-1 - afid, ori_f2,
+                     SlaveFaceInfo{FlipIndexSign(afid), ori_f2,
                         {fki(vMinID,0), fki(vMinID,1)},
                         {
                            fki((vMinID + 2) % 4,0) - fki(vMinID,0),
@@ -509,7 +509,7 @@ void NCNURBSExtension::FindAdditionalFacesSA(
                      auxFaces.push_back(auxFace);
 
                      facePairs.emplace_back(FacePairInfo{fverts[vMinID], f,
-                     SlaveFaceInfo{-1 - auxFaceId, ori_f,
+                     SlaveFaceInfo{FlipIndexSign(auxFaceId), ori_f,
                         {fki(vMinID,0), fki(vMinID,1)},
                         {
                            fki((vMinID + 2) % 4,0) - fki(vMinID,0),
@@ -622,7 +622,7 @@ void NCNURBSExtension::GetAuxFaceEdges(int auxFace, Array<int> &edges) const
       }
       else // Auxiliary edge
       {
-         edges[i] = -1 - auxv2e.at(edge_v);
+         edges[i] = FlipIndexSign(auxv2e.at(edge_v));
       }
    }
 }
@@ -633,7 +633,7 @@ int OffsetHelper(int i, int j, const Array<int> &a, const Array<int> &b)
 {
    if (i < 0)
    {
-      return b[-1 - i + j];
+      return b[FlipIndexSign(i) + j];
    }
    else if (i + j < a.Size())
    {
@@ -679,7 +679,7 @@ void NCNURBSExtension::GetMasterEdgeDofs(bool dof, int me,
       }
       else // Auxiliary edge
       {
-         GetAuxEdgeVertices(-1 - slaveId, svert);
+         GetAuxEdgeVertices(FlipIndexSign(slaveId), svert);
       }
 
       bool reverse = false;
@@ -872,7 +872,7 @@ void ReorderArray2D(int i0, int j0, const Array2D<int> &a,
 // Set a quadrilateral vertex index permutation for a given orientation.
 void GetVertexOrdering(int ori, std::array<int, 4> &perm)
 {
-   const int oriAbs = ori < 0 ? -1 - ori : ori;
+   const int oriAbs = UnsignIndex(ori);
 
    for (int i=0; i<4; ++i)
    {
@@ -1094,7 +1094,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
       if (slaveId < 0)
       {
          // Auxiliary face
-         const int auxFace = -1 - slaveId;
+         const int auxFace = FlipIndexSign(slaveId);
 
          // Set slave face entity dimensions.
          if (dof)
@@ -1171,7 +1171,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
             }
             else
             {
-               const int auxEdge = -1 - edge;
+               const int auxEdge = FlipIndexSign(edge);
                GetAuxEdgeVertices(auxEdge, evert);
             }
             MFEM_ASSERT(evert[0] == vstart || evert[1] == vstart, "");
@@ -1184,7 +1184,7 @@ void NCNURBSExtension::GetMasterFaceDofs(bool dof, int mf,
          // dimensions of the master face, by using ori.
          int e1 = -1, e2 = -1;
          {
-            const int aori = ori < 0 ? -1 - ori : ori;
+            const int aori = UnsignIndex(ori);
             if (aori % 2 == 0)
             {
                e1 = 0;
@@ -1416,14 +1416,15 @@ void NCNURBSExtension::ProcessVertexToKnot2D(const VertexToKnotSpan &v2k,
          {
             // Create a new auxiliary edge
             auxv2e[childPair] = auxEdges.size();
-            auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ?
-                                                parentEdge : -1 - parentEdge,
+            auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ? parentEdge :
+                                                FlipIndexSign(parentEdge),
             {childPair.first, childPair.second},
             {newParentEdge ? 0 : prevKI, ks}});
          }
       }
 
-      const int childEdge = childPairTopo ? v2e[childPair] : -1 - auxv2e[childPair];
+      const int childEdge = childPairTopo ? v2e[childPair] :
+                            FlipIndexSign(auxv2e[childPair]);
 
       // Check whether this is the final vertex in this parent edge. Note that
       // the logic for comparing (pv[0],pv[1]) to the next parents assumes the
@@ -1460,14 +1461,15 @@ void NCNURBSExtension::ProcessVertexToKnot2D(const VertexToKnotSpan &v2k,
 
                // -1 denotes `ne` at endpoint
                auxEdges.emplace_back(AuxiliaryEdge{pv[0] < pv[1] ?
-                                                   -1 - parentEdge : parentEdge,
+                                                   FlipIndexSign(parentEdge) :
+                                                   parentEdge,
                {finalChildPair.first, finalChildPair.second},
                {ks, -1}});
             }
          }
 
          const int finalChildEdge = finalChildPairTopo ? v2e[finalChildPair] :
-                                    -1 - auxv2e[finalChildPair];
+                                    FlipIndexSign(auxv2e[finalChildPair]);
          edgePairs.emplace_back(-1, -1, finalChildEdge, parentEdge);
       }
 
@@ -1805,7 +1807,7 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
                   auxFaces.push_back(auxFace);
                   facePairs.emplace_back(
                      FacePairInfo{cv[0], parentFace,
-                                  SlaveFaceInfo{-1 - auxv2f[childPair],
+                                  SlaveFaceInfo{FlipIndexSign(auxv2f[childPair]),
                                                 0, {i0, j0}, {d0, d1}}});
                }
             }
@@ -2111,7 +2113,7 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
                      auxv2e[childPair] = auxEdges.size();
                      auxEdges.emplace_back(AuxiliaryEdge{pv0 < pv1 ?
                                                          parentEdge :
-                                                         -1 - parentEdge,
+                                                         FlipIndexSign(parentEdge),
                      {childPair.first, childPair.second},
                      {knotIndex0, knotIndex1}});
                   }
@@ -2131,7 +2133,8 @@ void NCNURBSExtension::ProcessVertexToKnot3D(
 
                   const EdgePairInfo ep_e((e_idx == n_d - de) ? -1 : tv,
                                           (e_idx == n_d - de) ? -1 : tvki,
-                                          -1 - auxv2e[childPair], parentEdge);
+                                          FlipIndexSign(auxv2e[childPair]),
+                                          parentEdge);
 
                   const bool unset = !edgePairs[edgePairOS[parentEdge] + e_idx].isSet;
                   if (unset)
@@ -2226,7 +2229,7 @@ void NCNURBSExtension::GetAuxFaceToPatchTable(Array2D<int> &auxface2patch)
                if (s < 0)
                {
                   // Auxiliary face.
-                  const int aux = -1 - s;
+                  const int aux = FlipIndexSign(s);
                   if (auxface2patch(aux, 0) >= 0)
                   {
                      if (auxface2patch(aux, 1) != -1) { consistent = false; }
@@ -2316,7 +2319,7 @@ void NCNURBSExtension::UpdateAuxiliaryKnotSpans(const Array<int> &rf)
    for (auto auxEdge : auxEdges)
    {
       const int p = auxEdge.parent;
-      const int parent = p < 0 ? -1 - p : p;
+      const int parent = UnsignIndex(p);
       const int kv = KnotInd(parent);
       for (int i=0; i<2; ++i)
       {
@@ -2382,14 +2385,8 @@ int NCNURBSExtension::AuxiliaryEdgeNE(int aux_edge)
    const int signedParentEdge = auxEdges[aux_edge].parent;
    const int ki0 = auxEdges[aux_edge].ksi[0];
    const int ki1raw = auxEdges[aux_edge].ksi[1];
-   int ki1 = ki1raw;
-   if (ki1raw == -1)
-   {
-      const bool rev = signedParentEdge < 0;
-      const int parentEdge = rev ? -1 - signedParentEdge : signedParentEdge;
-      ki1 = KnotVec(parentEdge)->GetNE();
-   }
-
+   const int ki1 = ki1raw == -1 ? KnotVec(UnsignIndex(signedParentEdge))->GetNE()
+                   : ki1raw;
    return ki1 - ki0;
 }
 
@@ -2403,7 +2400,7 @@ void NCNURBSExtension::SlaveEdgeToParent(int se, int parent,
    Array<int> sev(2);
    if (se < 0) // Auxiliary edge
    {
-      for (int i=0; i<2; ++i) { sev[i] = auxEdges[-1 - se].v[i]; }
+      for (int i=0; i<2; ++i) { sev[i] = auxEdges[FlipIndexSign(se)].v[i]; }
    }
    else
    {
@@ -2459,7 +2456,7 @@ void NCNURBSExtension::GetMasterEdgePieceOffsets(int mid, Array<int> &os)
       }
       else
       {
-         nes = AuxiliaryEdgeNE(-1 - s);
+         nes = AuxiliaryEdgeNE(FlipIndexSign(s));
       }
 
       os[i+1] = os[i] + nes;
@@ -2565,7 +2562,7 @@ int NCNURBSExtension::SetPatchFactors(int p)
                }
                else // Aux edge
                {
-                  const int aux_edge = -1 - s;
+                  const int aux_edge = FlipIndexSign(s);
                   if (auxef[aux_edge].Size() == 0)
                   {
                      auxef[aux_edge].SetSize(AuxiliaryEdgeNE(aux_edge));
@@ -2611,7 +2608,7 @@ int NCNURBSExtension::SetPatchFactors(int p)
    }
 
    MFEM_VERIFY(consistent, "");
-   return partialChange ? -1 - dirSet : dirSet;
+   return partialChange ? FlipIndexSign(dirSet) : dirSet;
 }
 
 void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
@@ -2715,7 +2712,7 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
                if (s < 0)
                {
                   // Auxiliary face.
-                  const int aux = -1 - s;
+                  const int aux = FlipIndexSign(s);
                   for (int i=0; i<2; ++i)
                   {
                      const int patch = auxface2patch(aux, i);
@@ -2763,7 +2760,7 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
 
          const int dirSetSigned = SetPatchFactors(p);
          const bool partialChange = dirSetSigned < 0;
-         const int dirSet = partialChange ? -1 - dirSetSigned : dirSetSigned;
+         const int dirSet = UnsignIndex(dirSetSigned);
          const bool changed = (patchState[p] != dirSet) || partialChange;
          patchState[p] = dirSet;
 
@@ -2806,8 +2803,8 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
             {
                const int dirSetSigned_i = SetPatchFactors(i);
                const bool partialChange_i = dirSetSigned_i < 0;
-               const int dirSet_i = partialChange_i ? -1 - dirSetSigned_i :
-                                    dirSetSigned_i;
+               const int dirSet_i = partialChange_i ?
+                                    FlipIndexSign(dirSetSigned_i) : dirSetSigned_i;
                const bool changed_i = (patchState[i] != dirSet_i) ||
                                       partialChange_i;
                patchState[p] = dirSet_i;
@@ -3027,7 +3024,7 @@ int GetFaceOrientation(const Mesh *mesh, const int face,
 
    // Check whether ordering is reversed.
    const bool rev = verts[(s + 1) % 4] != fverts[1];
-   if (rev) { s = -1 - s; } // Reversed order is encoded by the sign.
+   if (rev) { s = FlipIndexSign(s); } // Reversed order is encoded by the sign.
    return s;
 }
 
@@ -3040,7 +3037,7 @@ int GetFaceOrientation(const Mesh *mesh, const int face,
 // see GetFaceOrientation.
 bool Reorder2D(int ori, std::array<int, 2> &s0)
 {
-   const int shift = ori < 0 ? -1 - ori : ori;
+   const int shift = UnsignIndex(ori);
 
    // Shift is an F1 index in the counter-clockwise ordering of 4 quad vertices.
    // Now find the (i,j) indices of this index, with i,j in {0,1}.
@@ -3064,7 +3061,7 @@ void GetInverseShiftedDimensions2D(int signedShift, int sm, int sn, int &m,
                                    int &n)
 {
    const bool rev = (signedShift < 0);
-   const int shift = rev ? -1 - signedShift : signedShift;
+   const int shift = UnsignIndex(signedShift);
    MFEM_ASSERT(0 <= shift && shift < 4, "");
 
    // We consider 8 cases for the possible values of rev and shift.
@@ -3136,7 +3133,7 @@ void GetShiftedGridPoints2D(int m, int n, int i, int j, int signedShift,
                             int& sm, int& sn, int& si, int& sj)
 {
    const bool rev = (signedShift < 0);
-   const int shift = rev ? -1 - signedShift : signedShift;
+   const int shift = UnsignIndex(signedShift);
    MFEM_ASSERT(0 <= shift && shift < 4, "");
 
    // (0,0) <= (i,j) < (m,n) are old indices, and old vertex [shift] maps
@@ -3798,8 +3795,7 @@ void NCNURBSExtension::GenerateOffsets()
       const int signedParentEdge = auxEdges[e].parent;
       const int ki0 = auxEdges[e].ksi[0];
       const int ki1raw = auxEdges[e].ksi[1];
-      const bool rev = signedParentEdge < 0;
-      const int parentEdge = rev ? -1 - signedParentEdge : signedParentEdge;
+      const int parentEdge = UnsignIndex(signedParentEdge);
       const int masterNE = KnotVec(parentEdge)->GetNE();
       const int ki1 = ki1raw == -1 ? masterNE : ki1raw;
       const int perEdgeCP = GetNCPperEdge(KnotVec(e));

mesh/nurbs.cpp

@@ -3723,10 +3723,7 @@ bool NURBSExtension::CheckPatches()
       for (int i = 0; i < edges.Size(); i++)
       {
          edges[i] = edge_to_ukv[edges[i]];
-         if (oedge[i] < 0)
-         {
-            edges[i] = -1 - edges[i];
-         }
+         if (oedge[i] < 0) { edges[i] = FlipIndexSign(edges[i]); }
       }
 
       // In 2d - opposite edges must be same knotvector with opposite sign.
@@ -3740,7 +3737,7 @@ bool NURBSExtension::CheckPatches()
       // {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7} for Geometry::CUBE in 3D
       // See fem/geom.cpp for these definitions.
       if ((dim == 2 &&
-           (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
+           (edges[0] != FlipIndexSign(edges[2]) || edges[1] != FlipIndexSign(edges[3]))) ||
 
           (dim == 3 &&
            (edges[0] != edges[2] || edges[0] != edges[4] ||
@@ -3769,7 +3766,7 @@ void NURBSExtension::CheckBdrPatches()
          edges[i] = edge_to_ukv[edges[i]];
          if (oedge[i] < 0)
          {
-            edges[i] = -1 - edges[i];
+            edges[i] = FlipIndexSign(edges[i]);
          }
       }
 
@@ -4766,14 +4763,13 @@ void NURBSExtension::GenerateBdrElementDofTable()
    SetPatchToBdrElements();
 
    int *dof = bel_dof->GetJ();
-   int ndof = bel_dof->Size_of_connections();
+   const int ndof = bel_dof->Size_of_connections();
    for (int i = 0; i < ndof; i++)
    {
-      int idx = dof[i];
+      const int idx = dof[i];
       if (idx < 0)
       {
-         dof[i] = -1 - (activeDof[-1-idx] - 1);
-         dof[i] = -activeDof[-1-idx];
+         dof[i] = -activeDof[FlipIndexSign(idx)];
       }
       else
       {
@@ -4858,12 +4854,12 @@ void NURBSExtension::Generate2DBdrElementDofTable()
                   for (int ii = 0; ii <= ord0; ii++)
                   {
                      conn.to = DofMap(p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]);
-                     if (s == -1) { conn.to = -1 -conn.to; }
+                     if (s == -1) { conn.to = FlipIndexSign(conn.to); }
                      bel_dof_list.Append(conn);
                   }
                }
                bel_to_patch[lbe] = b;
-               bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
+               bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : FlipIndexSign(i);
                lbe++;
             }
             gbe++;
@@ -4936,14 +4932,14 @@ void NURBSExtension::Generate3DBdrElementDofTable()
                            {
                               const int ii_ = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
                               conn.to = DofMap(p2g(ii_, jj_));
-                              if (s == -1) { conn.to = -1 -conn.to; }
+                              if (s == -1) { conn.to = FlipIndexSign(conn.to); }
                               bel_dof_list.Append(conn);
                            }
                         }
                      }
                      bel_to_patch[lbe] = b;
-                     bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
-                     bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
+                     bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : FlipIndexSign(i);
+                     bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : FlipIndexSign(j);
                      lbe++;
                   }
                   gbe++;

mesh/nurbs.hpp

@@ -220,7 +220,7 @@ public:
        @a u.
        The main purpose of this function is its use in FindInterpolant.
        Use GetBotella instead for each shape function separately, perhaps in
-       conjuction with GetSpan and GetRefPoint.*/
+       conjunction with GetSpan and GetRefPoint.*/
    MFEM_DEPRECATED void FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const;
 
    /** @brief Global curve interpolation through the points @a x (overwritten).
@@ -1398,8 +1398,7 @@ inline const real_t &NURBSPatch::operator()(int i, int j, int k, int l) const
 
 inline int NURBSExtension::KnotInd(int edge) const
 {
-   const int kv = edge_to_ukv[edge];
-   return kv >= 0 ? kv : -1 - kv;
+   return UnsignIndex(edge_to_ukv[edge]);
 }
 
 inline int NURBSExtension::KnotSign(int edge) const
@@ -1429,7 +1428,7 @@ const
    else
    {
       *okv = -oedge;
-      return knotVectors[-1-kv];
+      return knotVectors[FlipIndexSign(kv)];
    }
 }
 

mesh/pmesh.cpp

@@ -5639,6 +5639,12 @@ Mesh ParMesh::GetSerialMesh(int save_rank) const
       }
    }
 
+   if (MyRank == save_rank)
+   {
+      attribute_sets.Copy(serialmesh.attribute_sets);
+      bdr_attribute_sets.Copy(serialmesh.bdr_attribute_sets);
+   }
+
    MPI_Barrier(MyComm);
    return serialmesh;
 }

mesh/pncmesh.cpp

@@ -400,7 +400,7 @@ void ParNCMesh::MakeSharedList(const NCList &list, NCList &shared)
          }
          else // special case: prism edge-face constraint
          {
-            if (entity_owner[1][-1-si] != MyRank)
+            if (entity_owner[1][FlipIndexSign(si)] != MyRank)
             {
                master_flag |= 0x2;
             }
@@ -571,9 +571,10 @@ void ParNCMesh::CalculatePMatrixGroups()
       ranks.SetSize(0);
       for (int j = master_face.slaves_begin; j < master_face.slaves_end; j++)
       {
-         int si = face_list.slaves[j].index;
-         int owner = (si >= 0) ? entity_owner[2][si] // standard face dependency
-                     /*     */ : entity_owner[1][-1 - si]; // prism edge-face dep
+         const int si = face_list.slaves[j].index;
+         const int owner =
+            (si >= 0) ? entity_owner[2][si] : // standard face dependency
+            entity_owner[1][FlipIndexSign(si)]; // prism edge-face dep
          ranks.Append(groups[owner][0]);
       }
       ranks.Sort();
@@ -1181,7 +1182,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
       if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }
 
       Mesh::FaceInfo &fi = pmesh.faces_info[cf.index];
-      fi.Elem2No = -1 - fnbr_index[e[0]->index - NElements];
+      fi.Elem2No = FlipIndexSign(fnbr_index[e[0]->index - NElements]);
 
       if (Dim == 3)
       {
@@ -1270,7 +1271,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
                // In other words, side 1 IS the side that generated the face.
             }
             MFEM_ASSERT(fi.Elem2No >= NElements, "");
-            fi.Elem2No = -1 - fnbr_index[fi.Elem2No - NElements];
+            fi.Elem2No = FlipIndexSign(fnbr_index[fi.Elem2No - NElements]);
 
             const DenseMatrix* pm = full_list.point_matrices[sf.geom][sf.matrix];
             if (!sloc && Dim == 3)
@@ -2286,7 +2287,7 @@ void ParNCMesh::Derefine(const Array<int> &derefs)
       if (element_type[index] == 0)
       {
          // this coarse element will get pruned, encode who owns it now
-         index = -1 - elements[coarse[i]].rank;
+         index = FlipIndexSign(elements[coarse[i]].rank);
       }
       transforms.embeddings[i].parent = index;
    }

mesh/submesh/psubmesh.hpp

@@ -227,15 +227,29 @@ public:
                                            const ParGridFunction &dst);
 
    /**
-   * @brief Check if ParMesh @a m is a ParSubMesh.
+   * @brief Check if Mesh @a m is a ParSubMesh.
    *
-   * @param m The input ParMesh
+   * @param m The input Mesh
    */
-   static bool IsParSubMesh(const ParMesh *m)
+   static bool IsParSubMesh(const Mesh *m)
    {
       return dynamic_cast<const ParSubMesh *>(m) != nullptr;
    }
 
+   /**
+   * @brief Check if Mesh @a sub is a ParSubMesh of Mesh @a parent.
+   *
+   * @param sub The potential submesh Mesh
+   * @param parent The potential parent Mesh
+   */
+   static bool IsParSubMesh(const Mesh* sub, const Mesh* parent)
+   {
+      while (IsParSubMesh(sub) &&
+             (sub = static_cast<const ParSubMesh *>(sub)->GetParent()) &&
+             sub != parent);
+      return sub == parent;
+   }
+
 private:
    ParSubMesh(const ParMesh &parent, SubMesh::From from,
               const Array<int> &attributes);

mesh/submesh/submesh.hpp

@@ -225,6 +225,20 @@ public:
       return dynamic_cast<const SubMesh *>(m) != nullptr;
    }
 
+   /**
+   * @brief Check if Mesh @a sub is a SubMesh of Mesh @a parent.
+   *
+   * @param sub The potential submesh Mesh
+   * @param parent The potential parent Mesh
+   */
+   static bool IsSubMesh(const Mesh* sub, const Mesh* parent)
+   {
+      while (IsSubMesh(sub) &&
+             (sub = static_cast<const SubMesh *>(sub)->GetParent()) &&
+             sub != parent);
+      return sub == parent;
+   }
+
 private:
    /// Private constructor
    SubMesh(const Mesh &parent, From from, const Array<int> &attributes);

miniapps/CMakeLists.txt

@@ -35,6 +35,7 @@ add_subdirectory(multidomain)
 add_subdirectory(nurbs)
 add_subdirectory(parelag)
 add_subdirectory(performance)
+add_subdirectory(plasma)
 add_subdirectory(shifted)
 add_subdirectory(solvers)
 add_subdirectory(spde)

miniapps/electromagnetics/CMakeLists.txt

@@ -43,19 +43,39 @@ endif()
   
   # Add the corresponding tests to the "test" target
   if (MFEM_ENABLE_TESTING)
-    add_test(NAME tesla_np=4
+    add_test(NAME tesla_1_np=${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
       $<TARGET_FILE:tesla> -no-vis -maxit 2 -cr "0 0 -0.2 0 0 0.2 0.2 0.4 1"
       ${MPIEXEC_POSTFLAGS})
-
-    add_test(NAME volta_np=4
+    add_test(NAME tesla_2_np=${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:volta> -no-vis -maxit 2 -dbcs 1 -dbcg -ds "0.0 0.0 0.0 0.2 8.0"
+      $<TARGET_FILE:tesla>
+      -no-vis -maxit 2 -m ../../data/inline-hex.mesh -ubbc "0 0 1"
       ${MPIEXEC_POSTFLAGS})
 
-    add_test(NAME joule_np=4
+    add_test(NAME volta_1_np=${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:volta>
+      -no-vis -maxit 2 -dbcs 1 -dbcg -ds "0.0 0.0 0.0 0.2 8.0"
+      ${MPIEXEC_POSTFLAGS})
+    add_test(NAME volta_2_np=${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:volta>
+      -no-vis -maxit 2 -m ../../data/square-disc.mesh -dbcs "1 2 3 4 5 6 7 8"
+      -dbcv "0 0 0 0 1 1 1 1"
+      ${MPIEXEC_POSTFLAGS})
+    add_test(NAME volta_3_np=${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:volta>
+      -no-vis -maxit 2 -m ../../data/inline-hex.mesh -dbcs "1 6" -dbcv "0 1"
+      ${MPIEXEC_POSTFLAGS})
+
+    add_test(NAME joule_np=${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
       $<TARGET_FILE:joule>
@@ -63,12 +83,41 @@ endif()
       ${MPIEXEC_POSTFLAGS})
 
     if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME maxwell_np=4
+      add_test(NAME maxwell_np=${MFEM_MPI_NP}
         COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
         ${MPIEXEC_PREFLAGS}
         $<TARGET_FILE:maxwell>
         -no-vis -abcs "-1" -dp "-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5"
         ${MPIEXEC_POSTFLAGS})
     endif()
+
+    if (MFEM_USE_GSLIB)
+      add_test(NAME lorentz_1_np=${MFEM_MPI_NP}
+        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+        ${MPIEXEC_PREFLAGS}
+        $<TARGET_FILE:lorentz>
+        -no-vis -er Volta-AMR-Parallel -ec 2 -npt 100 -xmin "0.0 0.0 0.0"
+        -xmax "1.0 1.0 1.0" -pmin "1 0 0" -pmax "1 0 0" -rdf 0 -vt 0 -nt 100
+        ${MPIEXEC_POSTFLAGS})
+      # Setup dependency on volta_3_np=<np>
+      set_tests_properties(volta_3_np=${MFEM_MPI_NP}
+        PROPERTIES FIXTURES_SETUP Volta3)
+      set_tests_properties(lorentz_1_np=${MFEM_MPI_NP}
+        PROPERTIES FIXTURES_REQUIRED Volta3)
+
+      add_test(NAME lorentz_2_np=${MFEM_MPI_NP}
+        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+        ${MPIEXEC_PREFLAGS}
+        $<TARGET_FILE:lorentz>
+        -no-vis -br Tesla-AMR-Parallel -bc 2 -npt 10 -xmin "0.0 0.0 0.0"
+        -xmax "1.0 1.0 1.0" -pmin "0 0.1 0.05" -pmax "0 0.4 0.1" -nt 1000 -rdf 0
+        -vt 0
+        ${MPIEXEC_POSTFLAGS})
+      # Setup dependency on tesla_2_np=<np>
+      set_tests_properties(tesla_2_np=${MFEM_MPI_NP}
+        PROPERTIES FIXTURES_SETUP Tesla2)
+      set_tests_properties(lorentz_2_np=${MFEM_MPI_NP}
+        PROPERTIES FIXTURES_REQUIRED Tesla2)
+    endif()
   endif()
 endif()

miniapps/electromagnetics/makefile

@@ -117,10 +117,10 @@ joule-test-par: joule
 lorentz-test-par: lorentz-test-1 lorentz-test-2
 lorentz-test-1: lorentz volta-test-3
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-er Volta-AMR-Parallel -ec 2 -npt 100 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '1 0 0' -pmax '1 0 0' -rdf 0 -vt 0 -nt 100')
+	-er Volta-AMR-Parallel -ec 2 -npt 100 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '1 0 0' -pmax '1 0 0' -rdf 0 -vt 0 -nt 100)
 lorentz-test-2: lorentz tesla-test-2
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-br Tesla-AMR-Parallel -bc 2 -br Tesla-AMR-Parallel -npt 10 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '0 0.1 0.05' -pmax '0 0.4 0.1' -nt 1000 -rdf 0 -vt 0)
+	-br Tesla-AMR-Parallel -bc 2 -npt 10 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '0 0.1 0.05' -pmax '0 0.4 0.1' -nt 1000 -rdf 0 -vt 0)
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk