Fix signed char issue in Device::GetUUID

Fixed missing parenthesis in the comment
Adjust seed values in sample runs in ex12p to ensure LOBPCG convergence in
2026-04-17 11:38:22 -04:00 · 2026-04-17 11:38:22 -04:00 · 2026-04-17 11:38:22 -04:00 · 2026-04-17 11:38:22 -04:00 · 2026-04-17 11:38:22 -04:00 · 2026-04-17 11:38:22 -04:00
62 changed files with 3242 additions and 649 deletions
@@ -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:
@@ -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:
@@ -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}}
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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

@@ -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
@@ -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}}
@@ -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
@@ -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
@@ -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:
@@ -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)
 ===========================

@@ -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()
@@ -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
@@ -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
@@ -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)
   {
@@ -97,13 +97,7 @@ int main(int argc, char *argv[])
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   args.PrintOptions(cout);
+   args.ParseCheck();

   // 2. Read the mesh from the mesh file.
   const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -128,7 +122,7 @@ int main(int argc, char *argv[])
   *nodes /= scale;

   // 4. Define the necessary finite element spaces on the mesh.
-   H1_FECollection H1fec(order+1, dim);
+   H1Bubble_FECollection H1fec(order, order - 1, dim);
   FiniteElementSpace H1fes(&mesh, &H1fec);

   L2_FECollection L2fec(order-1, dim);
@@ -103,19 +103,7 @@ int main(int argc, char *argv[])
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
+   args.ParseCheck();

   // 2. Read the mesh from the mesh file.
   const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -143,7 +131,7 @@ int main(int argc, char *argv[])
   mesh.Clear();

   // 4. Define the necessary finite element spaces on the mesh.
-   H1_FECollection H1fec(order+1, dim);
+   H1Bubble_FECollection H1fec(order, order - 1, dim);
   ParFiniteElementSpace H1fes(&pmesh, &H1fec);

   L2_FECollection L2fec(order-1, dim);
@@ -73,6 +73,7 @@ set(SRCS
  fe/fe_base.cpp
  fe/fe_fixed_order.cpp
  fe/fe_h1.cpp
+  fe/fe_h1_bubble.cpp
  fe/fe_l2.cpp
  fe/fe_nd.cpp
  fe/fe_nurbs.cpp
@@ -133,7 +134,7 @@ set(SRCS
  tmop/assemble/diag2.cpp
  tmop/assemble/grad2_limit.cpp
  tmop/assemble/grad2.cpp
-  tmop/assemble/diag3_limit.cpp  
+  tmop/assemble/diag3_limit.cpp
  tmop/assemble/diag3.cpp
  tmop/assemble/grad3_limit.cpp
  tmop/assemble/grad3.cpp
@@ -221,6 +222,7 @@ set(HDRS
  fe/fe_base.hpp
  fe/fe_fixed_order.hpp
  fe/fe_h1.hpp
+  fe/fe_h1_bubble.hpp
  fe/fe_l2.hpp
  fe/fe_nd.hpp
  fe/fe_nurbs.hpp
@@ -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);
            }
         }
@@ -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:
@@ -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; }
@@ -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");
   }
 }

@@ -20,6 +20,7 @@
 #include "fe/fe_base.hpp"
 #include "fe/fe_fixed_order.hpp"
 #include "fe/fe_h1.hpp"
+#include "fe/fe_h1_bubble.hpp"
 #include "fe/fe_nd.hpp"
 #include "fe/fe_rt.hpp"
 #include "fe/fe_l2.hpp"
@@ -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);
   }
@@ -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). */
+   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).
+   */
+   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 { return space_dim; }
 };

 /// @brief Class for computing 1D special polynomials and their associated basis
@@ -0,0 +1,973 @@
+// 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.
+
+// H1 Finite Element classes
+
+#include "fe_h1_bubble.hpp"
+
+namespace mfem
+{
+
+using namespace std;
+
+H1Bubble_TriangleElement::H1Bubble_TriangleElement(int p, int q, int btype)
+   : NodalFiniteElement(2, Geometry::TRIANGLE, 3*p + ((q+1)*(q+2))/2,
+                        max(p, 3 + q), FunctionSpace::Pk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 3, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   shape_l.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+   dshape_l.SetSize(n1d);
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set2(cp[0], cp[0]);
+   Nodes.IntPoint(1).Set2(cp[p], cp[0]);
+   Nodes.IntPoint(2).Set2(cp[0], cp[p]);
+
+   // edges
+   int o = 3;
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
+   }
+
+   // Interior P_{q+3} nodes
+   for (int j = 1; j < q + 3; j++)
+   {
+      for (int i = 1; i + j < q + 3; i++)
+      {
+         const real_t w = cp2[i] + cp2[j] + cp2[q+3-i-j];
+         Nodes.IntPoint(o++).Set2(cp2[i]/w, cp2[j]/w);
+      }
+   }
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(n1d), shape_y(n1d), shape_l(n1d);
+#endif
+
+   DenseMatrix Tt(dof, npq);
+   for (int k = 0; k < dof; ++k)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(k);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
+      o = 0;
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i + j <= p; i++)
+         {
+            Tt(k, o++) = shape_x[i]*shape_y[j]*shape_l[p-i-j];
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      poly1d.CalcBasis(q, 1. - ip.x - ip.y, shape_l);
+      const real_t b_T = ip.x * ip.y * (1 - ip.x - ip.y);
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i + j <= q; i++)
+         {
+            Tt(k, o++) = b_T*shape_x[i]*shape_y[j]*shape_l[q-i-j];
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_TriangleElement::CalcShape(const IntegrationPoint &ip,
+                                         Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
+   Vector shape_x(n1d), shape_y(n1d), shape_l(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i + j <= p; i++)
+      {
+         u(o++) = shape_x[i]*shape_y[j]*shape_l[p-i-j];
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   poly1d.CalcBasis(q, 1. - ip.x - ip.y, shape_l);
+   const real_t b_T = ip.x * ip.y * (1 - ip.x - ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i + j <= q; i++)
+      {
+         u(o++) = b_T*shape_x[i]*shape_y[j]*shape_l[q-i-j];
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_TriangleElement::CalcDShape(const IntegrationPoint &ip,
+                                          DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
+   Vector  shape_x(n1d),  shape_y(n1d),  shape_l(n1d);
+   Vector dshape_x(n1d), dshape_y(n1d), dshape_l(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   const real_t lambda = 1.0 - ip.x - ip.y;
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, lambda, shape_l, dshape_l);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i + j <= p; i++)
+      {
+         int k = p - i - j;
+         du(o,0) = (dshape_x[i]*shape_l[k] - shape_x[i]*dshape_l[k])*shape_y[j];
+         du(o,1) = (dshape_y[j]* shape_l[k] - shape_y[j]*dshape_l[k])*shape_x[i];
+         o++;
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(q, lambda, shape_l, dshape_l);
+   const real_t b_T = ip.x * ip.y * lambda;
+   const real_t dxb_T = ip.y * (lambda - ip.x);
+   const real_t dyb_T = ip.x * (lambda - ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i + j <= q; i++)
+      {
+         int k = q - i - j;
+         du(o,0) = shape_y[j]*(dxb_T*shape_x[i]*shape_l[k]
+                               + b_T*dshape_x[i]*shape_l[k]
+                               - b_T*shape_x[i]*dshape_l[k]);
+         du(o,1) = shape_x[i]*(dyb_T*shape_y[j]*shape_l[k]
+                               + b_T*dshape_y[j]*shape_l[k]
+                               - b_T*shape_y[j]*dshape_l[k]);
+         o++;
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+H1Bubble_QuadrilateralElement::H1Bubble_QuadrilateralElement(
+   int p, int q, int btype)
+   : NodalFiniteElement(2, Geometry::SQUARE, 4*p + (q+1)*(q+1),
+                        max(p, 2 + q), FunctionSpace::Qk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 2, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set2(cp[0], cp[0]);
+   Nodes.IntPoint(1).Set2(cp[p], cp[0]);
+   Nodes.IntPoint(2).Set2(cp[p], cp[p]);
+   Nodes.IntPoint(3).Set2(cp[0], cp[p]);
+
+   // edges
+   int o = 4;
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[p], cp[i]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[p-i], cp[p]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
+   }
+
+   // interior P_{q+2} nodes
+   for (int j = 1; j < q+2; j++)
+   {
+      for (int i = 1; i < q+2; i++)
+      {
+         Nodes.IntPoint(o++).Set2(cp2[i], cp2[j]);
+      }
+   }
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(n1d), shape_y(n1d);
+#endif
+
+   DenseMatrix Tt(dof, npq);
+   for (int k = 0; k < dof; ++k)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(k);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+
+      o = 0;
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i <= p; i++)
+         {
+            Tt(k, o++) = shape_x[i]*shape_y[j];
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i <= q; i++)
+         {
+            Tt(k, o++) = b_T*shape_x[i]*shape_y[j];
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
+                                              Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i <= p; i++)
+      {
+         u(o++) = shape_x[i]*shape_y[j];
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i <= q; i++)
+      {
+         u(o++) = b_T*shape_x[i]*shape_y[j];
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
+                                               DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), dshape_x(n1d), dshape_y(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i <= p; i++)
+      {
+         du(o,0) = dshape_x[i]*shape_y[j];
+         du(o,1) = shape_x[i]*dshape_y[j];
+         o += 1;
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
+   const real_t dxb_T = (1.0 - 2*ip.x)*ip.y*(1.0 - ip.y);
+   const real_t dyb_T = ip.x*(1.0 - ip.x)*(1.0 - 2*ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i <= q; i++)
+      {
+         du(o,0) = (dxb_T*shape_x[i] + b_T*dshape_x[i])*shape_y[j];
+         du(o,1) = (dyb_T*shape_y[j] + b_T*dshape_y[j])*shape_x[i];
+         o += 1;
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+H1Bubble_TetrahedronElement::H1Bubble_TetrahedronElement(
+   int p, int q, int btype)
+   : NodalFiniteElement(3, Geometry::TETRAHEDRON,
+                        2*(p*p + 1) + ((q+1)*(q+2)*(q+3))/6,
+                        max(p, 4 + q), FunctionSpace::Pk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 4, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p+1, q+1);
+   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   shape_z.SetSize(n1d);
+   shape_l.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+   dshape_z.SetSize(n1d);
+   dshape_l.SetSize(n1d);
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#else
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
+   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
+   Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
+   Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
+
+   // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
+   int o = 4;
+   for (int i = 1; i < p; i++)  // (0,1)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (1,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (1,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (2,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
+   }
+
+   // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (1,2,3)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (0,3,2)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (0,1,3)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (0,2,1)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
+      }
+   }
+
+   // Interior P_{q+4} nodes
+   for (int k = 1; k < q + 4; k++)
+   {
+      for (int j = 1; j + k < q + 4; j++)
+      {
+         for (int i = 1; i + j + k < q + 4; i++)
+         {
+            real_t w = cp2[i] + cp2[j] + cp2[k] + cp2[q+4-i-j-k];
+            Nodes.IntPoint(o++).Set3(cp2[i]/w, cp2[j]/w, cp2[k]/w);
+         }
+      }
+   }
+
+   DenseMatrix Tt(dof, npq);
+   for (int m = 0; m < dof; ++m)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(m);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, ip.z, shape_z);
+      poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
+
+      o = 0;
+      for (int k = 0; k <= p; k++)
+      {
+         for (int j = 0; j + k <= p; j++)
+         {
+            for (int i = 0; i + j + k <= p; i++)
+            {
+               Tt(m, o++) = shape_x[i]*shape_y[j]*shape_z[k]*shape_l[p-i-j-k];
+            }
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      poly1d.CalcBasis(q, ip.z, shape_z);
+      poly1d.CalcBasis(q, 1. - ip.x - ip.y - ip.z, shape_l);
+      const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
+
+      for (int k = 0; k <= q; k++)
+      {
+         for (int j = 0; j + k <= q; j++)
+         {
+            for (int i = 0; i + j + k <= q; i++)
+            {
+               Tt(m, o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k]*shape_l[q-i-j-k];
+            }
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_TetrahedronElement::CalcShape(const IntegrationPoint &ip,
+                                            Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j + k <= p; j++)
+      {
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            u[o++] = shape_x[i]*shape_y[j]*shape_z[k]*shape_l[p-i-j-k];
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z);
+   poly1d.CalcBasis(q, 1. - ip.x - ip.y - ip.z, shape_l);
+   const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j + k <= q; j++)
+      {
+         for (int i = 0; i + j + k <= q; i++)
+         {
+            u(o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k]*shape_l[q-i-j-k];
+         }
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
+                                             DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p+1, q+1);
+   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
+
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d);
+   Vector dshape_x(n1d), dshape_y(n1d), dshape_z(n1d), dshape_l(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   const real_t lambda = 1.0 - ip.x - ip.y - ip.z;
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
+   poly1d.CalcBasis(p, lambda, shape_l, dshape_l);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j + k <= p; j++)
+      {
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            int l = p - i - j - k;
+            du(o,0) = (dshape_x[i]*shape_l[l] - shape_x[i]*dshape_l[l])
+                      *shape_y[j]*shape_z[k];
+            du(o,1) = (dshape_y[j]*shape_l[l] - shape_y[j]*dshape_l[l])
+                      *shape_x[i]*shape_z[k];
+            du(o,2) = (dshape_z[k]*shape_l[l] - shape_z[k]*dshape_l[l])
+                      *shape_x[i]*shape_y[j];
+            o++;
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z, dshape_z);
+   poly1d.CalcBasis(q, lambda, shape_l, dshape_l);
+   const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
+   const real_t dxb_T = ip.y * ip.z * (lambda - ip.x);
+   const real_t dyb_T = ip.x * ip.z * (lambda - ip.y);
+   const real_t dzb_T = ip.x * ip.y * (lambda - ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j + k <= q; j++)
+      {
+         for (int i = 0; i + j + k <= q; i++)
+         {
+            int l = q - i - j - k;
+            du(o,0) = shape_y[j]*shape_z[k]*(dxb_T*shape_x[i]*shape_l[l]
+                                             + b_T*dshape_x[i]*shape_l[l]
+                                             - b_T*shape_x[i]*dshape_l[l]);
+            du(o,1) = shape_x[i]*shape_z[k]*(dyb_T*shape_y[j]*shape_l[l]
+                                             + b_T*dshape_y[j]*shape_l[l]
+                                             - b_T*shape_y[j]*dshape_l[l]);
+            du(o,2) = shape_x[i]*shape_y[j]*(dzb_T*shape_z[k]*shape_l[l]
+                                             + b_T*dshape_z[k]*shape_l[l]
+                                             - b_T*shape_z[k]*dshape_l[l]);
+            o++;
+         }
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+H1Bubble_HexahedronElement::H1Bubble_HexahedronElement(
+   int p, int q, int btype)
+   : NodalFiniteElement(3, Geometry::CUBE, (2 + 6*p*p) + (q+1)*(q+1)*(q+1),
+                        max(p, 2 + q), FunctionSpace::Qk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 2, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   shape_z.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+   dshape_z.SetSize(n1d);
+
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
+   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
+   Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
+   Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
+
+   Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
+   Nodes.IntPoint(5).Set3(cp[p], cp[0], cp[p]);
+   Nodes.IntPoint(6).Set3(cp[p], cp[p], cp[p]);
+   Nodes.IntPoint(7).Set3(cp[0], cp[p], cp[p]);
+
+   int o = 8;
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);   // (0,1)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);   // (1,2)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);   // (3,2)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);   // (0,3)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[p]);   // (4,5)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[p]);   // (5,6)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[p]);   // (7,6)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[p]);   // (4,7)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);   // (0,4)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[0], cp[i]);   // (1,5)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[p], cp[i]);   // (2,6)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[p], cp[i]);   // (3,7)
+   }
+
+   // faces
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[p-j], cp[0]);   // (3,2,1,0)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[j]);   // (0,1,5,4)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[j]);   // (1,2,6,5)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[p-i], cp[p], cp[j]);   // (2,3,7,6)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[j]);   // (3,0,4,7)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[p]);   // (4,5,6,7)
+      }
+   }
+
+   // interior P_{q+2} nodes
+   for (int k = 1; k < q+2; k++)
+   {
+      for (int j = 1; j < q+2; j++)
+      {
+         for (int i = 1; i < q+2; i++)
+         {
+            Nodes.IntPoint(o++).Set3(cp2[i], cp2[j], cp2[k]);
+         }
+      }
+   }
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d);
+#endif
+
+   DenseMatrix Tt(dof, npq);
+   for (int m = 0; m < dof; ++m)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(m);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, ip.z, shape_z);
+
+      o = 0;
+      for (int k = 0; k <= p; k++)
+      {
+         for (int j = 0; j <= p; j++)
+         {
+            for (int i = 0; i <= p; i++)
+            {
+               Tt(m, o++) = shape_x[i]*shape_y[j]*shape_z[k];
+            }
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      poly1d.CalcBasis(q, ip.z, shape_z);
+      const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+      for (int k = 0; k <= q; k++)
+      {
+         for (int j = 0; j <= q; j++)
+         {
+            for (int i = 0; i <= q; i++)
+            {
+               Tt(m, o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k];
+            }
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_HexahedronElement::CalcShape(const IntegrationPoint &ip,
+                                           Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i <= p; i++)
+         {
+            u(o++) = shape_x[i]*shape_y[j]*shape_z[k];
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i <= q; i++)
+         {
+            u(o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k];
+         }
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
+                                            DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), dshape_x(n1d),
+          dshape_y(n1d), dshape_z(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i <= p; i++)
+         {
+            du(o,0) = dshape_x[i]*shape_y[j]*shape_z[k];
+            du(o,1) = shape_x[i]*dshape_y[j]*shape_z[k];
+            du(o,2) = shape_x[i]*shape_y[j]*dshape_z[k];
+            o += 1;
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z, dshape_z);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+   const real_t dxb_T = (1.0 - 2*ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+   const real_t dyb_T = ip.x*(1.0 - ip.x)*(1.0 - 2*ip.y)*ip.z*(1.0 - ip.z);
+   const real_t dzb_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*(1.0 - 2*ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i <= q; i++)
+         {
+            du(o,0) = (dxb_T*shape_x[i] + b_T*dshape_x[i])*shape_y[j]*shape_z[k];
+            du(o,1) = (dyb_T*shape_y[j] + b_T*dshape_y[j])*shape_x[i]*shape_z[k];
+            du(o,2) = (dzb_T*shape_z[k] + b_T*dshape_z[k])*shape_x[i]*shape_y[j];
+            o += 1;
+         }
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+}
@@ -0,0 +1,109 @@
+// 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_FE_H1_BUBBLE
+#define MFEM_FE_H1_BUBBLE
+
+#include "fe_base.hpp"
+
+namespace mfem
+{
+
+/// Arbitrary order H1 plus bubble elements in 2D on a triangle
+class H1Bubble_TriangleElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the triangular bubble element with degree-p polynomials,
+   /// enriched with cubic bubble times degree q polynomial.
+   H1Bubble_TriangleElement(int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+/// Arbitrary order H1 plus bubble elements in 2D on a quadrilateral
+class H1Bubble_QuadrilateralElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, dshape_x, dshape_y, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the quadrilateral bubble element with degree-p
+   /// polynomials, enriched with biquadratic bubble times degree q polynomial.
+   H1Bubble_QuadrilateralElement(
+      int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+/// Arbitrary order H1 plus bubble elements in 3D on a tetrahedron
+class H1Bubble_TetrahedronElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z, shape_l;
+   mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the tetrahedral bubble element with degree-p
+   /// polynomials, enriched with quartic bubble times degree q polynomial.
+   H1Bubble_TetrahedronElement(int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+/// Arbitrary order H1 plus bubble elements in 3D on a hexahedron
+class H1Bubble_HexahedronElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z;
+   mutable Vector dshape_x, dshape_y, dshape_z, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the hexahedral bubble element with degree-p polynomials,
+   /// enriched with triquadratic bubble times degree q polynomial.
+   H1Bubble_HexahedronElement(int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+} // namespace mfem
+
+#endif
@@ -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)
@@ -663,6 +663,9 @@ public:
                         const int cb_type = BasisType::GaussLobatto,
                         const int ob_type = BasisType::GaussLegendre);

+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { 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 { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 3; }
+
   using FiniteElement::CalcVShape;
   using FiniteElement::CalcPhysCurlShape;

@@ -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)
@@ -510,6 +510,9 @@ public:
   RT_R2D_SegmentElement(const int p,
                         const int ob_type = BasisType::GaussLegendre);

+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { 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 { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 0; }
+
   using FiniteElement::CalcVShape;

   void CalcVShape(ElementTransformation &Trans,
@@ -243,11 +243,21 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
   {
      fec = new H1Ser_FECollection(atoi(name + 10), atoi(name + 6));
   }
+   else if (!strncmp(name, "H1Bubble_", 9))
+   {
+      fec = new H1Bubble_FECollection(atoi(name + 13), atoi(name + 16),
+                                      atoi(name + 9));
+   }
   else if (!strncmp(name, "H1@", 3))
   {
      fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
                                BasisType::GetType(name[3]));
   }
+   else if (!strncmp(name, "H1Bubble@", 9))
+   {
+      fec = new H1Bubble_FECollection(atoi(name + 15), atoi(name + 18),
+                                      atoi(name + 11), BasisType::GetType(name[9]));
+   }
   else if (!strncmp(name, "L2_T", 4))
      fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
                                atoi(name + 4));
@@ -2122,6 +2132,173 @@ H1_FECollection::~H1_FECollection()
   }
 }

+static int GetBubbleSpaceOrder(int p, int q, int dim)
+{
+   switch (dim)
+   {
+      case 0: return 0;
+      case 1: return std::max(p, q + 2);
+      case 2: return std::max(p, q + 3);
+      case 3: return std::max(p, q + 4);
+   }
+   MFEM_ABORT("Unsupported dimension.");
+}
+
+H1Bubble_FECollection::H1Bubble_FECollection(const int p, const int q,
+                                             const int dim, const int btype)
+   : FiniteElementCollection(GetBubbleSpaceOrder(p, q, dim)),
+     dim(dim),
+     b_type(BasisType::Check(btype)),
+     h1_order(p),
+     bubble_order(q)
+{
+   MFEM_VERIFY(p >= 1, "H1Bubble_FECollection requires order >= 1.");
+   MFEM_VERIFY(dim >= 0 && dim <= 3, "Unsupported dimension.");
+
+   switch (btype)
+   {
+      case BasisType::GaussLobatto:
+      {
+         snprintf(fec_name, 32, "H1Bubble_%dD_P%d_P%d", dim, p, q);
+         break;
+      }
+      default:
+      {
+         const int pt_type = BasisType::GetQuadrature1D(btype);
+         MFEM_VERIFY(Quadrature1D::CheckClosed(pt_type) != Quadrature1D::Invalid,
+                     "unsupported BasisType: " << BasisType::Name(btype));
+         snprintf(fec_name, 32, "H1Bubble@%c_%dD_P%d_P%d",
+                  (int)BasisType::GetChar(btype), dim, p, q);
+      }
+   }
+
+   dofs[Geometry::POINT] = 1;
+   elements[Geometry::POINT] = make_unique<PointFiniteElement>();
+
+   if (dim >= 1)
+   {
+      dofs[Geometry::SEGMENT] = p - 1;
+      elements[Geometry::SEGMENT] = make_unique<H1_SegmentElement>(p, btype);
+   }
+
+   if (dim == 2)
+   {
+      dofs[Geometry::TRIANGLE] = ((q+1)*(q+2))/2;
+      dofs[Geometry::SQUARE] = (q+1)*(q+1);
+
+      elements[Geometry::TRIANGLE] =
+         make_unique<H1Bubble_TriangleElement>(p, q, btype);
+      elements[Geometry::SQUARE] =
+         make_unique<H1Bubble_QuadrilateralElement>(p, q, btype);
+   }
+
+   if (dim == 3)
+   {
+      dofs[Geometry::TRIANGLE] = ((p-1)*(p-2))/2;
+      dofs[Geometry::SQUARE] = (p-1)*(p-1);
+      dofs[Geometry::TETRAHEDRON] = ((q+1)*(q+2)*(q+3))/6;
+      dofs[Geometry::CUBE] = (q+1)*(q+1)*(q+1);
+
+      elements[Geometry::TRIANGLE] = make_unique<H1_TriangleElement>(p, btype);
+      elements[Geometry::SQUARE] = make_unique<H1_QuadrilateralElement>(p, btype);
+
+      elements[Geometry::TETRAHEDRON] =
+         make_unique<H1Bubble_TetrahedronElement>(p, q, btype);
+      elements[Geometry::CUBE] =
+         make_unique<H1Bubble_HexahedronElement>(p, q, btype);
+   }
+
+   // DOF orderings. Need only for lower-dimensional entities.
+   // Segment DOF orderings in 2D.
+   if (dim >= 2)
+   {
+      seg_dof_ord[0].resize(p - 1);
+      seg_dof_ord[1].resize(p - 1);
+      for (int i = 0; i < p - 1; i++)
+      {
+         seg_dof_ord[0][i] = i;
+         seg_dof_ord[1][i] = p - 2 - i;
+      }
+   }
+
+   // Face (triangle or quadrilateral) DOF orderings in 3D.
+   if (dim == 3)
+   {
+      const int n_tri_dof = dofs[Geometry::TRIANGLE];
+      for (int i = 0; i < 6; i++)
+      {
+         tri_dof_ord[i].resize(n_tri_dof);
+      }
+      // see Mesh::GetTriOrientation in mesh/mesh.cpp
+      const int pm1 = p - 1;
+      const int pm2 = p - 2;
+      for (int j = 0; j < pm2; j++)
+      {
+         for (int i = 0; i + j < pm2; i++)
+         {
+            int o = n_tri_dof - ((pm1 - j)*(pm2 - j))/2 + i;
+            int k = (p - 3) - j - i;
+            tri_dof_ord[0][o] = o;  // (0,1,2)
+            tri_dof_ord[1][o] = n_tri_dof - ((pm1-j)*(pm2-j))/2 + k;  // (1,0,2)
+            tri_dof_ord[2][o] = n_tri_dof - ((pm1-i)*(pm2-i))/2 + k;  // (2,0,1)
+            tri_dof_ord[3][o] = n_tri_dof - ((pm1-k)*(pm2-k))/2 + i;  // (2,1,0)
+            tri_dof_ord[4][o] = n_tri_dof - ((pm1-k)*(pm2-k))/2 + j;  // (1,2,0)
+            tri_dof_ord[5][o] = n_tri_dof - ((pm1-i)*(pm2-i))/2 + j;  // (0,2,1)
+         }
+      }
+
+      const int n_quad_dof = dofs[Geometry::SQUARE];
+      for (int i = 0; i < 8; i++)
+      {
+         quad_dof_ord[i].resize(n_quad_dof);
+      }
+      for (int j = 0; j < pm1; j++)
+      {
+         for (int i = 0; i < pm1; i++)
+         {
+            int o = i + j*pm1;
+            quad_dof_ord[0][o] = i + j*pm1;  // (0,1,2,3)
+            quad_dof_ord[1][o] = j + i*pm1;  // (0,3,2,1)
+            quad_dof_ord[2][o] = j + (pm2 - i)*pm1;  // (1,2,3,0)
+            quad_dof_ord[3][o] = (pm2 - i) + j*pm1;  // (1,0,3,2)
+            quad_dof_ord[4][o] = (pm2 - i) + (pm2 - j)*pm1;  // (2,3,0,1)
+            quad_dof_ord[5][o] = (pm2 - j) + (pm2 - i)*pm1;  // (2,1,0,3)
+            quad_dof_ord[6][o] = (pm2 - j) + i*pm1;  // (3,0,1,2)
+            quad_dof_ord[7][o] = i + (pm2 - j)*pm1;  // (3,2,1,0)
+         }
+      }
+   }
+}
+
+const FiniteElement *
+H1Bubble_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   return elements[GeomType].get();
+}
+
+const int *H1Bubble_FECollection::DofOrderForOrientation(
+   Geometry::Type GeomType, int Or) const
+{
+   if (GeomType == Geometry::SEGMENT)
+   {
+      return (Or > 0) ? seg_dof_ord[0].data() : seg_dof_ord[1].data();
+   }
+   else if (GeomType == Geometry::TRIANGLE)
+   {
+      return tri_dof_ord[Or%6].data();
+   }
+   else if (GeomType == Geometry::SQUARE)
+   {
+      return quad_dof_ord[Or%8].data();
+   }
+   return nullptr;
+}
+
+FiniteElementCollection *H1Bubble_FECollection::GetTraceCollection() const
+{
+   return (dim < 0) ? NULL : new H1_Trace_FECollection(h1_order, dim, b_type);
+}
+

 H1_Trace_FECollection::H1_Trace_FECollection(const int p, const int dim,
                                             const int btype)
@@ -111,6 +111,8 @@ public:
   | :------: | :---: | :---: | :-------: | :-----: | :---: |
   | H1_[DIM]_[ORDER] | H1 | * | 1 | VALUE | H1 nodal elements |
   | H1@[BTYPE]_[DIM]_[ORDER] | H1 | * | * | VALUE | H1 nodal elements |
+   | H1Bubble_[DIM]_[ORDER]_[BUBBLE_ORDER] | H1 | * | 1 | VALUE | H1 nodal elements enriched with bubble functions |
+   | H1Bubble@[BTYPE]_[DIM]_[ORDER]_[BUBBLE_ORDER] | H1 | * | 1 | VALUE | H1 nodal elements enriched with bubble functions |
   | H1Pos_[DIM]_[ORDER] | H1 | * | 2 | VALUE | H1 nodal elements |
   | H1Pos_Trace_[DIM]_[ORDER] | H^{1/2} | * | 2 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
   | H1_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
@@ -317,6 +319,59 @@ public:
   virtual ~H1_FECollection();
 };

+/// @brief Arbitrary order $H^1$-conforming (continuous) finite elements
+/// enriched with bubble functions.
+///
+/// The bubble space consists of the standard $P_p$ or $Q_p$ space, enriched
+/// with bubble functions, which are degree-$q$ polynomials times $b$, where $b$
+/// is the lowest-order bubble function.
+///
+/// The traces are the same as the standard $H^1$ traces.
+class H1Bubble_FECollection : public FiniteElementCollection
+{
+protected:
+   int dim;
+   int b_type;
+   int h1_order;
+   int bubble_order;
+
+   char fec_name[32];
+   std::array<int, Geometry::NumGeom> dofs{}; // zero initialize
+   std::array<std::unique_ptr<FiniteElement>, Geometry::NumGeom> elements;
+
+   std::array<std::vector<int>, 2> seg_dof_ord;
+   std::array<std::vector<int>, 6> tri_dof_ord;
+   std::array<std::vector<int>, 8> quad_dof_ord;
+   std::array<std::vector<int>, 24> tet_dof_ord;
+
+public:
+   /// Construct the $H^1$ bubble collection consisting of degree-$p$
+   /// polynomials enriched with the bubble function times degree-$q$
+   /// polynomials.
+   explicit H1Bubble_FECollection(const int p, const int q, const int dim = 3,
+                                  const int btype = BasisType::GaussLobatto);
+
+   const FiniteElement *
+   FiniteElementForGeometry(Geometry::Type GeomType) const override;
+
+   int DofForGeometry(Geometry::Type GeomType) const override
+   { return dofs[GeomType]; }
+
+   const int *DofOrderForOrientation(Geometry::Type GeomType,
+                                     int Or) const override;
+
+   const char *Name() const override { return fec_name; }
+
+   int GetContType() const override { return CONTINUOUS; }
+
+   int GetBasisType() const { return b_type; }
+
+   FiniteElementCollection *GetTraceCollection() const override;
+
+   FiniteElementCollection *Clone(int p) const override
+   { return new H1Bubble_FECollection(p, bubble_order, dim, b_type); }
+};
+
 /** @brief Arbitrary order H1-conforming (continuous) finite elements with
    positive basis functions. */
 class H1Pos_FECollection : public H1_FECollection
@@ -3877,9 +3877,12 @@ const FiniteElement *FiniteElementSpace::GetFE(int i) const
   else
   {
 #ifdef MFEM_DEBUG
-      // consistency check: fec->GetOrder() and FE->GetOrder() should return
-      // the same value (for standard, constant-order spaces)
-      if (!IsVariableOrder() && FE->GetDim() > 0)
+      // Consistency check: fec->GetOrder() and FE->GetOrder() should return
+      // the same value (for standard, constant-order spaces). Skip this check
+      // even for constant-order bubble spaces, since the bubble functions on
+      // different geometries have different orders.
+      if (!IsVariableOrder() && FE->GetDim() > 0 &&
+          dynamic_cast<const H1Bubble_FECollection*>(fec) == nullptr)
      {
         MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
                     "internal error: " <<
@@ -3137,52 +3137,29 @@ 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++)
+   if (fes->GetNBE() > 0)
   {
-      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);
+      // TODO: Replace this by GetTypicalBdrElement() once implemented
+      const FiniteElement *be = fes->GetBE(0);
+      MFEM_VERIFY(be->GetRangeType() == FiniteElement::SCALAR &&
+                  be->GetMapType() == FiniteElement::INTEGRAL, "Not an RT FE space!");
   }
-#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 +3175,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 +4991,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 +5049,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);
@@ -532,6 +532,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 +659,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 +1928,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 +1936,52 @@ 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;
+
+   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);

@@ -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;
@@ -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:
@@ -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);
            }
         }
@@ -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
         {
@@ -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
 }
@@ -278,9 +278,18 @@ void ArraysByName<T>::Load(std::istream &in)
         q1 = ArrayLine.find(' ');
         ArrayName = ArrayLine.substr(0,q1-1);
      }
-
-      // Ignore the remainder of the line which may contain explanatory comments
-      data[ArrayName].Load(in, 0);
+      if (q1+2 < ArrayLine.size())
+      {
+         // Read the remainder of the line which contains the array data
+         std::istringstream ArrayDataStream(ArrayLine.substr(q1+2,
+                                                             ArrayLine.size()));
+         data[ArrayName].Load(ArrayDataStream, 0);
+      }
+      else
+      {
+         // Read the array data starting on the next line
+         data[ArrayName].Load(in, 0);
+      }
   }

 }
@@ -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();
@@ -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();
@@ -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);
      }
   }
@@ -5329,21 +5322,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 +5355,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 +5609,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 +5703,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);
@@ -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
 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)

 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)

@@ -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
@@ -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))
   {
@@ -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)
@@ -0,0 +1,25 @@
+# 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()
@@ -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 ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/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
+
+.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)
+
+# 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)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+
@@ -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
+
@@ -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

 ////////////////////////////////////////////////////////////////////////////////

@@ -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
@@ -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))
@@ -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);
+         }
+      }
+   }
+}
+
+
@@ -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++)
@@ -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);
+}
Author	SHA1	Message	Date
Will Pazner	cc7ebccc54	Fix signed char issue in Device::GetUUID	2026-04-17 11:38:22 -04:00
adam-sim-dev	9e423f2f8e	Fixed missing parenthesis in the comment	2026-04-17 11:38:22 -04:00
Veselin Dobrev	c14938cd1c	Adjust seed values in sample runs in ex12p to ensure LOBPCG convergence in older hypre versions.	2026-04-17 11:38:22 -04:00
Hugh Carson	419758be1e	Address PR feedback - Use [IntegrationRules] test tag instead of [PositiveWeightRules] - Remove redundant case 21: (default branch handles it via the overwrite guard) - Remove trailing blank line	2026-04-17 11:38:22 -04:00
Hugh Carson	4444d8ed70	Remove unused private helper methods from IntegrationRule AddTriPoints3R, AddTetPoints4b, and AddTetPoints12bc are no longer called after the legacy simplex rules were removed.	2026-04-17 11:38:22 -04:00
Hugh Carson	48670b9f87	Use exact fractions for trivial quadrature weights and coordinates For rules where the mathematical value is an exact simple fraction (midpoint weights, equal-weight symmetric rules), use the fraction directly rather than the Polyquad decimal expansion. Cleaner to read and avoids any rounding from decimal-to-double conversion.	2026-04-17 11:38:22 -04:00
Hugh Carson	bddb52ace5	Restore original function order in intrules.cpp Move TriangleIntegrationRule before SquareIntegrationRule to match the original file layout, reducing diff noise against master.	2026-04-17 11:38:22 -04:00
Hugh Carson	2d99e1e2de	Remove legacy simplex rules; positive-weight rules are now the default The positive-weight rules now cover the full tabulated range for both triangles (0-25) and tetrahedra (0-20), so the old rules with negative weights are no longer needed. Remove the SimplexQuadrature enum, simplex_type member, and legacy rule functions — all simplex quadrature now uses positive-weight rules by default, with Grundmann-Moller fallback for higher orders.	2026-04-17 11:38:22 -04:00
Hugh Carson	fdad993654	Add existing order 21-25 triangle rule to positive-weight rules The 126-point degree-25 rule already has all positive weights. Copy it into TrianglePositiveIntegrationRule so the positive-weight path covers orders 0-25.	2026-04-17 11:38:22 -04:00
Hugh Carson	7f51024345	Fix memory leak	2026-04-17 11:38:22 -04:00
Hugh Carson	0d00e79cd5	Add positive-weight simplex quadrature rules for orders 0-20 Triangle rules from Witherden & Vincent (2015), tet rules d=0-13 from Witherden & Vincent, tet rules d=14-20 from Chuluunbaatar et al. (2022). All rules have strictly positive weights and interior points, replacing the legacy rules which use negative weights at several orders and fall back to Grundmann-Moller (negative weights, high point counts) for tets at d>=9.	2026-04-17 11:38:22 -04:00
Gabriele Bozzola	b406cdaf79	Improve error message for gmsh versions != 2.2 I am a new user of [palace](https://github.com/awslabs/palace). As I was trying to set a simple mesh up (with gmsh), I kept getting indexing errors that I could not decipher. I eventually [learned](https://mfem.org/mesh-formats/) that supported version for gmsh meshes is 2.2. This commit catches this and adds an informative error.	2026-04-17 11:38:22 -04:00
chapman39	8adba4e1bb	add comments showing each modulus replacement	2026-04-17 11:38:22 -04:00
chapman39	2b86c7300b	added comment	2026-04-17 11:38:22 -04:00
chapman39	7357a9b4bf	eliminate usage of modulus to avoid llvm backend bug	2026-04-17 11:38:22 -04:00
Veselin Dobrev	b44728af9e	Update the action actions/cache/restore to v5	2026-04-17 11:38:22 -04:00
Veselin Dobrev	8d002d09c8	Updated the github/codeql-action/* actions to the latest, v4	2026-04-17 11:38:22 -04:00
Veselin Dobrev	fed8e6bc1b	Updated actions/checkout to the latest major version, v6	2026-04-17 11:38:22 -04:00
Veselin Dobrev	0aa0ac0637	Update actions/{checkout,cache} to v5 Update github/codeql-action/* to v3	2026-04-17 11:38:22 -04:00
Jan Nikl	77ef843c2e	Added scalar unit test of ProjectBdrCoefficientNormal().	2026-04-17 11:38:22 -04:00
Jan Nikl	51bc8037d4	Added a unit test for vector ProjectBdrCoefficientNormal().	2026-04-17 11:38:22 -04:00
Jan Nikl	523c208d87	Made the ProjectBdrCoefficientNormal check non-debug.	2026-04-17 11:38:22 -04:00
Jan Nikl	90e0e8e289	Minor unification of docstrings.	2026-04-17 11:38:22 -04:00
Jan Nikl	3f41665e4f	Generalized RT normal projection.	2026-04-17 11:38:22 -04:00
Jan Nikl	6f07de9114	Removed unused code.	2026-04-17 11:38:22 -04:00
Jan Nikl	9f74ee130a	Fixed vis of the initial exact solution.	2026-04-17 11:38:22 -04:00
Jan Nikl	0a7eb2c39e	Fixed visulization in ex22p.	2026-04-17 11:38:22 -04:00
Jan Nikl	d85723ce29	Fixed spelling of transverse.	2026-04-17 11:38:22 -04:00
Jan Nikl	b33340edab	Added documentation and checks to the extrusion classes.	2026-04-17 11:38:22 -04:00
Jan Nikl	9c1bf9704d	Added extrusion of vector 1D grid functions.	2026-04-17 11:38:22 -04:00
Stowell, Mark L.	c58816905d	Adding bugfix and unit test which would have caught the bug	2026-04-17 11:38:22 -04:00
Wouter Tonnon	ded65cf364	extended to MixedBilinearForm	2026-04-17 11:38:22 -04:00
Wouter Tonnon	422f42ec0d	added missing face orientation	2026-04-17 11:38:22 -04:00
Veselin Dobrev	acac245260	Small change in error message + formatting.	2026-04-17 11:38:22 -04:00
thartland	ed0b39b732	VERIFY instead of ASSERT	2026-04-17 11:38:22 -04:00
Tucker Hartland	d609bee2cc	style	2026-04-17 11:38:22 -04:00
thartland	421f3f03ba	adding a check to make sure that each process owns at least one entry of the HypreParVector prior to calling GlobalVector	2026-04-17 11:38:22 -04:00
Will Pazner	842f88a7b3	Use constexpr in unit test	2026-04-17 11:38:22 -04:00
Will Pazner	abf97587d4	Add comment about the shape of FaceNbrData	2026-04-17 11:38:22 -04:00
Will Pazner	355a2cd570	Add unit test for parallel L2 face restriction with vdim > 1	2026-04-17 11:38:22 -04:00
Will Pazner	e7e184a24d	Fix bug in ParL2FaceRestriction with vdim > 1 The layout of the FaceNbrData vector was not handled properly	2026-04-17 11:38:22 -04:00
Stowell, Mark L.	616eaec18e	Updating unit tests	2026-04-17 11:38:22 -04:00
Stowell, Mark L.	5c0711f334	Using new MapType entries and implementing new GetPhys*Dim methods	2026-04-17 11:38:22 -04:00
Stowell, Mark L.	999e4c4f46	Adding new MapType entries for R2D and R1D classes	2026-04-17 11:38:22 -04:00
Barry Smith Nuno Nobre Satish Balay	214750edc8	Update to use PetscCtxRt from (3,25,0), and cleanup duplicate code Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk> Co-authored-by: Satish Balay <balay@mcs.anl.gov>	2026-04-17 11:38:22 -04:00
Satish Balay	4cf708d4fa	update KSPMonitorFn usage for < (3,24,0)	2026-04-17 11:38:22 -04:00
Satish Balay	91a40a1d1a	update PetscCtxDestroyFn usage for < (3,23,0)	2026-04-17 11:38:22 -04:00
Satish Balay	dffe36f382	rework PetscContainerSetCtxDestroy() usage for < (3,23,0)	2026-04-17 11:38:22 -04:00
chapman39	49363859ee	80 chars/ line	2026-04-17 11:38:22 -04:00
chapman39	9f47892f62	dfem integrate: use mfem abort kernel in device code	2026-04-17 11:38:22 -04:00
Stowell, Mark L.	1c286184be	Changing copyright date to pass CI checks	2026-04-17 11:38:22 -04:00
Stowell, Mark L.	6657cf2760	Adding miniapps/plasma subdirectory to build system	2026-04-17 11:38:22 -04:00
Stowell, Mark L.	731224d5e8	Adding plasma miniapp directory	2026-04-17 11:38:22 -04:00
Jan Nikl	a8a85c68fb	Minor docstring correction.	2026-04-17 11:38:22 -04:00
Jan Nikl	591cc1ca41	Fixed complex grid function copy assignment.	2026-04-17 11:38:22 -04:00
Andrew Ho	6b1c2644e6	comment on why TPL_LIBRARIES is reversed twice	2026-04-17 11:38:22 -04:00
Andrew HoandNuno Nobre	7742ad8355	Update config/cmake/modules/MfemCmakeUtilities.cmake Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2026-04-17 11:38:22 -04:00
Andrew HoandNuno Nobre	275e98264c	Update config/cmake/modules/MfemCmakeUtilities.cmake Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2026-04-17 11:38:21 -04:00
Andrew Ho	cd4d7c292f	move cudart to MFEM_EXT_LIBS	2026-04-17 11:38:21 -04:00
Andrew Ho	e130ae7dd8	fixed wrong dir being marked as system	2026-04-17 11:38:21 -04:00
Andrew HoandNuno Nobre	e026e15fe3	Update config/cmake/modules/MfemCmakeUtilities.cmake Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2026-04-17 11:38:21 -04:00
Andrew Ho	493b5a942e	MFEM_EXPORT_GPU_CONFIG should export CPU config.mk when set to off	2026-04-17 11:38:21 -04:00
Andrew Ho	f27a13cbad	revert change, updated comment to why libdl gets special treatment	2026-04-17 11:38:21 -04:00
Andrew Ho	559d0e42c7	suggestions from Veselin	2026-04-17 11:38:21 -04:00
Andrew Ho	5ba64a774e	missed one old unsetting of shared_link_flag	2026-04-17 11:38:21 -04:00
Andrew Ho	896e251d3a	review suggestions	2026-04-17 11:38:21 -04:00
Andrew Ho	9b0c9d3f6e	fixes for hip	2026-04-17 11:38:21 -04:00
Andrew Ho	c90d6f9d60	remove debug printout	2026-04-17 11:38:21 -04:00
Andrew Ho	6f2b8b82d1	seems to be building external laghos now	2026-04-17 11:38:21 -04:00
Andrew Ho	90cf6af2bb	improving config.mk file generated by cmake to work with hip/cuda Still need to export compiler flags	2026-04-17 11:38:21 -04:00
jdongg	cd6bfc7de8	fix clang compiler warnings from origin/catch-tests	2026-03-06 14:07:34 -08:00
Will Pazner	b42ad0a57c	Merge remote-tracking branch 'origin/master' into bubble # Conflicts: # fem/fe_coll.hpp	2026-03-01 16:30:46 -08:00
Dohyun Kim	8e67185297	Merge branch 'master' into bubble	2026-01-03 01:58:22 +09:00
Will Pazner	456c236cc5	Small fixes Add local variables in thread-safe mode Fix MFEM_VERIFY message Fix trace collection order	2025-12-05 11:10:57 -08:00
Will Pazner	f80902b776	Re-add assertion; skip check for bubble spaces	2025-12-05 10:14:06 -08:00
Will Pazner	4c9f6edef0	Improve Doxygen	2025-12-05 10:14:06 -08:00
Will Pazner	45ec9d451d	Support "H1Bubble@" in FiniteElementCollection::New	2025-12-05 10:14:06 -08:00
Dohyun Kim	988ab81e5f	FEColl::New	2025-12-05 10:14:06 -08:00
Will Pazner	3a7b1d7c67	Use bubble elements in ex36 and ex36p	2025-12-05 10:14:06 -08:00
Will Pazner	411ffcfef6	Fix DOF orderings in bubble elements	2025-12-05 10:14:06 -08:00
Will Pazner	27d79fc463	Revert "Return nullptr for H1Bubble_FECollection::DofOrderForOrientation" This reverts commit 58e7bb6e6e213eb90839feb67de7d0e03b5799da.	2025-12-05 10:14:06 -08:00
Will Pazner	898337f772	Return nullptr for H1Bubble_FECollection::DofOrderForOrientation Some features (e.g. node reordering) won't be supported; this could be added later.	2025-12-05 10:14:06 -08:00
Will Pazner	d3063a0982	Add bubble tets and hexes	2025-12-05 10:14:06 -08:00
Will Pazner	15853215ff	Disable check that FE and FEC orders are the same With enriched bubble elements, the orders could be different. For example, linear triangle enriched with bubble has max total degree 3, but the linear quadrilateral enriched with bubble has max degree 2 in each variable (and max total degree 4).	2025-12-05 10:14:06 -08:00
Will Pazner	8b0e9ff064	Move bubble elements to their own file	2025-12-05 10:14:06 -08:00
Will Pazner	03973ad244	Add quad bubble element, change meaning of q	2025-12-05 09:57:09 -08:00
Will Pazner	88c04a6e45	Add H1 bubble triangle element and collection	2025-12-05 09:57:09 -08:00