documentation

* 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.

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

@@ -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@20 enzyme
-          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
-          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV
+          brew install 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

.gitlab-ci.yml

@@ -22,7 +22,7 @@ include:
 # the "needs" keyword and express the DAG of jobs for more efficiency.
 # - We use setup and setup_baseline phases to download content outside of mfem
 #   directory.
-# - Allocate/Release is where 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:

.gitlab/README.md

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

.gitlab/configs/common.yml

@@ -24,7 +24,7 @@ variables:
 # TODO: add a clean-up mechanism
   BUILD_ROOT: ${USER_CI_TOP_DIR}/${CI_PROJECT_NAME}-${MACHINE_NAME}-pipeline-${CI_PIPELINE_ID}
 
-# On LLNL's 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}

.gitlab/configs/dane-config.yml

@@ -9,17 +9,17 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-# GitLab pipelines configurations for the 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}

.gitlab/configs/setup-baseline.yml

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

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

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

.gitlab/dane-baseline.yml

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

.gitlab/dane-build-and-test.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

.gitlab/scripts/baseline

@@ -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}"

.gitlab/scripts/rebaseline

@@ -11,7 +11,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-# There will be collision between corona and 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

INSTALL

@@ -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).

config/config.hpp

@@ -23,6 +23,9 @@
 #include "_config.hpp"
 #endif
 
+#include <cstdint>
+#include <climits>
+
 namespace mfem
 {
 

fem/bilinearform_ext.cpp

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

fem/bilinearform_ext.hpp

@@ -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;

fem/dfem/util.hpp

@@ -327,8 +327,8 @@ void print_mpi_sync(const std::string& msg)
    // First gather string lengths
    size_t msg_len = msg.length();
    std::vector<size_t> lengths(nranks);
-   MPI_Gather(&msg_len, 1, 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)

fem/fespace.hpp

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

fem/integ/bilininteg_mass_pa.cpp

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

fem/linearform_ext.cpp

@@ -51,7 +51,7 @@ void LinearFormExtension::Assemble()
       {
          // scan the attributes to set the markers to 0 or 1
          const int NE = fes.GetNE();
-         const auto attr = attributes.Read();
+         const auto attr = attributes->Read();
          const auto dimk = domain_integs_marker_k->Read();
          auto markers_w = markers.Write();
          mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
@@ -93,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,

fem/linearform_ext.hpp

@@ -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;

fem/qfunction.hpp

@@ -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;

fem/qspace.cpp

@@ -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++)

fem/qspace.hpp

@@ -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; }

general/communication.cpp

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

general/communication.hpp

@@ -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).

general/forall.hpp

@@ -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
 
 
@@ -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)
 {

linalg/CMakeLists.txt

@@ -23,6 +23,7 @@ list(APPEND SRCS
   complex_operator.cpp
   constraints.cpp
   densemat.cpp
+  functional.cpp
   symmat.cpp
   handle.cpp
   matrix.cpp
@@ -49,6 +50,7 @@ list(APPEND HDRS
   complex_operator.hpp
   constraints.hpp
   densemat.hpp
+  functional.hpp
   dinvariants.hpp
   symmat.hpp
   dtensor.hpp

linalg/functional.cpp

@@ -0,0 +1,68 @@
+#include "functional.hpp"
+
+namespace mfem
+{
+QuadraticFunctional::QuadraticFunctional(const Operator *A_,
+                                         const Vector *b_, const real_t beta_, const real_t c_)
+   : Functional(A_ ? A_->Width() : 0)
+   , A(A_), beta(beta_), b(b_), c(c_)
+   , aux(A_ ? A_->Width() : 0)
+{
+   MFEM_VERIFY(A_ != nullptr,
+               "QuadraticFunctional: A must not be nullptr. "
+               << "Use QuadraticFunctional() constructor to create an empty Quadratic functional.");
+   MFEM_VERIFY(A_->Width() == A_->Height(),
+               "QuadraticFunctional: A must be a square operator.");
+   MFEM_VERIFY(b_ == nullptr || A_->Width() == b_->Size(),
+               "QuadraticFunctional: A and b must have compatible sizes");
+}
+
+void QuadraticFunctional::SetOperator(const Operator &A_)
+{
+   MFEM_VERIFY(A_.Width() == A_.Height(),
+               "QuadraticFunctional: A must be a square operator.");
+   A = &A_;
+   width = A_.Width();
+   aux.SetSize(width);
+}
+
+void QuadraticFunctional::SetVector(const Vector &b_, const real_t beta_)
+{
+   MFEM_VERIFY(A != nullptr && A->Width() == b_.Size(),
+               "QuadraticFunctional: A and b must have compatible sizes.");
+   b = &b_;
+   beta = beta_;
+}
+void QuadraticFunctional::SetConstant(real_t c_) { c = c_; }
+
+void QuadraticFunctional::Mult(const Vector &x, Vector &y) const
+{
+   MFEM_ASSERT(A != nullptr, "QuadraticFunctional: A must not be nullptr");
+   A->Mult(x, aux); // aux = A(x)
+   if (b != nullptr) { aux.Add(2.0*beta, *b); } // aux = A(x) + 2*beta*b
+   real_t result = 0.0;
+#ifdef MFEM_USE_MPI
+   if (IsParallel())
+   {
+      result = InnerProduct(GetComm(), x, aux); // result = <A(x) + 2*beta*b, x>
+   }
+   else
+   {
+      result = InnerProduct(x, aux); // result = <A(x) + 2*beta*b, x>
+   }
+#else
+   result = InnerProduct(x, aux); // result  = <A(x) + 2*beta*b, x>
+#endif
+   y.SetSize(1);
+   y[0] = result*0.5 + c; // y = 0.5 * <A(x) + 2*beta*b, x> + c
+}
+
+void QuadraticFunctional::EvalGradient(const Vector &x,
+                                       Vector &y) const
+{
+   y.SetSize(A->Width());
+   A->Mult(x, y);
+   if (b) { y.Add(beta,*b); }
+}
+
+} // namespace mfem

linalg/functional.hpp

@@ -0,0 +1,744 @@
+#ifndef MFEM_FUNCTIONAL_HPP
+#define MFEM_FUNCTIONAL_HPP
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+#include "../general/communication.hpp"
+#endif
+
+#include "operator.hpp"
+#include "blockvector.hpp"
+#include "solvers.hpp"
+#include <cxxabi.h>
+#include <vector>
+
+namespace mfem
+{
+/// @brief A base class for functionals F:R^n->R
+///
+/// This class provides an interface for evaluating
+/// $ F:R^n->R, \nabla F:R^n->R^n $, and $ \nabla^2 F:R^n x R^n ->R^n $.
+/// F.Mult(x, y) evaluates the functional at a point x, and stores the result in y[0]
+/// F.GetGradient() returns an operator that evaluates the gradient
+/// F.GetGradient().GetGradient(x) returns an Hessian action operator.
+///
+/// The usual Operator::GetGradient(const Vector &x) method for this method
+/// is deprecated as $ \nabla F $ only takes a single argument x.
+/// It is redundant to use F.GetGradient(x).Mult(x, y) to evaluate the gradient.
+/// Instead, use F.GetGradient().Mult(x, y) to evaluate the gradient at x.
+///
+/// The gradient and Hessian can be defined in two ways:
+/// 1. If the gradient is available as a sperate operator,
+///    then override Functional::GetGradient().
+///    In this case, Functional::HessianMult() will not be called.
+///
+/// 2. Otherwise, override Functional::EvalGradient() and Functional::HessianMult()
+///    to evaluate the gradient and Hessian action, respectively.
+///    The helper classes, GradientOperator and HessianActionOperator,
+///    will call these methods to evaluate the gradient and Hessian action.
+///    If Hessian is a seperate operator, then you can override
+///    The GradientOperator::GetGradient(x) will call Functional::GetHessian(x)
+///
+class Functional : public Operator
+{
+   Operator * riesz_map = nullptr; ///< Riesz map operator, if available
+public:
+   /// @brief Create a Functional with optional gradient and hessian
+   /// @param n number of variables
+   Functional(int n=0)
+      : Operator(1, n)
+      , grad_operator(*this)
+      , hessian_action_operator(*this)
+   { }
+
+#ifdef MFEM_USE_MPI
+   Functional(MPI_Comm comm, int n=0)
+      : Functional(n)
+   { SetComm(comm); }
+
+   void SetComm(MPI_Comm comm_)
+   {
+      parallel = comm_ != MPI_COMM_NULL;
+      comm = comm_;
+   }
+   MPI_Comm GetComm() const { return comm; }
+   bool IsParallel() const { return parallel; }
+#else
+   constexpr bool IsParallel() const { return false; }
+#endif
+
+   void SetRieszMap(Operator &op) { riesz_map = &op; }
+   /// @brief return the GradientOperator that evaluates the gradient
+   ///        input x is not used. Use GetGradient().Mult(x,y) to evaluate the gradient
+   ///        we recommend using GetGradient() instead of GetGradient(x)
+   /// Deprecated. See Functional::GetGradient()
+   MFEM_DEPRECATED
+   Operator &GetGradient(const Vector &dummy) const override final { return GetGradient(); }
+
+   /// @brief Return the GradientOperator that wraps Functional::EvalGradient() for Mult().
+   /// @note If the functional has a corresponding standalone gradient operator,
+   /// override this method to return the gradient operator.
+   virtual Operator &GetGradient() const { return grad_operator; }
+
+   /// @brief Evaluate the functional at a point x that will be called by GradientOperator::Mult()
+   /// @note This method is not meant to be called directly. See, GradientOperator
+   virtual void EvalGradient(const Vector &x, Vector &y) const
+   {
+      MFEM_ABORT("Functional::EvalGradient() not implemented");
+   }
+   /// @brief Evaluate the Hessian action at a point x and direction d
+   /// that will be called by Functional::GetGradient().GetHessian(x).Mult(d,y)
+   /// @note This method is not meant to be called directly. See, HessianActionOperator
+   virtual void HessianMult(const Vector &x, const Vector &d, Vector &y) const
+   {
+      MFEM_ABORT("Functional::HessianMult() not implemented.");
+   }
+   /// @brief Return the HessianActionOperator at evaluation point x
+   ///        that wraps Functional::HessianMult() for Mult().
+   ///        See, HessianActionOperator and Functional::HessianMult().
+   ///
+   /// @note If the Hessian is available as a seperate operator, override this method.
+   ///
+   /// @warning If GetGradient() is overridden, this method will not be used.
+   virtual Operator &GetHessian(const Vector &x) const
+   {
+      hessian_action_operator.SetX(x);
+      return hessian_action_operator;
+   }
+private:
+#ifdef MFEM_USE_MPI
+   bool parallel=false;
+#else
+   const static bool parallel=false;
+#endif
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+   /// @brief A helper class to return an operator that evaluates the gradient
+   ///        using Functional::EvalGradient() method.
+   class GradientOperator : public Operator
+   {
+   private: const Functional &f; mutable Vector der;
+   public:
+      GradientOperator(const Functional &f) : Operator(f.Width()), f(f) {}
+      /// @brief Evaluate the gradient of Functional at a point x
+      void Mult(const Vector &x, Vector &y) const override final
+      {
+         if (f.riesz_map)
+         {
+            der.SetSize(f.Width());
+            f.EvalGradient(x, der);
+            f.riesz_map->Mult(der, y);
+         }
+         else
+         {
+            f.EvalGradient(x, y);
+         }
+      }
+      /// @brief Evaluate the Hessian of Functional at a point x
+      Operator &GetGradient(const Vector &x) const override final { return f.GetHessian(x); }
+   };
+   friend class GradientOperator;
+
+   /// @brief A helper class to return an operator that applies the Hessian action
+   ///        using Functional::HessianMult() method.
+   class HessianActionOperator : public Operator
+   {
+   private:
+      const Functional &f;
+      const Vector *x;
+   public:
+      HessianActionOperator(const Functional &f) : Operator(f.Width()), f(f) {}
+      void SetX(const Vector &new_x) { x = &new_x; }
+      void Mult(const Vector &d, Vector &y) const override { f.HessianMult(*x, d, y); }
+   };
+   friend class HessianActionOperator;
+
+   mutable GradientOperator grad_operator;
+   mutable HessianActionOperator hessian_action_operator;
+};
+
+/// @brief Stacked functioanl operator, [f1, ..., fk] where fi:R^n->R are functionals
+/*
+   Typical usage of this class is to provide a single operator for multiple constraints.
+   For example, consider a minimization problem with k constraints,
+   min f0(u) s.t. fi(u)=0, i=1,...,k.
+   The Lagrangian functional is
+   $ L(u, lambda) = F0(u) + sum_i lambda_i * fi(u) $
+   The first-order optimality conditions are
+   $ \nabla f0(u) + \sum_i lambda_i * grad fi(u) = 0 $
+   $ fi(u) = 0 $
+   where lambda_i are the Lagrange multipliers.
+   The StackedFunctional class can be used to represent the list of constraints fi(u).
+
+   StackedFunctional::Mult(u, y) will evaluate each functional y[i]=fi(u)
+
+   StackedFunctional::GetGradient(u) represents an operator, column-stacked gradient
+   That is, [grad f0(u), ..., grad fk(u)] in R^{n x k}
+   If you want to extract the gradient as a matrix, use
+   StackedFunctional::GetGradientMatrix(const Vector &x, DenseMatrix &grad)
+   As functionals are not assumed to return a sparse vector, the gradient is dense.
+
+   StackedFunctional::GetGradient(u).Mult(lambda, y) contract the gradients with the Lagrange multipliers
+   y = sum lambda_i * grad fi(u)
+   StackedFunctional::GetGradient(u).MultTranspose(d, y) return the directional derivative for each k
+   y[i] = <grad fi(u), d>
+
+   StackedFunctional::GetHessian(u, lambda).Mult(d, y) will return the contracted Hessian action
+   $ y = \sum_i \lambda_i * H_{fi}(u, d) $
+*/
+/// @warning Functionals should be all serial or all parallel.
+///
+class StackedFunctional : public Operator
+{
+public:
+   StackedFunctional(int n=0)
+      : Operator(0, n)
+      , funcs(0)
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   {}
+   StackedFunctional(Functional &f)
+      : Operator(0, f.Width())
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   { AddFunctional(f); }
+
+   StackedFunctional(const std::vector<Functional*> &funcs)
+      : Operator((int)funcs.size(), funcs[0]->Width())
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   { for (auto &f : funcs) { AddFunctional(*f); } }
+
+   void AddFunctional(Functional &f)
+   {
+#ifdef MFEM_USE_MPI
+      if (funcs.empty()) { if (f.IsParallel()) { SetComm(f.GetComm()); } }
+#endif
+      MFEM_VERIFY(f.Width() == Width(),
+                  "StackedFunctional::AddFunctional: Functional width does not match with the operator.");
+      MFEM_VERIFY(parallel == f.IsParallel(),
+                  "StackedFunctional::AddFunctional: Parallelism mismatch.");
+      funcs.push_back(&f);
+      height++;
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(Height());
+      Vector yview;
+      for (int i=0; i<Height(); i++)
+      {
+         yview.MakeRef(y, i, 1);
+         funcs[i]->Mult(x, yview);
+      }
+   }
+
+   Operator &GetGradient(const Vector &x) const override
+   {
+      grad_helper_op.SetX(x);
+      return grad_helper_op;
+   }
+
+   void GetGradientMatrix(const Vector &x, DenseMatrix &grads) const
+   {
+      grads.SetSize(Width(), Height());
+      Vector grad;
+      for (int i=0; i<Height(); i++)
+      {
+         grads.GetColumnReference(i, grad);
+         funcs[i]->GetGradient().Mult(x, grad);
+      }
+   }
+   Functional &GetFunctional(int i) const
+   {
+      MFEM_VERIFY(i >= 0 && i < Height(),
+                  "StackedFunctional::GetFunctional: Index out of bounds.");
+      return *funcs[i];
+   }
+
+   Operator &GetHessian(const Vector &x, const Vector &lambda) const
+   {
+      hessian_helper_op.SetX(x, lambda);
+      return hessian_helper_op;
+   }
+
+   bool parallel;
+   bool IsParallel() const { return parallel; }
+#ifdef MFEM_USE_MPI
+   void SetComm(MPI_Comm comm_)
+   {
+      parallel = comm != MPI_COMM_NULL;
+      comm = comm_;
+   }
+   MPI_Comm GetComm() const { return comm; }
+#endif
+
+protected:
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+   std::vector<Functional*> funcs;
+
+   class GradientOperator : public Operator
+   {
+   public:
+      GradientOperator(const StackedFunctional &op)
+         : Operator(op.Width(), op.Height())
+         , op(op)
+         , tmp_grad(op.Width())
+      {}
+      void SetX(const Vector &x) const { x_curr = &x; }
+      Operator &GetGradient(const Vector &lambda) const override
+      {
+         op.hessian_helper_op.SetX(*x_curr, lambda);
+         return op.hessian_helper_op;
+      }
+
+      void Mult(const Vector &lambda, Vector &y) const override
+      {
+         y.SetSize(op.Width());
+         y = 0.0;
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().Mult(*x_curr, tmp_grad);
+            y.Add(lambda[i], tmp_grad);
+         }
+      }
+      void MultTranspose(const Vector &x, Vector &y) const override
+      {
+         y.SetSize(op.Height());
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().Mult(x, tmp_grad);
+            y[i] = InnerProduct(tmp_grad, *x_curr);
+         }
+#ifdef MFEM_USE_MPI
+         if (op.IsParallel())
+         {
+            MPI_Allreduce(MPI_IN_PLACE, y.GetData(), op.Height(),
+                          MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                          op.GetComm());
+         }
+#endif
+      }
+
+   private:
+      const StackedFunctional &op;
+      mutable const Vector *x_curr;
+      mutable Vector tmp_grad;
+
+   };
+   class HessianActionOperator : public Operator
+   {
+   public:
+      HessianActionOperator(const StackedFunctional &op)
+         : Operator(op.Width()), op(op)
+      {}
+      void SetX(const Vector &x, const Vector &lambda) const { x_curr = &x; lambda_curr = &lambda; }
+
+      void Mult(const Vector &d, Vector &y) const override
+      {
+         y.SetSize(op.Width());
+         y = 0.0;
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().GetGradient(*x_curr).Mult(d, tmp_hessian);
+            y.Add((*lambda_curr)[i], tmp_hessian);
+         }
+      }
+   private:
+      const StackedFunctional &op;
+      mutable Vector tmp_hessian;
+      mutable const Vector *x_curr;
+      mutable const Vector *lambda_curr;
+   };
+   friend class GradientOperator;
+   friend class HessianActionOperator;
+   mutable GradientOperator grad_helper_op;
+   mutable HessianActionOperator hessian_helper_op;
+private:
+};
+
+class ConstrainedOptimizationProblem : public Functional
+{
+public:
+   ConstrainedOptimizationProblem(Functional &objective_,
+                                  Operator *eq_constraints_=nullptr,
+                                  Operator *ineq_constraints_=nullptr)
+      : Functional(objective_.Width())
+      , objective(objective_)
+      , eq_constraints(eq_constraints_)
+      , ineq_constraints(ineq_constraints_)
+   {
+      // Check Size
+      MFEM_VERIFY((eq_constraints == nullptr ||
+                   eq_constraints->Width() == objective.Width()),
+                  "ConstrainedFunctional: Equality constraints width does not match with the objective.");
+      MFEM_VERIFY((ineq_constraints == nullptr ||
+                   ineq_constraints->Width() == objective.Width()),
+                  "ConstrainedFunctional: Inequality constraints width does not match with the objective.");
+#ifdef MFEM_USE_MPI
+      if (objective.IsParallel()) { SetComm(objective.GetComm()); }
+#endif
+   }
+
+   Functional &GetObjective() { return objective; }
+   const Functional &GetObjective() const { return objective; }
+   Operator *GetEqualityConstraints() { return eq_constraints; }
+   const Operator *GetEqualityConstraints() const { return eq_constraints; }
+   Operator *GetInequalityConstraints() { return ineq_constraints; }
+   const Operator *GetInequalityConstraints() const { return ineq_constraints; }
+protected:
+   Functional &objective;
+   Operator *eq_constraints;
+   Operator *ineq_constraints;
+};
+
+/// @brief A Lagrangian functional for
+/// min F(u)
+/// subject to C(u) = 0
+/// That is, L(u, lambda) = F(u) + <lambda, C(u)>
+///
+/// We assume that $ F:R^n -> R $ is a functional,
+///                $ C:R^n -> R^k $ is an equality constraint operator,
+/// C should return a residual. That is,
+/// C(u) = c, then C.Mult(u, y) should return y[i] = C_i(u) - c_i.
+///
+/// C.GetGradient(u):R^k -> R^n that takes lambda and returns the contracted gradient at x
+/// C.GetGradient(u).Mult(lambda, y) returns y = sum lambda_i * grad C_i(u)
+///
+/// C's gradient should support MultTranspose method
+/// That is, C.GetGradient(u).MultTranspose(d, y) returns y[i] = <grad C_i(u), d>
+///
+class LagrangianFunctional : public ConstrainedOptimizationProblem
+{
+private:
+   mutable Vector eq_residual;
+public:
+   LagrangianFunctional(Functional &objective,
+                        Operator &eq_constraints)
+      : ConstrainedOptimizationProblem(objective, &eq_constraints)
+      , eq_residual(eq_constraints.Height())
+   {
+      width = objective.Width() + eq_constraints.Height();
+
+      offsets.SetSize(3);
+      offsets[0] = 0;
+      offsets[1] = objective.Width();
+      offsets[2] = eq_constraints.Height();
+      offsets.PartialSum();
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      y.SetSize(1);
+      y[0] = 0.0;
+      objective.Mult(u, y);
+      eq_constraints->Mult(u, eq_residual);
+      y[0] += InnerProduct(lambda, eq_residual);
+   }
+
+   void EvalGradient(const Vector &x, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      y.SetSize(Width());
+      BlockVector output_block(y, offsets);
+      Vector &opt_residual = output_block.GetBlock(0);
+      eq_residual = output_block.GetBlock(1);
+      y = 0.0;
+      // grad F(u) + \sum_i lambda_i grad C_i(u)
+      objective.GetGradient().Mult(u, opt_residual);
+      eq_constraints->GetGradient(u).AddMult(lambda, opt_residual);
+      // grad C_i(u)^T
+      eq_constraints->GetGradient(u).MultTranspose(u, eq_residual);
+   }
+
+   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
+   /// and direction d=[v, mu]
+   /// $ [H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
+   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      const BlockVector direction_block(const_cast<Vector&>(x), offsets);
+      const Vector &v = direction_block.GetBlock(0);
+      const Vector &mu = direction_block.GetBlock(1);
+
+      y.SetSize(Width());
+      BlockVector output_block(y, offsets);
+      Vector &opt_H = output_block.GetBlock(0); // Optimality Hessian
+      Vector &eq_H = output_block.GetBlock(1); // Equality Hessian
+
+      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u)
+      objective.GetGradient().GetGradient(u).Mult(v, opt_H);
+      eq_constraints->GetGradient(u).GetGradient(lambda).AddMult(v, opt_H);
+      eq_constraints->GetGradient(u).Mult(mu, eq_H);
+      // <grad C_i(u), d>
+      eq_constraints->GetGradient(u).MultTranspose(d, eq_H);
+   }
+
+protected:
+   Array<int> offsets; // offsets for [x, lambda, mu]
+};
+
+/// @brief An augmented Lagrangian functional of the form
+/// F(u) + 0.5 mu * ||C(u)||^2 + <lambda, C(u)>
+/// where F is the objective functional,
+/// C is the equality constraint operator,
+/// lambda is the Lagrange multiplier vector (initialized to zero),
+/// mu is the penalty parameter (defaults to 1.0)
+///
+/// Currently, only equality constraints are supported.
+///
+/// AugLagrangianFunctional::Update() will update the penalty and Lagrange multiplier vectors
+/// By default, lambda <- lambda + mu * C(u)
+///             mu <- mu (no update)
+class AugLagrangianFunctional : public ConstrainedOptimizationProblem
+{
+public:
+   AugLagrangianFunctional(Functional &objective_,
+                           Operator &eq_constraints_)
+      : ConstrainedOptimizationProblem(objective_, &eq_constraints_)
+      , lambda(eq_constraints_.Height())
+      , mu(1.0)
+      , eq_residual(eq_constraints_.Height())
+      , eq_dir(eq_constraints_.Height())
+   {
+      lambda = 0.0;
+   }
+
+   void SetLambda(const Vector &lambda_)
+   {
+      MFEM_VERIFY(lambda_.Size() == eq_constraints->Height(),
+                  "AugLagrangianFunctional: Lambda size does not match with the equality constraints.");
+      lambda = lambda_;
+   }
+
+   void SetPenalty(real_t mu_)
+   {
+      MFEM_VERIFY(mu_ >= 0.0,
+                  "AugLagrangianFunctional: Penalty parameter mu must be non-negative.");
+      mu = mu_;
+   }
+
+   virtual void Update(const Vector &x)
+   {
+      // Update the Lagrange multipliers
+      eq_constraints->AddMult(x, lambda, mu);
+      // Update the penalty parameter
+      // Do nothing
+   }
+   const Vector &GetLambda() const { return lambda; }
+   real_t GetPenalty() const { return mu; }
+
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(1);
+      objective.Mult(x, y);
+      eq_constraints->Mult(x, eq_residual);
+      y[0] += lambda*eq_residual;
+      y[0] += 0.5 * mu * (eq_residual*eq_residual);
+   }
+
+   void EvalGradient(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(Width());
+      // grad F(x) + \sum_i (lambda_i + mu * C_i(x)) grad C_i(x)
+      Vector curr_lambda = lambda; // store lambda + mu * C(x)
+      objective.GetGradient().Mult(x, y);
+      eq_constraints->Mult(x, eq_residual);
+      curr_lambda.Add(mu, eq_residual);
+      eq_constraints->GetGradient(x).AddMult(curr_lambda, y);
+   }
+
+   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
+   /// and direction d=[v, mu]
+   /// $ H_F(u,d) + \sum_i \lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
+   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
+   {
+      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u) <grad C_i(u), d>
+      objective.GetGradient().GetGradient(x).Mult(d, y);
+
+      Vector curr_lambda = lambda;
+      eq_constraints->Mult(x, eq_residual);
+      curr_lambda.Add(mu, eq_residual);
+
+      eq_constraints->GetGradient(x).GetGradient(curr_lambda).AddMult(d, y);
+      // eq_dir = <grad C_i(u), d>
+      eq_constraints->GetGradient(x).MultTranspose(d, eq_dir);
+      // mu_i <grad C_i(u), eq_dir>
+      eq_constraints->GetGradient(x).AddMult(eq_dir, y, mu);
+   }
+
+protected:
+   Vector lambda;
+   real_t mu;
+   mutable Vector eq_residual; // residual of the equality constraints, R^k
+   // directional derivative of the equality constraints, R^k
+   mutable Vector eq_dir;
+};
+
+/// @brief Quadratic functional of the form
+///        f(u) = 0.5 * <A u, u> + beta<b, u> + c
+///        where A is a square (possibly nonlinear) operator,
+///        beta is a scalar (defaults to 1.0, not used when b is nullptr),
+///        b is a vector (independent of u, optional),
+///        c is a constant (independent of u, optional).
+///        GetHessian() returns the operator A.
+///
+class QuadraticFunctional : public Functional
+{
+public:
+   QuadraticFunctional()
+      : Functional(0)
+      , A(nullptr), b(nullptr), c(0.0)
+   {}
+   QuadraticFunctional(const Operator *A_, const Vector *b_=nullptr,
+                       const real_t beta_=1.0,
+                       const real_t c_=0.0);
+#ifdef MFEM_USE_MPI
+   QuadraticFunctional(MPI_Comm comm_)
+      : QuadraticFunctional()
+   { SetComm(comm_); }
+
+   QuadraticFunctional(MPI_Comm comm_, const Operator *A_,
+                       const Vector *b_=nullptr,
+                       const real_t beta_=1.0, const real_t c_=0.0)
+      : QuadraticFunctional(A_, b_, beta_, c_)
+   { SetComm(comm_); }
+#endif
+
+   void SetOperator(const Operator &A_);
+   void SetVector(const Vector &b_, const real_t beta_=1.0);
+   void SetConstant(real_t c_);
+
+   void Mult(const Vector &x, Vector &y) const override;
+   void EvalGradient(const Vector &x, Vector &y) const override;
+protected:
+   const Operator *A;
+   real_t beta;
+   const Vector *b;
+   real_t c;
+   mutable Vector aux;
+protected:
+   /// @brief return the underlying Operator A
+   /// @warning Modifying the returned operator leads to undefined behavior.
+   Operator& GetHessian(const Vector &dummy) const override
+   {
+      return const_cast<Operator&>(*A);
+   }
+};
+
+class Optimizer : public IterativeSolver
+{
+public:
+   Optimizer() : IterativeSolver(), f(nullptr) { }
+#ifdef MFEM_USE_MPI
+   Optimizer(MPI_Comm comm) : IterativeSolver(comm), f(nullptr) { }
+#endif
+   // @brief Set the subproblem functional operator
+   // @param op the functional operator
+   // @note The functional will be stored in subproblem, and oper will be set to the gradient of the functional.
+   void SetOperator(const Functional &f_)
+   {
+      f = &f_;
+      IterativeSolver::SetOperator(f_.GetGradient());
+   }
+   virtual void SetLinearSolver(Solver &prec) { IterativeSolver::SetPreconditioner(prec); }
+
+   /// @brief This will abort. Should be called only with a Functional operator.
+   void SetOperator(const Operator &op) override
+   {
+      MFEM_ABORT("OptSolver::SetOperator() should not be called directly. Use SetFunctional() instead.");
+   }
+protected:
+   const Functional * f;
+};
+
+class NewtonOptimizer : public Optimizer
+{
+private:
+   real_t step_size = 1.0; // default step size
+public:
+   NewtonOptimizer() : Optimizer() { }
+#ifdef MFEM_USE_MPI
+   NewtonOptimizer(MPI_Comm comm) : Optimizer(comm) { }
+#endif
+   void SetStepSize(real_t step_size_) { step_size = step_size_; }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      dx.SetSize(x.Size());
+      y.SetSize(x.Size());
+      y = x;
+      MFEM_ASSERT(f != nullptr,
+                  "NewtonOptimizer::Mult() called without a functional operator.");
+      MFEM_ASSERT(prec != nullptr,
+                  "NewtonOptimizer::Mult() called without a linear solver.");
+      for (int i=0; i<max_iter; i++)
+      {
+         oper->Mult(y, grad);
+         Operator &hess = oper->GetGradient(y);
+         prec->SetOperator(hess);
+         prec->Mult(grad, dx);
+         y.Add(-step_size, dx);
+         if (Dot(dx, dx) < abs_tol*abs_tol)
+         {
+            break;
+         }
+      }
+   }
+private:
+   mutable Vector grad;
+   mutable Vector dx;
+};
+
+class GradientDescentOptimizer : public Optimizer
+{
+private:
+   real_t step_size = 1.0; // default step size
+public:
+   GradientDescentOptimizer() : Optimizer() { }
+#ifdef MFEM_USE_MPI
+   GradientDescentOptimizer(MPI_Comm comm) : Optimizer(comm) { }
+#endif
+   void SetStepSize(real_t step_size_) { step_size = step_size_; }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      grad.SetSize(x.Size());
+      y.SetSize(x.Size());
+      y = x;
+      MFEM_ASSERT(f != nullptr,
+                  "NewtonOptimizer::Mult() called without a functional operator.");
+      MFEM_ASSERT(prec != nullptr,
+                  "NewtonOptimizer::Mult() called without a linear solver.");
+      for (int i=0; i<max_iter; i++)
+      {
+         oper->Mult(y, grad);
+         y.Add(-step_size, grad);
+         if (Dot(grad, grad) < abs_tol*abs_tol)
+         {
+            break;
+         }
+      }
+   }
+private:
+   mutable Vector grad;
+};
+
+
+} // namespace mfem
+#endif // MFEM_FUNCTIONAL_HPP

linalg/vector.cpp

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

mesh/mesh.cpp

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

mesh/mesh.hpp

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

mesh/pmesh.cpp

@@ -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!");
+      }
    }
 }
 

mesh/pmesh.hpp

@@ -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;