Merge remote-tracking branch 'upstream/feature/partial_ess_bcs' into exaconstit-csr

fix build issue noted on github ci
make style
2022-06-07 11:03:02 -07:00 · 2022-06-07 10:45:36 -07:00 · 2022-06-07 10:34:25 -07:00 · 2022-06-07 10:14:43 -07:00 · 2022-06-07 10:10:31 -07:00 · 2022-06-02 09:49:38 -07:00
235 changed files with 18971 additions and 2211 deletions
@@ -1,10 +1,13 @@
 version: '{build}'

 # https://www.appveyor.com/docs/build-environment/#build-worker-images
-image: Visual Studio 2017
+image: Visual Studio 2019

 install:

+# Start from outside clone directory
+- cd ..
+
 # Install MS-MPI
 - ps: Start-FileDownload 'https://download.microsoft.com/download/B/2/E/B2EB83FE-98C2-4156-834A-E1711E6884FB/MSMpiSetup.exe'
 - MSMpiSetup.exe -unattend
@@ -15,6 +18,11 @@ install:
 - msmpisdk.msi /passive
 - set PATH=C:\Program Files\Microsoft MPI\Bin;%PATH%

+# Set MSMPI environment variables needed for CMake detection
+- set MSMPI_LIB32=C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x86
+- set MSMPI_LIB64=C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x64
+- set MSMPI_INC=C:\Program Files (x86)\Microsoft SDKs\MPI\Include
+
 # Install METIS, use a mirror because the original source server is not always
 # up. Original url:
 #    http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz
@@ -25,20 +33,24 @@ install:
 - cmake -H. -Bbuild
 # -DCMAKE_BUILD_TYPE=Release
 - cmake --build build
+- set METIS_PATH=%cd%
 - cd ..

 # Install hypre
 - ps: Start-FileDownload 'https://github.com/hypre-space/hypre/archive/v2.19.0.tar.gz'
 - 7z x v2.19.0.tar.gz -so | 7z x -si -ttar > nul
 - cd hypre-2.19.0/src
- cmake -H. -Bbuild -DMPI_C_INCLUDE_PATH="C:\Program Files (x86)\Microsoft SDKs\MPI\Include" -DMPI_C_LIBRARIES="C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x86\msmpi.lib" -DMPI_CXX_LIBRARIES="C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x86\msmpi.lib" -DMPI_CXX_INCLUDE_PATH="C:\Program Files (x86)\Microsoft SDKs\MPI\Include"
+- cmake -H. -Bbuild
 - cmake --build build
 - cmake --build build --target install
 - cd ../..

+# Return to clone directory
+- cd %APPVEYOR_BUILD_FOLDER%
+
 # MFEM
 before_build:
- cmake -H. -DCMAKE_INSTALL_PREFIX=install -Bbuild_parallel -DMFEM_USE_MPI=TRUE -DMFEM_USE_METIS_5=TRUE -DMPI_CXX_LIBRARIES="C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x86\msmpi.lib" -DMPI_CXX_INCLUDE_PATH="C:\Program Files (x86)\Microsoft SDKs\MPI\Include" -DHYPRE_DIR=%cd%\hypre-2.19.0\src\hypre -DMETIS_LIBRARIES=%cd%\metis-5.1.0\build\libmetis\Debug\metis.lib -DMETIS_INCLUDE_DIRS=%cd%\metis-5.1.0\include
+- cmake -H. -DCMAKE_INSTALL_PREFIX=install -Bbuild_parallel -DMFEM_USE_MPI=TRUE -DMFEM_USE_METIS_5=TRUE -DHYPRE_DIR=%cd%\..\hypre-2.19.0\src\hypre -DMETIS_LIBRARIES=%METIS_PATH%\build\libmetis\Debug\metis.lib -DMETIS_INCLUDE_DIRS=%METIS_PATH%\include
 - cmake -H. -DCMAKE_INSTALL_PREFIX=install -Bbuild_serial -DMFEM_USE_MPI=FALSE

 build_script:
@@ -0,0 +1,80 @@
+name: Build Deploy Container
+
+on:
+
+  # Always have a base image ready to go - this is a nightly build
+  schedule:
+    - cron: 0 3 * * *
+
+  # Allow manual trigger of a build
+  workflow_dispatch:
+
+  # On push to main we build and deploy images
+  push:
+    branches:
+      - master
+
+  # Publish packages on release
+  release:
+    types: [published]
+
+jobs:
+  build:
+    permissions:
+      packages: write
+    strategy:
+      fail-fast: false
+      matrix:
+
+        # Dockerfiles to build, a matrix supports future expanded builds
+        container: [["config/docker/Dockerfile", "ghcr.io/mfem/mfem-ubuntu-base"]]
+
+    runs-on: ubuntu-latest
+    name: Build
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v3
+
+      - name: Make Space For Build
+        run: |
+          sudo rm -rf /usr/share/dotnet
+          sudo rm -rf /opt/ghc
+
+       # It's easier to reference named variables than indexes of the matrix
+      - name: Set Environment
+        env:
+          dockerfile: ${{ matrix.container[0] }}
+          uri: ${{ matrix.container[1] }}
+        run: |
+          echo "dockerfile=$dockerfile" >> $GITHUB_ENV
+          echo "uri=$uri" >> $GITHUB_ENV
+
+      - name: Pull previous layers for cache
+        run: docker pull ${uri}:latest || echo "No container to pull"
+
+      - name: Build Container
+        run: |
+           container=$uri:latest
+           docker build -f ${dockerfile} -t ${container} .
+           echo "container=$container" >> $GITHUB_ENV
+
+      - name: GHCR Login
+        if: (github.event_name != 'pull_request')
+        uses: docker/login-action@v1
+        with:
+          registry: ghcr.io
+          username: ${{ github.actor }}
+          password: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Deploy
+        if: (github.event_name != 'pull_request')
+        run: |
+            docker push ${container}
+
+      - name: Tag and Push Release
+        if: (github.event_name == 'release')
+        run: |
+            tag=${GITHUB_REF#refs/tags/}
+            echo "Tagging and releasing ${uri}:${tag}"
+            docker tag ${uri}:latest ${uri}:${tag}
+            docker push ${uri}:${tag}
@@ -27,6 +27,7 @@ on:
      - master
      - next
  pull_request:
+  workflow_dispatch:

 env:
  HYPRE_ARCHIVE: v2.19.0.tar.gz
@@ -46,11 +47,18 @@ jobs:
  builds-and-tests:
    strategy:
      matrix:
-        os: [ubuntu-18.04, macos-10.15]
+        os: [ubuntu-20.04, macos-10.15, windows-2022]
        target: [dbg, opt]
        mpi: [seq, par]
-        build-system: [make]
+        build-system: [make, cmake]
        hypre-target: [int32]
+        exclude:
+          - os: ubuntu-20.04
+            build-system: cmake
+          - os: macos-10.15
+            build-system: cmake
+          - os: windows-2022
+            build-system: make
        # 'include' allows us to:
        # - Add a variable to all jobs without creating a new matrix dimension.
        #   Codecov is defined that way.
@@ -64,13 +72,15 @@ jobs:
            codecov: NO
          - target: opt
            codecov: YES
-          - os: ubuntu-18.04
+          - os: windows-2022
+            codecov: NO
+          - os: ubuntu-20.04
            target: opt
            codecov: NO
            mpi: par
            build-system: cmake
            hypre-target: int32
-          - os: ubuntu-18.04
+          - os: ubuntu-20.04
            target: opt
            codecov: NO
            mpi: par
@@ -102,13 +112,13 @@ jobs:
    # TODO: It would be nice to have only one step, e.g. with a dedicated
    # action, but I (@adrienbernede) don't see how at the moment.
    - name: get MPI (Linux)
-      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-18.04'
+      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-20.04'
      run: |
        sudo apt-get install mpich libmpich-dev
        export MAKE_CXX_FLAG="MPICXX=mpic++"

    - name: get lcov (Linux)
-      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-18.04'
+      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-20.04'
      run: |
        sudo apt-get install lcov

@@ -129,6 +139,10 @@ jobs:
        export HOMEBREW_NO_INSTALL_CLEANUP=1
        brew install lcov

+    - name: get MPI (Windows)
+      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      uses: mpi4py/setup-mpi@v1.0.3
+
    # Get Hypre through cache, or build it.
    # Install will only run on cache miss.
    - name: cache hypre
@@ -137,36 +151,67 @@ jobs:
      uses: actions/cache@v2
      with:
        path: ${{ env.HYPRE_TOP_DIR }}
-        key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.0
+        key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2

    - name: get hypre
-      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true'
-      uses: mfem/github-actions/build-hypre@v2.0
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-2022'
+      uses: mfem/github-actions/build-hypre@v2.2
      with:
        archive: ${{ env.HYPRE_ARCHIVE }}
        dir: ${{ env.HYPRE_TOP_DIR }}
        target: ${{ matrix.hypre-target }}
+        build-system: make
+
+    - name: get hypre (Windows)
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-2022'
+      uses: mfem/github-actions/build-hypre@v2.2
+      with:
+        archive: ${{ env.HYPRE_ARCHIVE }}
+        dir: ${{ env.HYPRE_TOP_DIR }}
+        target: ${{ matrix.hypre-target }}
+        build-system: cmake

    # Get Metis through cache, or build it.
    # Install will only run on cache miss.
    - name: cache metis
      id: metis-cache
-      if: matrix.mpi == 'par'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-2022'
      uses: actions/cache@v2
      with:
        path: ${{ env.METIS_TOP_DIR }}
-        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.0
+        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2

    - name: install metis
-      if: matrix.mpi == 'par' && steps.metis-cache.outputs.cache-hit != 'true'
-      uses: mfem/github-actions/build-metis@v2.0
+      if: matrix.mpi == 'par' && matrix.os != 'windows-2022' && steps.metis-cache.outputs.cache-hit != 'true'
+      uses: mfem/github-actions/build-metis@v2.2
      with:
        archive: ${{ env.METIS_ARCHIVE }}
        dir: ${{ env.METIS_TOP_DIR }}

+    - name: cache vcpkg (Windows)
+      id: vcpkg-cache
+      uses: actions/cache@v3
+      with:
+        path: vcpkg_cache
+        key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
+
+    - name: prepare binary cache location
+      if: matrix.os == 'windows-2022' && steps.vcpkg-cache.outputs.cache-hit != 'true'
+      run: |
+        mkdir -p vcpkg_cache
+
+    - name: install metis (Windows)
+      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      env:
+        VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
+      run: |
+        vcpkg install metis --triplet=x64-windows-static
+
    # MFEM build and test
    - name: build
-      uses: mfem/github-actions/build-mfem@v2.1
+      uses: mfem/github-actions/build-mfem@v2.2
+      env:
+        VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
      with:
        os: ${{ matrix.os }}
        target: ${{ matrix.target }}
@@ -176,6 +221,8 @@ jobs:
        hypre-dir: ${{ env.HYPRE_TOP_DIR }}
        metis-dir: ${{ env.METIS_TOP_DIR }}
        mfem-dir: ${{ env.MFEM_TOP_DIR }}
+        config-options: ${{ env.MFEM_EXTRA_CONFIG }}
+        library-only: ${{ matrix.target == 'dbg' }}

    # Run checks (and only checks) on debug targets
    - name: checks
@@ -193,10 +240,27 @@ jobs:
      run: |
        cd ${{ env.MFEM_TOP_DIR }} && make test

-    - name: cmake unit tests
-      if: matrix.build-system == 'cmake'
+    - name: cmake checks
+      if: matrix.build-system == 'cmake' && matrix.target == 'dbg'
      run: |
-        cd ${{ env.MFEM_TOP_DIR }}/build/tests/unit && ctest --output-on-failure
+        CTEST_CONFIG="Debug"
+        cd ${{ env.MFEM_TOP_DIR }} && cmake --build build --target check --config ${CTEST_CONFIG}
+      shell: bash
+
+    - name: cmake unit tests (Ubuntu 20.04)
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-20.04'
+      run: |
+        CTEST_CONFIG="Release"
+        [[ ${{ matrix.target }} == 'dbg' ]] && CTEST_CONFIG="Debug"
+        cd ${{ env.MFEM_TOP_DIR }}/build/tests/unit && ctest --output-on-failure -C ${CTEST_CONFIG}
+      shell: bash
+
+    - name: cmake tests
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-20.04'
+      run: |
+        CTEST_CONFIG="Release"
+        cd ${{ env.MFEM_TOP_DIR }}/build && ctest --output-on-failure -C ${CTEST_CONFIG}
+      shell: bash

    # Code coverage (process and upload reports)
    - name: codecov
@@ -20,6 +20,7 @@ on:
      - master
      - next
  pull_request:
+  workflow_dispatch:

 env:
  HYPRE_ARCHIVE: v2.19.0.tar.gz
@@ -53,11 +54,11 @@ jobs:
      uses: actions/cache@v2
      with:
        path: ${{ env.HYPRE_TOP_DIR }}
-        key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2.0
+        key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2.2

    - name: Get Hypre
      if: steps.hypre-cache.outputs.cache-hit != 'true'
-      uses: mfem/github-actions/build-hypre@v2.0
+      uses: mfem/github-actions/build-hypre@v2.2
      with:
        archive: ${{ env.HYPRE_ARCHIVE }}
        dir: ${{ env.HYPRE_TOP_DIR }}
@@ -68,18 +69,18 @@ jobs:
      uses: actions/cache@v2
      with:
        path: ${{ env.METIS_TOP_DIR }}
-        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.0
+        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2

    - name: Install Metis
      if: steps.metis-cache.outputs.cache-hit != 'true'
-      uses: mfem/github-actions/build-metis@v2.0
+      uses: mfem/github-actions/build-metis@v2.2
      with:
        archive: ${{ env.METIS_ARCHIVE }}
        dir: ${{ env.METIS_TOP_DIR }}

    # MFEM build and test
    - name: build-mfem
-      uses: mfem/github-actions/build-mfem@v2.0
+      uses: mfem/github-actions/build-mfem@v2.2
      with:
        os: ${{ runner.os }}
        target: opt
@@ -16,16 +16,28 @@ permissions:

 on:
  push:
+  pull_request:
+  workflow_dispatch:
+
+# This workflow is run on pushes to any branch in the MFEM repo (with or without
+# PRs), as well as on updates to PRs from forks. In particular, we do not
+# duplicate work by running on both pushes and updates to local PRs. We do that
+# by checking if the workflow trigger is 'push' ("github.event_name == 'push'")
+# and if we are in a fork ("github.event.pull_request.head.repo.full_name !=
+# github.repository").

 jobs:
  file-headers-check:
    runs-on: ubuntu-18.04
-
+    if: |
+      (github.event_name == 'push' ||
+       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: Cancel Previous Runs
      uses: styfle/cancel-workflow-action@0.9.0
      with:
        access_token: ${{ github.token }}
+
    - name: checkout mfem
      uses: actions/checkout@v2

@@ -49,7 +61,10 @@ jobs:
      continue-on-error: true

    - name: wrap-up
-      if: steps.copyright.outcome != 'success' || steps.license.outcome != 'success' || steps.release.outcome != 'success'
+      if: |
+        steps.copyright.outcome != 'success' ||
+        steps.license.outcome != 'success' ||
+        steps.release.outcome != 'success'
      run: |
        if [[ "${{ steps.copyright.outcome }}" != "success" ]]; then
          echo "copyright check failed, unroll log for details"
@@ -64,7 +79,9 @@ jobs:

  code-style:
    runs-on: ubuntu-18.04
-
+    if: |
+      (github.event_name == 'push' ||
+       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: checkout mfem
      uses: actions/checkout@v2
@@ -79,7 +96,9 @@ jobs:

  documentation:
    runs-on: ubuntu-18.04
-
+    if: |
+      (github.event_name == 'push' ||
+       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: checkout mfem
      uses: actions/checkout@v2
@@ -94,9 +113,12 @@ jobs:
        ./runtest documentation

  branch-history:
-    if: github.ref != 'refs/heads/next' && github.ref != 'refs/heads/master'
+    if: |
+      github.ref != 'refs/heads/next' &&
+      github.ref != 'refs/heads/master' &&
+      (github.event_name == 'push' ||
+       github.event.pull_request.head.repo.full_name != github.repository)
    runs-on: ubuntu-18.04
-
    steps:
    - name: checkout mfem
      uses: actions/checkout@v2
@@ -240,6 +240,7 @@ miniapps/navier/navier_kovasznay_vs
 miniapps/navier/navier_tgv
 miniapps/navier/navier_shear
 miniapps/navier/navier_3dfoc
+miniapps/navier/navier_turbchan
 miniapps/navier/tgv_out*.txt
 miniapps/navier/*_output

@@ -260,12 +261,14 @@ miniapps/performance/sol.*

 miniapps/shifted/distance
 miniapps/shifted/ParaViewDistance
+miniapps/shifted/ParaViewLSF
 miniapps/shifted/extrapolate
 miniapps/shifted/ParaViewExtrapolate
 miniapps/shifted/diffusion
 miniapps/shifted/diffusion.mesh
 miniapps/shifted/diffusion.gf
 miniapps/shifted/ParaViewDiffusion
+miniapps/shifted/lsf_integral

 miniapps/tools/display-basis
 miniapps/tools/load-dc
@@ -45,5 +45,5 @@ variables:
    - echo ${MFEM_DATA_DIR}
    - echo ${SPEC}
    # Next script uses 'THREADS': leaving it empty --> it uses 'make all -j'
-    - lalloc 1 -W 30 -q pdebug --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - lalloc 1 -W 45 -q pdebug --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
  needs: [setup]
@@ -52,4 +52,4 @@ variables:
    - echo ${JOBID}
    - echo ${MFEM_DATA_DIR}
    - echo ${SPEC}
-    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -t 30 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
@@ -23,7 +23,7 @@ allocate_resource:
  stage: allocate_resource
  script:
    - echo ${ALLOC_NAME}
-    - salloc --exclusive --nodes=1 --partition=mi60 --time=30 --no-shell --job-name=${ALLOC_NAME}
+    - salloc --exclusive --nodes=1 --partition=mi60 --time=45 --no-shell --job-name=${ALLOC_NAME}
  timeout: 6h
  needs: [setup]

@@ -23,7 +23,7 @@ allocate_resource:
  stage: allocate_resource
  script:
    - echo ${ALLOC_NAME}
-    - salloc --exclusive --nodes=1 --partition=pdebug --time=30 --no-shell --job-name=${ALLOC_NAME}
+    - salloc --exclusive --nodes=1 --partition=pdebug --time=45 --no-shell --job-name=${ALLOC_NAME}
  timeout: 6h

 # GitLab jobs for the Quartz machine at LLNL
@@ -8,6 +8,42 @@
                               https://mfem.org


+Version 4.4.1 (development)
+===========================
+- Added example for body-fitted volumetric and shape integration using the 
+  Algoim library.
+
+- Added WhiteGaussianNoiseDomainLFIntegrator: a LinearFormIntegrator class for
+  spatial Gaussian white noise.
+
+- Added a new Zienkiewicz-Zhu patch recovery-based a posteriori error estimator.
+  See fem/estimators.hpp.
+
+- Added support for ParMoonolith, https://bitbucket.org/zulianp/par_moonolith,
+  which provides parallel non-conforming, non-matching, variational, volumetric
+  mesh information transfer. With ParMortarAssember, fields can be exchanged
+  between arbitrarily distributed and unrelated finite element meshes in a
+  variationally consistent way.
+
+- Added full assembly and device support for several LinearForm integrators:
+  * DomainLF: (f, v)
+  * VectorDomainLF: ((f1,...,fn), (v1,...,vn))
+  * DomainLFGrad: (f, grad(v))
+  * VectorDomainLFGrad: ((f1x,f1y,f1z,...,fnx,fny,fnz), grad(v1,...,vn))
+
+- Add a new example code, Example 33/33p, to demonstrate the solution of
+  spectral fractional PDEs with MFEM.
+
+- Added a Dockerfile for a simple MFEM container, see config/docker/README.md.
+
+- Added support for assembling low-order-refined matrices using a GPU-enabled
+  "batched" algorithm. The lor_solvers and plor_solvers now fully support GPU
+  acceleration.
+
+- Added Windows 2022 CI testing with GitHub actions.
+
+- Added support for mixed meshes and pyramids in GSLIB-FindPoints.
+
 Version 4.4, released on March 21, 2022
 =======================================

@@ -128,6 +164,9 @@ Integrations, testing and documentation
 - Switched from Artistic Style (astyle) version 2.05.1 to version 3.1 for code
  formatting. See the "make style" target.

+- New benchmark for the different assembly levels inspired by the CEED
+  Bake-Off Problems, see tests/benchmarks/bench_assembly_levels.cpp.
+
 Miscellaneous
 -------------
 - Added a simple singleton class, Mpi, as a replacement for MPI_Session. New
@@ -51,7 +51,7 @@ project(mfem NONE)
 # Current version of MFEM, see also `makefile`.
 #   mfem_VERSION = (string)
 #   MFEM_VERSION = (int)   [automatically derived from mfem_VERSION]
-set(${PROJECT_NAME}_VERSION 4.4.0)
+set(${PROJECT_NAME}_VERSION 4.4.1)

 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -268,7 +268,7 @@ if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
  if(APPLE)
    # On macOS, the compiler needs additional help to find the <omp.h> header.
    # See issue #2642 for more information.
-    include_directories(${OpenMP_CXX_INCLUDE_DIRS})
+    set(OPENMP_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS})
  endif(APPLE)
 endif()

@@ -379,6 +379,21 @@ if (MFEM_USE_PUMI)
  endif()
 endif()

+# Moonolith
+if(MFEM_USE_MOONOLITH)
+  find_package(ParMoonolith REQUIRED)
+  if(ParMoonolith_FOUND)
+    get_target_property(
+      MOONOLITH_INCLUDE_DIRS ParMoonolith::par_moonolith
+      INTERFACE_INCLUDE_DIRECTORIES)
+    set(MOONOLITH_FOUND TRUE)
+    set(MOONOLITH_LIBRARIES ParMoonolith::par_moonolith)
+    message(
+      STATUS
+      "MOONOLITH_LIBRARIES=${MOONOLITH_LIBRARIES}, MOONOLITH_INCLUDE_DIRS=${MOONOLITH_INCLUDE_DIRS}")
+  endif()
+endif()
+
 # HiOp optimizer
 if (MFEM_USE_HIOP)
  find_package(HIOP REQUIRED)
@@ -416,6 +431,11 @@ if (MFEM_USE_CALIPER)
  find_package(Caliper REQUIRED)
 endif()

+# Algoim
+if (MFEM_USE_ALGOIM)
+  find_package(Algoim REQUIRED)
+endif()
+
 # ADIOS2 for parallel I/O
 if (MFEM_USE_ADIOS2)
   find_package(ADIOS2 REQUIRED)
@@ -459,7 +479,7 @@ set(MFEM_TPLS OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS
    PETSC SLEPC MESQUITE MUMPS STRUMPACK AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
    ADIOS2 CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
-    MPI_CXX HIP HIPSPARSE)
+    MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM)

 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
@@ -474,7 +494,6 @@ endforeach(TPL)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_INCLUDE_DIRS)
 # message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")
-include_directories(${TPL_INCLUDE_DIRS})

 if (OPENMP_FOUND)
  message(STATUS "MFEM: using package OpenMP")
@@ -496,6 +515,11 @@ message(STATUS "MFEM git string: ${MFEM_GIT_STRING}")
 set(SOURCES "")
 set(HEADERS "")
 set(MFEM_SOURCE_DIRS general linalg mesh fem)
+
+if(MFEM_USE_MOONOLITH)
+  set(MFEM_SOURCE_DIRS ${MFEM_SOURCE_DIRS} fem/moonolith)
+endif()
+
 foreach(DIR IN LISTS MFEM_SOURCE_DIRS)
  add_subdirectory(${DIR})
 endforeach()
@@ -526,6 +550,15 @@ target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES})
 if (MINGW)
  target_link_libraries(mfem PRIVATE ws2_32)
 endif()
+if (MSVC)
+  target_compile_options(mfem PUBLIC "/wd4819")
+endif()
+message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")
+target_include_directories(mfem
+  PUBLIC
+    ${TPL_INCLUDE_DIRS}
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>)
 set_target_properties(mfem PROPERTIES VERSION "${mfem_VERSION}")
 set_target_properties(mfem PROPERTIES SOVERSION "${mfem_VERSION}")

@@ -654,6 +687,10 @@ set(INSTALL_LIB_DIR lib
 set(INSTALL_CMAKE_DIR lib/cmake/mfem
  CACHE PATH "Relative path for installing cmake config files.")

+target_include_directories(mfem
+  PUBLIC
+    $<INSTALL_INTERFACE:${INSTALL_INCLUDE_DIR}>)
+
 # The 'install' target will not depend on 'all'.
 # set(CMAKE_SKIP_INSTALL_ALL_DEPENDENCY TRUE)

@@ -687,7 +724,7 @@ endif()
 if (MFEM_USE_CEED)
  install(DIRECTORY ${MFEM_SOURCE_DIRS}
    DESTINATION ${INSTALL_INCLUDE_DIR}/mfem
-    FILES_MATCHING PATTERN "fem/ceed/*.h")
+    FILES_MATCHING PATTERN "fem/ceed/integrators/*/*.h")
 endif()

 # Install ${HEADERS}
@@ -718,14 +755,8 @@ export(TARGETS ${PROJECT_NAME}
 # TODO: How do we register the install-tree? Replacing the build-tree?
 export(PACKAGE ${PROJECT_NAME})

-# Extract the include directories required to use MFEM
-get_target_property(MFEM_TPL_INCLUDE_DIRS mfem INCLUDE_DIRECTORIES)
-if (NOT MFEM_TPL_INCLUDE_DIRS)
-  set(MFEM_TPL_INCLUDE_DIRS "")
-endif()
-
 # This is the build-tree version
-set(INCLUDE_INSTALL_DIRS ${PROJECT_BINARY_DIR} ${MFEM_TPL_INCLUDE_DIRS})
+set(INCLUDE_INSTALL_DIRS ${PROJECT_BINARY_DIR} ${TPL_INCLUDE_DIRS})
 set(LIB_INSTALL_DIR ${PROJECT_BINARY_DIR})
 configure_package_config_file(config/cmake/MFEMConfig.cmake.in
  ${CMAKE_CURRENT_BINARY_DIR}/MFEMConfig.cmake
@@ -733,7 +764,7 @@ configure_package_config_file(config/cmake/MFEMConfig.cmake.in
  PATH_VARS INCLUDE_INSTALL_DIRS LIB_INSTALL_DIR)

 # This is the version that will be installed
-set(INCLUDE_INSTALL_DIRS ${INSTALL_INCLUDE_DIR}  ${MFEM_TPL_INCLUDE_DIRS})
+set(INCLUDE_INSTALL_DIRS ${INSTALL_INCLUDE_DIR}  ${TPL_INCLUDE_DIRS})
 set(LIB_INSTALL_DIR ${INSTALL_LIB_DIR})
 configure_package_config_file(config/cmake/MFEMConfig.cmake.in
  ${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/MFEMConfig.cmake
@@ -119,6 +119,7 @@ The MFEM source code has the following structure:
  │   ├── ceed
  │   ├── fe
  │   ├── qinterp
+  │   ├── moonolith
  │   └── tmop
  ├── general
  ├── linalg
@@ -471,6 +471,14 @@ MFEM_USE_CODIPACK = YES/NO
  Enable automatic differentiation using the CoDiPack library.
  www.scicomp.uni-kl.de/codi/

+MFEM_USE_ALGOIM = YES/NO
+  Enable the usage of Algoim - a collection of high-order accurate numerical 
+  methods and C++ algorithms for working with implicitly-defined geometry and 
+  level set methods. The Algoim library requires the Blitz++ library. The MFEM 
+  provides interface to Algoim v1. Thus, to check out the specific state use: 
+  git checkout 9c9ca0ef094d8ab0390ed36367a1151b459bbe0a
+  https://algoim.github.io
+
 MFEM_USE_ADFORWARD = YES/NO
  Enable forward mode for AD packages. This option is valid
  only if the AD package supports two modes (backward/forward).
@@ -520,6 +528,16 @@ MFEM_USE_MKL_CPARDISO = YES/NO
   MFEM_USE_LAPACK=YES, verify that the MKL LAPACK libraries are used. The
   OpenMP capabilities are disabled at link time.

+MFEM_USE_MOONOLITH = YES/NO
+   Enables the ParMoonolith interface for parallel non-conforming, non-matching,
+   variational, volumetric mesh information transfer. It requires the variable
+   MOONOLITH_DIR=<path to installation> to be defined in the environment in
+   order to be used with the Makefile. Makefile users are also required to
+   install moonolith using the command `make install_all`, see
+   https://bitbucket.org/zulianp/par_moonolith for details.
+   Although Moonolith is an MPI-based library, both serial (MFEM_USE_MPI=NO) and
+   parallel (MFEM_USE_MPI=YES) versions of MFEM are supported.
+
 MFEM_USE_CALIPER = YES/NO
   Enables the interface to Caliper. Caliper is a library to integrate
   performance profiling capabilities into applications. To use Caliper,
@@ -728,6 +746,22 @@ The specific libraries and their options are:
  Options: GSLIB_OPT, GSLIB_LIB.
  Versions: GSLIB >= 1.0.7.

+- ALGOIM (optional), used when MFE_USE_ALGOIM=YES. The library provides only
+  headers so it just needs to be downloaded at the same level as MFEM. Download
+  the specific version we use as:
+  "git clone https://github.com/algoim/algoim.git;
+   git checkout 9c9ca0ef094d8ab0390ed36367a1151b459bbe0a"
+  ALGOIM depends on BLITZ and rhe library must be built prior to the MFEM build.
+  Download v1.0.2, untar it at the same level as MFEM and create a symbolic link:
+  "ln -s blitz-1.0.2 blitz".
+  Build Blitz using CMake as:
+  "cmake . -DCMAKE_INSTALL_PREFIX=.; make lib; make install"
+  URL: https://github.com/blitzpp/blitz/archive/refs/tags/1.0.2.tar.gz
+  Options: BLITZ_OPT, BLITZ_LIB
+  Versions: BLITZ = 1.0.2
+
+
+
 - MKL CPardiso (optional), used when MFEM_USE_MKL_CPARDISO = YES.
  URL: https://software.intel.com/content/www/us/en/develop/tools/math-kernel-library.html
  Options: MKL_CPARDISO_OPT, MKL_CPARDISO_LIB.
@@ -759,6 +793,11 @@ The specific libraries and their options are:
  Options: RAJA_DIR, RAJA_OPT, RAJA_LIB.
  Versions: RAJA >= 0.14.0.

+- Moonolith (optional), use when MFEM_USE_MOONOLITH = YES.
+  URL: https://bitbucket.org/zulianp/par_moonolith
+  Options: MOONOLITH_DIR
+  Versions: MOONOLITH >= 1.1.0.
+
 - Caliper (optional), used when MFEM_USE_CALIPER = YES.
  URL: https://github.com/LLNL/Caliper
  Options: CALIPER_DIR
@@ -932,6 +971,7 @@ MFEM_USE_CEED
 MFEM_USE_RAJA
 MFEM_USE_UMPIRE
 MFEM_USE_SIDRE
+MFEM_USE_MOONOLITH
 MFEM_USE_CALIPER
 MFEM_USE_FMS
 MFEM_USE_BENCHMARK
@@ -990,6 +1030,7 @@ The CMake build system adds auto-detection for the following packages/libraries:
 - RAJA
 - UMPIRE
 - AXOM - Used when MFEM_USE_SIDRE is enabled
+ - MOONOLITH
 - CALIPER
 - FMS
 - BENCHMARK
@@ -7,21 +7,24 @@

                               https://mfem.org

-MFEM is a modular parallel C++ library for finite element methods. Its goal is
-to enable high-performance scalable finite element discretization research and
-application development on a wide variety of platforms, ranging from laptops to
-supercomputers.
+[MFEM](https://mfem.org) is a modular parallel C++ library for finite element
+methods. Its goal is to enable high-performance scalable finite element
+discretization research and application development on a wide variety of
+platforms, ranging from laptops to supercomputers.

 We welcome contributions and feedback from the community. Please see the file
-CONTRIBUTING.md for additional details about our development process.
+[CONTRIBUTING.md](CONTRIBUTING.md) for additional details about our development
+process.

-* For building instructions, see the file INSTALL, or type "make help".
+* For building instructions, see the file [INSTALL](INSTALL), or type "make help".

-* Copyright and licensing information can be found in files LICENSE and NOTICE.
+* Copyright and licensing information can be found in files [LICENSE](LICENSE) and [NOTICE](NOTICE).

 * The best starting point for new users interested in MFEM's features is to
  review the examples and miniapps at https://mfem.org/examples.

+* Instructions for learning with Docker are in [config/docker](config/docker).
+
 Conceptually, MFEM can be viewed as a finite element toolbox that provides the
 building blocks for developing finite element algorithms in a manner similar to
 that of MATLAB for linear algebra methods. In particular, MFEM provides support
@@ -58,12 +61,16 @@ solvers from the hypre library. Comprehensive support for other external
 packages, e.g. PETSc, SUNDIALS and libCEED is also included, giving access to
 additional linear and nonlinear solvers, preconditioners, time integrators, etc.

-For examples of using MFEM, see the examples/ and miniapps/ directories, as well
-as the OpenGL visualization tool GLVis which is available at https://glvis.org.
+For examples of using MFEM, see the [examples/](examples) and [miniapps/](miniapps)
+directories, as well as the OpenGL visualization tool GLVis which is available
+at https://glvis.org.
+
+## License

 MFEM is distributed under the terms of the BSD-3 license. All new contributions
-must be made under this license. See LICENSE and NOTICE for details.
+must be made under this license. See [LICENSE](LICENSE) and [NOTICE](NOTICE) for
+details.

-SPDX-License-Identifier: BSD-3-Clause
-LLNL Release Number: LLNL-CODE-806117
+SPDX-License-Identifier: BSD-3-Clause <br>
+LLNL Release Number: LLNL-CODE-806117 <br>
 DOI: 10.11578/dc.20171025.1248
@@ -54,9 +54,11 @@ set(MFEM_USE_CEED @MFEM_USE_CEED@)
 set(MFEM_USE_UMPIRE @MFEM_USE_UMPIRE@)
 set(MFEM_USE_SIMD @MFEM_USE_SIMD@)
 set(MFEM_USE_ADIOS2 @MFEM_USE_ADIOS2@)
+set(MFEM_USE_MOONOLITH @MFEM_USE_MOONOLITH@)
 set(MFEM_USE_CODIPACK @MFEM_USE_CODIPACK@)
 set(MFEM_USE_ADFORWARD @MFEM_USE_ADFORWARD@)
 set(MFEM_USE_CALIPER @MFEM_USE_CALIPER@)
+set(MFEM_USE_ALGOIM @MFEM_USE_ALGOIM@)
 set(MFEM_USE_BENCHMARK @MFEM_USE_BENCHMARK@)
 set(MFEM_USE_PARELAG @MFEM_USE_PARELAG@)

@@ -128,6 +128,9 @@
 // Enable MFEM functionality based on the PUMI library
 #cmakedefine MFEM_USE_PUMI

+// Enable MFEM functionality based on the Moonolith library
+#cmakedefine MFEM_USE_MOONOLITH
+
 // Enable MFEM functionality based on the HiOp library
 #cmakedefine MFEM_USE_HIOP

@@ -157,6 +160,9 @@
 // Enable MFEM functionality based on the Caliper library
 #cmakedefine MFEM_USE_CALIPER

+// Enable MFEM functionality based on the Algoim library
+#cmakedefine MFEM_USE_ALGOIM
+
 // Which library functions to use in class StopWatch for measuring time.
 // For a list of the available options, see INSTALL.
 // If not defined, an option is selected automatically.
@@ -0,0 +1,22 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - ALGOIM_FOUND
+#   - ALGOIM_LIBRARIES
+#   - ALGOIM_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(Algoim ALGOIM ALGOIM_DIR
+       	"include" "algoim_quad.hpp"
+       	"" ""
+	"Paths to headers required by Algoim."
+	"Libraries required by Algoim.")
@@ -0,0 +1,22 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - BLITZ_FOUND
+#   - BLITZ_LIBRARIES
+#   - BLITZ_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(Blitz BLITZ BLITZ_DIR
+       	"include" "blitz/blitz.h"
+       	"lib" "blitz"
+  	"Paths to headers required by Blitz."
+       	"Libraries required by Blitz.")
@@ -15,5 +15,5 @@
 #   - GSLIB_INCLUDE_DIRS

 include(MfemCmakeUtilities)
-mfem_find_package(gslib GSLIB GSLIB_DIR "include" gslib.h "lib" gs
+mfem_find_package(GSLIB GSLIB GSLIB_DIR "include" gslib.h "lib" gs
  "Paths to headers required by GSLIB." "Libraries required by GSLIB.")
@@ -893,7 +893,8 @@ function(mfem_export_mk_files)
      MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA
      MFEM_USE_HIP MFEM_USE_RAJA MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER
      MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO
-      MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG)
+      MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG
+      MFEM_USE_MOONOLITH MFEM_USE_ALGOIM)
  foreach(var ${CONFIG_MK_BOOL_VARS})
    if (${var})
      set(${var} YES)
@@ -977,9 +978,11 @@ function(mfem_export_mk_files)
      string(REGEX REPLACE "^SCOREC::" "" libname ${pumilib})
      string(FIND "${pumilib}" ".a" staticlib)
      string(FIND "${pumilib}" ".so" sharedlib)
+      string(FIND "${pumilib}" ".dylib" dynamiclib)
      find_library(lib ${libname} PATHS ${PUMI_DIR}/lib NO_DEFAULT_PATH)
      if (NOT "${sharedlib}" MATCHES "-1" OR
-          NOT "${staticlib}" MATCHES "-1"   )
+          NOT "${staticlib}" MATCHES "-1" OR
+          NOT "${dynamiclib}" MATCHES "-1"  )
        set(MFEM_EXT_LIBS "${pumilib} ${MFEM_EXT_LIBS}")
      elseif (NOT "${lib}" MATCHES "lib-NOTFOUND")
        set(MFEM_EXT_LIBS "${lib} ${MFEM_EXT_LIBS}")
@@ -994,7 +997,7 @@ function(mfem_export_mk_files)
  foreach(lib ${TPL_LIBRARIES})
    get_filename_component(suffix ${lib} EXT)
    # handle interfaces (e.g., SCOREC::apf)
-    if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*")
+    if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*" OR "${lib}" MATCHES "ParMoonolith::.*")
    elseif (TARGET "${lib}")
      mfem_get_target_options(${lib} CompileOpts LinkOpts)
      # Removing duplicates may lead to issues:
@@ -138,6 +138,9 @@
 // Enable MFEM functionality based on the PUMI library
 // #define MFEM_USE_PUMI

+// Enable Moonolith-based general interpolation between finite element spaces.
+// #define MFEM_USE_MOONOLITH
+
 // Enable MFEM functionality based on the HIOP library.
 // #define MFEM_USE_HIOP

@@ -164,6 +167,9 @@
 // Enable functionality based on the Caliper library.
 // #define MFEM_USE_CALIPER

+// Enable functionality based on the Algoim library.
+// #define MFEM_USE_ALGOIM
+
 // Enable functionality based on the Umpire library.
 // #define MFEM_USE_UMPIRE

@@ -58,6 +58,7 @@ MFEM_USE_UMPIRE        = @MFEM_USE_UMPIRE@
 MFEM_USE_SIMD          = @MFEM_USE_SIMD@
 MFEM_USE_ADIOS2        = @MFEM_USE_ADIOS2@
 MFEM_USE_MKL_CPARDISO  = @MFEM_USE_MKL_CPARDISO@
+MFEM_USE_MOONOLITH     = @MFEM_USE_MOONOLITH@
 MFEM_USE_ADFORWARD     = @MFEM_USE_ADFORWARD@
 MFEM_USE_CODIPACK      = @MFEM_USE_CODIPACK@
 MFEM_USE_BENCHMARK     = @MFEM_USE_BENCHMARK@
@@ -58,6 +58,7 @@ option(MFEM_USE_UMPIRE "Enable Umpire" OFF)
 option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" OFF)
 option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)
 option(MFEM_USE_CALIPER "Enable Caliper support" OFF)
+option(MFEM_USE_ALGOIM "Enable Algoim support" OFF)
 option(MFEM_USE_MKL_CPARDISO "Enable MKL CPardiso" OFF)
 option(MFEM_USE_ADFORWARD "Enable forward mode for AD" OFF)
 option(MFEM_USE_CODIPACK "Enable automatic differentiation (AD) using CoDiPack" OFF)
@@ -236,6 +237,11 @@ set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
 set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
 set(UMPIRE_DIR "${MFEM_DIR}/../umpire" CACHE PATH "Path to Umpire")
 set(CALIPER_DIR "${MFEM_DIR}/../caliper" CACHE PATH "Path to Caliper")
+set(BLITZ_DIR "${MFEM_DIR}/../blitz" CACHE PATH "Path to Blitz")
+set(ALGOIM_DIR "${MFEM_DIR}/../algoim" CACHE PATH "Path to Algoim")
+set(ALGOIM_REQUIRED_PACKAGES "BLITZ" CACHE STRING
+    "Packages that ALGOIM depends on.")
+
 set(BENCHMARK_DIR "${MFEM_DIR}/../google-benchmark" CACHE PATH
    "Path to Google Benchmark")

@@ -153,10 +153,12 @@ MFEM_USE_RAJA          = NO
 MFEM_USE_OCCA          = NO
 MFEM_USE_CEED          = NO
 MFEM_USE_CALIPER       = NO
+MFEM_USE_ALGOIM        = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
 MFEM_USE_MKL_CPARDISO  = NO
+MFEM_USE_MOONOLITH     = NO
 MFEM_USE_ADFORWARD     = NO
 MFEM_USE_CODIPACK      = NO
 MFEM_USE_BENCHMARK     = NO
@@ -387,6 +389,11 @@ ifeq ($(SLEPC_FOUND),YES)
      $(subst $(CXX_XLINKER),$(XLINKER),$(SLEPC_DEP))
 endif

+ifeq ($(MFEM_USE_MOONOLITH),YES)
+  include $(MOONOLITH_DIR)/config/moonolith-config.makefile
+  MOONOLITH_LIB=$(MOONOLITH_LIBRARIES)
+endif
+
 # MPFR library configuration
 MPFR_OPT =
 MPFR_LIB = -lmpfr
@@ -463,6 +470,16 @@ CALIPER_DIR = @MFEM_DIR@/../caliper
 CALIPER_OPT = -I$(CALIPER_DIR)/include
 CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib64 -lcaliper

+# BLITZ library configuration
+BLITZ_DIR = @MFEM_DIR@/../blitz
+BLITZ_OPT = -I$(BLITZ_DIR)/include
+BLITZ_LIB = $(XLINKER)-rpath,$(BLITZ_DIR)/lib -L$(BLITZ_DIR)/lib -lblitz
+
+# ALGOIM library configuration
+ALGOIM_DIR = @MFEM_DIR@/../algoim
+ALGOIM_OPT = -I$(ALGOIM_DIR)/src $(BLITZ_OPT)
+ALGOIM_LIB = $(BLITZ_LIB)
+
 # BENCHMARK library configuration
 BENCHMARK_DIR = @MFEM_DIR@/../google-benchmark
 BENCHMARK_OPT = -I$(BENCHMARK_DIR)/include
@@ -0,0 +1,30 @@
+FROM ghcr.io/rse-ops/cuda-ubuntu-20.04:cuda-11.0.3
+
+# docker build -t ghcr.io/mfem/mfem-ubuntu-base .
+
+RUN apt-get update && \
+    apt-get install -y unzip gfortran && \
+    spack compiler find && \
+    apt-get install -y libcurl4-openssl-dev libssl-dev
+
+# /code is the working directory for code
+WORKDIR /code
+COPY . /code
+
+# This is for a spack environment/view to install from there
+WORKDIR /opt/mfem-env
+RUN . /opt/spack/share/spack/setup-env.sh && \
+    spack env create -d . && \
+    echo "  concretization: together" >> spack.yaml && \
+    spack env activate . && \
+    spack develop --path /code mfem@master+examples+miniapps && \
+    spack add mfem@master+examples+miniapps && \
+    spack install
+
+# ensure mfem always on various paths
+RUN cd /opt/mfem-env && \
+    spack env activate --sh -d . >> /etc/profile.d/z10_spack_environment.sh
+
+# The user will see the view on shell into the container
+WORKDIR /opt/mfem-env/.spack-env/view/
+ENTRYPOINT ["/bin/bash", "--rcfile", "/etc/profile", "-l", "-c"]
@@ -0,0 +1,130 @@
+# mfem Docker
+
+We provide a [Dockerfile](Dockerfile) to build an ubuntu base image. You can use
+this image for a demo of using mfem! 🎉️
+
+Updated containers are built and deployed on merges to the main branch and releases.
+If you want to request a build on demand, you can [manually run the workflow](https://docs.github.com/en/actions/managing-workflow-runs/manually-running-a-workflow) thanks to the workflow dispatch event.
+
+### Usage
+
+Here is how to build the container. Note that we build so it belongs to the same
+namespace as the repository here. "ghcr.io" means "GitHub Container Registry" and
+is the [GitHub packages](https://github.com/features/packages) registry that supports
+ Docker images and other OCI artifacts. From the root of the repository:
+
+```bash
+$ docker build -f config/docker/Dockerfile -t ghcr.io/mfem/mfem-ubuntu-base .
+```
+
+or this directory:
+
+```bash
+$ docker build -f Dockerfile -t ghcr.io/mfem/mfem-ubuntu-base ../../
+```
+
+### Shell
+
+To shell into a container (here is an example with ubuntu):
+
+```bash
+$ docker run -it ghcr.io/mfem/mfem-ubuntu-base bash
+```
+
+Off the bat, you can see mfem libraries are in your path so you can jump into development:
+
+```bash
+env | grep mfem
+```
+```bash
+PKG_CONFIG_PATH=/opt/mfem-env/.spack-env/view/lib/pkgconfig:/opt/mfem-env/.spack-env/view/share/pkgconfig:/opt/mfem-env/.spack-env/view/lib64/pkgconfig
+PWD=/opt/mfem-env
+MANPATH=/opt/mfem-env/.spack-env/view/share/man:/opt/mfem-env/.spack-env/view/man:
+CMAKE_PREFIX_PATH=/opt/mfem-env/.spack-env/view
+SPACK_ENV=/opt/mfem-env
+ACLOCAL_PATH=/opt/mfem-env/.spack-env/view/share/aclocal
+LD_LIBRARY_PATH=/opt/mfem-env/.spack-env/view/lib:/opt/mfem-env/.spack-env/view/lib64
+PATH=/opt/mfem-env/.spack-env/view/bin:/opt/view/bin:/opt/spack/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin
+```
+
+#### Examples and MiniApps
+
+If you want to develop a tool that _uses_ mfem, you can find the built libraries in:
+
+```
+$ ls /opt/mfem-env/.spack-env/view/
+bin  etc  include  lib  libexec  sbin  share  var
+```
+
+And yes, this is the working directory when you shell into the container!
+You can find the examples here:
+
+
+```bash
+cd share/mfem/examples
+```
+```bash
+$ ./ex0
+Options used:
+   --mesh ../data/star.mesh
+   --order 1
+Number of unknowns: 101
+   Iteration :   0  (B r, r) = 0.184259
+   Iteration :   1  (B r, r) = 0.102754
+   Iteration :   2  (B r, r) = 0.00558141
+   Iteration :   3  (B r, r) = 1.5247e-05
+   Iteration :   4  (B r, r) = 1.13807e-07
+   Iteration :   5  (B r, r) = 6.27231e-09
+   Iteration :   6  (B r, r) = 3.76268e-11
+   Iteration :   7  (B r, r) = 6.07423e-13
+   Iteration :   8  (B r, r) = 4.10615e-15
+Average reduction factor = 0.140201
+```
+
+Try running a few, and look at the associated .cpp file for the source code!
+You can also explore the "mini apps," also in share/mfem, but under miniapps.
+
+```bash
+# This is run from the examples directory
+$ cd ../miniapps
+```
+```bash
+$ ls
+CMakeLists.txt  common            meshing  nurbs        shifted  toys
+adjoint         electromagnetics  mtop     parelag      solvers
+autodiff        gslib             navier   performance  tools
+```
+
+And an example in "toys"
+
+```bash
+cd toys
+```
+```bash
+$ ./automata -no-vis
+Options used:
+   --num-steps 16
+   --rule 90
+   --no-visualization
+
+Rule:
+ 111 110 101 100 011 010 001 000
+  0   1   0   1   1   0   1   0
+
+Applying rule...done.
+```
+
+Have fun!
+
+
+#### Your own App
+If you want to develop with your own code base
+(and mfem as is in the container) you can bind to somewhere else in the container (e.g., src)
+
+```bash
+$ docker run -it ghcr.io/mfem/mfem-ubuntu-base -v $PWD:/src bash
+```
+
+In the above, we can pretend your project is in the present working directory (PWD) and we are
+binding to source. You can then use the mfem in the container for development, and if you
+want to distribute your library or app in a container, you can use the mfem container as the base.
@@ -0,0 +1,86 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+8
+1 5  0   1   10  9   3   4   13  12 
+1 5  1   2   11  10  4   5   14  13 
+1 5  3   4   13  12  6   7   16  15
+1 5  4   5   14  13  7   8   17  16
+1 5  9   10  19  18  12  13  22  21
+1 5  10  11  20  19  13  14  23  22
+1 5  12  13  22  21  15  16  25  24
+1 5  13  14  23  22  16  17  26  25
+
+#
+
+boundary
+24
+1 3 1  0  9  10 
+1 3 2  1  10 11 
+1 3 10 9  18 19 
+1 3 11 10 19 20 
+2 3 0  3  12 9
+2 3 9  12 21 18
+2 3 3  6  15 12
+2 3 12 15 24 21
+3 3 0  1  4  3
+3 3 1  2  5  4
+3 3 3  4  7  6
+3 3 4  5  8  7 
+4 3 6  7  16 15
+4 3 7  8  17 16
+4 3 15 16 25 24
+4 3 16 17 26 25 
+5 3 18 21 22 19 
+5 3 19 22 23 20 
+5 3 21 24 25 22 
+5 3 22 25 26 23 
+6 3 2  11 14 5
+6 3 11 20 23 14
+6 3 5  14 17 8 
+6 3 14 23 26 17
+
+vertices
+27
+3
+0.0 0.0 0.0
+0.5 0.0 0.0
+1.0 0.0 0.0
+0.0 0.0 0.5
+0.5 0.0 0.5
+1.0 0.0 0.5
+0.0 0.0 1.0
+0.5 0.0 1.0
+1.0 0.0 1.0
+0.0 0.5 0.0
+0.5 0.5 0.0
+1.0 0.5 0.0
+0.0 0.5 0.5
+0.5 0.5 0.5
+1.0 0.5 0.5
+0.0 0.5 1.0
+0.5 0.5 1.0
+1.0 0.5 1.0
+0.0 1.0 0.0
+0.5 1.0 0.0
+1.0 1.0 0.0
+0.0 1.0 0.5
+0.5 1.0 0.5
+1.0 1.0 0.5
+0.0 1.0 1.0
+0.5 1.0 1.0
+1.0 1.0 1.0
@@ -0,0 +1,84 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+8
+1 5 0  1  4  3  9  10 13 12 
+1 5 1  2  5  4  10 11 14 13
+1 5 9  10 13 12 18 19 22 21
+1 5 10 11 14 13 19 20 23 22
+1 5 3  4  7  6  12 13 16 15
+1 5 4  5  8  7  13 14 17 16
+1 5 12 13 16 15 21 22 25 24
+1 5 13 14 17 16 22 23 26 25
+
+boundary
+24
+1 3 0  1  10 9 
+1 3 1  2  11 10 
+1 3 9  10 19 18 
+1 3 10 11 20 19 
+3 3 0  3  4  1 
+3 3 1  4  5  2 
+3 3 3  6  7  4 
+3 3 4  7  8  5 
+3 3 18 19 22 21 
+3 3 19 20 23 22 
+3 3 21 22 25 24 
+3 3 22 23 26 25 
+3 3 2  5  14 11 
+3 3 11 14 23 20 
+3 3 5  8  17 14 
+3 3 14 17 26 23 
+3 3 0  9  12 3 
+3 3 9  18 21 12 
+3 3 3  12 15 6 
+3 3 12 21 24 15 
+2 3 6  15 16 7 
+2 3 7  16 17 8 
+2 3 15 24 25 16 
+2 3 16 25 26 17 
+
+vertices
+27
+3
+0.0 0.0 0.0
+0.5 0.0 0.0
+1.0 0.0 0.0
+0.0 0.0 0.5
+0.5 0.0 0.5
+1.0 0.0 0.5
+0.0 0.0 1.0
+0.5 0.0 1.0
+1.0 0.0 1.0
+0.0 0.5 0.0
+0.5 0.5 0.0
+1.0 0.5 0.0
+0.0 0.5 0.5
+0.5 0.5 0.5
+1.0 0.5 0.5
+0.0 0.5 1.0
+0.5 0.5 1.0
+1.0 0.5 1.0
+0.0 1.0 0.0
+0.5 1.0 0.0
+1.0 1.0 0.0
+0.0 1.0 0.5
+0.5 1.0 0.5
+1.0 1.0 0.5
+0.0 1.0 1.0
+0.5 1.0 1.0
+1.0 1.0 1.0
@@ -0,0 +1,56 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geomety Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+
+dimension
+3
+
+elements
+4
+1 5 0 1 4 3 6 7 10 9
+1 6 4 1 5 10 7 11
+1 7 11 7 1 5 8
+1 4 2 5 1 8
+
+boundary
+14
+1 3 0 3 4 1
+1 3 6 7 10 9
+1 3 0 6 9 3
+1 3 0 1 7 6
+1 3 3 9 10 4
+1 2 1 5 2
+1 2 1 4 5
+1 2 1 8 7
+1 2 1 2 8
+1 2 2 5 8
+1 2 5 11 8
+1 2 7 8 11
+1 2 7 11 10
+1 3 4 10 11 5
+
+vertices
+12
+3
+0 0 0
+1 0 0
+2 0 0
+0 1 0
+1 1 0
+2 1 0
+0 0 1
+1 0 1
+2 0 1
+0 1 1
+1 1 1
+2 1 1
@@ -38,7 +38,7 @@ PROJECT_NAME           = "MFEM"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.

-PROJECT_NUMBER         = v4.4.0
+PROJECT_NUMBER         = v4.4.1

 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
@@ -765,11 +765,14 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                         @MFEM_SOURCE_DIR@/linalg \
                         @MFEM_SOURCE_DIR@/mesh \
                         @MFEM_SOURCE_DIR@/fem \
+                         @MFEM_SOURCE_DIR@/fem/moonolith \
                         @MFEM_SOURCE_DIR@/fem/fe \
+                         @MFEM_SOURCE_DIR@/fem/lor \
                         @MFEM_SOURCE_DIR@/examples \
                         @MFEM_SOURCE_DIR@/examples/caliper \
                         @MFEM_SOURCE_DIR@/examples/amgx \
                         @MFEM_SOURCE_DIR@/examples/ginkgo \
+                         @MFEM_SOURCE_DIR@/examples/moonolith \
                         @MFEM_SOURCE_DIR@/examples/hiop \
                         @MFEM_SOURCE_DIR@/examples/petsc \
                         @MFEM_SOURCE_DIR@/examples/pumi \
@@ -103,6 +103,8 @@ namespace mfem {
 * - <a class="el" href="ex31_8cpp_source.html">Example 31</a>: Nedelec H(curl) FEM for the definite anisotropic Maxwell problem
 * - <a class="el" href="ex31p_8cpp_source.html">Example 31p</a>: parallel Nedelec H(curl) FEM for the definite anisotropic Maxwell problem
 * - <a class="el" href="ex32p_8cpp_source.html">Example 32p</a>: parallel anisotropic Maxwell eigensolver
+ * - <a class="el" href="ex33_8cpp_source.html">Example 33</a>: nodal H1 FEM for the fractional Laplacian problem
+ * - <a class="el" href="ex33p_8cpp_source.html">Example 33p</a>: parallel nodal H1 FEM for the fractional Laplacian problem
 *
 * <H4>AmgX Examples</H4>
 * - Variants of Examples
@@ -39,6 +39,7 @@ list(APPEND ALL_EXE_SRCS
  ex29.cpp
  ex30.cpp
  ex31.cpp
+  ex33.cpp
  )

 if (MFEM_USE_MPI)
@@ -75,6 +76,7 @@ if (MFEM_USE_MPI)
    ex30p.cpp
    ex31p.cpp
    ex32p.cpp
+    ex33p.cpp
    )
 endif()

@@ -199,3 +201,7 @@ endif()
 if (MFEM_USE_SUPERLU)
  add_subdirectory(superlu)
 endif()
+
+if(MFEM_USE_MOONOLITH)
+  add_subdirectory(moonolith)
+endif()
@@ -305,6 +305,8 @@ int main(int argc, char *argv[])
         vis_w.precision(8);
         visualize(vis_w, mesh, &x, &w, "Elastic energy density", true);
      }
+      cout << "GLVis visualization paused."
+           << " Press space (in the GLVis window) to resume it.\n";
   }

   double ee0 = oper.ElasticEnergy(x.GetTrueVector());
@@ -351,6 +351,11 @@ int main(int argc, char *argv[])
         vis_w.precision(8);
         visualize(vis_w, pmesh, &x_gf, &w_gf, "Elastic energy density", true);
      }
+      if (myid == 0)
+      {
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
   }

   double ee0 = oper.ElasticEnergy(x_gf);
@@ -1,4 +1,4 @@
-//                    MFEM Example 17 - Parallel Version
+//                       MFEM Example 17 - Parallel Version
 //
 // Compile with: make ex17p
 //
@@ -43,7 +43,7 @@
 #include <sstream>
 #include <iostream>

-// Classes FE_Evolution, RiemannSolver, DomainIntegrator and FaceIntegrator
+// Classes FE_Evolution, RiemannSolver, and FaceIntegrator
 // shared between the serial and parallel version of the example.
 #include "ex18.hpp"

@@ -189,7 +189,7 @@ int main(int argc, char *argv[])
   // 7. Set up the nonlinear form corresponding to the DG discretization of the
   //    flux divergence, and assemble the corresponding mass matrix.
   MixedBilinearForm Aflux(&dfes, &fes);
-   Aflux.AddDomainIntegrator(new DomainIntegrator(dim));
+   Aflux.AddDomainIntegrator(new TransposeIntegrator(new GradientIntegrator()));
   Aflux.Assemble();

   NonlinearForm A(&vfes);
@@ -56,27 +56,6 @@ public:
               const Vector &nor, Vector &flux);
 };

-
-// Constant (in time) mixed bilinear form multiplying the flux grid function.
-// The form is (vec(v), grad(w)) where the trial space = vector L2 space (mesh
-// dim) and test space = scalar L2 space.
-class DomainIntegrator : public BilinearFormIntegrator
-{
-private:
-   Vector shape;
-   DenseMatrix flux;
-   DenseMatrix dshapedr;
-   DenseMatrix dshapedx;
-
-public:
-   DomainIntegrator(const int dim);
-
-   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
-                                       const FiniteElement &test_fe,
-                                       ElementTransformation &Tr,
-                                       DenseMatrix &elmat);
-};
-
 // Interior face term: <F.n(u),[w]>
 class FaceIntegrator : public NonlinearFormIntegrator
 {
@@ -318,60 +297,6 @@ double RiemannSolver::Eval(const Vector &state1, const Vector &state2,
   return maxE;
 }

-// Implementation of class DomainIntegrator
-DomainIntegrator::DomainIntegrator(const int dim) : flux(num_equation, dim) { }
-
-void DomainIntegrator::AssembleElementMatrix2(const FiniteElement &trial_fe,
-                                              const FiniteElement &test_fe,
-                                              ElementTransformation &Tr,
-                                              DenseMatrix &elmat)
-{
-   // Assemble the form (vec(v), grad(w))
-
-   // Trial space = vector L2 space (mesh dim)
-   // Test space  = scalar L2 space
-
-   const int dof_trial = trial_fe.GetDof();
-   const int dof_test = test_fe.GetDof();
-   const int dim = trial_fe.GetDim();
-
-   shape.SetSize(dof_trial);
-   dshapedr.SetSize(dof_test, dim);
-   dshapedx.SetSize(dof_test, dim);
-
-   elmat.SetSize(dof_test, dof_trial * dim);
-   elmat = 0.0;
-
-   const int maxorder = max(trial_fe.GetOrder(), test_fe.GetOrder());
-   const int intorder = 2 * maxorder;
-   const IntegrationRule *ir = &IntRules.Get(trial_fe.GetGeomType(), intorder);
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-
-      // Calculate the shape functions
-      trial_fe.CalcShape(ip, shape);
-      shape *= ip.weight;
-
-      // Compute the physical gradients of the test functions
-      Tr.SetIntPoint(&ip);
-      test_fe.CalcDShape(ip, dshapedr);
-      Mult(dshapedr, Tr.AdjugateJacobian(), dshapedx);
-
-      for (int d = 0; d < dim; d++)
-      {
-         for (int j = 0; j < dof_test; j++)
-         {
-            for (int k = 0; k < dof_trial; k++)
-            {
-               elmat(j, k + d * dof_trial) += shape(k) * dshapedx(j, d);
-            }
-         }
-      }
-   }
-}
-
 // Implementation of class FaceIntegrator
 FaceIntegrator::FaceIntegrator(RiemannSolver &rsolver_, const int dim) :
   rsolver(rsolver_),
@@ -1,4 +1,4 @@
-//                         MFEM Example 18 - Parallel Version
+//                       MFEM Example 18 - Parallel Version
 //
 // Compile with: make ex18
 //
@@ -43,7 +43,7 @@
 #include <sstream>
 #include <iostream>

-// Classes FE_Evolution, RiemannSolver, DomainIntegrator and FaceIntegrator
+// Classes FE_Evolution, RiemannSolver, and FaceIntegrator
 // shared between the serial and parallel version of the example.
 #include "ex18.hpp"

@@ -219,7 +219,7 @@ int main(int argc, char *argv[])
   // 9. Set up the nonlinear form corresponding to the DG discretization of the
   //    flux divergence, and assemble the corresponding mass matrix.
   MixedBilinearForm Aflux(&dfes, &fes);
-   Aflux.AddDomainIntegrator(new DomainIntegrator(dim));
+   Aflux.AddDomainIntegrator(new TransposeIntegrator(new GradientIntegrator()));
   Aflux.Assemble();

   ParNonlinearForm A(&vfes);
@@ -1,4 +1,4 @@
-//                                MFEM Example 21
+//                       MFEM Example 21 - Parallel Version
 //
 // Compile with: make ex21p
 //
@@ -1,4 +1,4 @@
-//                               MFEM Example 22
+//                                MFEM Example 22
 //
 // Compile with: make ex22
 //
@@ -1,4 +1,4 @@
-//                               MFEM Example 24
+//                                MFEM Example 24
 //
 // Compile with: make ex24
 //
@@ -1,4 +1,4 @@
-//                       MFEM Example 27 - Serial Version
+//                                MFEM Example 27
 //
 // Compile with: make ex27
 //
@@ -1,4 +1,4 @@
-//                      MFEM Example 29 - Parallel Version
+//                       MFEM Example 29 - Parallel Version
 //
 // Compile with: make ex29p
 //
@@ -1,4 +1,4 @@
-//                            MFEM Example 30
+//                                MFEM Example 30
 //
 // Compile with: make ex30
 //
@@ -1,4 +1,4 @@
-//                        MFEM Example 30 - Parallel Version
+//                       MFEM Example 30 - Parallel Version
 //
 // Compile with: make ex30p
 //
@@ -45,9 +45,10 @@ int dim;
 int main(int argc, char *argv[])
 {
   // 1. Initialize MPI.
-   MPI_Session mpi;
-   int num_procs = mpi.WorldSize();
-   int myid = mpi.WorldRank();
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();

   // 2. Parse command-line options.
   const char *mesh_file = "../data/inline-quad.mesh";
@@ -119,7 +120,7 @@ int main(int argc, char *argv[])
   }
   ParFiniteElementSpace fespace(&pmesh, fec);
   HYPRE_Int size = fespace.GlobalTrueVSize();
-   if (mpi.Root()) { cout << "Number of H(Curl) unknowns: " << size << endl; }
+   if (Mpi::Root()) { cout << "Number of H(Curl) unknowns: " << size << endl; }

   // 7. Determine the list of true (i.e. parallel conforming) essential
   //    boundary dofs. In this example, the boundary conditions are defined
@@ -181,7 +182,7 @@ int main(int argc, char *argv[])
   // 12. Solve the system AX=B using PCG with the AMS preconditioner from hypre
   if (use_ams)
   {
-      if (mpi.Root())
+      if (Mpi::Root())
      {
         cout << "Size of linear system: "
              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
@@ -199,7 +200,7 @@ int main(int argc, char *argv[])
   else
 #ifdef MFEM_USE_SUPERLU
   {
-      if (mpi.Root())
+      if (Mpi::Root())
      {
         cout << "Size of linear system: "
              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
@@ -212,7 +213,7 @@ int main(int argc, char *argv[])
   }
 #else
   {
-      if (mpi.Root()) { cout << "No solvers available." << endl; }
+      if (Mpi::Root()) { cout << "No solvers available." << endl; }
      return 1;
   }
 #endif
@@ -224,7 +225,7 @@ int main(int argc, char *argv[])
   // 14. Compute and print the H(Curl) norm of the error.
   {
      double error = sol.ComputeHCurlError(&E, &CurlE);
-      if (mpi.Root())
+      if (Mpi::Root())
      {
         cout << "\n|| E_h - E ||_{H(Curl)} = " << error << '\n' << endl;
      }
@@ -41,9 +41,10 @@ double GetScalarMax(const ParGridFunction &x);
 int main(int argc, char *argv[])
 {
   // 1. Initialize MPI.
-   MPI_Session mpi;
-   int num_procs = mpi.WorldSize();
-   int myid = mpi.WorldRank();
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();

   // 2. Parse command-line options.
   const char *mesh_file = "../data/inline-quad.mesh";
@@ -122,7 +123,7 @@ int main(int argc, char *argv[])
   ParFiniteElementSpace fespace_rt(&pmesh, fec_rt);
   HYPRE_Int size_nd = fespace_nd.GlobalTrueVSize();
   HYPRE_Int size_rt = fespace_rt.GlobalTrueVSize();
-   if (mpi.Root())
+   if (Mpi::Root())
   {
      cout << "Number of H(Curl) unknowns: " << size_nd << endl;
      cout << "Number of H(Div) unknowns: " << size_rt << endl;
@@ -164,7 +165,7 @@ int main(int argc, char *argv[])
         // closed surface.
         a.AddDomainIntegrator(new VectorFEMassIntegrator(epsilon));
         shift = 1.0;
-         if (mpi.Root())
+         if (Mpi::Root())
         {
            cout << "Computing eigenvalues shifted by " << shift << endl;
         }
@@ -287,7 +288,7 @@ int main(int argc, char *argv[])

         for (int i=0; i<nev; i++)
         {
-            if (mpi.Root())
+            if (Mpi::Root())
            {
               cout << "Eigenmode " << i+1 << '/' << nev
                    << ", Lambda = " << eigenvalues[i] - shift << endl;
@@ -408,7 +409,7 @@ int main(int argc, char *argv[])
               MPI_Barrier(MPI_COMM_WORLD);
            }
            char c;
-            if (mpi.Root())
+            if (Mpi::Root())
            {
               cout << "press (q)uit or (c)ontinue --> " << flush;
               cin >> c;
@@ -461,7 +462,7 @@ int main(int argc, char *argv[])

         for (int i=0; i<nev; i++)
         {
-            if (mpi.Root())
+            if (Mpi::Root())
            {
               cout << "Eigenmode " << i+1 << '/' << nev
                    << ", Lambda = " << eigenvalues[i] - shift << endl;
@@ -565,7 +566,7 @@ int main(int argc, char *argv[])
               MPI_Barrier(MPI_COMM_WORLD);
            }
            char c;
-            if (mpi.Root())
+            if (Mpi::Root())
            {
               cout << "press (q)uit or (c)ontinue --> " << flush;
               cin >> c;
@@ -591,7 +592,7 @@ int main(int argc, char *argv[])

         for (int i=0; i<nev; i++)
         {
-            if (mpi.Root())
+            if (Mpi::Root())
            {
               cout << "Eigenmode " << i+1 << '/' << nev
                    << ", Lambda = " << eigenvalues[i] - shift << endl;
@@ -620,7 +621,7 @@ int main(int argc, char *argv[])
            MPI_Barrier(MPI_COMM_WORLD);

            char c;
-            if (mpi.Root())
+            if (Mpi::Root())
            {
               cout << "press (q)uit or (c)ontinue --> " << flush;
               cin >> c;
@@ -0,0 +1,192 @@
+//                                MFEM Example 33
+//
+// Compile with: make ex33
+//
+// Sample runs:  ex33 -m ../data/square-disc.mesh -alpha 0.33 -o 2
+//               ex33 -m ../data/star.mesh -alpha 0.99 -o 3
+//               ex33 -m ../data/inline-quad.mesh -alpha 0.5 -o 3
+//               ex33 -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3
+//               ex33 -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
+//
+// Description:
+//
+//  In this example we solve the following fractional PDE with MFEM:
+//
+//    ( - Δ )^α u = f  in Ω,      u = 0  on ∂Ω,      0 < α < 1,
+//
+//  To solve this FPDE, we rely on a rational approximation [2] of the normal
+//  linear operator A^{-α}, where A = - Δ (with associated homogeneous
+//  boundary conditions). Namely, we first approximate the operator
+//
+//    A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1},      d_0 = 0,   d_i > 0,
+//
+//  where I is the L2-identity operator and the coefficients c_i and d_i
+//  are generated offline to a prescribed accuracy in a pre-processing step.
+//  We use the triple-A algorithm [1] to generate the rational approximation
+//  that this partial fractional expansion derives from. We then solve N+1
+//  independent integer-order PDEs,
+//
+//    A u_i + d_i u_i = c_i f  in Ω,      u_i = 0  on ∂Ω,      i=0,...,N,
+//
+//  using MFEM and sum u_i to arrive at an approximate solution of the FPDE
+//
+//    u ≈ Σ_{i=0}^N u_i.
+//
+// References:
+//
+// [1] Nakatsukasa, Y., Sète, O., & Trefethen, L. N. (2018). The AAA algorithm
+//     for rational approximation. SIAM Journal on Scientific Computing, 40(3),
+//     A1494-A1522.
+//
+// [2] Harizanov, S., Lazarov, R., Margenov, S., Marinov, P., & Pasciak, J.
+//     (2020). Analysis of numerical methods for spectral fractional elliptic
+//     equations based on the best uniform rational approximation. Journal of
+//     Computational Physics, 408, 109285.
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "ex33.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int num_refs = 3;
+   bool visualization = true;
+   double alpha = 0.5;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&num_refs, "-r", "--refs",
+                  "Number of uniform refinements");
+   args.AddOption(&alpha, "-alpha", "--alpha",
+                  "Fractional exponent");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Array<double> coeffs, poles;
+
+   // 2. Compute the coefficients that define the integer-order PDEs.
+   ComputePartialFractionApproximation(alpha,coeffs,poles);
+
+   // 3. Read the mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. Refine the mesh to increase the resolution.
+   for (int i = 0; i < num_refs; i++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 5. Define a finite element space on the mesh.
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   FiniteElementSpace fespace(&mesh, fec);
+   cout << "Number of finite element unknowns: "
+        << fespace.GetTrueVSize() << endl;
+
+   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
+   Array<int> ess_tdof_list;
+   if (mesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 7. Define diffusion coefficient, load, and solution GridFunction.
+   ConstantCoefficient f(1.0);
+   ConstantCoefficient one(1.0);
+   GridFunction u(&fespace);
+   u = 0.;
+
+   // 8. Prepare for visualization.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream xout, uout;
+   ostringstream oss_x, oss_u;
+   if (visualization)
+   {
+      xout.open(vishost, visport);
+      xout.precision(8);
+      uout.open(vishost, visport);
+      uout.precision(8);
+   }
+
+   for (int i = 0; i < coeffs.Size(); i++)
+   {
+      // 9. Set up the linear form b(.) for integer-order PDE solve.
+      LinearForm b(&fespace);
+      ProductCoefficient cf(coeffs[i], f);
+      b.AddDomainIntegrator(new DomainLFIntegrator(cf));
+      b.Assemble();
+
+      // 10. Define GridFunction for integer-order PDE solve.
+      GridFunction x(&fespace);
+      x = 0.0;
+
+      // 11. Set up the bilinear form a(.,.) for integer-order PDE solve.
+      BilinearForm a(&fespace);
+      a.AddDomainIntegrator(new DiffusionIntegrator(one));
+      ConstantCoefficient c2(-poles[i]);
+      a.AddDomainIntegrator(new MassIntegrator(c2));
+      a.Assemble();
+
+      // 12. Assemble the bilinear form and the corresponding linear system.
+      OperatorPtr A;
+      Vector B, X;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+      // 13. Solve the linear system A X = B.
+      GSSmoother M((SparseMatrix&)(*A));
+
+      mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
+                << " u = " << coeffs[i] << " f " << endl;
+      PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
+
+      // 14. Recover the solution as a finite element grid function.
+      a.RecoverFEMSolution(X, b, x);
+
+      // 15. Accumulate integer-order PDE solutions.
+      u+=x;
+
+      // 16. Send the solutions by socket to a GLVis server.
+      if (visualization)
+      {
+         oss_x.str(""); oss_x.clear();
+         oss_x << "Solution of PDE -Δ u + " << -poles[i]
+               << " u = " << coeffs[i] << " f";
+         xout << "solution\n" << mesh << x
+              << "window_title '" << oss_x.str() << "'" << flush;
+
+         oss_u.str(""); oss_u.clear();
+         oss_u << "Solution of fractional PDE -Δ^" << alpha
+               << " u = f";
+         uout << "solution\n" << mesh << u
+              << "window_title '" << oss_u.str() << "'" << flush;
+      }
+   }
+
+   // 17. Free the used memory.
+   delete fec;
+   return 0;
+}
@@ -0,0 +1,390 @@
+//                  MFEM Example 33 - Serial/Parallel Shared Code
+//                      (Implementation of the AAA algorithm)
+//
+//  Here, we implement the triple-A algorithm [1] for the rational approximation
+//  of complex-valued functions,
+//
+//          p(z)/q(z) ≈ f(z).
+//
+//  In this file, we always assume f(z) = z^{-α}. The triple-A algorithm
+//  provides a robust, accurate approximation in rational barycentric form.
+//  This representation must be transformed into a partial fraction
+//  representation in order to be used to solve a spectral FPDE.
+//
+//  More specifically, we first expand the numerator in terms of the zeros of
+//  the rational approximation,
+//
+//          p(z) ∝ Π_i (z - z_i),
+//
+//  and expand the denominator in terms of the poles of the rational
+//  approximation,
+//
+//          q(z) ∝ Π_i (z - p_i).
+//
+//  We then use these zeros and poles to derive the partial fraction expansion
+//
+//          f(z) ≈ p(z)/q(z) = Σ_i c_i / (z - p_i).
+//
+//  [1] Nakatsukasa, Y., Sète, O., & Trefethen, L. N. (2018). The AAA algorithm
+//      for rational approximation. SIAM Journal on Scientific Computing, 40(3),
+//      A1494-A1522.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+/** RationalApproximation_AAA: compute the rational approximation (RA) of data
+    @a val [in] at the set of points @a pt [in].
+
+    @param[in]  val        Vector of data values
+    @param[in]  pt         Vector of sample points
+    @param[in]  tol        Relative tolerance
+    @param[in]  max_order  Maximum number of terms (order) of the RA
+    @param[out] z          Support points of the RA in rational barycentric form
+    @param[out] f          Data values at support points @a z
+    @param[out] w          Weights of the RA in rational barycentric form
+
+    See pg. A1501 of Nakatsukasa et al. [1]. */
+void RationalApproximation_AAA(const Vector &val, const Vector &pt,
+                               Array<double> &z, Array<double> &f, Vector &w,
+                               double tol, int max_order)
+{
+
+   // number of sample points
+   int size = val.Size();
+   MFEM_VERIFY(pt.Size() == size, "size mismatch");
+
+   // Initializations
+   Array<int> J(size);
+   for (int i = 0; i < size; i++) { J[i] = i; }
+   z.SetSize(0);
+   f.SetSize(0);
+
+   DenseMatrix C, Ctemp, A, Am;
+   // auxiliary arrays and vectors
+   Vector f_vec;
+   Array<double> c_i;
+
+   // mean of the value vector
+   Vector R(val.Size());
+   double mean_val = val.Sum()/size;
+
+   for (int i = 0; i<R.Size(); i++) { R(i) = mean_val; }
+
+   for (int k = 0; k < max_order; k++)
+   {
+      // select next support point
+      int idx = 0;
+      double tmp_max = 0;
+      for (int j = 0; j < size; j++)
+      {
+         double tmp = abs(val(j)-R(j));
+         if (tmp > tmp_max)
+         {
+            tmp_max = tmp;
+            idx = j;
+         }
+      }
+
+      // Append support points and data values
+      z.Append(pt(idx));
+      f.Append(val(idx));
+
+      // Update index vector
+      J.DeleteFirst(idx);
+
+      // next column in Cauchy matrix
+      Array<double> C_tmp(size);
+      for (int j = 0; j < size; j++)
+      {
+         C_tmp[j] = 1.0/(pt(j)-pt(idx));
+      }
+      c_i.Append(C_tmp);
+      int h_C = C_tmp.Size();
+      int w_C = k+1;
+      C.UseExternalData(c_i.GetData(),h_C,w_C);
+
+      Ctemp = C;
+
+      f_vec.SetDataAndSize(f.GetData(),f.Size());
+      Ctemp.InvLeftScaling(val);
+      Ctemp.RightScaling(f_vec);
+
+      A.SetSize(C.Height(), C.Width());
+      Add(C,Ctemp,-1.0,A);
+      A.LeftScaling(val);
+
+      int h_Am = J.Size();
+      int w_Am = A.Width();
+      Am.SetSize(h_Am,w_Am);
+      for (int i = 0; i<h_Am; i++)
+      {
+         int ii = J[i];
+         for (int j = 0; j<w_Am; j++)
+         {
+            Am(i,j) = A(ii,j);
+         }
+      }
+
+#ifdef MFEM_USE_LAPACK
+      DenseMatrixSVD svd(Am,false,true);
+      svd.Eval(Am);
+      DenseMatrix &v = svd.RightSingularvectors();
+      v.GetRow(k,w);
+#else
+      mfem_error("Compiled without LAPACK");
+#endif
+
+      // N = C*(w.*f); D = C*w; % numerator and denominator
+      Vector aux(w);
+      aux *= f_vec;
+      Vector N(C.Height()); // Numerator
+      C.Mult(aux,N);
+      Vector D(C.Height()); // Denominator
+      C.Mult(w,D);
+
+      R = val;
+      for (int i = 0; i<J.Size(); i++)
+      {
+         int ii = J[i];
+         R(ii) = N(ii)/D(ii);
+      }
+
+      Vector verr(val);
+      verr-=R;
+
+      if (verr.Normlinf() <= tol*val.Normlinf()) { break; }
+   }
+}
+
+/** ComputePolesAndZeros: compute the @a poles [out] and @a zeros [out] of the
+    rational function f(z) = C p(z)/q(z) from its ration barycentric form.
+
+    @param[in]  z      Support points in rational barycentric form
+    @param[in]  f      Data values at support points @a z
+    @param[in]  w      Weights in rational barycentric form
+    @param[out] poles  Array of poles (roots of p(z))
+    @param[out] zeros  Array of zeros (roots of q(z))
+    @param[out] scale  Scaling constant in f(z) = C p(z)/q(z)
+
+    See pg. A1501 of Nakatsukasa et al. [1]. */
+void ComputePolesAndZeros(const Vector &z, const Vector &f, const Vector &w,
+                          Array<double> & poles, Array<double> & zeros, double &scale)
+{
+   // Initialization
+   poles.SetSize(0);
+   zeros.SetSize(0);
+
+   // Compute the poles
+   int m = w.Size();
+   DenseMatrix B(m+1); B = 0.;
+   DenseMatrix E(m+1); E = 0.;
+   for (int i = 1; i<=m; i++)
+   {
+      B(i,i) = 1.;
+      E(0,i) = w(i-1);
+      E(i,0) = 1.;
+      E(i,i) = z(i-1);
+   }
+
+#ifdef MFEM_USE_LAPACK
+   DenseMatrixGeneralizedEigensystem eig1(E,B);
+   eig1.Eval();
+   Vector & evalues = eig1.EigenvaluesRealPart();
+   for (int i = 0; i<evalues.Size(); i++)
+   {
+      if (IsFinite(evalues(i)))
+      {
+         poles.Append(evalues(i));
+      }
+   }
+#else
+   mfem_error("Compiled without LAPACK");
+#endif
+   // compute the zeros
+   B = 0.;
+   E = 0.;
+   for (int i = 1; i<=m; i++)
+   {
+      B(i,i) = 1.;
+      E(0,i) = w(i-1) * f(i-1);
+      E(i,0) = 1.;
+      E(i,i) = z(i-1);
+   }
+
+#ifdef MFEM_USE_LAPACK
+   DenseMatrixGeneralizedEigensystem eig2(E,B);
+   eig2.Eval();
+   evalues = eig2.EigenvaluesRealPart();
+   for (int i = 0; i<evalues.Size(); i++)
+   {
+      if (IsFinite(evalues(i)))
+      {
+         zeros.Append(evalues(i));
+      }
+   }
+#else
+   mfem_error("Compiled without LAPACK");
+#endif
+
+   scale = w * f / w.Sum();
+}
+
+/** PartialFractionExpansion: compute the partial fraction expansion of the
+    rational function f(z) = Σ_i c_i / (z - p_i) from its @a poles [in] and
+    @a zeros [in].
+
+    @param[in]  poles   Array of poles (same as p_i above)
+    @param[in]  zeros   Array of zeros
+    @param[in]  scale   Scaling constant
+    @param[out] coeffs  Coefficients c_i */
+void PartialFractionExpansion(double scale, Array<double> & poles,
+                              Array<double> & zeros, Array<double> & coeffs)
+{
+   int psize = poles.Size();
+   int zsize = zeros.Size();
+   coeffs.SetSize(psize);
+   coeffs = scale;
+
+   for (int i=0; i<psize; i++)
+   {
+      double tmp_numer=1.0;
+      for (int j=0; j<zsize; j++)
+      {
+         tmp_numer *= poles[i]-zeros[j];
+      }
+
+      double tmp_denom=1.0;
+      for (int k=0; k<psize; k++)
+      {
+         if (k != i) { tmp_denom *= poles[i]-poles[k]; }
+      }
+      coeffs[i] *= tmp_numer / tmp_denom;
+   }
+}
+
+
+/** ComputePartialFractionApproximation: compute a rational approximation (RA)
+    in partial fraction form, e.g., f(z) ≈ Σ_i c_i / (z - p_i), from sampled
+    values of the function f(z) = z^{-a}, 0 < a < 1.
+
+    @param[in]  alpha         Exponent a in f(z) = z^-a
+    @param[in] lmax, npoints  f(z) is uniformly sampled @a npoints times in the
+                              interval [ 0, @a lmax ]
+    @param[in]  tol           Relative tolerance
+    @param[in]  max_order     Maximum number of terms (order) of the RA
+    @param[out] coeffs        Coefficients c_i
+    @param[out] poles         Poles p_i
+
+    NOTES: When MFEM is not built with LAPACK support, only @a alpha = 0.33,
+           0.5, and 0.99 are possible. In this case, if @a alpha != 0.33 and
+           @a alpha != 0.99, then @a alpha = 0.5 is used by default.
+
+   See pg. A1501 of Nakatsukasa et al. [1]. */
+void ComputePartialFractionApproximation(double & alpha,
+                                         Array<double> & coeffs, Array<double> & poles,
+                                         double lmax = 1000.,
+                                         double tol=1e-10, int npoints = 1000,
+                                         int max_order = 100)
+{
+   MFEM_VERIFY(alpha < 1., "alpha must be less than 1");
+   MFEM_VERIFY(alpha > 0., "alpha must be greater than 0");
+   MFEM_VERIFY(npoints > 2, "npoints must be greater than 2");
+   MFEM_VERIFY(lmax > 0,  "lmin must be greater than 0");
+   MFEM_VERIFY(tol > 0,  "tol must be greater than 0");
+
+   bool print_warning = true;
+#ifdef MFEM_USE_MPI
+   if ((Mpi::IsInitialized() && !Mpi::Root())) { print_warning = false; }
+#endif
+
+#ifndef MFEM_USE_LAPACK
+   if (print_warning)
+   {
+      mfem::out
+            << "\nMFEM is compiled without LAPACK.\nUsing precomputed values for PartialFractionApproximation. \n"
+            << "Only alpha = 0.33, 0.5, and 0.99 are available.\nThe default is alpha = 0.5."
+            << std::endl;
+   }
+   const double eps = std::numeric_limits<double>::epsilon();
+
+   if (abs(alpha - 0.33) < eps)
+   {
+      coeffs = Array<double> ({1.821898e+03, 9.101221e+01, 2.650611e+01,
+                               1.174937e+01, 6.140444e+00, 3.441713e+00,
+                               1.985735e+00, 1.162634e+00, 6.891560e-01,
+                               4.111574e-01, 2.298736e-01});
+      poles = Array<double> ({-4.155583e+04, -2.956285e+03, -8.331715e+02,
+                              -3.139332e+02, -1.303448e+02, -5.563385e+01,
+                              -2.356255e+01, -9.595516e+00, -3.552160e+00,
+                              -1.032136e+00, -1.241480e-01});
+   }
+   else if (abs(alpha - 0.99) < eps)
+   {
+      coeffs = Array<double>({2.919591e-02, 1.419750e-02, 1.065798e-02,
+                              9.395094e-03, 8.915329e-03, 8.822991e-03,
+                              9.058247e-03, 9.814521e-03, 1.180396e-02,
+                              1.834554e-02, 9.840482e-01});
+      poles = Array<double> ({-1.069683e+04, -1.769370e+03, -5.718374e+02,
+                              -2.242095e+02, -9.419132e+01, -4.031012e+01,
+                              -1.701525e+01, -6.810088e+00, -2.382810e+00,
+                              -5.700059e-01, -1.384324e-03});
+   }
+   else
+   {
+      if (abs(alpha - 0.5) > eps && print_warning)
+      {
+         alpha = 0.5;
+      }
+      coeffs = Array<double>({2.290262e+02, 2.641819e+01, 1.005566e+01,
+                              5.390411e+00, 3.340725e+00, 2.211205e+00,
+                              1.508883e+00, 1.049474e+00, 7.462709e-01,
+                              5.482686e-01, 4.232510e-01, 3.578967e-01});
+      poles = Array<double>({-3.168211e+04, -3.236077e+03, -9.868287e+02,
+                             -3.945597e+02, -1.738889e+02, -7.925178e+01,
+                             -3.624992e+01, -1.629196e+01, -6.982956e+00,
+                             -2.679984e+00, -7.782607e-01, -7.649166e-02});
+   }
+
+   if (print_warning)
+   {
+      mfem::out << "Using precomputed values for alpha = "
+                << alpha << "\n" << std::endl;
+   }
+
+
+   return;
+#endif
+
+   Vector x(npoints);
+   Vector val(npoints);
+   double dx = lmax / (double)(npoints-1);
+   for (int i = 0; i<npoints; i++)
+   {
+      x(i) = dx * (double)i;
+      val(i) = pow(x(i),1.-alpha);
+   }
+
+   // Apply triple-A algorithm to f(x) = x^{1-a}
+   Array<double> z, f;
+   Vector w;
+   RationalApproximation_AAA(val,x,z,f,w,tol,max_order);
+
+   Vector vecz, vecf;
+   vecz.SetDataAndSize(z.GetData(), z.Size());
+   vecf.SetDataAndSize(f.GetData(), f.Size());
+
+   // Compute poles and zeros for RA of f(x) = x^{1-a}
+   double scale;
+   Array<double> zeros;
+   ComputePolesAndZeros(vecz, vecf, w, poles, zeros, scale);
+
+   // Remove the zero at x=0, thus, delivering a RA for f(x) = x^{-a}
+   zeros.DeleteFirst(0.0);
+
+   // Compute partial fraction approximation of f(x) = x^{-a}
+   PartialFractionExpansion(scale, poles, zeros, coeffs);
+}
@@ -0,0 +1,290 @@
+//                       MFEM Example 33 - Parallel Version
+//
+// Compile with: make ex33p
+//
+// Sample runs:  mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 0.33 -o 2
+//               mpirun -np 4 ex33p -m ../data/star.mesh -alpha 0.99 -o 3
+//               mpirun -np 4 ex33p -m ../data/inline-quad.mesh -alpha 0.5 -o 3
+//               mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3
+//               mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
+//
+// Description:
+//
+//  In this example we solve the following fractional PDE with MFEM:
+//
+//    ( - Δ )^α u = f  in Ω,      u = 0  on ∂Ω,      0 < α < 1,
+//
+//  To solve this FPDE, we rely on a rational approximation [2] of the normal
+//  linear operator A^{-α}, where A = - Δ (with associated homogeneous
+//  boundary conditions). Namely, we first approximate the operator
+//
+//    A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1},      d_0 = 0,   d_i > 0,
+//
+//  where I is the L2-identity operator and the coefficients c_i and d_i
+//  are generated offline to a prescribed accuracy in a pre-processing step.
+//  We use the triple-A algorithm [1] to generate the rational approximation
+//  that this partial fractional expansion derives from. We then solve N+1
+//  independent integer-order PDEs,
+//
+//    A u_i + d_i u_i = c_i f  in Ω,      u_i = 0  on ∂Ω,      i=0,...,N,
+//
+//  using MFEM and sum u_i to arrive at an approximate solution of the FPDE
+//
+//    u ≈ Σ_{i=0}^N u_i.
+//
+// References:
+//
+// [1] Nakatsukasa, Y., Sète, O., & Trefethen, L. N. (2018). The AAA algorithm
+//     for rational approximation. SIAM Journal on Scientific Computing, 40(3),
+//     A1494-A1522.
+//
+// [2] Harizanov, S., Lazarov, R., Margenov, S., Marinov, P., & Pasciak, J.
+//     (2020). Analysis of numerical methods for spectral fractional elliptic
+//     equations based on the best uniform rational approximation. Journal of
+//     Computational Physics, 408, 109285.
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "ex33.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int num_refs = 3;
+   bool visualization = true;
+   bool visualize_x = false;
+   double alpha = 0.5;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&num_refs, "-r", "--refs",
+                  "Number of uniform refinements");
+   args.AddOption(&alpha, "-alpha", "--alpha",
+                  "Fractional exponent");
+   args.AddOption(&visualize_x, "-vis_x", "--visualize_x", "-no-vis_x",
+                  "--no-visualization_x",
+                  "Enable or disable GLVis visualization of each integer-order PDE solution.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization of the fractional PDE solution.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (Mpi::Root())
+   {
+      args.PrintOptions(cout);
+   }
+
+   Array<double> coeffs, poles;
+
+   // 2. Compute the coefficients that define the integer-order PDEs.
+   ComputePartialFractionApproximation(alpha,coeffs,poles);
+
+   int num_par_solves;
+   int max_par_solves = max(1,num_procs/2);
+   for (num_par_solves=max_par_solves; num_par_solves>0; num_par_solves--)
+   {
+      if (num_procs%num_par_solves==0 && num_par_solves<coeffs.Size())
+      {
+         break;
+      }
+   }
+   if (num_par_solves == 1) {num_par_solves = num_procs;}
+
+   int solver_ranks = num_procs/num_par_solves;
+
+   // 3. Split the MPI communicator:
+   //    row_comm is used for parallel partition of the mesh
+   //    col_comm is used for independent integer-order solves
+   int row_color = myid / solver_ranks; // Determine color based on row
+   int col_color = myid % solver_ranks; // Determine color based on col
+
+   MPI_Comm row_comm, col_comm;
+   MPI_Comm_split(MPI_COMM_WORLD, row_color, myid, &row_comm);
+   MPI_Comm_split(MPI_COMM_WORLD, col_color, myid, &col_comm);
+
+   int row_rank, row_size, col_rank, col_size;
+   MPI_Comm_rank(row_comm, &row_rank);
+   MPI_Comm_size(row_comm, &row_size);
+   MPI_Comm_rank(col_comm, &col_rank);
+   MPI_Comm_size(col_comm, &col_size);
+
+   if (Mpi::Root())
+   {
+      mfem::out << "\nTotal number of MPI ranks = " << num_procs << endl;
+      mfem::out << "Number of independent parallel solves = " << col_size << endl;
+      mfem::out << "Number of MPI ranks within each solve = " << row_size
+                <<"\n" << endl;
+   }
+
+   // 4. Read the mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the mesh to increase the resolution.
+   for (int i = 0; i < num_refs; i++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   ParMesh pmesh(row_comm, mesh);
+   mesh.Clear();
+
+   // 6. Define a finite element space on the mesh.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+   if (Mpi::Root())
+   {
+      cout << "Number of finite element unknowns: "
+           << fespace.GetTrueVSize() << endl;
+   }
+
+   // 7. Determine the list of true (i.e. conforming) essential boundary dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 8. Define diffusion coefficient, load, and solution GridFunction.
+   ConstantCoefficient f(1.0);
+   ConstantCoefficient one(1.0);
+   ParGridFunction u(&fespace);
+   ParGridFunction x(&fespace);
+   u = 0.0;
+
+   // 9. Set up the linear form b(.) for integer-order PDE solves.
+   ParLinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(f));
+   b.Assemble();
+
+   int my_coeff_size = max(coeffs.Size()/col_size,1);
+   int ibeg = col_rank*my_coeff_size;
+   if (ibeg + 2*my_coeff_size > coeffs.Size())
+   {
+      my_coeff_size = coeffs.Size()-col_rank*my_coeff_size;
+   }
+   else if (ibeg > coeffs.Size() - 1)
+   {
+      my_coeff_size = 0;
+   }
+
+   int iend = ibeg+my_coeff_size;
+
+
+   for (int i = ibeg; i < iend; i++)
+   {
+      // 10. Reset GridFunction for integer-order PDE solve.
+      x = 0.0;
+
+      // 11. Set up the bilinear form a(.,.) for integer-order PDE solve.
+      ParBilinearForm a(&fespace);
+      a.AddDomainIntegrator(new DiffusionIntegrator(one));
+      ConstantCoefficient d_i(-poles[i]);
+      a.AddDomainIntegrator(new MassIntegrator(d_i));
+      a.Assemble();
+
+      // 12. Assemble the bilinear form and the corresponding linear system.
+      OperatorPtr A;
+      Vector B, X;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+      // 13. Solve the linear system A X = B.
+      HypreBoomerAMG * prec = new HypreBoomerAMG;
+      prec->SetPrintLevel(-1);
+
+      int print_level = (col_rank==0) ? 3 : 0;
+      if (Mpi::Root())
+      {
+         mfem::out << "\nMPI rank " <<  myid
+                   << ": Solving PDE -Δ u + " << -poles[i]
+                   << " u = " << coeffs[i] << " f " << endl;
+      }
+      CGSolver cg(row_comm);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(2000);
+      cg.SetPrintLevel(print_level);
+      cg.SetPreconditioner(*prec);
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+      delete prec;
+
+      // 14. Recover the solution as a finite element grid function.
+      a.RecoverFEMSolution(X, b, x);
+
+      // 15. Accumulate integer-order PDE solutions.
+      x *= coeffs[i];
+      u += x;
+
+      // 16. Send integer-order PDE solutions to a GLVis server.
+      if (visualize_x)
+      {
+         if (col_rank > 0 && i < iend-1)
+         {
+            MPI_Status status;
+            MPI_Recv(nullptr,0,MPI_INT, col_rank-1,0,col_comm,&status);
+         }
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream xout(vishost, visport);
+         xout.precision(8);
+         ostringstream oss;
+         oss << "Solution of PDE -Δ u + " << -poles[i]
+             << " u = " << coeffs[i] << " f" ;
+         xout << "parallel " << row_size << " " << row_rank << "\n";
+         xout << "solution\n" << pmesh << x
+              << "window_title '" << oss.str() << "'" << flush;
+         if (col_rank < col_size-1)
+         {
+            MPI_Send(nullptr,0,MPI_INT,col_rank+1,0,col_comm);
+         }
+      }
+   }
+
+   // 17. Accumulate for the fractional PDE solution
+   MPI_Allreduce(MPI_IN_PLACE, u.GetData(), u.Size(),
+                 MPI_DOUBLE, MPI_SUM,col_comm);
+
+   // 18. Send fractional PDE solution to a GLVis server.
+   if (visualization)
+   {
+      if (col_rank == 0)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream uout(vishost, visport);
+         uout.precision(8);
+         ostringstream oss;
+         oss << "Solution of fractional PDE -Δ^" << alpha
+             << " u = f" ;
+         uout << "parallel " << row_size << " " << row_rank << "\n";
+         uout << "solution\n" << pmesh << u
+              << "window_title '" << oss.str() << "'" << flush;
+      }
+   }
+
+   return 0;
+}
@@ -55,6 +55,7 @@ int main(int argc, char *argv[])
   bool pa = false;
   const char *device_config = "cpu";
   int max_dofs = 50000;
+   bool LSZZ = false;
   bool visualization = true;

   OptionsParser args(argc, argv);
@@ -68,6 +69,9 @@ int main(int argc, char *argv[])
                  "Device configuration string, see Device::Configure().");
   args.AddOption(&max_dofs, "-md", "--max-dofs",
                  "Stop after reaching this many degrees of freedom.");
+   args.AddOption(&LSZZ, "-ls", "--ls-zz", "-no-ls",
+                  "--no-ls-zz",
+                  "Switch to least-squares ZZ estimator.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
@@ -151,15 +155,29 @@ int main(int argc, char *argv[])
   //     recover a smoothed flux (gradient) that is subtracted from the element
   //     flux to get an error indicator. We need to supply the space for the
   //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
-   FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
-   ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
-   estimator.SetAnisotropic();
+   ErrorEstimator *estimator{nullptr};
+
+   if (LSZZ)
+   {
+      estimator = new LSZienkiewiczZhuEstimator(*integ, x);
+      if (dim == 3 && mesh.GetElementType(0) != Element::HEXAHEDRON)
+      {
+         dynamic_cast<LSZienkiewiczZhuEstimator *>
+         (estimator)->SetTichonovRegularization();
+      }
+   }
+   else
+   {
+      auto flux_fes = new FiniteElementSpace(&mesh, &fec, sdim);
+      estimator = new ZienkiewiczZhuEstimator(*integ, x, flux_fes);
+      dynamic_cast<ZienkiewiczZhuEstimator *>(estimator)->SetAnisotropic();
+   }

   // 11. A refiner selects and refines elements based on a refinement strategy.
   //     The strategy here is to refine elements with errors larger than a
   //     fraction of the maximum element error. Other strategies are possible.
   //     The refiner will call the given error estimator.
-   ThresholdRefiner refiner(estimator);
+   ThresholdRefiner refiner(*estimator);
   refiner.SetTotalErrorFraction(0.7);

   // 12. The main AMR loop. In each iteration we solve the problem on the
@@ -256,5 +274,6 @@ int main(int argc, char *argv[])
      b.Update();
   }

+   delete estimator;
   return 0;
 }
@@ -23,10 +23,10 @@ MFEM_LIB_FILE = mfem_is_not_built

 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
  ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-  ex31
+  ex31 ex33
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
  ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
-  ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p
+  ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p \
  ex24p ex25p ex26p
@@ -58,6 +58,9 @@ endif
 ifeq ($(MFEM_USE_SUPERLU),YES)
   SUBDIRS += superlu
 endif
+ifeq ($(MFEM_USE_MOONOLITH),YES)
+   SUBDIRS += moonolith
+endif
 ifeq ($(MFEM_USE_CALIPER),YES)
   SUBDIRS += caliper
 endif
@@ -88,8 +91,11 @@ $(SUBDIRS_TPRINT):

 # Additional dependencies
 ex18: $(SRC)ex18.hpp
+ex33: $(SRC)ex33.hpp
+
 ifeq ($(MFEM_USE_MPI),YES)
 ex18p: $(SRC)ex18.hpp
+ex33p: $(SRC)ex33.hpp
 endif

 MFEM_TESTS = EXAMPLES
@@ -0,0 +1,57 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+set(MOONOLITH_EXAMPLES_SRCS)
+list(APPEND MOONOLITH_EXAMPLES_SRCS ex1.cpp)
+
+if (MFEM_USE_MPI)
+  list(APPEND MOONOLITH_EXAMPLES_SRCS ex1p.cpp ex2p.cpp)
+endif()
+
+# Include the source directory where mfem.hpp and mfem-performance.hpp are.
+include_directories(BEFORE ${PROJECT_BINARY_DIR})
+
+# Add "test_moonolith" target, see below.
+add_custom_target(test_moonolith
+  ${CMAKE_CTEST_COMMAND} -R moonolith USES_TERMINAL)
+
+# Add one executable per cpp file, adding "moonolith_" as prefix. Sets
+# "test_moonolith" as a target that depends on the given examples.
+set(PFX moonolith_)
+add_mfem_examples(MOONOLITH_EXAMPLES_SRCS ${PFX} "" test_moonolith)
+
+# Testing.
+# The MOONOLITH tests can be run separately using the target "test_moonolith"
+# which builds the examples and runs:
+#   ctest -R moonolith
+
+# Add the tests: one test per source file.
+foreach(SRC_FILE ${MOONOLITH_EXAMPLES_SRCS})
+  get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
+  string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
+  string(TOUPPER ${TEST_NAME} UP_TEST_NAME)
+  set(TEST_NAME ${PFX}${TEST_NAME})
+
+  set(THIS_TEST_OPTIONS "-no-vis")
+  list(APPEND THIS_TEST_OPTIONS ${${UP_TEST_NAME}_TEST_OPTS})
+  # message(STATUS "Test ${TEST_NAME} options: ${THIS_TEST_OPTIONS}")
+
+  if (NOT (${TEST_NAME} MATCHES ".*p$"))
+    add_test(NAME ${TEST_NAME}_ser
+      COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
+  else()
+    add_test(NAME ${TEST_NAME}_np=4
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS}
+      ${MPIEXEC_POSTFLAGS})
+  endif()
+endforeach()
@@ -0,0 +1,18 @@
+                    Finite Element Discretization Library
+                                   __
+                       _ __ ___   / _|  ___  _ __ ___
+                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
+                      | | | | | ||  _||  __/| | | | | |
+                      |_| |_| |_||_|   \___||_| |_| |_|
+
+                               https://mfem.org
+
+This directory contains modifications of the example codes that illustrate the
+use of MFEM features based on ParMoonolith library for variational information
+transfer between non-matching meshes.
+
+To build these examples, make sure that MFEM is configured with the option
+"MFEM_USE_MOONOLITH = YES".
+
+For the parallel version add also the options MFEM_USE_MPI = YES", and
+MFEM_USE_METIS = YES". See the top-level INSTALL file for details.
@@ -0,0 +1,229 @@
+//                         MFEM + Moonolith Example 1
+//
+// Compile with: make ex1
+//
+// Moonolith sample runs:
+//   ex1
+//   ex1 --source_refinements 1 --dest_refinements 2
+//   ex1 --source_refinements 1 --dest_refinements 2 --use_vector_fe
+//   ex1 -s ../../data/inline-hex.mesh -d ../../data/inline-tet.mesh
+//
+// Description:  This example code demonstrates the use of MFEM for transferring
+//               discrete fields from one finite element mesh to another. The
+//               meshes can be of arbitrary shape and completely unrelated with
+//               each other. This feature can be used for implementing immersed
+//               domain methods for fluid-structure interaction or general
+//               multi-physics applications.
+//
+//               This particular example is only for serial runtimes.
+
+#include "example_utils.hpp"
+#include "mfem.hpp"
+
+using namespace mfem;
+using namespace std;
+
+int main(int argc, char *argv[])
+{
+   // Init transfer library context
+   InitTransfer(argc, argv);
+
+   const char *source_mesh_file = "../../data/inline-tri.mesh";
+   const char *destination_mesh_file = "../../data/inline-quad.mesh";
+
+   int src_n_refinements = 0;
+   int dest_n_refinements = 0;
+   int source_fe_order = 1;
+   int dest_fe_order = 1;
+   bool visualization = true;
+   bool use_vector_fe = false;
+   bool verbose = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&source_mesh_file, "-s", "--source_mesh",
+                  "Mesh file to use for src.");
+   args.AddOption(&destination_mesh_file, "-d", "--destination_mesh",
+                  "Mesh file to use for dest.");
+   args.AddOption(&src_n_refinements, "-sr", "--source_refinements",
+                  "Number of src refinements");
+   args.AddOption(&dest_n_refinements, "-dr", "--dest_refinements",
+                  "Number of dest refinements");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&source_fe_order, "-so", "--source_fe_order",
+                  "Order of the src finite elements");
+   args.AddOption(&dest_fe_order, "-do", "--dest_fe_order",
+                  "Order of the dest finite elements");
+   args.AddOption(&verbose, "-verb", "--verbose", "--no-verb", "--no-verbose",
+                  "Enable/Disable verbose output");
+   args.AddOption(&use_vector_fe, "-vfe", "--use_vector_fe", "-no-vfe",
+                  "--no-vector_fe", "Use vector finite elements");
+   args.Parse();
+   check_options(args);
+
+   shared_ptr<Mesh> src_mesh, dest_mesh;
+
+   ifstream imesh;
+
+   imesh.open(destination_mesh_file);
+   if (imesh)
+   {
+      dest_mesh = make_shared<Mesh>(imesh, 1, 1);
+      imesh.close();
+   }
+   else
+   {
+      mfem::err << "WARNING: Destination mesh file not found: "
+                << destination_mesh_file << "\n"
+                << "Using default 2D quad mesh.";
+
+      dest_mesh = make_shared<Mesh>(4, 4, Element::QUADRILATERAL);
+   }
+
+   const int dim = dest_mesh->Dimension();
+
+   Vector box_min(dim), box_max(dim), range(dim);
+   dest_mesh->GetBoundingBox(box_min, box_max);
+   range = box_max;
+   range -= box_min;
+
+   imesh.open(source_mesh_file);
+
+   if (imesh)
+   {
+      src_mesh = make_shared<Mesh>(imesh, 1, 1);
+      imesh.close();
+   }
+   else
+   {
+      mfem::err << "WARNING: Source mesh file not found: " << source_mesh_file
+                << "\n"
+                << "Using default box mesh.\n";
+
+      if (dim == 2)
+      {
+         src_mesh =
+            make_shared<Mesh>(4, 4, Element::TRIANGLE, 1, range[0], range[1]);
+      }
+      else if (dim == 3)
+      {
+         src_mesh = make_shared<Mesh>(4, 4, 4, Element::TETRAHEDRON, 1, range[0],
+                                      range[1], range[2]);
+      }
+
+      for (int i = 0; i < src_mesh->GetNV(); ++i)
+      {
+         double *v = src_mesh->GetVertex(i);
+
+         for (int d = 0; d < dim; ++d)
+         {
+            v[d] += box_min[d];
+         }
+      }
+   }
+
+   for (int i = 0; i < src_n_refinements; ++i)
+   {
+      src_mesh->UniformRefinement();
+   }
+
+   for (int i = 0; i < dest_n_refinements; ++i)
+   {
+      dest_mesh->UniformRefinement();
+   }
+
+   shared_ptr<FiniteElementCollection> src_fe_coll, dest_fe_coll;
+
+   if (use_vector_fe)
+   {
+      src_fe_coll =
+         make_shared<RT_FECollection>(source_fe_order, src_mesh->Dimension());
+      dest_fe_coll =
+         make_shared<RT_FECollection>(dest_fe_order, dest_mesh->Dimension());
+   }
+   else
+   {
+      src_fe_coll =
+         make_shared<L2_FECollection>(source_fe_order, src_mesh->Dimension());
+      dest_fe_coll =
+         make_shared<L2_FECollection>(dest_fe_order, dest_mesh->Dimension());
+   }
+
+   auto src_fe =
+      make_shared<FiniteElementSpace>(src_mesh.get(), src_fe_coll.get());
+
+   auto dest_fe =
+      make_shared<FiniteElementSpace>(dest_mesh.get(), dest_fe_coll.get());
+
+   GridFunction src_fun(src_fe.get());
+   GridFunction dest_fun(dest_fe.get());
+   src_fun = 1.0;
+
+   // To be used with standard fe
+   FunctionCoefficient coeff(example_fun);
+
+   // To be used with vector fe
+   VectorFunctionCoefficient vector_coeff(dim, &vector_fun);
+
+   if (use_vector_fe)
+   {
+      src_fun.ProjectCoefficient(vector_coeff);
+      src_fun.Update();
+   }
+   else
+   {
+      src_fun.ProjectCoefficient(coeff);
+      src_fun.Update();
+   }
+
+   dest_fun = 0.0;
+   dest_fun.Update();
+
+   MortarAssembler assembler(src_fe, dest_fe);
+   assembler.SetVerbose(verbose);
+
+   if (use_vector_fe)
+   {
+      assembler.AddMortarIntegrator(make_shared<VectorL2MortarIntegrator>());
+   }
+   else
+   {
+      assembler.AddMortarIntegrator(make_shared<L2MortarIntegrator>());
+   }
+
+   if (assembler.Transfer(src_fun, dest_fun))
+   {
+      if (visualization)
+      {
+         dest_fun.Update();
+
+         double src_err = 0;
+         double dest_err = 0;
+
+         if (use_vector_fe)
+         {
+            src_err = src_fun.ComputeL2Error(vector_coeff);
+            dest_err = dest_fun.ComputeL2Error(vector_coeff);
+         }
+         else
+         {
+            src_err = src_fun.ComputeL2Error(coeff);
+            dest_err = dest_fun.ComputeL2Error(coeff);
+         }
+
+         mfem::out << "l2 error: src: " << src_err << ", dest: " << dest_err
+                   << std::endl;
+
+         plot(*src_mesh, src_fun, "source");
+         plot(*dest_mesh, dest_fun, "destination");
+      }
+   }
+   else
+   {
+      mfem::out << "No intersection -> no transfer!" << std::endl;
+   }
+
+   // Finalize transfer library context
+   return FinalizeTransfer();
+}
@@ -0,0 +1,258 @@
+//                   MFEM + Moonolith Example 1 (parallel version)
+//
+// Compile with: make ex1p
+//
+// Moonolith sample runs:
+//   mpirun -np 4 ex1p
+//   mpirun -np 4 ex1p --source_refinements 1 --dest_refinements 2
+//   mpirun -np 4 ex1p -s ../../data/inline-hex.mesh -d ../../data/inline-tet.mesh
+//
+// Description:  This example code demonstrates the use of MFEM for transferring
+//               discrete fields from one finite element mesh to another. The
+//               meshes can be of arbitrary shape and completely unrelated with
+//               each other. This feature can be used for implementing immersed
+//               domain methods for fluid-structure interaction or general
+//               multi-physics applications.
+//
+//               This particular example is for parallel runtimes. Vector FE is
+//               an experimental feature in parallel.
+
+#include "example_utils.hpp"
+#include "mfem.hpp"
+
+using namespace mfem;
+using namespace std;
+
+void destination_transform(const Vector &x, Vector &x_new)
+{
+   x_new = x;
+   // x_new *= 0.5;
+}
+
+int main(int argc, char *argv[])
+{
+   MPI_Init(&argc, &argv);
+
+   int num_procs, rank;
+
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+   // Init transfer library context, with MPI handled outside the library
+   InitTransfer(argc, argv, MPI_COMM_WORLD);
+
+   const char *source_mesh_file = "../../data/inline-tri.mesh";
+   const char *destination_mesh_file = "../../data/inline-quad.mesh";
+
+   int src_n_refinements = 0;
+   int dest_n_refinements = 0;
+   int source_fe_order = 1;
+   int dest_fe_order = 1;
+   bool visualization = true;
+   bool use_vector_fe = false;
+   bool verbose = false;
+   bool assemble_mass_and_coupling_together = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&source_mesh_file, "-s", "--source_mesh",
+                  "Mesh file to use for src.");
+   args.AddOption(&destination_mesh_file, "-d", "--destination_mesh",
+                  "Mesh file to use for dest.");
+   args.AddOption(&src_n_refinements, "-sr", "--source_refinements",
+                  "Number of src refinements");
+   args.AddOption(&dest_n_refinements, "-dr", "--dest_refinements",
+                  "Number of dest refinements");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&source_fe_order, "-so", "--source_fe_order",
+                  "Order of the src finite elements");
+   args.AddOption(&dest_fe_order, "-do", "--dest_fe_order",
+                  "Order of the dest finite elements");
+   args.AddOption(&verbose, "-verb", "--verbose", "--no-verb", "--no-verbose",
+                  "Enable/Disable verbose output");
+   args.AddOption(&use_vector_fe, "-vfe", "--use_vector_fe", "-no-vfe",
+                  "--no-vector_fe", "Use vector finite elements (Experimental)");
+   args.AddOption(&assemble_mass_and_coupling_together, "-act",
+                  "--assemble_mass_and_coupling_together", "-no-act",
+                  "--no-assemble_mass_and_coupling_together",
+                  "Assemble mass and coupling operators together (better for non-affine elements)");
+   args.Parse();
+   check_options(args);
+
+   shared_ptr<Mesh> src_mesh, dest_mesh;
+
+   ifstream imesh;
+
+   imesh.open(destination_mesh_file);
+   if (imesh)
+   {
+      dest_mesh = make_shared<Mesh>(imesh, 1, 1);
+      imesh.close();
+   }
+   else
+   {
+      if (rank == 0)
+         mfem::err << "WARNING: Destination mesh file not found: "
+                   << destination_mesh_file << "\n"
+                   << "Using default 2D quad mesh.";
+
+      dest_mesh = make_shared<Mesh>(4, 4, Element::QUADRILATERAL);
+   }
+
+   const int dim = dest_mesh->Dimension();
+
+   dest_mesh->Transform(&destination_transform);
+
+   Vector box_min(dim), box_max(dim), range(dim);
+   dest_mesh->GetBoundingBox(box_min, box_max);
+   range = box_max;
+   range -= box_min;
+
+   imesh.open(source_mesh_file);
+
+   if (imesh)
+   {
+      src_mesh = make_shared<Mesh>(imesh, 1, 1);
+      imesh.close();
+   }
+   else
+   {
+      if (rank == 0)
+         mfem::err << "WARNING: Source mesh file not found: " << source_mesh_file
+                   << "\n"
+                   << "Using default box mesh.\n";
+
+      if (dim == 2)
+      {
+         src_mesh =
+            make_shared<Mesh>(4, 4, Element::TRIANGLE, 1, range[0], range[1]);
+      }
+      else if (dim == 3)
+      {
+         src_mesh = make_shared<Mesh>(4, 4, 4, Element::TETRAHEDRON, 1, range[0],
+                                      range[1], range[2]);
+      }
+
+      for (int i = 0; i < src_mesh->GetNV(); ++i)
+      {
+         double *v = src_mesh->GetVertex(i);
+
+         for (int d = 0; d < dim; ++d)
+         {
+            v[d] += box_min[d];
+         }
+      }
+   }
+
+   for (int i = 0; i < src_n_refinements; ++i)
+   {
+      src_mesh->UniformRefinement();
+   }
+
+   for (int i = 0; i < dest_n_refinements; ++i)
+   {
+      dest_mesh->UniformRefinement();
+   }
+
+   auto p_src_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *src_mesh);
+   auto p_dest_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *dest_mesh);
+
+   shared_ptr<FiniteElementCollection> src_fe_coll, dest_fe_coll;
+
+   if (use_vector_fe)
+   {
+      src_fe_coll =
+         make_shared<RT_FECollection>(source_fe_order, src_mesh->Dimension());
+      dest_fe_coll =
+         make_shared<RT_FECollection>(dest_fe_order, dest_mesh->Dimension());
+   }
+   else
+   {
+      src_fe_coll =
+         make_shared<L2_FECollection>(source_fe_order, src_mesh->Dimension());
+      dest_fe_coll =
+         make_shared<L2_FECollection>(dest_fe_order, dest_mesh->Dimension());
+   }
+
+   auto src_fe =
+      make_shared<ParFiniteElementSpace>(p_src_mesh.get(), src_fe_coll.get());
+
+   auto dest_fe =
+      make_shared<ParFiniteElementSpace>(p_dest_mesh.get(), dest_fe_coll.get());
+
+   ParGridFunction src_fun(src_fe.get());
+
+   // To be used with standard fe
+   FunctionCoefficient coeff(example_fun);
+
+   // To be used with vector fe
+   VectorFunctionCoefficient vector_coeff(dim, &vector_fun);
+
+   if (use_vector_fe)
+   {
+      src_fun.ProjectCoefficient(vector_coeff);
+      src_fun.Update();
+   }
+   else
+   {
+      src_fun.ProjectCoefficient(coeff);
+      src_fun.Update();
+   }
+
+   ParGridFunction dest_fun(dest_fe.get());
+   dest_fun = 0.0;
+   dest_fun.Update();
+
+   ParMortarAssembler assembler(src_fe, dest_fe);
+   assembler.SetAssembleMassAndCouplingTogether(
+      assemble_mass_and_coupling_together);
+   assembler.SetVerbose(verbose);
+
+   if (use_vector_fe)
+   {
+      assembler.AddMortarIntegrator(make_shared<VectorL2MortarIntegrator>());
+   }
+   else
+   {
+      assembler.AddMortarIntegrator(make_shared<L2MortarIntegrator>());
+   }
+
+   if (assembler.Transfer(src_fun, dest_fun))
+   {
+
+      if (visualization)
+      {
+         double src_err = 0;
+         double dest_err = 0;
+
+         if (use_vector_fe)
+         {
+            src_err = src_fun.ComputeL2Error(vector_coeff);
+            dest_err = dest_fun.ComputeL2Error(vector_coeff);
+         }
+         else
+         {
+            src_err = src_fun.ComputeL2Error(coeff);
+            dest_err = dest_fun.ComputeL2Error(coeff);
+         }
+
+         if (rank == 0)
+         {
+            mfem::out << "l2 error: src: " << src_err << ", dest: " << dest_err
+                      << std::endl;
+         }
+
+         plot(*p_src_mesh, src_fun, "source");
+         plot(*p_dest_mesh, dest_fun, "destination");
+      }
+   }
+   else
+   {
+      mfem::out << "No intersection no transfer!" << std::endl;
+   }
+
+   // Finalize transfer library context
+   FinalizeTransfer();
+   return MPI_Finalize();
+}
@@ -0,0 +1,202 @@
+//                  MFEM + Moonolith Example 2 (parallel version)
+//
+// Compile with: make ex2p
+//
+// Moonolith sample runs:
+//   mpirun -np 4 ex2p
+//   mpirun -np 4 ex2p --source_refinements 1 --dest_refinements 2
+//   mpirun -np 4 ex2p -s ../../data/inline-hex.mesh -d ../../data/inline-tet.mesh
+//
+// Description:  This example code demonstrates the use of MFEM for transferring
+//               discrete fields from one finite element mesh to another. The
+//               meshes can be of arbitrary shape and completely unrelated with
+//               each other. This feature can be used for implementing immersed
+//               domain methods for fluid-structure interaction or general
+//               multi-physics applications.
+//
+//               This particular example concerns discontinuous Galerkin FEM with
+//               adaptive mesh refinement for parallel runtimes.
+
+#include "example_utils.hpp"
+#include "mfem.hpp"
+
+using namespace mfem;
+using namespace std;
+
+void destination_transform(const Vector &x, Vector &x_new)
+{
+   x_new = x;
+   // x_new *= .5;
+}
+
+int main(int argc, char *argv[])
+{
+   MPI_Init(&argc, &argv);
+
+   int num_procs, rank;
+
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+   // Init transfer library context, with MPI handled outside the library
+   InitTransfer(argc, argv, MPI_COMM_WORLD);
+
+   const char *source_mesh_file = "../../data/inline-tri.mesh";
+   const char *destination_mesh_file = "../../data/inline-quad.mesh";
+
+   int src_n_refinements = 0;
+   int dest_n_refinements = 0;
+
+   // Source fe order has to be greater or equal than destination order
+   int source_fe_order = 1;
+   int dest_fe_order = 0;
+   bool visualization = true;
+   bool verbose = false;
+   int max_iterations = 30000;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&source_mesh_file, "-s", "--source_mesh",
+                  "Mesh file to use for src.");
+   args.AddOption(&destination_mesh_file, "-d", "--destination_mesh",
+                  "Mesh file to use for dest.");
+   args.AddOption(&src_n_refinements, "-sr", "--source_refinements",
+                  "Number of src refinements");
+   args.AddOption(&dest_n_refinements, "-dr", "--dest_refinements",
+                  "Number of dest refinements");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&source_fe_order, "-so", "--source_fe_order",
+                  "Order of the src finite elements");
+   args.AddOption(&dest_fe_order, "-do", "--dest_fe_order",
+                  "Order of the dest finite elements");
+   args.AddOption(&verbose, "-verb", "--verbose", "--no-verb", "--no-verbose",
+                  "Enable/Disable verbose output");
+   args.AddOption(&max_iterations, "-m", "--max_iterations",
+                  "Max number of solver iterations");
+   args.Parse();
+   check_options(args);
+
+   if (source_fe_order == 0 &&  dest_fe_order != 0)
+   {
+      mfem::out <<
+                "Source fe order should not be 0 unless destination fe order is also 0!\n";
+
+      FinalizeTransfer();
+      return MPI_Finalize();
+   }
+
+   ifstream imesh(source_mesh_file);
+   shared_ptr<Mesh> src_mesh, dest_mesh;
+   if (imesh)
+   {
+      src_mesh = make_shared<Mesh>(imesh, 1, 1);
+      imesh.close();
+   }
+   else
+   {
+      if (rank == 0)
+         mfem::err << "WARNING: Source mesh file not found: " << source_mesh_file
+                   << "\n"
+                   << "Using default 2D triangle mesh.";
+      src_mesh = make_shared<Mesh>(4, 4, Element::TRIANGLE);
+   }
+
+   imesh.open(destination_mesh_file);
+   if (imesh)
+   {
+      dest_mesh = make_shared<Mesh>(imesh, 1, 1);
+      imesh.close();
+   }
+   else
+   {
+      if (rank == 0)
+         mfem::err << "WARNING: Destination mesh file not found: "
+                   << destination_mesh_file << "\n"
+                   << "Using default 2D quad mesh.";
+      dest_mesh = make_shared<Mesh>(4, 4, Element::QUADRILATERAL);
+   }
+
+   dest_mesh->Transform(&destination_transform);
+
+   for (int i = 0; i < src_n_refinements; ++i)
+   {
+      src_mesh->UniformRefinement();
+   }
+
+   for (int i = 0; i < dest_n_refinements; ++i)
+   {
+      dest_mesh->UniformRefinement();
+   }
+
+   src_mesh->EnsureNCMesh();
+   dest_mesh->EnsureNCMesh();
+
+   {
+      for (int l = 0; l < 4; l++)
+      {
+         src_mesh->RandomRefinement(0.1); // 10% probability
+      }
+   }
+
+   {
+      for (int l = 0; l < 4; l++)
+      {
+         dest_mesh->RandomRefinement(0.1); // 10% probability
+      }
+   }
+
+   auto p_src_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *src_mesh);
+   auto p_dest_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *dest_mesh);
+
+   auto src_fe_coll =
+      make_shared<DG_FECollection>(source_fe_order, p_src_mesh->Dimension());
+   auto src_fe =
+      make_shared<ParFiniteElementSpace>(p_src_mesh.get(), src_fe_coll.get());
+
+   auto dest_fe_coll =
+      make_shared<DG_FECollection>(dest_fe_order, p_dest_mesh->Dimension());
+   auto dest_fe =
+      make_shared<ParFiniteElementSpace>(p_dest_mesh.get(), dest_fe_coll.get());
+
+   ParGridFunction src_fun(src_fe.get());
+   FunctionCoefficient coeff(example_fun);
+   make_fun(*src_fe, coeff, src_fun);
+
+   ParGridFunction dest_fun(dest_fe.get());
+   dest_fun = 0.0;
+   dest_fun.Update();
+
+   ParMortarAssembler assembler(src_fe, dest_fe);
+   assembler.SetVerbose(verbose);
+   assembler.SetMaxSolverIterations(max_iterations);
+
+   assembler.AddMortarIntegrator(make_shared<L2MortarIntegrator>());
+   if (assembler.Transfer(src_fun, dest_fun))
+   {
+      if (visualization)
+      {
+
+         const double src_err = src_fun.ComputeL2Error(coeff);
+         const double dest_err = dest_fun.ComputeL2Error(coeff);
+
+         if (rank == 0)
+         {
+            mfem::out << "l2 error: src: " << src_err << ", dest: " << dest_err
+                      << std::endl;
+         }
+
+         plot(*p_src_mesh, src_fun, "source");
+         plot(*p_dest_mesh, dest_fun, "destination");
+      }
+   }
+   else
+   {
+      mfem::out << "Transfer failed! Use --verbose option for diagnostic!" <<
+                std::endl;
+   }
+
+   // Finalize transfer library context
+   FinalizeTransfer();
+   return MPI_Finalize();
+}
@@ -0,0 +1,96 @@
+#include <algorithm>
+#include <assert.h>
+#include <cstdlib>
+#include <memory>
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif // MFEM_USE_MPI
+
+#include "mfem.hpp"
+
+inline void check_options(mfem::OptionsParser &args)
+{
+   using namespace std;
+   using namespace mfem;
+
+   int rank = 0;
+#ifdef MFEM_USE_MPI
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+#endif // MFEM_USE_MPI
+
+   if (!args.Good())
+   {
+      if (rank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+
+#ifdef MFEM_USE_MPI
+      MPI_Finalize();
+      MPI_Abort(MPI_COMM_WORLD, 1);
+#else
+      abort();
+#endif // MFEM_USE_MPI
+   }
+
+   if (rank == 0)
+   {
+      args.PrintOptions(cout);
+   }
+}
+
+inline void make_fun(mfem::FiniteElementSpace &fe, mfem::Coefficient &c,
+                     mfem::GridFunction &f)
+{
+   using namespace std;
+   using namespace mfem;
+
+   f.SetSpace(&fe);
+   f.ProjectCoefficient(c);
+   f.Update();
+}
+
+inline double example_fun(const mfem::Vector &x)
+{
+   using namespace std;
+   using namespace mfem;
+
+   const int n = x.Size();
+   double ret = 0;
+   for (int k = 0; k < n; ++k)
+   {
+      ret += x(k) * x(k);
+   }
+
+   return sqrt(ret);
+}
+
+void vector_fun(const mfem::Vector &x, mfem::Vector &f)
+{
+   const double n = x.Norml2();
+   f.SetSize(x.Size());
+   f = n;
+}
+
+inline void plot(mfem::Mesh &mesh, mfem::GridFunction &x, std::string title)
+{
+   using namespace std;
+   using namespace mfem;
+
+   int num_procs = 1, rank = 0;
+
+#ifdef MFEM_USE_MPI
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+#endif // MFEM_USE_MPI
+
+   char vishost[] = "localhost";
+   int visport = 19916;
+   socketstream sol_sock(vishost, visport);
+   sol_sock << "parallel " << num_procs << " " << rank << "\n";
+   sol_sock.precision(8);
+   sol_sock << "solution\n" << mesh << x
+            << "window_title '"<< title << "'\n" << flush;
+   sol_sock << flush;
+}
@@ -0,0 +1,78 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/moonolith/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex1
+PAR_EXAMPLES = ex1p ex2p
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+ifeq ($(MFEM_USE_MOONOLITH),NO)
+$(EXAMPLES):
+	$(error MFEM is not configured with MOONOLITH)
+endif
+
+MFEM_TESTS = EXAMPLES
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI_NP = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP)
+RUN_MPI = $(RUN_MPI_NP) $(MFEM_MPI_NP)
+SERIAL_NAME := Serial MOONOLITH example
+PARALLEL_NAME := Parallel MOONOLITH example
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_NAME))
+%-test-seq: %
+	@$(call mfem-test,$<,, $(SERIAL_NAME))
+
+# Testing: Example-specific execution options:
+ex1-test-par: ex1
+	@$(call mfem-test,$<, $(RUN_MPI_NP) 1, $(PARALLEL_NAME))
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
@@ -58,17 +58,29 @@ set(SRCS
  gridfunc.cpp
  hybridization.cpp
  intrules.cpp
-  ceed/algebraic.cpp
-  ceed/full-assembly.cpp
-  ceed/solvers-atpmg.cpp
-  ceed/convection.cpp
-  ceed/diffusion.cpp
-  ceed/nlconvection.cpp
-  ceed/mass.cpp
-  ceed/operator.cpp
-  ceed/util.cpp
+  ceed/interface/basis.cpp
+  ceed/interface/restriction.cpp
+  ceed/interface/operator.cpp
+  ceed/interface/util.cpp
+  ceed/integrators/convection/convection.cpp
+  ceed/integrators/diffusion/diffusion.cpp
+  ceed/integrators/nlconvection/nlconvection.cpp
+  ceed/integrators/mass/mass.cpp
+  ceed/solvers/algebraic.cpp
+  ceed/solvers/full-assembly.cpp
+  ceed/solvers/solvers-atpmg.cpp
  linearform.cpp
+  linearform_ext.cpp
  lininteg.cpp
+  lininteg_domain.cpp
+  lininteg_domain_grad.cpp
+  lor/lor.cpp
+  lor/lor_ads.cpp
+  lor/lor_ams.cpp
+  lor/lor_batched.cpp
+  lor/lor_h1.cpp
+  lor/lor_nd.cpp
+  lor/lor_rt.cpp
  multigrid.cpp
  nonlinearform.cpp
  nonlinearform_ext.cpp
@@ -118,7 +130,6 @@ set(SRCS
  tmop_amr.cpp
  gslib.cpp
  transfer.cpp
-  lor.cpp
  )

 set(HDRS
@@ -149,19 +160,31 @@ set(HDRS
  gridfunc.hpp
  hybridization.hpp
  intrules.hpp
-  ceed/algebraic.hpp
-  ceed/full-assembly.hpp
-  ceed/solvers-atpmg.hpp
-  ceed/coefficient.hpp
-  ceed/convection.hpp
-  ceed/diffusion.hpp
-  ceed/integrator.hpp
-  ceed/mass.hpp
-  ceed/nlconvection.hpp
-  ceed/operator.hpp
-  ceed/util.hpp
+  ceed/interface/basis.hpp
+  ceed/interface/integrator.hpp
+  ceed/interface/interface.hpp
+  ceed/interface/operator.hpp
+  ceed/interface/restriction.hpp
+  ceed/interface/util.hpp
+  ceed/integrators/convection/convection.hpp
+  ceed/integrators/diffusion/diffusion.hpp
+  ceed/integrators/mass/mass.hpp
+  ceed/integrators/nlconvection/nlconvection.hpp
+  ceed/interface/coefficient.hpp
+  ceed/solvers/algebraic.hpp
+  ceed/solvers/full-assembly.hpp
+  ceed/solvers/solvers-atpmg.hpp
  linearform.hpp
+  linearform_ext.hpp
  lininteg.hpp
+  lor/lor.hpp
+  lor/lor_ads.hpp
+  lor/lor_ams.hpp
+  lor/lor_batched.hpp
+  lor/lor_h1.hpp
+  lor/lor_nd.hpp
+  lor/lor_rt.hpp
+  lor/lor_util.hpp
  multigrid.hpp
  nonlinearform.hpp
  nonlinearform_ext.hpp
@@ -188,7 +211,6 @@ set(HDRS
  tmop_amr.hpp
  gslib.hpp
  transfer.hpp
-  lor.hpp
  )

 if (MFEM_USE_SIDRE)
@@ -1763,12 +1763,24 @@ void MixedBilinearForm::FormRectangularSystemMatrix(

   mat->Finalize();

-   if (test_P) // TODO: Must actually check for trial_P too
+   if (test_P && trial_P)
   {
      SparseMatrix *m = RAP(*test_P, *mat, *trial_P);
      delete mat;
      mat = m;
   }
+   else if (test_P)
+   {
+      SparseMatrix *m = TransposeMult(*test_P, *mat);
+      delete mat;
+      mat = m;
+   }
+   else if (trial_P)
+   {
+      SparseMatrix *m = mfem::Mult(*mat, *trial_P);
+      delete mat;
+      mat = m;
+   }

   Array<int> ess_trial_tdof_marker, ess_test_tdof_marker;
   FiniteElementSpace::ListToMarker(trial_tdof_list, trial_fes->GetTrueVSize(),
@@ -16,7 +16,7 @@
 #include "bilinearform.hpp"
 #include "pbilinearform.hpp"
 #include "pgridfunc.hpp"
-#include "ceed/util.hpp"
+#include "ceed/interface/util.hpp"

 namespace mfem
 {
@@ -903,7 +903,8 @@ void FABilinearFormExtension::Assemble()
   }
   else // We create, compute the sparsity, and fill the sparse matrix
   {
-      mat = new SparseMatrix(height, width, 0);
+      mat = new SparseMatrix;
+      mat->OverrideSize(height, width);
      if (fes.IsDGSpace())
      {
         const L2ElementRestriction *restE =
@@ -1036,7 +1036,8 @@ void DiffusionIntegrator::AssembleElementVector(

 void DiffusionIntegrator::ComputeElementFlux
 ( const FiniteElement &el, ElementTransformation &Trans,
-  Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef )
+  Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef,
+  const IntegrationRule *ir)
 {
   int nd, spaceDim, fnd;

@@ -1057,8 +1058,6 @@ void DiffusionIntegrator::ComputeElementFlux
                  "Unexpected height for MatrixCoefficient");
   }

-   MFEM_VERIFY(!SMQ, "SymmetricMatrixCoefficient not supported here");
-
 #ifdef MFEM_THREAD_SAFE
   DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim);
   DenseMatrix M(MQ ? spaceDim : 0);
@@ -1073,13 +1072,16 @@ void DiffusionIntegrator::ComputeElementFlux
   vecdxt.SetSize(spaceDim);
   pointflux.SetSize(MQ || VQ ? spaceDim : 0);

-   const IntegrationRule &ir = fluxelem.GetNodes();
-   fnd = ir.GetNPoints();
+   if (!ir)
+   {
+      ir = &fluxelem.GetNodes();
+   }
+   fnd = ir->GetNPoints();
   flux.SetSize( fnd * spaceDim );

   for (int i = 0; i < fnd; i++)
   {
-      const IntegrationPoint &ip = ir.IntPoint(i);
+      const IntegrationPoint &ip = ir->IntPoint(i);
      el.CalcDShape(ip, dshape);
      dshape.MultTranspose(u, vec);

@@ -1146,8 +1148,6 @@ double DiffusionIntegrator::ComputeFluxEnergy
   D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif

-   MFEM_VERIFY(!SMQ, "SymmetricMatrixCoefficient not supported here");
-
   shape.SetSize(nd);
   pointflux.SetSize(spaceDim);
   if (d_energy) { vec.SetSize(spaceDim); }
@@ -1672,7 +1672,8 @@ void VectorFEDivergenceIntegrator::AssembleElementMatrix2(
   {
      const IntegrationPoint &ip = ir->IntPoint(i);
      trial_fe.CalcDivShape(ip, divshape);
-      test_fe.CalcShape(ip, shape);
+      Trans.SetIntPoint(&ip);
+      test_fe.CalcPhysShape(Trans, shape);
      double w = ip.weight;
      if (Q)
      {
@@ -2005,12 +2006,14 @@ void CurlCurlIntegrator::AssembleElementMatrix
 void CurlCurlIntegrator
 ::ComputeElementFlux(const FiniteElement &el, ElementTransformation &Trans,
                     Vector &u, const FiniteElement &fluxelem, Vector &flux,
-                     bool with_coef)
+                     bool with_coef, const IntegrationRule *ir)
 {
 #ifdef MFEM_THREAD_SAFE
   DenseMatrix projcurl;
 #endif

+   MFEM_VERIFY(ir == NULL, "Integration rule (ir) must be NULL")
+
   fluxelem.ProjectCurl(el, Trans, projcurl);

   flux.SetSize(projcurl.Height());
@@ -2842,7 +2845,7 @@ void ElasticityIntegrator::AssembleElementMatrix(
 void ElasticityIntegrator::ComputeElementFlux(
   const mfem::FiniteElement &el, ElementTransformation &Trans,
   Vector &u, const mfem::FiniteElement &fluxelem, Vector &flux,
-   bool with_coef)
+   bool with_coef, const IntegrationRule *ir)
 {
   const int dof = el.GetDof();
   const int dim = el.GetDim();
@@ -2865,14 +2868,17 @@ void ElasticityIntegrator::ComputeElementFlux(
   DenseMatrix gh(gh_data, dim, dim);
   DenseMatrix grad(grad_data, dim, dim);

-   const IntegrationRule &ir = fluxelem.GetNodes();
-   const int fnd = ir.GetNPoints();
+   if (!ir)
+   {
+      ir = &fluxelem.GetNodes();
+   }
+   const int fnd = ir->GetNPoints();
   flux.SetSize(fnd * tdim);

   DenseMatrix loc_data_mat(u.GetData(), dof, dim);
   for (int i = 0; i < fnd; i++)
   {
-      const IntegrationPoint &ip = ir.IntPoint(i);
+      const IntegrationPoint &ip = ir->IntPoint(i);
      el.CalcDShape(ip, dshape);
      MultAtB(loc_data_mat, dshape, gh);

@@ -15,6 +15,7 @@
 #include "../config/config.hpp"
 #include "nonlininteg.hpp"
 #include "fespace.hpp"
+#include "ceed/interface/util.hpp"

 namespace mfem
 {
@@ -213,12 +214,18 @@ public:
                             of the method may choose not to scale the "flux"
                             function by any coefficients describing the
                             integrator.
+       @param[in] ir  If passed (the default value is NULL), the implementation
+                        of the method will ignore the integration rule provided
+                        by the @a fluxelem parameter and, instead, compute the
+                        discrete flux at the points specified by the integration
+                        rule @a ir.
    */
   virtual void ComputeElementFlux(const FiniteElement &el,
                                   ElementTransformation &Trans,
                                   Vector &u,
                                   const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef = true) { }
+                                   Vector &flux, bool with_coef = true,
+                                   const IntegrationRule *ir = NULL) { }

   /** @brief Virtual method required for Zienkiewicz-Zhu type error estimators.

@@ -1899,6 +1906,7 @@ protected:
                           const FiniteElementSpace &test_fes);

   virtual void AddMultPA(const Vector&, Vector&) const;
+   virtual void AddMultTransposePA(const Vector&, Vector&) const;

 private:
   // PA extension
@@ -1957,6 +1965,7 @@ protected:
                           const FiniteElementSpace &test_fes);

   virtual void AddMultPA(const Vector&, Vector&) const;
+   virtual void AddMultTransposePA(const Vector&, Vector&) const;

 private:
   // PA extension
@@ -2068,7 +2077,6 @@ protected:
   Coefficient *Q;
   VectorCoefficient *VQ;
   MatrixCoefficient *MQ;
-   SymmetricMatrixCoefficient *SMQ;

 private:
   Vector vec, vecdxt, pointflux, shape;
@@ -2090,30 +2098,24 @@ public:
   /// Construct a diffusion integrator with coefficient Q = 1
   DiffusionIntegrator(const IntegrationRule *ir = nullptr)
      : BilinearFormIntegrator(ir),
-        Q(NULL), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
+        Q(NULL), VQ(NULL), MQ(NULL), maps(NULL), geom(NULL) { }

   /// Construct a diffusion integrator with a scalar coefficient q
   DiffusionIntegrator(Coefficient &q, const IntegrationRule *ir = nullptr)
      : BilinearFormIntegrator(ir),
-        Q(&q), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
+        Q(&q), VQ(NULL), MQ(NULL), maps(NULL), geom(NULL) { }

   /// Construct a diffusion integrator with a vector coefficient q
   DiffusionIntegrator(VectorCoefficient &q,
                       const IntegrationRule *ir = nullptr)
      : BilinearFormIntegrator(ir),
-        Q(NULL), VQ(&q), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
+        Q(NULL), VQ(&q), MQ(NULL), maps(NULL), geom(NULL) { }

   /// Construct a diffusion integrator with a matrix coefficient q
   DiffusionIntegrator(MatrixCoefficient &q,
                       const IntegrationRule *ir = nullptr)
      : BilinearFormIntegrator(ir),
-        Q(NULL), VQ(NULL), MQ(&q), SMQ(NULL), maps(NULL), geom(NULL) { }
-
-   /// Construct a diffusion integrator with a symmetric matrix coefficient q
-   DiffusionIntegrator(SymmetricMatrixCoefficient &q,
-                       const IntegrationRule *ir = nullptr)
-      : BilinearFormIntegrator(ir),
-        Q(NULL), VQ(NULL), MQ(NULL), SMQ(&q), maps(NULL), geom(NULL) { }
+        Q(NULL), VQ(NULL), MQ(&q), maps(NULL), geom(NULL) { }

   /** Given a particular Finite Element computes the element stiffness matrix
       elmat. */
@@ -2135,7 +2137,8 @@ public:
   virtual void ComputeElementFlux(const FiniteElement &el,
                                   ElementTransformation &Trans,
                                   Vector &u, const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef = true);
+                                   Vector &flux, bool with_coef = true,
+                                   const IntegrationRule *ir = NULL);

   virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                    ElementTransformation &Trans,
@@ -2164,6 +2167,8 @@ public:
                                         const FiniteElement &test_fe);

   bool SupportsCeed() const { return DeviceCanUseCeed(); }
+
+   Coefficient *GetCoefficient() const { return Q; }
 };

 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */
@@ -2223,6 +2228,8 @@ public:
                                         ElementTransformation &Trans);

   bool SupportsCeed() const { return DeviceCanUseCeed(); }
+
+   const Coefficient *GetCoefficient() const { return Q; }
 };

 /** Mass integrator (u, v) restricted to the boundary of a domain */
@@ -2391,8 +2398,9 @@ public:
    scalar function given by FiniteElement through standard transformation.
    Here, u is the trial function and p is the test function.

-    Note: the element matrix returned by AssembleElementMatrix2 does NOT depend
-    on the ElementTransformation Trans. */
+    Note: if the test space does not have map type INTEGRAL, then the element
+    matrix returned by AssembleElementMatrix2 will not depend on the
+    ElementTransformation Trans. */
 class VectorFEDivergenceIntegrator : public BilinearFormIntegrator
 {
 protected:
@@ -2523,7 +2531,6 @@ protected:
   Coefficient *Q;
   DiagonalMatrixCoefficient *DQ;
   MatrixCoefficient *MQ;
-   SymmetricMatrixCoefficient *SMQ;

   // PA extension
   Vector pa_data;
@@ -2534,18 +2541,15 @@ protected:
   bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient

 public:
-   CurlCurlIntegrator() { Q = NULL; DQ = NULL; MQ = NULL; SMQ = NULL; }
+   CurlCurlIntegrator() { Q = NULL; DQ = NULL; MQ = NULL; }
   /// Construct a bilinear form integrator for Nedelec elements
   CurlCurlIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) :
-      BilinearFormIntegrator(ir), Q(&q), DQ(NULL), MQ(NULL), SMQ(NULL) { }
+      BilinearFormIntegrator(ir), Q(&q), DQ(NULL), MQ(NULL) { }
   CurlCurlIntegrator(DiagonalMatrixCoefficient &dq,
                      const IntegrationRule *ir = NULL) :
-      BilinearFormIntegrator(ir), Q(NULL), DQ(&dq), MQ(NULL), SMQ(NULL) { }
+      BilinearFormIntegrator(ir), Q(NULL), DQ(&dq), MQ(NULL) { }
   CurlCurlIntegrator(MatrixCoefficient &mq, const IntegrationRule *ir = NULL) :
-      BilinearFormIntegrator(ir), Q(NULL), DQ(NULL), MQ(&mq), SMQ(NULL) { }
-   CurlCurlIntegrator(SymmetricMatrixCoefficient &smq,
-                      const IntegrationRule *ir = NULL) :
-      BilinearFormIntegrator(ir), Q(NULL), DQ(NULL), MQ(NULL), SMQ(&smq) { }
+      BilinearFormIntegrator(ir), Q(NULL), DQ(NULL), MQ(&mq) { }

   /* Given a particular Finite Element, compute the
      element curl-curl matrix elmat */
@@ -2556,7 +2560,8 @@ public:
   virtual void ComputeElementFlux(const FiniteElement &el,
                                   ElementTransformation &Trans,
                                   Vector &u, const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef);
+                                   Vector &flux, bool with_coef,
+                                   const IntegrationRule *ir = NULL);

   virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                    ElementTransformation &Trans,
@@ -2566,6 +2571,8 @@ public:
   virtual void AssemblePA(const FiniteElementSpace &fes);
   virtual void AddMultPA(const Vector &x, Vector &y) const;
   virtual void AssembleDiagonalPA(Vector& diag);
+
+   const Coefficient *GetCoefficient() const { return Q; }
 };

 /** Integrator for (curl u, curl v) for FE spaces defined by 'dim' copies of a
@@ -2602,9 +2609,8 @@ public:
 class VectorFEMassIntegrator: public BilinearFormIntegrator
 {
 private:
-   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq,
-             SymmetricMatrixCoefficient *smq)
-   { Q = q; DQ = dq; MQ = mq; SMQ = smq; }
+   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq)
+   { Q = q; DQ = dq; MQ = mq; }

 #ifndef MFEM_THREAD_SAFE
   Vector shape;
@@ -2619,7 +2625,6 @@ protected:
   Coefficient *Q;
   DiagonalMatrixCoefficient *DQ;
   MatrixCoefficient *MQ;
-   SymmetricMatrixCoefficient *SMQ;

   // PA extension
   Vector pa_data;
@@ -2632,15 +2637,13 @@ protected:
   bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient

 public:
-   VectorFEMassIntegrator() { Init(NULL, NULL, NULL, NULL); }
-   VectorFEMassIntegrator(Coefficient *q_) { Init(q_, NULL, NULL, NULL); }
-   VectorFEMassIntegrator(Coefficient &q) { Init(&q, NULL, NULL, NULL); }
-   VectorFEMassIntegrator(DiagonalMatrixCoefficient *dq_) { Init(NULL, dq_, NULL, NULL); }
-   VectorFEMassIntegrator(DiagonalMatrixCoefficient &dq) { Init(NULL, &dq, NULL, NULL); }
-   VectorFEMassIntegrator(MatrixCoefficient *mq_) { Init(NULL, NULL, mq_, NULL); }
-   VectorFEMassIntegrator(MatrixCoefficient &mq) { Init(NULL, NULL, &mq, NULL); }
-   VectorFEMassIntegrator(SymmetricMatrixCoefficient &smq) { Init(NULL, NULL, NULL, &smq); }
-   VectorFEMassIntegrator(SymmetricMatrixCoefficient *smq) { Init(NULL, NULL, NULL, smq); }
+   VectorFEMassIntegrator() { Init(NULL, NULL, NULL); }
+   VectorFEMassIntegrator(Coefficient *q_) { Init(q_, NULL, NULL); }
+   VectorFEMassIntegrator(Coefficient &q) { Init(&q, NULL, NULL); }
+   VectorFEMassIntegrator(DiagonalMatrixCoefficient *dq_) { Init(NULL, dq_, NULL); }
+   VectorFEMassIntegrator(DiagonalMatrixCoefficient &dq) { Init(NULL, &dq, NULL); }
+   VectorFEMassIntegrator(MatrixCoefficient *mq_) { Init(NULL, NULL, mq_); }
+   VectorFEMassIntegrator(MatrixCoefficient &mq) { Init(NULL, NULL, &mq); }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
@@ -2655,7 +2658,10 @@ public:
   virtual void AssemblePA(const FiniteElementSpace &trial_fes,
                           const FiniteElementSpace &test_fes);
   virtual void AddMultPA(const Vector &x, Vector &y) const;
+   virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
   virtual void AssembleDiagonalPA(Vector& diag);
+
+   const Coefficient *GetCoefficient() const { return Q; }
 };

 /** Integrator for (Q div u, p) where u=(v1,...,vn) and all vi are in the same
@@ -2736,6 +2742,7 @@ public:
   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
+   const Coefficient *GetCoefficient() const { return Q; }
 };

 /** Integrator for
@@ -2875,12 +2882,14 @@ public:
       of the stress components is: s_xx, s_yy, s_xy. In 3D, it is: s_xx, s_yy,
       s_zz, s_xy, s_xz, s_yz. In other words, @a flux is the local vector for
       a FE space with dim*(dim+1)/2 vector components, based on the finite
-       element @a fluxelem. */
+       element @a fluxelem. The integration rule is taken from @a fluxelem.
+       @a ir exists to specific an alternative integration rule. */
   virtual void ComputeElementFlux(const FiniteElement &el,
                                   ElementTransformation &Trans,
                                   Vector &u,
                                   const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef = true);
+                                   Vector &flux, bool with_coef = true,
+                                   const IntegrationRule *ir = NULL);

   /** Compute the element energy (integral of the strain energy density)
       corresponding to the stress represented by @a flux which is a vector of
@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/convection.hpp"
+#include "ceed/integrators/convection/convection.hpp"

 using namespace std;

@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/convection.hpp"
+#include "ceed/integrators/convection/convection.hpp"
 #include "quadinterpolator.hpp"

 namespace mfem
@@ -43,9 +43,10 @@ static void PADGTraceSetup2D(const int Q1D,
   auto W = w.Read();
   auto qd = Reshape(op.Write(), Q1D, 2, 2, NF);

-   MFEM_FORALL(f, NF, // can be optimized with Q1D thread for NF blocks
+   MFEM_FORALL(tid, Q1D*NF,
   {
-      for (int q = 0; q < Q1D; ++q)
+      const int f = tid / Q1D;
+      const int q = tid % Q1D;
      {
         const double r = const_r ? R(0,0) : R(q,f);
         const double v0 = const_v ? V(0,0,0) : V(0,q,f);
@@ -85,11 +86,12 @@ static void PADGTraceSetup3D(const int Q1D,
   auto W = w.Read();
   auto qd = Reshape(op.Write(), Q1D, Q1D, 2, 2, NF);

-   MFEM_FORALL(f, NF, // can be optimized with Q1D*Q1D threads for NF blocks
+   MFEM_FORALL(tid, Q1D*Q1D*NF,
   {
-      for (int q1 = 0; q1 < Q1D; ++q1)
+      int f = tid / (Q1D * Q1D);
+      int q2 = (tid / Q1D) % Q1D;
+      int q1 = tid % Q1D;
      {
-         for (int q2 = 0; q2 < Q1D; ++q2)
         {
            const double r = const_r ? R(0,0,0) : R(q1,q2,f);
            const double v0 = const_v ? V(0,0,0,0) : V(0,q1,q2,f);
@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/diffusion.hpp"
+#include "ceed/integrators/diffusion/diffusion.hpp"

 using namespace std;

@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/diffusion.hpp"
+#include "ceed/integrators/diffusion/diffusion.hpp"

 using namespace std;

@@ -362,7 +362,7 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
   if (DeviceCanUseCeed())
   {
      delete ceedOp;
-      MFEM_VERIFY(!VQ && !MQ && !SMQ,
+      MFEM_VERIFY(!VQ && !MQ,
                  "Only scalar coefficient supported for DiffusionIntegrator"
                  " with libCEED");
      ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
@@ -381,41 +381,14 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
   int coeffDim = 1;
   Vector coeff;
   const int MQfullDim = MQ ? MQ->GetHeight() * MQ->GetWidth() : 0;
-   if (MQ)
-   {
-      symmetric = false;
-      MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
-
-      coeffDim = MQfullDim;
-
-      coeff.SetSize(MQfullDim * nq * ne);
-
-      DenseMatrix GM;
-      GM.SetSize(dim);
-
-      auto C = Reshape(coeff.HostWrite(), MQfullDim, nq, ne);
-      for (int e=0; e<ne; ++e)
-      {
-         ElementTransformation *tr = mesh->GetElementTransformation(e);
-         for (int p=0; p<nq; ++p)
-         {
-            MQ->Eval(GM, *tr, ir->IntPoint(p));
-            for (int i=0; i<dim; ++i)
-               for (int j=0; j<dim; ++j)
-               {
-                  C(j+(i*dim), p, e) = GM(i,j);
-               }
-         }
-      }
-   }
-   else if (SMQ)
+   if (auto *SMQ = dynamic_cast<SymmetricMatrixCoefficient *>(MQ))
   {
      MFEM_VERIFY(SMQ->GetSize() == dim, "");
      coeffDim = symmDims;
      coeff.SetSize(symmDims * nq * ne);

-      DenseSymmetricMatrix SM;
-      SM.SetSize(dim);
+      DenseSymmetricMatrix sym_mat;
+      sym_mat.SetSize(dim);

      auto C = Reshape(coeff.HostWrite(), symmDims, nq, ne);

@@ -424,12 +397,39 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
         ElementTransformation *tr = mesh->GetElementTransformation(e);
         for (int p=0; p<nq; ++p)
         {
-            SMQ->Eval(SM, *tr, ir->IntPoint(p));
+            SMQ->Eval(sym_mat, *tr, ir->IntPoint(p));
            int cnt = 0;
            for (int i=0; i<dim; ++i)
               for (int j=i; j<dim; ++j, ++cnt)
               {
-                  C(cnt, p, e) = SM(i,j);
+                  C(cnt, p, e) = sym_mat(i,j);
+               }
+         }
+      }
+   }
+   else if (MQ)
+   {
+      symmetric = false;
+      MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
+
+      coeffDim = MQfullDim;
+
+      coeff.SetSize(MQfullDim * nq * ne);
+
+      DenseMatrix mat;
+      mat.SetSize(dim);
+
+      auto C = Reshape(coeff.HostWrite(), MQfullDim, nq, ne);
+      for (int e=0; e<ne; ++e)
+      {
+         ElementTransformation *tr = mesh->GetElementTransformation(e);
+         for (int p=0; p<nq; ++p)
+         {
+            MQ->Eval(mat, *tr, ir->IntPoint(p));
+            for (int i=0; i<dim; ++i)
+               for (int j=0; j<dim; ++j)
+               {
+                  C(j+(i*dim), p, e) = mat(i,j);
               }
         }
      }
@@ -978,12 +978,14 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)

   MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");

+   auto SMQ = dynamic_cast<SymmetricMatrixCoefficient *>(MQ);
+
   const int MQsymmDim = SMQ ? (SMQ->GetSize() * (SMQ->GetSize() + 1)) / 2 : 0;
   const int MQfullDim = MQ ? (MQ->GetHeight() * MQ->GetWidth()) : 0;
-   const int MQdim = MQ ? MQfullDim : MQsymmDim;
-   const int coeffDim = (MQ || SMQ) ? MQdim : (DQ ? DQ->GetVDim() : 1);
+   const int MQdim = SMQ ? MQsymmDim : MQfullDim;
+   const int coeffDim = MQ ? MQdim : (DQ ? DQ->GetVDim() : 1);

-   symmetric = (MQ == NULL);
+   symmetric = (SMQ || MQ == NULL);

   const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
   const int ndata = (dim == 2) ? 1 : (symmetric ? symmDims : MQfullDim);
@@ -992,7 +994,7 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
   Vector coeff(coeffDim * ne * nq);
   coeff = 1.0;
   auto coeffh = Reshape(coeff.HostWrite(), coeffDim, nq, ne);
-   if (Q || DQ || MQ || SMQ)
+   if (Q || DQ || MQ)
   {
      Vector DM(DQ ? coeffDim : 0);
      DenseMatrix GM;
@@ -1002,35 +1004,24 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
      {
         MFEM_VERIFY(coeffDim == dimc, "");
      }
-      if (MQ)
-      {
-         GM.SetSize(dimc);
-         MFEM_VERIFY(coeffDim == MQdim, "");
-         MFEM_VERIFY(MQ->GetHeight() == dimc && MQ->GetWidth() == dimc, "");
-      }
      if (SMQ)
      {
         SM.SetSize(dimc);
         MFEM_VERIFY(SMQ->GetSize() == dimc, "");
      }
+      else if (MQ)
+      {
+         GM.SetSize(dimc);
+         MFEM_VERIFY(coeffDim == MQdim, "");
+         MFEM_VERIFY(MQ->GetHeight() == dimc && MQ->GetWidth() == dimc, "");
+      }

      for (int e=0; e<ne; ++e)
      {
         ElementTransformation *tr = mesh->GetElementTransformation(e);
         for (int p=0; p<nq; ++p)
         {
-            if (MQ)
-            {
-               MQ->Eval(GM, *tr, ir->IntPoint(p));
-
-               for (int i=0; i<dimc; ++i)
-                  for (int j=0; j<dimc; ++j)
-                  {
-                     coeffh(j+(i*dimc), p, e) = GM(i,j);
-                  }
-
-            }
-            else if (SMQ)
+            if (SMQ)
            {
               SMQ->Eval(SM, *tr, ir->IntPoint(p));

@@ -1041,6 +1032,17 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
                     coeffh(cnt, p, e) = SM(i,j);
                  }

+            }
+            else if (MQ)
+            {
+               MQ->Eval(GM, *tr, ir->IntPoint(p));
+
+               for (int i=0; i<dimc; ++i)
+                  for (int j=0; j<dimc; ++j)
+                  {
+                     coeffh(j+(i*dimc), p, e) = GM(i,j);
+                  }
+
            }
            else if (DQ)
            {
@@ -3633,8 +3635,8 @@ void MixedVectorCurlIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
      }
      else
      {
-         PAHcurlHdivSetup3D(quad1D, coeffDim, ne, false, ir->GetWeights(), geom->J,
-                            coeff, pa_data);
+         PAHcurlHdivSetup3D(quad1D, coeffDim, ne, false, ir->GetWeights(),
+                            geom->J, coeff, pa_data);
      }
   }
   else if (testType == mfem::FiniteElement::DIV &&
@@ -4679,6 +4681,362 @@ static void PAHcurlHdivApply3D(const int D1D,
   }); // end of element loop
 }

+// Apply to x corresponding to DOF's in H(div) (test), integrated against the
+// curl of H(curl) trial functions corresponding to y.
+template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
+static void PAHcurlHdivApply3DTranspose(const int D1D,
+                                        const int D1Dtest,
+                                        const int Q1D,
+                                        const int NE,
+                                        const Array<double> &bo,
+                                        const Array<double> &bc,
+                                        const Array<double> &bot,
+                                        const Array<double> &bct,
+                                        const Array<double> &gct,
+                                        const Vector &pa_data,
+                                        const Vector &x,
+                                        Vector &y)
+{
+   MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
+   MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
+   // Using Piola transformations (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u}
+   // for u in H(curl) and w = (1 / det (dF)) dF \hat{w} for w in H(div), we get
+   // (\nabla\times u) \cdot w = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{w}
+   // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
+   // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
+   // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
+
+   constexpr static int VDIM = 3;
+
+   auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
+   auto Bc = Reshape(bc.Read(), Q1D, D1D);
+   auto Bot = Reshape(bot.Read(), D1Dtest-1, Q1D);
+   auto Bct = Reshape(bct.Read(), D1Dtest, Q1D);
+   auto Gct = Reshape(gct.Read(), D1D, Q1D);
+   auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, 6, NE);
+   auto X = Reshape(x.Read(), 3*(D1Dtest-1)*(D1Dtest-1)*D1D, NE);
+   auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
+
+   MFEM_FORALL(e, NE,
+   {
+      double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];  // Assuming HDIV_MAX_D1D <= HCURL_MAX_D1D
+
+      for (int qz = 0; qz < Q1D; ++qz)
+      {
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int c = 0; c < VDIM; ++c)
+               {
+                  mass[qz][qy][qx][c] = 0.0;
+               }
+            }
+         }
+      }
+
+      int osc = 0;
+
+      for (int c = 0; c < VDIM; ++c)  // loop over x, y, z components
+      {
+         const int D1Dz = (c == 2) ? D1D : D1D - 1;
+         const int D1Dy = (c == 1) ? D1D : D1D - 1;
+         const int D1Dx = (c == 0) ? D1D : D1D - 1;
+
+         for (int dz = 0; dz < D1Dz; ++dz)
+         {
+            double massXY[HDIV_MAX_Q1D][HDIV_MAX_Q1D];
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  massXY[qy][qx] = 0.0;
+               }
+            }
+
+            for (int dy = 0; dy < D1Dy; ++dy)
+            {
+               double massX[HDIV_MAX_Q1D];
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  massX[qx] = 0.0;
+               }
+
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
+                  for (int qx = 0; qx < Q1D; ++qx)
+                  {
+                     massX[qx] += t * ((c == 0) ? Bc(qx,dx) : Bo(qx,dx));
+                  }
+               }
+
+               for (int qy = 0; qy < Q1D; ++qy)
+               {
+                  const double wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
+                  for (int qx = 0; qx < Q1D; ++qx)
+                  {
+                     const double wx = massX[qx];
+                     massXY[qy][qx] += wx * wy;
+                  }
+               }
+            }
+
+            for (int qz = 0; qz < Q1D; ++qz)
+            {
+               const double wz = (c == 2) ? Bc(qz,dz) : Bo(qz,dz);
+               for (int qy = 0; qy < Q1D; ++qy)
+               {
+                  for (int qx = 0; qx < Q1D; ++qx)
+                  {
+                     mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
+                  }
+               }
+            }
+         }
+
+         osc += D1Dx * D1Dy * D1Dz;
+      }  // loop (c) over components
+
+      // Apply D operator.
+      for (int qz = 0; qz < Q1D; ++qz)
+      {
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               const double O11 = op(qx,qy,qz,0,e);
+               const double O12 = op(qx,qy,qz,1,e);
+               const double O13 = op(qx,qy,qz,2,e);
+               const double O22 = op(qx,qy,qz,3,e);
+               const double O23 = op(qx,qy,qz,4,e);
+               const double O33 = op(qx,qy,qz,5,e);
+               const double massX = mass[qz][qy][qx][0];
+               const double massY = mass[qz][qy][qx][1];
+               const double massZ = mass[qz][qy][qx][2];
+               mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
+               mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
+               mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
+            }
+         }
+      }
+
+      // x component
+      osc = 0;
+      {
+         const int D1Dz = D1D;
+         const int D1Dy = D1D;
+         const int D1Dx = D1D - 1;
+
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            double gradXY12[MAX_D1D][MAX_D1D];
+            double gradXY21[MAX_D1D][MAX_D1D];
+
+            for (int dy = 0; dy < D1Dy; ++dy)
+            {
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  gradXY12[dy][dx] = 0.0;
+                  gradXY21[dy][dx] = 0.0;
+               }
+            }
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               double massX[MAX_D1D][2];
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  for (int n = 0; n < 2; ++n)
+                  {
+                     massX[dx][n] = 0.0;
+                  }
+               }
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  for (int dx = 0; dx < D1Dx; ++dx)
+                  {
+                     const double wx = Bot(dx,qx);
+
+                     massX[dx][0] += wx * mass[qz][qy][qx][1];
+                     massX[dx][1] += wx * mass[qz][qy][qx][2];
+                  }
+               }
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  const double wy = Bct(dy,qy);
+                  const double wDy = Gct(dy,qy);
+
+                  for (int dx = 0; dx < D1Dx; ++dx)
+                  {
+                     gradXY21[dy][dx] += massX[dx][0] * wy;
+                     gradXY12[dy][dx] += massX[dx][1] * wDy;
+                  }
+               }
+            }
+
+            for (int dz = 0; dz < D1Dz; ++dz)
+            {
+               const double wz = Bct(dz,qz);
+               const double wDz = Gct(dz,qz);
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  for (int dx = 0; dx < D1Dx; ++dx)
+                  {
+                     // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
+                     // (u_0)_{x_2} * (op * curl)_1 - (u_0)_{x_1} * (op * curl)_2
+                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                       e) += (gradXY21[dy][dx] * wDz) - (gradXY12[dy][dx] * wz);
+                  }
+               }
+            }
+         }  // loop qz
+
+         osc += D1Dx * D1Dy * D1Dz;
+      }
+
+      // y component
+      {
+         const int D1Dz = D1D;
+         const int D1Dy = D1D - 1;
+         const int D1Dx = D1D;
+
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            double gradXY02[MAX_D1D][MAX_D1D];
+            double gradXY20[MAX_D1D][MAX_D1D];
+
+            for (int dy = 0; dy < D1Dy; ++dy)
+            {
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  gradXY02[dy][dx] = 0.0;
+                  gradXY20[dy][dx] = 0.0;
+               }
+            }
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               double massY[MAX_D1D][2];
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  massY[dy][0] = 0.0;
+                  massY[dy][1] = 0.0;
+               }
+               for (int qy = 0; qy < Q1D; ++qy)
+               {
+                  for (int dy = 0; dy < D1Dy; ++dy)
+                  {
+                     const double wy = Bot(dy,qy);
+
+                     massY[dy][0] += wy * mass[qz][qy][qx][2];
+                     massY[dy][1] += wy * mass[qz][qy][qx][0];
+                  }
+               }
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  const double wx = Bct(dx,qx);
+                  const double wDx = Gct(dx,qx);
+
+                  for (int dy = 0; dy < D1Dy; ++dy)
+                  {
+                     gradXY02[dy][dx] += massY[dy][0] * wDx;
+                     gradXY20[dy][dx] += massY[dy][1] * wx;
+                  }
+               }
+            }
+
+            for (int dz = 0; dz < D1Dz; ++dz)
+            {
+               const double wz = Bct(dz,qz);
+               const double wDz = Gct(dz,qz);
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  for (int dx = 0; dx < D1Dx; ++dx)
+                  {
+                     // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
+                     // -(u_1)_{x_2} * (op * curl)_0 + (u_1)_{x_0} * (op * curl)_2
+                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                       e) += (-gradXY20[dy][dx] * wDz) + (gradXY02[dy][dx] * wz);
+                  }
+               }
+            }
+         }  // loop qz
+
+         osc += D1Dx * D1Dy * D1Dz;
+      }
+
+      // z component
+      {
+         const int D1Dz = D1D - 1;
+         const int D1Dy = D1D;
+         const int D1Dx = D1D;
+
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            double gradYZ01[MAX_D1D][MAX_D1D];
+            double gradYZ10[MAX_D1D][MAX_D1D];
+
+            for (int dy = 0; dy < D1Dy; ++dy)
+            {
+               for (int dz = 0; dz < D1Dz; ++dz)
+               {
+                  gradYZ01[dz][dy] = 0.0;
+                  gradYZ10[dz][dy] = 0.0;
+               }
+            }
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               double massZ[MAX_D1D][2];
+               for (int dz = 0; dz < D1Dz; ++dz)
+               {
+                  for (int n = 0; n < 2; ++n)
+                  {
+                     massZ[dz][n] = 0.0;
+                  }
+               }
+               for (int qz = 0; qz < Q1D; ++qz)
+               {
+                  for (int dz = 0; dz < D1Dz; ++dz)
+                  {
+                     const double wz = Bot(dz,qz);
+
+                     massZ[dz][0] += wz * mass[qz][qy][qx][0];
+                     massZ[dz][1] += wz * mass[qz][qy][qx][1];
+                  }
+               }
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  const double wy = Bct(dy,qy);
+                  const double wDy = Gct(dy,qy);
+
+                  for (int dz = 0; dz < D1Dz; ++dz)
+                  {
+                     gradYZ01[dz][dy] += wy * massZ[dz][1];
+                     gradYZ10[dz][dy] += wDy * massZ[dz][0];
+                  }
+               }
+            }
+
+            for (int dx = 0; dx < D1Dx; ++dx)
+            {
+               const double wx = Bct(dx,qx);
+               const double wDx = Gct(dx,qx);
+
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  for (int dz = 0; dz < D1Dz; ++dz)
+                  {
+                     // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
+                     // (u_2)_{x_1} * (op * curl)_0 - (u_2)_{x_0} * (op * curl)_1
+                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                       e) += (gradYZ10[dz][dy] * wx) - (gradYZ01[dz][dy] * wDx);
+                  }
+               }
+            }
+         }  // loop qx
+      }
+   }); // end of element loop
+}
+
 void MixedVectorCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
   if (testType == mfem::FiniteElement::CURL &&
@@ -4723,6 +5081,20 @@ void MixedVectorCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
   }
 }

+void MixedVectorCurlIntegrator::AddMultTransposePA(const Vector &x,
+                                                   Vector &y) const
+{
+   if (testType == mfem::FiniteElement::DIV &&
+       trialType == mfem::FiniteElement::CURL && dim == 3)
+      PAHcurlHdivApply3DTranspose(dofs1D, dofs1Dtest, quad1D, ne, mapsO->B,
+                                  mapsC->B, mapsOtest->Bt, mapsCtest->Bt,
+                                  mapsC->Gt, pa_data, x, y);
+   else
+   {
+      MFEM_ABORT("Unsupported dimension or space!");
+   }
+}
+
 void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
                                               &trial_fes,
                                               const FiniteElementSpace &test_fes)
@@ -4761,8 +5133,16 @@ void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace

   MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");

+   testType = test_el->GetDerivType();
+   trialType = trial_el->GetDerivType();
+
+   const bool curlSpaces = (testType == mfem::FiniteElement::CURL &&
+                            trialType == mfem::FiniteElement::CURL);
+
+   const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
+
   coeffDim = DQ ? 3 : 1;
-   const int ndata = DQ ? 9 : 1;
+   const int ndata = curlSpaces ? (DQ ? 9 : 1) : symmDims;

   pa_data.SetSize(ndata * nq * ne, Device::GetMemoryType());

@@ -4799,9 +5179,6 @@ void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
      }
   }

-   testType = test_el->GetDerivType();
-   trialType = trial_el->GetDerivType();
-
   if (trialType == mfem::FiniteElement::CURL && dim == 3)
   {
      if (coeffDim == 1)
@@ -4810,10 +5187,16 @@ void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
      }
      else
      {
-         PAHcurlHdivSetup3D(quad1D, coeffDim, ne, false, ir->GetWeights(), geom->J,
-                            coeff, pa_data);
+         PAHcurlHdivSetup3D(quad1D, coeffDim, ne, false, ir->GetWeights(),
+                            geom->J, coeff, pa_data);
      }
   }
+   else if (trialType == mfem::FiniteElement::DIV && dim == 3 &&
+            test_el->GetOrder() == trial_el->GetOrder())
+   {
+      PACurlCurlSetup3D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J, coeff,
+                        pa_data);
+   }
   else
   {
      MFEM_ABORT("Unknown kernel.");
@@ -5463,6 +5846,29 @@ void MixedVectorWeakCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
         PAHcurlL2Apply3DTranspose(dofs1D, quad1D, ndata, ne, mapsO->B,
                                   mapsC->B, mapsO->Bt, mapsC->Bt, mapsC->Gt, pa_data, x, y);
   }
+   else if (testType == mfem::FiniteElement::CURL &&
+            trialType == mfem::FiniteElement::DIV && dim == 3)
+   {
+      PAHcurlHdivApply3DTranspose(dofs1D, dofs1D, quad1D, ne, mapsO->B,
+                                  mapsC->B, mapsO->Bt, mapsC->Bt,
+                                  mapsC->Gt, pa_data, x, y);
+   }
+   else
+   {
+      MFEM_ABORT("Unsupported dimension or space!");
+   }
+}
+
+void MixedVectorWeakCurlIntegrator::AddMultTransposePA(const Vector &x,
+                                                       Vector &y) const
+{
+   if (testType == mfem::FiniteElement::CURL &&
+       trialType == mfem::FiniteElement::DIV && dim == 3)
+   {
+      PAHcurlHdivApply3D(dofs1D, dofs1D, quad1D, ne, mapsO->B,
+                         mapsC->B, mapsO->Bt, mapsC->Bt, mapsC->G,
+                         pa_data, x, y);
+   }
   else
   {
      MFEM_ABORT("Unsupported dimension or space!");
@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/mass.hpp"
+#include "ceed/integrators/mass/mass.hpp"

 using namespace std;

@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/mass.hpp"
+#include "ceed/integrators/mass/mass.hpp"

 using namespace std;

@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/diffusion.hpp"
+#include "ceed/integrators/diffusion/diffusion.hpp"

 using namespace std;

@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/diffusion.hpp"
+#include "ceed/integrators/diffusion/diffusion.hpp"

 using namespace std;

@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/mass.hpp"
+#include "ceed/integrators/mass/mass.hpp"

 using namespace std;

@@ -12,7 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
-#include "ceed/mass.hpp"
+#include "ceed/integrators/mass/mass.hpp"

 using namespace std;

@@ -372,8 +372,8 @@ void PAHcurlHdivSetup2D(const int Q1D,
               const double R21 = D2*J21;
               const double R22 = D2*J22;
               y(i11,qx,qy,e) = w_detJ * ( J22*R11 - J12*R21); // 1,1
-               y(i12,qx,qy,e) = w_detJ * ( J22*R12 - J12*R22); // 1,2
-               y(i21,qx,qy,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1
+               y(i21,qx,qy,e) = w_detJ * ( J22*R12 - J12*R22); // 1,2 (transpose)
+               y(i12,qx,qy,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1 (transpose)
               y(i22,qx,qy,e) = w_detJ * (-J21*R12 + J11*R22); // 2,2
            }
         }
@@ -389,6 +389,7 @@ void PAHcurlHdivMassApply3D(const int D1D,
                            const int NE,
                            const bool scalarCoeff,
                            const bool trialHcurl,
+                            const bool transpose,
                            const Array<double> &Bo_,
                            const Array<double> &Bc_,
                            const Array<double> &Bot_,
@@ -413,6 +414,13 @@ void PAHcurlHdivMassApply3D(const int D1D,
   auto y = Reshape(y_.ReadWrite(), 3*(D1Dtest-1)*D1Dtest*
                    (trialHcurl ? D1Dtest-1 : D1Dtest), NE);

+   const int i12 = transpose ? 3 : 1;
+   const int i13 = transpose ? 6 : 2;
+   const int i21 = transpose ? 1 : 3;
+   const int i23 = transpose ? 7 : 5;
+   const int i31 = transpose ? 2 : 6;
+   const int i32 = transpose ? 5 : 7;
+
   MFEM_FORALL(e, NE,
   {
      double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
@@ -507,13 +515,13 @@ void PAHcurlHdivMassApply3D(const int D1D,
            for (int qx = 0; qx < Q1D; ++qx)
            {
               const double O11 = op(0,qx,qy,qz,e);
-               const double O12 = scalarCoeff ? 0.0 : op(1,qx,qy,qz,e);
-               const double O13 = scalarCoeff ? 0.0 : op(2,qx,qy,qz,e);
-               const double O21 = scalarCoeff ? 0.0 : op(3,qx,qy,qz,e);
+               const double O12 = scalarCoeff ? 0.0 : op(i12,qx,qy,qz,e);
+               const double O13 = scalarCoeff ? 0.0 : op(i13,qx,qy,qz,e);
+               const double O21 = scalarCoeff ? 0.0 : op(i21,qx,qy,qz,e);
               const double O22 = scalarCoeff ? O11 : op(4,qx,qy,qz,e);
-               const double O23 = scalarCoeff ? 0.0 : op(5,qx,qy,qz,e);
-               const double O31 = scalarCoeff ? 0.0 : op(6,qx,qy,qz,e);
-               const double O32 = scalarCoeff ? 0.0 : op(7,qx,qy,qz,e);
+               const double O23 = scalarCoeff ? 0.0 : op(i23,qx,qy,qz,e);
+               const double O31 = scalarCoeff ? 0.0 : op(i31,qx,qy,qz,e);
+               const double O32 = scalarCoeff ? 0.0 : op(i32,qx,qy,qz,e);
               const double O33 = scalarCoeff ? O11 : op(8,qx,qy,qz,e);
               const double massX = mass[qz][qy][qx][0];
               const double massY = mass[qz][qy][qx][1];
@@ -601,6 +609,7 @@ void PAHcurlHdivMassApply2D(const int D1D,
                            const int NE,
                            const bool scalarCoeff,
                            const bool trialHcurl,
+                            const bool transpose,
                            const Array<double> &Bo_,
                            const Array<double> &Bc_,
                            const Array<double> &Bot_,
@@ -624,6 +633,9 @@ void PAHcurlHdivMassApply2D(const int D1D,
   auto x = Reshape(x_.Read(), 2*(D1D-1)*D1D, NE);
   auto y = Reshape(y_.ReadWrite(), 2*(D1Dtest-1)*D1Dtest, NE);

+   const int i12 = transpose ? 2 : 1;
+   const int i21 = transpose ? 1 : 2;
+
   MFEM_FORALL(e, NE,
   {
      double mass[MAX_Q1D][MAX_Q1D][VDIM];
@@ -685,8 +697,8 @@ void PAHcurlHdivMassApply2D(const int D1D,
         for (int qx = 0; qx < Q1D; ++qx)
         {
            const double O11 = op(0,qx,qy,e);
-            const double O12 = scalarCoeff ? 0.0 : op(1,qx,qy,e);
-            const double O21 = scalarCoeff ? 0.0 : op(2,qx,qy,e);
+            const double O12 = scalarCoeff ? 0.0 : op(i12,qx,qy,e);
+            const double O21 = scalarCoeff ? 0.0 : op(i21,qx,qy,e);
            const double O22 = scalarCoeff ? O11 : op(3,qx,qy,e);
            const double massX = mass[qy][qx][0];
            const double massY = mass[qy][qx][1];
@@ -785,12 +797,14 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
   trial_fetype = trial_el->GetDerivType();
   test_fetype = test_el->GetDerivType();

+   auto SMQ = dynamic_cast<SymmetricMatrixCoefficient *>(MQ);
+
   const int MQsymmDim = SMQ ? (SMQ->GetSize() * (SMQ->GetSize() + 1)) / 2 : 0;
   const int MQfullDim = MQ ? (MQ->GetHeight() * MQ->GetWidth()) : 0;
-   const int MQdim = MQ ? MQfullDim : MQsymmDim;
-   const int coeffDim = (MQ || SMQ) ? MQdim : (DQ ? DQ->GetVDim() : 1);
+   const int MQdim = SMQ ? MQsymmDim : MQfullDim;
+   const int coeffDim = MQ ? MQdim : (DQ ? DQ->GetVDim() : 1);

-   symmetric = (MQ == NULL);
+   symmetric = (SMQ || MQ == NULL);

   const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
   const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
@@ -807,7 +821,7 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
   Vector coeff(coeffDim * ne * nq);
   coeff = 1.0;
   auto coeffh = Reshape(coeff.HostWrite(), coeffDim, nq, ne);
-   if (Q || DQ || MQ || SMQ)
+   if (Q || DQ || MQ)
   {
      Vector DM(DQ ? coeffDim : 0);
      DenseMatrix M;
@@ -817,34 +831,25 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
      {
         MFEM_VERIFY(coeffDim == dim, "");
      }
-      if (MQ)
-      {
-         MFEM_VERIFY(coeffDim == MQdim, "");
-         MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
-         M.SetSize(dim);
-      }
      if (SMQ)
      {
         MFEM_VERIFY(SMQ->GetSize() == dim, "");
         SM.SetSize(dim);
      }
+      else if (MQ)
+      {
+         MFEM_VERIFY(coeffDim == MQdim, "");
+         MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
+         M.SetSize(dim);
+      }
+

      for (int e=0; e<ne; ++e)
      {
         ElementTransformation *tr = mesh->GetElementTransformation(e);
         for (int p=0; p<nq; ++p)
         {
-            if (MQ)
-            {
-               MQ->Eval(M, *tr, ir->IntPoint(p));
-
-               for (int i=0; i<dim; ++i)
-                  for (int j=0; j<dim; ++j)
-                  {
-                     coeffh(j+(i*dim), p, e) = M(i,j);
-                  }
-            }
-            else if (SMQ)
+            if (SMQ)
            {
               SMQ->Eval(SM, *tr, ir->IntPoint(p));
               int cnt = 0;
@@ -854,6 +859,16 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
                     coeffh(cnt, p, e) = SM(i,j);
                  }
            }
+            else if (MQ)
+            {
+               MQ->Eval(M, *tr, ir->IntPoint(p));
+
+               for (int i=0; i<dim; ++i)
+                  for (int j=0; j<dim; ++j)
+                  {
+                     coeffh(j+(i*dim), p, e) = M(i,j);
+                  }
+            }
            else if (DQ)
            {
               DQ->Eval(DM, *tr, ir->IntPoint(p));
@@ -1024,16 +1039,16 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
      }
      else if (trial_curl && test_div)
      {
-         const bool scalarCoeff = !(DQ || MQ || SMQ);
+         const bool scalarCoeff = !(DQ || MQ);
         PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                true, mapsO->B, mapsC->B, mapsOtest->Bt,
+                                true, false, mapsO->B, mapsC->B, mapsOtest->Bt,
                                mapsCtest->Bt, pa_data, x, y);
      }
      else if (trial_div && test_curl)
      {
-         const bool scalarCoeff = !(DQ || MQ || SMQ);
+         const bool scalarCoeff = !(DQ || MQ);
         PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                false, mapsO->B, mapsC->B, mapsOtest->Bt,
+                                false, false, mapsO->B, mapsC->B, mapsOtest->Bt,
                                mapsCtest->Bt, pa_data, x, y);
      }
      else
@@ -1055,10 +1070,10 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
      }
      else if ((trial_curl && test_div) || (trial_div && test_curl))
      {
-         const bool scalarCoeff = !(DQ || MQ || SMQ);
+         const bool scalarCoeff = !(DQ || MQ);
         PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                trial_curl, mapsO->B, mapsC->B, mapsOtest->Bt,
-                                mapsCtest->Bt, pa_data, x, y);
+                                trial_curl, false, mapsO->B, mapsC->B,
+                                mapsOtest->Bt, mapsCtest->Bt, pa_data, x, y);
      }
      else
      {
@@ -1067,6 +1082,39 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
   }
 }

+void VectorFEMassIntegrator::AddMultTransposePA(const Vector &x,
+                                                Vector &y) const
+{
+   const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
+   const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
+   const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
+   const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
+
+   bool symmetricSpaces = true;
+
+   if (dim == 3 && ((trial_div && test_curl) || (trial_curl && test_div)))
+   {
+      const bool scalarCoeff = !(DQ || MQ);
+      PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                             trial_div, true, mapsO->B, mapsC->B, mapsOtest->Bt,
+                             mapsCtest->Bt, pa_data, x, y);
+      symmetricSpaces = false;
+   }
+   else if (dim == 2 && ((trial_curl && test_div) || (trial_div && test_curl)))
+   {
+      const bool scalarCoeff = !(DQ || MQ);
+      PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                             !trial_curl, true, mapsO->B, mapsC->B, mapsOtest->Bt,
+                             mapsCtest->Bt, pa_data, x, y);
+      symmetricSpaces = false;
+   }
+
+   if (symmetricSpaces)
+   {
+      this->AddMultPA(x, y);
+   }
+}
+
 void MixedVectorGradientIntegrator::AssemblePA(const FiniteElementSpace
                                               &trial_fes,
                                               const FiniteElementSpace &test_fes)
@@ -11,7 +11,7 @@

 #include "convection.hpp"

-#include "../../config/config.hpp"
+#include "../../../../config/config.hpp"
 #ifdef MFEM_USE_CEED
 #include "convection_qf.h"
 #endif
@@ -28,7 +28,7 @@ struct ConvectionOperatorInfo : public OperatorInfo
   ConvectionContext ctx;
   ConvectionOperatorInfo(int dim, double alpha)
   {
-      header = "/convection_qf.h";
+      header = "/integrators/convection/convection_qf.h";
      build_func_const = ":f_build_conv_const";
      build_qf_const = &f_build_conv_const;
      build_func_quad = ":f_build_conv_quad";
@@ -12,8 +12,8 @@
 #ifndef MFEM_LIBCEED_CONV_HPP
 #define MFEM_LIBCEED_CONV_HPP

-#include "integrator.hpp"
-#include "../fespace.hpp"
+#include "../../interface/integrator.hpp"
+#include "../../../fespace.hpp"

 namespace mfem
 {
@@ -11,7 +11,7 @@

 #include "diffusion.hpp"

-#include "../../config/config.hpp"
+#include "../../../../config/config.hpp"
 #ifdef MFEM_USE_CEED
 #include "diffusion_qf.h"
 #endif
@@ -28,7 +28,7 @@ struct DiffusionOperatorInfo : public OperatorInfo
   DiffusionContext ctx;
   DiffusionOperatorInfo(int dim)
   {
-      header = "/diffusion_qf.h";
+      header = "/integrators/diffusion/diffusion_qf.h";
      build_func_const = ":f_build_diff_const";
      build_qf_const = &f_build_diff_const;
      build_func_quad = ":f_build_diff_quad";
@@ -12,8 +12,8 @@
 #ifndef MFEM_LIBCEED_DIFF_HPP
 #define MFEM_LIBCEED_DIFF_HPP

-#include "integrator.hpp"
-#include "../fespace.hpp"
+#include "../../interface/integrator.hpp"
+#include "../../../fespace.hpp"

 namespace mfem
 {
@@ -11,7 +11,7 @@

 #include "mass.hpp"

-#include "../../config/config.hpp"
+#include "../../../../config/config.hpp"
 #ifdef MFEM_USE_CEED
 #include "mass_qf.h"
 #endif
@@ -28,7 +28,7 @@ struct MassOperatorInfo : public OperatorInfo
   MassContext ctx;
   MassOperatorInfo()
   {
-      header = "/mass_qf.h";
+      header = "/integrators/mass/mass_qf.h";
      build_func_const = ":f_build_mass_const";
      build_qf_const = &f_build_mass_const;
      build_func_quad = ":f_build_mass_quad";
@@ -12,8 +12,8 @@
 #ifndef MFEM_LIBCEED_MASS_HPP
 #define MFEM_LIBCEED_MASS_HPP

-#include "integrator.hpp"
-#include "../fespace.hpp"
+#include "../../interface/integrator.hpp"
+#include "../../../fespace.hpp"

 namespace mfem
 {
@@ -11,7 +11,7 @@

 #include "nlconvection.hpp"

-#include "../../config/config.hpp"
+#include "../../../../config/config.hpp"
 #ifdef MFEM_USE_CEED
 #include "nlconvection_qf.h"
 #endif
@@ -28,7 +28,7 @@ struct NLConvectionOperatorInfo : public OperatorInfo
   NLConvectionContext ctx;
   NLConvectionOperatorInfo(int dim)
   {
-      header = "/nlconvection_qf.h";
+      header = "/integrators/nlconvection/nlconvection_qf.h";
      build_func_const = ":f_build_conv_const";
      build_qf_const = &f_build_conv_const;
      build_func_quad = ":f_build_conv_quad";
@@ -12,8 +12,8 @@
 #ifndef MFEM_LIBCEED_NLCONV_HPP
 #define MFEM_LIBCEED_NLCONV_HPP

-#include "integrator.hpp"
-#include "../fespace.hpp"
+#include "../../interface/integrator.hpp"
+#include "../../../fespace.hpp"

 namespace mfem
 {
@@ -0,0 +1,134 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../../gridfunc.hpp"
+#include "util.hpp"
+
+namespace mfem
+{
+
+namespace ceed
+{
+
+#ifdef MFEM_USE_CEED
+
+static CeedElemTopology GetCeedTopology(Geometry::Type geom)
+{
+   switch (geom)
+   {
+      case Geometry::SEGMENT:
+         return CEED_TOPOLOGY_LINE;
+      case Geometry::TRIANGLE:
+         return CEED_TOPOLOGY_TRIANGLE;
+      case Geometry::SQUARE:
+         return CEED_TOPOLOGY_QUAD;
+      case Geometry::TETRAHEDRON:
+         return CEED_TOPOLOGY_TET;
+      case Geometry::CUBE:
+         return CEED_TOPOLOGY_HEX;
+      case Geometry::PRISM:
+         return CEED_TOPOLOGY_PRISM;
+      default:
+         MFEM_ABORT("This type of element is not supported");
+         return CEED_TOPOLOGY_PRISM; // Silence warning
+   }
+}
+
+static void InitNonTensorBasis(const mfem::FiniteElementSpace &fes,
+                               const mfem::IntegrationRule &ir,
+                               Ceed ceed, CeedBasis *basis)
+{
+   const mfem::DofToQuad &maps = fes.GetFE(0)->
+                                 GetDofToQuad(ir,mfem::DofToQuad::FULL);
+   mfem::Mesh *mesh = fes.GetMesh();
+   const int dim = mesh->Dimension();
+   const int ndofs = maps.ndof;
+   const int nqpts = maps.nqpt;
+   mfem::DenseMatrix qX(dim,nqpts);
+   mfem::Vector qW(nqpts);
+   for (int i = 0; i < nqpts; i++)
+   {
+      const mfem::IntegrationPoint &ip = ir.IntPoint(i);
+      qX(0,i) = ip.x;
+      if (dim>1) { qX(1,i) = ip.y; }
+      if (dim>2) { qX(2,i) = ip.z; }
+      qW(i) = ip.weight;
+   }
+   CeedBasisCreateH1(ceed, GetCeedTopology(fes.GetFE(0)->GetGeomType()),
+                     fes.GetVDim(), ndofs, nqpts,
+                     maps.Bt.GetData(), maps.Gt.GetData(),
+                     qX.GetData(), qW.GetData(), basis);
+}
+
+static void InitTensorBasis(const mfem::FiniteElementSpace &fes,
+                            const mfem::IntegrationRule &ir,
+                            Ceed ceed, CeedBasis *basis)
+{
+   const mfem::DofToQuad &maps =
+      fes.GetFE(0)->GetDofToQuad(ir, mfem::DofToQuad::TENSOR);
+   mfem::Mesh *mesh = fes.GetMesh();
+   const int ndofs = maps.ndof;
+   const int nqpts = maps.nqpt;
+   mfem::Vector qX(nqpts), qW(nqpts);
+   // The x-coordinates of the first `nqpts` points of the integration rule are
+   // the points of the corresponding 1D rule. We also scale the weights
+   // accordingly.
+   double w_sum = 0.0;
+   for (int i = 0; i < nqpts; i++)
+   {
+      const mfem::IntegrationPoint &ip = ir.IntPoint(i);
+      qX(i) = ip.x;
+      qW(i) = ip.weight;
+      w_sum += ip.weight;
+   }
+   qW *= 1.0/w_sum;
+   CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes.GetVDim(), ndofs,
+                           nqpts, maps.Bt.GetData(),
+                           maps.Gt.GetData(), qX.GetData(),
+                           qW.GetData(), basis);
+}
+
+void InitBasis(const FiniteElementSpace &fes,
+               const IntegrationRule &irm,
+               Ceed ceed, CeedBasis *basis)
+{
+   // Check for FES -> basis, restriction in hash tables
+   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const int P = fe->GetDof();
+   const int Q = irm.GetNPoints();
+   const int ncomp = fes.GetVDim();
+   BasisKey basis_key(&fes, &irm, ncomp, P, Q);
+   auto basis_itr = mfem::internal::ceed_basis_map.find(basis_key);
+
+   // Init or retreive key values
+   if (basis_itr == mfem::internal::ceed_basis_map.end())
+   {
+      if (UsesTensorBasis(fes))
+      {
+         InitTensorBasis(fes, irm, ceed, basis);
+      }
+      else
+      {
+         InitNonTensorBasis(fes, irm, ceed, basis);
+      }
+      mfem::internal::ceed_basis_map[basis_key] = *basis;
+   }
+   else
+   {
+      *basis = basis_itr->second;
+   }
+}
+
+#endif
+
+} // namespace ceed
+
+} // namespace mfem
@@ -0,0 +1,42 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_LIBCEED_BASIS
+#define MFEM_LIBCEED_BASIS
+
+#include "ceed.hpp"
+
+namespace mfem
+{
+
+namespace ceed
+{
+
+#ifdef MFEM_USE_CEED
+
+/** @brief Initialize a CeedBasis.
+
+   @param[in] fes Input finite element space.
+   @param[in] irm Input integration rule.
+   @param[in] ceed Input Ceed object.
+   @param[out] basis The address of the initialized CeedBasis object.
+*/
+void InitBasis(const FiniteElementSpace &fes,
+               const IntegrationRule &irm,
+               Ceed ceed, CeedBasis *basis);
+
+#endif
+
+} // namespace ceed
+
+} // namespace mfem
+
+#endif // MFEM_LIBCEED_BASIS
@@ -0,0 +1,35 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_LIBCEED_CEED
+#define MFEM_LIBCEED_CEED
+
+#ifdef MFEM_USE_CEED
+#include <ceed.h>
+#if !CEED_VERSION_GE(0,10,0)
+#error MFEM requires a libCEED version >= 0.10.0
+#endif
+
+namespace mfem
+{
+
+namespace internal
+{
+
+extern Ceed ceed;
+
+} // namespace internal
+
+} // namespace mfem
+
+#endif // MFEM_USE_CEED
+
+#endif // MFEM_LIBCEED_CEED
@@ -12,14 +12,15 @@
 #ifndef MFEM_LIBCEED_COEFF
 #define MFEM_LIBCEED_COEFF

-#include "../../config/config.hpp"
 #ifdef MFEM_USE_CEED
+
+#include "../../../config/config.hpp"
+#include "../../../linalg/vector.hpp"
+#include "../../../linalg/dtensor.hpp"
+#include "../../../mesh/mesh.hpp"
+#include "../../gridfunc.hpp"
 #include "util.hpp"
-#include "../../linalg/vector.hpp"
-#include "../../linalg/dtensor.hpp"
-#include "../../mesh/mesh.hpp"
-#include "../gridfunc.hpp"
-#include "../../linalg/dtensor.hpp"
+#include "ceed.hpp"

 namespace mfem
 {
@@ -9,14 +9,16 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.

-#ifndef MFEM_LIBCEED_PAINTEG
-#define MFEM_LIBCEED_PAINTEG
+#ifndef MFEM_LIBCEED_INTEG
+#define MFEM_LIBCEED_INTEG

+#include "../../../config/config.hpp"
+#include "../../fespace.hpp"
+#include "../../gridfunc.hpp"
 #include "operator.hpp"
-#include "../../config/config.hpp"
 #include "coefficient.hpp"
-#include "../fespace.hpp"
-#include "../gridfunc.hpp"
+#include "restriction.hpp"
+#include "ceed.hpp"

 namespace mfem
 {
@@ -76,8 +78,8 @@ class PAIntegrator : public ceed::Operator
 {
 #ifdef MFEM_USE_CEED
 protected:
-   CeedBasis basis, mesh_basis;
-   CeedElemRestriction restr, mesh_restr, restr_i;
+   CeedBasis  trial_basis, test_basis, mesh_basis;
+   CeedElemRestriction trial_restr, test_restr, mesh_restr, restr_i;
   CeedQFunction build_qfunc, apply_qfunc;
   CeedVector node_coords, qdata;
   Coefficient *coeff;
@@ -85,27 +87,52 @@ protected:
   CeedOperator build_oper;

   PAIntegrator()
-      : Operator(), basis(nullptr), mesh_basis(nullptr),
-        restr(nullptr), mesh_restr(nullptr),
+      : Operator(),
+        trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
+        trial_restr(nullptr), test_restr(nullptr), mesh_restr(nullptr),
        restr_i(nullptr),
        build_qfunc(nullptr), apply_qfunc(nullptr), node_coords(nullptr),
        qdata(nullptr), coeff(nullptr), build_ctx(nullptr), build_oper(nullptr)
   { }

 public:
-   /** This method assembles the PAIntegrator with the given CeedOperatorInfo
-       @a info, an mfem::FiniteElementSpace @a fes, an mfem::IntegrationRule
-       @a ir, and mfem::Coefficient or mfem::VectorCoefficient @a Q.
-       The CeedOperatorInfo type is expected to inherit from OperatorInfo and
-       contain a Context type relevant to the qFunctions. */
+   /** This method assembles the PAIntegrator.
+
+       @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
+                       the `CeedOperatorInfo` type is expected to inherit from
+                       `OperatorInfo` and contain a `Context` type relevant to
+                       the qFunctions.
+       @param[in] fes the `FiniteElementSpace` for the form,
+       @param[in] ir the `IntegrationRule` for the numerical integration,
+       @param[in] Q `Coefficient` or `VectorCoefficient`. */
   template <typename CeedOperatorInfo, typename CoeffType>
   void Assemble(CeedOperatorInfo &info,
                 const mfem::FiniteElementSpace &fes,
                 const mfem::IntegrationRule &irm,
                 CoeffType *Q)
+   {
+      Assemble(info, fes, fes, irm, Q);
+   }
+
+   /** This method assembles the PAIntegrator for mixed forms.
+
+       @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
+                       the `CeedOperatorInfo` type is expected to inherit from
+                       `OperatorInfo` and contain a `Context` type relevant to
+                       the qFunctions.
+       @param[in] trial_fes the trial `FiniteElementSpace` for the form,
+       @param[in] test_fes the test `FiniteElementSpace` for the form,
+       @param[in] ir the `IntegrationRule` for the numerical integration,
+       @param[in] Q `Coefficient` or `VectorCoefficient`. */
+   template <typename CeedOperatorInfo, typename CoeffType>
+   void Assemble(CeedOperatorInfo &info,
+                 const mfem::FiniteElementSpace &trial_fes,
+                 const mfem::FiniteElementSpace &test_fes,
+                 const mfem::IntegrationRule &irm,
+                 CoeffType *Q)
   {
      Ceed ceed(internal::ceed);
-      mfem::Mesh &mesh = *fes.GetMesh();
+      mfem::Mesh &mesh = *trial_fes.GetMesh();
      InitCoefficient(Q, mesh, irm, coeff, info.ctx);
      bool const_coeff = coeff->IsConstant();
      std::string build_func = const_coeff ? info.build_func_const
@@ -119,17 +146,38 @@ public:
                     info.test_op
                    };
      CeedInt nqpts, nelem = mesh.GetNE();
-      CeedInt dim = mesh.SpaceDimension(), vdim = fes.GetVDim();
+      CeedInt dim = mesh.SpaceDimension();
+      CeedInt trial_vdim = trial_fes.GetVDim();
+      CeedInt test_vdim = test_fes.GetVDim();

      mesh.EnsureNodes();
-      InitBasisAndRestriction(fes, irm, ceed, &basis, &restr);
+      if ( &trial_fes == &test_fes )
+      {
+         InitBasisAndRestriction(trial_fes, irm, ceed,
+                                 &trial_basis, &trial_restr);
+         test_basis = trial_basis;
+         test_restr = trial_restr;
+      }
+      else
+      {
+         InitBasisAndRestriction(trial_fes, irm, ceed,
+                                 &trial_basis, &trial_restr);
+         InitBasisAndRestriction(test_fes, irm, ceed,
+                                 &test_basis, &test_restr);
+      }

      const mfem::FiniteElementSpace *mesh_fes = mesh.GetNodalFESpace();
      MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
      InitBasisAndRestriction(*mesh_fes, irm, ceed, &mesh_basis,
                              &mesh_restr);

-      CeedBasisGetNumQuadraturePoints(basis, &nqpts);
+      CeedInt trial_nqpts, test_nqpts;
+      CeedBasisGetNumQuadraturePoints(trial_basis, &trial_nqpts);
+      CeedBasisGetNumQuadraturePoints(test_basis, &test_nqpts);
+      MFEM_VERIFY(trial_nqpts == test_nqpts,
+                  "Trial and test basis must have the same number of quadrature"
+                  " points.");
+      nqpts = trial_nqpts;

      const int qdatasize = op.qdatasize;
      InitStridedRestriction(*mesh_fes, nelem, nqpts, qdatasize,
@@ -143,7 +191,7 @@ public:
      // Context data to be passed to the Q-function.
      info.ctx.dim = mesh.Dimension();
      info.ctx.space_dim = mesh.SpaceDimension();
-      info.ctx.vdim = fes.GetVDim();
+      info.ctx.vdim = trial_fes.GetVDim();

      std::string qf_file = GetCeedPath() + op.header;
      std::string qf = qf_file + op.build_func;
@@ -206,17 +254,22 @@ public:
      switch (op.trial_op)
      {
         case EvalMode::None:
-            CeedQFunctionAddInput(apply_qfunc, "u", vdim, CEED_EVAL_NONE);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
+                                  CEED_EVAL_NONE);
            break;
         case EvalMode::Interp:
-            CeedQFunctionAddInput(apply_qfunc, "u", vdim, CEED_EVAL_INTERP);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
+                                  CEED_EVAL_INTERP);
            break;
         case EvalMode::Grad:
-            CeedQFunctionAddInput(apply_qfunc, "gu", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
+                                  CEED_EVAL_GRAD);
            break;
         case EvalMode::InterpAndGrad:
-            CeedQFunctionAddInput(apply_qfunc, "u", vdim, CEED_EVAL_INTERP);
-            CeedQFunctionAddInput(apply_qfunc, "gu", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
+                                  CEED_EVAL_INTERP);
+            CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
+                                  CEED_EVAL_GRAD);
            break;
      }
      // qdata
@@ -225,17 +278,22 @@ public:
      switch (op.test_op)
      {
         case EvalMode::None:
-            CeedQFunctionAddOutput(apply_qfunc, "v", vdim, CEED_EVAL_NONE);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
+                                   CEED_EVAL_NONE);
            break;
         case EvalMode::Interp:
-            CeedQFunctionAddOutput(apply_qfunc, "v", vdim, CEED_EVAL_INTERP);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
+                                   CEED_EVAL_INTERP);
            break;
         case EvalMode::Grad:
-            CeedQFunctionAddOutput(apply_qfunc, "gv", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
+                                   CEED_EVAL_GRAD);
            break;
         case EvalMode::InterpAndGrad:
-            CeedQFunctionAddOutput(apply_qfunc, "v", vdim, CEED_EVAL_INTERP);
-            CeedQFunctionAddOutput(apply_qfunc, "gv", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
+                                   CEED_EVAL_INTERP);
+            CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
+                                   CEED_EVAL_GRAD);
            break;
      }
      CeedQFunctionSetContext(apply_qfunc, build_ctx);
@@ -246,18 +304,22 @@ public:
      switch (op.trial_op)
      {
         case EvalMode::None:
-            CeedOperatorSetField(oper, "u", restr,
+            CeedOperatorSetField(oper, "u", trial_restr,
                                 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Interp:
-            CeedOperatorSetField(oper, "u", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Grad:
-            CeedOperatorSetField(oper, "gu", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::InterpAndGrad:
-            CeedOperatorSetField(oper, "u", restr, basis, CEED_VECTOR_ACTIVE);
-            CeedOperatorSetField(oper, "gu", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
      }
      // qdata
@@ -267,23 +329,27 @@ public:
      switch (op.test_op)
      {
         case EvalMode::None:
-            CeedOperatorSetField(oper, "v", restr,
+            CeedOperatorSetField(oper, "v", test_restr,
                                 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Interp:
-            CeedOperatorSetField(oper, "v", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "v", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Grad:
-            CeedOperatorSetField(oper, "gv", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gv", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::InterpAndGrad:
-            CeedOperatorSetField(oper, "v", restr, basis, CEED_VECTOR_ACTIVE);
-            CeedOperatorSetField(oper, "gv", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "v", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gv", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
      }

-      CeedVectorCreate(ceed, vdim*fes.GetNDofs(), &u);
-      CeedVectorCreate(ceed, vdim*fes.GetNDofs(), &v);
+      CeedVectorCreate(ceed, trial_vdim*trial_fes.GetNDofs(), &u);
+      CeedVectorCreate(ceed, test_vdim*test_fes.GetNDofs(), &v);
   }

   virtual ~PAIntegrator()
@@ -329,34 +395,59 @@ class MFIntegrator : public ceed::Operator
 {
 #ifdef MFEM_USE_CEED
 protected:
-   CeedBasis basis, mesh_basis;
-   CeedElemRestriction restr, mesh_restr, restr_i;
+   CeedBasis trial_basis, test_basis, mesh_basis;
+   CeedElemRestriction trial_restr, test_restr, mesh_restr, restr_i;
   CeedQFunction apply_qfunc;
   CeedVector node_coords, qdata;
   Coefficient *coeff;
   CeedQFunctionContext build_ctx;

   MFIntegrator()
-      : Operator(), basis(nullptr), mesh_basis(nullptr),
-        restr(nullptr), mesh_restr(nullptr),
+      : Operator(),
+        trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
+        trial_restr(nullptr), test_restr(nullptr), mesh_restr(nullptr),
        restr_i(nullptr),
        apply_qfunc(nullptr), node_coords(nullptr),
        qdata(nullptr), coeff(nullptr), build_ctx(nullptr) { }

 public:
-   /** This method assembles the MFIntegrator with the given CeedOperatorInfo
-       @a info, an mfem::FiniteElementSpace @a fes, an mfem::IntegrationRule
-       @a ir, and mfem::Coefficient or mfem::VectorCoefficient @a Q.
-       The CeedOperatorInfo type is expected to inherit from OperatorInfo and
-       contain a Context type relevant to the qFunctions. */
+   /** This method assembles the MFIntegrator.
+
+       @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
+                       the `CeedOperatorInfo` type is expected to inherit from
+                       `OperatorInfo` and contain a `Context` type relevant to
+                       the qFunctions.
+       @param[in] fes the `FiniteElementSpace` for the form,
+       @param[in] ir the `IntegrationRule` for the numerical integration,
+       @param[in] Q `Coefficient` or `VectorCoefficient`. */
   template <typename CeedOperatorInfo, typename CoeffType>
   void Assemble(CeedOperatorInfo &info,
                 const mfem::FiniteElementSpace &fes,
                 const mfem::IntegrationRule &irm,
                 CoeffType *Q)
+   {
+      Assemble(info, fes, fes, irm, Q);
+   }
+
+   /** This method assembles the MFIntegrator for mixed forms.
+
+       @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
+                       the `CeedOperatorInfo` type is expected to inherit from
+                       `OperatorInfo` and contain a `Context` type relevant to
+                       the qFunctions.
+       @param[in] trial_fes the trial `FiniteElementSpace` for the form,
+       @param[in] test_fes the test `FiniteElementSpace` for the form,
+       @param[in] ir the `IntegrationRule` for the numerical integration,
+       @param[in] Q `Coefficient` or `VectorCoefficient`. */
+   template <typename CeedOperatorInfo, typename CoeffType>
+   void Assemble(CeedOperatorInfo &info,
+                 const mfem::FiniteElementSpace &trial_fes,
+                 const mfem::FiniteElementSpace &test_fes,
+                 const mfem::IntegrationRule &irm,
+                 CoeffType *Q)
   {
      Ceed ceed(internal::ceed);
-      Mesh &mesh = *fes.GetMesh();
+      Mesh &mesh = *trial_fes.GetMesh();
      InitCoefficient(Q, mesh, irm, coeff, info.ctx);
      bool const_coeff = coeff->IsConstant();
      std::string apply_func = const_coeff ? info.apply_func_mf_const
@@ -369,24 +460,45 @@ public:
                     info.test_op
                    };
      CeedInt nqpts, nelem = mesh.GetNE();
-      CeedInt dim = mesh.SpaceDimension(), vdim = fes.GetVDim();
+      CeedInt dim = mesh.SpaceDimension();
+      CeedInt trial_vdim = trial_fes.GetVDim();
+      CeedInt test_vdim = test_fes.GetVDim();

      mesh.EnsureNodes();
-      InitBasisAndRestriction(fes, irm, ceed, &basis, &restr);
+      if ( &trial_fes == &test_fes )
+      {
+         InitBasisAndRestriction(trial_fes, irm, ceed,
+                                 &trial_basis, &trial_restr);
+         test_basis = trial_basis;
+         test_restr = trial_restr;
+      }
+      else
+      {
+         InitBasisAndRestriction(trial_fes, irm, ceed,
+                                 &trial_basis, &trial_restr);
+         InitBasisAndRestriction(test_fes, irm, ceed,
+                                 &test_basis, &test_restr);
+      }

      const mfem::FiniteElementSpace *mesh_fes = mesh.GetNodalFESpace();
      MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
      InitBasisAndRestriction(*mesh_fes, irm, ceed, &mesh_basis,
                              &mesh_restr);

-      CeedBasisGetNumQuadraturePoints(basis, &nqpts);
+      CeedInt trial_nqpts, test_nqpts;
+      CeedBasisGetNumQuadraturePoints(trial_basis, &trial_nqpts);
+      CeedBasisGetNumQuadraturePoints(trial_basis, &test_nqpts);
+      MFEM_VERIFY(trial_nqpts == test_nqpts,
+                  "Trial and test basis must have the same number of quadrature"
+                  " points.");
+      nqpts = trial_nqpts;

      InitVector(*mesh.GetNodes(), node_coords);

      // Context data to be passed to the Q-function.
      info.ctx.dim = mesh.Dimension();
      info.ctx.space_dim = mesh.SpaceDimension();
-      info.ctx.vdim = fes.GetVDim();
+      info.ctx.vdim = trial_fes.GetVDim();

      std::string qf_file = GetCeedPath() + op.header;
      std::string qf = qf_file + op.apply_func;
@@ -405,17 +517,22 @@ public:
      switch (op.trial_op)
      {
         case EvalMode::None:
-            CeedQFunctionAddInput(apply_qfunc, "u", vdim, CEED_EVAL_NONE);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
+                                  CEED_EVAL_NONE);
            break;
         case EvalMode::Interp:
-            CeedQFunctionAddInput(apply_qfunc, "u", vdim, CEED_EVAL_INTERP);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
+                                  CEED_EVAL_INTERP);
            break;
         case EvalMode::Grad:
-            CeedQFunctionAddInput(apply_qfunc, "gu", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
+                                  CEED_EVAL_GRAD);
            break;
         case EvalMode::InterpAndGrad:
-            CeedQFunctionAddInput(apply_qfunc, "u", vdim, CEED_EVAL_INTERP);
-            CeedQFunctionAddInput(apply_qfunc, "gu", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
+                                  CEED_EVAL_INTERP);
+            CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
+                                  CEED_EVAL_GRAD);
            break;
      }
      CeedQFunctionAddInput(apply_qfunc, "dx", dim * dim, CEED_EVAL_GRAD);
@@ -424,17 +541,22 @@ public:
      switch (op.test_op)
      {
         case EvalMode::None:
-            CeedQFunctionAddOutput(apply_qfunc, "v", vdim, CEED_EVAL_NONE);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
+                                   CEED_EVAL_NONE);
            break;
         case EvalMode::Interp:
-            CeedQFunctionAddOutput(apply_qfunc, "v", vdim, CEED_EVAL_INTERP);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
+                                   CEED_EVAL_INTERP);
            break;
         case EvalMode::Grad:
-            CeedQFunctionAddOutput(apply_qfunc, "gv", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
+                                   CEED_EVAL_GRAD);
            break;
         case EvalMode::InterpAndGrad:
-            CeedQFunctionAddOutput(apply_qfunc, "v", vdim, CEED_EVAL_INTERP);
-            CeedQFunctionAddOutput(apply_qfunc, "gv", vdim*dim, CEED_EVAL_GRAD);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
+                                   CEED_EVAL_INTERP);
+            CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
+                                   CEED_EVAL_GRAD);
            break;
      }

@@ -470,18 +592,22 @@ public:
      switch (op.trial_op)
      {
         case EvalMode::None:
-            CeedOperatorSetField(oper, "u", restr,
+            CeedOperatorSetField(oper, "u", trial_restr,
                                 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Interp:
-            CeedOperatorSetField(oper, "u", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Grad:
-            CeedOperatorSetField(oper, "gu", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::InterpAndGrad:
-            CeedOperatorSetField(oper, "u", restr, basis, CEED_VECTOR_ACTIVE);
-            CeedOperatorSetField(oper, "gu", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
      }
      CeedOperatorSetField(oper, "dx", mesh_restr,
@@ -492,23 +618,27 @@ public:
      switch (op.test_op)
      {
         case EvalMode::None:
-            CeedOperatorSetField(oper, "v", restr,
+            CeedOperatorSetField(oper, "v", test_restr,
                                 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Interp:
-            CeedOperatorSetField(oper, "v", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "v", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::Grad:
-            CeedOperatorSetField(oper, "gv", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gv", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
         case EvalMode::InterpAndGrad:
-            CeedOperatorSetField(oper, "v", restr, basis, CEED_VECTOR_ACTIVE);
-            CeedOperatorSetField(oper, "gv", restr, basis, CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "v", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gv", test_restr, test_basis,
+                                 CEED_VECTOR_ACTIVE);
            break;
      }

-      CeedVectorCreate(ceed, vdim*fes.GetNDofs(), &u);
-      CeedVectorCreate(ceed, vdim*fes.GetNDofs(), &v);
+      CeedVectorCreate(ceed, trial_vdim*trial_fes.GetNDofs(), &u);
+      CeedVectorCreate(ceed, test_vdim*test_fes.GetNDofs(), &v);
   }

   virtual ~MFIntegrator()
@@ -545,4 +675,4 @@ private:

 } // namespace mfem

-#endif // MFEM_LIBCEED_PAINTEG
+#endif // MFEM_LIBCEED_INTEG
@@ -0,0 +1,30 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_LIBCEED_INTERFACE
+#define MFEM_LIBCEED_INTERFACE
+
+// Object wrapping a CeedOperator in a mfem::Operator.
+#include "operator.hpp"
+// Functions to initialize CeedBasis objects.
+#include "basis.hpp"
+// Functions to initialize CeedRestriction objects.
+#include "restriction.hpp"
+// Functions to initialize coefficients.
+#include "coefficient.hpp"
+// PA or MF Operator using libCEED.
+#include "integrator.hpp"
+// Utility functions
+#include "util.hpp"
+// Wrapper to include <ceed.h>
+#include "ceed.hpp"
+
+#endif // MFEM_LIBCEED_INTERFACE
@@ -11,9 +11,11 @@

 #include "operator.hpp"

-#include "../../config/config.hpp"
-#include "../../linalg/vector.hpp"
-#include "../fespace.hpp"
+#include "../../../config/config.hpp"
+#include "../../../linalg/vector.hpp"
+#include "../../fespace.hpp"
+#include "util.hpp"
+#include "ceed.hpp"

 namespace mfem
 {
@@ -12,8 +12,8 @@
 #ifndef MFEM_LIBCEED_OPERATOR
 #define MFEM_LIBCEED_OPERATOR

-#include "util.hpp"
-#include "../../linalg/operator.hpp"
+#include "../../../linalg/operator.hpp"
+#include "ceed.hpp"

 namespace mfem
 {
@@ -41,7 +41,6 @@ public:
   void AddMult(const mfem::Vector &x, mfem::Vector &y) const;
   void GetDiagonal(mfem::Vector &diag) const;
   using mfem::Operator::SetupRAP;
-
   virtual ~Operator()
   {
 #ifdef MFEM_USE_CEED
@@ -0,0 +1,162 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../../../fem/gridfunc.hpp"
+#include "ceed.hpp"
+
+namespace mfem
+{
+
+namespace ceed
+{
+
+#ifdef MFEM_USE_CEED
+
+static void InitNonTensorRestriction(const mfem::FiniteElementSpace &fes,
+                                     Ceed ceed, CeedElemRestriction *restr)
+{
+   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const int P = fe->GetDof();
+   CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
+   const mfem::Table &el_dof = fes.GetElementToDofTable();
+   mfem::Array<int> tp_el_dof(el_dof.Size_of_connections());
+   const mfem::TensorBasisElement * tfe =
+      dynamic_cast<const mfem::TensorBasisElement *>(fe);
+   const int stride = compstride == 1 ? fes.GetVDim() : 1;
+   if (tfe) // Lexicographic ordering using dof_map
+   {
+      const mfem::Array<int>& dof_map = tfe->GetDofMap();
+      for (int i = 0; i < fes.GetNE(); i++)
+      {
+         const int el_offset = P * i;
+         for (int j = 0; j < P; j++)
+         {
+            tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[dof_map[j]+el_offset];
+         }
+      }
+   }
+   else  // Native ordering
+   {
+      for (int e = 0; e < fes.GetNE(); e++)
+      {
+         for (int i = 0; i < P; i++)
+         {
+            tp_el_dof[i + e*P] = stride*el_dof.GetJ()[i + e*P];
+         }
+      }
+   }
+   CeedElemRestrictionCreate(ceed, fes.GetNE(), P, fes.GetVDim(),
+                             compstride, (fes.GetVDim())*(fes.GetNDofs()),
+                             CEED_MEM_HOST, CEED_COPY_VALUES,
+                             tp_el_dof.GetData(), restr);
+}
+
+// TODO fuse Tensor and NonTensor Restriction
+void InitTensorRestriction(const mfem::FiniteElementSpace &fes,
+                           Ceed ceed, CeedElemRestriction *restr)
+{
+   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const mfem::TensorBasisElement * tfe =
+      dynamic_cast<const mfem::TensorBasisElement *>(fe);
+   MFEM_VERIFY(tfe, "invalid FE");
+   const mfem::Array<int>& dof_map = tfe->GetDofMap();
+
+   CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
+   const mfem::Table &el_dof = fes.GetElementToDofTable();
+   mfem::Array<int> tp_el_dof(el_dof.Size_of_connections());
+   const int dof = fe->GetDof();
+   const int stride = compstride == 1 ? fes.GetVDim() : 1;
+   if (dof_map.Size()>0)
+   {
+      for (int i = 0; i < fes.GetNE(); i++)
+      {
+         const int el_offset = dof * i;
+         for (int j = 0; j < dof; j++)
+         {
+            tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[dof_map[j]+el_offset];
+         }
+      }
+   }
+   else // dof_map.Size == 0, means dof_map[j]==j;
+   {
+      for (int i = 0; i < fes.GetNE(); i++)
+      {
+         const int el_offset = dof * i;
+         for (int j = 0; j < dof; j++)
+         {
+            tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[j+el_offset];
+         }
+      }
+   }
+   CeedElemRestrictionCreate(ceed, fes.GetNE(), dof, fes.GetVDim(),
+                             compstride, (fes.GetVDim())*(fes.GetNDofs()),
+                             CEED_MEM_HOST, CEED_COPY_VALUES,
+                             tp_el_dof.GetData(), restr);
+}
+
+void InitStridedRestriction(const mfem::FiniteElementSpace &fes,
+                            CeedInt nelem, CeedInt nqpts, CeedInt qdatasize,
+                            const CeedInt *strides,
+                            CeedElemRestriction *restr)
+{
+   RestrKey restr_key(&fes, nelem, nqpts, qdatasize, restr_type::Strided);
+   auto restr_itr = mfem::internal::ceed_restr_map.find(restr_key);
+   if (restr_itr == mfem::internal::ceed_restr_map.end())
+   {
+      CeedElemRestrictionCreateStrided(mfem::internal::ceed, nelem, nqpts, qdatasize,
+                                       nelem*nqpts*qdatasize,
+                                       strides,
+                                       restr);
+      // Will be automatically destroyed when @a fes gets destroyed.
+      mfem::internal::ceed_restr_map[restr_key] = *restr;
+   }
+   else
+   {
+      *restr = restr_itr->second;
+   }
+}
+
+void InitRestriction(const FiniteElementSpace &fes,
+                     Ceed ceed,
+                     CeedElemRestriction *restr)
+{
+   // Check for FES -> basis, restriction in hash tables
+   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const int P = fe->GetDof();
+   const int nelem = fes.GetNE();
+   const int ncomp = fes.GetVDim();
+   RestrKey restr_key(&fes, nelem, P, ncomp, restr_type::Standard);
+   auto restr_itr = mfem::internal::ceed_restr_map.find(restr_key);
+
+   // Init or retreive key values
+   if (restr_itr == mfem::internal::ceed_restr_map.end())
+   {
+      if (UsesTensorBasis(fes))
+      {
+         InitTensorRestriction(fes, ceed, restr);
+      }
+      else
+      {
+         InitNonTensorRestriction(fes, ceed, restr);
+      }
+      mfem::internal::ceed_restr_map[restr_key] = *restr;
+   }
+   else
+   {
+      *restr = restr_itr->second;
+   }
+}
+
+#endif
+
+} // namespace ceed
+
+} // namespace mfem
@@ -0,0 +1,62 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_LIBCEED_RESTR
+#define MFEM_LIBCEED_RESTR
+
+#include "ceed.hpp"
+
+namespace mfem
+{
+
+namespace ceed
+{
+
+#ifdef MFEM_USE_CEED
+/// @brief Initialize a strided CeedElemRestriction
+/** @a nelem is the number of elements,
+    @a nqpts is the total number of quadrature points
+    @a qdatasize is the number of data per quadrature point
+    @a strides Array for strides between [nodes, components, elements].
+    Data for node i, component j, element k can be found in the L-vector at
+    index i*strides[0] + j*strides[1] + k*strides[2]. CEED_STRIDES_BACKEND may
+    be used with vectors created by a Ceed backend. */
+void InitStridedRestriction(const mfem::FiniteElementSpace &fes,
+                            CeedInt nelem, CeedInt nqpts, CeedInt qdatasize,
+                            const CeedInt *strides,
+                            CeedElemRestriction *restr);
+
+/** @brief Initialize a CeedElemRestriction.
+ *
+ * @param[in] fes Input finite element space.
+ * @param[in] ceed Input Ceed object.
+   @param[out] restr The address of the initialized CeedElemRestriction object.
+*/
+void InitRestriction(const FiniteElementSpace &fes,
+                     Ceed ceed,
+                     CeedElemRestriction *restr);
+
+/** @brief Initialize a CeedElemRestriction.
+ *
+ * @param[in] fes Input finite element space.
+ * @param[in] ceed Input Ceed object.
+   @param[out] restr The address of the initialized CeedElemRestriction object.
+*/
+void InitTensorRestriction(const FiniteElementSpace &fes,
+                           Ceed ceed, CeedElemRestriction *restr);
+
+#endif
+
+} // namespace ceed
+
+} // namespace mfem
+
+#endif // MFEM_LIBCEED_RESTR
--- a/Show More
+++ b/Show More