Merge branch 'master' into feature/partial_ess_bcs

fix CMake umpire build and install

.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

CHANGELOG

@@ -8,6 +8,22 @@
                                https://mfem.org
 
 
+Version 4.10 (development)
+==========================
+
+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.
+
+
 Version 4.9.1 (development)
 ===========================
 

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()

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/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/makefile

@@ -157,6 +157,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/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/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/fe/fe_base.cpp

@@ -1044,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)
@@ -1064,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_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

@@ -547,11 +547,22 @@ void MarkDofs(const Array<int> &dofs, Array<int> &mark_array)
 
 void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
                                            Array<int> &ess_vdofs,
-                                           int component) const
+                                           int component,
+                                           bool overwrite) const
 {
    Array<int> dofs;
-   ess_vdofs.SetSize(GetVSize());
-   ess_vdofs = 0;
+
+   if (overwrite)
+   {
+      ess_vdofs.SetSize(GetVSize());
+      ess_vdofs = 0;
+   }
+   else
+   {
+      MFEM_ASSERT(ess_vdofs.Size() == GetVSize(),
+                  "ess_vdofs size is not equal to FESpaces GetVSize().");
+   }
+
    for (int i = 0; i < GetNBE(); i++)
    {
       if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
@@ -663,6 +674,54 @@ void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
    MarkerToList(ess_tdofs, ess_tdof_list);
 }
 
+void FiniteElementSpace::GetEssentialVDofsFromComponent(
+   const Array<int> &bdr_attr_is_ess,
+   const Array2D<bool> &component,
+   Array<int> &ess_vdofs) const
+{
+   MFEM_ASSERT(component.NumCols() == vdim,
+               "Number of columns of component was not equal to FESpace vdim");
+   MFEM_ASSERT(component.NumRows() == bdr_attr_is_ess.Size(),
+               "Number of rows of component was not equal to bdr_attr_is_ess.Size()");
+
+   Array<int> bdr_attr_is_ess_single_comp;
+   bdr_attr_is_ess_single_comp.SetSize(bdr_attr_is_ess.Size());
+
+   for (int i = 0; i < vdim; i++)
+   {
+      const bool overwrite = (i == 0);
+      bdr_attr_is_ess_single_comp = 0;
+      for (int j = 0; j < bdr_attr_is_ess.Size(); j++)
+      {
+         if (bdr_attr_is_ess[j] && component(j, i))
+         {
+            bdr_attr_is_ess_single_comp[j] = bdr_attr_is_ess[j];
+         }
+      }
+      GetEssentialVDofs(bdr_attr_is_ess_single_comp, ess_vdofs, i, overwrite);
+   }
+}
+
+void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                              Array<int> &ess_tdof_list,
+                                              const Array2D<bool> &component)
+{
+   Array<int> ess_vdofs, ess_tdofs;
+
+   GetEssentialVDofsFromComponent(bdr_attr_is_ess, component, ess_vdofs);
+
+   const SparseMatrix *R = GetConformingRestriction();
+   if (!R)
+   {
+      ess_tdofs.MakeRef(ess_vdofs);
+   }
+   else
+   {
+      R->BooleanMult(ess_vdofs, ess_tdofs);
+   }
+   MarkerToList(ess_tdofs, ess_tdof_list);
+}
+
 void FiniteElementSpace::GetBoundaryTrueDofs(Array<int> &boundary_dofs,
                                              int component)
 {
@@ -3934,6 +3993,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");
@@ -3964,6 +4033,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

@@ -595,6 +595,29 @@ protected:
    virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes,
                                                const Array<int> *perm);
 
+   /** @brief Helper function to mark essential VDofs based on a component
+       matrix specifying which vector components are essential on each boundary
+       attribute.
+
+       This method loops over all vector dimensions (vdim) and calls
+       GetEssentialVDofs() for each component where the corresponding entry
+       in the @a component matrix is true.
+
+       @param[in] bdr_attr_is_ess Array marking which boundary attributes are
+                                  essential (1 for essential, 0 otherwise).
+       @param[in] component 2D boolean array of size (num_bdr_attributes x vdim)
+                           indicating which vector components are essential for
+                           each boundary attribute. component(j,i) == true means
+                           component i is essential on boundary attribute j.
+       @param[out] ess_vdofs Marker array for essential VDofs. On exit,
+                            ess_vdofs[i] != 0 if VDof i is essential.
+
+       @note This is a helper method used by GetEssentialTrueDofs() when
+             component-wise boundary conditions are specified. */
+   void GetEssentialVDofsFromComponent(const Array<int> &bdr_attr_is_ess,
+                                       const Array2D<bool> &component,
+                                       Array<int> &ess_vdofs) const;
+
 public:
 
 
@@ -839,7 +862,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;
 
@@ -1323,12 +1346,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;
@@ -1344,11 +1379,17 @@ public:
 
    /** @brief Mark degrees of freedom associated with boundary elements with
        the specified boundary attributes (marked in 'bdr_attr_is_ess').
+
        For spaces with 'vdim' > 1, the 'component' parameter can be used
-       to restricts the marked vDOFs to the specified component. */
+       to restricts the marked vDOFs to the specified component.
+       If overwrite is set to false then values in ess_vdofs are preserved
+       and not reset which allows the accumulation of multiple DOFs into a single array.
+       However, the assumption here is that ess_vdofs is set to
+       the correct size already.*/
    virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
                                   Array<int> &ess_vdofs,
-                                  int component = -1) const;
+                                  int component = -1,
+                                  bool overwrite = true) const;
 
    /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess.
@@ -1358,6 +1399,20 @@ public:
                                      Array<int> &ess_tdof_list,
                                      int component = -1) const;
 
+
+   /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       boundary attributes marked in the array bdr_attr_is_ess.
+
+       For spaces with 'vdim' > 1, the 'component' array can be used
+       to restricts the marked tDOFs per boundary to the specified components.
+       If vdim > 1 then one can specify per boundary attribute which components
+       on a boundary are essential by assigning a value of true to its location
+       in the component array.
+       The component has dimensions number of boundary attributes x vdim. */
+   virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                     Array<int> &ess_tdof_list,
+                                     const Array2D<bool> &component);
+
    /** @brief Get a list of all boundary true dofs, @a boundary_dofs. For spaces
        with 'vdim' > 1, the 'component' parameter can be used to restricts the
        marked tDOFs to the specified component. Equivalent to

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;
@@ -3033,6 +3033,7 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
                                           Array<int> &dof_attr)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> vdofs;
    Vector vals;
 
@@ -3064,6 +3065,7 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
 
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> dof_attr;
    ProjectDiscCoefficient(coeff, dof_attr);
 }
@@ -3073,6 +3075,10 @@ 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);
@@ -5758,4 +5768,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
    return std::make_pair(global_max_lower, global_max_upper);
 }
 
-}
+}

fem/gridfunc.hpp

@@ -150,11 +150,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
@@ -532,6 +534,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 +661,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 +1930,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 +1938,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.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/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/intrules.hpp

@@ -125,18 +125,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 +171,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 +189,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/nonlinearform.cpp

@@ -53,6 +53,22 @@ void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
    }
 }
 
+void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
+                                   const Array2D<bool> &bdr_component,
+                                   Vector *rhs)
+{
+   // virtual call, works in parallel too
+   fes->GetEssentialTrueDofs(bdr_attr_is_ess, ess_tdof_list, bdr_component);
+
+   if (rhs)
+   {
+      for (int i = 0; i < ess_tdof_list.Size(); i++)
+      {
+         (*rhs)(ess_tdof_list[i]) = 0.0;
+      }
+   }
+}
+
 void NonlinearForm::SetEssentialVDofs(const Array<int> &ess_vdofs_list)
 {
    if (!P)

fem/nonlinearform.hpp

@@ -165,14 +165,39 @@ public:
    const Array<NonlinearFormIntegrator*> &GetBdrFaceIntegrators() const
    { return bfnfi; }
 
-   /// Specify essential boundary conditions.
-   /** This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
+   /** @brief Specify essential boundary conditions.
+
+       This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
        the result internally for use by other methods. If the @a rhs pointer is
        not NULL, its essential true dofs will be set to zero. This makes it
        "compatible" with the output vectors from the Mult() method which also
-       have zero entries at the essential true dofs. */
+       have zero entries at the essential true dofs.
+
+       @note The values in the essential vdofs have to come from the initial guess.
+       */
    void SetEssentialBC(const Array<int> &bdr_attr_is_ess, Vector *rhs = NULL);
 
+   /** @brief Specify essential boundary conditions.
+
+       For spaces with 'vdim' > 1, the 'bdr_component' array can be used
+       to restricts the marked tDOFs per boundary to the specified components.
+       If vdim > 1 then one can specify per boundary attribute which components
+       on a boundary are essential by assigning a value of true to its location
+       in the bdr_component array.
+       The bdr_component has dimensions number of boundary attributes x vdim
+
+       This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
+       the result internally for use by other methods. If the @a rhs pointer is
+       not NULL, its essential true dofs will be set to zero. This makes it
+       "compatible" with the output vectors from the Mult() method which also
+       have zero entries at the essential true dofs.
+
+       @note The values in the essential vdofs have to come from the initial guess.
+       */
+   void SetEssentialBC(const Array<int> &bdr_attr_is_ess,
+                       const Array2D<bool> &bdr_component,
+                       Vector *rhs);
+
    /// Specify essential boundary conditions.
    /** Use either SetEssentialBC() or SetEssentialTrueDofs() if possible. */
    void SetEssentialVDofs(const Array<int> &ess_vdofs_list);

fem/pfespace.cpp

@@ -1122,9 +1122,11 @@ void ParFiniteElementSpace::Synchronize(Array<int> &ldof_marker) const
 
 void ParFiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
                                               Array<int> &ess_dofs,
-                                              int component) const
+                                              int component,
+                                              bool overwrite) const
 {
-   FiniteElementSpace::GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component);
+   FiniteElementSpace::GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component,
+                                         overwrite);
 
    // Make sure that processors without boundary elements mark
    // their boundary dofs (if they have any).
@@ -1184,6 +1186,39 @@ void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
    MarkerToList(true_ess_dofs, ess_tdof_list);
 }
 
+void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
+                                                 &bdr_attr_is_ess,
+                                                 Array<int> &ess_tdof_list,
+                                                 const Array2D<bool> &component)
+{
+   MFEM_VERIFY(!IsVariableOrderH1(),
+               "Variable order H1 spaces are currently not supported with this feature");
+
+   Array<int> ess_dofs, true_ess_dofs;
+
+   GetEssentialVDofsFromComponent(bdr_attr_is_ess, component, ess_dofs);
+   GetRestrictionMatrix()->BooleanMult(ess_dofs, true_ess_dofs);
+
+#ifdef MFEM_DEBUG
+   // Verify that in boolean arithmetic: P^T ess_dofs = R ess_dofs.
+   Array<int> true_ess_dofs2(true_ess_dofs.Size());
+   HypreParMatrix *Pt = Dof_TrueDof_Matrix()->Transpose();
+   const int *ess_dofs_data = ess_dofs.HostRead();
+   Pt->BooleanMult(1, ess_dofs_data, 0, true_ess_dofs2);
+   delete Pt;
+   int counter = 0;
+   const int *ted = true_ess_dofs.HostRead();
+   for (int i = 0; i < true_ess_dofs.Size(); i++)
+   {
+      if (bool(ted[i]) != bool(true_ess_dofs2[i])) { counter++; }
+   }
+   MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter
+               << ", rank = " << MyRank);
+#endif
+
+   MarkerToList(true_ess_dofs, ess_tdof_list);
+}
+
 void ParFiniteElementSpace::GetEssentialTrueDofsVar(const Array<int>
                                                     &bdr_attr_is_ess,
                                                     const Array<int> &ess_dofs,

fem/pfespace.hpp

@@ -420,10 +420,18 @@ public:
        "partially conforming") space. */
    void Synchronize(Array<int> &ldof_marker) const;
 
-   /// Determine the boundary degrees of freedom
-   void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
-                          Array<int> &ess_dofs,
-                          int component = -1) const override;
+   /** @brief Mark degrees of freedom associated with boundary elements with
+       the specified boundary attributes (marked in 'bdr_attr_is_ess').
+
+       For spaces with 'vdim' > 1, the 'component' parameter can be used
+       to restricts the marked vDOFs to the specified component.
+       If overwrite is set to false then values in ess_vdofs are preserved
+       and not reset. However, the assumption here is that ess_vdofs is set to
+       the correct size already.*/
+   virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
+                                  Array<int> &ess_dofs,
+                                  int component = -1,
+                                  bool overwrite = true) const override;
 
    /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess. */
@@ -448,6 +456,19 @@ public:
    void GetExteriorTrueDofs(Array<int> &ext_tdof_list,
                             int component = -1) const override;
 
+   /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       boundary attributes marked in the array bdr_attr_is_ess.
+
+       For spaces with 'vdim' > 1, the 'component' array can be used
+       to restricts the marked tDOFs per boundary to the specified components.
+       If vdim > 1 then one can specify per boundary attribute which components
+       on a boundary are essential by assigning a value of true to its location
+       in the component array.
+       The component has dimensions number of boundary attributes x vdim. */
+   virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                     Array<int> &ess_tdof_list,
+                                     const Array2D<bool> &component) override;
+
    /** If the given ldof is owned by the current processor, return its local
        tdof number, otherwise return -1 */
    int GetLocalTDofNumber(int ldof) const;

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)
@@ -717,6 +719,7 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 
 void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    // 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

@@ -280,8 +280,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/restriction.cpp

@@ -1398,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
 }

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/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/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

@@ -126,11 +126,11 @@ EXAMPLE_TEST_DIRS := examples
 MINIAPP_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_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 shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
 
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
 

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_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/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;
 }

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/hdiv-linear-solver/change_basis.cpp

@@ -22,7 +22,7 @@ void ComputeInverse(const Array<real_t> &A, Array<real_t> &Ainv)
 {
    Array<real_t> A2 = A;
    const int n2 = A.Size();
-   const int n = static_cast<const int>(sqrt(n2));
+   const int n = static_cast<int>(sqrt(n2));
    Array<int> ipiv(n);
    LUFactors lu(A2.GetData(), ipiv.GetData());
    lu.Factor(n);
@@ -58,7 +58,7 @@ void SubcellIntegrals(int n, const Poly_1D::Basis &basis, Array<real_t> &B)
 
 void Transpose(const Array<real_t> &B, Array<real_t> &Bt)
 {
-   const int n = static_cast<const int>(sqrt(B.Size()));
+   const int n = static_cast<int>(sqrt(B.Size()));
    Bt.SetSize(n*n);
    for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) { Bt[i+j*n] = B[j+i*n]; }
 }

miniapps/plasma/CMakeLists.txt

@@ -0,0 +1,27 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+if (MFEM_USE_MPI)
+  list(APPEND PLASMA_COMMON_SOURCES)
+
+  list(APPEND PLASMA_COMMON_HEADERS
+  plasma.hpp)
+
+  convert_filenames_to_full_paths(PLASMA_COMMON_SOURCES)
+  convert_filenames_to_full_paths(PLASMA_COMMON_HEADERS)
+
+  set(PLASMA_COMMON_FILES
+      EXTRA_SOURCES ${PLASMA_COMMON_SOURCES}
+      EXTRA_HEADERS ${PLASMA_COMMON_HEADERS})
+
+endif()
+
+add_subdirectory(pic)

miniapps/plasma/makefile

@@ -0,0 +1,93 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS = 
+PAR_MINIAPPS = 
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+PLASMA_SUBDIRS = pic
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all lib-common clean clean-build clean-exec
+.PRECIOUS: %.o
+
+COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
+
+# If MFEM_SHARED is set, add the ../common rpath
+COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+COMMON_O=
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(MINIAPPS) subdirs
+
+.PHONY: subdirs $(PLASMA_SUBDIRS)
+subdirs: $(PLASMA_SUBDIRS)
+$(PLASMA_SUBDIRS): lib-common
+	$(MAKE) -C $(BLD)$(@)
+
+# Rules for building the miniapps
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_O) $(COMMON_LIB) \
+        $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+# Rule for building lib-common
+lib-common:
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Specific execution options
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+ALL_CLEAN_SUBDIRS = $(addsuffix /clean,$(PLASMA_SUBDIRS))
+.PHONY: $(ALL_CLEAN_SUBDIRS)
+$(ALL_CLEAN_SUBDIRS):
+	$(MAKE) -C $(BLD)$(@D) $(@F)
+
+clean: clean-build clean-exec
+
+clean-build: $(addsuffix /clean,$(PLASMA_SUBDIRS))
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:

miniapps/plasma/pic/CMakeLists.txt

@@ -0,0 +1,28 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+if (MFEM_USE_MPI AND MFEM_USE_GSLIB)
+  add_mfem_miniapp(electrostatic-pic
+    MAIN electrostatic-pic.cpp
+    EXTRA_HEADERS ${MFEM_MINIAPPS_COMMON_HEADERS}
+    LIBRARIES mfem-common)
+
+  # Add the corresponding tests to the "test" target
+  if (MFEM_ENABLE_TESTING)
+    add_test(NAME electrostatic-pic_np=${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:electrostatic-pic> -rdi 2 -npt 40960 -k 0.2855993321 -a 0.05
+        -nt 200 -nx 16 -ny 16 -O 1 -q 0.01181640625 -m 0.01181640625 -oci 1000
+        -dt 0.1
+      ${MPIEXEC_POSTFLAGS})
+  endif()
+endif()

miniapps/plasma/pic/electrostatic-pic.cpp

@@ -0,0 +1,788 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+//           -----------------------------------------------------
+//           Particle-In-Cell (PIC) Simulation (2D/3D)
+//           -----------------------------------------------------
+//
+// This miniapp performs a Particle-In-Cell simulation (supports 2D or 3D
+// spatial dimensions) of multiple charged particles subject to electric
+// field forces.
+//
+//                           dp/dt = q E
+//
+// The method used is explicit time integration with a leap-frog scheme.
+//
+// The electric field is computed from the particle charge distribution using
+// a Poisson solver. The particle trajectories are computed within a periodic
+// domain (2D or 3D).
+//
+// Solution process (per timestep, repeating steps 1-6):
+//   (1) Deposit charge from particles to grid via Dirac delta function
+//       to form the RHS of the Poisson equation
+//   (2) Solve Poisson equation (-Δφ = ρ - ρ_0) to compute potential φ, where
+//       ρ_0 is a constant neutralizing term that enforces global charge
+//       neutrality.
+//   (3) Compute electric field E = -∇φ from the potential
+//   (4) Interpolate E-field to particle positions
+//   (5) Push particles using leap-frog scheme (update momentum and position)
+//   (6) Redistribute particles across processors
+//
+// Compile with: make electrostatic-pic
+//
+// Sample runs:
+//
+//   2D2V Linear Landau damping test case (Ricketson & Hu, 2025):
+//      mpirun -n 4 ./electrostatic-pic -rdi 1 -npt 409600 -k 0.2855993321 -a 0.05 -nt 200 -nx 32 -ny 32 -O 1 -q 0.001181640625 -m 0.001181640625 -oci 1000 -dt 0.1
+//   3D3V Linear Landau damping test case (Zheng et al., 2025):
+//    * mpirun -n 128 ./electrostatic-pic -dim 3 -rdi 1 -npt 40960000 -k 0.5 -a 0.01 -nt 100 -nx 32 -ny 32 -nz 32 -O 1 -q 0.00004844730731 -m 0.00004844730731 -oci 1000 -dt 0.02 -no-vis
+
+#include "mfem.hpp"
+#include "../../../general/text.hpp"
+#include "../../common/fem_extras.hpp"
+#include "../../common/particles_extras.hpp"
+#include "../../common/pfem_extras.hpp"
+
+#include <ctime>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <random>
+#include <string>
+#include <vector>
+
+#define EPSILON 1  // ε_0
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::common;
+
+struct PICContext
+{
+   int dim = 2;     ///< Spatial dimension.
+   int order = 1;   ///< FE order for spatial discretization.
+   int nx = 100;    ///< Number of grid cells in x-direction.
+   int ny = 100;    ///< Number of grid cells in y-direction.
+   int nz = 100;    ///< Number of grid cells in z-direction.
+   real_t L = 1.0;  ///< Domain length.
+
+   int ordering = 1;  ///< Ordering of particles.
+   int npt = 1000;    ///< Number of particles.
+   real_t q = 1.0;    ///< Particle charge.
+   real_t m = 1.0;    ///< Particle mass.
+
+   real_t k = 1.0;      ///< Wave number (Landau damping init).
+   real_t alpha = 0.1;  ///< Perturbation amplitude (Landau damping init).
+
+   real_t dt = 1e-2;  ///< Time step size.
+
+   int nt = 1000;            ///< Number of time steps to run.
+   int redist_interval = 5;  ///< Redistribution and update E_gf interval.
+   int output_csv_interval = 1000;  ///< Interval for outputting CSV data files.
+
+   bool visualization = true;  ///< Enable visualization.
+   int visport = 19916;        ///< Port number for visualization server.
+   bool reproduce = true;      ///< Enable reproducible results.
+} ctx;
+
+/** This class implements explicit time integration for charged particles
+    in an electric field using ParticleSet. */
+class ParticleMover
+{
+public:
+   enum Fields
+   {
+      MASS,    // vdim = 1
+      CHARGE,  // vdim = 1
+      MOM,     // vdim = dim
+      EFIELD   // vdim = dim
+   };
+
+protected:
+   /// Pointers to E field GridFunctions
+   ParGridFunction* E_gf;
+
+   /// FindPointsGSLIB object for E field mesh
+   FindPointsGSLIB& E_finder;
+
+   /// ParticleSet of charged particles
+   std::unique_ptr<ParticleSet> charged_particles;
+
+   /// Temporary vectors for particle computation
+   mutable Vector pm_, pp_;
+
+public:
+   ParticleMover(MPI_Comm comm, ParGridFunction* E_gf_,
+                 FindPointsGSLIB& E_finder_, int num_particles,
+                 Ordering::Type pdata_ordering);
+
+   /// Initialize charged particles with given parameters
+   void InitializeChargedParticles(const real_t& k, const real_t& alpha,
+                                   real_t m, real_t q, real_t L,
+                                   bool reproduce = false);
+
+   /// Find Particles in mesh corresponding to E and field
+   void FindParticles();
+
+   /// Advance particles one time step using Boris algorithm
+   void Step(real_t& t, real_t dt, real_t L, bool first_step = false);
+
+   /// Redistribute particles across processors
+   void Redistribute();
+
+   /// Get reference to ParticleSet
+   ParticleSet& GetParticles() { return *charged_particles; }
+
+   /// Compute (global) kinetic energy from particles
+   /** Optionally, advance the particle momenta by time step @a dt. */
+   real_t ComputeKineticEnergy(real_t dt = 0.) const;
+};
+
+/** Field solver responsible for updating the electrostatic potential and field
+    from the particle charge density. Assembles and solves the periodic Poisson
+    problem, computes the electric field via a discrete gradient operator, and
+    provides utilities for field diagnostics (e.g. global field energy). */
+class FieldSolver
+{
+private:
+   real_t domain_volume;
+   real_t neutralizing_const;
+   ParLinearForm* precomputed_neutralizing_lf = nullptr;
+   bool precompute_neutralizing_const = false;
+   // Diffusion matrix
+   HypreParMatrix* diffusion_matrix;
+   // Gradient operator for computing E = -∇φ
+   ParDiscreteLinearOperator* grad_interpolator;
+   FindPointsGSLIB& E_finder;
+   ParLinearForm b;
+
+protected:
+   /** Compute neutralizing constant and initialize with the constant.
+       Returns a reference to the precomputed neutralizing ParLinearForm. */
+   const ParLinearForm& ComputeNeutralizingRHS(ParFiniteElementSpace* pfes,
+                                               const ParticleVector& Q,
+                                               MPI_Comm comm);
+
+   /** Deposit charge from particles into a ParLinearForm (RHS b).
+       b_i = sum_p q_p * φ_i(x_p) */
+   void DepositCharge(ParFiniteElementSpace* pfes, const ParticleVector& Q);
+
+public:
+   FieldSolver(ParFiniteElementSpace* phi_fes, ParFiniteElementSpace* E_fes,
+               FindPointsGSLIB& E_finder_,
+               bool precompute_neutralizing_const_ = false);
+
+   ~FieldSolver();
+
+   /** Update the phi_gf grid function from the particles.
+       Solve periodic Poisson: diffusion_matrix * phi = (rho - <rho>)
+       with zero-mean enforcement via OrthoSolver. */
+   void UpdatePhiGridFunction(ParticleSet& particles, ParGridFunction& phi_gf);
+
+   /** Update E_gf grid function from phi_gf grid function.
+       Compute the gradient: E = -∇φ. */
+   void UpdateEGridFunction(ParGridFunction& phi_gf, ParGridFunction& E_gf);
+
+   /// Compute (global) field energy: 0.5 * ∫ ||E||^2 dx
+   real_t ComputeFieldEnergy(const ParGridFunction& E_gf) const;
+};
+
+/// Prints the program's logo to the given output stream
+void display_banner(ostream& os);
+
+int main(int argc, char* argv[])
+{
+   Mpi::Init(argc, argv);
+   int num_ranks = Mpi::WorldSize();
+   int rank = Mpi::WorldRank();
+   Hypre::Init();
+
+   if (Mpi::Root()) { display_banner(cout); }
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ctx.dim, "-dim", "--dimension",
+                  "Spatial dimension (2 or 3)");
+   args.AddOption(&ctx.order, "-O", "--order",
+                  "Finite element polynomial degree");
+   args.AddOption(&ctx.nx, "-nx", "--num-x",
+                  "Number of elements in the x direction.");
+   args.AddOption(&ctx.ny, "-ny", "--num-y",
+                  "Number of elements in the y direction.");
+   args.AddOption(&ctx.nz, "-nz", "--num-z",
+                  "Number of elements in the z direction.");
+   args.AddOption(&ctx.q, "-q", "--charge", "Particle charge.");
+   args.AddOption(&ctx.m, "-m", "--mass", "Particle mass.");
+   args.AddOption(&ctx.dt, "-dt", "--time-step", "Time Step.");
+   args.AddOption(&ctx.nt, "-nt", "--num-timesteps", "Number of timesteps.");
+   args.AddOption(&ctx.npt, "-npt", "--num-particles",
+                  "Total number of particles.");
+   args.AddOption(&ctx.k, "-k", "--k", "Wave number for initial distribution.");
+   args.AddOption(&ctx.alpha, "-a", "--alpha",
+                  "Perturbation amplitude for initial distribution.");
+   args.AddOption(&ctx.ordering, "-o", "--ordering",
+                  "Ordering of particle data. 0 = byNODES, 1 = byVDIM.");
+   args.AddOption(&ctx.redist_interval, "-rdi", "--redist-interval",
+                  "Redistribution and update E_gf interval. Disabled if < 0.");
+   args.AddOption(&ctx.output_csv_interval, "-oci", "--output-csv-interval",
+                  "Output CSV interval. Disabled if < 0.");
+   args.AddOption(&ctx.visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&ctx.visport, "-p", "--send-port", "Socket for GLVis.");
+   args.AddOption(&ctx.reproduce, "-rep", "--reproduce", "-no-rep",
+                  "--no-reproduce",
+                  "Enable or disable reproducible random seed.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (Mpi::Root()) { args.PrintUsage(cout); }
+      return 1;
+   }
+   if (Mpi::Root()) { args.PrintOptions(cout); }
+
+   // Assert that dimension is 2 or 3
+   MFEM_VERIFY(ctx.dim == 2 || ctx.dim == 3,
+               "Dimension must be 2 or 3, got " << ctx.dim);
+   MFEM_VERIFY(ctx.alpha >= -1.0 && ctx.alpha < 1.0,
+               "Alpha should be in range [-1, 1).");
+   MFEM_VERIFY(ctx.k > 0.0,
+               "k must be nonzero for displacement initialization.");
+
+   ctx.L = 2.0 * M_PI / ctx.k;
+
+   // 1. make a Cartesian Mesh (2D or 3D)
+   Mesh serial_mesh;
+   std::vector<Vector> translations;
+
+   if (ctx.dim == 2)
+   {
+      serial_mesh = Mesh(Mesh::MakeCartesian2D(
+                            ctx.nx, ctx.ny, Element::QUADRILATERAL, false, ctx.L, ctx.L));
+      translations = {Vector({ctx.L, 0.0}), Vector({0.0, ctx.L})};
+   }
+   else  // ctx.dim == 3
+   {
+      serial_mesh = Mesh(Mesh::MakeCartesian3D(
+                            ctx.nx, ctx.ny, ctx.nz, Element::HEXAHEDRON, ctx.L, ctx.L, ctx.L));
+      translations = {Vector({ctx.L, 0.0, 0.0}), Vector({0.0, ctx.L, 0.0}),
+                      Vector({0.0, 0.0, ctx.L})
+                     };
+   }
+
+   Mesh periodic_mesh(Mesh::MakePeriodic(
+                         serial_mesh, serial_mesh.CreatePeriodicVertexMapping(translations)));
+   // 2. Partition and distribute the mesh
+   ParMesh mesh(MPI_COMM_WORLD, periodic_mesh);
+   serial_mesh.Clear();    // the serial mesh is no longer needed
+   periodic_mesh.Clear();  // the periodic mesh is no longer needed
+
+   // 3. Build the interpolator of E field
+   mesh.EnsureNodes();
+   FindPointsGSLIB E_finder(mesh);
+
+   // 4. Define finite element spaces on the parallel mesh
+   H1_FECollection phi_fec(ctx.order, ctx.dim);
+   ParFiniteElementSpace phi_fespace(&mesh, &phi_fec);
+   ND_FECollection E_fec(ctx.order, ctx.dim);
+   ParFiniteElementSpace E_fespace(&mesh, &E_fec);
+
+   // 5. Initialize the grid functions for the electric field and potential
+   ParGridFunction phi_gf(&phi_fespace);
+   ParGridFunction E_gf(&E_fespace);
+   phi_gf = 0.0;  // Initialize phi_gf to zero
+   E_gf = 0.0;    // Initialize E_gf to zero
+
+   // 6. Construct the field solver
+   FieldSolver field_solver(&phi_fespace, &E_fespace, E_finder, true);
+
+   // 7. Initialize ParticleMover
+   Ordering::Type ordering_type =
+      ctx.ordering == 0 ? Ordering::byNODES : Ordering::byVDIM;
+   int num_particles =
+      ctx.npt / num_ranks + (rank < (ctx.npt % num_ranks) ? 1 : 0);
+   ParticleMover particle_mover(MPI_COMM_WORLD, &E_gf, E_finder, num_particles,
+                                ordering_type);
+   particle_mover.InitializeChargedParticles(ctx.k, ctx.alpha, ctx.m, ctx.q,
+                                             ctx.L, ctx.reproduce);
+
+   // 8. Start the main loop
+   real_t t = 0;
+   real_t dt = ctx.dt;
+
+   mfem::StopWatch sw;
+   sw.Start();
+   for (int step = 1; step <= ctx.nt; step++)
+   {
+      // Step the FieldSolver
+      if (ctx.redist_interval > 0 &&
+          (step % ctx.redist_interval == 0 || step == 1) &&
+          particle_mover.GetParticles().GetGlobalNParticles() > 0)
+      {
+         // Redistribute
+         particle_mover.Redistribute();
+
+         // Update phi_gf from particles
+         field_solver.UpdatePhiGridFunction(particle_mover.GetParticles(),
+                                            phi_gf);
+         // Update E_gf from phi_gf
+         field_solver.UpdateEGridFunction(phi_gf, E_gf);
+
+         // Visualize fields if requested
+         if (ctx.visualization)
+         {
+            static socketstream vis_e, vis_phi;
+            common::VisualizeField(vis_e, "localhost", ctx.visport, E_gf,
+                                   "E_field", 0, 0, 500, 500);
+            common::VisualizeField(vis_phi, "localhost", ctx.visport, phi_gf,
+                                   "Potential", 500, 0, 500, 500);
+         }
+      }
+
+      // Step the ParticleMover
+      particle_mover.Step(t, dt, ctx.L, step == 1);
+      if (Mpi::Root())
+      {
+         mfem::out << "Step: " << step << " | Time: " << t;
+         mfem::out << " | Time per step: " << sw.RealTime() / step;
+         mfem::out << endl;
+      }
+      // Output particle data to CSV
+      if (ctx.output_csv_interval > 0 &&
+          (step % ctx.output_csv_interval == 0 || step == 1))
+      {
+         std::string csv_prefix = "PIC_Part_";
+         Array<int> field_idx{2}, tag_idx;
+         std::string file_name =
+            csv_prefix + mfem::to_padded_string(step, 6) + ".csv";
+         particle_mover.GetParticles().PrintCSV(file_name.c_str(), field_idx,
+                                                tag_idx);
+      }
+
+      if (ctx.redist_interval > 0 &&
+          (step % ctx.redist_interval == 0 || step == 1) &&
+          particle_mover.GetParticles().GetGlobalNParticles() > 0)
+      {
+         // Compute energies
+         // Note that particle momenta are a half time step ahead of the field
+         // after particle_mover.Step(). Therefore they are returned to the
+         // time level of the field for calculation of kinetic energy.
+         real_t kinetic_energy = particle_mover.ComputeKineticEnergy(-dt/2.);
+         real_t field_energy = field_solver.ComputeFieldEnergy(E_gf);
+
+         // Output energies
+         if (Mpi::Root())
+         {
+            cout << "Kinetic energy: " << kinetic_energy << "\t"
+                 << "Field energy: " << field_energy << "\t"
+                 << "Total energy: " << kinetic_energy + field_energy
+                 << endl;
+         }
+         // Write energies to a CSV file
+         if (Mpi::Root())
+         {
+            std::ofstream energy_file("energy.csv", std::ios::app);
+            energy_file << setprecision(10) << kinetic_energy << ","
+                        << field_energy << "," << kinetic_energy + field_energy
+                        << "\n";
+         }
+      }
+   }
+}
+
+ParticleMover::ParticleMover(MPI_Comm comm, ParGridFunction* E_gf_,
+                             FindPointsGSLIB& E_finder_, int num_particles,
+                             Ordering::Type pdata_ordering)
+   : E_gf(E_gf_), E_finder(E_finder_)
+{
+   MFEM_ASSERT(E_gf, "Must pass an E field to ParticleMover.");
+
+   int dim = E_gf->ParFESpace()->GetMesh()->SpaceDimension();
+
+   pm_.SetSize(dim);
+   pp_.SetSize(dim);
+
+   // Create particle set: 2 scalars of mass and charge,
+   // 2 vectors of size space dim for momentum and e field
+   Array<int> field_vdims({1, 1, dim, dim});
+   charged_particles = std::make_unique<ParticleSet>(
+                          comm, num_particles, dim, field_vdims, 1, pdata_ordering);
+}
+
+void ParticleMover::InitializeChargedParticles(const real_t& k,
+                                               const real_t& alpha, real_t m,
+                                               real_t q, real_t L,
+                                               bool reproduce)
+{
+   int rank;
+   MPI_Comm_rank(charged_particles->GetComm(), &rank);
+   // use time-based seed for randomness
+   std::mt19937 gen(
+      reproduce ? rank : (rank + static_cast<unsigned int>(time(nullptr))));
+   std::uniform_real_distribution<> real_dist(0.0, 1.0);
+   std::normal_distribution<> norm_dist(0.0, 1.0);
+
+   int dim = charged_particles->Coords().GetVDim();
+
+   ParticleVector& X = charged_particles->Coords();
+   ParticleVector& P = charged_particles->Field(ParticleMover::MOM);
+   ParticleVector& M = charged_particles->Field(ParticleMover::MASS);
+   ParticleVector& Q = charged_particles->Field(ParticleMover::CHARGE);
+
+   for (int i = 0; i < charged_particles->GetNParticles(); i++)
+   {
+      // Initialize momentum
+      for (int d = 0; d < dim; d++) { P(i, d) = m * norm_dist(gen); }
+
+      // Uniform positions (no accept-reject)
+      for (int d = 0; d < dim; d++) { X(i, d) = real_dist(gen) * L; }
+
+      // Displacement along x for perturbation ~ cos(k x)
+      for (int d = 0; d < dim; d++)
+      {
+         real_t x = X(i, d);
+         x -= (alpha / k) * std::sin(k * x);
+
+         // periodic wrap to [0, L)
+         x = std::fmod(x, L);
+         if (x < 0) { x += L; }
+
+         X(i, d) = x;
+      }
+
+      // Initialize mass + charge
+      M(i) = m;
+      Q(i) = q;
+   }
+   FindParticles();
+}
+
+void ParticleMover::FindParticles()
+{
+   E_finder.FindPoints(charged_particles->Coords());
+}
+
+void ParticleMover::Step(real_t& t, real_t dt, real_t L, bool first_step)
+{
+   // Update E field at particles
+   ParticleVector& E = charged_particles->Field(EFIELD);
+   E_finder.Interpolate(*E_gf, E, E.GetOrdering());
+
+   // Extract particle data
+   ParticleVector& X = charged_particles->Coords();
+   ParticleVector& P = charged_particles->Field(MOM);
+   ParticleVector& M = charged_particles->Field(MASS);
+   ParticleVector& Q = charged_particles->Field(CHARGE);
+
+   // Accelerate the particles by the electric field
+   const int npt = charged_particles->GetNParticles();
+   const int dim = X.GetVDim();
+
+   for (int particle = 0; particle < npt; ++particle)
+   {
+      for (int d = 0; d < dim; ++d)
+      {
+         P(particle, d) +=
+            (first_step ? dt / 2.0 : dt) * Q(particle) * E(particle, d);
+      }
+   }
+
+   // Periodic boundary: wrap coordinates to [0, L)
+   for (int particle = 0; particle < npt; ++particle)
+   {
+      for (int d = 0; d < dim; ++d)
+      {
+         X(particle, d) += dt / M(particle) * P(particle, d);
+         while (X(particle, d) > L) { X(particle, d) -= L; }
+         while (X(particle, d) < 0.0) { X(particle, d) += L; }
+      }
+   }
+
+   FindParticles();
+
+   // Update time
+   t += dt;
+}
+
+void ParticleMover::Redistribute()
+{
+   charged_particles->Redistribute(E_finder.GetProc());
+   FindParticles();
+}
+
+real_t ParticleMover::ComputeKineticEnergy(real_t dt) const
+{
+   const ParticleVector& P = charged_particles->Field(MOM);
+   const ParticleVector& M = charged_particles->Field(MASS);
+   const ParticleVector& Q = charged_particles->Field(CHARGE);
+   const ParticleVector& E = charged_particles->Field(EFIELD);
+
+   // Note the electric field is not reinterpolated here and the last
+   // update from Step() is used directly.
+
+   real_t kinetic_energy = 0.0;
+   for (int p = 0; p < charged_particles->GetNParticles(); ++p)
+   {
+      real_t p_square_p = 0.0;
+      for (int d = 0; d < P.GetVDim(); ++d)
+      {
+         const real_t P_m = P(p, d) + dt * Q(p) * E(p, d);
+         p_square_p += P_m * P_m;
+      }
+      kinetic_energy += 0.5 * p_square_p / M(p);
+   }
+
+   real_t global_kinetic_energy = 0.0;
+   MPI_Allreduce(&kinetic_energy, &global_kinetic_energy, 1, MPI_DOUBLE,
+                 MPI_SUM, charged_particles->GetComm());
+   return global_kinetic_energy;
+}
+
+FieldSolver::FieldSolver(ParFiniteElementSpace* phi_fes,
+                         ParFiniteElementSpace* E_fes,
+                         FindPointsGSLIB& E_finder_,
+                         bool precompute_neutralizing_const_)
+   : precompute_neutralizing_const(precompute_neutralizing_const_),
+     E_finder(E_finder_),
+     b(phi_fes)
+{
+   // compute domain volume
+   ParMesh* pmesh = phi_fes->GetParMesh();
+   real_t local_domain_volume = 0.0;
+   for (int i = 0; i < pmesh->GetNE(); i++)
+   {
+      local_domain_volume += pmesh->GetElementVolume(i);
+   }
+   MPI_Allreduce(&local_domain_volume, &domain_volume, 1, MPI_DOUBLE, MPI_SUM,
+                 phi_fes->GetParMesh()->GetComm());
+
+   {
+      // Par bilinear form for the gradgrad matrix
+      ParBilinearForm dm(phi_fes);
+      ConstantCoefficient epsilon(EPSILON);  // ε_0
+      dm.AddDomainIntegrator(
+         new DiffusionIntegrator(epsilon));  // ∫ ∇φ_i · ∇φ_j
+
+      dm.Assemble();
+      dm.Finalize();
+
+      diffusion_matrix = dm.ParallelAssemble();  // global gradgrad matrix
+   }
+
+   {
+      // Compute E = -∇φ using DiscreteLinearOperator
+      grad_interpolator = new ParDiscreteLinearOperator(phi_fes, E_fes);
+      grad_interpolator->AddDomainInterpolator(new GradientInterpolator);
+      grad_interpolator->Assemble();
+   }
+}
+
+FieldSolver::~FieldSolver()
+{
+   delete diffusion_matrix;
+   delete precomputed_neutralizing_lf;
+   delete grad_interpolator;
+}
+
+const ParLinearForm& FieldSolver::ComputeNeutralizingRHS(
+   ParFiniteElementSpace* pfes, const ParticleVector& Q, MPI_Comm comm)
+{
+   int npt = Q.Size();
+   // Get E_finder references
+   const Array<unsigned int>& code = E_finder.GetCode();
+
+   if (!precompute_neutralizing_const || precomputed_neutralizing_lf == nullptr)
+   {
+      // compute neutralizing constant
+      real_t local_sum = 0.0;
+      for (int p = 0; p < npt; ++p)
+      {
+         // Skip particles not successfully found
+         MFEM_ASSERT(code[p] != 2, "Particle " << p << " not found.");
+         local_sum += Q(p);
+      }
+
+      real_t global_sum = 0.0;
+      MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, comm);
+
+      neutralizing_const = -global_sum / domain_volume;
+      if (Mpi::Root())
+      {
+         cout << "Total charge: " << global_sum
+              << ", Domain volume: " << domain_volume
+              << ", Neutralizing constant: " << neutralizing_const << endl;
+         if (precompute_neutralizing_const)
+         {
+            cout << "Further updates will use this precomputed neutralizing "
+                 "constant."
+                 << endl;
+         }
+      }
+      delete precomputed_neutralizing_lf;
+      precomputed_neutralizing_lf = new ParLinearForm(pfes);
+      *precomputed_neutralizing_lf = 0.0;
+      ConstantCoefficient neutralizing_coeff(neutralizing_const);
+      precomputed_neutralizing_lf->AddDomainIntegrator(
+         new DomainLFIntegrator(neutralizing_coeff));
+      precomputed_neutralizing_lf->Assemble();
+   }
+   return *precomputed_neutralizing_lf;
+}
+
+void FieldSolver::DepositCharge(ParFiniteElementSpace* pfes,
+                                const ParticleVector& Q)
+{
+   int npt = Q.Size();
+   ParMesh* pmesh = pfes->GetParMesh();
+   int dim = pmesh->SpaceDimension();
+   int curr_rank;
+   MPI_Comm_rank(pmesh->GetComm(), &curr_rank);
+
+   // Get E_finder references
+   // 0: inside, 1: boundary, 2: not found
+   const Array<unsigned int>& code = E_finder.GetCode();
+   const Array<unsigned int>& proc = E_finder.GetProc();  // owning MPI rank
+   const Array<unsigned int>& elem = E_finder.GetElem();  // local element id
+   const Vector& rref = E_finder.GetReferencePosition();  // (r,s,t) byVDIM
+
+   Array<int> dofs;
+
+   for (int p = 0; p < npt; ++p)
+   {
+      // Skip particles not successfully found
+      MFEM_ASSERT(code[p] != 2, "Particle " << p << " not found.");
+
+      // Assert particle is on the current rank
+      MFEM_ASSERT((int)proc[p] == curr_rank,
+                  "Particle " << p << " found in element owned by rank "
+                  << proc[p] << " but current rank is " << curr_rank
+                  << "." << endl
+                  << "You must call redistribute everytime before "
+                  "updating the density grid function.");
+      const int e = elem[p];
+
+      // Reference coordinates for this particle (r,s[,t]) with byVDIM layout
+      IntegrationPoint ip;
+      ip.Set(rref.GetData() + dim * p, dim);
+
+      const FiniteElement& fe = *pfes->GetFE(e);
+      const int ldofs = fe.GetDof();
+
+      Vector shape(ldofs);
+      fe.CalcShape(ip, shape);  // φ_i(x_p) in this element
+
+      pfes->GetElementDofs(e, dofs);  // local dof indices
+
+      const real_t q_p = Q(p);
+
+      // Add q_p * φ_i(x_p) to b_i
+      b.AddElementVector(dofs, q_p, shape);
+   }
+}
+
+void FieldSolver::UpdatePhiGridFunction(ParticleSet& particles,
+                                        ParGridFunction& phi_gf)
+{
+   // FE space / mesh
+   ParFiniteElementSpace* pfes = phi_gf.ParFESpace();
+
+   // Particle data: Q - charges (npt x 1)
+   ParticleVector& Q = particles.Field(ParticleMover::CHARGE);
+
+   // --------------------------------------------------------
+   // 1) Make RHS and pre-subtract averaged charge density for zero-mean RHS
+   // --------------------------------------------------------
+   MPI_Comm comm = pfes->GetComm();
+   b = ComputeNeutralizingRHS(pfes, Q, comm);
+
+   // --------------------------------------------------------
+   // 2) Deposit q_p * phi_i(x_p) into a ParLinearForm (RHS b)
+   //      b_i = sum_p q_p * φ_i(x_p)
+   // --------------------------------------------------------
+   DepositCharge(pfes, Q);
+
+   // Assemble to a global true-dof RHS vector compatible with MassMatrix
+   HypreParVector B(pfes);
+   b.ParallelAssemble(B);
+
+   // ------------------------------------------------------------------
+   // 3) Solve A * phi = B with zero-mean enforcement via OrthoSolver
+   // ------------------------------------------------------------------
+   phi_gf = 0.0;
+   HypreParVector Phi_true(pfes);
+   Phi_true = 0.0;
+
+   HyprePCG solver(diffusion_matrix->GetComm());
+   solver.SetOperator(*diffusion_matrix);
+   solver.SetTol(1e-12);
+   solver.SetMaxIter(200);
+   solver.SetPrintLevel(0);
+
+   HypreBoomerAMG prec(*diffusion_matrix);
+   prec.SetPrintLevel(0);
+   solver.SetPreconditioner(prec);
+
+   OrthoSolver ortho(comm);
+   ortho.SetSolver(solver);
+   ortho.Mult(B, Phi_true);
+
+   // Map true-dof solution back to the ParGridFunction
+   phi_gf.Distribute(Phi_true);
+}
+
+void FieldSolver::UpdateEGridFunction(ParGridFunction& phi_gf,
+                                      ParGridFunction& E_gf)
+{
+   // Compute ∇φ using precomputed gradient operator
+   grad_interpolator->Mult(phi_gf, E_gf);
+   // Scale by -1 to get E = -∇φ
+   E_gf.Neg();
+}
+
+real_t FieldSolver::ComputeFieldEnergy(const ParGridFunction& E_gf) const
+{
+   // ---- Field energy: 0.5 * ∫ ||E||^2 dx ----
+   const ParFiniteElementSpace* fes = E_gf.ParFESpace();
+   const ParMesh* pmesh = fes->GetParMesh();
+
+   const int order = fes->GetMaxElementOrder();
+   const int qorder = std::max(2, 2 * order + 1);
+
+   const IntegrationRule* irs[Geometry::NumGeom];
+   for (int g = 0; g < Geometry::NumGeom; g++)
+   {
+      irs[g] = &IntRules.Get(g, qorder);
+   }
+
+   real_t field_energy = 0.0;
+
+   Vector zero(pmesh->Dimension());
+   zero = 0.0;
+   VectorConstantCoefficient zero_vec(zero);
+
+   const real_t E_l2 = E_gf.ComputeL2Error(zero_vec, irs);
+   field_energy = 0.5 * EPSILON * E_l2 * E_l2;
+
+   return field_energy;
+}
+
+void display_banner(ostream& os)
+{
+   os << R"(
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         )"
+      << endl
+      << flush;
+}

miniapps/plasma/pic/makefile

@@ -0,0 +1,85 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../../..
+MFEM_BUILD_DIR ?= ../../..
+MFEM_INSTALL_DIR ?= ../../../mfem
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/pic/,)
+CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
+   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+PAR_MINIAPPS =
+
+ifeq ($(MFEM_USE_GSLIB),YES)
+   PAR_MINIAPPS += electrostatic-pic
+endif
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS =
+else
+   MINIAPPS = $(PAR_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all lib-common clean clean-build clean-exec
+.PRECIOUS: %.o
+
+COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
+
+# If MFEM_SHARED is set, add the ../common rpath
+COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(MINIAPPS)
+
+# Rules for building the miniapps
+electrostatic-pic: electrostatic-pic.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
+
+# Rule for building lib-common
+lib-common:
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
+
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Testing: Specific execution options
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+electrostatic-pic-test-par: electrostatic-pic
+	@$(call mfem-test,$<, $(RUN_MPI), PIC miniapp,\
+		-rdi 2 -npt 40960 -k 0.2855993321 -a 0.05 -nt 200 -nx 16 -ny 16\
+		-O 1 -q 0.01181640625 -m 0.01181640625 -oci 1000 -dt 0.1)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf electrostatic-pic_* *.csv energy.csv

miniapps/plasma/plasma.hpp

@@ -0,0 +1,62 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_PLASMA_HPP
+#define MFEM_PLASMA_HPP
+
+#include <cmath>
+#include <complex>
+
+namespace mfem
+{
+
+namespace plasma
+{
+
+// Physical Constants
+
+// Permittivity of Free Space (units F/m)
+static const real_t epsilon0_ = 8.8541878176e-12;
+
+// Permeability of Free Space (units H/m)
+static const real_t mu0_ = 4.0e-7 * M_PI;
+
+// Speed of light in Free Space (units m/s)
+static const real_t c0_ = 1.0 / sqrt(epsilon0_ * mu0_);
+
+// Impedance of Free Space (units Ohm)
+static const real_t Z0_ = sqrt(mu0_ / epsilon0_);
+
+static const real_t q_     = 1.602176634e-19; // Elementary charge in coulombs
+static const real_t eV_    = 1.602176634e-19; // 1 eV in Joules
+static const real_t amu_   = 1.660539040e-27; // Atomic mass unit in kilograms
+static const real_t me_kg_ = 9.10938356e-31;  // Mass of electron in kilograms
+static const real_t me_u_  = 5.4857990907e-4; // Mass of electron in a.m.u
+
+/**
+   Returns the cyclotron frequency in radians/second
+   m is the mass in a.m.u
+   q is the charge in units of elementary electric charge
+   B is the magnetic field magnitude in tesla
+ */
+inline real_t cyclotronFrequency(real_t B, real_t m, real_t q)
+{
+   return fabs(q * q_ * B / (m * amu_));
+}
+
+typedef std::complex<real_t> complex_t;
+
+} // namespace plasma
+
+} // namespace mfem
+
+#endif // MFEM_PLASMA_HPP
+

tests/unit/CMakeLists.txt

@@ -71,10 +71,12 @@ set(UNIT_TESTS_SRCS
   linalg/test_ode2.cpp
   linalg/test_operator.cpp
   linalg/test_particlevector.cpp
+  linalg/test_petsc_nonlinear.cpp
   linalg/test_sparsesmoothers.cpp
   linalg/test_vector.cpp
   mesh/mesh_test_utils.cpp
   mesh/test_exodus_reader.cpp
+  mesh/test_mfem_mesh_reader.cpp
   mesh/test_exodus_writer.cpp
   mesh/test_face_orientations.cpp
   mesh/test_fms.cpp
@@ -120,6 +122,7 @@ set(UNIT_TESTS_SRCS
   fem/test_fe_pos.cpp
   fem/test_fe_symmetry.cpp
   fem/test_fe.cpp
+  fem/test_fespace_get_ess_true_dofs.cpp
   fem/test_get_value.cpp
   fem/test_getderivative.cpp
   fem/test_getgradient.cpp

tests/unit/catch.hpp

@@ -295,8 +295,13 @@ namespace Catch {
 // Otherwise all supported compilers support COUNTER macro,
 // but user still might want to turn it off
 #if ( !defined(__JETBRAINS_IDE__) || __JETBRAINS_IDE__ >= 20170300L )
+#if ( !(defined(__clang__) && __clang_major__ >= 22 ) )
+// don't use __COUNTER__ if compiling with clang 22+ to avoid compiler warning
+// https://github.com/llvm/llvm-project/pull/162662
+// TODO: can enable if building with C2y
     #define CATCH_INTERNAL_CONFIG_COUNTER
 #endif
+#endif
 
 ////////////////////////////////////////////////////////////////////////////////
 

tests/unit/data/compass-testing.mesh

@@ -0,0 +1,118 @@
+MFEM mesh v1.3
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+12
+10 2 7 0 1
+11 2 0 7 2
+12 2 9 0 2
+13 2 0 9 3
+14 2 11 0 3
+15 2 0 11 4
+16 2 5 0 4
+17 2 0 5 1
+9 3 1 5 6 7
+9 3 2 7 8 9
+9 3 3 9 10 11
+9 3 4 11 12 5
+
+attribute_sets
+16
+"Base" 1 9
+"E Even" 1 16
+"E Odd" 1 17
+"East"
+2
+16
+17
+"N Even" 1 10
+"N Odd" 1 11
+"North" 2 10     11
+"Rose" 8 10 11     12
+13    14 
+15    16     17
+"Rose Even" 4
+10
+12
+14
+16
+"Rose Odd"
+4
+11
+13
+15
+17
+"S Even" 1 14
+"S Odd" 1 15
+South 2
+14
+15
+"W Even" 1 12
+"W Odd" 1 13
+West 2 12 13
+
+boundary
+8
+1 1 5 6
+2 1 6 7
+3 1 7 8
+4 1 8 9
+5 1 9 10
+6 1 10 11
+7 1 11 12
+8 1 12 5
+
+bdr_attribute_sets
+13
+"Boundary" 8 1 2 3 4 5 6 7 8
+"ENE" 1 1
+"ESE" 1 8
+"Eastern Boundary" 2 1 8
+"NNE" 1 2
+"NNW" 1 3
+"Northern Boundary"
+2
+2
+3
+"SSE" 1 7
+"SSW" 1 6
+"Southern Boundary" 2
+6
+7
+"WNW" 1 4
+"WSW" 1 5
+"Western Boundary" 2 4
+5
+
+vertices
+13
+2
+0 0
+0.14142136 0.14142136
+-0.14142136 0.14142136
+-0.14142136 -0.14142136
+0.14142136 -0.14142136
+1 0
+0.70710678 0.70710678
+0 1
+-0.70710678 0.70710678
+-1 0
+-0.70710678 -0.70710678
+0 -1
+0.70710678 -0.70710678
+
+mfem_mesh_end

tests/unit/fem/test_face_restriction.cpp

@@ -117,3 +117,52 @@ TEST_CASE("Vector FE Face Restriction", "[FaceRestriction]")
    gf2 -= gf;
    REQUIRE(gf2.Normlinf() == MFEM_Approx(0.0));
 }
+
+#ifdef MFEM_USE_MPI
+
+TEST_CASE("L2 Face Restriction", "[FaceRestriction][Parallel]")
+{
+   const int dim = GENERATE(2, 3);
+   constexpr int nx = 3;
+   constexpr int order = 2;
+   constexpr int vdim = 2;
+   const Ordering::Type ordering = GENERATE(Ordering::byNODES, Ordering::byVDIM);
+
+   Mesh serial_mesh = MakeCartesianMesh(nx, dim);
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+
+   L2_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace fes(&mesh, &fec, vdim, ordering);
+
+   auto *R = fes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,
+                                    FaceType::Interior);
+
+   Vector vals({1.0, 2.0});
+   VectorConstantCoefficient coeff(vals);
+
+   ParGridFunction gf(&fes);
+   gf.ProjectCoefficient(coeff);
+
+   Vector face_vec(R->Height());
+   R->Mult(gf, face_vec);
+
+   const int nf = mesh.GetNFbyType(FaceType::Interior);
+   const int face_dofs = fes.GetTypicalTraceElement()->GetDof();
+   auto h_face_vec = Reshape(face_vec.HostRead(), face_dofs, vdim, 2, nf);
+
+   for (int f = 0; f < nf; ++f)
+   {
+      for (int m = 0; m < 2; ++m)
+      {
+         for (int c = 0; c < vdim; ++c)
+         {
+            for (int i = 0; i < face_dofs; ++i)
+            {
+               REQUIRE(h_face_vec(i, c, m, f) == vals[c]);
+            }
+         }
+      }
+   }
+}
+
+#endif

tests/unit/fem/test_fe.cpp

@@ -281,8 +281,10 @@ TEST_CASE("Nedelec Segment Finite Element",
             REQUIRE( fe.GetRangeType()      == (int) FiniteElement::VECTOR  );
             REQUIRE( fe.GetMapType()        == (int) FiniteElement::H_CURL  );
             REQUIRE( fe.GetDerivType()      == (int) FiniteElement::NONE    );
-            REQUIRE( fe.GetDerivRangeType() == (int) FiniteElement::SCALAR  );
-            REQUIRE( fe.GetDerivMapType()   == (int) FiniteElement::INTEGRAL);
+            REQUIRE( fe.GetDerivRangeType() ==
+                     (int) FiniteElement::UNKNOWN_RANGE_TYPE);
+            REQUIRE( fe.GetDerivMapType()   ==
+                     (int) FiniteElement::UNKNOWN_MAP_TYPE);
          }
       }
       SECTION("Sizes for p = " + std::to_string(p))

tests/unit/fem/test_fespace_get_ess_true_dofs.cpp

@@ -0,0 +1,74 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+TEST_CASE("FESpace Get Essential True DOFs",
+          "[FESpace Get Essential True DOFs]")
+{
+   std::cout << "Testing get essential true dofs" << std::endl;
+   int order_h1 = 3, n = 2, dim = 3;
+
+   Mesh mesh = Mesh::MakeCartesian3D(
+                  n, n, n, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
+   mesh.SetCurvature(order_h1);
+
+   H1_FECollection fec(order_h1, dim);
+   FiniteElementSpace fe_space(&mesh, &fec, dim);
+
+
+   const int num_bdr_attr = fe_space.GetMesh()->bdr_attributes.Max();
+
+
+   Array<int> ess_tdofs_2d, ess_tdofs_1d, ess_tdofs_tmp, ess_bdrs;
+   Array2D<bool> comps(num_bdr_attr, dim);
+   ess_bdrs.SetSize(num_bdr_attr);
+   comps = false;
+   ess_bdrs = 0;
+
+   // simple xy boundary condition on all surfaces
+   // could do something more complex but don't really want to...
+   for (int i = 0; i < num_bdr_attr; i++)
+   {
+      ess_bdrs[i] = 1;
+      comps(i, 0) = true;
+      comps(i, 2) = true;
+   }
+
+   fe_space.GetEssentialTrueDofs(ess_bdrs, ess_tdofs_2d, comps);
+
+   // Now for the old way
+   fe_space.GetEssentialTrueDofs(ess_bdrs, ess_tdofs_tmp, 0);
+   ess_tdofs_1d.Append(ess_tdofs_tmp);
+   ess_tdofs_tmp.DeleteAll();
+   fe_space.GetEssentialTrueDofs(ess_bdrs, ess_tdofs_tmp, 2);
+   ess_tdofs_1d.Append(ess_tdofs_tmp);
+   // Sort the 2 arrays in order to compare them
+   ess_tdofs_2d.Sort();
+   ess_tdofs_1d.Sort();
+
+   int diff = 0;
+
+   for (int i = 0; i < ess_tdofs_2d.Size(); i++)
+   {
+      diff += std::abs(ess_tdofs_2d[i] - ess_tdofs_1d[i]);
+   }
+
+   std::cout << "Difference in essential tdofs approaches is: " << diff <<
+             std::endl;
+
+   REQUIRE(diff == 0);
+
+}
+

tests/unit/fem/test_intrules.cpp

@@ -105,39 +105,39 @@ TEST_CASE("Integration rule order initialization", "[IntegrationRules]")
    SECTION("Segment rule constructed by accessing square rule")
    {
       auto &quad5_ir = intrules.Get(Geometry::SQUARE, 5);
-      REQUIRE(quad5_ir.GetOrder() == 5);
+      REQUIRE(quad5_ir.GetOrder() >= 5);
       // The segment integration rule of order 5 is lazy constructed when we get
       // the square integration rule of order 5. Make sure its order was
       // properly set:
       auto &line5_ir = intrules.Get(Geometry::SEGMENT, 5);
-      REQUIRE(line5_ir.GetOrder() == 5);
+      REQUIRE(line5_ir.GetOrder() >= 5);
    }
 
    SECTION("Segment rule constructed by accessing cube rule")
    {
       auto &hex7_ir = intrules.Get(Geometry::CUBE, 7);
-      REQUIRE(hex7_ir.GetOrder() == 7);
+      REQUIRE(hex7_ir.GetOrder() >= 7);
       // The segment integration rule of order 7 is lazy constructed when we get
       // the cube integration rule of order 7. Make sure its order was properly
       // set:
       auto &line7_ir = intrules.Get(Geometry::SEGMENT, 7);
-      REQUIRE(line7_ir.GetOrder() == 7);
+      REQUIRE(line7_ir.GetOrder() >= 7);
    }
 
    SECTION("Segment and triangle rules constructed by accessing prism rule")
    {
       auto &prism3_ir = intrules.Get(Geometry::PRISM, 3);
-      REQUIRE(prism3_ir.GetOrder() == 3);
+      REQUIRE(prism3_ir.GetOrder() >= 3);
       // The segment integration rule of order 3 is lazy constructed when we get
       // the prism integration rule of order 3. Make sure its order was properly
       // set:
       auto &line3_ir = intrules.Get(Geometry::SEGMENT, 3);
-      REQUIRE(line3_ir.GetOrder() == 3);
+      REQUIRE(line3_ir.GetOrder() >= 3);
       // The triangle integration rule of order 3 is lazy constructed when we
       // get the prism integration rule of order 3. Make sure its order was
       // properly set:
       auto &tri3_ir = intrules.Get(Geometry::TRIANGLE, 3);
-      REQUIRE(tri3_ir.GetOrder() == 3);
+      REQUIRE(tri3_ir.GetOrder() >= 3);
    }
 }
 
@@ -271,3 +271,43 @@ TEST_CASE("Simplex integration rules", "[SimplexRules]")
       }
    }
 }
+
+
+// Monomial exactness is tested by [SimplexRules] above, which now uses
+// positive-weight rules by default. The tests below verify properties
+// specific to the positive-weight rules: weight positivity, stability,
+// and interior point placement.
+
+TEST_CASE("Simplex rule positivity", "[IntegrationRules]")
+{
+   IntegrationRules rules;
+
+   SECTION("triangle rules have all positive weights for orders 0-25")
+   {
+      for (int order = 0; order <= 25; order++)
+      {
+         const IntegrationRule &ir = rules.Get(Geometry::TRIANGLE, order);
+         for (int i = 0; i < ir.GetNPoints(); i++)
+         {
+            INFO("order=" << order << ", point=" << i);
+            REQUIRE(ir.IntPoint(i).weight > 0.0);
+         }
+      }
+   }
+
+   SECTION("tet rules have all positive weights for orders 0-20")
+   {
+      for (int order = 0; order <= 20; order++)
+      {
+         const IntegrationRule &ir =
+            rules.Get(Geometry::TETRAHEDRON, order);
+         for (int i = 0; i < ir.GetNPoints(); i++)
+         {
+            INFO("order=" << order << ", point=" << i);
+            REQUIRE(ir.IntPoint(i).weight > 0.0);
+         }
+      }
+   }
+}
+
+

tests/unit/fem/test_project_bdr.cpp

@@ -25,15 +25,201 @@ void Func_3D_lin(const Vector &x, Vector &v)
    v[2] = -2.572 * x[0] + 1.321 * x[1] + 3.234 * x[2];
 }
 
+TEST_CASE("3D ProjectBdrCoefficientNormal Vector",
+          "[GridFunction]"
+          "[VectorGridFunctionCoefficient]")
+{
+   const int n = 1;
+   const int dim = 3;
+   const int order = 1;
+
+   const double tol = 1e-6;
+
+   for (int type = (int)Element::TETRAHEDRON;
+        type <= (int)Element::HEXAHEDRON; type++)
+   {
+      Mesh mesh = Mesh::MakeCartesian3D(
+                     n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
+
+      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
+
+      SECTION("3D GetVectorValue tests for element type " +
+              std::to_string(type))
+      {
+         RT_FECollection  rt_fec(order+1, dim);
+
+         FiniteElementSpace  rt_fespace(&mesh,  &rt_fec);
+
+         GridFunction  rt_x( &rt_fespace);
+
+         VectorGridFunctionCoefficient  rt_xCoef( &rt_x);
+
+         Array<int> bdr_marker(6);
+
+         Vector normal(dim);
+         Vector  f_val(dim);
+         Vector rt_val(dim);
+
+         for (int b = 1; b<=6; b++)
+         {
+            bdr_marker = 0;
+            bdr_marker[b-1] = 1;
+
+            rt_x = 0.0;
+            rt_x.ProjectBdrCoefficientNormal(funcCoef, bdr_marker);
+
+            for (int be = 0; be < mesh.GetNBE(); be++)
+            {
+               Element *e = mesh.GetBdrElement(be);
+               if (e->GetAttribute() != b) { continue; }
+
+               ElementTransformation *T = mesh.GetBdrElementTransformation(be);
+               const FiniteElement   *fe = rt_fespace.GetBE(be);
+               const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
+                                                        2*order + 2);
+
+               double rt_err = 0.0;
+
+               for (int j=0; j<ir.GetNPoints(); j++)
+               {
+                  const IntegrationPoint &ip = ir.IntPoint(j);
+                  T->SetIntPoint(&ip);
+
+                  CalcOrtho(T->Jacobian(), normal);
+
+                  funcCoef.Eval(f_val, *T, ip);
+                  rt_xCoef.Eval(rt_val, *T, ip);
+
+                  rt_val -= f_val;
+
+                  double rt_dist = rt_val * normal;
+
+                  rt_err += rt_dist;
+
+                  if (verbose_tests && rt_dist > tol)
+                  {
+                     mfem::out << be << ":" << j << " rt ("
+                               << f_val[0] << "," << f_val[1] << "," << f_val[2]
+                               << ") vs. ("
+                               << rt_val[0] << "," << rt_val[1] << ","
+                               << rt_val[2] << ") " << rt_dist << std::endl;
+                  }
+               }
+               rt_err  /= ir.GetNPoints();
+
+               REQUIRE( rt_err == MFEM_Approx(0.0));
+            }
+         }
+      }
+   }
+}
+
+TEST_CASE("3D ProjectBdrCoefficientNormal Scalar",
+          "[GridFunction]"
+          "[VectorGridFunctionCoefficient]")
+{
+   const int n = 1;
+   const int dim = 3;
+   const int order = 1;
+
+   const double tol = 1e-6;
+
+   const char bdrs_axis[] = {2, 1, 0, 1, 0, 2};
+   const char bdrs_sign[] = {-1, -1, +1, +1, -1, +1};
+
+   for (int type = (int)Element::TETRAHEDRON;
+        type <= (int)Element::HEXAHEDRON; type++)
+   {
+      Mesh mesh = Mesh::MakeCartesian3D(
+                     n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
+
+      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
+
+      SECTION("3D GetVectorValue tests for element type " +
+              std::to_string(type))
+      {
+         RT_FECollection  rt_fec(order+1, dim);
+
+         FiniteElementSpace  rt_fespace(&mesh,  &rt_fec);
+
+         GridFunction  rt_x( &rt_fespace);
+
+         VectorGridFunctionCoefficient  rt_xCoef( &rt_x);
+
+         Array<int> bdr_marker(6);
+
+         Vector normal(dim);
+         Vector  f_val(dim);
+         Vector rt_val(dim);
+
+         for (int b = 1; b<=6; b++)
+         {
+            bdr_marker = 0;
+            bdr_marker[b-1] = 1;
+
+            rt_x = 0.0;
+
+            normal = 0.;
+            normal(bdrs_axis[b-1]) = (bdrs_sign[b-1] > 0)?(+1.):(-1.);
+            VectorConstantCoefficient normCoef(normal);
+            InnerProductCoefficient prodCoef(funcCoef, normCoef);
+            rt_x.ProjectBdrCoefficientNormal(prodCoef, bdr_marker);
+
+            for (int be = 0; be < mesh.GetNBE(); be++)
+            {
+               Element *e = mesh.GetBdrElement(be);
+               if (e->GetAttribute() != b) { continue; }
+
+               ElementTransformation *T = mesh.GetBdrElementTransformation(be);
+               const FiniteElement   *fe = rt_fespace.GetBE(be);
+               const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
+                                                        2*order + 2);
+
+               double rt_err = 0.0;
+
+               for (int j=0; j<ir.GetNPoints(); j++)
+               {
+                  const IntegrationPoint &ip = ir.IntPoint(j);
+                  T->SetIntPoint(&ip);
+
+                  CalcOrtho(T->Jacobian(), normal);
+
+                  funcCoef.Eval(f_val, *T, ip);
+                  rt_xCoef.Eval(rt_val, *T, ip);
+
+                  rt_val -= f_val;
+
+                  double rt_dist = rt_val * normal;
+
+                  rt_err += rt_dist;
+
+                  if (verbose_tests && rt_dist > tol)
+                  {
+                     mfem::out << be << ":" << j << " rt ("
+                               << f_val[0] << "," << f_val[1] << "," << f_val[2]
+                               << ") vs. ("
+                               << rt_val[0] << "," << rt_val[1] << ","
+                               << rt_val[2] << ") " << rt_dist << std::endl;
+                  }
+               }
+               rt_err  /= ir.GetNPoints();
+
+               REQUIRE( rt_err == MFEM_Approx(0.0));
+            }
+         }
+      }
+   }
+}
+
 TEST_CASE("3D ProjectBdrCoefficientTangent",
           "[GridFunction]"
           "[VectorGridFunctionCoefficient]")
 {
-   int n = 1;
-   int dim = 3;
-   int order = 1;
+   const int n = 1;
+   const int dim = 3;
+   const int order = 1;
 
-   double tol = 1e-6;
+   const double tol = 1e-6;
 
    for (int type = (int)Element::TETRAHEDRON;
         type <= (int)Element::HEXAHEDRON; type++)

tests/unit/fem/test_var_order.cpp

@@ -271,6 +271,8 @@ TEST_CASE("Variable Order FiniteElementSpace",
 
       const auto space_type = GENERATE(SpaceType::RT, SpaceType::ND);
       const int dim = GENERATE(2, 3);
+      CAPTURE(space_type);
+      CAPTURE(dim);
 
       Mesh mesh = MakeCartesianMesh(dim == 2 ? 4 : 2, dim);
       mesh.EnsureNCMesh();
@@ -698,7 +700,14 @@ static void TestSolveVec(FiniteElementSpace &fespace)
 
    GridFunction x(&fespace);
    x = 0.0;
-   x.ProjectBdrCoefficient(exsol, ess_attr);
+   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
+   {
+      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
+   }
+   else
+   {
+      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
+   }
 
    // Assemble the linear form
    LinearForm lf(&fespace);
@@ -1082,7 +1091,14 @@ static void TestSolveParVec(ParFiniteElementSpace &fespace)
 
    ParGridFunction x(&fespace);
    x = 0.0;
-   x.ProjectBdrCoefficient(exsol, ess_attr);
+   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
+   {
+      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
+   }
+   else
+   {
+      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
+   }
 
    // Assemble the linear form
    ParLinearForm lf(&fespace);

tests/unit/general/test_arrays_by_name.cpp

@@ -200,3 +200,39 @@ TEST_CASE("ArraysByName Sort/Unique Methods", "[ArraysByName]")
       }
    }
 }
+
+TEST_CASE("ArraysByName Print/Load Methods", "[ArraysByName]")
+{
+   ArraysByName<int> abn;
+
+   FillArraysByName(abn);
+
+   // Print object to string using default format
+   std::ostringstream oss1;
+   abn.Print(oss1);
+
+   // Load new object from printed output
+   ArraysByName<int> abn_load1;
+   std::istringstream iss1(oss1.str());
+   abn_load1.Load(iss1);
+   REQUIRE(abn == abn_load1);
+
+   // Print object to string using one line per array
+   std::ostringstream oss2;
+   oss2 << abn.Size() << '\n';
+   for (auto a : abn)
+   {
+      oss2 << '"' << a.first << "\" " << a.second.Size();
+      for (auto d : a.second)
+      {
+         oss2 << ' ' << d;
+      }
+      oss2 << '\n';
+   }
+
+   // Load new object from printed output
+   ArraysByName<int> abn_load2;
+   std::istringstream iss2(oss2.str());
+   abn_load2.Load(iss2);
+   REQUIRE(abn == abn_load2);
+}

tests/unit/general/test_text.cpp

@@ -30,3 +30,29 @@ TEST_CASE("String Manipulation", "[General]")
       }
    }
 }
+
+TEST_CASE("Quoted String Input", "[General]")
+{
+   const auto test_strings =
+   {
+      "Test",
+      "Test with spaces",
+      "Test with \"quoted text\"",
+      "Test string ending with \\",
+      "\nTest with\tvarious white\v\rspace characters.",
+      "Test with some unicode characters: ∆, ∉, ∑, 🍎."
+   };
+
+   for (const auto c_str : test_strings)
+   {
+      CAPTURE(c_str);
+      const std::string str(c_str);
+      std::stringstream ss;
+      ss << std::quoted(str);
+
+      std::string read_str;
+      int error = parse_quoted_string(read_str, ss);
+      CHECK(error == 0);
+      CHECK(read_str == str);
+   }
+}

tests/unit/linalg/test_petsc_nonlinear.cpp

@@ -0,0 +1,74 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+#if defined(MFEM_USE_MPI) && defined(MFEM_USE_PETSC)
+
+namespace
+{
+struct PetscSession
+{
+   PetscSession() { MFEMInitializePetsc(); }
+   ~PetscSession() { MFEMFinalizePetsc(); }
+};
+
+class IdentityGradientOperator : public IdentityOperator
+{
+public:
+   IdentityGradientOperator() : IdentityOperator(1), _jac(1)
+   {
+      _jac.Add(0, 0, 1.0);
+      _jac.Finalize();
+   }
+
+   Operator &GetGradient(const Vector &) const override
+   {
+      return const_cast<SparseMatrix &>(_jac);
+   }
+
+private:
+   SparseMatrix _jac;
+};
+}
+
+TEST_CASE("PetscNonlinearSolver accepts non-empty rhs", "[Parallel][PETSc]")
+{
+   static PetscSession petsc_session;
+
+   IdentityGradientOperator oper;
+   PetscNonlinearSolver solver(MPI_COMM_WORLD, "nl_");
+   solver.SetRelTol(1.0e-12);
+   solver.SetAbsTol(1.0e-12);
+   solver.SetMaxIter(5);
+   solver.SetPrintLevel(0);
+   solver.SetJacobianType(Operator::PETSC_MATAIJ);
+   solver.SetOperator(oper);
+
+   Vector x(1);
+
+   Vector empty_rhs;
+   x = 0.0;
+   solver.Mult(empty_rhs, x);
+   REQUIRE(x(0) == MFEM_Approx(0.0));
+
+   Vector nonempty_rhs(1);
+   nonempty_rhs(0) = 2.5;
+   x = 0.0;
+   solver.Mult(nonempty_rhs, x);
+   REQUIRE(x.Size() == 1);
+   REQUIRE(x(0) == MFEM_Approx(nonempty_rhs(0)));
+}
+
+#endif

tests/unit/mesh/test_mfem_mesh_reader.cpp

@@ -0,0 +1,108 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include <algorithm>
+#include <string>
+#include <utility>
+#include <vector>
+
+using namespace mfem;
+
+TEST_CASE("MFEM Mesh Named Attributes", "[Mesh]")
+{
+   // Path relative to the directory tests/unit
+   Mesh mesh("data/compass-testing.mesh");
+
+   REQUIRE(mesh.Dimension() == 2);
+   REQUIRE(mesh.GetNE() == 12);
+   REQUIRE(mesh.GetNV() == 13);
+
+   REQUIRE(mesh.attribute_sets.attr_sets.Size() == 16);
+   REQUIRE(mesh.bdr_attribute_sets.attr_sets.Size() == 13);
+
+   std::vector<std::pair<std::string, std::vector<int>>> expected_attr_sets =
+   {
+      {"Base",  {9}},
+      {"E Even",  {16}},
+      {"E Odd",  {17}},
+      {"East",  {16, 17}},
+      {"N Even",  {10}},
+      {"N Odd",  {11}},
+      {"North",  {10, 11}},
+      {"Rose",  {10, 11, 12, 13, 14, 15, 16, 17}},
+      {"Rose Even",  {10, 12, 14, 16}},
+      {"Rose Odd",  {11, 13, 15, 17}},
+      {"S Even",  {14}},
+      {"S Odd",  {15}},
+      {"South",  {14, 15}},
+      {"W Even",  {12}},
+      {"W Odd",  {13}},
+      {"West",  {12, 13}}
+   };
+
+   for (auto const &attr_name_index_pair: expected_attr_sets )
+   {
+      REQUIRE(mesh.attribute_sets.AttributeSetExists(
+                 attr_name_index_pair.first));
+
+      auto const &attr_set = mesh.attribute_sets.GetAttributeSet(
+                                attr_name_index_pair.first);
+      auto const &expected_attr_set = attr_name_index_pair.second;
+
+      REQUIRE(static_cast<std::size_t>(attr_set.Size()) ==
+              expected_attr_set.size());
+
+      bool const elements_equal = std::equal(attr_set.begin(), attr_set.end(),
+                                             expected_attr_set.begin());
+
+      REQUIRE(elements_equal);
+   }
+
+   std::vector<std::pair<std::string, std::vector<int>>> expected_bdr_attr_sets
+   =
+   {
+      {"Boundary",  {1, 2, 3, 4, 5, 6, 7, 8}},
+      {"ENE",  { 1}},
+      {"ESE",  { 8}},
+      {"Eastern Boundary",  {1, 8}},
+      {"NNE",  { 2}},
+      {"NNW",  { 3}},
+      {"Northern Boundary", {2, 3}},
+      {"SSE",  { 7}},
+      {"SSW",  { 6}},
+      {"Southern Boundary", {6,7}},
+      {"WNW",  { 4}},
+      {"WSW",  { 5}},
+      {"Western Boundary", {4,5}}
+   };
+
+   for (auto const &attr_bdr_name_index_pair: expected_bdr_attr_sets )
+   {
+      REQUIRE(mesh.bdr_attribute_sets.AttributeSetExists(
+                 attr_bdr_name_index_pair.first));
+
+      auto const &bdr_attr_set = mesh.bdr_attribute_sets.GetAttributeSet(
+                                    attr_bdr_name_index_pair.first);
+      auto const &expected_bdr_attr_set = attr_bdr_name_index_pair.second;
+
+      REQUIRE(static_cast<std::size_t>(bdr_attr_set.Size()) ==
+              expected_bdr_attr_set.size());
+
+      bool const elements_equal = std::equal(bdr_attr_set.begin(),
+                                             bdr_attr_set.end(),
+                                             expected_bdr_attr_set.begin());
+
+      REQUIRE(elements_equal);
+   }
+}

tests/unit/mesh/test_psubmesh.cpp

@@ -486,8 +486,14 @@ void multidomain_test_3d(FECType fec_type)
          {
             cylinder_gf.ProjectCoefficient(vcoeff);
             outer_gf.ProjectCoefficient(vcoeff);
-            outer_gf.ProjectBdrCoefficient(vzerocoeff,
-                                           outer_cyl_surf_marker);
+            if (fec_type == FECType::RT)
+            {
+               outer_gf.ProjectBdrCoefficientNormal(vzerocoeff, outer_cyl_surf_marker);
+            }
+            else
+            {
+               outer_gf.ProjectBdrCoefficientTangent(vzerocoeff, outer_cyl_surf_marker);
+            }
             outer_gf_ex.ProjectCoefficient(vcoeff);
          }
          ParSubMesh::Transfer(cylinder_gf, outer_gf);
@@ -507,8 +513,14 @@ void multidomain_test_3d(FECType fec_type)
          {
             outer_gf.ProjectCoefficient(vcoeff);
             cylinder_gf.ProjectCoefficient(vcoeff);
-            cylinder_gf.ProjectBdrCoefficient(vzerocoeff,
-                                              cylinder_cyl_surf_marker);
+            if (fec_type == FECType::RT)
+            {
+               cylinder_gf.ProjectBdrCoefficientNormal(vzerocoeff, cylinder_cyl_surf_marker);
+            }
+            else
+            {
+               cylinder_gf.ProjectBdrCoefficientTangent(vzerocoeff, cylinder_cyl_surf_marker);
+            }
             cylinder_gf_ex.ProjectCoefficient(vcoeff);
          }
          ParSubMesh::Transfer(outer_gf, cylinder_gf);