fix year

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

@@ -168,13 +168,10 @@ jobs:
           env
         shell: bash
 
-      # For info on Xcode see:
-      #  - https://github.com/actions/runner-images/issues/12541
-      #  - https://github.com/actions/runner-images/blob/releases/macos-15-arm64/20250811/images/macos/macos-15-arm64-Readme.md#xcode
       - name: Xcode version setup (MacOS)
         if: matrix.os == 'macos-latest'
         run: |
-          XCODE_PATH="/Applications/Xcode_16.4.app"
+          XCODE_PATH="/Applications/Xcode_15.3.app"
           echo "> sudo xcode-select -s ${XCODE_PATH}"
           sudo xcode-select -s ${XCODE_PATH}
           echo "> g++ -v"
@@ -292,10 +289,10 @@ jobs:
         run: |
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
-          brew install llvm@20 enzyme
-          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
-          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV
+          brew install llvm@19 enzyme
+          echo "LLVM_PREFIX=$(brew --prefix llvm@19)" >> $GITHUB_ENV
+          echo "OMPI_CC=$(brew --prefix llvm@19)/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$(brew --prefix llvm@19)/bin/clang++" >> $GITHUB_ENV
 
       # MFEM build and test
       - name: build

.gitlab-ci.yml

@@ -22,7 +22,7 @@ include:
 # the "needs" keyword and express the DAG of jobs for more efficiency.
 # - We use setup and setup_baseline phases to download content outside of mfem
 #   directory.
-# - Allocate/Release is where Dane resource are allocated/released once for all.
+# - Allocate/Release is where ruby resource are allocated/released once for all.
 # - Build and Test is where we build and MFEM for multiple toolchains.
 # - Baseline_checks gathers baseline-type test suites execution
 # - Baseline_publish, only available on master, allows to update baseline
@@ -53,7 +53,7 @@ variables:
   AUTOTEST_COMMIT: "YES"
 
 # Trigger subpipelines:
-dane-build-and-test:
+ruby-build-and-test:
   stage: sub-pipelines
   variables:
     # Explicitly pass down values that we want to be able to set when triggering
@@ -61,10 +61,10 @@ dane-build-and-test:
     AUTOTEST: "${AUTOTEST}"
     AUTOTEST_COMMIT: "${AUTOTEST_COMMIT}"
   trigger:
-    include: .gitlab/dane-build-and-test.yml
+    include: .gitlab/ruby-build-and-test.yml
     strategy: depend
 
-dane-baseline:
+ruby-baseline:
   stage: sub-pipelines
   variables:
     # Explicitly pass down values that we want to be able to set when triggering
@@ -73,7 +73,7 @@ dane-baseline:
     AUTOTEST: "${AUTOTEST}"
     AUTOTEST_COMMIT: "${AUTOTEST_COMMIT}"
   trigger:
-    include: .gitlab/dane-baseline.yml
+    include: .gitlab/ruby-baseline.yml
     strategy: depend
 
 lassen-build-and-test:

.gitlab/README.md

@@ -24,7 +24,7 @@ and `test type`.
 
 Machines typically include:
 
-* Dane: Intel Sapphire Rapids
+* Ruby: 2nd Gen Intel Xeon (Cascade Lake)
 * Lassen: Power9 + Nvidia GPU
 * Corona: AMD GPU
 
@@ -76,13 +76,13 @@ with a spack spec of MFEM, within the limits permitted by the MFEM spack
 package.
 
 In any build-and-test sub-pipeline a job basically consists in defining the
-spack spec to use. Adding a job on Dane for example resumes to:
+spack spec to use. Adding a job on ruby for example resumes to:
 
 ```yaml
 <job_name>:
   variables:
     SPEC: "<spack_spec>"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 ```
 
 The remaining and non trivial work is to make sure this spec is working. To

.gitlab/configs/common.yml

@@ -24,7 +24,7 @@ variables:
 # TODO: add a clean-up mechanism
   BUILD_ROOT: ${USER_CI_TOP_DIR}/${CI_PROJECT_NAME}-${MACHINE_NAME}-pipeline-${CI_PIPELINE_ID}
 
-# On LLNL's Dane, there is only one allocation shared among jobs in order to
+# On LLNL's ruby, there is only one allocation shared among jobs in order to
 # save time and resource. This allocation has to be uniquely named so that we
 # are sure to retrieve it.
   ALLOC_NAME: ${CI_PROJECT_NAME}_ci_${CI_PIPELINE_ID}

.gitlab/configs/ruby-config.yml

@@ -9,17 +9,17 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-# GitLab pipelines configurations for the Dane machine at LLNL
+# GitLab pipelines configurations for the Ruby machine at LLNL
 variables:
-  MACHINE_NAME: dane
+  MACHINE_NAME: ruby
 
-.on_dane:
+.on_ruby:
   tags:
     - shell
-    - dane
+    - ruby
   rules:
-    # Don't run dane jobs if...
-    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_DANE == "OFF"'
+    # Don't run ruby jobs if...
+    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_RUBY == "OFF"'
       when: never
     # Don't run autotest update if...
     - if: '$CI_JOB_NAME =~ /report/ && $AUTOTEST != "YES"'
@@ -40,14 +40,13 @@ variables:
     - when: on_success
 
 # Spack helped builds
-# Generic dane build job, extending build script
-.build_and_test_on_dane:
-  extends: [.on_dane]
+# Generic ruby build job, extending build script
+.build_and_test_on_ruby:
+  extends: [.on_ruby]
   stage: build_and_test
   script:
     # THREADS is used by 'tests/gitlab/build_and_test', run below
-    # Dane has 224 threads/node and we run 7 separate jobs: 224=7*32
-    - export THREADS=28
+    - export THREADS=16
     - echo ${ALLOC_NAME}
     - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
     - echo ${JOBID}

.gitlab/configs/setup-baseline.yml

@@ -18,7 +18,7 @@
 setup_baseline:
   tags:
     - shell
-    - dane
+    - ruby
   stage: setup
   variables:
     GIT_STRATEGY: none

.gitlab/configs/setup-build-and-test.yml

@@ -16,7 +16,7 @@
 setup:
   tags:
     - shell
-    - dane
+    - ruby
   stage: setup
   variables:
     GIT_STRATEGY: none

.gitlab/ruby-baseline.yml

@@ -19,8 +19,8 @@ stages:
   - cleanup
   - baseline_publish
 
-baselinecheck_mfem_intel_dane:
-  extends: [.on_dane]
+baselinecheck_mfem_intel_ruby:
+  extends: [.on_ruby]
   stage: baseline_check
   variables:
     # TPLS_DIR is used in .gitlab/scripts/baseline to provide the tpls location
@@ -31,8 +31,8 @@ baselinecheck_mfem_intel_dane:
   script:
     - echo ${BUILD_ROOT}
     - echo ${TPLS_DIR}
-    # Used by the tests in MFEM/tests, dane has 224 threads/node:
-    - export MFEM_TEST_NP=192
+    # Used by the tests in MFEM/tests:
+    - export MFEM_TEST_NP=48
     # The next script uses the following environment variables:
     # * BASELINE_TEST, SYS_TYPE, CI_PROJECT_DIR, ARTIFACTS_DIR,
     # * BUILD_ROOT, TPLS_DIR, MACHINE_NAME
@@ -44,7 +44,7 @@ baselinecheck_mfem_intel_dane:
   allow_failure: true
 
 cleanup:
-  extends: .on_dane
+  extends: .on_ruby
   stage: cleanup
   variables:
     GIT_STRATEGY: none
@@ -53,7 +53,7 @@ cleanup:
     - rm -rf "${BUILD_ROOT}" || true
 
 report_baseline:
-  extends: [.on_dane]
+  extends: [.on_ruby]
   stage: baseline_report
   script:
     - echo ${MACHINE_NAME}
@@ -113,8 +113,8 @@ report_baseline:
         exit $err
       ) 9> autotest.lock
 
-baselinepublish_mfem_dane:
-  extends: [.on_dane]
+baselinepublish_mfem_ruby:
+  extends: [.on_ruby]
   stage: baseline_publish
   rules:
     # - if: '$CI_COMMIT_BRANCH == "master" || $REBASELINE == "YES"'
@@ -129,5 +129,5 @@ baselinepublish_mfem_dane:
 
 include:
   - local: .gitlab/configs/common.yml
-  - local: .gitlab/configs/dane-config.yml
+  - local: .gitlab/configs/ruby-config.yml
   - local: .gitlab/configs/setup-baseline.yml

.gitlab/ruby-build-and-test.yml

@@ -19,54 +19,54 @@ stages:
 allocate_resource:
   variables:
     GIT_STRATEGY: none
-  extends: .on_dane
+  extends: .on_ruby
   stage: allocate_resource
   script:
     - echo ${ALLOC_NAME}
     - salloc --exclusive --nodes=1 --reservation=ci --time=60 --no-shell --job-name=${ALLOC_NAME}
   timeout: 6h
 
-# GitLab jobs for the Dane machine at LLNL
+# GitLab jobs for the Ruby machine at LLNL
 debug_ser_gcc_10:
   variables:
     SPEC: "%gcc@10.3.1 +debug~mpi"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 
 debug_par_gcc_10:
   variables:
     SPEC: "%gcc@10.3.1 +debug+mpi"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 
 opt_ser_gcc_10:
   variables:
     SPEC: "%gcc@10.3.1 ~mpi"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 
 opt_par_gcc_10:
   variables:
     SPEC: "%gcc@10.3.1"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 
 opt_par_gcc_10_sundials:
   variables:
     SPEC: "%gcc@10.3.1 +sundials"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 
 opt_par_gcc_10_petsc:
   variables:
     SPEC: "%gcc@10.3.1 +petsc ^petsc+mumps~superlu-dist"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 
 opt_par_gcc_10_pumi:
   variables:
     SPEC: "%gcc@10.3.1 +pumi"
-  extends: .build_and_test_on_dane
+  extends: .build_and_test_on_ruby
 
 # Release
 release_resource:
   variables:
     GIT_STRATEGY: none
-  extends: .on_dane
+  extends: .on_ruby
   stage: release_resource_and_report
   script:
     - echo ${ALLOC_NAME}
@@ -78,17 +78,17 @@ release_resource:
 report_job_success:
   stage: release_resource_and_report
   extends:
-    - .on_dane
+    - .on_ruby
     - .report_job_success
 
 report_job_failure:
   stage: release_resource_and_report
   extends:
-    - .on_dane
+    - .on_ruby
     - .report_job_failure
 
 include:
   - local: .gitlab/configs/common.yml
-  - local: .gitlab/configs/dane-config.yml
+  - local: .gitlab/configs/ruby-config.yml
   - local: .gitlab/configs/setup-build-and-test.yml
   - local: .gitlab/configs/report-build-and-test.yml

.gitlab/scripts/baseline

@@ -14,7 +14,7 @@
 # locals
 glob_err=${BASELINE_TEST}.err
 base=${BASELINE_TEST}-${SYS_TYPE}
-if [[ "${MACHINE_NAME}" == "dane" ]]; then
+if [[ "${MACHINE_NAME}" == "ruby" ]]; then
   base="${BASELINE_TEST}-${MACHINE_NAME}"
 fi
 base_diff=${base}.diff
@@ -31,7 +31,7 @@ cd tests
 mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
 
 # run
-if [[ "${MACHINE_NAME}" == "dane" ]]; then
+if [[ "${MACHINE_NAME}" == "ruby" ]]; then
   salloc --nodes=1 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
   salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"

.gitlab/scripts/rebaseline

@@ -11,7 +11,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-# There will be collision between corona and dane baselines.
+# There will be collision between corona and ruby baselines.
 # Once the corresponding files have been generated, we can switch to machine
 # specific ref.
 ARTIFACT_PATH=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/baseline-${SYS_TYPE}
@@ -21,7 +21,7 @@ PATCH_FILE=${ARTIFACT_PATH}.patch
 FULL_FILE=${ARTIFACT_PATH}.out
 DIFF_FILE=${ARTIFACT_PATH}.diff
 
-# There will be collision between corona and dane baselines.
+# There will be collision between corona and ruby baselines.
 # Once the corresponding files have been generated, we can switch to machine
 # specific ref.
 SAVED_NAME=baseline-${SYS_TYPE}.saved

INSTALL

@@ -668,7 +668,6 @@ The specific libraries and their options are:
 - OpenMP (optional), usually part of compiler, used when either MFEM_USE_OPENMP
   or MFEM_USE_LEGACY_OPENMP is set to YES.
   Options: OPENMP_OPT, OPENMP_LIB.
-  Versions: OpenMP >= 3.1 when MFEM_USE_OPENMP=YES.
 
 - High-resolution POSIX clocks: when using MFEM_TIMER_TYPE = 2, it may be
   necessary to link with a system library (e.g. librt.so).
@@ -848,7 +847,6 @@ The specific libraries and their options are:
 - HIP (optional), used when MFEM_USE_HIP = YES.
   URL: https://rocmdocs.amd.com
   Options: HIP_CXX, HIP_ARCH, HIP_OPT, HIP_LIB.
-  Versions: ROCm >= 5.6.1.
 
 - OCCA (optional), used when MFEM_USE_OCCA = YES.
   URL: https://libocca.org

config/config.hpp

@@ -23,14 +23,11 @@
 #include "_config.hpp"
 #endif
 
-#include <cstdint>
-#include <climits>
-
 namespace mfem
 {
 
 #if (defined(MFEM_USE_CUDA) && defined(__CUDACC__)) || \
-    (defined(MFEM_USE_HIP) && defined(__HIP__))
+    (defined(MFEM_USE_HIP) && defined(__HIPCC__))
 #define MFEM_HOST_DEVICE __host__ __device__
 #else
 #define MFEM_HOST_DEVICE

fem/CMakeLists.txt

@@ -165,6 +165,7 @@ set(SRCS
   hyperbolic.cpp
   integrator.cpp
   bounds.cpp
+  blockform.cpp
   )
 
 set(HDRS
@@ -277,6 +278,7 @@ set(HDRS
   hyperbolic.hpp
   integrator.hpp
   bounds.hpp
+  blockform.hpp
   )
 
 if (MFEM_USE_SIDRE)

fem/bilinearform.cpp

@@ -1275,22 +1275,6 @@ void BilinearForm::Update(FiniteElementSpace *nfes)
    height = width = fes->GetVSize();
 
    if (ext) { ext->Update(); }
-   for (int k = 0; k < domain_integs.Size(); ++k)
-   {
-      domain_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_integs[k]->Update();
-   }
-   for (int k = 0; k < interior_face_integs.Size(); ++k)
-   {
-      interior_face_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_face_integs[k]->Update();
-   }
 }
 
 void BilinearForm::SetDiagonalPolicy(DiagonalPolicy policy)
@@ -2353,31 +2337,6 @@ void MixedBilinearForm::Update()
    height = test_fes->GetVSize();
    width = trial_fes->GetVSize();
    if (ext) { ext->Update(); }
-
-   for (int k = 0; k < domain_integs.Size(); ++k)
-   {
-      domain_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_integs[k]->Update();
-   }
-   for (int k = 0; k < interior_face_integs.Size(); ++k)
-   {
-      interior_face_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_face_integs[k]->Update();
-   }
-   for (int k = 0; k < trace_face_integs.Size(); ++k)
-   {
-      trace_face_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_trace_face_integs.Size(); ++k)
-   {
-      boundary_trace_face_integs[k]->Update();
-   }
 }
 
 MixedBilinearForm::~MixedBilinearForm()

fem/bilinearform_ext.cpp

@@ -266,7 +266,11 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
 
       // Gather the attributes on the host from all the elements
       const Mesh &mesh = *trial_fes->GetMesh();
-      elem_attributes = &mesh.GetElementAttributes();
+      elem_attributes.SetSize(mesh.GetNE());
+      for (int i = 0; i < mesh.GetNE(); ++i)
+      {
+         elem_attributes[i] = mesh.GetAttribute(i);
+      }
    }
 
    // Construct face restriction operators only if the bilinear form has
@@ -325,7 +329,45 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
          bdr_face_dYdn.SetSize(bdr_face_restrict_lex->Height());
       }
 
-      bdr_face_attributes = &trial_fes->GetMesh()->GetBdrFaceAttributes();
+      const Mesh &mesh = *trial_fes->GetMesh();
+      // See LinearFormExtension::Update for explanation of f_to_be logic.
+      std::unordered_map<int,int> f_to_be;
+      for (int i = 0; i < mesh.GetNBE(); ++i)
+      {
+         const int f = mesh.GetBdrElementFaceIndex(i);
+         f_to_be[f] = i;
+      }
+      const int nf_bdr = trial_fes->GetNFbyType(FaceType::Boundary);
+      bdr_attributes.SetSize(nf_bdr);
+      int f_ind = 0;
+      int missing_bdr_elems = 0;
+      for (int f = 0; f < mesh.GetNumFaces(); ++f)
+      {
+         if (!mesh.GetFaceInformation(f).IsOfFaceType(FaceType::Boundary))
+         {
+            continue;
+         }
+         int attribute = 1; // default value
+         if (f_to_be.find(f) != f_to_be.end())
+         {
+            const int be = f_to_be[f];
+            attribute = mesh.GetBdrAttribute(be);
+         }
+         else
+         {
+            // If a boundary face does not correspond to the a boundary element,
+            // we assign it the default attribute of 1. We also generate a
+            // warning at runtime with the number of such missing elements.
+            ++missing_bdr_elems;
+         }
+         bdr_attributes[f_ind] = attribute;
+         ++f_ind;
+      }
+      if (missing_bdr_elems)
+      {
+         MFEM_WARNING("Missing " << missing_bdr_elems << " boundary elements "
+                      "for boundary faces.");
+      }
    }
 }
 
@@ -387,7 +429,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
          mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
          {
             const int attr = d_attr[e];
-            if (attr <= 0 || d_m[attr - 1] == 0)
+            if (d_m[attr - 1] == 0)
             {
                for (int i = 0; i < nd; ++i)
                {
@@ -408,7 +450,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
          for (int i = 0; i < iSz; ++i)
          {
             assemble_diagonal_with_markers(*integrators[i], elem_markers[i],
-                                           *elem_attributes, localY);
+                                           elem_attributes, localY);
          }
          const ElementRestriction* H1elem_restrict =
             dynamic_cast<const ElementRestriction*>(elem_restrict);
@@ -434,7 +476,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
       for (int i = 0; i < iSz; ++i)
       {
          assemble_diagonal_with_markers(*integrators[i], elem_markers[i],
-                                        *elem_attributes, y);
+                                        elem_attributes, y);
       }
    }
 
@@ -447,7 +489,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
       for (int i = 0; i < n_bdr_integs; ++i)
       {
          assemble_diagonal_with_markers(*bdr_integs[i], bdr_markers[i],
-                                        *bdr_face_attributes, bdr_face_Y);
+                                        bdr_attributes, bdr_face_Y);
       }
       bdr_face_restrict_lex->AddAbsMultTranspose(bdr_face_Y, y);
    }
@@ -546,7 +588,7 @@ void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
          for (int i = 0; i < iSz; ++i)
          {
             AddMultWithMarkers(*integrators[i], localX, elem_markers[i],
-                               *elem_attributes, false, localY, useAbs);
+                               elem_attributes, false, localY, useAbs);
          }
          if (H1elem_restrict && useAbs)
          {
@@ -648,8 +690,8 @@ void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
          }
          for (int i = 0; i < n_bdr_integs; ++i)
          {
-            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i],
-                               *bdr_face_attributes, false, bdr_face_Y);
+            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i], bdr_attributes,
+                               false, bdr_face_Y);
          }
          for (int i = 0; i < n_bdr_face_integs; ++i)
          {
@@ -657,14 +699,12 @@ void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
             {
                AddMultNormalDerivativesWithMarkers(
                   *bdr_face_integs[i], bdr_face_X, bdr_face_dXdn,
-                  bdr_face_markers[i], *bdr_face_attributes, bdr_face_Y,
-                  bdr_face_dYdn);
+                  bdr_face_markers[i], bdr_attributes, bdr_face_Y, bdr_face_dYdn);
             }
             else
             {
-               AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X,
-                                  bdr_face_markers[i], *bdr_face_attributes, false,
-                                  bdr_face_Y);
+               AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X, bdr_face_markers[i],
+                                  bdr_attributes, false, bdr_face_Y);
             }
          }
          bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
@@ -687,7 +727,7 @@ void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
       localY = 0.0;
       for (int i = 0; i < iSz; ++i)
       {
-         AddMultWithMarkers(*integrators[i], localX, elem_markers[i], *elem_attributes,
+         AddMultWithMarkers(*integrators[i], localX, elem_markers[i], elem_attributes,
                             true, localY);
       }
       elem_restrict->MultTranspose(localY, y);
@@ -734,14 +774,13 @@ void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
          bdr_face_Y = 0.0;
          for (int i = 0; i < n_bdr_integs; ++i)
          {
-            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i],
-                               *bdr_face_attributes, true, bdr_face_Y);
+            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i], bdr_attributes,
+                               true, bdr_face_Y);
          }
          for (int i = 0; i < n_bdr_face_integs; ++i)
          {
-            AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X,
-                               bdr_face_markers[i], *bdr_face_attributes, true,
-                               bdr_face_Y);
+            AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X, bdr_face_markers[i],
+                               bdr_attributes, true, bdr_face_Y);
          }
          bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
       }
@@ -765,7 +804,7 @@ static void AddWithMarkers_(
    mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
    {
       const int attr = d_attr[e];
-      if (attr <= 0 || d_m[attr - 1] == 0) { return; }
+      if (d_m[attr - 1] == 0) { return; }
       for (int i = 0; i < nd; ++i)
       {
          d_y(i, e) += d_x(i, e);
@@ -881,8 +920,7 @@ void EABilinearFormExtension::Assemble()
       {
          const int i = idx % sz;
          const int e = idx / sz;
-         const real_t val =
-            d_a[e] > 0 ? (d_m[d_a[e] - 1] ? d_ea_1(i, e) : 0) : 0;
+         const real_t val = d_m[d_a[e] - 1] ? d_ea_1(i, e) : 0.0;
          if (add)
          {
             d_ea_2(i, e) += val;
@@ -915,7 +953,7 @@ void EABilinearFormExtension::Assemble()
             ea_data_tmp.SetSize(ea_data.Size());
             integrators[i]->AssembleEA(*a->FESpace(), ea_data_tmp, false);
             add_with_markers(ea_data_tmp, ea_data, ne, *markers,
-                             *elem_attributes, add);
+                             elem_attributes, add);
          }
       }
    }
@@ -944,7 +982,7 @@ void EABilinearFormExtension::Assemble()
             ea_data_tmp.SetSize(ea_data_bdr.Size());
             bdr_integs[i]->AssembleEABoundary(*a->FESpace(), ea_data_tmp, add);
             add_with_markers(ea_data_tmp, ea_data_bdr, nf_bdr, *markers,
-                             *bdr_face_attributes, add);
+                             bdr_attributes, add);
          }
       }
    }
@@ -993,7 +1031,7 @@ void EABilinearFormExtension::Assemble()
                                                         ea_data_tmp,
                                                         add);
             add_with_markers(ea_data_tmp, ea_data_bdr, nf_bdr, *markers,
-                             *bdr_face_attributes, add);
+                             bdr_attributes, add);
          }
       }
    }

fem/bilinearform_ext.hpp

@@ -69,8 +69,7 @@ class PABilinearFormExtension : public BilinearFormExtension
 protected:
    const FiniteElementSpace *trial_fes, *test_fes; // Not owned
    /// Attributes of all mesh elements.
-   const Array<int> *elem_attributes; // Not owned
-   const Array<int> *bdr_face_attributes; // Not owned
+   Array<int> elem_attributes, bdr_attributes;
    mutable Vector tmp_evec; // Work array
    mutable Vector localX, localY;
    mutable Vector int_face_X, int_face_Y;

fem/bilininteg.cpp

@@ -21,11 +21,6 @@ using namespace std;
 namespace mfem
 {
 
-void BilinearFormIntegrator::Update()
-{
-   // default no-op
-}
-
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&)
 {
    MFEM_ABORT("BilinearFormIntegrator::AssemblePA(fes)\n"
@@ -3465,12 +3460,6 @@ real_t ElasticityIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
    return energy;
 }
 
-void DGTraceIntegrator::Update()
-{
-   qspace[0].reset();
-   qspace[1].reset();
-}
-
 void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
                                            const FiniteElement &el2,
                                            FaceElementTransformations &Trans,

fem/bilininteg.hpp

@@ -23,8 +23,6 @@
 
 namespace mfem
 {
-class QuadratureSpace;
-class FaceQuadratureSpace;
 
 /// Abstract base class BilinearFormIntegrator
 class BilinearFormIntegrator : public NonlinearFormIntegrator
@@ -46,10 +44,6 @@ public:
    // make sense for the action of the nonlinear operator (but they all make
    // sense for its Jacobian).
 
-   /// Signal this integrator that something about either the trial or test space has changed.
-   virtual void Update();
-
-
    /// Method defining partial assembly.
    /** The result of the partial assembly is stored internally so that it can be
        used later in the methods AddMultPA() and AddMultTransposePA(). */
@@ -818,7 +812,7 @@ protected:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetDim() == 1 && test_fe.GetDim() == 1 &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -890,7 +884,7 @@ protected:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
+      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -925,7 +919,7 @@ protected:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
+      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 
@@ -1606,7 +1600,7 @@ public:
    {
       return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL   &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 
@@ -1641,7 +1635,7 @@ public:
    {
       return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL   &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 
@@ -1675,7 +1669,7 @@ public:
    {
       return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -1766,7 +1760,7 @@ public:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -1799,7 +1793,7 @@ public:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR &&
               test_fe.GetDerivType()  == mfem::FiniteElement::DIV   );
    }
@@ -1838,7 +1832,7 @@ public:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::DIV   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::DIV    &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR &&
               test_fe.GetDerivType()  == mfem::FiniteElement::GRAD
              );
@@ -1979,7 +1973,7 @@ protected:
       const FiniteElement & test_fe) const override
    {
       return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 
@@ -3317,7 +3311,6 @@ protected:
    VectorCoefficient *u;
    real_t alpha, beta;
    // PA extension
-   std::unique_ptr<FaceQuadratureSpace> qspace[2];
    Vector pa_data;
    const DofToQuad *maps;             ///< Not owned
    const FaceGeometricFactors *geom;  ///< Not owned
@@ -3340,8 +3333,6 @@ public:
                      real_t a, real_t b)
    { rho = &rho_; u = &u_; alpha = a; beta = b; }
 
-   void Update() override;
-
    using BilinearFormIntegrator::AssembleFaceMatrix;
    void AssembleFaceMatrix(const FiniteElement &el1,
                            const FiniteElement &el2,

fem/blockform.cpp

@@ -0,0 +1,353 @@
+// 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.
+
+// Implementation of class BlockForm
+
+#include "../config/config.hpp"
+
+#include "fem.hpp"
+
+namespace mfem
+{
+
+void BlockForm::BuildProlongation()
+{
+   P = new BlockMatrix(dof_offsets, tdof_offsets);
+   R = new BlockMatrix(tdof_offsets, dof_offsets);
+   P->owns_blocks = 0;
+   R->owns_blocks = 0;
+
+   for (int i = 0; i<nblocks; i++)
+   {
+      const SparseMatrix * P_ = fes[i]->GetConformingProlongation();
+      P->SetBlock(i,i,const_cast<SparseMatrix*>(P_));
+      const SparseMatrix * R_ = fes[i]->GetRestrictionMatrix();
+      R->SetBlock(i,i,const_cast<SparseMatrix*>(R_));
+   }
+}
+
+void BlockForm::Finalize(int skip_zeros)
+{
+   if (mat) { mat->Finalize(skip_zeros); }
+   if (mat_e) { mat_e->Finalize(skip_zeros); }
+}
+
+void BlockForm::ConformingAssemble()
+{
+   Finalize(0);
+   if (!P) { BuildProlongation(); }
+
+   BlockMatrix * Pt = Transpose(*P);
+   BlockMatrix * PtA = mfem::Mult(*Pt, *mat);
+   // mat->owns_blocks = 0;
+   for (int i = 0; i<nblocks; i++)
+   {
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (mat->IsZeroBlock(i,j)) { continue; }
+         if (Pt->IsZeroBlock(i,i))
+         {
+            PtA->SetBlock(i,j,&mat->GetBlock(i,j));
+         }
+      }
+   }
+   delete mat;
+   if (mat_e)
+   {
+      BlockMatrix *PtAe = mfem::Mult(*Pt, *mat_e);
+      mat_e->owns_blocks = 0;
+      for (int i = 0; i<nblocks; i++)
+      {
+         for (int j = 0; j<nblocks; j++)
+         {
+            if (mat_e->IsZeroBlock(i,j)) { continue; }
+            SparseMatrix * tmp = &mat_e->GetBlock(i,j);
+            if (Pt->IsZeroBlock(i,i))
+            {
+               PtAe->SetBlock(i,j,tmp);
+            }
+            else
+            {
+               delete tmp;
+            }
+         }
+      }
+      delete mat_e;
+      mat_e = PtAe;
+      mat_e->owns_blocks = 1;
+   }
+   delete Pt;
+
+   mat = mfem::Mult(*PtA, *P);
+   PtA->owns_blocks = 0;
+   for (int i = 0; i<nblocks; i++)
+   {
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (PtA->IsZeroBlock(j,i)) { continue; }
+         SparseMatrix * tmp = &PtA->GetBlock(j,i);
+         if (P->IsZeroBlock(i,i))
+         {
+            mat->SetBlock(j,i,tmp);
+         }
+         else
+         {
+            delete tmp;
+         }
+      }
+   }
+   delete PtA;
+
+   if (mat_e)
+   {
+      BlockMatrix *PtAeP = mfem::Mult(*mat_e, *P);
+      mat_e->owns_blocks = 0;
+      for (int i = 0; i<nblocks; i++)
+      {
+         for (int j = 0; j<nblocks; j++)
+         {
+            if (mat_e->IsZeroBlock(j,i)) { continue; }
+            SparseMatrix * tmp = &mat_e->GetBlock(j,i);
+            if (P->IsZeroBlock(i,i))
+            {
+               PtAeP->SetBlock(j,i,tmp);
+            }
+            else
+            {
+               delete tmp;
+            }
+         }
+      }
+
+      delete mat_e;
+      mat_e = PtAeP;
+   }
+   height = mat->Height();
+   width = mat->Width();
+}
+
+
+BlockForm::BlockForm(const Array<FiniteElementSpace*> fes_ ): fes(
+      fes_)
+{
+   nblocks = fes.Size();
+   bforms.SetSize(nblocks,nblocks);
+   mforms.SetSize(nblocks,nblocks);
+   dof_offsets.Append(0);
+   tdof_offsets.Append(0);
+   for (int i = 0; i<nblocks; i++)
+   {
+      dof_offsets.Append(fes[i]->GetVSize());
+      tdof_offsets.Append(fes[i]->GetTrueVSize());
+      for (int j = 0; j<nblocks; j++)
+      {
+         bforms(i,j) = nullptr;
+         mforms(i,j) = nullptr;
+      }
+   }
+   dof_offsets.PartialSum();
+   tdof_offsets.PartialSum();
+   diag_policy = mfem::Operator::DIAG_ONE;
+
+}
+
+void BlockForm::SetBlock(BilinearForm * bform, int row_idx, int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   bforms(row_idx,col_idx) = bform;
+}
+void BlockForm::SetBlock(MixedBilinearForm * mform, int row_idx,
+                         int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   mforms(row_idx,col_idx) = mform;
+}
+
+/// Assemble the local matrix
+void BlockForm::Assemble(int skip_zeros)
+{
+   mat = new BlockMatrix(dof_offsets);
+   for (int i = 0; i<nblocks; i++)
+   {
+      int h = dof_offsets[i+1]-dof_offsets[i];
+      for (int j = 0; j<nblocks; j++)
+      {
+         int w = dof_offsets[j+1]-dof_offsets[j];
+         if (bforms(i,j))
+         {
+            bforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = bforms(i,j)->Height(), "inconsistent height of bilinear form");
+            MFEM_VERIFY(w = bforms(i,j)->Width(), "inconsistent width of bilinear form");
+            mat->SetBlock(i,j,&bforms(i,j)->SpMat());
+         }
+         else if (mforms(i,j))
+         {
+            mforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = mforms(i,j)->Height(),
+                        "inconsistent height of MixedBilinear form");
+            MFEM_VERIFY(w = mforms(i,j)->Width(),
+                        "inconsistent width of Mixedbilinear form");
+            mat->SetBlock(i,j,&mforms(i,j)->SpMat());
+         }
+         else
+         {
+            mat->SetBlock(i,j,nullptr);
+         }
+      }
+   }
+}
+
+void BlockForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
+                                 Vector & b,
+                                 OperatorHandle &A, Vector &X,
+                                 Vector &B, int copy_interior)
+{
+   FormSystemMatrix(ess_tdof_list, A);
+
+   if (!P)
+   {
+      EliminateVDofsInRHS(ess_tdof_list, x, b);
+      X.MakeRef(x, 0, x.Size());
+      B.MakeRef(b, 0, b.Size());
+      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+   }
+   else // non conforming space
+   {
+      B.SetSize(P->Width());
+
+      P->MultTranspose(b, B);
+      real_t *data = b.GetData();
+      Vector tmp;
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (P->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            B.SetVector(tmp,offset);
+         }
+      }
+
+      X.SetSize(R->Height());
+
+      R->Mult(x, X);
+      data = x.GetData();
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (R->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            X.SetVector(tmp,offset);
+         }
+      }
+
+      EliminateVDofsInRHS(ess_tdof_list, X, B);
+      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+   }
+}
+
+void BlockForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
+                                 OperatorHandle &A)
+{
+   if (!mat_e)
+   {
+      bool conforming = true;
+      for (int i = 0; i<nblocks; i++)
+      {
+         const SparseMatrix *P_ = fes[i]->GetConformingProlongation();
+         if (P_)
+         {
+            conforming = false;
+            break;
+         }
+      }
+      if (!conforming) { ConformingAssemble(); }
+      const int remove_zeros = 0;
+      EliminateVDofs(ess_tdof_list, diag_policy);
+      Finalize(remove_zeros);
+   }
+   A.Reset(mat, false);
+}
+
+void BlockForm::RecoverFEMSolution(const Vector &X, Vector &x)
+{
+   if (!P)
+   {
+      x.SyncMemory(X);
+   }
+   else
+   {
+      x.SetSize(P->Height());
+      P->Mult(X, x);
+      real_t *data = X.GetData();
+      Vector tmp;
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (P->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            x.SetVector(tmp,offset);
+         }
+      }
+   }
+}
+
+void BlockForm::EliminateVDofs(const Array<int> &vdofs,
+                               Operator::DiagonalPolicy dpolicy)
+{
+   if (mat_e == NULL)
+   {
+      Array<int> offsets;
+
+      offsets.MakeRef( (P) ? tdof_offsets : dof_offsets);
+
+      mat_e = new BlockMatrix(offsets);
+      mat_e->owns_blocks = 1;
+      for (int i = 0; i<mat_e->NumRowBlocks(); i++)
+      {
+         int h = offsets[i+1] - offsets[i];
+         for (int j = 0; j<mat_e->NumColBlocks(); j++)
+         {
+            int w = offsets[j+1] - offsets[j];
+            mat_e->SetBlock(i,j,new SparseMatrix(h, w));
+         }
+      }
+   }
+
+   mat->EliminateRowCols(vdofs,mat_e,diag_policy);
+}
+
+void BlockForm::EliminateVDofsInRHS(
+   const Array<int> &vdofs, const Vector &x, Vector &b)
+{
+   mat_e->AddMult(x,b,-1.);
+   mat->PartMult(vdofs,x,b);
+}
+
+BlockForm::~BlockForm()
+{
+   delete mat_e;
+   mat_e = nullptr;
+   delete mat;
+   mat = nullptr;
+   delete P;
+   delete R;
+}
+
+};

fem/blockform.hpp

@@ -0,0 +1,88 @@
+// 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_BLOCKFORM
+#define MFEM_BLOCKFORM
+
+#include "../config/config.hpp"
+
+#include "fespace.hpp"
+#include "gridfunc.hpp"
+#include "bilinearform.hpp"
+
+namespace mfem
+{
+// square block forms
+class BlockForm
+{
+
+private:
+   int height, width;
+   int nblocks;
+   Array<int> dof_offsets;
+   Array<int> tdof_offsets;
+   // BilinearForms
+   Array2D<BilinearForm * > bforms;
+   Array2D<MixedBilinearForm * > mforms;
+   Array<FiniteElementSpace *> fes;
+
+   // Block Prolongation
+   BlockMatrix * P = nullptr;
+   // Block Restriction
+   BlockMatrix * R = nullptr;
+
+   BlockMatrix * mat = nullptr;
+   BlockMatrix * mat_e = nullptr;
+
+   void BuildProlongation();
+   void ConformingAssemble();
+
+   mfem::Operator::DiagonalPolicy diag_policy;
+
+public:
+
+   BlockForm(const Array<FiniteElementSpace*> pfes_ );
+
+   void SetBlock(BilinearForm * bform, int row_idx, int col_idx);
+   void SetBlock(MixedBilinearForm * mform, int row_idx, int col_idx);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector & b,
+                         OperatorHandle &A, Vector &X,
+                         Vector &B, int copy_interior = 0);
+
+   void FormSystemMatrix(const Array<int> &ess_tdof_list,
+                         OperatorHandle &A);
+
+   void RecoverFEMSolution(const Vector &X, Vector &x);
+
+   ///  Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   void EliminateVDofs(const Array<int> &vdofs,
+                       Operator::DiagonalPolicy dpolicy = Operator::DIAG_ONE);
+
+   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
+
+   void SetDiagonalPolicy(Operator::DiagonalPolicy policy)
+   {
+      diag_policy = policy;
+   }
+
+   /// Destroys bilinear form.
+   ~BlockForm();
+};
+
+}
+
+#endif

fem/dfem/util.hpp

@@ -327,8 +327,8 @@ void print_mpi_sync(const std::string& msg)
    // First gather string lengths
    size_t msg_len = msg.length();
    std::vector<size_t> lengths(nranks);
-   MPI_Gather(&msg_len, 1, MPITypeMap<size_t>::mpi_type,
-              lengths.data(), 1, MPITypeMap<size_t>::mpi_type,
+   MPI_Gather(&msg_len, 1, MPI_INT,
+              lengths.data(), 1, MPI_INT,
               0, MPI_COMM_WORLD);
 
    if (myrank == 0)
@@ -568,7 +568,7 @@ struct ThreadBlocks
    int z = 1;
 };
 
-#if defined(MFEM_USE_CUDA_OR_HIP)
+#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
 template <typename func_t>
 __global__ void forall_kernel_shmem(func_t f, int n)
 {
@@ -591,7 +591,7 @@ void forall(func_t f,
    if (Device::Allows(Backend::CUDA_MASK) ||
        Device::Allows(Backend::HIP_MASK))
    {
-#if defined(MFEM_USE_CUDA_OR_HIP)
+#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
       // int gridsize = (N + Z - 1) / Z;
       int num_bytes = num_shmem * sizeof(decltype(shmem));
       dim3 block_size(blocks.x, blocks.y, blocks.z);

fem/fem.hpp

@@ -50,8 +50,8 @@
 #include "dgmassinv.hpp"
 #include "hyperbolic.hpp"
 #include "bounds.hpp"
-
 #include "dfem/doperator.hpp"
+#include "blockform.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"

fem/integ/bilininteg_dgtrace_pa.cpp

@@ -139,6 +139,8 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                          Device::GetDeviceMemoryType() : pa_mt;
 
+   nf = fes.GetNFbyType(type);
+   if (nf==0) { return; }
    // Assumes tensor-product elements
    Mesh *mesh = fes.GetMesh();
    const FiniteElement &el = *fes.GetTypicalTraceElement();
@@ -146,17 +148,6 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
                                IntRule:
                                &GetRule(el.GetGeomType(), el.GetOrder(),
                                         *mesh->GetTypicalElementTransformation());
-
-   if (!qspace[static_cast<int>(type)])
-   {
-      qspace[static_cast<int>(type)].reset(
-         new FaceQuadratureSpace(*mesh, *ir, type));
-   }
-
-   FaceQuadratureSpace& qs = *qspace[static_cast<int>(type)];
-   nf = qs.GetNumFaces();
-   if (nf==0) { return; }
-
    const int symmDims = 4;
    nq = ir->GetNPoints();
    dim = mesh->Dimension();
@@ -168,6 +159,8 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
+
+   FaceQuadratureSpace qs(*mesh, *ir, type);
    CoefficientVector vel(*u, qs, CoefficientStorage::COMPRESSED);
 
    CoefficientVector r(qs, CoefficientStorage::COMPRESSED);

fem/integ/bilininteg_mass_pa.cpp

@@ -59,23 +59,26 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 
    QuadratureSpace qs(*mesh, *ir);
    CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);
+
+   const int NE = ne;
+   const int Q1D = quad1D;
+   const int NQ = static_cast<int>(std::pow(Q1D, dim));
+   const bool const_c = coeff.Size() == 1;
+   const bool by_val = map_type == FiniteElement::VALUE;
+   const auto W = Reshape(ir->GetWeights().Read(), NQ);
+   const auto J = Reshape(geom->detJ.Read(), NQ, NE);
+   const auto C = const_c ? Reshape(coeff.Read(), 1, 1) :
+                  Reshape(coeff.Read(), NQ,NE);
+   auto v = Reshape(pa_data.Write(), NQ, NE);
+   mfem::forall_2D(NE, NQ, 1, [=] MFEM_HOST_DEVICE (int e)
    {
-      const int NE = ne;
-      const int NQ = nq;
-      const bool const_c = coeff.Size() == 1;
-      const bool by_val = map_type == FiniteElement::VALUE;
-      const auto W = Reshape(ir->GetWeights().Read(), NQ);
-      const auto J = Reshape(geom->detJ.Read(), NQ, NE);
-      const auto C =
-         const_c ? Reshape(coeff.Read(), 1, 1) : Reshape(coeff.Read(), NQ, NE);
-      auto v = Reshape(pa_data.Write(), NQ, NE);
-      mfem::forall(NQ, NE, [=] MFEM_HOST_DEVICE(int q, int e)
+      MFEM_FOREACH_THREAD(i, x, NQ)
       {
-         const real_t detJ = J(q, e);
-         const real_t coeff = const_c ? C(0, 0) : C(q, e);
-         v(q, e) = W(q) * coeff * (by_val ? detJ : 1.0 / detJ);
-      });
-   }
+         const real_t detJ = J(i,e);
+         const real_t coeff = const_c ? C(0,0) : C(i,e);
+         v(i,e) =  W(i) * coeff * (by_val ? detJ : 1.0/detJ);
+      }
+   });
 }
 
 void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
@@ -106,22 +109,50 @@ void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
    CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);
 
    const int NE = ne;
-   const int NQ = nq;
+   const int Q1D = quad1D;
    const bool const_c = coeff.Size() == 1;
    const bool by_val = map_type == FiniteElement::VALUE;
+   if (dim==1)
    {
-      const auto W = Reshape(ir->GetWeights().Read(), NQ);
-      const auto J = Reshape(face_geom->detJ.Read(), NQ, NE);
-      const auto C = const_c ? Reshape(coeff.Read(), 1, 1)
-                     : Reshape(coeff.Read(), NQ, NE);
-      auto v = Reshape(pa_data.Write(), NQ, NE);
-      mfem::forall(NQ, NE, [=] MFEM_HOST_DEVICE(int q, int e)
+      const auto W = Reshape(ir->GetWeights().Read(), Q1D);
+      const auto J = Reshape(face_geom->detJ.Read(), Q1D, NE);
+      const auto C = const_c ? Reshape(coeff.Read(), 1, 1) :
+                     Reshape(coeff.Read(), Q1D, NE);
+      auto v = Reshape(pa_data.Write(), Q1D, NE);
+      mfem::forall_2D(NE, Q1D, 1, [=] MFEM_HOST_DEVICE (int e)
       {
-         const real_t detJ = J(q, e);
-         const real_t coeff = const_c ? C(0, 0) : C(q, e);
-         v(q, e) = W(q) * coeff * (by_val ? detJ : 1.0 / detJ);
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            const real_t detJ = J(qx,e);
+            const real_t coeff = const_c ? C(0,0) : C(qx,e);
+            v(qx,e) =  W(qx) * coeff * (by_val ? detJ : 1.0/detJ);
+         }
       });
    }
+   else if (dim==2)
+   {
+      const auto W = Reshape(ir->GetWeights().Read(), Q1D,Q1D);
+      const auto J = Reshape(face_geom->detJ.Read(), Q1D,Q1D,NE);
+      const auto C = const_c ? Reshape(coeff.Read(), 1,1,1) :
+                     Reshape(coeff.Read(), Q1D,Q1D,NE);
+      auto v = Reshape(pa_data.Write(), Q1D,Q1D, NE);
+      mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+      {
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            MFEM_FOREACH_THREAD(qy,y,Q1D)
+            {
+               const real_t detJ = J(qx,qy,e);
+               const real_t coeff = const_c ? C(0,0,0) : C(qx,qy,e);
+               v(qx,qy,e) =  W(qx,qy) * coeff * (by_val ? detJ : 1.0/detJ);
+            }
+         }
+      });
+   }
+   else
+   {
+      MFEM_ABORT("Not supported.");
+   }
 }
 
 void MassIntegrator::AssembleDiagonalPA(Vector &diag)

fem/linearform_ext.cpp

@@ -51,7 +51,7 @@ void LinearFormExtension::Assemble()
       {
          // scan the attributes to set the markers to 0 or 1
          const int NE = fes.GetNE();
-         const auto attr = attributes->Read();
+         const auto attr = attributes.Read();
          const auto dimk = domain_integs_marker_k->Read();
          auto markers_w = markers.Write();
          mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
@@ -93,14 +93,13 @@ void LinearFormExtension::Assemble()
       else
       {
          // scan the attributes to set the markers to 0 or 1
-         const int NBE = bdr_face_attributes->Size();
-         const auto attr = bdr_face_attributes->Read();
+         const int NBE = bdr_attributes.Size();
+         const auto attr = bdr_attributes.Read();
          const auto attr_markers = boundary_integs_marker_k->Read();
          auto markers_w = bdr_markers.Write();
-         mfem::forall(NBE, [=] MFEM_HOST_DEVICE(int e)
+         mfem::forall(NBE, [=] MFEM_HOST_DEVICE (int e)
          {
-            markers_w[e] =
-               attr[e] > 0 ? (attr_markers[attr[e] - 1] == 1) : false;
+            markers_w[e] = attr_markers[attr[e]-1] == 1;
          });
       }
 
@@ -126,7 +125,8 @@ void LinearFormExtension::Update()
       //markers.UseDevice(true);
 
       // Gather the attributes on the host from all the elements
-      attributes = &mesh.GetElementAttributes();
+      attributes.SetSize(NE);
+      for (int i = 0; i < NE; ++i) { attributes[i] = mesh.GetAttribute(i); }
 
       elem_restrict_lex = fes.GetElementRestriction(ordering);
       MFEM_VERIFY(elem_restrict_lex, "Element restriction not available");
@@ -136,12 +136,34 @@ void LinearFormExtension::Update()
 
    if (lf->boundary_integs.Size() > 0)
    {
-      bdr_face_attributes = &mesh.GetBdrFaceAttributes();
-
-      const int nf_bdr = bdr_face_attributes->Size();
+      const int nf_bdr = fes.GetNFbyType(FaceType::Boundary);
       bdr_markers.SetSize(nf_bdr);
       // bdr_markers.UseDevice(true);
 
+      // The face restriction will give us "face E-vectors" on the boundary that
+      // are numbered in the order of the faces of mesh. This numbering will be
+      // different than the numbering of the boundary elements. We compute
+      // mappings so that the array `bdr_attributes[i]` gives the boundary
+      // attribute of the `i`th boundary face in the mesh face order.
+      std::unordered_map<int,int> f_to_be;
+      for (int i = 0; i < mesh.GetNBE(); ++i)
+      {
+         const int f = mesh.GetBdrElementFaceIndex(i);
+         f_to_be[f] = i;
+      }
+      MFEM_VERIFY(size_t(nf_bdr) == f_to_be.size(), "Incompatible sizes");
+      bdr_attributes.SetSize(nf_bdr);
+      int f_ind = 0;
+      for (int f = 0; f < mesh.GetNumFaces(); ++f)
+      {
+         if (f_to_be.find(f) != f_to_be.end())
+         {
+            const int be = f_to_be[f];
+            bdr_attributes[f_ind] = mesh.GetBdrAttribute(be);
+            ++f_ind;
+         }
+      }
+
       bdr_restrict_lex =
          dynamic_cast<const FaceRestriction*>(
             fes.GetFaceRestriction(ordering, FaceType::Boundary,

fem/linearform_ext.hpp

@@ -25,8 +25,7 @@ class LinearForm;
 class LinearFormExtension
 {
    /// Attributes of all mesh elements.
-   const Array<int> *attributes; // Not owned
-   const Array<int> *bdr_face_attributes; // Not owned
+   Array<int> attributes, bdr_attributes;
 
    /// Temporary markers for device kernels.
    Array<int> markers, bdr_markers;

fem/qfunction.hpp

@@ -214,16 +214,16 @@ public:
 inline void QuadratureFunction::GetValues(
    int idx, Vector &values)
 {
-   const int s_offset = qspace->Offset(idx);
-   const int sl_size = qspace->Offset(idx + 1) - s_offset;
+   const int s_offset = qspace->offsets[idx];
+   const int sl_size = qspace->offsets[idx+1] - s_offset;
    values.MakeRef(*this, vdim*s_offset, vdim*sl_size);
 }
 
 inline void QuadratureFunction::GetValues(
    int idx, Vector &values) const
 {
-   const int s_offset = qspace->Offset(idx);
-   const int sl_size = qspace->Offset(idx + 1) - s_offset;
+   const int s_offset = qspace->offsets[idx];
+   const int sl_size = qspace->offsets[idx+1] - s_offset;
    values.SetSize(vdim*sl_size);
    values.HostWrite();
    const real_t *q = HostRead() + vdim*s_offset;
@@ -236,14 +236,14 @@ inline void QuadratureFunction::GetValues(
 inline void QuadratureFunction::GetValues(
    int idx, const int ip_num, Vector &values)
 {
-   const int s_offset = qspace->Offset(idx) * vdim + ip_num * vdim;
+   const int s_offset = qspace->offsets[idx] * vdim + ip_num * vdim;
    values.MakeRef(*this, s_offset, vdim);
 }
 
 inline void QuadratureFunction::GetValues(
    int idx, const int ip_num, Vector &values) const
 {
-   const int s_offset = qspace->Offset(idx) * vdim + ip_num * vdim;
+   const int s_offset = qspace->offsets[idx] * vdim + ip_num * vdim;
    values.SetSize(vdim);
    values.HostWrite();
    const real_t *q = HostRead() + s_offset;
@@ -256,8 +256,8 @@ inline void QuadratureFunction::GetValues(
 inline void QuadratureFunction::GetValues(
    int idx, DenseMatrix &values)
 {
-   const int s_offset = qspace->Offset(idx);
-   const int sl_size = qspace->Offset(idx + 1) - s_offset;
+   const int s_offset = qspace->offsets[idx];
+   const int sl_size = qspace->offsets[idx+1] - s_offset;
    // Make the values matrix memory an alias of the quadrature function memory
    Memory<real_t> &values_mem = values.GetMemory();
    values_mem.Delete();
@@ -268,8 +268,8 @@ inline void QuadratureFunction::GetValues(
 inline void QuadratureFunction::GetValues(
    int idx, DenseMatrix &values) const
 {
-   const int s_offset = qspace->Offset(idx);
-   const int sl_size = qspace->Offset(idx + 1) - s_offset;
+   const int s_offset = qspace->offsets[idx];
+   const int sl_size = qspace->offsets[idx+1] - s_offset;
    values.SetSize(vdim, sl_size);
    values.HostWrite();
    const real_t *q = HostRead() + vdim*s_offset;

fem/qspace.cpp

@@ -17,9 +17,8 @@ namespace mfem
 {
 
 QuadratureSpaceBase::QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                                         const IntegrationRule &ir,
-                                         QSpaceStorage storage)
-   : mesh(mesh_), order(ir.GetOrder()), storage(storage)
+                                         const IntegrationRule &ir)
+   : mesh(mesh_), order(ir.GetOrder())
 {
    for (int g = 0; g < Geometry::NumGeom; g++)
    {
@@ -97,37 +96,27 @@ void QuadratureSpaceBase::Integrate(VectorCoefficient &coeff,
 
 void QuadratureSpace::ConstructOffsets()
 {
-   const int num_elem = ne;
-
-   if (storage == QSpaceStorage::COMPRESSED &&
-       mesh.GetNumGeometries(mesh.Dimension()) == 1)
+   const int num_elem = mesh.GetNE();
+   offsets.SetSize(num_elem + 1);
+   int offset = 0;
+   for (int i = 0; i < num_elem; i++)
    {
-      Array<Geometry::Type> geoms;
-      mesh.GetGeometries(mesh.Dimension(), geoms);
-      offsets.SetSize(1);
-      offsets.HostWrite();
-      offsets[0] = int_rule[geoms[0]]->GetNPoints();
-      size = num_elem * offsets[0];
-   }
-   else
-   {
-      offsets.SetSize(num_elem + 1);
-      int offset = 0;
-      for (int i = 0; i < num_elem; i++)
-      {
-         offsets[i] = offset;
-         const Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule.");
-         offset += int_rule[geom]->GetNPoints();
-      }
-      offsets[num_elem] = offset;
-      size = offsets.Last();
+      offsets[i] = offset;
+      int geom = mesh.GetElementBaseGeometry(i);
+      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule.");
+      offset += int_rule[geom]->GetNPoints();
    }
+   offsets[num_elem] = size = offset;
 }
 
-QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in,
-                                 QSpaceStorage storage)
-   : QuadratureSpaceBase(*mesh_, 0, storage)
+void QuadratureSpace::Construct()
+{
+   ConstructIntRules(mesh.Dimension());
+   ConstructOffsets();
+}
+
+QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
+   : QuadratureSpaceBase(*mesh_)
 {
    const char *msg = "invalid input stream";
    std::string ident;
@@ -146,24 +135,15 @@ QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in,
       return;
    }
 
-   ne = mesh.GetNE();
-   ConstructIntRules(mesh.Dimension());
+   Construct();
 }
 
-QuadratureSpace::QuadratureSpace(Mesh *mesh_, int order_, QSpaceStorage storage)
-   : QuadratureSpaceBase(*mesh_, order_, storage)
-{
-   ne = mesh.GetNE();
-   ConstructIntRules(mesh.Dimension());
-}
-
-QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                                 QSpaceStorage storage)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir, storage)
+QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir)
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) <= 1,
                "Constructor not valid for mixed meshes");
-   ne = mesh.GetNE();
+   ConstructOffsets();
 }
 
 void QuadratureSpace::Save(std::ostream &os) const
@@ -185,53 +165,54 @@ const Vector &QuadratureSpace::GetGeometricFactorWeights() const
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, int order_,
-                                         FaceType face_type_,
-                                         QSpaceStorage storage)
-   : QuadratureSpaceBase(mesh_, order_, storage), face_type(face_type_),
-     face_indices(mesh.GetFaceIndices(face_type_)),
-     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
+                                         FaceType face_type_)
+   : QuadratureSpaceBase(mesh_, order_),
+     face_type(face_type_),
+     num_faces(mesh.GetNFbyType(face_type))
 {
-   ne = face_indices.Size();
-   ConstructIntRules(mesh.Dimension() - 1);
+   Construct();
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                                         FaceType face_type_,
-                                         QSpaceStorage storage)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir, storage),
-     face_type(face_type_), face_indices(mesh.GetFaceIndices(face_type_)),
-     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
+                                         FaceType face_type_)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir),
+     face_type(face_type_),
+     num_faces(mesh.GetNFbyType(face_type))
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension() - 1) <= 1,
                "Constructor not valid for mixed meshes");
-   ne = face_indices.Size();
+   ConstructOffsets();
 }
 
 void FaceQuadratureSpace::ConstructOffsets()
 {
-   if (storage == QSpaceStorage::COMPRESSED &&
-       mesh.GetNumGeometries(mesh.Dimension() - 1) == 1)
+   face_indices.SetSize(num_faces);
+   offsets.SetSize(num_faces + 1);
+   int offset = 0;
+   int f_idx = 0;
+   for (int i = 0; i < mesh.GetNumFacesWithGhost(); i++)
    {
-      Array<Geometry::Type> geoms;
-      mesh.GetGeometries(mesh.Dimension() - 1, geoms);
-      offsets.SetSize(1);
-      offsets.HostWrite();
-      offsets[0] = int_rule[geoms[0]]->GetNPoints();
-      size = ne * offsets[0];
-   }
-   else
-   {
-      offsets.SetSize(face_indices.Size() + 1);
-      int offset = 0;
-      for (int i = 0; i < mesh.GetNFbyType(face_type); ++i)
+      const Mesh::FaceInformation face = mesh.GetFaceInformation(i);
+      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(face_type))
       {
-         offsets[i] = offset;
-         Geometry::Type geom = mesh.GetFaceGeometry(face_indices[i]);
-         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
-         offset += int_rule[geom]->GetNPoints();
+         continue;
       }
-      offsets[face_indices.Size()] = size = offset;
+      face_indices[f_idx] = i;
+      face_indices_inv[i] = f_idx;
+      offsets[f_idx] = offset;
+      Geometry::Type geom = mesh.GetFaceGeometry(i);
+      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
+      offset += int_rule[geom]->GetNPoints();
+
+      f_idx++;
    }
+   offsets[num_faces] = size = offset;
+}
+
+void FaceQuadratureSpace::Construct()
+{
+   ConstructIntRules(mesh.Dimension() - 1);
+   ConstructOffsets();
 }
 
 int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const

fem/qspace.hpp

@@ -19,49 +19,35 @@
 namespace mfem
 {
 
-enum class QSpaceStorage
-{
-   FULL,
-   COMPRESSED
-};
-
 /// Abstract base class for QuadratureSpace and FaceQuadratureSpace.
 /** This class represents the storage layout for QuadratureFunction%s, that may
     be defined either on mesh elements or mesh faces. */
 class QuadratureSpaceBase
 {
 protected:
+   friend class QuadratureFunction; // Uses the offsets.
+
    Mesh &mesh; ///< The underlying mesh.
    int order; ///< The order of integration rule.
-   int size = -1; ///< Total number of quadrature points. -1 indicates
-   ///< offsets/size not computed yet.
-   int ne; ///< Actual number of entities
+   int size; ///< Total number of quadrature points.
    mutable Vector weights; ///< Integration weights.
    mutable long nodes_sequence = 0; ///< Nodes counter for cache invalidation.
 
-   QSpaceStorage storage;
-
-   /// @brief Entity quadrature point offset array.
-   ///
-   /// Supports a constant compression scheme for meshes which have a single
-   /// geometry type. When compressed, will have a single value. The true offset
-   /// can be computed as i * offsets[0], where i is the entity index. Otherwise
-   /// has size num_entities + 1. Lazily constructed.
+   /// @brief Entity quadrature point offset array, of size num_entities + 1.
    ///
+   /// The quadrature point values for entity i are stored in the indices between
+   /// offsets[i] and offsets[i+1].
    Array<int> offsets;
    /// The quadrature rules used for each geometry type.
    const IntegrationRule *int_rule[Geometry::NumGeom];
 
    /// Protected constructor. Used by derived classes.
-   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED)
-      : mesh(mesh_), order(order_), storage(storage)
-   {}
+   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0)
+      : mesh(mesh_), order(order_) { }
 
    /// Protected constructor. Used by derived classes.
    QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                       const IntegrationRule &ir,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+                       const IntegrationRule &ir);
 
    /// Fill the @ref int_rule array for each geometry type using @ref order.
    void ConstructIntRules(int dim);
@@ -72,55 +58,15 @@ protected:
    /// Compute the integration weights.
    void ConstructWeights() const;
 
-   virtual void ConstructOffsets() = 0;
-
 public:
-   QSpaceStorage StorageType() const { return storage; }
-
-   /// @brief Gets the offset for a given entity @a idx.
-   ///
-   /// The quadrature point values for entity i are stored in the indices
-   /// between Offset(i) and Offset(i+1)
-   int Offset(int idx) const
-   {
-      if (size < 0)
-      {
-         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
-      }
-      return (offsets.Size() == 1) ? (idx * offsets[0]) : offsets[idx];
-   }
-
-   /// @brief Entity quadrature point offset array.
-   ///
-   /// Supports a constant compression scheme for meshes which have a single
-   /// geometry type. When compressed, will have a single value. The true offset
-   /// can be computed as i * offsets[0], where i is the entity index. Otherwise
-   /// has size num_entities + 1.
-   ///
-   const Array<int> &Offsets() const
-   {
-      if (size < 0)
-      {
-         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
-      }
-      return offsets;
-   }
-
    /// Return the total number of quadrature points.
-   int GetSize() const
-   {
-      if (size < 0)
-      {
-         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
-      }
-      return size;
-   }
+   int GetSize() const { return size; }
 
    /// Return the order of the quadrature rule(s) used by all elements.
    int GetOrder() const { return order; }
 
    /// Return the number of entities.
-   int GetNE() const { return ne; }
+   int GetNE() const { return offsets.Size() - 1; }
 
    /// Returns the mesh.
    inline Mesh *GetMesh() const { return &mesh; }
@@ -174,20 +120,19 @@ class QuadratureSpace : public QuadratureSpaceBase
 {
 protected:
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets() override;
+   void ConstructOffsets();
+   void Construct();
 public:
    /// Create a QuadratureSpace based on the global rules from #IntRules.
-   QuadratureSpace(Mesh *mesh_, int order_,
-                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   QuadratureSpace(Mesh *mesh_, int order_)
+      : QuadratureSpaceBase(*mesh_, order_) { Construct(); }
 
    /// @brief Create a QuadratureSpace with an IntegrationRule, valid only when
    /// the mesh has one element type.
-   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir);
 
    /// Read a QuadratureSpace from the stream @a in.
-   QuadratureSpace(Mesh *mesh_, std::istream &in,
-                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   QuadratureSpace(Mesh *mesh_, std::istream &in);
 
    /// Returns number of elements in the mesh.
    inline int GetNE() const { return mesh.GetNE(); }
@@ -224,29 +169,29 @@ public:
 class FaceQuadratureSpace : public QuadratureSpaceBase
 {
    FaceType face_type; ///< Is the space defined on interior or boundary faces?
+   const int num_faces; ///< Number of faces.
 
    /// Map from boundary or interior face indices to mesh face indices.
-   const Array<int> &face_indices;
+   Array<int> face_indices;
 
    /// Inverse of the map @a face_indices.
-   const std::unordered_map<int,int> &face_indices_inv;
+   std::unordered_map<int,int> face_indices_inv;
 
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets() override;
+   void ConstructOffsets();
+   void Construct();
 
 public:
    /// Create a FaceQuadratureSpace based on the global rules from #IntRules.
-   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_);
 
    /// @brief Create a FaceQuadratureSpace with an IntegrationRule, valid only
    /// when the mesh has one type of face geometry.
    FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                       FaceType face_type_,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+                       FaceType face_type_);
 
    /// Returns number of faces in the mesh.
-   inline int GetNumFaces() const { return face_indices.Size(); }
+   inline int GetNumFaces() const { return num_faces; }
 
    /// Returns the face type (boundary or interior).
    FaceType GetFaceType() const { return face_type; }

general/array.hpp

@@ -211,9 +211,6 @@ public:
    /// Delete the first entry with value == 'el'.
    inline void DeleteFirst(const T &el);
 
-   /// Delete entries at @a indices, and resize.
-   inline void DeleteAt(const Array<int> &indices);
-
    /// Delete the whole array.
    inline void DeleteAll();
 
@@ -938,30 +935,6 @@ inline void Array<T>::DeleteFirst(const T &el)
    }
 }
 
-template <class T>
-inline void Array<T>::DeleteAt(const Array<int> &indices)
-{
-   // Make a copy of the indices, sorted.
-   Array<int> sorted_indices(indices);
-   sorted_indices.Sort();
-
-   int rm_count = 0;
-   for (int i = 0; i < size; i++)
-   {
-      if (rm_count < sorted_indices.Size() && i == sorted_indices[rm_count])
-      {
-         rm_count++;
-      }
-      else
-      {
-         data[i-rm_count] = data[i]; // shift data rm_count
-      }
-   }
-
-   // Resize to remove tail
-   size -= rm_count;
-}
-
 template <class T>
 inline void Array<T>::DeleteAll()
 {

general/backends.hpp

@@ -14,7 +14,7 @@
 
 #include "../config/config.hpp"
 
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#ifdef MFEM_USE_CUDA
 #include <cusparse.h>
 #include <library_types.h>
 #include <cuda_runtime.h>
@@ -22,7 +22,7 @@
 #endif
 #include "cuda.hpp"
 
-#if defined(MFEM_USE_HIP) && defined(__HIP__)
+#ifdef MFEM_USE_HIP
 #include <hip/hip_runtime.h>
 #endif
 #include "hip.hpp"
@@ -43,7 +43,7 @@
 #endif
 #endif
 
-#if !defined(MFEM_USE_CUDA_OR_HIP)
+#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
 #define MFEM_DEVICE
 #define MFEM_HOST
 #define MFEM_LAMBDA
@@ -55,7 +55,7 @@
 #endif
 
 #if !((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \
-      (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
+      (defined(MFEM_USE_HIP)  && defined(__HIP_DEVICE_COMPILE__)))
 #define MFEM_SHARED
 #define MFEM_SYNC_THREAD
 #define MFEM_BLOCK_ID(k) 0
@@ -66,7 +66,7 @@
 #endif
 
 // 'double' and 'float' atomicAdd implementation for previous versions of CUDA
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && (__CUDA_ARCH__ < 600)
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600
 MFEM_DEVICE inline mfem::real_t atomicAdd(mfem::real_t *add, mfem::real_t val)
 {
    unsigned long long int *ptr = (unsigned long long int *) add;
@@ -94,7 +94,7 @@ template <typename T>
 MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
 {
 #if ((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \
-     (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
+     (defined(MFEM_USE_HIP)  && defined(__HIP_DEVICE_COMPILE__)))
    return atomicAdd(&add,val);
 #else
    T old = add;

general/communication.cpp

@@ -347,23 +347,13 @@ void GroupTopology::Swap(GroupTopology &other)
    mfem::Swap(group_mgroup, other.group_mgroup);
 }
 
-/// \cond DO_NOT_DOCUMENT
 // Initialize the static mpi_type for the specializations of MPITypeMap:
-const MPI_Datatype MPITypeMap<bool>::mpi_type = MFEM_MPI_CXX_BOOL;
-const MPI_Datatype MPITypeMap<char>::mpi_type = MPI_CHAR;
-const MPI_Datatype MPITypeMap<unsigned char>::mpi_type = MPI_UNSIGNED_CHAR;
-const MPI_Datatype MPITypeMap<short>::mpi_type = MPI_SHORT;
-const MPI_Datatype MPITypeMap<unsigned short>::mpi_type = MPI_UNSIGNED_SHORT;
 const MPI_Datatype MPITypeMap<int>::mpi_type = MPI_INT;
-const MPI_Datatype MPITypeMap<unsigned int>::mpi_type = MPI_UNSIGNED;
-const MPI_Datatype MPITypeMap<long>::mpi_type = MPI_LONG;
-const MPI_Datatype MPITypeMap<unsigned long>::mpi_type = MPI_UNSIGNED_LONG;
-const MPI_Datatype MPITypeMap<long long>::mpi_type = MPI_LONG_LONG;
-const MPI_Datatype MPITypeMap<unsigned long long>::mpi_type =
-   MPI_UNSIGNED_LONG_LONG;
 const MPI_Datatype MPITypeMap<float>::mpi_type = MPI_FLOAT;
 const MPI_Datatype MPITypeMap<double>::mpi_type = MPI_DOUBLE;
-/// \endcond DO_NOT_DOCUMENT
+const MPI_Datatype MPITypeMap<int64_t>::mpi_type = MPI_INT64_T;
+const MPI_Datatype MPITypeMap<uint64_t>::mpi_type = MPI_UINT64_T;
+
 
 GroupCommunicator::GroupCommunicator(const GroupTopology &gt, Mode m)
    : gtopo(gt), mode(m)

general/communication.hpp

@@ -607,50 +607,10 @@ protected:
 template <typename Type> struct MPITypeMap;
 
 // Specializations of MPITypeMap; mpi_type initialized in communication.cpp:
-template<> struct MPITypeMap<bool>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<char>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<unsigned char>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<short>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<unsigned short>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
 template<> struct MPITypeMap<int>
 {
    static MFEM_EXPORT const MPI_Datatype mpi_type;
 };
-template<> struct MPITypeMap<unsigned int>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<long>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<unsigned long>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<long long>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<unsigned long long>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
 template<> struct MPITypeMap<double>
 {
    static MFEM_EXPORT const MPI_Datatype mpi_type;
@@ -659,6 +619,15 @@ template<> struct MPITypeMap<float>
 {
    static MFEM_EXPORT const MPI_Datatype mpi_type;
 };
+template<> struct MPITypeMap<int64_t>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<uint64_t>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+
 
 /** Reorder MPI ranks to follow the Z-curve within the physical machine topology
     (provided that functions to query physical node coordinates are available).

general/cuda.cpp

@@ -24,7 +24,7 @@ void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
                      const char *file, int line)
 {
    mfem::err << "\n\nCUDA error: (" << expr << ") failed with error:\n --> "
-             << cudaGetErrorString(err) << " [code: " << (int)err << ']'
+             << cudaGetErrorString(err)
              << "\n ... in function: " << func
              << "\n ... in file: " << file << ':' << line << '\n';
    mfem_error();

general/cuda.hpp

@@ -18,7 +18,7 @@
 // CUDA block size used by MFEM.
 #define MFEM_CUDA_BLOCKS 256
 
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#ifdef MFEM_USE_CUDA
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,9 +37,10 @@
                               __FILE__, __LINE__);                             \
     }                                                                          \
   } while (0)
+#endif // MFEM_USE_CUDA
 
 // Define the MFEM inner threading macros
-#if defined(__CUDA_ARCH__)
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) blockIdx.k
@@ -47,13 +48,12 @@
 #define MFEM_THREAD_SIZE(k) blockDim.k
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=threadIdx.k; i<N; i+=blockDim.k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) if(const int i=threadIdx.k; i<N)
-#endif // defined(__CUDA_ARCH__)
-#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#endif
 
 namespace mfem
 {
 
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#ifdef MFEM_USE_CUDA
 // Function used by the macro MFEM_GPU_CHECK.
 void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
                      const char *file, int line);

general/error.hpp

@@ -176,7 +176,7 @@ __device__ void abort_msg(T & msg)
       printf(__VA_ARGS__);     \
       asm("trap;");            \
    }
-#elif defined(__HIP_DEVICE_COMPILE__)
+#elif defined(MFEM_USE_HIP)
 #define MFEM_ABORT_KERNEL(...) \
    {                           \
       printf(__VA_ARGS__);     \

general/forall.hpp

@@ -158,8 +158,8 @@ private:
 #define MFEM_PRAGMA(X) _Pragma(#X)
 
 // MFEM_UNROLL pragma macro that can be used inside MFEM_FORALL macros.
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) // Clang cuda or nvcc
-#ifdef __NVCC__ // nvcc specifically
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
+#ifdef __NVCC__
 #define MFEM_UNROLL(N) MFEM_PRAGMA(unroll(N))
 #else // Assuming Clang CUDA
 #define MFEM_UNROLL(N) MFEM_PRAGMA(unroll N)
@@ -169,12 +169,12 @@ private:
 #endif
 
 // MFEM_GPU_FORALL: "parallel for" executed with CUDA or HIP based on the MFEM
-// build-time configuration (MFEM_USE_CUDA or MFEM_USE_HIP), and if compiling
-// with CUDA/HIP language. Otherwise, this macro is a no-op.
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+// build-time configuration (MFEM_USE_CUDA or MFEM_USE_HIP). If neither CUDA nor
+// HIP is enabled, this macro is a no-op.
+#if defined(MFEM_USE_CUDA)
 #define MFEM_GPU_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
                                        (int i) {__VA_ARGS__})
-#elif defined(MFEM_USE_HIP) && defined(__HIP__)
+#elif defined(MFEM_USE_HIP)
 #define MFEM_GPU_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
                                         (int i) {__VA_ARGS__})
 #else
@@ -225,52 +225,6 @@ void OmpWrap(const int N, HBODY &&h_body)
 #endif
 }
 
-template <typename HBODY>
-void OmpWrap2D(const int Nx, const int Ny, HBODY &&h_body)
-{
-#ifdef MFEM_USE_OPENMP
-   // requires OpenMP 3.1
-   #pragma omp parallel for collapse(2)
-   for (int j = 0; j < Ny; j++)
-   {
-      for (int i = 0; i < Nx; i++)
-      {
-         h_body(i, j);
-      }
-   }
-#else
-   MFEM_CONTRACT_VAR(Nx);
-   MFEM_CONTRACT_VAR(Ny);
-   MFEM_CONTRACT_VAR(h_body);
-   MFEM_ABORT("OpenMP requested for MFEM but OpenMP is not enabled!");
-#endif
-}
-
-template <typename HBODY>
-void OmpWrap3D(const int Nx, const int Ny, const int Nz, HBODY &&h_body)
-{
-#ifdef MFEM_USE_OPENMP
-   // requires OpenMP 3.1
-   #pragma omp parallel for collapse(3)
-   for (int k = 0; k < Nz; k++)
-   {
-      for (int j = 0; j < Ny; j++)
-      {
-         for (int i = 0; i < Nx; i++)
-         {
-            h_body(i, j, k);
-         }
-      }
-   }
-#else
-   MFEM_CONTRACT_VAR(Nx);
-   MFEM_CONTRACT_VAR(Ny);
-   MFEM_CONTRACT_VAR(Nz);
-   MFEM_CONTRACT_VAR(h_body);
-   MFEM_ABORT("OpenMP requested for MFEM but OpenMP is not enabled!");
-#endif
-}
-
 
 /// RAJA Cuda and Hip backends
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_CUDA)
@@ -500,42 +454,6 @@ void RajaOmpWrap(const int N, HBODY &&h_body)
    RAJA::forall<RAJA::omp_parallel_for_exec>(RAJA::RangeSegment(0,N), h_body);
 }
 
-template <typename HBODY>
-void RajaOmpWrap2D(const int Nx, const int Ny, HBODY &&h_body)
-{
-   using omp_launch_policy = RAJA::LaunchPolicy<RAJA::omp_launch_t>;
-   using global_thread_xy = RAJA::LoopPolicy<RAJA::omp_for_exec>;
-   RAJA::RangeSegment xrange(0, Nx);
-   RAJA::RangeSegment yrange(0, Ny);
-   RAJA::launch<omp_launch_policy>(RAJA::ExecPlace::HOST, RAJA::LaunchParams(),
-                                   [=](RAJA::LaunchContext ctx)
-   {
-      // contiguous in x
-      RAJA::expt::loop<global_thread_xy>(ctx, xrange, yrange, [&](int i, int j)
-      {
-         h_body(i, j);
-      });
-   });
-}
-
-template <typename HBODY>
-void RajaOmpWrap3D(const int Nx, const int Ny, const int Nz, HBODY &&h_body)
-{
-   using omp_launch_policy = RAJA::LaunchPolicy<RAJA::omp_launch_t>;
-   using global_thread_xyz = RAJA::LoopPolicy<RAJA::omp_for_exec>;
-   RAJA::RangeSegment xrange(0, Nx);
-   RAJA::RangeSegment yrange(0, Ny);
-   RAJA::RangeSegment zrange(0, Nz);
-   RAJA::launch<omp_launch_policy>(RAJA::ExecPlace::HOST, RAJA::LaunchParams(),
-                                   [=](RAJA::LaunchContext ctx)
-   {
-      // contiguous in x
-      RAJA::expt::loop<global_thread_xyz>(ctx, xrange, yrange, zrange,
-                                          [&](int i, int j, int k)
-      { h_body(i, j, k); });
-   });
-}
-
 #endif
 
 
@@ -563,7 +481,7 @@ void RajaSeqWrap(const int N, HBODY &&h_body)
 
 
 /// CUDA backend
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#ifdef MFEM_USE_CUDA
 
 template <typename BODY> __global__ static
 void CuKernel1D(const int N, BODY body)
@@ -655,11 +573,11 @@ struct CuWrap<3>
    }
 };
 
-#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#endif // MFEM_USE_CUDA
 
 
 /// HIP backend
-#if defined(MFEM_USE_HIP) && defined(__HIP__)
+#ifdef MFEM_USE_HIP
 
 template <typename BODY> __global__ static
 void HipKernel1D(const int N, BODY body)
@@ -750,7 +668,7 @@ struct HipWrap<3>
    }
 };
 
-#endif // defined(MFEM_USE_HIP) && defined(__HIP__)
+#endif // MFEM_USE_HIP
 
 
 /// The forall kernel body wrapper
@@ -783,7 +701,7 @@ inline void ForallWrap(const bool use_dev, const int N,
    }
 #endif
 
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#ifdef MFEM_USE_CUDA
    // If Backend::CUDA is allowed, use it
    if (Device::Allows(Backend::CUDA))
    {
@@ -791,7 +709,7 @@ inline void ForallWrap(const bool use_dev, const int N,
    }
 #endif
 
-#if defined(MFEM_USE_HIP) && defined(__HIP__)
+#ifdef MFEM_USE_HIP
    // If Backend::HIP is allowed, use it
    if (Device::Allows(Backend::HIP))
    {
@@ -835,83 +753,6 @@ inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
 template<typename lambda>
 inline void forall(int N, lambda &&body) { ForallWrap<1>(true, N, body); }
 
-template<typename lambda>
-inline void forall(int Nx, int Ny, lambda &&body)
-{
-   if (Device::Allows(Backend::DEVICE_MASK))
-   {
-      forall(Nx * Ny, [=] MFEM_HOST_DEVICE(int idx)
-      {
-         int j = idx / Nx;
-         int i = idx % Nx;
-         body(i, j);
-      });
-   }
-#if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_OPENMP)
-   else if (Device::Allows(Backend::RAJA_OMP))
-   {
-      return RajaOmpWrap2D(Nx, Ny, body);
-   }
-#endif
-#ifdef MFEM_USE_OPENMP
-   else if (Device::Allows(Backend::OMP))
-   {
-      return OmpWrap2D(Nx, Ny, body);
-   }
-#endif
-   else
-   {
-      for (int j = 0; j < Ny; ++j)
-      {
-         for (int i = 0; i < Nx; ++i)
-         {
-            body(i, j);
-         }
-      }
-   }
-}
-
-template<typename lambda>
-inline void forall(int Nx, int Ny, int Nz, lambda &&body)
-{
-   if (Device::Allows(Backend::DEVICE_MASK))
-   {
-      forall(Nx * Ny * Nz, [=] MFEM_HOST_DEVICE(int idx)
-      {
-         int i = idx % Nx;
-         int j = idx / Nx;
-         int k = j / Ny;
-         j = j % Ny;
-         body(i, j, k);
-      });
-   }
-#if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_OPENMP)
-   else if (Device::Allows(Backend::RAJA_OMP))
-   {
-      return RajaOmpWrap3D(Nx, Ny, Nz, body);
-   }
-#endif
-#ifdef MFEM_USE_OPENMP
-   else if (Device::Allows(Backend::OMP))
-   {
-      return OmpWrap3D(Nx, Ny, Nz, body);
-   }
-#endif
-   else
-   {
-      for (int k = 0; k < Nz; ++k)
-      {
-         for (int j = 0; j < Ny; ++j)
-         {
-            for (int i = 0; i < Nx; ++i)
-            {
-               body(i, j, k);
-            }
-         }
-      }
-   }
-}
-
 template<typename lambda>
 inline void forall_switch(bool use_dev, int N, lambda &&body)
 {

general/hip.cpp

@@ -24,7 +24,7 @@ void mfem_hip_error(hipError_t err, const char *expr, const char *func,
                     const char *file, int line)
 {
    mfem::err << "\n\nHIP error: (" << expr << ") failed with error:\n --> "
-             << hipGetErrorString(err) << " [code: " << (int)err << ']'
+             << hipGetErrorString(err)
              << "\n ... in function: " << func
              << "\n ... in file: " << file << ':' << line << '\n';
    mfem_error();

general/hip.hpp

@@ -18,7 +18,7 @@
 // HIP block size used by MFEM.
 #define MFEM_HIP_BLOCKS 256
 
-#if defined(MFEM_USE_HIP) && defined(__HIP__)
+#ifdef MFEM_USE_HIP
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,9 +37,10 @@
                              __FILE__, __LINE__);                              \
     }                                                                          \
   } while (0)
+#endif // MFEM_USE_HIP
 
 // Define the MFEM inner threading macros
-#if defined(__HIP_DEVICE_COMPILE__)
+#if defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) hipBlockIdx_ ##k
@@ -49,8 +50,7 @@
    for(int i=hipThreadIdx_ ##k; i<N; i+=hipBlockDim_ ##k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) \
    if(const int i=hipThreadIdx_ ##k; i<N)
-#endif // defined(__HIP_DEVICE_COMPILE__)
-#endif // defined(MFEM_USE_HIP) && defined(__HIP__)
+#endif
 
 namespace mfem
 {

general/mem_manager.cpp

@@ -513,7 +513,10 @@ public:
    void *HtoD(void *dst, const void *src, size_t bytes) override
    { return HipMemcpyHtoD(dst, src, bytes); }
    void *DtoD(void* dst, const void* src, size_t bytes) override
-   { return HipMemcpyDtoD(dst, src, bytes); }
+   // Unlike cudaMemcpy(DtoD), hipMemcpy(DtoD) causes a host-side synchronization so
+   // instead we use hipMemcpyAsync to get similar behavior.
+   // for more info see: https://github.com/mfem/mfem/pull/2780
+   { return HipMemcpyDtoDAsync(dst, src, bytes); }
    void *DtoH(void *dst, const void *src, size_t bytes) override
    { return HipMemcpyDtoH(dst, src, bytes); }
 };
@@ -655,7 +658,10 @@ public:
       return CuMemcpyDtoD(dst, src, bytes);
 #endif
 #ifdef MFEM_USE_HIP
-      return HipMemcpyDtoD(dst, src, bytes);
+      // Unlike cudaMemcpy(DtoD), hipMemcpy(DtoD) causes a host-side synchronization so
+      // instead we use hipMemcpyAsync to get similar behavior.
+      // for more info see: https://github.com/mfem/mfem/pull/2780
+      return HipMemcpyDtoDAsync(dst, src, bytes);
 #endif
       // rm.copy(dst, const_cast<void*>(src), bytes); return dst;
    }

general/mem_manager.hpp

@@ -896,7 +896,6 @@ inline HYPRE_MemoryLocation GetHypreMemoryLocation()
 #elif MFEM_HYPRE_VERSION < 23100
    return HYPRE_MEMORY_DEVICE;
 #else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
-   if (!HYPRE_Initialized()) { return HYPRE_MEMORY_HOST; }
    HYPRE_MemoryLocation loc;
    HYPRE_GetMemoryLocation(&loc);
    return loc;
@@ -1058,8 +1057,7 @@ inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
          // register the 'base' if the MemoryManager::Exists():
          MemoryManager::Exists()
 #else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
-         IsDeviceMemory(MemoryManager::GetDeviceMemoryType()) ||
-         (MemoryManager::Exists() && HypreUsingGPU())
+         MemoryManager::Exists() && HypreUsingGPU()
 #endif
       )
       {

general/reducers.hpp

@@ -537,7 +537,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
       return;
    }
 
-#if defined(MFEM_USE_CUDA_OR_HIP)
+#if defined(MFEM_USE_HIP) || defined(MFEM_USE_CUDA)
    if (use_dev &&
        mfem::Device::Allows(Backend::CUDA | Backend::HIP | Backend::RAJA_CUDA |
                             Backend::RAJA_HIP))

general/scan.hpp

@@ -1,176 +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_SCAN_HPP
-#define MFEM_SCAN_HPP
-
-#ifdef MFEM_USE_CUDA
-#include <cub/device/device_scan.cuh>
-#define MFEM_CUB_NAMESPACE cub
-#elif MFEM_USE_HIP
-#include <hipcub/device/device_scan.hpp>
-#define MFEM_CUB_NAMESPACE hipcub
-#endif
-
-#include <functional>
-#include <numeric>
-
-namespace mfem
-{
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
-/// std::plus<>{})
-template <class InputIt, class OutputIt>
-void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   Array<char> &workspace)
-{
-   // forward to InclusiveSum for potentially faster kernels
-#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
-   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
-   {
-      size_t bytes = workspace.Size();
-      if (bytes)
-      {
-         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       workspace.Write(), bytes, d_in, d_out, num_items);
-#if defined(MFEM_USE_CUDA)
-         if (err == cudaSuccess)
-         {
-            return;
-         }
-#elif defined(MFEM_USE_HIP)
-         if (err == hipSuccess)
-         {
-            return;
-         }
-#endif
-      }
-      // try allocating a larger buffer
-      bytes = 0;
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        nullptr, bytes, d_in, d_out, num_items));
-      workspace.SetSize(bytes);
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        workspace.Write(), bytes, d_in, d_out, num_items));
-      return;
-   }
-#endif
-   std::inclusive_scan(d_in, d_in + num_items, d_out);
-}
-
-/// Performs an inclusive scan of [d_in, d_in+num_items) -> [d_out,
-/// d_out+num_items). This call is potentially asynchronous on the device.
-/// @a d_in input start.
-/// @a d_out output start. Can perform in-place scans with d_out = d_in
-/// @a workspace temporary workspace used for device scans. TODO: replace with
-/// internal temporary workspace once that's added to the memory manager.
-/// @a scan_op binary scan functor. Must be associative. If only weakly
-/// associative (i.e. floating point addition) results are not deterministic. On
-/// device this must also be commutative.
-template <class InputIt, class OutputIt, class ScanOp>
-void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   Array<char> &workspace, ScanOp scan_op)
-{
-#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
-   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
-   {
-      size_t bytes = workspace.Size();
-      if (bytes)
-      {
-         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
-#if defined(MFEM_USE_CUDA)
-         if (err == cudaSuccess)
-         {
-            return;
-         }
-#elif defined(MFEM_USE_HIP)
-         if (err == hipSuccess)
-         {
-            return;
-         }
-#endif
-      }
-      // try allocating a larger buffer
-      bytes = 0;
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        nullptr, bytes, d_in, d_out, scan_op, num_items));
-      workspace.SetSize(bytes);
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
-      return;
-   }
-#endif
-   std::inclusive_scan(d_in, d_in + num_items, d_out, scan_op);
-}
-
-/// Performs an exclusive scan of [d_in, d_in+num_items) -> [d_out,
-/// d_out+num_items). This call is potentially asynchronous on the device.
-/// @a d_in input start.
-/// @a d_out output start. Can perform in-place scans with d_out = d_in
-/// @a workspace temporary workspace used for device scans. TODO: replace with
-/// internal temporary workspace once that's added to the memory manager.
-/// @a scan_op binary scan functor. Must be associative. If only weakly
-/// associative (i.e. floating point addition) results are not deterministic. On
-/// device this must also be commutative.
-template <class InputIt, class OutputIt, class T, class ScanOp>
-void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   T init_value, Array<char> &workspace, ScanOp scan_op)
-{
-#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
-   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
-   {
-      size_t bytes = workspace.Size();
-      if (bytes)
-      {
-         auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                       num_items);
-#if defined(MFEM_USE_CUDA)
-         if (err == cudaSuccess)
-         {
-            return;
-         }
-#elif defined(MFEM_USE_HIP)
-         if (err == hipSuccess)
-         {
-            return;
-         }
-#endif
-      }
-      // try allocating a larger buffer
-      bytes = 0;
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
-      workspace.SetSize(bytes);
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                        num_items));
-      return;
-   }
-#endif
-   std::exclusive_scan(d_in, d_in + num_items, d_out, init_value, scan_op);
-}
-
-/// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// workspace, std::plus<>{})
-template <class InputIt, class OutputIt, class T>
-void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   T init_value, Array<char> &workspace)
-{
-   ExclusiveScan(use_dev, d_in, d_out, num_items, init_value, workspace,
-                 std::plus<> {});
-}
-
-} // namespace mfem
-
-#undef MFEM_CUB_NAMESPACE
-
-#endif

linalg/vector.cpp

@@ -14,7 +14,6 @@
 #include "../general/forall.hpp"
 #include "../general/reducers.hpp"
 #include "../general/hash.hpp"
-#include "../general/scan.hpp"
 #include "vector.hpp"
 
 #ifdef MFEM_USE_OPENMP
@@ -118,7 +117,7 @@ Vector::Vector(const Vector &v)
    UseDevice(v.UseDevice());
 }
 
-Vector::Vector(Vector &&v) : Vector()
+Vector::Vector(Vector &&v)
 {
    *this = std::move(v);
 }
@@ -1253,44 +1252,4 @@ real_t Vector::Sum() const
    return res;
 }
 
-void Vector::DeleteAt(const Array<int> &indices)
-{
-   const bool use_dev = UseDevice();
-
-   Array<int> flag(size);
-   const auto d_flag = flag.Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
-   {
-      d_flag[i] = true;
-   });
-   const auto d_indices = indices.Read(use_dev);
-   mfem::forall_switch(use_dev, indices.Size(), [=] MFEM_HOST_DEVICE (int i)
-   {
-      d_flag[d_indices[i]] = false;
-   });
-
-   Array<int> out_idx(size);
-   auto d_out_idx = out_idx.Write(use_dev);
-   Array<char> workspace;
-   // Perform inclusive scan so that the last entry is the new size.
-   InclusiveScan(use_dev, d_flag, d_out_idx, size, workspace);
-
-   Vector copy(*this);
-   auto d_in = copy.Read(use_dev);
-   auto d_out = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
-   {
-      if (d_flag[i])
-      {
-         // Transform inclusive scan to exclusive by shifting.
-         const int j = (i > 0) ? d_out_idx[i - 1] : 0;
-         d_out[j] = d_in[i];
-      }
-   });
-
-   // Get the new size of the vector. Copy only the last entry.
-   Memory<int> submem(out_idx.GetMemory(), out_idx.Size() - 1, 1);
-   size = submem.Read(MemoryClass::HOST, 1)[0];
-}
-
 } // namespace mfem

linalg/vector.hpp

@@ -171,12 +171,6 @@ public:
    /// Resize the vector to size @a s using the MemoryType of @a v.
    void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
 
-   /// Update \ref Capacity() to @a res (if less than current), keeping existing entries.
-   void Reserve(int res);
-
-   /// Delete entries at @a indices and resize vector accordingly.
-   void DeleteAt(const Array<int> &indices);
-
    /// Set the Vector data.
    /// @warning This method should be called only when OwnsData() is false.
    void SetData(real_t *d) { data.Wrap(d, data.Capacity(), false); }
@@ -627,18 +621,6 @@ inline void Vector::SetSize(int s, MemoryType mt)
    data.UseDevice(use_dev);
 }
 
-inline void Vector::Reserve(int res)
-{
-   if (res > Capacity())
-   {
-      Memory<real_t> p(res, data.GetMemoryType());
-      p.CopyFrom(data, size);
-      p.UseDevice(data.UseDevice());
-      data.Delete();
-      data = p;
-   }
-}
-
 inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
                                      bool own_mem)
 {

mesh/mesh.cpp

@@ -32,7 +32,6 @@
 #include <cstring>
 #include <ctime>
 #include <functional>
-#include <set>
 #include <numeric>
 #include <unordered_map>
 #include <unordered_set>
@@ -922,103 +921,6 @@ const FaceGeometricFactors* Mesh::GetFaceGeometricFactors(
    return gf;
 }
 
-const Array<int>& Mesh::GetBdrFaceAttributes() const
-{
-   if (bdr_face_attrs_cache.Size() == 0)
-   {
-      std::unordered_map<int, int> f_to_be;
-      for (int i = 0; i < GetNBE(); ++i)
-      {
-         const int f = GetBdrElementFaceIndex(i);
-         f_to_be[f] = i;
-      }
-      const int nf_bdr = GetNFbyType(FaceType::Boundary);
-      // MFEM_VERIFY(size_t(nf_bdr) == f_to_be.size(), "Incompatible sizes");
-      bdr_face_attrs_cache.SetSize(nf_bdr);
-      int f_ind = 0;
-      const int nf = GetNumFaces();
-      for (int f = 0; f < nf; ++f)
-      {
-         if (!GetFaceInformation(f).IsOfFaceType(FaceType::Boundary))
-         {
-            continue;
-         }
-         int attribute = -1; // default value
-         auto iter = f_to_be.find(f);
-         if (iter != f_to_be.end())
-         {
-            const int be = iter->second;
-            attribute = GetBdrAttribute(be);
-         }
-         else
-         {
-            // If a boundary face does not correspond to the a boundary element,
-            // we assign it the default attribute of -1.
-         }
-         bdr_face_attrs_cache[f_ind] = attribute;
-         ++f_ind;
-      }
-   }
-   return bdr_face_attrs_cache;
-}
-
-const Array<int>& Mesh::GetElementAttributes() const
-{
-   if (elem_attrs_cache.Size() == 0)
-   {
-      // re-compute cache
-      elem_attrs_cache.SetSize(GetNE());
-      elem_attrs_cache.HostWrite();
-      for (int i = 0; i < GetNE(); ++i)
-      {
-         elem_attrs_cache[i] = GetAttribute(i);
-         MFEM_ASSERT(elem_attrs_cache[i] > 0,
-                     "Negative attribute on element " << i);
-      }
-   }
-   return elem_attrs_cache;
-}
-
-void Mesh::ComputeFaceInfo(FaceType ftype) const
-{
-   auto &fidcs = face_indices[static_cast<int>(ftype)];
-   auto &ifidcs = inv_face_indices[static_cast<int>(ftype)];
-   fidcs.SetSize(GetNFbyType(ftype));
-   fidcs.HostWrite();
-   ifidcs.reserve(fidcs.Size());
-   int f_idx = 0;
-   for (int i = 0; i < GetNumFacesWithGhost(); ++i)
-   {
-      const FaceInformation face = GetFaceInformation(i);
-      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(ftype))
-      {
-         continue;
-      }
-      fidcs[f_idx] = i;
-      ifidcs[i] = f_idx;
-      ++f_idx;
-   }
-}
-
-const Array<int> &Mesh::GetFaceIndices(FaceType ftype) const
-{
-   if (face_indices[static_cast<int>(ftype)].Size() == 0)
-   {
-      ComputeFaceInfo(ftype);
-   }
-   return face_indices[static_cast<int>(ftype)];
-}
-
-const std::unordered_map<int, int> &
-Mesh::GetInvFaceIndices(FaceType ftype) const
-{
-   if (inv_face_indices[static_cast<int>(ftype)].empty())
-   {
-      ComputeFaceInfo(ftype);
-   }
-   return inv_face_indices[static_cast<int>(ftype)];
-}
-
 void Mesh::DeleteGeometricFactors()
 {
    for (int i = 0; i < geom_factors.Size(); i++)
@@ -1902,15 +1804,8 @@ void Mesh::Destroy()
    TetMemory.Clear();
 #endif
 
-   elem_attrs_cache.DeleteAll();
-   bdr_face_attrs_cache.DeleteAll();
    attributes.DeleteAll();
    bdr_attributes.DeleteAll();
-
-   face_indices[0].DeleteAll();
-   face_indices[1].DeleteAll();
-   inv_face_indices[0] = std::unordered_map<int, int>();
-   inv_face_indices[1] = std::unordered_map<int, int>();
 }
 
 void Mesh::ResetLazyData()
@@ -1922,48 +1817,36 @@ void Mesh::ResetLazyData()
    DeleteGeometricFactors();
    nbInteriorFaces = -1;
    nbBoundaryFaces = -1;
-   // set size to 0 so re-computations can potentially avoid a new allocation
-   bdr_face_attrs_cache.SetSize(0);
-   elem_attrs_cache.SetSize(0);
 }
 
-void Mesh::SetAttributes(bool elem_attrs_changed, bool bdr_face_attrs_changed)
+void Mesh::SetAttributes()
 {
-   if (bdr_face_attrs_changed)
+   Array<int> attribs;
+
+   attribs.SetSize(GetNBE());
+   for (int i = 0; i < attribs.Size(); i++)
    {
-      bdr_face_attrs_cache.SetSize(0); // Invalidate the cache
-
-      // Get sorted list of unique boundary element attributes
-      std::set<int> attribs;
-      for (int i = 0; i < GetNBE(); i++)
-      {
-         attribs.emplace(GetBdrAttribute(i));
-      }
-
-      bdr_attributes.SetSize(attribs.size());
-      bdr_attributes.HostWrite();
-      std::copy(attribs.begin(), attribs.end(), bdr_attributes.begin());
-      if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
-      {
-         MFEM_WARNING("Non-positive attributes on the boundary!");
-      }
+      attribs[i] = GetBdrAttribute(i);
+   }
+   attribs.Sort();
+   attribs.Unique();
+   attribs.Copy(bdr_attributes);
+   if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
+   {
+      MFEM_WARNING("Non-positive attributes on the boundary!");
    }
 
-   if (elem_attrs_changed)
+   attribs.SetSize(GetNE());
+   for (int i = 0; i < attribs.Size(); i++)
    {
-      // Re-compute the attributes cache
-      elem_attrs_cache.SetSize(0);
-      GetElementAttributes();
-      // Get sorted list of unique element attributes
-      std::set<int> attribs(elem_attrs_cache.begin(), elem_attrs_cache.end());
-      attributes.SetSize(attribs.size());
-      attributes.HostWrite();
-      std::copy(attribs.begin(), attribs.end(), attributes.begin());
-
-      if (attributes.Size() > 0 && attributes[0] <= 0)
-      {
-         MFEM_WARNING("Non-positive attributes in the domain!");
-      }
+      attribs[i] = GetAttribute(i);
+   }
+   attribs.Sort();
+   attribs.Unique();
+   attribs.Copy(attributes);
+   if (attributes.Size() > 0 && attributes[0] <= 0)
+   {
+      MFEM_WARNING("Non-positive attributes in the domain!");
    }
 }
 
@@ -4571,10 +4454,6 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
       Nodes = mesh.Nodes;
       own_nodes = 0;
    }
-
-   // copy attribute caches
-   elem_attrs_cache = mesh.elem_attrs_cache;
-   bdr_face_attrs_cache = mesh.bdr_face_attrs_cache;
 }
 
 Mesh::Mesh(Mesh &&mesh) : Mesh()
@@ -7842,12 +7721,6 @@ void Mesh::GetBdrElementAdjacentElement2(
 void Mesh::SetAttribute(int i, int attr)
 {
   elements[i]->SetAttribute(attr);
-  if (elem_attrs_cache.Size() == GetNE())
-  {
-     // update the existing cache instead of deleting it
-     elem_attrs_cache.HostReadWrite();
-     elem_attrs_cache[i] = attr;
-  }
   if (ncmesh) ncmesh->SetAttribute(i, attr);
 }
 
@@ -8182,12 +8055,6 @@ void Mesh::GenerateFaces()
       FreeElement(f);
    }
 
-   // delete caches
-   face_indices[0].SetSize(0);
-   face_indices[1].SetSize(0);
-   inv_face_indices[0].clear();
-   inv_face_indices[1].clear();
-
    // (re)generate the interior faces and the info for them
    faces.SetSize(nfaces);
    faces_info.SetSize(nfaces);
@@ -10983,15 +10850,6 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
       mfem::Swap(nodes_sequence, other.nodes_sequence);
       mfem::Swap(last_operation, other.last_operation);
    }
-
-   // copy attribute caches
-   mfem::Swap(elem_attrs_cache, other.elem_attrs_cache);
-   mfem::Swap(bdr_face_attrs_cache, other.bdr_face_attrs_cache);
-
-   mfem::Swap(face_indices[0], other.face_indices[0]);
-   mfem::Swap(face_indices[1], other.face_indices[1]);
-   inv_face_indices[0].swap(other.inv_face_indices[0]);
-   inv_face_indices[1].swap(other.inv_face_indices[1]);
 }
 
 void Mesh::GetElementData(const Array<Element*> &elem_array, int geom,

mesh/mesh.hpp

@@ -107,11 +107,6 @@ protected:
    Array<Element *> boundary;
    Array<Element *> faces;
 
-   /// internal cache for element attributes
-   mutable Array<int> elem_attrs_cache;
-   /// internal cache for boundary element attributes
-   mutable Array<int> bdr_face_attrs_cache;
-
    /** @brief This structure stores the low level information necessary to
        interpret the configuration of elements on a specific face. This
        information can be accessed using methods like GetFaceElements(),
@@ -278,13 +273,6 @@ protected:
 
    // used during NC mesh initialization only
    Array<Triple<int, int, int> > tmp_vertex_parents;
-   /// cache for FaceIndices(ftype)
-   mutable Array<int> face_indices[2];
-   /// cache for FaceIndices(ftype)
-   mutable std::unordered_map<int, int> inv_face_indices[2];
-
-   /// compute face_indices[ftype] and inv_face_indices[type]
-   void ComputeFaceInfo(FaceType ftype) const;
 
 public:
    typedef Geometry::Constants<Geometry::SEGMENT>     seg_t;
@@ -319,11 +307,6 @@ public:
    // (true) is set in mesh_readers.cpp.
    static bool remove_unused_vertices;
 
-   /// Map from boundary or interior face indices to mesh face indices.
-   const Array<int>& GetFaceIndices(FaceType ftype) const;
-   /// Inverse of the map FaceIndices(ftype)
-   const std::unordered_map<int, int>& GetInvFaceIndices(FaceType ftype) const;
-
 protected:
    Operation last_operation;
 
@@ -1139,14 +1122,13 @@ public:
        Mesh vertices or nodes are set. */
    virtual void Finalize(bool refine = false, bool fix_orientation = false);
 
-   /// @brief Determine the sets of unique attribute values in domain if @a
-   /// elem_attrs_changed and boundary elements if @a bdr_face_attrs_changed.
+   /// @brief Determine the sets of unique attribute values in domain and
+   /// boundary elements.
    ///
    /// Separately scan the domain and boundary elements to generate unique,
    /// sorted sets of the element attribute values present in the mesh and
    /// store these in the Mesh::attributes and Mesh::bdr_attributes arrays.
-   virtual void SetAttributes(bool elem_attrs_changed = true,
-                              bool bdr_face_attrs_changed = true);
+   virtual void SetAttributes();
 
    /// Check (and optionally attempt to fix) the orientation of the elements
    /** @param[in] fix_it  If `true`, attempt to fix the orientations of some
@@ -2285,35 +2267,6 @@ public:
        method does not modify the nodes. */
    void NodesUpdated() { DeleteGeometricFactors(); }
 
-   /// @brief Returns the attributes for all elements in this mesh. The i'th
-   /// entry of the array is the attribute of the i'th element of the mesh.
-   ///
-   /// The returned array points to an internal object that may be invalidated
-   /// by mesh operations such as refinement or any element attributes are
-   /// modified. Since not all such modifications can be tracked by the Mesh
-   /// class (e.g. if a user calls GetElement() then changes the element
-   /// attribute directly), one needs to account for such changes by calling the
-   /// method SetAttributes().
-   const Array<int>& GetElementAttributes() const;
-
-   /// @brief Returns the attributes for all boundary elements in this mesh.
-   ///
-   /// The face restriction will give "face E-vectors" on the boundary that
-   /// are numbered in the order of the faces of mesh. This numbering will be
-   /// different than the numbering of the boundary elements. We compute
-   /// mappings so that the array `bdr_attributes[i]` gives the boundary
-   /// attribute of the `i`th boundary face in the mesh face order.
-   /// Attributes <= 0 indicate there is no boundary element and should be
-   /// skipped.
-   ///
-   /// The returned array points to an internal object that may be invalidated
-   /// by mesh operations such as refinement or any element attributes are
-   /// modified. Since not all such modifications can be tracked by the Mesh
-   /// class (e.g. if a user calls GetElement() then changes the element
-   /// attribute directly), one needs to account for such changes by calling the
-   /// method SetAttributes().
-   const Array<int>& GetBdrFaceAttributes() const;
-
    /// @}
 
    /// @anchor mfem_Mesh_gf_nodes

mesh/pmesh.cpp

@@ -1590,27 +1590,21 @@ void ParMesh::DistributeAttributes(Array<int> &attr)
    delete [] glb_attr_marker;
 }
 
-void ParMesh::SetAttributes(bool elem_attrs_changed, bool bdr_attrs_changed)
+void ParMesh::SetAttributes()
 {
    // Determine the attributes occurring in local interior and boundary elements
-   Mesh::SetAttributes(elem_attrs_changed, bdr_attrs_changed);
+   Mesh::SetAttributes();
 
-   if (bdr_attrs_changed)
+   DistributeAttributes(bdr_attributes);
+   if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
    {
-      DistributeAttributes(bdr_attributes);
-      if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
-      {
-         MFEM_WARNING("Non-positive boundary element attributes found!");
-      }
+      MFEM_WARNING("Non-positive boundary element attributes found!");
    }
 
-   if (elem_attrs_changed)
+   DistributeAttributes(attributes);
+   if (attributes.Size() > 0 && attributes[0] <= 0)
    {
-      DistributeAttributes(attributes);
-      if (attributes.Size() > 0 && attributes[0] <= 0)
-      {
-         MFEM_WARNING("Non-positive element attributes found!");
-      }
+      MFEM_WARNING("Non-positive element attributes found!");
    }
 }
 

mesh/pmesh.hpp

@@ -394,8 +394,7 @@ public:
 
    void Finalize(bool refine = false, bool fix_orientation = false) override;
 
-   void SetAttributes(bool elem_attrs_changed = true,
-                      bool bdr_attrs_changed = true) override;
+   void SetAttributes() override;
 
    /// Checks if any rank in the mesh has boundary elements
    bool HasBoundaryElements() const override;

tests/unit/CMakeLists.txt

@@ -22,7 +22,6 @@ set(UNIT_TESTS_SRCS
   dfem/test_mass.cpp
   general/test_array.cpp
   general/test_reduction.cpp
-  general/test_scan.cpp
   general/test_arrays_by_name.cpp
   general/test_error.cpp
   general/test_mem.cpp

tests/unit/general/test_array.cpp

@@ -124,19 +124,3 @@ TEST_CASE("Array stl-interactions", "[Array]")
       CHECK(x[i] == y[i]);
    }
 }
-
-TEST_CASE("Array delete at indices", "[Array]")
-{
-   Array<int>       test({0,1,2,3,4,5,6,7,8});
-   Array<int> rm_indices({0,    3,4,  6,  8});
-   Array<int>     result({  1,2,    5,  7  });
-
-   test.DeleteAt(rm_indices);
-
-   REQUIRE(test.Size() == result.Size());
-
-   for (int i = 0; i < test.Size(); i++)
-   {
-      CHECK(test[i] == result[i]);
-   }
-}

tests/unit/general/test_scan.cpp

@@ -1,102 +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.
-
-#include <algorithm>
-#include <limits>
-
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-// must be included after mfem.hpp
-#include "general/scan.hpp"
-
-using namespace mfem;
-
-TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
-{
-   Array<char> workspace;
-   Array<int> a(10);
-
-   for (int use_dev = 0; use_dev < 2; ++use_dev)
-   {
-      CAPTURE(use_dev);
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i;
-      }
-      auto dptr = a.ReadWrite(use_dev);
-      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace);
-      a.HostRead();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         int expected = (i + 1) * i / 2;
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-      }
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i + 1;
-      }
-      a.ReadWrite(use_dev);
-      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace, std::multiplies<> {});
-      a.HostRead();
-      int expected = 1;
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         expected *= i + 1;
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-      }
-   }
-}
-
-TEST_CASE("Exclusive Scan", "[Scan],[GPU]")
-{
-   Array<char> workspace;
-   Array<int> a(10);
-
-   for (int use_dev = 0; use_dev < 2; ++use_dev)
-   {
-      CAPTURE(use_dev);
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i;
-      }
-      auto dptr = a.ReadWrite(use_dev);
-      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace);
-      a.HostRead();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         int expected = (i + 1) * i / 2 - i + 5;
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-      }
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i + 1;
-      }
-      a.ReadWrite(use_dev);
-      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace,
-                    std::multiplies<> {});
-      a.HostRead();
-      int expected = 5;
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-         expected *= i + 1;
-      }
-   }
-}

tests/unit/linalg/test_vector.cpp

@@ -247,21 +247,3 @@ TEST_CASE("Vector Sum", "[Vector],[GPU]")
 
    REQUIRE(sum_1 == MFEM_Approx(sum_2));
 }
-
-TEST_CASE("Vector delete at indices", "[Vector][GPU]")
-{
-   Vector           test({0,1,2,3,4,5,6,7,8});
-   Array<int> rm_indices({0,    3,4,  6,  8});
-   Vector         result({  1,2,    5,  7  });
-
-   test.UseDevice(true);
-   test.DeleteAt(rm_indices);
-
-   REQUIRE(test.Size() == result.Size());
-
-   test.HostReadWrite();
-   for (int i = 0; i < test.Size(); i++)
-   {
-      CHECK(test[i] == result[i]);
-   }
-}