style

This has been superseded by the MFEMSidreDataCollection in Axom.

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

@@ -295,7 +295,8 @@ jobs:
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
           brew install enzyme
-          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
+          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required.*:.*\(llvm[^ ]*\).*/\1/p')
+          echo "ENZYME_LLVM=$ENZYME_LLVM"
           LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
           echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
           echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV

CHANGELOG

@@ -23,7 +23,6 @@ Discretization improvements
   Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
   2022.
 
-
 Version 4.9.1 (development)
 ===========================
 
@@ -46,6 +45,11 @@ New and updated examples and miniapps
 - Electromagnetics/lorentz miniapp has been updated to leverage the ParticleSet
   capability.
 
+Miscellaneous
+-------------
+- Removes the SidreDataCollection class from MFEM in favor of the 
+  MFEMSidreDataCollection class in the Axom library (https://github.com/llnl/axom).
+
 
 Version 4.9, released on Dec 11, 2025
 =====================================

CMakeLists.txt

@@ -75,12 +75,10 @@ set(XSDK_ENABLE_Fortran OFF)
 
 # Check if we need to enable C or Fortran.
 if (MFEM_USE_CONDUIT OR
-    MFEM_USE_SIDRE OR
     MFEM_USE_PETSC)
    # This seems to be needed by:
    #  * find_package(BLAS REQUIRED) and
-   #  * find_package(HDF5 REQUIRED) needed, in turn, by:
-   #    - find_package(AXOM REQUIRED)
+   #  * find_package(HDF5 REQUIRED) and
    #  * find_package(PETSc REQUIRED)
    set(XSDK_ENABLE_C ON)
 endif()
@@ -478,11 +476,6 @@ if (MFEM_USE_FMS)
     find_package(FMS REQUIRED fms)
 endif()
 
-# Axom/Sidre
-if (MFEM_USE_SIDRE)
-  find_package(Axom REQUIRED Axom)
-endif()
-
 # PUMI
 if (MFEM_USE_PUMI)
   # If PUMI_DIR was specified, only link to that directory,
@@ -629,7 +622,7 @@ find_package(Threads REQUIRED)
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
 set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
-    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
+    SUNDIALS PETSC SLEPC MUMPS FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
     NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
     ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CALIPER CODIPACK
     BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ

INSTALL

@@ -452,13 +452,6 @@ MFEM_USE_MPFR = YES/NO
    quadrature rules. When enabled, this option uses the MPFR_* library options,
    see below.
 
-MFEM_USE_SIDRE = YES/NO
-   Sidre is a component of LLNL's axom project, https://github.com/LLNL/axom,
-   that provides an HDF5-based file format for visualization or restart
-   capability following the Conduit (https://github.com/LLNL/conduit) mesh
-   blueprint specification. When enabled, this option requires installation of
-   HDF5 (see also MFEM_USE_NETCDF), Conduit and LLNL's axom project.
-
 MFEM_USE_SIMD = YES/NO
    Enables the high performance templated classes to use architecture dependent
    SIMD intrinsics instead of the generic implementation of class AutoSIMD in
@@ -778,14 +771,6 @@ The specific libraries and their options are:
   Options: SLEPC_OPT, SLEPC_LIB.
   Versions: SLEPc >= 3.8.0.
 
-- Sidre (optional), part of LLNL's axom project, used when MFEM_USE_SIDRE = YES.
-  Starting with MFEM v4.1, Axom version 0.3.1 or later is required.
-  URL: https://github.com/LLNL/axom
-       https://github.com/LLNL/conduit (Conduit)
-       https://support.hdfgroup.org/HDF5 (HDF5)
-  Options: SIDRE_OPT, SIDRE_LIB.
-  Versions: Axom >= 0.3.1.
-
 - Conduit (optional), used when MFEM_USE_CONDUIT = YES. Conduit Mesh Blueprint
   support requires Conduit >= v0.3.1 and VisIt >= v2.13.1 to read the output.
   URL: https://github.com/LLNL/conduit (Conduit)
@@ -1069,7 +1054,6 @@ MFEM_USE_OCCA
 MFEM_USE_CEED
 MFEM_USE_RAJA
 MFEM_USE_UMPIRE
-MFEM_USE_SIDRE
 MFEM_USE_MOONOLITH
 MFEM_USE_CALIPER
 MFEM_USE_FMS
@@ -1133,7 +1117,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - OCCA
  - RAJA
  - UMPIRE
- - AXOM - Used when MFEM_USE_SIDRE is enabled
  - MOONOLITH
  - CALIPER
  - FMS

config/XSDKDefaults.cmake

@@ -248,10 +248,6 @@ IF (DEFINED TPL_ENABLE_MPFR)
   SET(MFEM_USE_MPFR ${TPL_ENABLE_MPFR} CACHE BOOL "Enable MPFR usage." FORCE)
 ENDIF()
 
-IF (DEFINED TPL_ENABLE_SIDRE)
-  SET(MFEM_USE_SIDRE ${TPL_ENABLE_SIDRE} CACHE BOOL "Enable Axom/Sidre usage" FORCE)
-ENDIF()
-
 IF (DEFINED TPL_ENABLE_FMS)
   SET(MFEM_USE_FMS ${TPL_ENABLE_FMS} CACHE BOOL "Enable FMS usage" FORCE)
 ENDIF()

config/cmake/MFEMConfig.cmake.in

@@ -46,7 +46,6 @@ set(MFEM_USE_NETCDF @MFEM_USE_NETCDF@)
 set(MFEM_USE_PETSC @MFEM_USE_PETSC@)
 set(MFEM_USE_SLEPC @MFEM_USE_SLEPC@)
 set(MFEM_USE_MPFR @MFEM_USE_MPFR@)
-set(MFEM_USE_SIDRE @MFEM_USE_SIDRE@)
 set(MFEM_USE_FMS @MFEM_USE_FMS@)
 set(MFEM_USE_CONDUIT @MFEM_USE_CONDUIT@)
 set(MFEM_USE_PUMI @MFEM_USE_PUMI@)

config/cmake/config.hpp.in

@@ -120,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 #cmakedefine MFEM_USE_GNUTLS
 
-// Enable Sidre support.
-#cmakedefine MFEM_USE_SIDRE
-
 // Enable the use of SIMD in the high performance templated classes.
 #cmakedefine MFEM_USE_SIMD
 

config/cmake/modules/FindAdiak.cmake

@@ -0,0 +1,24 @@
+# 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.
+
+# Defines the following variables:
+#   - ADIAK_FOUND
+#   - ADIAK_LIBRARIES
+#   - ADIAK_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+
+mfem_find_package(Adiak ADIAK ADIAK_DIR
+      "include" "adiak.h"
+      "lib" "adiak"
+      "Paths to headers required by Adiak."
+      "Libraries required by Adiak.")
+

config/cmake/modules/FindAxom.cmake

@@ -13,6 +13,9 @@
 #   - AXOM_FOUND
 #   - AXOM_LIBRARIES
 #   - AXOM_INCLUDE_DIRS
+#
+# MFEM itself does not depend on Axom, however Tribol does. This module exists
+# to support MFEM's Tribol integration (e.g. the contact miniapp).
 
 include(MfemCmakeUtilities)
 # Note: components are enabled based on the find_package() parameters.

config/cmake/modules/FindCAMP.cmake

@@ -0,0 +1,36 @@
+# 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.
+
+# Defines the following variables:
+#   - CAMP_FOUND
+#   - CAMP_LIBRARIES
+#   - CAMP_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+
+mfem_find_package(CAMP CAMP CAMP_DIR
+      "include" "camp/camp.hpp"
+      "lib" "camp"
+      "Paths to headers required by CAMP."
+      "Libraries required by CAMP.")
+
+# RAJA commonly lists "camp" in INTERFACE_LINK_LIBRARIES. If there is no CMake
+# target named "camp", CMake treats it as a bare library name (-lcamp).
+if (CAMP_FOUND AND NOT TARGET camp)
+   list(GET CAMP_LIBRARIES 0 _camp_lib0)
+   add_library(camp UNKNOWN IMPORTED)
+   set_target_properties(camp PROPERTIES
+      IMPORTED_LOCATION "${_camp_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CAMP_INCLUDE_DIRS}")
+   set(CAMP_LIBRARIES "camp" CACHE STRING "CAMP imported target." FORCE)
+   unset(_camp_lib0)
+endif()
+

config/cmake/modules/FindCaliper.cmake

@@ -21,6 +21,21 @@ mfem_find_package(Caliper CALIPER CALIPER_DIR
       "Paths to headers required by Caliper."
       "Libraries required by Caliper.")
 
+# Some downstream CMake packages (notably RAJA) may list "caliper" in their
+# INTERFACE_LINK_LIBRARIES. If there is no CMake target named "caliper", CMake
+# treats it as a bare library name and will pass -lcaliper to the linker.
+# Create a minimal imported target when we only located the library by path.
+if (CALIPER_FOUND AND NOT TARGET caliper)
+   list(GET CALIPER_LIBRARIES 0 _caliper_lib0)
+   add_library(caliper UNKNOWN IMPORTED)
+   set_target_properties(caliper PROPERTIES
+      IMPORTED_LOCATION "${_caliper_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CALIPER_INCLUDE_DIRS}")
+   # Prefer linking via the target.
+   set(CALIPER_LIBRARIES "caliper" CACHE STRING "Caliper imported target." FORCE)
+   unset(_caliper_lib0)
+endif()
+
 # Append adiak path/lib if the user provided ADIAK_DIR
 if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
     find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -878,7 +878,7 @@ function(mfem_export_mk_files)
       MFEM_USE_SUITESPARSE MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS
       MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_MAGMA
       MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC
-      MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
+      MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
       MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA
       MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD
       MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO

config/config.hpp.in

@@ -120,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 // #define MFEM_USE_GNUTLS
 
-// Enable Sidre support.
-// #define MFEM_USE_SIDRE
-
 // Enable the use of SIMD in the high performance templated classes.
 // #define MFEM_USE_SIMD
 

config/config.mk.in

@@ -45,7 +45,6 @@ MFEM_USE_NETCDF        = @MFEM_USE_NETCDF@
 MFEM_USE_PETSC         = @MFEM_USE_PETSC@
 MFEM_USE_SLEPC         = @MFEM_USE_SLEPC@
 MFEM_USE_MPFR          = @MFEM_USE_MPFR@
-MFEM_USE_SIDRE         = @MFEM_USE_SIDRE@
 MFEM_USE_FMS           = @MFEM_USE_FMS@
 MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@

config/defaults.cmake

@@ -48,7 +48,6 @@ option(MFEM_USE_NETCDF "Enable NETCDF usage" OFF)
 option(MFEM_USE_PETSC "Enable PETSc support." OFF)
 option(MFEM_USE_SLEPC "Enable SLEPc support." OFF)
 option(MFEM_USE_MPFR "Enable MPFR usage." OFF)
-option(MFEM_USE_SIDRE "Enable Axom/Sidre usage" OFF)
 option(MFEM_USE_FMS "Enable FMS usage" OFF)
 option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
 option(MFEM_USE_PUMI "Enable PUMI" OFF)
@@ -224,17 +223,8 @@ set(FMS_DIR "${MFEM_DIR}/../fms" CACHE PATH
 set(CONDUIT_DIR "${MFEM_DIR}/../conduit" CACHE PATH
     "Path to the Conduit library.")
 
-set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH "Path to the Axom library.")
-# May need to add "Boost" as requirement.
-if (MFEM_USE_SIDRE)
-    if (MFEM_USE_MPI)
-        set(Axom_REQUIRED_PACKAGES "Conduit/blueprint/blueprint_mpi/relay/relay_mpi" CACHE STRING
-            "Additional packages required by Axom.")
-    elseif()
-        set(Axom_REQUIRED_PACKAGES "Conduit/blueprint/relay" CACHE STRING
-            "Additional packages required by Axom.")
-    endif()
-endif()
+set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH
+    "Path to the Axom library (required by Tribol for the contact mini-app).")
 
 set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
     "Directory where PUMI is installed")
@@ -252,6 +242,7 @@ set(MKL_PARDISO_DIR "" CACHE STRING "MKL installation path.")
 
 set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
 set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
+set(CAMP_DIR "${MFEM_DIR}/../camp" CACHE PATH "Path to CAMP (required by RAJA/Umpire)")
 set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
 set(UMPIRE_DIR "${MFEM_DIR}/../umpire" CACHE PATH "Path to Umpire")
 set(CALIPER_DIR "${MFEM_DIR}/../caliper" CACHE PATH "Path to Caliper")
@@ -272,8 +263,12 @@ set(PARELAG_LIBRARIES "${PARELAG_DIR}/build/src/libParELAG.a" CACHE STRING
     "The ParELAG library.")
 
 set(TRIBOL_DIR "${MFEM_DIR}/../tribol" CACHE PATH "Path to Tribol")
-set(Tribol_REQUIRED_PACKAGES "Axom/core/mint/slam/slic" CACHE STRING
-    "Additional packages required by Tribol")
+# Tribol requires Axom. Many Tribol builds also enable optional TPLs like
+# RAJA/UMPIRE/Caliper, and may pull additional Axom components (e.g. quest,
+# lumberjack) via its exported targets.
+set(Tribol_REQUIRED_PACKAGES
+    "REQUIRED:;Axom/core/primal/mint/slam/slic/quest/lumberjack;OPTIONAL:;Adiak;CAMP;RAJA;UMPIRE;Caliper"
+    CACHE STRING "Additional packages required by Tribol")
 
 set(ENZYME_DIR "${MFEM_DIR}/../enzyme" CACHE PATH "Path to Enzyme")
 

config/defaults.mk

@@ -162,7 +162,6 @@ MFEM_USE_NETCDF        = NO
 MFEM_USE_PETSC         = NO
 MFEM_USE_SLEPC         = NO
 MFEM_USE_MPFR          = NO
-MFEM_USE_SIDRE         = NO
 MFEM_USE_FMS           = NO
 MFEM_USE_CONDUIT       = NO
 MFEM_USE_PUMI          = NO
@@ -249,6 +248,15 @@ endif
 
 # METIS library configuration
 ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
+   # MFEM_USE_METIS_5: when the user supplies METIS_DIR, try to auto-detect
+   # METIS 5 installs that follow the common <prefix>/{include,lib,lib64} layout.
+   ifeq ($(MFEM_USE_METIS_5),NO)
+      ifneq ($(wildcard $(METIS_DIR)/include/metis.h),)
+         ifneq ($(wildcard $(METIS_DIR)/lib/libmetis.* $(METIS_DIR)/lib64/libmetis.*),)
+            MFEM_USE_METIS_5 = YES
+         endif
+      endif
+   endif
    ifeq ($(MFEM_USE_METIS_5),NO)
      METIS_DIR = @MFEM_DIR@/../metis-4.0
      METIS_OPT =
@@ -487,17 +495,6 @@ ifneq (,$(wildcard $(CONDUIT_HDF5_HEADER)))
                   -lhdf5 $(ZLIB_LIB)
 endif
 
-# Sidre and required libraries configuration
-# Be sure to check the HDF5_DIR (set above) is correct
-SIDRE_DIR = @MFEM_DIR@/../axom
-SIDRE_OPT = -I$(SIDRE_DIR)/include -I$(CONDUIT_DIR)/include/conduit\
- -I$(HDF5_DIR)/include
-SIDRE_LIB = \
-   $(XLINKER)-rpath,$(SIDRE_DIR)/lib -L$(SIDRE_DIR)/lib \
-   $(XLINKER)-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
-   $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
-   -laxom -lconduit -lconduit_relay -lconduit_blueprint -lhdf5 $(ZLIB_LIB) -ldl
-
 # PUMI
 # Note that PUMI_DIR is needed -- it is used to check for gmi_sim.h
 PUMI_DIR = @MFEM_DIR@/../pumi-2.1.0
@@ -579,7 +576,13 @@ ifdef CUB_DIR
    RAJA_OPT += -I$(CUB_DIR)
 endif
 
+# CAMP library configuration (required by RAJA/Umpire for most installs)
 CAMP_LIB = -lcamp
+# If the common sibling layout exists, use it as a default (handles versioned
+# directories like camp-<hash>).
+ifneq ($(wildcard $(RAJA_DIR)/../camp*/include/camp/camp.hpp),)
+   CAMP_DIR ?= $(patsubst %/include/camp/camp.hpp,%,$(firstword $(wildcard $(RAJA_DIR)/../camp*/include/camp/camp.hpp)))
+endif
 ifdef CAMP_DIR
    RAJA_OPT += -I$(CAMP_DIR)/include
    CAMP_LIB = $(XLINKER)-rpath,$(CAMP_DIR)/lib -L$(CAMP_DIR)/lib -lcamp
@@ -589,7 +592,12 @@ RAJA_LIB = $(XLINKER)-rpath,$(RAJA_DIR)/lib -L$(RAJA_DIR)/lib -lRAJA $(CAMP_LIB)
 # UMPIRE library configuration
 UMPIRE_DIR = @MFEM_DIR@/../umpire
 UMPIRE_OPT = -I$(UMPIRE_DIR)/include $(if $(CAMP_DIR), -I$(CAMP_DIR)/include)
-UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -L$(UMPIRE_DIR)/lib64 -lumpire $(CAMP_LIB)
+UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -L$(UMPIRE_DIR)/lib64 -lumpire $(CAMP_LIB) -lpthread
+# If the common sibling layout exists, use it as a default (handles versioned
+# directories like fmt-<hash>).
+ifneq ($(wildcard $(UMPIRE_DIR)/../fmt*/include/fmt/format.h),)
+   FMT_DIR ?= $(patsubst %/include/fmt/format.h,%,$(firstword $(wildcard $(UMPIRE_DIR)/../fmt*/include/fmt/format.h)))
+endif
 ifdef FMT_DIR
    UMPIRE_OPT += -I$(FMT_DIR)/include
    UMPIRE_LIB += -L$(FMT_DIR)/lib -L$(FMT_DIR)/lib64 -lfmt
@@ -621,8 +629,63 @@ PARELAG_LIB = -L$(PARELAG_DIR)/build/src -lParELAG
 AXOM_DIR = @MFEM_DIR@/../axom
 TRIBOL_DIR = @MFEM_DIR@/../tribol
 TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
-TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
-   -laxom_slam -laxom_slic -laxom_core
+# Tribol may be built with optional dependencies (e.g. RAJA/UMPIRE/CALIPER).
+# Add those options only when the corresponding headers/libraries exist.
+ifneq ($(wildcard $(RAJA_DIR)/include/RAJA/RAJA.hpp),)
+   TRIBOL_OPT += $(RAJA_OPT)
+endif
+ifneq ($(wildcard $(UMPIRE_DIR)/include/umpire/Umpire.hpp),)
+   TRIBOL_OPT += $(UMPIRE_OPT)
+endif
+ifneq ($(wildcard $(CALIPER_DIR)/include/caliper/cali.h),)
+   TRIBOL_OPT += $(CALIPER_OPT)
+endif
+
+TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -L$(TRIBOL_DIR)/lib64
+ifneq ($(wildcard $(TRIBOL_DIR)/lib/libtribol.* $(TRIBOL_DIR)/lib64/libtribol.*),)
+   TRIBOL_LIB += -ltribol
+endif
+ifneq ($(wildcard $(TRIBOL_DIR)/lib/libtribol_shared.* $(TRIBOL_DIR)/lib64/libtribol_shared.*),)
+   TRIBOL_LIB += -ltribol_shared
+endif
+ifneq ($(wildcard $(TRIBOL_DIR)/lib/libredecomp.* $(TRIBOL_DIR)/lib64/libredecomp.*),)
+   TRIBOL_LIB += -lredecomp
+endif
+
+TRIBOL_LIB += -L$(AXOM_DIR)/lib -L$(AXOM_DIR)/lib64
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_quest.* $(AXOM_DIR)/lib64/libaxom_quest.*),)
+   TRIBOL_LIB += -laxom_quest
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_mint.* $(AXOM_DIR)/lib64/libaxom_mint.*),)
+   TRIBOL_LIB += -laxom_mint
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_slam.* $(AXOM_DIR)/lib64/libaxom_slam.*),)
+   TRIBOL_LIB += -laxom_slam
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_slic.* $(AXOM_DIR)/lib64/libaxom_slic.*),)
+   TRIBOL_LIB += -laxom_slic
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_lumberjack.* $(AXOM_DIR)/lib64/libaxom_lumberjack.*),)
+   TRIBOL_LIB += -laxom_lumberjack
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_core.* $(AXOM_DIR)/lib64/libaxom_core.*),)
+   TRIBOL_LIB += -laxom_core
+endif
+
+# Add common optional Tribol TPLs when their libraries are present.
+ifneq ($(wildcard $(ADIAK_DIR)/lib/libadiak.* $(ADIAK_DIR)/lib64/libadiak.*),)
+   TRIBOL_LIB += $(XLINKER)-rpath,$(ADIAK_DIR)/lib64 $(XLINKER)-rpath,$(ADIAK_DIR)/lib \
+      -L$(ADIAK_DIR)/lib64 -L$(ADIAK_DIR)/lib -ladiak -ldl
+endif
+ifneq ($(wildcard $(UMPIRE_DIR)/lib/libumpire.* $(UMPIRE_DIR)/lib64/libumpire.*),)
+   TRIBOL_LIB += $(UMPIRE_LIB)
+endif
+ifneq ($(wildcard $(RAJA_DIR)/lib/libRAJA.* $(RAJA_DIR)/lib64/libRAJA.*),)
+   TRIBOL_LIB += $(RAJA_LIB)
+endif
+ifneq ($(wildcard $(CALIPER_DIR)/lib/libcaliper.* $(CALIPER_DIR)/lib64/libcaliper.*),)
+   TRIBOL_LIB += $(CALIPER_LIB)
+endif
 
 # Enzyme configuration
 ENZYME_DIR = @MFEM_DIR@/../enzyme

examples/ex36.cpp

@@ -97,7 +97,13 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.ParseCheck();
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
 
    // 2. Read the mesh from the mesh file.
    const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -122,7 +128,7 @@ int main(int argc, char *argv[])
    *nodes /= scale;
 
    // 4. Define the necessary finite element spaces on the mesh.
-   H1Bubble_FECollection H1fec(order, order - 1, dim);
+   H1_FECollection H1fec(order+1, dim);
    FiniteElementSpace H1fes(&mesh, &H1fec);
 
    L2_FECollection L2fec(order-1, dim);

examples/ex36p.cpp

@@ -103,7 +103,19 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.ParseCheck();
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
 
    // 2. Read the mesh from the mesh file.
    const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -131,7 +143,7 @@ int main(int argc, char *argv[])
    mesh.Clear();
 
    // 4. Define the necessary finite element spaces on the mesh.
-   H1Bubble_FECollection H1fec(order, order - 1, dim);
+   H1_FECollection H1fec(order+1, dim);
    ParFiniteElementSpace H1fes(&pmesh, &H1fec);
 
    L2_FECollection L2fec(order-1, dim);

examples/ex41.cpp

@@ -433,16 +433,18 @@ int main(int argc, char *argv[])
    u.ProjectCoefficient(*u0);
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example41", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example41", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/ex41p.cpp

@@ -518,10 +518,12 @@ int main(int argc, char *argv[])
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example41-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example41-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/ex9.cpp

@@ -305,16 +305,18 @@ int main(int argc, char *argv[])
    }
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/ex9p.cpp

@@ -441,16 +441,18 @@ int main(int argc, char *argv[])
    }
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/hiop/ex9.cpp

@@ -354,16 +354,18 @@ int main(int argc, char *argv[])
    }
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9", mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9", mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/hiop/ex9p.cpp

@@ -414,16 +414,18 @@ int main(int argc, char *argv[])
    }
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/petsc/ex9p.cpp

@@ -368,16 +368,18 @@ int main(int argc, char *argv[])
    }
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/sundials/ex9.cpp

@@ -316,16 +316,18 @@ int main(int argc, char *argv[])
    }
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

examples/sundials/ex9p.cpp

@@ -453,16 +453,18 @@ int main(int argc, char *argv[])
    }
 
    // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
    DataCollection *dc = NULL;
    if (visit)
    {
       if (binary)
       {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
       }
       else

fem/CMakeLists.txt

@@ -73,7 +73,6 @@ 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
@@ -134,7 +133,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
@@ -222,7 +221,6 @@ 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
@@ -314,11 +312,6 @@ set(HDRS
   particleset.hpp
   )
 
-if (MFEM_USE_SIDRE)
-  list(APPEND SRCS sidredatacollection.cpp)
-  list(APPEND HDRS sidredatacollection.hpp)
-endif()
-
 if (MFEM_USE_CONDUIT)
   list(APPEND SRCS conduitdatacollection.cpp)
   list(APPEND HDRS conduitdatacollection.hpp)

fem/conduitdatacollection.cpp

@@ -1453,8 +1453,6 @@ ConduitDataCollection::LoadMeshAndFields(int domain_id,
 std::string
 ConduitDataCollection::ElementTypeToShapeName(Element::Type element_type)
 {
-   // Adapted from SidreDataCollection
-
    // Note -- the mapping from Element::Type to string is based on
    //   enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
    //                        TETRAHEDRON, HEXAHEDRON };

fem/conduitdatacollection.hpp

@@ -34,10 +34,10 @@ namespace mfem
     - HDF5 library, https://support.hdfgroup.org/HDF5
 
     @note The ConduitDataCollection only wraps the mfem objects to save them and
-    creates them on load, Conduit does not own any of the data. The
-    SidreDataCollection provides more features, for example the
-    SidreDataCollection allocates and will own the data backing the mfem objects
-    in the data collection.
+    creates them on load, Conduit does not own any of the data.
+    The MFEMSidreDataCollection in the Axom package (https://github.com/LLNL/axom)
+    derives from mfem::DataCollection and provides more features, for example
+    it allocates and will own the data backing the mfem objects in the data collection.
 
     This class also provides public static methods that convert between MFEM
     Meshes and GridFunctions and Conduit Mesh Blueprint descriptions.

fem/fe.hpp

@@ -20,7 +20,6 @@
 #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"

fem/fe/fe_base.hpp

@@ -349,7 +349,7 @@ public:
        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; }
+   virtual int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
 
    /** Returns the dimension of the curl for vector-valued finite elements,
        which is also the width of the DenseMatrix argument in
@@ -360,7 +360,7 @@ public:
        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; }
+   virtual int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
 
    /// Returns the Geometry::Type of the reference element.
    Geometry::Type GetGeomType() const { return geom_type; }
@@ -1017,7 +1017,7 @@ 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; }
+   int GetPhysRangeDim(int space_dim) const override { return space_dim; }
 };
 
 /// @brief Class for computing 1D special polynomials and their associated basis

fem/fe/fe_h1_bubble.cpp

@@ -1,973 +0,0 @@
-// 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);
-}
-
-}

fem/fe/fe_h1_bubble.hpp

@@ -1,109 +0,0 @@
-// 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

fem/fe/fe_nd.hpp

@@ -663,8 +663,8 @@ 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; }
+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { return 1; }
 
    void CalcVShape(const IntegrationPoint &ip,
                    DenseMatrix &shape) const override;
@@ -708,8 +708,8 @@ private:
                            DenseMatrix &I) const;
 
 public:
-   int GetPhysRangeDim(int space_dim) const { return 3; }
-   int GetPhysCurlDim(int space_dim) const { return 3; }
+   int GetPhysRangeDim(int space_dim) const override { return 3; }
+   int GetPhysCurlDim(int space_dim) const override { return 3; }
 
    using FiniteElement::CalcVShape;
    using FiniteElement::CalcPhysCurlShape;

fem/fe/fe_rt.hpp

@@ -510,8 +510,8 @@ 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; }
+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { return 0; }
 
    void CalcVShape(const IntegrationPoint &ip,
                    DenseMatrix &shape) const override;
@@ -550,8 +550,8 @@ private:
                            DenseMatrix &I) const;
 
 public:
-   int GetPhysRangeDim(int space_dim) const { return 3; }
-   int GetPhysCurlDim(int space_dim) const { return 0; }
+   int GetPhysRangeDim(int space_dim) const override { return 3; }
+   int GetPhysCurlDim(int space_dim) const override { return 0; }
 
    using FiniteElement::CalcVShape;
 

fem/fe_coll.cpp

@@ -243,21 +243,11 @@ 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));
@@ -2132,173 +2122,6 @@ 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)

fem/fe_coll.hpp

@@ -111,8 +111,6 @@ 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) |
@@ -319,59 +317,6 @@ 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

fem/fem.hpp

@@ -62,10 +62,6 @@
 #include "pnonlinearform.hpp"
 #endif
 
-#ifdef MFEM_USE_SIDRE
-#include "sidredatacollection.hpp"
-#endif
-
 #ifdef MFEM_USE_CONDUIT
 #include "conduitdatacollection.hpp"
 #endif

fem/fespace.cpp

@@ -3877,12 +3877,9 @@ 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). 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)
+      // consistency check: fec->GetOrder() and FE->GetOrder() should return
+      // the same value (for standard, constant-order spaces)
+      if (!IsVariableOrder() && FE->GetDim() > 0)
       {
          MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
                      "internal error: " <<
@@ -3937,6 +3934,16 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const
    return BE;
 }
 
+const FiniteElement *FiniteElementSpace::GetTypicalBE() const
+{
+   if (mesh->GetNBE() > 0) { return GetBE(0); }
+
+   Geometry::Type geom = mesh->GetTypicalFaceGeometry();
+   const FiniteElement *be = fec->FiniteElementForGeometry(geom);
+   MFEM_VERIFY(be != nullptr, "Could not determine a typical BE!");
+   return be;
+}
+
 const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
 {
    MFEM_VERIFY(!IsVariableOrder(), "not implemented");
@@ -3967,6 +3974,11 @@ const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
    return fe;
 }
 
+const FiniteElement *FiniteElementSpace::GetTypicalFaceElement() const
+{
+   return fec->FiniteElementForGeometry(mesh->GetTypicalFaceGeometry());
+}
+
 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
                                                         int variant) const
 {

fem/fespace.hpp

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

fem/gridfunc.cpp

@@ -345,27 +345,6 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi,
    }
 }
 
-int GridFunction::VectorDim() const
-{
-   const FiniteElement *fe = fes->GetTypicalFE();
-   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
-   {
-      return fes->GetVDim();
-   }
-   return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
-                                  fe->GetRangeDim());
-}
-
-int GridFunction::CurlDim() const
-{
-   const FiniteElement *fe = fes->GetTypicalFE();
-   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
-   {
-      return 2 * fes->GetMesh()->SpaceDimension() - 3;
-   }
-   return fes->GetVDim()*fe->GetCurlDim();
-}
-
 void GridFunction::GetTrueDofs(Vector &tv) const
 {
    const SparseMatrix *R = fes->GetRestrictionMatrix();
@@ -2050,6 +2029,18 @@ void GridFunction::AccumulateAndCountBdrValues(
    Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
    Array<int> &values_counter)
 {
+   if (vcoeff)
+   {
+      MFEM_VERIFY(fes->GetVDim() == vcoeff->GetVDim(),
+                  "vcoeff vdim != fes VDim");
+      MFEM_VERIFY(fes->GetTypicalBE()->GetMapType() == FiniteElement::VALUE &&
+                  fes->GetTypicalBE()->GetRangeType() ==
+                  FiniteElement::SCALAR,
+                  "Can only call ProjectBdrCoefficient on scalar value-type "
+                  "boundary elements. "
+                  "Did you intended to call ProjectBdrCoefficientNormal or "
+                  "ProjectBdrCoefficientTangent for vector finite elements?");
+   }
    Array<int> vdofs;
    Vector vc;
 
@@ -2202,6 +2193,9 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
    VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
    Array<int> &values_counter)
 {
+   MFEM_VERIFY(fes->GetTypicalBE()->GetPhysRangeDim(
+                  fes->GetMesh()->SpaceDimension()) == vcoeff.GetVDim(),
+               "vcoeff vdim != PhysRangeDim");
    const FiniteElement *fe;
    ElementTransformation *T;
    Array<int> dofs;
@@ -2355,6 +2349,9 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
 
 void GridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
+   MFEM_VERIFY(
+      VectorDim() == 1,
+      "Cannot project scalar Coefficient onto vector GridFunction");
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
    DofTransformation doftrans;
    Array<int> vdofs;
@@ -2630,6 +2627,7 @@ void GridFunction::ProjectCoefficient(
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff,
                                       ProjectType type)
 {
+   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    Array<int> vdofs;
    Vector vals;
    DofTransformation doftrans;
@@ -2945,6 +2943,7 @@ void GridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 void GridFunction::ProjectCoefficient(
    VectorCoefficient &vcoeff, Array<int> &dofs)
 {
+   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int el = -1;
    ElementTransformation *T = NULL;
    const FiniteElement *fe = NULL;
@@ -2974,6 +2973,7 @@ void GridFunction::ProjectCoefficient(
 
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff, int attribute)
 {
+   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int i;
    Array<int> vdofs;
    Vector vals;
@@ -3033,6 +3033,7 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
                                           Array<int> &dof_attr)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> vdofs;
    Vector vals;
 
@@ -3064,6 +3065,7 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
 
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> dof_attr;
    ProjectDiscCoefficient(coeff, dof_attr);
 }
@@ -3073,6 +3075,10 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
+   MFEM_VERIFY(
+      VectorDim() == 1,
+      "Cannot project a scalar coefficient onto a vector GridFunction");
+
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -3082,6 +3088,7 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
                                           AvgType type)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -3139,12 +3146,16 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
 void GridFunction::ProjectBdrCoefficientNormal(
    Coefficient *coeff, VectorCoefficient *vcoeff, const Array<int> &bdr_attr)
 {
-   if (fes->GetNBE() > 0)
+   MFEM_VERIFY(fes->GetVDim() == 1, "fespace VDim != 1");
+   MFEM_VERIFY(fes->GetTypicalBE()->GetRangeType() == FiniteElement::SCALAR &&
+               fes->GetTypicalBE()->GetMapType() == FiniteElement::INTEGRAL,
+               "Not an RT FE space!");
+   if (vcoeff)
    {
-      // 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!");
+      MFEM_VERIFY(vcoeff->GetVDim() == fes->GetMesh()->SpaceDimension(),
+                  "vcoeff vdim (" << vcoeff->GetVDim()
+                  << ") != SpaceDimension ("
+                  << fes->GetMesh()->SpaceDimension() << ")");
    }
 
    // implementation for the case when the face dofs are scaled point
@@ -5757,4 +5768,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
    return std::make_pair(global_max_lower, global_max_upper);
 }
 
-}
+}

fem/gridfunc.hpp

@@ -150,11 +150,13 @@ public:
 
    FiniteElementCollection *OwnFEC() { return fec_owned; }
 
-   /// Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying #fes
-   int VectorDim() const;
+   /** @brief Shortcut for calling FiniteElementSpace::GetVectorDim() on the
+       underlying #fes */
+   int VectorDim() const { return fes->GetVectorDim(); }
 
-   /// Shortcut for calling FiniteElementSpace::GetCurlDim() on the underlying #fes
-   int CurlDim() const;
+   /** @brief Shortcut for calling FiniteElementSpace::GetCurlDim() on the
+       underlying #fes */
+   int CurlDim() const { return fes->GetCurlDim(); }
 
    /// Read only access to the (optional) internal true-dof Vector.
    const Vector &GetTrueVector() const
@@ -1971,6 +1973,7 @@ public:
 
    void Eval(Vector &v, ElementTransformation &T,
              const IntegrationPoint &ip) override;
+   using VectorCoefficient::Eval;
 
    virtual ~VectorExtrudeCoefficient() { }
 };

fem/pgridfunc.cpp

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

fem/pgridfunc.hpp

@@ -280,8 +280,7 @@ public:
    using GridFunction::ProjectBdrCoefficient;
 
    void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
-                              const Array<int> &attr) override
-   { ProjectBdrCoefficient(NULL, &vcoeff, attr); }
+                              const Array<int> &attr) override;
 
    void ProjectBdrCoefficient(Coefficient *coeff[],
                               const Array<int> &attr) override

fem/sidredatacollection.hpp

@@ -1,539 +0,0 @@
-// 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_SIDREDATACOLLECTION
-#define MFEM_SIDREDATACOLLECTION
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_SIDRE
-
-#include "datacollection.hpp"
-
-// Ignore warnings from the axom/sidre header (GCC + Clang versions)
-#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
-# pragma GCC diagnostic push
-# if defined(__clang__)
-#  pragma GCC diagnostic ignored "-Wextra-semi"
-# else // real GCC?
-#  pragma GCC diagnostic ignored "-Wpedantic"
-# endif
-#endif
-#include <axom/sidre.hpp>
-#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
-# pragma GCC diagnostic pop
-#endif
-
-namespace mfem
-{
-
-/** @brief Data collection with Sidre routines following the Conduit mesh
-    blueprint specification. */
-/** SidreDataCollection provides an HDF5-based file format for visualization or
-    restart capability.  This functionality is aimed primarily at customers of
-    LLNL's axom project that run problems at extreme scales.
-
-    For more information, see:
-    - Sidre component of LLNL's axom project (to be open-sourced), http://goo.gl/cZyJdn
-    - LLNL conduit/blueprint library, https://github.com/LLNL/conduit
-    - HDF5 library, https://support.hdfgroup.org/HDF5
-
-    The layout created in the Sidre DataStore is: (`"──"` denote groups,
-    `"─•"` denote views, `"─>"` denote links, i.e. shallow-copy view)
-
-        <root>
-         ├── <collection-name>_global          (global group)
-         │    └── blueprint_index
-         │         └── <collection-name>       (bp_index group)
-         │              ├── state
-         │              │    ├─• cycle
-         │              │    ├─• time
-         │              │    └─• number_of_domains = <mesh-mpi-comm-size>
-         │              ├── coordsets
-         │              │    └── coords
-         │              │         ├─• path = "<bp-path>/coordsets/coords"
-         │              │         ├─• type ─> <bp-grp>/coordsets/coords/type = "explicit"
-         │              │         └─• coord_system = "x"|"xy"|"xyz"
-         │              ├── topologies
-         │              │    ├── mesh
-         │              │    │    ├─• path = "<bp-path>/topologies/mesh"
-         │              │    │    ├─• type              ─> <bp-grp>/topologies/mesh/type = "unstructured"
-         │              │    │    ├─• coordset          ─> <bp-grp>/topologies/mesh/coordset = "coords"
-         │              │    │    ├─• grid_function     ─> <bp-grp>/topologies/mesh/grid_function = "<nodes-field-name>"
-         │              │    │    └─• boundary_topology ─> <bp-grp>/topologies/mesh/boundary_topology = "boundary"
-         │              │    └── boundary
-         │              │         ├─• path = "<bp-path>/topologies/mesh"
-         │              │         ├─• type     ─> <bp-grp>/topologies/boundary/type = "unstructured"
-         │              │         └─• coordset ─> <bp-grp>/topologies/boundary/coordset = "coords"
-         │              └── fields
-         │                   ├── mesh_material_attribute
-         │                   │    ├─• path = "<bp-path>/fields/mesh_material_attribute"
-         │                   │    ├─• association ─> <bp-grp>/fields/mesh_material_attribute/association = "element"
-         │                   │    ├─• topology    ─> <bp-grp>/fields/mesh_material_attribute/topology = "mesh"
-         │                   │    └─• number_of_components = 1
-         │                   ├── boundary_material_attribute
-         │                   │    ├─• path = "<bp-path>/fields/boundary_material_attribute"
-         │                   │    ├─• association ─> <bp-grp>/fields/boundary_material_attribute/association = "element"
-         │                   │    ├─• topology    ─> <bp-grp>/fields/boundary_material_attribute/topology = "boundary"
-         │                   │    └─• number_of_components = 1
-         │                   ├── grid-function-1
-         │                   │    ├─• path = "<bp-path>/fields/grid-function-1"
-         │                   │    ├─• basis    ─> <bp-grp>/fields/grid-function-1/basis = "<fe-coll-name>"
-         │                   │    ├─• topology ─> <bp-grp>/fields/grid-function-1/topology = "mesh"
-         │                   │    └─• number_of_components = gf1->VectorDim()
-         │                   ├── grid-function-2
-         │                   │    ├─• path = "<bp-path>/fields/grid-function-2"
-         │                   │    ├─• basis    ─> <bp-grp>/fields/grid-function-2/basis = "<fe-coll-name>"
-         │                   │    ├─• topology ─> <bp-grp>/fields/grid-function-2/topology = "mesh"
-         │                   │    └─• number_of_components = gf2->VectorDim()
-         │                   ├── ...
-         │                  ...
-         └── <collection-name>                 (domain group)
-              ├── blueprint                    (blueprint group)
-              │    ├── state
-              │    │    ├─• cycle
-              │    │    ├─• time
-              │    │    ├─• domain = <mesh-mpi-rank>
-              │    │    └─• time_step
-              │    ├── coordsets
-              │    │    └── coords
-              │    │         ├─• type = "explicit"
-              │    │         └── values
-              │    │              ├─• x = view in <vertex-coords-buffer>/<ext-double-data>
-              │    │              ├─• y = view in <vertex-coords-buffer>/<ext-double-data>
-              │    │              └─• z = view in <vertex-coords-buffer>/<ext-double-data>
-              │    ├── topologies
-              │    │    ├── mesh
-              │    │    │    ├─• type = "unstructured"
-              │    │    │    ├── elements
-              │    │    │    │    ├─• shape = "points"|"lines"|...
-              │    │    │    │    └─• connectivity = <vert-idx-array>
-              │    │    │    ├─• coordset = "coords"
-              │    │    │    ├─• grid_function = "<nodes-field-name>"
-              │    │    │    └─• boundary_topology = "boundary"
-              │    │    └── boundary
-              │    │         ├─• type = "unstructured"
-              │    │         ├── elements
-              │    │         │    ├─• shape = "points"|"lines"|...
-              │    │         │    └─• connectivity = <vert-idx-array>
-              │    │         └─• coordset = "coords"
-              │    └── fields
-              │         ├── mesh_material_attribute
-              │         │    ├─• association = "element"
-              │         │    ├─• topology = "mesh"
-              │         │    └─• values = <attr-array>
-              │         ├── boundary_material_attribute
-              │         │    ├─• association = "element"
-              │         │    ├─• topology = "boundary"
-              │         │    └─• values = <attr-array>
-              │         ├── grid-function-1   (name can include path)
-              │         │    ├─• basis = "<fe-coll-name>"
-              │         │    ├─• topology = "mesh"
-              │         │    └─• values = <ext-double-array>/<named-buffer> (vdim == 1)
-              │         ├── grid-function-2   (name can include path)
-              │         │    ├─• basis = "<fe-coll-name>"
-              │         │    ├─• topology = "mesh"
-              │         │    └── values   (vdim > 1)
-              │         │         ├─• x0 = view into <ext-double-array>/<named-buffer>
-              │         │         ├─• x1 = view into <ext-double-array>/<named-buffer>
-              │         │         └─• x2 = view into <ext-double-array>/<named-buffer>
-              │         ├── ...
-              │        ...
-              └── named_buffers                (named_buffers group)
-                   ├─• vertex_coords = <double-array>
-                   ├─• grid-function-1 = <double-array>
-                   ├─• grid-function-2 = <double-array>
-                  ...
-
-    @note blueprint_index is used both in serial and in parallel. In parallel,
-    only rank 0 will add entries to the blueprint index.
-
-    @note QuadratureFunction%s (q-fields) are not supported.
-
-    @note SidreDataCollection does not manage the FiniteElementSpace%s and
-    FiniteElementCollection%s associated with registered GridFunction%s.
-    Therefore, field registration is left to the user of SidreDataCollection and
-    there are no methods that automatically register GridFunction%s using just
-    the content of the Sidre DataStore. Such capabilities can be implemented in
-    a derived class, adding any desired object management routines.
-
-    @warning This class is still _experimental_, meaning that in future
-    releases, it may not be backward compatible, and the output files generated
-    by the current version may become unreadable.
-*/
-class SidreDataCollection : public DataCollection
-{
-public:
-   typedef NamedFieldsMap< Array<int> > AttributeFieldMap;
-   AttributeFieldMap attr_map;
-
-public:
-
-   /// Constructor that allocates and initializes a Sidre DataStore.
-   /**
-       @param[in] collection_name  Name of the collection used as a file name
-                                   when saving
-       @param[in] the_mesh         Mesh shared by all grid functions in the
-                                   collection (can be NULL)
-       @param[in] owns_mesh_data   Does the SidreDC own the mesh vertices?
-
-       With this constructor, the SidreDataCollection owns the allocated Sidre
-       DataStore.
-    */
-   SidreDataCollection(const std::string& collection_name,
-                       Mesh *the_mesh = NULL,
-                       bool owns_mesh_data = false);
-
-   /// Constructor that links to an external Sidre DataStore.
-   /** Specifically, the global and domain groups can be at arbitrary paths.
-
-       @param[in] collection_name  Name of the collection used as a file name
-                                   when saving
-       @param[in] bp_index_grp     Pointer to the blueprint index group in the
-                                   datastore, see the above schematic
-       @param[in] domain_grp       Pointer to the domain group in the datastore,
-                                   see the above schematic
-       @param[in] owns_mesh_data   Does the SidreDC own the mesh vertices?
-
-       With this constructor, the SidreDataCollection does not own the Sidre
-       DataStore.
-       @note No mesh or fields are read from the given Groups. The mesh has
-       to be set with SetMesh() and fields registered with RegisterField().
-    */
-   SidreDataCollection(const std::string& collection_name,
-                       axom::sidre::Group * bp_index_grp,
-                       axom::sidre::Group * domain_grp,
-                       bool owns_mesh_data = false);
-
-#ifdef MFEM_USE_MPI
-   /// Associate an MPI communicator with the collection.
-   /** If no mesh was associated with the collection, this method should be
-       called before using any of the Load() methods to read parallel data. */
-   void SetComm(MPI_Comm comm);
-#endif
-
-   /// Register a GridFunction in the Sidre DataStore.
-   /** This method is a shortcut for the call
-       `RegisterField(field_name, gf, field_name, 0)`.
-    */
-   virtual void RegisterField(const std::string &field_name, GridFunction *gf)
-   {
-      RegisterField(field_name, gf, field_name, 0);
-   }
-
-   /// Register a GridFunction in the Sidre DataStore.
-   /** The registration procedure is as follows:
-       - if (@a gf's data is NULL), allocate named buffer with the name
-         @a buffer_name with size _offset + gf->FESpace()->GetVSize()_ and use
-         its data (plus the given @a offset) to set @a gf's data;
-       - else, if (DataStore has a named buffer @a buffer_name), replace @a gf's
-         data array with that named buffer plus the given @a offset;
-       - else, use @a gf's data as external data associated with @a field_name
-         in the DataStore;
-       - register @a field_name in #field_map.
-
-       Both the @a field_name and @a buffer_name can contain a path prefix.
-       @note If @a field_name or @a buffer_name is empty, the method does
-       nothing.
-       @note If the GridFunction pointer @a gf or it's FiniteElementSpace
-       pointer are NULL, the method does nothing.
-    */
-   void RegisterField(const std::string &field_name, GridFunction *gf,
-                      const std::string &buffer_name,
-                      axom::sidre::IndexType offset);
-
-   /// Registers an attribute field in the Sidre DataStore
-   /** The registration process is similar to that of RegisterField()
-       The attribute field is associated with the elements of the mesh
-       when @a is_bdry is false, and with the boundary elements, when
-       @a is_bdry is true.
-       @sa RegisterField()  */
-   void RegisterAttributeField(const std::string& name, bool is_bdry);
-   void DeregisterAttributeField(const std::string& name);
-
-   /** Returns a pointer to the attribute field associated with
-       @a field_name, or NULL when there is no associated field */
-   Array<int>* GetAttributeField(const std::string& field_name) const
-   { return attr_map.Get(field_name); }
-
-   /** Checks if there is an attribute field associated with @a field_name */
-   bool HasAttributeField(const std::string& field_name) const
-   { return attr_map.Has(field_name); }
-
-   /** Checks if any rank in the mesh has boundary elements */
-   bool HasBoundaryMesh() const;
-
-   /// Set the name of the mesh nodes field.
-   /** This name will be used by SetMesh() to register the mesh nodes, if not
-       already registered. Also, this method should be called if the mesh nodes
-       GridFunction was or will be registered directly by the user. The default
-       value for the name is "mesh_nodes". */
-   void SetMeshNodesName(const std::string &nodes_name)
-   {
-      if (!nodes_name.empty()) { m_meshNodesGFName = nodes_name; }
-   }
-
-   /// De-register @a field_name from the SidreDataCollection.
-   /** The field is removed from the #field_map and the DataStore, including
-       deleting it from the named_buffers group, if allocated. */
-   virtual void DeregisterField(const std::string& field_name);
-
-   /// Delete all owned data.
-   virtual ~SidreDataCollection();
-
-   /// Set/change the mesh associated with the collection
-   /** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
-       to register the mesh nodes GridFunction, if the mesh uses nodes. */
-   virtual void SetMesh(Mesh *new_mesh);
-
-#ifdef MFEM_USE_MPI
-   /// Set/change the mesh associated with the collection
-   /** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
-       to register the mesh nodes GridFunction, if the mesh uses nodes. */
-   virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh);
-#endif
-
-   /// Reset the domain and global datastore group pointers.
-   /** These are set in the constructor, but if a host code changes the
-       datastore contents ( such as wiping out the datastore and loading in new
-       contents from a file, i.e. a restart ) these pointers will need to be
-       reset to valid groups in the datastore.
-       @sa Load(const std::string &path, const std::string &protocol).
-    */
-   void SetGroupPointers(axom::sidre::Group * global_grp,
-                         axom::sidre::Group * domain_grp);
-
-   axom::sidre::Group * GetBPGroup() { return m_bp_grp; }
-   axom::sidre::Group * GetBPIndexGroup() { return m_bp_index_grp; }
-
-   /// Prepare the DataStore for writing
-   virtual void PrepareToSave();
-
-   /// Save the collection to file.
-   /** This method calls `Save(collection_name, "sidre_hdf5")`. */
-   virtual void Save();
-
-   /// Save the collection to @a filename.
-   /** The collection path prefix is prepended to the @a filename and the
-       current cycle is appended, if cycle >= 0. */
-   void Save(const std::string& filename, const std::string& protocol);
-
-   /// Load the Sidre DataStore from file.
-   /** No mesh or fields are read from the loaded DataStore.
-
-       If the data collection created the datastore, it knows the layout of
-       where the domain and global groups are, and can restore them after the
-       Load().
-
-       If, however, the data collection does not own the datastore (e.g. it did
-       not create the datastore), the host code must reset these pointers after
-       the load operation, using SetGroupPointers(), and also reset the state
-       variables, using UpdateStateFromDS().
-    */
-   void Load(const std::string& path, const std::string& protocol);
-
-   /// Load SidreDataCollection from file.
-   /** The used file path is based on the current prefix path, collection name,
-       and the given @a cycle_. The protocol is "sidre_hdf5".
-       @sa Load(const std::string &path, const std::string &protocol).
-    */
-   virtual void Load(int cycle_ = 0)
-   {
-      SetCycle(cycle_);
-      Load(get_file_path(name), "sidre_hdf5");
-   }
-
-   /// Load external data after registering externally owned fields.
-   void LoadExternalData(const std::string& path);
-
-   /** @brief Updates the DataCollection's cycle, time, and time-step variables
-       with the values from the data store. */
-   void UpdateStateFromDS();
-
-   /** @brief Updates the data store's cycle, time, and time-step variables with
-       the values from the SidreDataCollection. */
-   void UpdateStateToDS();
-
-   /** @name Methods for named buffer access and manipulation. */
-   ///@{
-
-   /** @brief Get a pointer to the sidre::View holding the named buffer for
-       @a buffer_name. */
-   /** If such named buffer is not allocated, the method returns NULL.
-       @note To access the underlying pointer, use View::getData().
-       @note To query the size of the buffer, use View::getNumElements().
-    */
-   axom::sidre::View *
-   GetNamedBuffer(const std::string& buffer_name) const
-   {
-      return named_buffers_grp()->hasView(buffer_name)
-             ? named_buffers_grp()->getView(buffer_name)
-             : NULL;
-   }
-
-   /// Return newly allocated or existing named buffer for @a buffer_name.
-   /** The buffer is stored in the named_buffers group. If the currently
-       allocated buffer size is smaller than @a sz, then the buffer is
-       reallocated with size @a sz, destroying its contents.
-       @note To access the underlying pointer, use View::getData().
-    */
-   axom::sidre::View *
-   AllocNamedBuffer(const std::string& buffer_name,
-                    axom::sidre::IndexType sz,
-                    axom::sidre::TypeID type =
-                       axom::sidre::DOUBLE_ID);
-
-   /// Deallocate the named buffer @a buffer_name.
-   void FreeNamedBuffer(const std::string& buffer_name)
-   { named_buffers_grp()->destroyViewAndData(buffer_name); }
-
-   ///@}
-
-private:
-   // Used if the Sidre data collection is providing the datastore itself.
-   const bool m_owns_datastore;
-
-   // TODO - Need to evaluate if this bool member can be combined with own_data
-   // in parent data collection class. m_owns_mesh_data indicates whether the
-   // Sidre dc owns the mesh element data and node positions gf. The DC base
-   // class own_data indicates if the dc owns the mesh object pointer itself and
-   // GF objects. Can we use one flag and just have DC own all objects vs none?
-   const bool m_owns_mesh_data;
-
-   // Name to be used for registering the mesh nodes in the SidreDataCollection.
-   // This name is used by SetMesh() and can be overwritten by the method
-   // SetMeshNodesName().
-   // Default value: "mesh_nodes".
-   std::string m_meshNodesGFName;
-
-   // If the data collection owns the datastore, it will store a pointer to it.
-   // Otherwise, this pointer is NULL.
-   axom::sidre::DataStore * m_datastore_ptr;
-
-protected:
-   axom::sidre::Group *named_buffers_grp() const;
-
-   axom::sidre::View *
-   alloc_view(axom::sidre::Group *grp,
-              const std::string &view_name);
-
-   axom::sidre::View *
-   alloc_view(axom::sidre::Group *grp,
-              const std::string &view_name,
-              const axom::sidre::DataType &dtype);
-
-   axom::sidre::Group *
-   alloc_group(axom::sidre::Group *grp,
-               const std::string &group_name);
-
-   // return the filename based on prefix_path, collection name and cycle.
-   std::string get_file_path(const std::string &filename) const;
-
-private:
-   // If the data collection does not own the datastore, it will need pointers
-   // to the blueprint and blueprint index group to use.
-   axom::sidre::Group * m_bp_grp;
-   axom::sidre::Group * m_bp_index_grp;
-
-   // This is stored for convenience.
-   axom::sidre::Group * m_named_bufs_grp;
-
-   // Private helper functions
-
-   void RegisterFieldInBPIndex(const std::string& field_name,
-                               GridFunction *gf);
-   void DeregisterFieldInBPIndex(const std::string & field_name);
-
-   void RegisterAttributeFieldInBPIndex(const std::string& attr_name);
-   void DeregisterAttributeFieldInBPIndex(const std::string& attr_name);
-
-   /** @brief Return a string with the conduit blueprint name for the given
-       Element::Type. */
-   std::string getElementName( Element::Type elementEnum );
-
-   /**
-    * \brief A private helper function to set up the views associated with the
-       data of a scalar valued grid function in the blueprint style.
-    * \pre gf is not null
-    * \note This function is expected to be called by RegisterField()
-    * \note Handles cases where hierarchy is already set up,
-    *      where the data was allocated by this data collection
-    *      and where the grid function data is external to Sidre
-    */
-   void addScalarBasedGridFunction(const std::string& field_name,
-                                   GridFunction* gf,
-                                   const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
-
-   /**
-    * \brief A private helper function to set up the views associated with the
-       data of a vector valued grid function in the blueprint style.
-    * \pre gf is not null
-    * \note This function is expected to be called by RegisterField()
-    * \note Handles cases where hierarchy is already set up,
-    *      where the data was allocated by this data collection
-    *      and where the grid function data is external to Sidre
-    */
-   void addVectorBasedGridFunction(const std::string& field_name,
-                                   GridFunction* gf,
-                                   const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
-
-   /** @brief A private helper function to set up the Views associated with
-       attribute field named @a field_name */
-   void addIntegerAttributeField(const std::string& field_name, bool is_bdry);
-
-   /// Sets up the four main mesh blueprint groups.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    */
-   void createMeshBlueprintStubs(bool hasBP);
-
-   /// Sets up the mesh blueprint 'state' group.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    */
-   void createMeshBlueprintState(bool hasBP);
-
-   /// Sets up the mesh blueprint 'coordsets' group.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    */
-   void createMeshBlueprintCoordset(bool hasBP);
-
-   /// Sets up the mesh blueprint 'topologies' group.
-   /**
-    * This method is called from SetMesh().
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    * \param mesh_name The name of the topology.
-    * \note Valid values for @a mesh_name are "mesh" and "boundary" and the
-            former has to be created with this method before the latter.
-    */
-   void createMeshBlueprintTopologies(bool hasBP, const std::string& mesh_name);
-
-#ifdef MFEM_USE_MPI
-   /// Sets up the mesh blueprint 'adjacencies' group.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    * \note Only valid when using parallel meshes
-    */
-   void createMeshBlueprintAdjacencies(bool hasBP);
-#endif
-
-   /// Verifies that the contents of the mesh blueprint data is valid.
-   void verifyMeshBlueprint();
-};
-
-} // end namespace mfem
-
-#endif
-
-#endif

general/version.cpp

@@ -160,9 +160,6 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_RAJA
       "MFEM_USE_RAJA\n"
 #endif
-#ifdef MFEM_USE_SIDRE
-      "MFEM_USE_SIDRE\n"
-#endif
 #ifdef MFEM_USE_SIMD
       "MFEM_USE_SIMD\n"
 #endif

linalg/petsc.cpp

@@ -4156,20 +4156,31 @@ void PetscNonlinearSolver::SetUpdate(void (*update)(Operator *,int,
 void PetscNonlinearSolver::Mult(const Vector &b, Vector &x) const
 {
    SNES snes = (SNES)obj;
+   MPI_Comm comm = PetscObjectComm(obj);
 
-   bool b_nonempty = b.Size();
-   if (!B) { B = new PetscParVector(PetscObjectComm(obj), *this, true); }
-   if (!X) { X = new PetscParVector(PetscObjectComm(obj), *this, false, false); }
+   // Reduction needed: some processes may have null local size while others don't,
+   // and VecPlaceArray (used by PlaceMemory) is a logically collective operation.
+   PetscBool b_nonempty = b.Size() ? PETSC_TRUE : PETSC_FALSE;
+#if PETSC_VERSION_LT(3,24,0)
+   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPIU_BOOL,MPI_LOR,comm);
+#else
+   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPI_C_BOOL,MPI_LOR,comm);
+#endif
+   CCHKERRQ(comm,mpiierr);
+
+   // Always create B with allocate=false so that PlaceMemory can be called on
+   // it regardless of whether b was empty on a previous call.
+   if (!B) { B = new PetscParVector(comm, *this, true, false); }
+   if (!X) { X = new PetscParVector(comm, *this, false, false); }
    X->PlaceMemory(x.GetMemory(),iterative_mode);
    if (b_nonempty) { B->PlaceMemory(b.GetMemory()); }
-   else { *B = 0.0; }
 
    Customize();
 
    if (!iterative_mode) { *X = 0.; }
 
-   // Solve the system.
-   ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
+   // Solve the system. Pass nullptr for b when empty (PETSc treats it as zero RHS).
+   ierr = SNESSolve(snes, b_nonempty ? B->x : nullptr, X->x); PCHKERRQ(snes, ierr);
    X->ResetMemory();
    if (b_nonempty) { B->ResetMemory(); }
 }

makefile

@@ -299,7 +299,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT LAPACK SUNDIALS\
  SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
  TRIBOL BENCHMARK MOONOLITH ALGOIM
@@ -365,7 +365,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS\
  MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_HDF5 MFEM_USE_NETCDF MFEM_USE_PETSC\
- MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
+ MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
@@ -746,7 +746,6 @@ status info:
 	$(info MFEM_USE_PETSC         = $(MFEM_USE_PETSC))
 	$(info MFEM_USE_SLEPC         = $(MFEM_USE_SLEPC))
 	$(info MFEM_USE_MPFR          = $(MFEM_USE_MPFR))
-	$(info MFEM_USE_SIDRE         = $(MFEM_USE_SIDRE))
 	$(info MFEM_USE_FMS           = $(MFEM_USE_FMS))
 	$(info MFEM_USE_CONDUIT       = $(MFEM_USE_CONDUIT))
 	$(info MFEM_USE_PUMI          = $(MFEM_USE_PUMI))

mesh/pmesh.cpp

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

miniapps/contact/README.md

@@ -82,6 +82,43 @@ Although Tribol can be built automatically via **uberenv** and **Spack**,
 for this miniapp it is simpler to build **Axom** and **MFEM** manually and 
 point Tribol to them. The steps are as follows:
 
+### Using pre-built Tribol/Axom installs
+
+If you already have compatible installs of Tribol and Axom, point MFEM to the install prefixes.
+
+- Hypre install prefix: `<path/to/hypre>`
+- METIS install prefix: `<path/to/metis>`
+- Axom install prefix: `<path/to/axom>`
+- Tribol install prefix: `<path/to/tribol>`
+
+**MFEM make build (configure):**
+```bash
+make config MFEM_USE_MPI=YES MFEM_USE_METIS=YES MFEM_USE_TRIBOL=YES \
+  HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
+  AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
+  UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
+```
+
+**MFEM CMake build (configure):**
+```bash
+cmake -S . -B <mfem-build-dir> -DMFEM_USE_MPI=YES -DMFEM_USE_METIS=YES -DMFEM_USE_TRIBOL=YES \
+  HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
+  AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
+  UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
+```
+
+Note: RAJA/UMPIRE/CALIPER are optional for MFEM itself, but many Tribol builds
+enable them. If your Tribol install does not depend on them, you can omit the
+corresponding `*_DIR` entries above.
+
+Note: `FMT_DIR` only needs to be added for the make-based build (and only when
+the Umpire install uses `fmt`). If `FMT_DIR` is not set and a sibling `fmt-*`
+directory exists next to your `UMPIRE_DIR`, MFEM's make configuration will try
+to pick it up automatically.
+
+Note: when using pre-built Tribol/Axom, you typically need to use a compatible
+compiler/MPI wrapper (same C++ standard library ABI).
+
 ### Manual Build Steps
 
 1. Pull axom and tribol (starting from the mfem folder):
@@ -99,7 +136,7 @@ point Tribol to them. The steps are as follows:
    TRIBOL_DIR = @MFEM_DIR@/../tribol-repo/tribol
    TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
    TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib \
-   -laxom_mint -laxom_slam -laxom_slic -laxom_core
+   -laxom_quest -laxom_mint -laxom_slam -laxom_slic -laxom_lumberjack -laxom_core
    ```
 3. [**Axom:**](https://github.com/LLNL/axom.git) Starting from the MFEM root
    directory (we assume this directory is named mfem):

miniapps/tools/convert-dc.cpp

@@ -18,7 +18,6 @@
 //
 // Currently supported data collection type options:
 //    visit:                VisItDataCollection (default)
-//    sidre or sidre_hdf5:  SidreDataCollection
 //    json:                 ConduitDataCollection w/ protocol json
 //    conduit_json:         ConduitDataCollection w/ protocol conduit_json
 //    conduit_bin:          ConduitDataCollection w/ protocol conduit_bin
@@ -52,14 +51,6 @@ DataCollection *create_data_collection(const std::string &dc_name,
       dc = new VisItDataCollection(MPI_COMM_WORLD, dc_name);
 #else
       dc = new VisItDataCollection(dc_name);
-#endif
-   }
-   else if ( dc_type == "sidre" || dc_type == "sidre_hdf5")
-   {
-#ifdef MFEM_USE_SIDRE
-      dc = new SidreDataCollection(dc_name);
-#else
-      MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for sidre support.");
 #endif
    }
    else if ( dc_type == "json" ||
@@ -140,7 +131,6 @@ int main(int argc, char *argv[])
    args.AddOption(&src_coll_type, "-st", "--source-type",
                   "Set the source data collection type. Options:\n"
                   "\t   visit:                VisItDataCollection (default)\n"
-                  "\t   sidre or sidre_hdf5:  SidreDataCollection\n"
                   "\t   json:                 ConduitDataCollection w/ protocol json\n"
                   "\t   conduit_json:         ConduitDataCollection w/ protocol conduit_json\n"
                   "\t   conduit_bin:          ConduitDataCollection w/ protocol conduit_bin\n"
@@ -152,7 +142,6 @@ int main(int argc, char *argv[])
    args.AddOption(&out_coll_type, "-ot", "--output-type",
                   "Set the output data collection type. Options:\n"
                   "\t   visit:                VisItDataCollection (default)\n"
-                  "\t   sidre or sidre_hdf5:  SidreDataCollection\n"
                   "\t   json:                 ConduitDataCollection w/ protocol json\n"
                   "\t   conduit_json:         ConduitDataCollection w/ protocol conduit_json\n"
                   "\t   conduit_bin:          ConduitDataCollection w/ protocol conduit_bin\n"

tests/unit/CMakeLists.txt

@@ -71,6 +71,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_ode2.cpp
   linalg/test_operator.cpp
   linalg/test_particlevector.cpp
+  linalg/test_petsc_nonlinear.cpp
   linalg/test_sparsesmoothers.cpp
   linalg/test_vector.cpp
   mesh/mesh_test_utils.cpp

tests/unit/fem/test_var_order.cpp

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

tests/unit/linalg/test_petsc_nonlinear.cpp

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

tests/unit/mesh/test_psubmesh.cpp

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