Add DG vector diffusion example

Merge pull request #4937 from mfem/fec-new-fix
bugfix for 4935
2025-09-04 21:57:02 -07:00 · 2025-09-04 12:19:26 -07:00 · 2025-08-29 08:16:41 -07:00 · 2025-08-26 12:16:19 -07:00 · 2025-08-19 16:54:27 -07:00 · 2025-08-19 16:32:39 -07:00
48 changed files with 1257 additions and 370 deletions
@@ -132,12 +132,14 @@ jobs:
            hypre-target: int32
            precision: fp64
            enzyme: true
-            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme)
+            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme) LDFLAGS=-L$LLVM_PREFIX/lib/c++

    name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}

    runs-on: ${{ matrix.os }}

+    continue-on-error: ${{ matrix.enzyme && true || false }}
+
    steps:
      # Fix 'No space left on device' errors for Ubuntu builds.
      - name: Run Actions Cleaner
@@ -168,10 +170,13 @@ 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_15.3.app"
+          XCODE_PATH="/Applications/Xcode_16.4.app"
          echo "> sudo xcode-select -s ${XCODE_PATH}"
          sudo xcode-select -s ${XCODE_PATH}
          echo "> g++ -v"
@@ -289,10 +294,12 @@ jobs:
        run: |
          export HOMEBREW_NO_INSTALL_CLEANUP=1
          brew update
-          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
+          brew install enzyme
+          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
+          LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
+          echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
+          echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV

      # MFEM build and test
      - name: build
@@ -63,6 +63,7 @@ examples/ex3[0-9]
 examples/ex3[0-9]p
 examples/ex4[0-9]
 examples/ex4[0-9]p
+examples/vector-dg-diffusion

 examples/refined.mesh
 examples/displaced.mesh
@@ -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 ruby resource are allocated/released once for all.
+# - Allocate/Release is where Dane 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:
-ruby-build-and-test:
+dane-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 @@ ruby-build-and-test:
    AUTOTEST: "${AUTOTEST}"
    AUTOTEST_COMMIT: "${AUTOTEST_COMMIT}"
  trigger:
-    include: .gitlab/ruby-build-and-test.yml
+    include: .gitlab/dane-build-and-test.yml
    strategy: depend

-ruby-baseline:
+dane-baseline:
  stage: sub-pipelines
  variables:
    # Explicitly pass down values that we want to be able to set when triggering
@@ -73,7 +73,7 @@ ruby-baseline:
    AUTOTEST: "${AUTOTEST}"
    AUTOTEST_COMMIT: "${AUTOTEST_COMMIT}"
  trigger:
-    include: .gitlab/ruby-baseline.yml
+    include: .gitlab/dane-baseline.yml
    strategy: depend

 lassen-build-and-test:
@@ -24,7 +24,7 @@ and `test type`.

 Machines typically include:

-* Ruby: 2nd Gen Intel Xeon (Cascade Lake)
+* Dane: Intel Sapphire Rapids
 * 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 ruby for example resumes to:
+spack spec to use. Adding a job on Dane for example resumes to:

 ```yaml
 <job_name>:
  variables:
    SPEC: "<spack_spec>"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane
 ```

 The remaining and non trivial work is to make sure this spec is working. To
@@ -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 ruby, there is only one allocation shared among jobs in order to
+# On LLNL's Dane, 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}
@@ -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 Ruby machine at LLNL
+# GitLab pipelines configurations for the Dane machine at LLNL
 variables:
-  MACHINE_NAME: ruby
+  MACHINE_NAME: dane

-.on_ruby:
+.on_dane:
  tags:
    - shell
-    - ruby
+    - dane
  rules:
-    # Don't run ruby jobs if...
-    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_RUBY == "OFF"'
+    # Don't run dane jobs if...
+    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_DANE == "OFF"'
      when: never
    # Don't run autotest update if...
    - if: '$CI_JOB_NAME =~ /report/ && $AUTOTEST != "YES"'
@@ -40,13 +40,14 @@ variables:
    - when: on_success

 # Spack helped builds
-# Generic ruby build job, extending build script
-.build_and_test_on_ruby:
-  extends: [.on_ruby]
+# Generic dane build job, extending build script
+.build_and_test_on_dane:
+  extends: [.on_dane]
  stage: build_and_test
  script:
    # THREADS is used by 'tests/gitlab/build_and_test', run below
-    - export THREADS=16
+    # Dane has 224 threads/node and we run 7 separate jobs: 224=7*32
+    - export THREADS=28
    - echo ${ALLOC_NAME}
    - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
    - echo ${JOBID}
@@ -18,7 +18,7 @@
 setup_baseline:
  tags:
    - shell
-    - ruby
+    - dane
  stage: setup
  variables:
    GIT_STRATEGY: none
@@ -16,7 +16,7 @@
 setup:
  tags:
    - shell
-    - ruby
+    - dane
  stage: setup
  variables:
    GIT_STRATEGY: none
@@ -19,8 +19,8 @@ stages:
  - cleanup
  - baseline_publish

-baselinecheck_mfem_intel_ruby:
-  extends: [.on_ruby]
+baselinecheck_mfem_intel_dane:
+  extends: [.on_dane]
  stage: baseline_check
  variables:
    # TPLS_DIR is used in .gitlab/scripts/baseline to provide the tpls location
@@ -31,8 +31,8 @@ baselinecheck_mfem_intel_ruby:
  script:
    - echo ${BUILD_ROOT}
    - echo ${TPLS_DIR}
-    # Used by the tests in MFEM/tests:
-    - export MFEM_TEST_NP=48
+    # Used by the tests in MFEM/tests, dane has 224 threads/node:
+    - export MFEM_TEST_NP=192
    # 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_ruby:
  allow_failure: true

 cleanup:
-  extends: .on_ruby
+  extends: .on_dane
  stage: cleanup
  variables:
    GIT_STRATEGY: none
@@ -53,7 +53,7 @@ cleanup:
    - rm -rf "${BUILD_ROOT}" || true

 report_baseline:
-  extends: [.on_ruby]
+  extends: [.on_dane]
  stage: baseline_report
  script:
    - echo ${MACHINE_NAME}
@@ -113,8 +113,8 @@ report_baseline:
        exit $err
      ) 9> autotest.lock

-baselinepublish_mfem_ruby:
-  extends: [.on_ruby]
+baselinepublish_mfem_dane:
+  extends: [.on_dane]
  stage: baseline_publish
  rules:
    # - if: '$CI_COMMIT_BRANCH == "master" || $REBASELINE == "YES"'
@@ -129,5 +129,5 @@ baselinepublish_mfem_ruby:

 include:
  - local: .gitlab/configs/common.yml
-  - local: .gitlab/configs/ruby-config.yml
+  - local: .gitlab/configs/dane-config.yml
  - local: .gitlab/configs/setup-baseline.yml
@@ -19,54 +19,54 @@ stages:
 allocate_resource:
  variables:
    GIT_STRATEGY: none
-  extends: .on_ruby
+  extends: .on_dane
  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 Ruby machine at LLNL
+# GitLab jobs for the Dane machine at LLNL
 debug_ser_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1 +debug~mpi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 debug_par_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1 +debug+mpi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_ser_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1 ~mpi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10_sundials:
  variables:
    SPEC: "%gcc@10.3.1 +sundials"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10_petsc:
  variables:
    SPEC: "%gcc@10.3.1 +petsc ^petsc+mumps~superlu-dist"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10_pumi:
  variables:
    SPEC: "%gcc@10.3.1 +pumi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 # Release
 release_resource:
  variables:
    GIT_STRATEGY: none
-  extends: .on_ruby
+  extends: .on_dane
  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_ruby
+    - .on_dane
    - .report_job_success

 report_job_failure:
  stage: release_resource_and_report
  extends:
-    - .on_ruby
+    - .on_dane
    - .report_job_failure

 include:
  - local: .gitlab/configs/common.yml
-  - local: .gitlab/configs/ruby-config.yml
+  - local: .gitlab/configs/dane-config.yml
  - local: .gitlab/configs/setup-build-and-test.yml
  - local: .gitlab/configs/report-build-and-test.yml
@@ -14,7 +14,7 @@
 # locals
 glob_err=${BASELINE_TEST}.err
 base=${BASELINE_TEST}-${SYS_TYPE}
-if [[ "${MACHINE_NAME}" == "ruby" ]]; then
+if [[ "${MACHINE_NAME}" == "dane" ]]; 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}" == "ruby" ]]; then
+if [[ "${MACHINE_NAME}" == "dane" ]]; 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}"
@@ -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 ruby baselines.
+# There will be collision between corona and dane 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 ruby baselines.
+# There will be collision between corona and dane baselines.
 # Once the corresponding files have been generated, we can switch to machine
 # specific ref.
 SAVED_NAME=baseline-${SYS_TYPE}.saved
@@ -668,6 +668,7 @@ 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).
@@ -847,6 +848,7 @@ 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
@@ -23,11 +23,14 @@
 #include "_config.hpp"
 #endif

+#include <cstdint>
+#include <climits>
+
 namespace mfem
 {

 #if (defined(MFEM_USE_CUDA) && defined(__CUDACC__)) || \
-    (defined(MFEM_USE_HIP) && defined(__HIPCC__))
+    (defined(MFEM_USE_HIP) && defined(__HIP__))
 #define MFEM_HOST_DEVICE __host__ __device__
 #else
 #define MFEM_HOST_DEVICE
@@ -0,0 +1,164 @@
+#include "mfem.hpp"
+#include "vector-dg-diffusion.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class RepeatedCoefficient : public VectorCoefficient
+{
+   Coefficient &coeff;
+public:
+   RepeatedCoefficient(int dim, Coefficient &coeff_)
+      : VectorCoefficient(dim), coeff(coeff_)
+   { }
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      V.SetSize(vdim);
+      V = coeff.Eval(T, ip);
+   }
+};
+
+real_t u_fn(const Vector &xvec);
+real_t f_fn(const Vector &xvec);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ref_levels = 0;
+   int order = 1;
+   real_t sigma = -1.0;
+   real_t kappa = -1.0;
+   const char *device_config = "cpu";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the three DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the three DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.ParseCheck();
+
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+
+   Device device(device_config);
+   device.Print();
+
+   Mesh mesh(mesh_file);
+   const int dim = mesh.Dimension();
+
+   {
+      if (ref_levels < 0)
+      {
+         ref_levels = (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
+      }
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec, dim);
+   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
+
+
+   FunctionCoefficient scalar_f_coeff(f_fn);
+   FunctionCoefficient scalar_u_coeff(u_fn);
+   RepeatedCoefficient f_coeff(dim, scalar_f_coeff);
+   RepeatedCoefficient u_coeff(dim, scalar_u_coeff);
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(5.0);
+   RepeatedCoefficient zero_vec(dim, zero);
+
+   LinearForm b(&fespace);
+   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f_coeff));
+   b.AddBdrFaceIntegrator(
+      new VectorDGDirichletLFIntegrator(u_coeff, one, sigma, kappa));
+   b.Assemble();
+
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new VectorDiffusionIntegrator(one));
+   a.AddInteriorFaceIntegrator(new VectorDGDiffusionIntegrator(
+                                  one, sigma, kappa, dim));
+   a.AddBdrFaceIntegrator(new VectorDGDiffusionIntegrator(
+                             one, sigma, kappa, dim));
+   a.Assemble();
+   a.Finalize();
+
+   const SparseMatrix &A = a.SpMat();
+#ifndef MFEM_USE_SUITESPARSE
+   GSSmoother M(A);
+   if (sigma == -1.0)
+   {
+      PCG(A, M, b, x, 1, 500, 1e-12, 0.0);
+   }
+   else
+   {
+      GMRES(A, M, b, x, 1, 500, 10, 1e-12, 0.0);
+   }
+#else
+   UMFPackSolver umf_solver;
+   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   umf_solver.SetOperator(A);
+   umf_solver.Mult(b, x);
+#endif
+
+   ParaViewDataCollection pv("DGDiffusion", &mesh);
+   pv.SetPrefixPath("ParaView");
+   pv.SetHighOrderOutput(true);
+   pv.SetLevelsOfDetail(order);
+   pv.RegisterField("u", &x);
+   pv.SetCycle(0);
+   pv.SetTime(0.0);
+   pv.Save();
+
+   cout << "L2 error: " << x.ComputeL2Error(u_coeff) << '\n';
+
+   return 0;
+}
+
+constexpr real_t pi = M_PI;
+constexpr real_t pi2 = pi*pi;
+
+real_t u_fn(const Vector &xvec)
+{
+   int dim = xvec.Size();
+   real_t x = pi*xvec[0], y = pi*xvec[1];
+   if (dim == 2) { return sin(x)*sin(y); }
+   else { real_t z = pi*xvec[2]; return sin(x)*sin(y)*sin(z); }
+}
+
+real_t f_fn(const Vector &xvec)
+{
+   int dim = xvec.Size();
+   real_t x = pi*xvec[0], y = pi*xvec[1];
+
+   if (dim == 2)
+   {
+      return 2*pi2*sin(x)*sin(y);
+   }
+   else // dim == 3
+   {
+      real_t z = pi*xvec[2];
+      return 3*pi2*sin(x)*sin(y)*sin(z);
+   }
+}
@@ -0,0 +1,433 @@
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+class VectorDGDiffusionIntegrator : public BilinearFormIntegrator
+{
+protected:
+   Coefficient *Q = nullptr;
+   MatrixCoefficient *MQ = nullptr;
+   real_t sigma, kappa;
+   int vdim;
+
+   // these are not thread-safe!
+   Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni;
+   DenseMatrix jmat, dshape1, dshape2, mq, adjJ;
+
+public:
+   VectorDGDiffusionIntegrator(real_t s, real_t k, int vd=-1)
+      : sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDiffusionIntegrator(Coefficient &q, real_t s, real_t k, int vd=-1)
+      : Q(&q), sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDiffusionIntegrator(MatrixCoefficient &mq, real_t s, real_t k,
+                               int vd=-1)
+      : MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
+
+   using BilinearFormIntegrator::AssembleFaceMatrix;
+   virtual void AssembleFaceMatrix(const FiniteElement &el1,
+                                   const FiniteElement &el2,
+                                   FaceElementTransformations &Trans,
+                                   DenseMatrix &full_elmat);
+};
+
+class VectorDGDirichletLFIntegrator : public LinearFormIntegrator
+{
+protected:
+   VectorCoefficient &uD;
+   Coefficient *Q = nullptr;
+   MatrixCoefficient *MQ = nullptr;
+   real_t sigma, kappa;
+   int vdim;
+
+   // these are not thread-safe!
+   Vector shape, dshape_dn, nor, nh, ni, uD_vec;
+   DenseMatrix dshape, mq, adjJ;
+
+public:
+   VectorDGDirichletLFIntegrator(VectorCoefficient &u, real_t s, real_t k,
+                                 int vd=-1)
+      : uD(u), sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDirichletLFIntegrator(VectorCoefficient &u, Coefficient &q, real_t s,
+                                 real_t k, int vd=-1)
+      : uD(u), Q(&q), sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDirichletLFIntegrator(VectorCoefficient &u, MatrixCoefficient &mq,
+                                 real_t s, real_t k, int vd=-1)
+      : uD(u), MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               ElementTransformation &Tr,
+                               Vector &elvect) override
+   { MFEM_ABORT("Not implemented."); }
+
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               FaceElementTransformations &Tr,
+                               Vector &elvect) override;
+};
+
+void VectorDGDiffusionIntegrator::AssembleFaceMatrix(
+   const FiniteElement &el1, const FiniteElement &el2,
+   FaceElementTransformations &Trans, DenseMatrix &full_elmat)
+{
+   int dim, ndof1, ndof2, ndofs;
+   bool kappa_is_nonzero = (kappa != 0.);
+   real_t w, wq = 0.0;
+
+   const int sdim = Trans.GetSpaceDim();
+   if (vdim < 0) { vdim = sdim; }
+
+   dim = el1.GetDim();
+   ndof1 = el1.GetDof();
+
+   nor.SetSize(dim);
+   nh.SetSize(dim);
+   ni.SetSize(dim);
+   adjJ.SetSize(dim);
+   if (MQ)
+   {
+      mq.SetSize(dim);
+   }
+
+   shape1.SetSize(ndof1);
+   dshape1.SetSize(ndof1, dim);
+   dshape1dn.SetSize(ndof1);
+   if (Trans.Elem2No >= 0)
+   {
+      ndof2 = el2.GetDof();
+      shape2.SetSize(ndof2);
+      dshape2.SetSize(ndof2, dim);
+      dshape2dn.SetSize(ndof2);
+   }
+   else
+   {
+      ndof2 = 0;
+   }
+
+   ndofs = ndof1 + ndof2;
+   DenseMatrix elmat;
+   elmat.SetSize(ndofs);
+   elmat = 0.0;
+   if (kappa_is_nonzero)
+   {
+      jmat.SetSize(ndofs);
+      jmat = 0.;
+   }
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == nullptr)
+   {
+      // a simple choice for the integration order
+      int order;
+      if (ndof2)
+      {
+         order = 2 * std::max(el1.GetOrder(), el2.GetOrder());
+      }
+      else
+      {
+         order = 2 * el1.GetOrder();
+      }
+      ir = &IntRules.Get(Trans.GetGeometryType(), order);
+   }
+
+   // assemble: < {(Q \nabla u).n},[v] >      --> elmat
+   //           kappa < {h^{-1} Q} [u],[v] >  --> jmat
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: eip2 will only contain valid data if Elem2 exists
+      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
+
+      if (dim == 1)
+      {
+         nor(0) = 2 * eip1.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Trans.Jacobian(), nor);
+      }
+
+      el1.CalcShape(eip1, shape1);
+      el1.CalcDShape(eip1, dshape1);
+      w = ip.weight / Trans.Elem1->Weight();
+      if (ndof2)
+      {
+         w /= 2;
+      }
+      if (!MQ)
+      {
+         if (Q)
+         {
+            w *= Q->Eval(*Trans.Elem1, eip1);
+         }
+         ni.Set(w, nor);
+      }
+      else
+      {
+         nh.Set(w, nor);
+         MQ->Eval(mq, *Trans.Elem1, eip1);
+         mq.MultTranspose(nh, ni);
+      }
+      CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
+      adjJ.Mult(ni, nh);
+      if (kappa_is_nonzero)
+      {
+         wq = ni * nor;
+      }
+
+      // Note: in the jump term, we use 1/h1 = |nor|/det(J1) which is
+      // independent of Loc1 and always gives the size of element 1 in
+      // direction perpendicular to the face. Indeed, for linear transformation
+      //
+      //      |nor|=measure(face)/measure(ref. face),
+      //
+      //      det(J1)=measure(element)/measure(ref. element),
+      //
+      //      and the ratios measure(ref. element)/measure(ref. face)
+      //      are compatible for all element/face pairs.
+      //
+      // For example: meas(ref. tetrahedron)/meas(ref. triangle) = 1/3, and
+      // for any tetrahedron vol(tet)=(1/3)*height*area(base).
+      //
+      // For interior faces: q_e/h_e=(q1/h1+q2/h2)/2.
+
+      dshape1.Mult(nh, dshape1dn);
+      for (int i = 0; i < ndof1; i++)
+         for (int j = 0; j < ndof1; j++)
+         {
+            elmat(i, j) += shape1(i) * dshape1dn(j);
+         }
+
+      if (ndof2)
+      {
+         el2.CalcShape(eip2, shape2);
+         el2.CalcDShape(eip2, dshape2);
+         w = ip.weight / 2 / Trans.Elem2->Weight();
+         if (!MQ)
+         {
+            if (Q)
+            {
+               w *= Q->Eval(*Trans.Elem2, eip2);
+            }
+            ni.Set(w, nor);
+         }
+         else
+         {
+            nh.Set(w, nor);
+            MQ->Eval(mq, *Trans.Elem2, eip2);
+            mq.MultTranspose(nh, ni);
+         }
+         CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
+         adjJ.Mult(ni, nh);
+         if (kappa_is_nonzero)
+         {
+            wq += ni * nor;
+         }
+
+         dshape2.Mult(nh, dshape2dn);
+
+         for (int i = 0; i < ndof1; i++)
+            for (int j = 0; j < ndof2; j++)
+            {
+               elmat(i, ndof1 + j) += shape1(i) * dshape2dn(j);
+            }
+
+         for (int i = 0; i < ndof2; i++)
+            for (int j = 0; j < ndof1; j++)
+            {
+               elmat(ndof1 + i, j) -= shape2(i) * dshape1dn(j);
+            }
+
+         for (int i = 0; i < ndof2; i++)
+            for (int j = 0; j < ndof2; j++)
+            {
+               elmat(ndof1 + i, ndof1 + j) -= shape2(i) * dshape2dn(j);
+            }
+      }
+
+      if (kappa_is_nonzero)
+      {
+         // only assemble the lower triangular part of jmat
+         wq *= kappa;
+         for (int i = 0; i < ndof1; i++)
+         {
+            const real_t wsi = wq * shape1(i);
+            for (int j = 0; j <= i; j++)
+            {
+               jmat(i, j) += wsi * shape1(j);
+            }
+         }
+         if (ndof2)
+         {
+            for (int i = 0; i < ndof2; i++)
+            {
+               const int i2 = ndof1 + i;
+               const real_t wsi = wq * shape2(i);
+               for (int j = 0; j < ndof1; j++)
+               {
+                  jmat(i2, j) -= wsi * shape1(j);
+               }
+               for (int j = 0; j <= i; j++)
+               {
+                  jmat(i2, ndof1 + j) += wsi * shape2(j);
+               }
+            }
+         }
+      }
+   }
+
+   // elmat := -elmat + sigma*elmat^t + jmat
+   if (kappa_is_nonzero)
+   {
+      for (int i = 0; i < ndofs; i++)
+      {
+         for (int j = 0; j < i; j++)
+         {
+            real_t aij = elmat(i, j), aji = elmat(j, i), mij = jmat(i, j);
+            elmat(i, j) = sigma * aji - aij + mij;
+            elmat(j, i) = sigma * aij - aji + mij;
+         }
+         elmat(i, i) = (sigma - 1.) * elmat(i, i) + jmat(i, i);
+      }
+   }
+   else
+   {
+      for (int i = 0; i < ndofs; i++)
+      {
+         for (int j = 0; j < i; j++)
+         {
+            real_t aij = elmat(i, j), aji = elmat(j, i);
+            elmat(i, j) = sigma * aji - aij;
+            elmat(j, i) = sigma * aij - aji;
+         }
+         elmat(i, i) *= (sigma - 1.);
+      }
+   }
+
+   // populate full matrix following github issue #2909
+   full_elmat.SetSize(vdim*(ndof1 + ndof2));
+   full_elmat = 0.0;
+   for (int d=0; d<vdim; ++d)
+   {
+      for (int j=0; j<ndofs; ++j)
+      {
+         int jj = (j < ndof1) ? j + d*ndof1 : j - ndof1 + d*ndof2 + vdim*ndof1;
+         for (int i=0; i<ndofs; ++i)
+         {
+            int ii = (i < ndof1) ? i + d*ndof1 : i - ndof1 + d*ndof2 + vdim*ndof1;
+            full_elmat(ii, jj) += elmat(i, j);
+         }
+      }
+   }
+};
+
+void VectorDGDirichletLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
+{
+   const int dim = el.GetDim();
+   const int sdim = Tr.GetSpaceDim();
+
+   if (vdim < 0) { vdim = sdim; }
+
+   const int ndof = el.GetDof();
+
+   bool kappa_is_nonzero = (kappa != 0.);
+   real_t w;
+
+   nor.SetSize(dim);
+   nh.SetSize(dim);
+   ni.SetSize(dim);
+   adjJ.SetSize(dim);
+   if (MQ)
+   {
+      mq.SetSize(dim);
+   }
+
+   shape.SetSize(ndof);
+   dshape.SetSize(ndof, dim);
+   dshape_dn.SetSize(ndof);
+
+   elvect.SetSize(vdim * ndof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      // a simple choice for the integration order; is this OK?
+      int order = 2*el.GetOrder();
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
+
+      uD.Eval(uD_vec, Tr, ip);
+
+      if (dim == 1)
+      {
+         nor(0) = 2*eip.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Tr.Jacobian(), nor);
+      }
+
+      el.CalcShape(eip, shape);
+      el.CalcDShape(eip, dshape);
+
+      // compute uD through the face transformation
+      w = ip.weight / Tr.Elem1->Weight();
+      if (!MQ)
+      {
+         if (Q)
+         {
+            w *= Q->Eval(*Tr.Elem1, eip);
+         }
+         ni.Set(w, nor);
+      }
+      else
+      {
+         nh.Set(w, nor);
+         MQ->Eval(mq, *Tr.Elem1, eip);
+         mq.MultTranspose(nh, ni);
+      }
+      CalcAdjugate(Tr.Elem1->Jacobian(), adjJ);
+      adjJ.Mult(ni, nh);
+
+      dshape.Mult(nh, dshape_dn);
+
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         for (int i = 0; i < ndof; ++i)
+         {
+            elvect[i + vd*ndof] += sigma * uD_vec[vd] * dshape_dn[i];
+         }
+      }
+      if (kappa_is_nonzero)
+      {
+         for (int vd = 0; vd < vdim; ++vd)
+         {
+            for (int i = 0; i < ndof; ++i)
+            {
+               elvect[i + vd*ndof] += kappa*(ni*nor) * uD_vec[vd] * shape[i];
+            }
+         }
+      }
+   }
+}
+
+} // namespace mfem
@@ -266,11 +266,7 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)

      // Gather the attributes on the host from all the elements
      const Mesh &mesh = *trial_fes->GetMesh();
-      elem_attributes.SetSize(mesh.GetNE());
-      for (int i = 0; i < mesh.GetNE(); ++i)
-      {
-         elem_attributes[i] = mesh.GetAttribute(i);
-      }
+      elem_attributes = &mesh.GetElementAttributes();
   }

   // Construct face restriction operators only if the bilinear form has
@@ -329,45 +325,7 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
         bdr_face_dYdn.SetSize(bdr_face_restrict_lex->Height());
      }

-      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.");
-      }
+      bdr_face_attributes = &trial_fes->GetMesh()->GetBdrFaceAttributes();
   }
 }

@@ -429,7 +387,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
         mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
         {
            const int attr = d_attr[e];
-            if (d_m[attr - 1] == 0)
+            if (attr <= 0 || d_m[attr - 1] == 0)
            {
               for (int i = 0; i < nd; ++i)
               {
@@ -450,7 +408,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);
@@ -476,7 +434,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);
      }
   }

@@ -489,7 +447,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_attributes, bdr_face_Y);
+                                        *bdr_face_attributes, bdr_face_Y);
      }
      bdr_face_restrict_lex->AddAbsMultTranspose(bdr_face_Y, y);
   }
@@ -588,7 +546,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)
         {
@@ -690,8 +648,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_attributes,
-                               false, bdr_face_Y);
+            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i],
+                               *bdr_face_attributes, false, bdr_face_Y);
         }
         for (int i = 0; i < n_bdr_face_integs; ++i)
         {
@@ -699,12 +657,14 @@ void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
            {
               AddMultNormalDerivativesWithMarkers(
                  *bdr_face_integs[i], bdr_face_X, bdr_face_dXdn,
-                  bdr_face_markers[i], bdr_attributes, bdr_face_Y, bdr_face_dYdn);
+                  bdr_face_markers[i], *bdr_face_attributes, bdr_face_Y,
+                  bdr_face_dYdn);
            }
            else
            {
-               AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X, bdr_face_markers[i],
-                                  bdr_attributes, false, bdr_face_Y);
+               AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X,
+                                  bdr_face_markers[i], *bdr_face_attributes, false,
+                                  bdr_face_Y);
            }
         }
         bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
@@ -727,7 +687,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);
@@ -774,13 +734,14 @@ 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_attributes,
-                               true, bdr_face_Y);
+            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i],
+                               *bdr_face_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_attributes, true, bdr_face_Y);
+            AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X,
+                               bdr_face_markers[i], *bdr_face_attributes, true,
+                               bdr_face_Y);
         }
         bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
      }
@@ -804,7 +765,7 @@ static void AddWithMarkers_(
   mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
   {
      const int attr = d_attr[e];
-      if (d_m[attr - 1] == 0) { return; }
+      if (attr <= 0 || d_m[attr - 1] == 0) { return; }
      for (int i = 0; i < nd; ++i)
      {
         d_y(i, e) += d_x(i, e);
@@ -920,7 +881,8 @@ void EABilinearFormExtension::Assemble()
      {
         const int i = idx % sz;
         const int e = idx / sz;
-         const real_t val = d_m[d_a[e] - 1] ? d_ea_1(i, e) : 0.0;
+         const real_t val =
+            d_a[e] > 0 ? (d_m[d_a[e] - 1] ? d_ea_1(i, e) : 0) : 0;
         if (add)
         {
            d_ea_2(i, e) += val;
@@ -953,7 +915,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);
         }
      }
   }
@@ -982,7 +944,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_attributes, add);
+                             *bdr_face_attributes, add);
         }
      }
   }
@@ -1031,7 +993,7 @@ void EABilinearFormExtension::Assemble()
                                                        ea_data_tmp,
                                                        add);
            add_with_markers(ea_data_tmp, ea_data_bdr, nf_bdr, *markers,
-                             bdr_attributes, add);
+                             *bdr_face_attributes, add);
         }
      }
   }
@@ -69,7 +69,8 @@ class PABilinearFormExtension : public BilinearFormExtension
 protected:
   const FiniteElementSpace *trial_fes, *test_fes; // Not owned
   /// Attributes of all mesh elements.
-   Array<int> elem_attributes, bdr_attributes;
+   const Array<int> *elem_attributes; // Not owned
+   const Array<int> *bdr_face_attributes; // Not owned
   mutable Vector tmp_evec; // Work array
   mutable Vector localX, localY;
   mutable Vector int_face_X, int_face_Y;
@@ -812,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 );
   }

@@ -884,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 );
   }

@@ -919,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 );
   }

@@ -1600,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 );
   }

@@ -1635,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 );
   }

@@ -1669,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 );
   }

@@ -1760,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 );
   }

@@ -1793,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   );
   }
@@ -1832,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
             );
@@ -1973,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 );
   }

@@ -912,7 +912,7 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,

   if (vdim == 1) // scalar case
   {
-      n_field["values"].set_external(gf->GetData(),
+      n_field["values"].set_external(const_cast<real_t *>(gf->HostRead()),
                                     ndofs);
   }
   else // vector case
@@ -925,18 +925,18 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,
      int vdim_stride  = (ordering == Ordering::byNODES ? ndofs : 1);

      index_t offset = 0;
-      index_t stride = sizeof(double) * entry_stride;
+      index_t stride = sizeof(real_t) * entry_stride;

      for (int d = 0;  d < vdim; d++)
      {
         std::ostringstream oss;
         oss << "v" << d;
         std::string comp_name = oss.str();
-         n_field["values"][comp_name].set_external(gf->GetData(),
+         n_field["values"][comp_name].set_external(const_cast<real_t *>(gf->HostRead()),
                                                   ndofs,
                                                   offset,
                                                   stride);
-         offset +=  sizeof(double) * vdim_stride;
+         offset +=  sizeof(real_t) * vdim_stride;
      }
   }

@@ -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, MPI_INT,
-              lengths.data(), 1, MPI_INT,
+   MPI_Gather(&msg_len, 1, MPITypeMap<size_t>::mpi_type,
+              lengths.data(), 1, MPITypeMap<size_t>::mpi_type,
              0, MPI_COMM_WORLD);

   if (myrank == 0)
@@ -568,7 +568,7 @@ struct ThreadBlocks
   int z = 1;
 };

-#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if defined(MFEM_USE_CUDA_OR_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) || defined(MFEM_USE_HIP))
+#if defined(MFEM_USE_CUDA_OR_HIP)
      // int gridsize = (N + Z - 1) / Z;
      int num_bytes = num_shmem * sizeof(decltype(shmem));
      dim3 block_size(blocks.x, blocks.y, blocks.z);
@@ -401,9 +401,6 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
   {
      MFEM_ABORT("unknown FiniteElementCollection: " << name);
   }
-   MFEM_VERIFY(!strcmp(fec->Name(), name), "input name: \"" << name
-               << "\" does not match the created collection name: \""
-               << fec->Name() << '"');

   return fec;
 }
@@ -922,6 +922,9 @@ public:
   { return mesh->GetBdrElementType(i); }

   /// Returns ElementTransformation for the @a i-th element.
+   /// @note The returned pointer references an object owned by the associated
+   /// @a Mesh that will be modified by other calls to `GetElementTransformation`.
+   /// As such, this pointer should @b not be deleted by the caller.
   ElementTransformation *GetElementTransformation(int i) const
   { return mesh->GetElementTransformation(i); }

@@ -59,26 +59,23 @@ 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)
   {
-      MFEM_FOREACH_THREAD(i, x, NQ)
+      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)
      {
-         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);
-      }
-   });
+         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);
+      });
+   }
 }

 void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
@@ -109,50 +106,22 @@ void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
   CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);

   const int NE = ne;
-   const int Q1D = quad1D;
+   const int NQ = nq;
   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(), 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 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)
      {
-         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);
-         }
+         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);
      });
   }
-   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)
@@ -346,13 +346,13 @@ private:
      template<typename T>
      T operator() (const blitz::TinyVector<T,3>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
-         TmplPoly_1D::CalcBernstein(el_order, x[2], u3);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[2], u3.data());

         const Array<int>& dof_map=el->GetDofMap();

@@ -370,17 +370,17 @@ private:
      template<typename T>
      blitz::TinyVector<T,3> grad(const blitz::TinyVector<T,3>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
-         T d3[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);
+         std::vector<T> d3(el_order+1);

-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
-         TmplPoly_1D::CalcBernstein(el_order,x[2], u3, d3);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[2], u3.data(), d3.data());

         blitz::TinyVector<T,3> res(T(0.0),T(0.0),T(0.0));

@@ -415,11 +415,11 @@ private:
      template<typename T>
      T operator() (const blitz::TinyVector<T,2>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());

         const Array<int>& dof_map=el->GetDofMap();

@@ -437,14 +437,14 @@ private:
      template<typename T>
      blitz::TinyVector<T,2> grad(const blitz::TinyVector<T,2>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);

-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());

         blitz::TinyVector<T,2> res(T(0.0),T(0.0));

@@ -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,13 +93,14 @@ void LinearFormExtension::Assemble()
      else
      {
         // scan the attributes to set the markers to 0 or 1
-         const int NBE = bdr_attributes.Size();
-         const auto attr = bdr_attributes.Read();
+         const int NBE = bdr_face_attributes->Size();
+         const auto attr = bdr_face_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_markers[attr[e]-1] == 1;
+            markers_w[e] =
+               attr[e] > 0 ? (attr_markers[attr[e] - 1] == 1) : false;
         });
      }

@@ -125,8 +126,7 @@ void LinearFormExtension::Update()
      //markers.UseDevice(true);

      // Gather the attributes on the host from all the elements
-      attributes.SetSize(NE);
-      for (int i = 0; i < NE; ++i) { attributes[i] = mesh.GetAttribute(i); }
+      attributes = &mesh.GetElementAttributes();

      elem_restrict_lex = fes.GetElementRestriction(ordering);
      MFEM_VERIFY(elem_restrict_lex, "Element restriction not available");
@@ -136,34 +136,12 @@ void LinearFormExtension::Update()

   if (lf->boundary_integs.Size() > 0)
   {
-      const int nf_bdr = fes.GetNFbyType(FaceType::Boundary);
+      bdr_face_attributes = &mesh.GetBdrFaceAttributes();
+
+      const int nf_bdr = bdr_face_attributes->Size();
      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,
@@ -25,7 +25,8 @@ class LinearForm;
 class LinearFormExtension
 {
   /// Attributes of all mesh elements.
-   Array<int> attributes, bdr_attributes;
+   const Array<int> *attributes; // Not owned
+   const Array<int> *bdr_face_attributes; // Not owned

   /// Temporary markers for device kernels.
   Array<int> markers, bdr_markers;
@@ -214,16 +214,16 @@ public:
 inline void QuadratureFunction::GetValues(
   int idx, Vector &values)
 {
-   const int s_offset = qspace->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(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->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(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->offsets[idx] * vdim + ip_num * vdim;
+   const int s_offset = qspace->Offset(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->offsets[idx] * vdim + ip_num * vdim;
+   const int s_offset = qspace->Offset(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->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(idx + 1) - s_offset;
   // Make the values matrix memory an alias of the quadrature function memory
   Memory<real_t> &values_mem = values.GetMemory();
   values_mem.Delete();
@@ -268,8 +268,8 @@ inline void QuadratureFunction::GetValues(
 inline void QuadratureFunction::GetValues(
   int idx, DenseMatrix &values) const
 {
-   const int s_offset = qspace->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(idx + 1) - s_offset;
   values.SetSize(vdim, sl_size);
   values.HostWrite();
   const real_t *q = HostRead() + vdim*s_offset;
@@ -97,16 +97,31 @@ void QuadratureSpaceBase::Integrate(VectorCoefficient &coeff,
 void QuadratureSpace::ConstructOffsets()
 {
   const int num_elem = mesh.GetNE();
-   offsets.SetSize(num_elem + 1);
-   int offset = 0;
-   for (int i = 0; i < num_elem; i++)
+   ne = num_elem;
+
+   if (mesh.GetNumGeometries(mesh.Dimension()) == 1)
   {
-      offsets[i] = offset;
-      int geom = mesh.GetElementBaseGeometry(i);
-      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule.");
-      offset += int_rule[geom]->GetNPoints();
+      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[num_elem] = size = offset;
 }

 void QuadratureSpace::Construct()
@@ -188,6 +203,7 @@ void FaceQuadratureSpace::ConstructOffsets()
 {
   face_indices.SetSize(num_faces);
   offsets.SetSize(num_faces + 1);
+   ne = num_faces;
   int offset = 0;
   int f_idx = 0;
   for (int i = 0; i < mesh.GetNumFacesWithGhost(); i++)
@@ -30,13 +30,17 @@ protected:
   Mesh &mesh; ///< The underlying mesh.
   int order; ///< The order of integration rule.
   int size; ///< Total number of quadrature points.
+   int ne; ///< Actual number of entities
   mutable Vector weights; ///< Integration weights.
   mutable long nodes_sequence = 0; ///< Nodes counter for cache invalidation.

-   /// @brief Entity quadrature point offset array, of size num_entities + 1.
+   /// @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.
   ///
-   /// 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];
@@ -59,6 +63,24 @@ protected:
   void ConstructWeights() const;

 public:
+   /// @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
+   {
+      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 { return offsets; }
+
   /// Return the total number of quadrature points.
   int GetSize() const { return size; }

@@ -66,7 +88,7 @@ public:
   int GetOrder() const { return order; }

   /// Return the number of entities.
-   int GetNE() const { return offsets.Size() - 1; }
+   int GetNE() const { return ne; }

   /// Returns the mesh.
   inline Mesh *GetMesh() const { return &mesh; }
@@ -14,7 +14,7 @@

 #include "../config/config.hpp"

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #include <cusparse.h>
 #include <library_types.h>
 #include <cuda_runtime.h>
@@ -22,7 +22,7 @@
 #endif
 #include "cuda.hpp"

-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
 #include <hip/hip_runtime.h>
 #endif
 #include "hip.hpp"
@@ -43,7 +43,7 @@
 #endif
 #endif

-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if !defined(MFEM_USE_CUDA_OR_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;
@@ -347,13 +347,23 @@ 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;
-const MPI_Datatype MPITypeMap<int64_t>::mpi_type = MPI_INT64_T;
-const MPI_Datatype MPITypeMap<uint64_t>::mpi_type = MPI_UINT64_T;
-
+/// \endcond DO_NOT_DOCUMENT

 GroupCommunicator::GroupCommunicator(const GroupTopology &gt, Mode m)
   : gtopo(gt), mode(m)
@@ -607,10 +607,50 @@ 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;
@@ -619,15 +659,6 @@ 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).
@@ -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)
+             << cudaGetErrorString(err) << " [code: " << (int)err << ']'
             << "\n ... in function: " << func
             << "\n ... in file: " << file << ':' << line << '\n';
   mfem_error();
@@ -18,7 +18,7 @@
 // CUDA block size used by MFEM.
 #define MFEM_CUDA_BLOCKS 256

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,10 +37,9 @@
                              __FILE__, __LINE__);                             \
    }                                                                          \
  } while (0)
-#endif // MFEM_USE_CUDA

 // Define the MFEM inner threading macros
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
+#if defined(__CUDA_ARCH__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) blockIdx.k
@@ -48,12 +47,13 @@
 #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
+#endif // defined(__CUDA_ARCH__)
+#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)

 namespace mfem
 {

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 // 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);
@@ -176,7 +176,7 @@ __device__ void abort_msg(T & msg)
      printf(__VA_ARGS__);     \
      asm("trap;");            \
   }
-#elif defined(MFEM_USE_HIP)
+#elif defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_ABORT_KERNEL(...) \
   {                           \
      printf(__VA_ARGS__);     \
@@ -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__)
-#ifdef __NVCC__
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) // Clang cuda or nvcc
+#ifdef __NVCC__ // nvcc specifically
 #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). If neither CUDA nor
-// HIP is enabled, this macro is a no-op.
-#if defined(MFEM_USE_CUDA)
+// 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__)
 #define MFEM_GPU_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
                                       (int i) {__VA_ARGS__})
-#elif defined(MFEM_USE_HIP)
+#elif defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_GPU_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
                                        (int i) {__VA_ARGS__})
 #else
@@ -225,6 +225,52 @@ 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)
@@ -454,6 +500,42 @@ 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


@@ -481,7 +563,7 @@ void RajaSeqWrap(const int N, HBODY &&h_body)


 /// CUDA backend
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)

 template <typename BODY> __global__ static
 void CuKernel1D(const int N, BODY body)
@@ -573,11 +655,11 @@ struct CuWrap<3>
   }
 };

-#endif // MFEM_USE_CUDA
+#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)


 /// HIP backend
-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)

 template <typename BODY> __global__ static
 void HipKernel1D(const int N, BODY body)
@@ -668,7 +750,7 @@ struct HipWrap<3>
   }
 };

-#endif // MFEM_USE_HIP
+#endif // defined(MFEM_USE_HIP) && defined(__HIP__)


 /// The forall kernel body wrapper
@@ -701,7 +783,7 @@ inline void ForallWrap(const bool use_dev, const int N,
   }
 #endif

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
   // If Backend::CUDA is allowed, use it
   if (Device::Allows(Backend::CUDA))
   {
@@ -709,7 +791,7 @@ inline void ForallWrap(const bool use_dev, const int N,
   }
 #endif

-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
   // If Backend::HIP is allowed, use it
   if (Device::Allows(Backend::HIP))
   {
@@ -753,6 +835,83 @@ 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)
 {
@@ -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)
+             << hipGetErrorString(err) << " [code: " << (int)err << ']'
             << "\n ... in function: " << func
             << "\n ... in file: " << file << ':' << line << '\n';
   mfem_error();
@@ -18,7 +18,7 @@
 // HIP block size used by MFEM.
 #define MFEM_HIP_BLOCKS 256

-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,10 +37,9 @@
                             __FILE__, __LINE__);                              \
    }                                                                          \
  } while (0)
-#endif // MFEM_USE_HIP

 // Define the MFEM inner threading macros
-#if defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)
+#if defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) hipBlockIdx_ ##k
@@ -50,7 +49,8 @@
   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
+#endif // defined(__HIP_DEVICE_COMPILE__)
+#endif // defined(MFEM_USE_HIP) && defined(__HIP__)

 namespace mfem
 {
@@ -513,10 +513,7 @@ 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
-   // 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); }
+   { return HipMemcpyDtoD(dst, src, bytes); }
   void *DtoH(void *dst, const void *src, size_t bytes) override
   { return HipMemcpyDtoH(dst, src, bytes); }
 };
@@ -658,10 +655,7 @@ public:
      return CuMemcpyDtoD(dst, src, bytes);
 #endif
 #ifdef MFEM_USE_HIP
-      // 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);
+      return HipMemcpyDtoD(dst, src, bytes);
 #endif
      // rm.copy(dst, const_cast<void*>(src), bytes); return dst;
   }
@@ -896,6 +896,7 @@ 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;
@@ -1057,7 +1058,8 @@ 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
-         MemoryManager::Exists() && HypreUsingGPU()
+         IsDeviceMemory(MemoryManager::GetDeviceMemoryType()) ||
+         (MemoryManager::Exists() && HypreUsingGPU())
 #endif
      )
      {
@@ -537,7 +537,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
      return;
   }

-#if defined(MFEM_USE_HIP) || defined(MFEM_USE_CUDA)
+#if defined(MFEM_USE_CUDA_OR_HIP)
   if (use_dev &&
       mfem::Device::Allows(Backend::CUDA | Backend::HIP | Backend::RAJA_CUDA |
                            Backend::RAJA_HIP))
@@ -117,7 +117,7 @@ Vector::Vector(const Vector &v)
   UseDevice(v.UseDevice());
 }

-Vector::Vector(Vector &&v)
+Vector::Vector(Vector &&v) : Vector()
 {
   *this = std::move(v);
 }
@@ -32,6 +32,7 @@
 #include <cstring>
 #include <ctime>
 #include <functional>
+#include <set>
 #include <numeric>
 #include <unordered_map>
 #include <unordered_set>
@@ -921,6 +922,63 @@ 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::DeleteGeometricFactors()
 {
   for (int i = 0; i < geom_factors.Size(); i++)
@@ -1804,6 +1862,8 @@ void Mesh::Destroy()
   TetMemory.Clear();
 #endif

+   elem_attrs_cache.DeleteAll();
+   bdr_face_attrs_cache.DeleteAll();
   attributes.DeleteAll();
   bdr_attributes.DeleteAll();
 }
@@ -1817,36 +1877,48 @@ 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()
+void Mesh::SetAttributes(bool elem_attrs_changed, bool bdr_face_attrs_changed)
 {
-   Array<int> attribs;
+   if (bdr_face_attrs_changed)
+   {
+      bdr_face_attrs_cache.SetSize(0); // Invalidate the cache

-   attribs.SetSize(GetNBE());
-   for (int i = 0; i < attribs.Size(); i++)
-   {
-      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!");
+      // 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.SetSize(GetNE());
-   for (int i = 0; i < attribs.Size(); i++)
+   if (elem_attrs_changed)
   {
-      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!");
+      // 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!");
+      }
   }
 }

@@ -4454,6 +4526,10 @@ 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()
@@ -7721,6 +7797,12 @@ 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);
 }

@@ -10850,6 +10932,10 @@ 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);
 }

 void Mesh::GetElementData(const Array<Element*> &elem_array, int geom,
@@ -107,6 +107,11 @@ 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(),
@@ -1122,13 +1127,14 @@ 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 and
-   /// boundary elements.
+   /// @brief Determine the sets of unique attribute values in domain if @a
+   /// elem_attrs_changed and boundary elements if @a bdr_face_attrs_changed.
   ///
   /// 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();
+   virtual void SetAttributes(bool elem_attrs_changed = true,
+                              bool bdr_face_attrs_changed = true);

   /// Check (and optionally attempt to fix) the orientation of the elements
   /** @param[in] fix_it  If `true`, attempt to fix the orientations of some
@@ -2267,6 +2273,35 @@ 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
@@ -1590,21 +1590,27 @@ void ParMesh::DistributeAttributes(Array<int> &attr)
   delete [] glb_attr_marker;
 }

-void ParMesh::SetAttributes()
+void ParMesh::SetAttributes(bool elem_attrs_changed, bool bdr_attrs_changed)
 {
   // Determine the attributes occurring in local interior and boundary elements
-   Mesh::SetAttributes();
+   Mesh::SetAttributes(elem_attrs_changed, bdr_attrs_changed);

-   DistributeAttributes(bdr_attributes);
-   if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
+   if (bdr_attrs_changed)
   {
-      MFEM_WARNING("Non-positive boundary element attributes found!");
+      DistributeAttributes(bdr_attributes);
+      if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
+      {
+         MFEM_WARNING("Non-positive boundary element attributes found!");
+      }
   }

-   DistributeAttributes(attributes);
-   if (attributes.Size() > 0 && attributes[0] <= 0)
+   if (elem_attrs_changed)
   {
-      MFEM_WARNING("Non-positive element attributes found!");
+      DistributeAttributes(attributes);
+      if (attributes.Size() > 0 && attributes[0] <= 0)
+      {
+         MFEM_WARNING("Non-positive element attributes found!");
+      }
   }
 }

@@ -394,7 +394,8 @@ public:

   void Finalize(bool refine = false, bool fix_orientation = false) override;

-   void SetAttributes() override;
+   void SetAttributes(bool elem_attrs_changed = true,
+                      bool bdr_attrs_changed = true) override;

   /// Checks if any rank in the mesh has boundary elements
   bool HasBoundaryElements() const override;
@@ -105,8 +105,11 @@ MFEM_PERF_CXXFLAGS_xlc = -mcpu=native
 # - Clang extra options:
 ifeq ($(MFEM_MACHINE),riscv64)
  MFEM_PERF_CXXFLAGS_clang += -march=rv64gc
-else ifneq ($(MFEM_MACHINE),arm64)
-  # -march=native is unavailable on clang/ARM64 as of 05/2021: support could be added later.
+else ifneq (,$(findstring ppc,$(MFEM_MACHINE)))
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else ifeq ($(MFEM_MACHINE),arm64)
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else
  MFEM_PERF_CXXFLAGS_clang += -march=native
 endif
 MFEM_PERF_CXXFLAGS_clang += $(PEDANTIC_FLAG) -Wall
Author	SHA1	Message	Date
Will Pazner	35d13bebce	Add DG vector diffusion example	2025-09-04 21:57:02 -07:00
Veselin Dobrev	cd377f5201	Merge pull request #4937 from mfem/fec-new-fix bugfix for 4935	2025-09-04 12:19:26 -07:00
Veselin Dobrev	2d57401575	Enzyme/LLVM fixes in GitHub CI (#4997 ) * Update GitHub CI to handle Enzyme/LLVM changes in Homebrew. * GitHub CI: try to allow the Enzyme job to fail * GitHub CI: fix a typo in last commit * GitHub CI: another try to allow the Enzyme job to fail. * GitHub CI: another try to allow the Enzyme job to fail. * GitHub CI: another try to allow the Enzyme job to fail. * GitHub CI: another try to allow the Enzyme job to fail. * GitHub CI: try without fail-fast == false. * GitHub CI: try to fix the LLVM link issue in the Enzyme job.	2025-08-29 08:16:41 -07:00
Veselin Dobrev	7f788e83b9	Merge pull request #4926 from mfem/dev/docs-add Add description to GetElementTransformation	2025-08-26 12:16:19 -07:00
Veselin Dobrev	7b7f77379e	Merge branch 'master' into fec-new-fix	2025-08-19 16:54:27 -07:00
Will Pazner	f0cb31088c	Merge pull request #4989 from mfem/update-ci-mac Update Xcode version in macos CI from 15.3 -> 16.4	2025-08-19 16:32:39 -07:00
Andrew Ho	41cccee855	Merge branch 'master' into fec-new-fix	2025-08-19 15:17:23 -07:00
Justin Laughlin	f407ca7756	Xcode 16.4	2025-08-19 13:58:39 -07:00
Justin Laughlin	ab00472c5d	Try removing xcode version specification	2025-08-19 13:46:38 -07:00
Justin Laughlin	c4a3d31289	Update Xcode version in macos CI from 15.3 -> 16.4	2025-08-19 13:38:07 -07:00
Tzanio Kolev	20e77be88a	Merge pull request #4987 from mfem/gitlab-ci--ruby-to-dane Move Gitlab CI from Ruby to Dane	2025-08-16 16:01:21 -07:00
Andrew Ho	d70b99c4f3	Merge branch 'master' into fec-new-fix	2025-08-16 00:29:08 -07:00
Veselin Dobrev	97a568af13	Move Gitlab CI from Ruby to Dane	2025-08-15 21:27:43 -07:00
Tzanio Kolev	cd2715680f	Merge pull request #4725 from mfem/qspace-device Bilinear Operator Assembly optimizations	2025-08-15 07:22:53 -07:00
Tzanio Kolev	62997e9dc7	Merge pull request #4928 from mfem/mpi-dfem-fix change MPI_INT to MFEM_MPI_SIZE_T	2025-08-15 07:22:37 -07:00
Veselin Dobrev	91cdbeee7d	Merge pull request #4919 from mfem/revert-hip-memcpy-async-for-d2d Use `HipMemcpyDtoD` instead of `HipMemcpyDtoDAsync`	2025-08-12 13:58:36 -07:00
Tzanio Kolev	0466ee74f2	Merge pull request #4622 from mfem/hypre-gpu-tweaks HYPRE and GPU related tweaks	2025-08-09 10:28:25 -07:00
Tzanio Kolev	50b8873f20	Merge pull request #4923 from mfem/hughcars/cuda-hip-compiler-checks Check for CUDA or HIP compilers	2025-08-09 10:27:54 -07:00
Veselin Dobrev	cd7dfe2188	Merge pull request #4973 from mfem/macos-llvm-ci switch to llvm 20 for CI	2025-08-07 17:55:50 -07:00
Andrew Ho	1bf00d3761	switch to llvm20 for CI	2025-08-07 15:18:36 -07:00
Andrew Ho	b2a0f48001	fixed comment	2025-08-06 09:19:47 -07:00
Veselin Dobrev	fb200e4b98	In the 3D version of 'forall', fix (Raja)OmpWrap2D calls to be 3D. In INSTALL, add version requirement: OpenMP >= 3.1. Fix a warning in the Vector move constructor about use of uninitialized values by default initializing the Vector vefore calling its move assignment operator.	2025-08-06 00:36:29 -07:00
Andrew Ho	da852fd0e4	Merge remote-tracking branch 'base/qspace-device' into qspace-device	2025-08-05 18:10:35 -07:00
Andrew Ho	c2657d20a8	Check for negative element attributes in debug mode when computing elem_attrs_cache	2025-08-05 18:08:42 -07:00
Veselin Dobrev	47e9b29c0f	Restore the MFEM_UNROLL definition -- it is only used in CUDA device code. Tweak a comment and restore some formatting.	2025-08-05 16:27:15 -07:00
Andrew Ho	cbcae09fe8	Merge branch 'master' into qspace-device	2025-08-05 14:30:11 -07:00
Andrew Ho	b75c6e96e6	fix doxygen	2025-08-05 14:19:51 -07:00
Andrew Ho	a54c5d2105	changed to use MPITypeMap extend MPITypeMap to all fundamental C++ types except long double	2025-08-05 14:09:11 -07:00
Andrew Ho	cd5d098c95	Merge branch 'master' into mpi-dfem-fix	2025-08-05 13:56:30 -07:00
Andrew Ho	30506b52d9	Documentation	2025-08-05 12:29:02 -07:00
Andrew Ho	3f78374d2a	ignore attributes <= 0	2025-08-05 12:12:54 -07:00
Veselin Dobrev	bd7281bdfa	Merge pull request #4960 from mfem/fix-warnings-algoim Algoim compilation warnings	2025-08-05 12:05:40 -07:00
Andrew Ho	1299d05d21	Merge branch 'master' into mpi-dfem-fix	2025-08-05 11:11:22 -07:00
Veselin Dobrev	39944d8dfb	In the INSTALL file, add requirement for ROCm >= 5.6.1 due to the fix in that version that made HipMemcpyDtoD asynchronous w.r.t. host.	2025-08-05 09:54:53 -07:00
Hugh Carson	6df8ae80cb	PR feedback	2025-08-05 11:17:12 -04:00
John Camier	0730fa304a	Merge branch 'master' into revert-hip-memcpy-async-for-d2d	2025-08-05 06:06:02 -07:00
Tzanio Kolev	6dfbe25c07	Merge pull request #4955 from mfem/clang-perf-flags In `miniapps/performance/makefile`, adjust the flags for Clang	2025-08-04 07:44:47 -07:00
Andrew Ho	a8e251c28e	2D and 3D forall for OpenMP	2025-08-02 14:39:24 -07:00
Andrew Ho	9df402b565	suggestions from Veselin	2025-08-02 13:34:13 -07:00
Andrew Ho	5ddbaed1f4	Merge branch 'master' into qspace-device	2025-08-02 12:52:12 -07:00
Tzanio Kolev	02edb95a31	Merge pull request #4957 from farscape-project/conduit Guarantee data on host when writing conduit files	2025-08-02 09:27:45 -07:00
Andrew Ho	2d6cb3f0e5	Merge remote-tracking branch 'base/qspace-device' into qspace-device	2025-07-29 17:09:37 -07:00
Andrew Ho	d06f52461e	missing bdr_attrs_cache	2025-07-29 17:09:10 -07:00
Andrew Ho	3853cb6860	fixed typo	2025-07-29 17:07:45 -07:00
Andrew HoandVeselin Dobrev	f5d0db2a33	Update mesh/mesh.cpp Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>	2025-07-29 17:07:04 -07:00
Andrew Ho	a102daf621	Merge branch 'master' into qspace-device	2025-07-29 17:04:33 -07:00
Veselin Dobrev	6d8a96a8fb	Merge branch 'master' into hypre-gpu-tweaks	2025-07-29 09:08:02 -07:00
Vladimir Z Tomov	8d9579057f	Simplified to std:vector.	2025-07-28 15:55:21 -07:00
Vladimir Z Tomov	76d65f5866	Fixed compilation warnings when Algoim is on.	2025-07-28 15:39:50 -07:00
Nuno Nobre	bbde976378	Guarantee data on host when writing conduit files	2025-07-26 01:31:22 +01:00
Veselin Dobrev	447cd0f500	In miniapps/performance/makefile, adjust the flags for Clang	2025-07-22 22:14:35 -07:00
Andrew Ho	4cc9b5c6f6	Merge branch 'master' into qspace-device	2025-07-19 17:07:43 -07:00
Andrew Ho	2d3ec4ca67	remove the output name check sometimes there are equivalent basis with different names created	2025-07-16 08:03:16 -07:00
Gabriel Esteban Pinochet Soto	769f672ac1	Fix style	2025-07-15 07:22:15 -07:00
Gabriel Pinochet-SotoandChris Vogl	dca2a24af2	Update fem/fespace.hpp Co-authored-by: Chris Vogl <vogl2@llnl.gov>	2025-07-14 21:40:40 -07:00
Gabriel Pinochet-Soto	019194d42d	Merge branch 'master' into dev/docs-add	2025-07-14 21:39:43 -07:00
Andrew Ho	267082d6f2	Merge remote-tracking branch 'base/master' into qspace-device	2025-07-12 23:07:56 -07:00
Andrew Ho	dd6c207534	formatting	2025-07-11 10:02:11 -07:00
Andrew Ho	71af00ab6d	change MPI_INT to MFEM_MPI_SIZE_T this detects the underlying size of size_t	2025-07-11 09:57:25 -07:00
Gabriel Esteban Pinochet Soto	03da0c870c	Add description to GetElementTransformation	2025-07-09 08:48:47 -07:00
Hugh Carson	3f447cb11c	Check if a cuda or hip compiler is being used before accessing cuda or hip specific variable	2025-07-07 15:47:32 -04:00
Tom Stitt	ce434a9c82	HipMemcpyDtoDAsync -> HipMemcpyDtoD since hipMemcpy for d2d is async wip to the host as of 5.6.1	2025-07-03 09:39:58 -07:00
Andrew Ho	1ac7837dc0	Merge branch 'master' into qspace-device	2025-06-16 12:20:48 -07:00
Andrew Ho	056152aaa6	Merge branch 'master' into qspace-device	2025-06-02 11:38:14 -07:00
Andrew Ho	02424a41a9	Merge branch 'master' into qspace-device	2025-05-27 12:05:07 -07:00
Andrew Ho	b289def58a	Merge branch 'master' into qspace-device	2025-05-12 14:04:22 -07:00
Andrew Ho	932508ad1b	fixed bugs	2025-05-07 18:49:36 -07:00
Andrew Ho	0e30ee1e15	Added forall wrappers which flatten 2D and 3D ranges achieves good performance on CPU and GPU	2025-05-07 18:25:27 -07:00
Andrew Ho	dd236b8796	Improved the performance of the 1D forall modulo assembly for CPUs	2025-05-06 22:44:39 -07:00
Andrew Ho	9a47f4391d	make sure CPU is also fast for mass matrix assembly	2025-05-06 01:15:25 -07:00
Andrew Ho	ccc006febf	Merge branch 'master' into qspace-device	2025-05-05 13:13:57 -06:00
Andrew Ho	f42334cb4d	Merge branch 'master' into qspace-device	2025-05-02 10:29:15 -06:00
Andrew Ho	cc86e59167	Made getting offsets public, added accessor to internal held offsets	2025-04-23 13:25:28 -07:00
Andrew Ho	52eaf0deb7	Merge branch 'master' into qspace-device	2025-04-23 13:16:59 -07:00
Andrew Ho	44f58a169c	Merge branch 'master' into qspace-device	2025-04-15 11:01:58 -07:00
Andrew Ho	4a09b7a1fc	fix merge compile errors	2025-04-04 16:28:42 -07:00
Andrew Ho	1ea8e59f69	Merge branch 'master' into qspace-device	2025-04-04 16:02:17 -07:00
Andrew Ho	b6aa248fc9	Merge branch 'master' into qspace-device	2025-03-24 06:27:44 -07:00
Andrew Ho	6857e4b277	use set instead of unordered_set to get ordered attributes	2025-02-26 07:14:55 -08:00
Andrew Ho	639a0347fb	Merge branch 'master' into qspace-device	2025-02-25 19:43:38 -08:00
Andrew Ho	a02f7165d8	Cache bdr_attributes in mesh (used by linearform_ext and bilinearform_ext)	2025-02-25 19:41:36 -08:00
Andrew Ho	4fb48b4e41	Get rid of the ElementAttributesUpdated function Have users call SetAttributes() instead	2025-02-25 17:16:43 -08:00
Andrew Ho	751dbd0ed9	Merge remote-tracking branch 'base/qspace-device' into qspace-device	2025-02-25 17:05:04 -08:00
Andrew Ho	c2c6a94d42	Merge remote-tracking branch 'base/qspace-device' into qspace-device	2025-02-25 17:04:49 -08:00
Andrew HoandWill Pazner	088ec0f461	Update fem/qspace.hpp Co-authored-by: Will Pazner <11493037+pazner@users.noreply.github.com>	2025-02-25 17:03:59 -08:00
Andrew Ho	7e1be9da75	Use unordered_set to get unique attributes	2025-02-25 17:00:05 -08:00
Andrew HoandWill Pazner	be29d4a3d4	Update mesh/mesh.hpp Co-authored-by: Will Pazner <11493037+pazner@users.noreply.github.com>	2025-02-25 16:58:32 -08:00
Andrew HoandWill Pazner	be34c5a6bd	Update mesh/mesh.cpp Co-authored-by: Will Pazner <11493037+pazner@users.noreply.github.com>	2025-02-25 16:57:35 -08:00
Andrew Ho	8bd3569f6d	unused variable	2025-02-25 02:31:54 -08:00
Andrew Ho	e8908ad646	cache element attributes in mesh	2025-02-24 19:38:29 -08:00
Andrew Ho	6c7673df1c	code style	2025-02-24 19:18:40 -08:00
Andrew Ho	1f7f394f4e	Merge remote-tracking branch 'base/master' into qspace-device	2025-02-24 18:47:58 -08:00
Veselin Dobrev	4f0f353720	Print error code (in addition to the error string) on CUDA/HIP errors	2024-12-09 18:20:32 -08:00
Veselin Dobrev	a3156023d0	Avoid HYPRE initialization via calls to HypreUsingGPU() and GetHypreMemoryLocation() when running serial codes.	2024-12-09 18:03:01 -08:00
Andrew Ho	7794c291be	Expose more parallelism in mass integrator assemblePA	2024-11-22 13:01:58 -08:00
Andrew Ho	7ca49cb200	FaceQuadratureSpace also needs to save ne	2024-11-22 12:01:20 -08:00
Andrew Ho	e48d41c9dc	return correct NE for compressed qspace	2024-11-22 11:51:20 -08:00
Andrew Ho	0ead1d9c77	Implemented offset compression in qspace	2024-11-22 11:40:33 -08:00
Andrew Ho	e2168f9348	code style	2024-11-20 15:56:17 -08:00
Andrew Ho	65516eca23	offsets need to be accessible to the host	2024-11-20 15:47:43 -08:00
Andrew Ho	79d9314560	Fix issue with using CUDA lambda functions in private/protected member functions	2024-11-19 17:06:07 -08:00
Will Pazner	d91d39bc00	Construct quadrature space offsets on device (when possible) This can be done easily when there is only one geometry type in the mesh	2024-11-19 12:52:22 -08:00