Try the 'lgtm' CodeQL queries

Windows CI METIS workaround

.github/workflows/build-container.yaml

@@ -1,4 +1,4 @@
-name: Build Deploy Container
+name: "Docker"
 
 on:
 

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

@@ -10,7 +10,7 @@
 # CONTRIBUTING.md for details.
 
 # In this CI section, we build different variants of mfem and run test on them.
-name: builds-and-tests
+name: "Tests"
 
 # Github actions can use the default "GITHUB_TOKEN". By default, this token
 # is set to have permissive access. However, this is not a good practice
@@ -47,17 +47,17 @@ jobs:
   builds-and-tests:
     strategy:
       matrix:
-        os: [ubuntu-20.04, macos-10.15, windows-2022]
+        os: [ubuntu-latest, macos-latest, windows-latest]
         target: [dbg, opt]
         mpi: [seq, par]
         build-system: [make, cmake]
         hypre-target: [int32]
         exclude:
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             build-system: cmake
-          - os: macos-10.15
+          - os: macos-latest
             build-system: cmake
-          - os: windows-2022
+          - os: windows-latest
             build-system: make
         # 'include' allows us to:
         # - Add a variable to all jobs without creating a new matrix dimension.
@@ -72,15 +72,15 @@ jobs:
             codecov: NO
           - target: opt
             codecov: YES
-          - os: windows-2022
+          - os: windows-latest
             codecov: NO
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             target: opt
             codecov: NO
             mpi: par
             build-system: cmake
             hypre-target: int32
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             target: opt
             codecov: NO
             mpi: par
@@ -112,35 +112,35 @@ 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-20.04'
+      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-latest'
       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-20.04'
+      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-latest'
       run: |
         sudo apt-get install lcov
 
     - name: Set up Homebrew
-      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-10.15'
+      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
       uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3
 
     - name: get MPI (MacOS)
-      if: matrix.mpi == 'par' && matrix.os == 'macos-10.15'
+      if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install openmpi
         export MAKE_CXX_FLAG="MPICXX=mpic++"
 
     - name: get MPI (MacOS)
-      if: matrix.codecov == 'YES' && matrix.os == 'macos-10.15'
+      if: matrix.codecov == 'YES' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install lcov
 
     - name: get MPI (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
       uses: mpi4py/setup-mpi@v1.0.3
 
     # Get Hypre through cache, or build it.
@@ -154,7 +154,7 @@ jobs:
         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' && matrix.os != 'windows-2022'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-latest'
       uses: mfem/github-actions/build-hypre@v2.2
       with:
         archive: ${{ env.HYPRE_ARCHIVE }}
@@ -163,7 +163,7 @@ jobs:
         build-system: make
 
     - name: get hypre (Windows)
-      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-latest'
       uses: mfem/github-actions/build-hypre@v2.2
       with:
         archive: ${{ env.HYPRE_ARCHIVE }}
@@ -175,14 +175,14 @@ jobs:
     # Install will only run on cache miss.
     - name: cache metis
       id: metis-cache
-      if: matrix.mpi == 'par' && matrix.os != 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
       uses: actions/cache@v2
       with:
         path: ${{ env.METIS_TOP_DIR }}
         key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
 
     - name: install metis
-      if: matrix.mpi == 'par' && matrix.os != 'windows-2022' && steps.metis-cache.outputs.cache-hit != 'true'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-latest' && steps.metis-cache.outputs.cache-hit != 'true'
       uses: mfem/github-actions/build-metis@v2.2
       with:
         archive: ${{ env.METIS_ARCHIVE }}
@@ -196,20 +196,16 @@ jobs:
         key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
 
     - name: prepare binary cache location
-      if: matrix.os == 'windows-2022' && steps.vcpkg-cache.outputs.cache-hit != 'true'
+      if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
       run: |
         mkdir -p vcpkg_cache
 
     - name: install metis (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
       env:
         VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
       run: |
-        $PortFile = 'C:\vcpkg\ports\metis\portfile.cmake'
-        $OriginalURL = 'http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-${METIS_VERSION}.tar.gz'
-        $NewURL = 'https://github.com/mfem/tpls/raw/gh-pages/metis-5.1.0.tar.gz'
-        (Get-Content $PortFile).replace($OriginalURL, $NewURL) | Set-Content $PortFile
-        vcpkg install metis --triplet=x64-windows-static
+        vcpkg install metis-mfem --triplet=x64-windows-static --overlay-ports=${{ env.MFEM_TOP_DIR }}/config/vcpkg/ports
 
     # MFEM build and test
     - name: build
@@ -252,7 +248,7 @@ jobs:
       shell: bash
 
     - name: cmake unit tests (Ubuntu 20.04)
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-20.04'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
       run: |
         CTEST_CONFIG="Release"
         [[ ${{ matrix.target }} == 'dbg' ]] && CTEST_CONFIG="Debug"
@@ -260,7 +256,7 @@ jobs:
       shell: bash
 
     - name: cmake tests
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-20.04'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-latest'
       run: |
         CTEST_CONFIG="Release"
         cd ${{ env.MFEM_TOP_DIR }}/build && ctest --output-on-failure -C ${CTEST_CONFIG}

.github/workflows/codeql-analysis.yml

@@ -0,0 +1,71 @@
+# For most projects, this workflow file will not need changing; you simply need
+# to commit it to your repository.
+#
+# You may wish to alter this file to override the set of languages analyzed,
+# or to provide custom queries or build logic.
+#
+# ******** NOTE ********
+# We have attempted to detect the languages in your repository. Please check
+# the `language` matrix defined below to confirm you have the correct set of
+# supported CodeQL languages.
+#
+name: "Static Analysis"
+
+on:
+  push:
+    branches: [ "master", "next"]
+  pull_request:
+    # The branches below must be a subset of the branches above
+    branches: [ "master" ]
+
+jobs:
+  analyze:
+    name: Analyze
+    runs-on: ubuntu-latest
+    permissions:
+      actions: read
+      contents: read
+      security-events: write
+
+    strategy:
+      fail-fast: false
+      matrix:
+        language: [ 'cpp' ]
+        # CodeQL supports [ 'cpp', 'csharp', 'go', 'java', 'javascript', 'python', 'ruby' ]
+        # Learn more about CodeQL language support at https://aka.ms/codeql-docs/language-support
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v3
+
+    # Initializes the CodeQL tools for scanning.
+    - name: Initialize CodeQL
+      uses: github/codeql-action/init@v2
+      with:
+        languages: ${{ matrix.language }}
+        # If you wish to specify custom queries, you can do so here or in a config file.
+        # By default, queries listed here will override any specified in a config file.
+        # Prefix the list here with "+" to use these queries and those in the config file.
+
+        # Details on CodeQL's query packs refer to : https://docs.github.com/en/code-security/code-scanning/automatically-scanning-your-code-for-vulnerabilities-and-errors/configuring-code-scanning#using-queries-in-ql-packs
+        # queries: security-extended,security-and-quality
+        queries: lgtm
+
+
+    # Autobuild attempts to build any compiled languages  (C/C++, C#, or Java).
+    # If this step fails, then you should remove it and run the build manually (see below)
+    - name: Autobuild
+      uses: github/codeql-action/autobuild@v2
+
+    # ℹ️ Command-line programs to run using the OS shell.
+    # 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun
+
+    #   If the Autobuild fails above, remove it and uncomment the following three lines.
+    #   modify them (or add more) to build your code if your project, please refer to the EXAMPLE below for guidance.
+
+    # - run: |
+    #   echo "Run, Build Application using script"
+    #   ./location_of_script_within_repo/buildscript.sh
+
+    - name: Perform CodeQL Analysis
+      uses: github/codeql-action/analyze@v2

.github/workflows/mfem-analysis.yml

@@ -9,7 +9,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-name: build-analysis
+name: "Build Analysis"
 
 permissions:
   actions: write
@@ -31,7 +31,7 @@ env:
 
 jobs:
   gitignore:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
 
     steps:
     - name: Cancel Previous Runs

.github/workflows/repo-check.yml

@@ -9,7 +9,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-name: repo-check
+name: "Checks"
 
 permissions:
   actions: write
@@ -28,7 +28,7 @@ on:
 
 jobs:
   file-headers-check:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -78,7 +78,7 @@ jobs:
         exit 1
 
   code-style:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -88,14 +88,14 @@ jobs:
 
     - name: get astyle
       run: |
-        sudo apt-get install astyle=3.1-1ubuntu2
+        sudo apt-get install astyle
 
     - name: style check
       run: |
         ./config/githooks/pre-push --style
 
   documentation:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -106,6 +106,8 @@ jobs:
     - name: get doxygen and graphviz
       run: |
         sudo apt-get install doxygen graphviz
+        cd doc
+        doxygen -u CodeDocumentation.conf.in 2>/dev/null
 
     - name: build documentation
       run: |
@@ -118,7 +120,7 @@ jobs:
       github.ref != 'refs/heads/master' &&
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     steps:
     - name: checkout mfem
       uses: actions/checkout@v2

CHANGELOG

@@ -11,6 +11,38 @@
 Version 4.4.1 (development)
 ===========================
 
+Meshing improvements
+--------------------
+- Added support for mixed meshes and pyramids in GSLIB-FindPoints.
+
+Discretization improvements
+---------------------------
+- 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 support for partial assembly and fully matrix-free operators on mixed
+  meshes (different element types and p-adaptivity) through libCEED, including
+  device acceleration, e.g. with NVIDIA and AMD GPUs. The p-adaptivity is
+  currently limited by MFEM capabilities, i.e. 2D serial meshes. All mixed
+  element topologies are supported in serial and parallel: segment, triangle,
+  square, tetrahedron, cube, prism, and pyramid.
+
+- 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))
+
+- 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.
+
+Linear and nonlinear solvers
+----------------------------
+
 New and updated examples and miniapps
 -------------------------------------
 - Added a new elasticity miniapp, Hooke, that showcases a low-level approach of
@@ -19,14 +51,12 @@ New and updated examples and miniapps
   automatic differentiation tools like a native dual number implementation or a
   third party library such as Enzyme. See miniapps/elasticity for more details.
 
-- Added example for body-fitted volumetric and shape integration using the 
-  Algoim library.
+- Add a new example code, Example 33/33p, to demonstrate the solution of
+  spectral fractional PDEs with MFEM.
 
-- 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.
+Integrations, testing and documentation
+---------------------------------------
+- Added a Dockerfile for a simple MFEM container, see config/docker/README.md.
 
 - Added support for ParMoonolith, https://bitbucket.org/zulianp/par_moonolith,
   which provides parallel non-conforming, non-matching, variational, volumetric
@@ -34,24 +64,20 @@ New and updated examples and miniapps
   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))
+- Added support for the LLVM-based automatic differentiation tool Enzyme, see
+  https://github.com/EnzymeAD/Enzyme. Build system flags and a convenience
+  header are provided. The functionality and interaction are demonstrated in a
+  new miniapp in miniapps/elasticity.
 
-- 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 example for body-fitted volumetric and shape integration using the
+  Algoim library.
 
 - Added Windows 2022 CI testing with GitHub actions.
 
-- Added support for mixed meshes and pyramids in GSLIB-FindPoints.
+Miscellaneous
+-------------
+- Various other simplifications, extensions, and bugfixes in the code.
+
 
 - Added boundary elimination with device support for `SparseMatrix` and
   `HypreParMatrix`.
@@ -65,14 +91,6 @@ New and updated examples and miniapps
 Version 4.4, released on March 21, 2022
 =======================================
 
-Integrations, testing and documentation
----------------------------------------
-- Added support for the LLVM-based automatic differentiation tool Enzyme, see
-  https://github.com/EnzymeAD/Enzyme. Build system flags and a convenience
-  header are provided. The functionality and interaction are demonstrated in a
-  new miniapp in miniapps/elasticity.
-
-
 Linear and nonlinear solvers
 ----------------------------
 - Added support for using the hypre library built with HIP support. Similar to
@@ -101,6 +119,11 @@ Meshing improvements
 - Added a simpler interface to access mesh face information, see FaceInformation
   and GetFaceInformation in the Mesh class.
 
+- Added the method ParMesh::GetSerialMesh() that reconstructs a partitioned
+  parallel mesh on a given single rank. Also, added the method
+  ParMesh::PrintAsSerial() that saves the reconstructed serial mesh to a C++
+  stream on rank 0.
+
 - Gmsh meshes where all elements have zero physical tag (the default Gmsh output
   format if no physical groups are defined) are now successfully loaded, and
   elements are reassigned attribute number 1.
@@ -203,6 +226,13 @@ Miscellaneous
 
 - Fixed several MinGW build issues on Windows.
 
+- In various places in the library, replace the use of 'long' with 'long long'
+  to better support Win64 builds where 'long' is 32-bit and 'long long' is
+  64-bit. On Linux and MacOS, both types are typically 64-bit.
+
+- Update various "MemoryUsage" methods to return 'std::size_t' instead of 'long'
+  since the latter is 32-bit in Win64 builds.
+
 - Added 'double' atomicAdd implementation for previous versions of CUDA.
 
 - HypreParVector and Vector now support C++ move semantics, and the copy

CMakeLists.txt

@@ -136,6 +136,8 @@ if (MFEM_USE_CUDA)
        "CUDA flags set for MFEM" FORCE)
    set(CUSPARSE_FOUND TRUE)
    set(CUSPARSE_LIBRARIES "cusparse")
+   set(CUBLAS_FOUND TRUE)
+   set(CUSBLAS_LIBRARIES "cublas")
 endif()
 
 if (XSDK_ENABLE_C)
@@ -483,7 +485,7 @@ endif()
 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
+    ADIOS2 CUBLAS CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
     MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)
 
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.

config/defaults.mk

@@ -207,7 +207,7 @@ HYPRE_OPT = -I$(HYPRE_DIR)/include
 HYPRE_LIB = -L$(HYPRE_DIR)/lib -lHYPRE
 ifeq (YES,$(MFEM_USE_CUDA))
    # This is only necessary when hypre is built with cuda:
-   HYPRE_LIB += -lcusparse -lcurand
+   HYPRE_LIB += -lcusparse -lcurand -lcublas
 endif
 ifeq (YES,$(MFEM_USE_HIP))
    # This is only necessary when hypre is built with hip:

config/vcpkg/ports/metis-mfem/disable-programs.patch

@@ -0,0 +1,8 @@
+--- a/CMakeLists.txt	Wed Dec 21 18:24:22 2016
++++ b/CMakeLists.txt	Wed Dec 21 18:24:26 2016
+@@ -20,4 +20,4 @@
+ # Recursively look for CMakeLists.txt in subdirs.
+ add_subdirectory("include")
+ add_subdirectory("libmetis")
+-add_subdirectory("programs")
++# add_subdirectory("programs")

config/vcpkg/ports/metis-mfem/enable-install.patch

@@ -0,0 +1,15 @@
+--- a/CMakeLists.txt	Sat Mar 30 17:24:45 2013
++++ b/CMakeLists.txt	Wed Dec 21 18:23:43 2016
+@@ -4,11 +4,7 @@
+ set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
+ set(SHARED FALSE CACHE BOOL "build a shared library")
+ 
+-if(MSVC)
+-  set(METIS_INSTALL FALSE)
+-else()
+-  set(METIS_INSTALL TRUE)
+-endif()
++set(METIS_INSTALL TRUE)
+ 
+ # Configure libmetis library.
+ if(SHARED)

config/vcpkg/ports/metis-mfem/fix-INT_MIN_define.patch

@@ -0,0 +1,34 @@
+diff --git a/include/metis.h b/include/metis.h
+index dc5406a..7732437 100644
+--- a/include/metis.h
++++ b/include/metis.h
+@@ -72,10 +72,14 @@ typedef __int64 int64_t;
+ #define PRId64       "I64d"
+ #define SCNd32       "ld"
+ #define SCNd64       "I64d"
++#ifdef _WIN32
++#include <stdint.h>
++#else
+ #define INT32_MIN    ((int32_t)_I32_MIN)
+ #define INT32_MAX    _I32_MAX
+ #define INT64_MIN    ((int64_t)_I64_MIN)
+ #define INT64_MAX    _I64_MAX
++#endif
+ #else
+ #include <inttypes.h>
+ #endif
+diff --git a/GKlib/gk_arch.h b/GKlib/gk_arch.h
+index 78b1431..7258763 100644
+--- a/GKlib/gk_arch.h
++++ b/GKlib/gk_arch.h
+@@ -32,8 +32,8 @@
+ 
+ 
+ #ifdef __MSC__ 
+-  #include "ms_stdint.h"
+-  #include "ms_inttypes.h"
++  #include <stdint.h>
++  #include <inttypes.h>
+   #include "ms_stat.h"
+ #else
+ #ifndef SUNOS

config/vcpkg/ports/metis-mfem/fix-gklib-vs14-math.patch

@@ -0,0 +1,11 @@
+--- a/GKlib/gk_arch.h	Wed Dec 21 18:34:18 2016
++++ b/GKlib/gk_arch.h	Wed Dec 21 18:30:49 2016
+@@ -58,7 +58,7 @@
+ #define PTRDIFF_MAX  INT64_MAX
+ #endif
+ 
+-#ifdef __MSC__
++#if defined(__MSC__) && (_MSC_VER < 1900)
+ /* MSC does not have rint() function */
+ #define rint(x) ((int)((x)+0.5))  
+ 

config/vcpkg/ports/metis-mfem/fix-linux-build-error.patch

@@ -0,0 +1,14 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index e94f050..b9613a7 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -1,7 +1,8 @@
+ cmake_minimum_required(VERSION 2.8)
+ project(METIS)
+ 
+-set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
++set(GKLIB_PATH "${CMAKE_SOURCE_DIR}/GKlib" CACHE PATH "path to GKlib")
++
+ set(SHARED FALSE CACHE BOOL "build a shared library")
+ 
+ set(METIS_INSTALL TRUE)

config/vcpkg/ports/metis-mfem/fix-metis-vs14-math.patch

@@ -0,0 +1,11 @@
+--- a/libmetis/metislib.h	Sat Mar 30 17:24:45 2013
++++ b/libmetis/metislib.h	Wed Dec 21 18:30:59 2016
+@@ -31,7 +31,7 @@
+ #include <proto.h>
+ 
+ 
+-#if defined(COMPILER_MSC)
++#if defined(COMPILER_MSC) && (_MSC_VER < 1900)
+ #if defined(rint)
+   #undef rint
+ #endif

config/vcpkg/ports/metis-mfem/fix-runtime-install-destination.patch

@@ -0,0 +1,10 @@
+--- a/libmetis/CMakeLists.txt	Sat Mar 30 17:24:45 2013
++++ b/libmetis/CMakeLists.txt	Wed Dec 21 17:41:37 2016
+@@ -11,6 +11,6 @@
+ if(METIS_INSTALL)
+   install(TARGETS metis
+     LIBRARY DESTINATION lib
+-    RUNTIME DESTINATION lib
++    RUNTIME DESTINATION bin
+     ARCHIVE DESTINATION lib)
+ endif()

config/vcpkg/ports/metis-mfem/install-metisConfig.patch

@@ -0,0 +1,44 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index b9613a7..e43ffee 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -22,3 +22,23 @@ include_directories(include)
+ add_subdirectory("include")
+ add_subdirectory("libmetis")
+ # add_subdirectory("programs")
++
++if(METIS_INSTALL)
++  set(PRJ_NAME metis)
++  set(PRJ_VER 5.1.0)
++  install(EXPORT metisTargets 
++    FILE ${PRJ_NAME}Targets.cmake
++    DESTINATION lib/cmake/${PRJ_NAME})
++  include(CMakePackageConfigHelpers)
++  write_basic_package_version_file(
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
++    VERSION ${PRJ_VER}
++    COMPATIBILITY SameMajorVersion)
++  file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
++    "include(\${CMAKE_CURRENT_LIST_DIR}/${PRJ_NAME}Targets.cmake)")
++  install(FILES 
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
++    DESTINATION lib/cmake/${PRJ_NAME})
++endif()
++
+diff --git a/libmetis/CMakeLists.txt b/libmetis/CMakeLists.txt
+index 7a5fc74..5a68cf0 100644
+--- a/libmetis/CMakeLists.txt
++++ b/libmetis/CMakeLists.txt
+@@ -9,8 +9,9 @@ if(UNIX)
+ endif()
+ 
+ if(METIS_INSTALL)
+-  install(TARGETS metis
++  install(TARGETS metis EXPORT metisTargets
+     LIBRARY DESTINATION lib
+     RUNTIME DESTINATION bin
+-    ARCHIVE DESTINATION lib)
++    ARCHIVE DESTINATION lib
++    INCLUDES DESTINATION include)
+ endif()

config/vcpkg/ports/metis-mfem/portfile.cmake

@@ -0,0 +1,41 @@
+vcpkg_check_linkage(ONLY_STATIC_LIBRARY)
+set(OPTIONS -DSHARED=OFF)
+
+set(METIS_VERSION 5.1.0)
+
+vcpkg_download_distfile(ARCHIVE
+    URLS "https://github.com/mfem/tpls/raw/gh-pages/metis-${METIS_VERSION}.tar.gz"
+    FILENAME "metis-${METIS_VERSION}.tar.gz"
+    SHA512 deea47749d13bd06fbeaf98a53c6c0b61603ddc17a43dae81d72c8015576f6495fd83c11b0ef68d024879ed5415c14ebdbd87ce49c181bdac680573bea8bdb25
+)
+
+vcpkg_extract_source_archive_ex(
+    OUT_SOURCE_PATH SOURCE_PATH
+    ARCHIVE ${ARCHIVE}
+    REF ${METIS_VERSION}
+    PATCHES
+        enable-install.patch
+        disable-programs.patch
+        fix-runtime-install-destination.patch
+        fix-metis-vs14-math.patch
+        fix-gklib-vs14-math.patch
+        fix-linux-build-error.patch
+        install-metisConfig.patch
+        fix-INT_MIN_define.patch
+)
+
+vcpkg_configure_cmake(
+    SOURCE_PATH ${SOURCE_PATH}
+    PREFER_NINJA
+    OPTIONS ${OPTIONS}
+)
+
+vcpkg_install_cmake()
+vcpkg_copy_pdbs()
+vcpkg_fixup_cmake_targets(CONFIG_PATH lib/cmake/metis)
+
+file(REMOVE_RECURSE ${CURRENT_PACKAGES_DIR}/debug/include)
+
+# Handle copyright
+file(COPY ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/metis)
+file(INSTALL ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/${PORT} RENAME copyright)

config/vcpkg/ports/metis-mfem/vcpkg.json

@@ -0,0 +1,7 @@
+{
+  "name": "metis-mfem",
+  "version-string": "5.1.0",
+  "port-version": 0,
+  "description": "Serial Graph Partitioning and Fill-reducing Matrix Ordering",
+  "homepage": "https://glaros.dtc.umn.edu/gkhome/metis/metis/overview"
+}

doc/CodeDocumentation.conf.in

@@ -2349,7 +2349,7 @@ PLANTUML_INCLUDE_PATH  =
 # Minimum value: 0, maximum value: 10000, default value: 50.
 # This tag requires that the tag HAVE_DOT is set to YES.
 
-DOT_GRAPH_MAX_NODES    = 50
+DOT_GRAPH_MAX_NODES    = 100
 
 # The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the graphs
 # generated by dot. A depth value of 3 means that only nodes reachable from the

examples/ex1.cpp

@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               ex1 -pa -d cuda
-//               ex1 -fa -d cuda
+//             * ex1 -fa -d cuda
 //               ex1 -pa -d raja-cuda
 //             * ex1 -pa -d raja-hip
 //               ex1 -pa -d occa-cuda
@@ -38,9 +38,13 @@
 //               ex1 -pa -d occa-omp
 //               ex1 -pa -d ceed-cpu
 //               ex1 -pa -d ceed-cpu -o 4 -a
+//               ex1 -pa -d ceed-cpu -m ../data/square-mixed.mesh
+//               ex1 -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
 //             * ex1 -pa -d ceed-cuda
 //             * ex1 -pa -d ceed-hip
 //               ex1 -pa -d ceed-cuda:/gpu/cuda/shared
+//               ex1 -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
+//               ex1 -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
 //               ex1 -m ../data/beam-hex.mesh -pa -d cuda
 //               ex1 -m ../data/beam-tet.mesh -pa -d ceed-cpu
 //               ex1 -m ../data/beam-tet.mesh -pa -d ceed-cuda:/gpu/cuda/ref

examples/ex1p.cpp

@@ -30,14 +30,18 @@
 //
 // Device sample runs:
 //               mpirun -np 4 ex1p -pa -d cuda
-//               mpirun -np 4 ex1p -fa -d cuda
+//             * mpirun -np 4 ex1p -fa -d cuda
 //               mpirun -np 4 ex1p -pa -d occa-cuda
 //               mpirun -np 4 ex1p -pa -d raja-omp
 //               mpirun -np 4 ex1p -pa -d ceed-cpu
 //               mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
 //             * mpirun -np 4 ex1p -pa -d ceed-cuda
 //             * mpirun -np 4 ex1p -pa -d ceed-hip
 //               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
 //               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
 //
 // Description:  This example code demonstrates the use of MFEM to define a

examples/ex26.cpp

@@ -182,7 +182,7 @@ int main(int argc, char *argv[])
    }
    for (int level = 0; level < order_refinements; ++level)
    {
-      collections.Append(new H1_FECollection(std::pow(2, level+1), dim));
+      collections.Append(new H1_FECollection((int)std::pow(2, level+1), dim));
       fespaces.AddOrderRefinedLevel(collections.Last());
    }
 

examples/ex26p.cpp

@@ -219,7 +219,7 @@ int main(int argc, char *argv[])
    }
    for (int level = 0; level < order_refinements; ++level)
    {
-      collections.Append(new H1_FECollection(std::pow(2, level+1), dim));
+      collections.Append(new H1_FECollection((int)std::pow(2, level+1), dim));
       fespaces->AddOrderRefinedLevel(collections.Last());
    }
 

fem/bilinearform_ext.cpp

@@ -251,6 +251,7 @@ PABilinearFormExtension::PABilinearFormExtension(BilinearForm *form)
 
 void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
 {
+   if ( Device::Allows(Backend::CEED_MASK) ) { return; }
    ElementDofOrdering ordering = UsesTensorBasis(*a->FESpace())?
                                  ElementDofOrdering::LEXICOGRAPHIC:
                                  ElementDofOrdering::NATIVE;

fem/bilininteg.hpp

@@ -2423,7 +2423,6 @@ private:
    const DofToQuad *mapsO;         ///< Not owned. DOF-to-quad map, open.
    const DofToQuad *L2mapsO;       ///< Not owned. DOF-to-quad map, open.
    const DofToQuad *mapsC;         ///< Not owned. DOF-to-quad map, closed.
-   const GeometricFactors *geom;   ///< Not owned.
    int dim, ne, dofs1D, L2dofs1D, quad1D;
 
 public:

fem/bilininteg_convection_mf.cpp

@@ -30,7 +30,16 @@ void ConvectionIntegrator::AssembleMF(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::MFConvectionIntegrator(fes, *ir, Q, alpha);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedMFConvectionIntegrator(*this, fes, Q, alpha);
+      }
+      else
+      {
+         ceedOp = new ceed::MFConvectionIntegrator(fes, *ir, Q, alpha);
+      }
       return;
    }
    MFEM_ABORT("Error: ConvectionIntegrator::AssembleMF only implemented with"

fem/bilininteg_convection_pa.cpp

@@ -1386,7 +1386,16 @@ void ConvectionIntegrator::AssemblePA(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::PAConvectionIntegrator(fes, *ir, Q, alpha);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedPAConvectionIntegrator(*this, fes, Q, alpha);
+      }
+      else
+      {
+         ceedOp = new ceed::PAConvectionIntegrator(fes, *ir, Q, alpha);
+      }
       return;
    }
    const int dims = el.GetDim();

fem/bilininteg_diffusion_mf.cpp

@@ -33,7 +33,16 @@ void DiffusionIntegrator::AssembleMF(const FiniteElementSpace &fes)
       MFEM_VERIFY(!VQ && !MQ,
                   "Only scalar coefficient supported for DiffusionIntegrator"
                   " with libCEED");
-      ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedMFDiffusionIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
+      }
       return;
    }
    MFEM_ABORT("Error: DiffusionIntegrator::AssembleMF only implemented with"

fem/bilininteg_diffusion_pa.cpp

@@ -368,7 +368,16 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
       MFEM_VERIFY(!VQ && !MQ,
                   "Only scalar coefficient supported for DiffusionIntegrator"
                   " with libCEED");
-      ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedPADiffusionIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
+      }
       return;
    }
    const int dims = el.GetDim();

fem/bilininteg_hdiv.cpp

@@ -1799,15 +1799,9 @@ VectorFEDivergenceIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
    if (test_el->GetMapType() == FiniteElement::INTEGRAL)
    {
-      geom = mesh->GetGeometricFactors(*ir, GeometricFactors::DETERMINANTS);
-      for (int i=0; i<ne*nq; ++i)
-      {
-         coeff[i] /= geom->detJ[i];
-      }
-   }
-   else
-   {
-      geom = nullptr;
+      const GeometricFactors *geom =
+         mesh->GetGeometricFactors(*ir, GeometricFactors::DETERMINANTS);
+      coeff /= geom->detJ;
    }
 
    if (trial_el->GetDerivType() == mfem::FiniteElement::DIV && dim == 3)

fem/bilininteg_mass_mf.cpp

@@ -31,7 +31,16 @@ void MassIntegrator::AssembleMF(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedMFMassIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
+      }
       return;
    }
    MFEM_ABORT("Error: MassIntegrator::AssembleMF only implemented with"

fem/bilininteg_mass_pa.cpp

@@ -38,7 +38,16 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedPAMassIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
+      }
       return;
    }
    int map_type = el.GetMapType();

fem/bilininteg_vecdiffusion.cpp

@@ -149,7 +149,16 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedPADiffusionIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
+      }
       return;
    }
    const int dims = el.GetDim();

fem/bilininteg_vecdiffusion_mf.cpp

@@ -30,7 +30,19 @@ void VectorDiffusionIntegrator::AssembleMF(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
+      MFEM_VERIFY(!VQ && !MQ,
+                  "Only scalar coefficient supported for DiffusionIntegrator"
+                  " with libCEED");
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedMFDiffusionIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
+      }
       return;
    }
    MFEM_ABORT("Error: VectorDiffusionIntegrator::AssembleMF only implemented"

fem/bilininteg_vecmass.cpp

@@ -34,7 +34,16 @@ void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedPAMassIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
+      }
       return;
    }
    dim = mesh->Dimension();

fem/bilininteg_vecmass_mf.cpp

@@ -34,7 +34,16 @@ void VectorMassIntegrator::AssembleMF(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedMFMassIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
+      }
       return;
    }
    MFEM_ABORT("Error: VectorMassIntegrator::AssembleMF only implemented with"

fem/ceed/integrators/convection/convection.cpp

@@ -62,6 +62,20 @@ PAConvectionIntegrator::PAConvectionIntegrator(
 #endif
 }
 
+MixedPAConvectionIntegrator::MixedPAConvectionIntegrator(
+   const ConvectionIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::VectorCoefficient *Q,
+   const double alpha)
+{
+#ifdef MFEM_USE_CEED
+   ConvectionOperatorInfo info(fes.GetMesh()->Dimension(), alpha);
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 MFConvectionIntegrator::MFConvectionIntegrator(
    const mfem::FiniteElementSpace &fes,
    const mfem::IntegrationRule &irm,
@@ -77,6 +91,20 @@ MFConvectionIntegrator::MFConvectionIntegrator(
 #endif
 }
 
+MixedMFConvectionIntegrator::MixedMFConvectionIntegrator(
+   const ConvectionIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::VectorCoefficient *Q,
+   const double alpha)
+{
+#ifdef MFEM_USE_CEED
+   ConvectionOperatorInfo info(fes.GetMesh()->Dimension(), alpha);
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 } // namespace ceed
 
 } // namespace mfem

fem/ceed/integrators/convection/convection.hpp

@@ -13,6 +13,7 @@
 #define MFEM_LIBCEED_CONV_HPP
 
 #include "../../interface/integrator.hpp"
+#include "../../interface/mixed_integrator.hpp"
 #include "../../../fespace.hpp"
 
 namespace mfem
@@ -26,21 +27,39 @@ class PAConvectionIntegrator : public PAIntegrator
 {
 public:
    PAConvectionIntegrator(const mfem::FiniteElementSpace &fes,
-                          const mfem::IntegrationRule &irm,
+                          const mfem::IntegrationRule &ir,
                           mfem::VectorCoefficient *Q,
                           const double alpha);
 };
 
+class MixedPAConvectionIntegrator : public MixedIntegrator<PAIntegrator>
+{
+public:
+   MixedPAConvectionIntegrator(const ConvectionIntegrator &integ,
+                               const mfem::FiniteElementSpace &fes,
+                               mfem::VectorCoefficient *Q,
+                               const double alpha);
+};
+
 /// Represent a ConvectionIntegrator with AssemblyLevel::None using libCEED.
 class MFConvectionIntegrator : public MFIntegrator
 {
 public:
    MFConvectionIntegrator(const mfem::FiniteElementSpace &fes,
-                          const mfem::IntegrationRule &irm,
+                          const mfem::IntegrationRule &ir,
                           mfem::VectorCoefficient *Q,
                           const double alpha);
 };
 
+class MixedMFConvectionIntegrator : public MixedIntegrator<MFIntegrator>
+{
+public:
+   MixedMFConvectionIntegrator(const ConvectionIntegrator &integ,
+                               const mfem::FiniteElementSpace &fes,
+                               mfem::VectorCoefficient *Q,
+                               const double alpha);
+};
+
 }
 
 }

fem/ceed/integrators/diffusion/diffusion.cpp

@@ -60,6 +60,32 @@ PADiffusionIntegrator::PADiffusionIntegrator(
 #endif
 }
 
+MixedPADiffusionIntegrator::MixedPADiffusionIntegrator(
+   const DiffusionIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
+MixedPADiffusionIntegrator::MixedPADiffusionIntegrator(
+   const VectorDiffusionIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 MFDiffusionIntegrator::MFDiffusionIntegrator(
    const mfem::FiniteElementSpace &fes,
    const mfem::IntegrationRule &irm,
@@ -74,6 +100,32 @@ MFDiffusionIntegrator::MFDiffusionIntegrator(
 #endif
 }
 
+MixedMFDiffusionIntegrator::MixedMFDiffusionIntegrator(
+   const DiffusionIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
+MixedMFDiffusionIntegrator::MixedMFDiffusionIntegrator(
+   const VectorDiffusionIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 } // namespace ceed
 
 } // namespace mfem

fem/ceed/integrators/diffusion/diffusion.hpp

@@ -13,6 +13,7 @@
 #define MFEM_LIBCEED_DIFF_HPP
 
 #include "../../interface/integrator.hpp"
+#include "../../interface/mixed_integrator.hpp"
 #include "../../../fespace.hpp"
 
 namespace mfem
@@ -26,19 +27,43 @@ class PADiffusionIntegrator : public PAIntegrator
 {
 public:
    PADiffusionIntegrator(const mfem::FiniteElementSpace &fes,
-                         const mfem::IntegrationRule &irm,
+                         const mfem::IntegrationRule &ir,
                          mfem::Coefficient *Q);
 };
 
+class MixedPADiffusionIntegrator : public MixedIntegrator<PAIntegrator>
+{
+public:
+   MixedPADiffusionIntegrator(const DiffusionIntegrator &integ,
+                              const mfem::FiniteElementSpace &fes,
+                              mfem::Coefficient *Q);
+
+   MixedPADiffusionIntegrator(const VectorDiffusionIntegrator &integ,
+                              const mfem::FiniteElementSpace &fes,
+                              mfem::Coefficient *Q);
+};
+
 /// Represent a DiffusionIntegrator with AssemblyLevel::None using libCEED.
 class MFDiffusionIntegrator : public MFIntegrator
 {
 public:
    MFDiffusionIntegrator(const mfem::FiniteElementSpace &fes,
-                         const mfem::IntegrationRule &irm,
+                         const mfem::IntegrationRule &ir,
                          mfem::Coefficient *Q);
 };
 
+class MixedMFDiffusionIntegrator : public MixedIntegrator<MFIntegrator>
+{
+public:
+   MixedMFDiffusionIntegrator(const DiffusionIntegrator &integ,
+                              const mfem::FiniteElementSpace &fes,
+                              mfem::Coefficient *Q);
+
+   MixedMFDiffusionIntegrator(const VectorDiffusionIntegrator &integ,
+                              const mfem::FiniteElementSpace &fes,
+                              mfem::Coefficient *Q);
+};
+
 }
 
 }

fem/ceed/integrators/mass/mass.cpp

@@ -59,6 +59,30 @@ PAMassIntegrator::PAMassIntegrator(const mfem::FiniteElementSpace &fes,
 #endif
 }
 
+MixedPAMassIntegrator::MixedPAMassIntegrator(const MassIntegrator &integ,
+                                             const mfem::FiniteElementSpace &fes,
+                                             mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   MassOperatorInfo info;
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
+MixedPAMassIntegrator::MixedPAMassIntegrator(const VectorMassIntegrator &integ,
+                                             const mfem::FiniteElementSpace &fes,
+                                             mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   MassOperatorInfo info;
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 MFMassIntegrator::MFMassIntegrator(const mfem::FiniteElementSpace &fes,
                                    const mfem::IntegrationRule &irm,
                                    mfem::Coefficient *Q)
@@ -72,6 +96,30 @@ MFMassIntegrator::MFMassIntegrator(const mfem::FiniteElementSpace &fes,
 #endif
 }
 
+MixedMFMassIntegrator::MixedMFMassIntegrator(const MassIntegrator &integ,
+                                             const mfem::FiniteElementSpace &fes,
+                                             mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   MassOperatorInfo info;
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
+MixedMFMassIntegrator::MixedMFMassIntegrator(const VectorMassIntegrator &integ,
+                                             const mfem::FiniteElementSpace &fes,
+                                             mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   MassOperatorInfo info;
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 } // namespace ceed
 
 } // namespace mfem

fem/ceed/integrators/mass/mass.hpp

@@ -13,6 +13,7 @@
 #define MFEM_LIBCEED_MASS_HPP
 
 #include "../../interface/integrator.hpp"
+#include "../../interface/mixed_integrator.hpp"
 #include "../../../fespace.hpp"
 
 namespace mfem
@@ -26,19 +27,43 @@ class PAMassIntegrator : public PAIntegrator
 {
 public:
    PAMassIntegrator(const mfem::FiniteElementSpace &fes,
-                    const mfem::IntegrationRule &irm,
+                    const mfem::IntegrationRule &ir,
                     mfem::Coefficient *Q);
 };
 
+class MixedPAMassIntegrator : public MixedIntegrator<PAIntegrator>
+{
+public:
+   MixedPAMassIntegrator(const MassIntegrator &integ,
+                         const mfem::FiniteElementSpace &fes,
+                         mfem::Coefficient *Q);
+
+   MixedPAMassIntegrator(const VectorMassIntegrator &integ,
+                         const mfem::FiniteElementSpace &fes,
+                         mfem::Coefficient *Q);
+};
+
 /// Represent a MassIntegrator with AssemblyLevel::None using libCEED.
 class MFMassIntegrator : public MFIntegrator
 {
 public:
    MFMassIntegrator(const mfem::FiniteElementSpace &fes,
-                    const mfem::IntegrationRule &irm,
+                    const mfem::IntegrationRule &ir,
                     mfem::Coefficient *Q);
 };
 
+class MixedMFMassIntegrator : public MixedIntegrator<MFIntegrator>
+{
+public:
+   MixedMFMassIntegrator(const MassIntegrator &integ,
+                         const mfem::FiniteElementSpace &fes,
+                         mfem::Coefficient *Q);
+
+   MixedMFMassIntegrator(const VectorMassIntegrator &integ,
+                         const mfem::FiniteElementSpace &fes,
+                         mfem::Coefficient *Q);
+};
+
 }
 
 }

fem/ceed/integrators/nlconvection/nlconvection.cpp

@@ -60,6 +60,19 @@ PAVectorConvectionNLFIntegrator::PAVectorConvectionNLFIntegrator(
 #endif
 }
 
+MixedPAVectorConvectionNLIntegrator::MixedPAVectorConvectionNLIntegrator(
+   const VectorConvectionNLFIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   NLConvectionOperatorInfo info(fes.GetMesh()->Dimension());
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 MFVectorConvectionNLFIntegrator::MFVectorConvectionNLFIntegrator(
    const mfem::FiniteElementSpace &fes,
    const mfem::IntegrationRule &irm,
@@ -74,6 +87,19 @@ MFVectorConvectionNLFIntegrator::MFVectorConvectionNLFIntegrator(
 #endif
 }
 
+MixedMFVectorConvectionNLIntegrator::MixedMFVectorConvectionNLIntegrator(
+   const VectorConvectionNLFIntegrator &integ,
+   const mfem::FiniteElementSpace &fes,
+   mfem::Coefficient *Q)
+{
+#ifdef MFEM_USE_CEED
+   NLConvectionOperatorInfo info(fes.GetMesh()->Dimension());
+   Assemble(integ, info, fes, Q);
+#else
+   MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
+#endif
+}
+
 } // namespace ceed
 
 } // namespace mfem

fem/ceed/integrators/nlconvection/nlconvection.hpp

@@ -13,6 +13,7 @@
 #define MFEM_LIBCEED_NLCONV_HPP
 
 #include "../../interface/integrator.hpp"
+#include "../../interface/mixed_integrator.hpp"
 #include "../../../fespace.hpp"
 
 namespace mfem
@@ -31,6 +32,15 @@ public:
                                    mfem::Coefficient *coeff);
 };
 
+class MixedPAVectorConvectionNLIntegrator : public MixedIntegrator<PAIntegrator>
+{
+public:
+   MixedPAVectorConvectionNLIntegrator(
+      const VectorConvectionNLFIntegrator &integ,
+      const mfem::FiniteElementSpace &fes,
+      mfem::Coefficient *Q);
+};
+
 /** Represent a VectorConvectionNLFIntegrator with AssemblyLevel::None
     using libCEED. */
 class MFVectorConvectionNLFIntegrator : public MFIntegrator
@@ -41,6 +51,15 @@ public:
                                    mfem::Coefficient *coeff);
 };
 
+class MixedMFVectorConvectionNLIntegrator : public MixedIntegrator<MFIntegrator>
+{
+public:
+   MixedMFVectorConvectionNLIntegrator(
+      const VectorConvectionNLFIntegrator &integ,
+      const mfem::FiniteElementSpace &fes,
+      mfem::Coefficient *Q);
+};
+
 }
 
 }

fem/ceed/integrators/nlconvection/nlconvection_qf.h

@@ -327,13 +327,13 @@ CEED_QFUNCTION(f_apply_conv_mf_const)(void *ctx, CeedInt Q,
             const CeedScalar A33 = J11 * J22 - J12 * J21;
             const CeedScalar w = qw[i] * coeff;
             const CeedScalar qd00 = w * A11;
-            const CeedScalar qd01 = w * A21;
-            const CeedScalar qd02 = w * A31;
-            const CeedScalar qd10 = w * A12;
+            const CeedScalar qd10 = w * A21;
+            const CeedScalar qd20 = w * A31;
+            const CeedScalar qd01 = w * A12;
             const CeedScalar qd11 = w * A22;
-            const CeedScalar qd12 = w * A32;
-            const CeedScalar qd20 = w * A13;
-            const CeedScalar qd21 = w * A23;
+            const CeedScalar qd21 = w * A32;
+            const CeedScalar qd02 = w * A13;
+            const CeedScalar qd12 = w * A23;
             const CeedScalar qd22 = w * A33;
             const CeedScalar u0   = u[i + Q * 0];
             const CeedScalar u1   = u[i + Q * 1];
@@ -440,13 +440,13 @@ CEED_QFUNCTION(f_apply_conv_mf_quad)(void *ctx, CeedInt Q,
             const CeedScalar A33 = J11 * J22 - J12 * J21;
             const CeedScalar w = qw[i] * c[i];
             const CeedScalar qd00 = w * A11;
-            const CeedScalar qd01 = w * A21;
-            const CeedScalar qd02 = w * A31;
-            const CeedScalar qd10 = w * A12;
+            const CeedScalar qd10 = w * A21;
+            const CeedScalar qd20 = w * A31;
+            const CeedScalar qd01 = w * A12;
             const CeedScalar qd11 = w * A22;
-            const CeedScalar qd12 = w * A32;
-            const CeedScalar qd20 = w * A13;
-            const CeedScalar qd21 = w * A23;
+            const CeedScalar qd21 = w * A32;
+            const CeedScalar qd02 = w * A13;
+            const CeedScalar qd12 = w * A23;
             const CeedScalar qd22 = w * A33;
             const CeedScalar u0   = u[i + Q * 0];
             const CeedScalar u1   = u[i + Q * 1];

fem/ceed/interface/basis.cpp

@@ -36,6 +36,8 @@ static CeedElemTopology GetCeedTopology(Geometry::Type geom)
          return CEED_TOPOLOGY_HEX;
       case Geometry::PRISM:
          return CEED_TOPOLOGY_PRISM;
+      case Geometry::PYRAMID:
+         return CEED_TOPOLOGY_PYRAMID;
       default:
          MFEM_ABORT("This type of element is not supported");
          return CEED_TOPOLOGY_PRISM; // Silence warning
@@ -43,11 +45,11 @@ static CeedElemTopology GetCeedTopology(Geometry::Type geom)
 }
 
 static void InitNonTensorBasis(const mfem::FiniteElementSpace &fes,
+                               const mfem::FiniteElement &fe,
                                const mfem::IntegrationRule &ir,
                                Ceed ceed, CeedBasis *basis)
 {
-   const mfem::DofToQuad &maps = fes.GetFE(0)->
-                                 GetDofToQuad(ir,mfem::DofToQuad::FULL);
+   const mfem::DofToQuad &maps = fe.GetDofToQuad(ir, mfem::DofToQuad::FULL);
    mfem::Mesh *mesh = fes.GetMesh();
    const int dim = mesh->Dimension();
    const int ndofs = maps.ndof;
@@ -62,18 +64,18 @@ static void InitNonTensorBasis(const mfem::FiniteElementSpace &fes,
       if (dim>2) { qX(2,i) = ip.z; }
       qW(i) = ip.weight;
    }
-   CeedBasisCreateH1(ceed, GetCeedTopology(fes.GetFE(0)->GetGeomType()),
+   CeedBasisCreateH1(ceed, GetCeedTopology(fe.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::FiniteElement &fe,
                             const mfem::IntegrationRule &ir,
                             Ceed ceed, CeedBasis *basis)
 {
-   const mfem::DofToQuad &maps =
-      fes.GetFE(0)->GetDofToQuad(ir, mfem::DofToQuad::TENSOR);
+   const mfem::DofToQuad &maps = fe.GetDofToQuad(ir, mfem::DofToQuad::TENSOR);
    mfem::Mesh *mesh = fes.GetMesh();
    const int ndofs = maps.ndof;
    const int nqpts = maps.nqpt;
@@ -96,28 +98,30 @@ static void InitTensorBasis(const mfem::FiniteElementSpace &fes,
                            qW.GetData(), basis);
 }
 
-void InitBasis(const FiniteElementSpace &fes,
-               const IntegrationRule &irm,
-               Ceed ceed, CeedBasis *basis)
+static void InitBasisImpl(const FiniteElementSpace &fes,
+                          const FiniteElement &fe,
+                          const IntegrationRule &ir,
+                          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 P = fe.GetDof();
+   const int Q = ir.GetNPoints();
    const int ncomp = fes.GetVDim();
-   BasisKey basis_key(&fes, &irm, ncomp, P, Q);
+   BasisKey basis_key(&fes, &ir, ncomp, P, Q);
    auto basis_itr = mfem::internal::ceed_basis_map.find(basis_key);
+   const bool tensor = dynamic_cast<const mfem::TensorBasisElement *>
+                       (&fe) != nullptr;
 
    // Init or retreive key values
    if (basis_itr == mfem::internal::ceed_basis_map.end())
    {
-      if (UsesTensorBasis(fes))
+      if ( tensor )
       {
-         InitTensorBasis(fes, irm, ceed, basis);
+         InitTensorBasis(fes, fe, ir, ceed, basis);
       }
       else
       {
-         InitNonTensorBasis(fes, irm, ceed, basis);
+         InitNonTensorBasis(fes, fe, ir, ceed, basis);
       }
       mfem::internal::ceed_basis_map[basis_key] = *basis;
    }
@@ -127,6 +131,24 @@ void InitBasis(const FiniteElementSpace &fes,
    }
 }
 
+void InitBasis(const FiniteElementSpace &fes,
+               const IntegrationRule &ir,
+               Ceed ceed, CeedBasis *basis)
+{
+   const mfem::FiniteElement &fe = *fes.GetFE(0);
+   InitBasisImpl(fes, fe, ir, ceed, basis);
+}
+
+void InitBasisWithIndices(const FiniteElementSpace &fes,
+                          const IntegrationRule &ir,
+                          int nelem,
+                          const int* indices,
+                          Ceed ceed, CeedBasis *basis)
+{
+   const mfem::FiniteElement &fe = *fes.GetFE(indices[0]);
+   InitBasisImpl(fes, fe, ir, ceed, basis);
+}
+
 #endif
 
 } // namespace ceed

fem/ceed/interface/basis.hpp

@@ -22,17 +22,32 @@ namespace ceed
 
 #ifdef MFEM_USE_CEED
 
-/** @brief Initialize a CeedBasis.
+/** @brief Initialize a CeedBasis for non-mixed meshes.
 
    @param[in] fes Input finite element space.
-   @param[in] irm Input integration rule.
+   @param[in] ir 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,
+               const IntegrationRule &ir,
                Ceed ceed, CeedBasis *basis);
 
+/** @brief Initialize a CeedBasis for mixed meshes.
+
+    @param[in] fes The finite element space.
+    @param[in] ir is the integration rule for the operator.
+    @param[in] nelem The number of elements.
+    @param[in] indices The indices of the elements of same type in the
+                       `FiniteElementSpace`.
+    @param[in] ceed The Ceed object.
+    @param[out] basis The `CeedBasis` to initialize. */
+void InitBasisWithIndices(const FiniteElementSpace &fes,
+                          const IntegrationRule &ir,
+                          int nelem,
+                          const int* indices,
+                          Ceed ceed, CeedBasis *basis);
+
 #endif
 
 } // namespace ceed

fem/ceed/interface/coefficient.hpp

@@ -14,6 +14,7 @@
 
 #ifdef MFEM_USE_CEED
 
+#include "../../../general/forall.hpp"
 #include "../../../config/config.hpp"
 #include "../../../linalg/vector.hpp"
 #include "../../../linalg/dtensor.hpp"
@@ -77,7 +78,14 @@ struct QuadCoefficient : VariableCoefficient
 
 /** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
     mfem::Coefficient @a Q, an mfem::Mesh @a mesh, and an mfem::IntegrationRule
-    @a ir. */
+    @a ir.
+
+    @param[in] Q is the coefficient from the `Integrator`.
+    @param[in] mesh is the mesh.
+    @param[in] ir is the integration rule.
+    @param[out] coeff_ptr is the structure to store the coefficient for the
+                          `CeedOperator`.
+    @param[out] ctx is the Context associated to the QFunction. */
 template <typename Context>
 void InitCoefficient(mfem::Coefficient *Q, mfem::Mesh &mesh,
                      const mfem::IntegrationRule &ir,
@@ -143,8 +151,15 @@ void InitCoefficient(mfem::Coefficient *Q, mfem::Mesh &mesh,
 
 
 /** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
-    mfem::VectorCoefficient @a Q, an mfem::Mesh @a mesh, and an
-    mfem::IntegrationRule @a ir. */
+    mfem::VectorCoefficient @a VQ, an mfem::Mesh @a mesh, and an
+    mfem::IntegrationRule @a ir.
+
+    @param[in] VQ is the vector coefficient from the `Integrator`.
+    @param[in] mesh is the mesh.
+    @param[in] ir is the integration rule.
+    @param[out] coeff_ptr is the structure to store the coefficient for the
+                          `CeedOperator`.
+    @param[out] ctx is the Context associated to the QFunction. */
 template <typename Context>
 void InitCoefficient(mfem::VectorCoefficient *VQ, mfem::Mesh &mesh,
                      const mfem::IntegrationRule &ir,
@@ -214,6 +229,209 @@ void InitCoefficient(mfem::VectorCoefficient *VQ, mfem::Mesh &mesh,
    }
 }
 
+/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
+    mfem::Coefficient @a Q, an mfem::Mesh @a mesh, and an mfem::IntegrationRule
+    @a ir for the elements given by the indices @a indices.
+
+    @param[in] Q is the coefficient from the `Integrator`.
+    @param[in] mesh is the mesh.
+    @param[in] ir is the integration rule.
+    @param[in] nelem The number of elements.
+    @param[in] indices The indices of the elements of same type in the
+                       `FiniteElementSpace`.
+    @param[out] coeff_ptr is the structure to store the coefficient for the
+                          `CeedOperator`.
+    @param[out] ctx is the Context associated to the QFunction. */
+template <typename Context>
+void InitCoefficientWithIndices(mfem::Coefficient *Q, mfem::Mesh &mesh,
+                                const mfem::IntegrationRule &ir,
+                                int nelem,
+                                const int* indices,
+                                Coefficient*& coeff_ptr, Context &ctx)
+{
+   if ( Q == nullptr )
+   {
+      Coefficient *ceedCoeff = new Coefficient(1);
+      ctx.coeff = 1.0;
+      coeff_ptr = ceedCoeff;
+   }
+   else if (ConstantCoefficient *const_coeff =
+               dynamic_cast<ConstantCoefficient*>(Q))
+   {
+      Coefficient *ceedCoeff = new Coefficient(1);
+      ctx.coeff = const_coeff->constant;
+      coeff_ptr = ceedCoeff;
+   }
+   else if (GridFunctionCoefficient* gf_coeff =
+               dynamic_cast<GridFunctionCoefficient*>(Q))
+   {
+      GridCoefficient *ceedCoeff =
+         new GridCoefficient(*gf_coeff->GetGridFunction());
+      coeff_ptr = ceedCoeff;
+   }
+   else if (QuadratureFunctionCoefficient *cQ =
+               dynamic_cast<QuadratureFunctionCoefficient*>(Q))
+   {
+      QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
+      const int ne = mesh.GetNE();
+      const int nq = ir.GetNPoints();
+      const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
+      MFEM_VERIFY(qFun.Size() == nq * ne,
+                  "Incompatible QuadratureFunction dimension \n");
+
+      MFEM_VERIFY(&ir == &qFun.GetSpace()->GetElementIntRule(0),
+                  "IntegrationRule used within integrator and in"
+                  " QuadratureFunction appear to be different");
+      ceedCoeff->coeff.SetSize(nq * nelem);
+      Memory<int> m_indices((int*)indices, nelem, false);
+      auto in = Reshape(qFun.Read(), nq, ne);
+      auto d_indices = Read(m_indices, nelem);
+      auto out = Reshape(ceedCoeff->coeff.Write(), nq, nelem);
+      MFEM_FORALL(i, nelem * nq,
+      {
+         const int q = i%nq;
+         const int sub_e = i/nq;
+         const int e = d_indices[sub_e];
+         out(q, sub_e) = in(q, e);
+      });
+      m_indices.DeleteDevice();
+      InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
+      coeff_ptr = ceedCoeff;
+   }
+   else
+   {
+      QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
+      const int nq = ir.GetNPoints();
+      ceedCoeff->coeff.SetSize(nq * nelem);
+      auto C = Reshape(ceedCoeff->coeff.HostWrite(), nq, nelem);
+      for (int i = 0; i < nelem; ++i)
+      {
+         const int e = indices[i];
+         mfem::ElementTransformation &T = *mesh.GetElementTransformation(e);
+         for (int q = 0; q < nq; ++q)
+         {
+            C(q, i) = Q->Eval(T, ir.IntPoint(q));
+         }
+      }
+      InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
+      coeff_ptr = ceedCoeff;
+   }
+}
+
+
+/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
+    mfem::VectorCoefficient @a Q, an mfem::Mesh @a mesh, and an
+    mfem::IntegrationRule @a ir for the elements given by the indices @a indices.
+
+    @param[in] VQ is the vector coefficient from the `Integrator`.
+    @param[in] mesh is the mesh.
+    @param[in] ir is the integration rule.
+    @param[in] nelem The number of elements.
+    @param[in] indices The indices of the elements of same type in the
+                       `FiniteElementSpace`.
+    @param[out] coeff_ptr is the structure to store the coefficient for the
+                          `CeedOperator`.
+    @param[out] ctx is the Context associated to the QFunction. */
+template <typename Context>
+void InitCoefficientWithIndices(mfem::VectorCoefficient *VQ, mfem::Mesh &mesh,
+                                const mfem::IntegrationRule &ir,
+                                int nelem,
+                                const int* indices,
+                                Coefficient *&coeff_ptr, Context &ctx)
+{
+   if (VectorConstantCoefficient *const_coeff =
+          dynamic_cast<VectorConstantCoefficient*>(VQ))
+   {
+      const int vdim = const_coeff->GetVDim();
+      const mfem::Vector &val = const_coeff->GetVec();
+      Coefficient *ceedCoeff = new Coefficient(vdim);
+      for (int i = 0; i < vdim; i++)
+      {
+         ctx.coeff[i] = val[i];
+      }
+      coeff_ptr = ceedCoeff;
+   }
+   else if (VectorGridFunctionCoefficient* vgf_coeff =
+               dynamic_cast<VectorGridFunctionCoefficient*>(VQ))
+   {
+      GridCoefficient *ceedCoeff =
+         new GridCoefficient(*vgf_coeff->GetGridFunction());
+      coeff_ptr = ceedCoeff;
+   }
+   else if (VectorQuadratureFunctionCoefficient *cQ =
+               dynamic_cast<VectorQuadratureFunctionCoefficient*>(VQ))
+   {
+      QuadCoefficient *ceedCoeff = new QuadCoefficient(cQ->GetVDim());
+      const int dim = mesh.Dimension();
+      const int ne = mesh.GetNE();
+      const int nq = ir.GetNPoints();
+      const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
+      MFEM_VERIFY(qFun.Size() == dim * nq * ne,
+                  "Incompatible QuadratureFunction dimension \n");
+
+      MFEM_VERIFY(&ir == &qFun.GetSpace()->GetElementIntRule(0),
+                  "IntegrationRule used within integrator and in"
+                  " QuadratureFunction appear to be different");
+      ceedCoeff->coeff.SetSize(dim * nq * nelem);
+      Memory<int> m_indices((int*)indices, nelem, false);
+      auto in = Reshape(qFun.Read(), dim, nq, ne);
+      auto d_indices = Read(m_indices, nelem);
+      auto out = Reshape(ceedCoeff->coeff.Write(), dim, nq, nelem);
+      MFEM_FORALL(i, nelem * nq,
+      {
+         const int q = i%nq;
+         const int sub_e = i/nq;
+         const int e = d_indices[sub_e];
+         for (int d = 0; d < dim; d++)
+         {
+            out(d, q, sub_e) = in(d, q, e);
+         }
+      });
+      m_indices.DeleteDevice();
+      InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
+      coeff_ptr = ceedCoeff;
+   }
+   else
+   {
+      const int dim = mesh.Dimension();
+      QuadCoefficient *ceedCoeff = new QuadCoefficient(dim);
+      const int nq = ir.GetNPoints();
+      ceedCoeff->coeff.SetSize(dim * nq * nelem);
+      auto C = Reshape(ceedCoeff->coeff.HostWrite(), dim, nq, nelem);
+      mfem::DenseMatrix Q_ir;
+      for (int i = 0; i < nelem; ++i)
+      {
+         const int e = indices[i];
+         mfem::ElementTransformation &T = *mesh.GetElementTransformation(e);
+         VQ->Eval(Q_ir, T, ir);
+         for (int q = 0; q < nq; ++q)
+         {
+            for (int d = 0; d < dim; ++d)
+            {
+               C(d, q, i) = Q_ir(d, q);
+            }
+         }
+      }
+      InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
+      coeff_ptr = ceedCoeff;
+   }
+}
+
+template <typename Coeff, typename Context>
+void InitCoefficient(Coeff *Q, mfem::Mesh &mesh,
+                     const mfem::IntegrationRule &ir, int nelem,
+                     const int* indices, Coefficient *&coeff_ptr, Context &ctx)
+{
+   if (indices)
+   {
+      InitCoefficientWithIndices(Q, mesh, ir, nelem, indices, coeff_ptr, ctx);
+   }
+   else
+   {
+      InitCoefficient(Q, mesh, ir, coeff_ptr, ctx);
+   }
+}
+
 } // namespace ceed
 
 } // namespace mfem

fem/ceed/interface/integrator.hpp

@@ -18,6 +18,7 @@
 #include "operator.hpp"
 #include "coefficient.hpp"
 #include "restriction.hpp"
+#include "util.hpp"
 #include "ceed.hpp"
 
 namespace mfem
@@ -86,6 +87,7 @@ protected:
    CeedQFunctionContext build_ctx;
    CeedOperator build_oper;
 
+public:
    PAIntegrator()
       : Operator(),
         trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
@@ -95,23 +97,51 @@ protected:
         qdata(nullptr), coeff(nullptr), build_ctx(nullptr), build_oper(nullptr)
    { }
 
-public:
-   /** This method assembles the PAIntegrator.
+   /** @brief 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.
 
-       @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`. */
+       @param[in] info is the structure describing the CeedOperator to assemble.
+       @param[in] fes is the finite element space.
+       @param[in] ir is the integration rule for the operator.
+       @param[in] Q is the coefficient from the `Integrator`. */
    template <typename CeedOperatorInfo, typename CoeffType>
    void Assemble(CeedOperatorInfo &info,
                  const mfem::FiniteElementSpace &fes,
-                 const mfem::IntegrationRule &irm,
+                 const mfem::IntegrationRule &ir,
                  CoeffType *Q)
    {
-      Assemble(info, fes, fes, irm, Q);
+      Assemble(info, fes, ir, fes.GetNE(), nullptr, Q);
+   }
+
+   /** @brief 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 for the elements given by the indices
+       @a indices.
+       The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
+       and contain a `Context` type relevant to the qFunctions.
+
+       @param[in] info is the structure describing the CeedOperator to assemble.
+       @param[in] fes is the finite element space.
+       @param[in] ir is the integration rule for the operator.
+       @param[in] nelem The number of elements.
+       @param[in] indices The indices of the elements of same type in the
+                          `FiniteElementSpace`. If `indices == nullptr`, assumes
+                          that the `FiniteElementSpace` is not mixed.
+       @param[in] Q is the coefficient from the `Integrator`. */
+   template <typename CeedOperatorInfo, typename CoeffType>
+   void Assemble(CeedOperatorInfo &info,
+                 const mfem::FiniteElementSpace &fes,
+                 const mfem::IntegrationRule &ir,
+                 int nelem,
+                 const int* indices,
+                 CoeffType *Q)
+   {
+      Assemble(info, fes, fes, ir, nelem, indices, Q);
    }
 
    /** This method assembles the PAIntegrator for mixed forms.
@@ -128,12 +158,40 @@ public:
    void Assemble(CeedOperatorInfo &info,
                  const mfem::FiniteElementSpace &trial_fes,
                  const mfem::FiniteElementSpace &test_fes,
-                 const mfem::IntegrationRule &irm,
+                 const mfem::IntegrationRule &ir,
+                 CoeffType *Q)
+   {
+      Assemble(info, trial_fes, test_fes, ir, trial_fes.GetNE(), nullptr, Q);
+   }
+
+   /** This method assembles the PAIntegrator for mixed forms on mixed meshes.
+
+       @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] nelem The number of elements,
+       @param[in] indices The indices of the elements of same type in the
+                          `FiniteElementSpace`. If `indices == nullptr`, assumes
+                          that the `FiniteElementSpace` is not mixed,
+       @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 &ir,
+                 int nelem,
+                 const int* indices,
                  CoeffType *Q)
    {
       Ceed ceed(internal::ceed);
       mfem::Mesh &mesh = *trial_fes.GetMesh();
-      InitCoefficient(Q, mesh, irm, coeff, info.ctx);
+      MFEM_VERIFY(!(!indices && mesh.GetNumGeometries(mesh.Dimension()) > 1),
+                  "Use ceed::MixedIntegrator on mixed meshes.");
+      InitCoefficient(Q, mesh, ir, nelem, indices, coeff, info.ctx);
       bool const_coeff = coeff->IsConstant();
       std::string build_func = const_coeff ? info.build_func_const
                                : info.build_func_quad;
@@ -145,7 +203,6 @@ public:
                      info.trial_op,
                      info.test_op
                     };
-      CeedInt nqpts, nelem = mesh.GetNE();
       CeedInt dim = mesh.SpaceDimension();
       CeedInt trial_vdim = trial_fes.GetVDim();
       CeedInt test_vdim = test_fes.GetVDim();
@@ -153,23 +210,23 @@ public:
       mesh.EnsureNodes();
       if ( &trial_fes == &test_fes )
       {
-         InitBasisAndRestriction(trial_fes, irm, ceed,
-                                 &trial_basis, &trial_restr);
+         InitBasisAndRestriction(trial_fes, ir, nelem, indices,
+                                 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);
+         InitBasisAndRestriction(trial_fes, ir, nelem, indices,
+                                 ceed, &trial_basis, &trial_restr);
+         InitBasisAndRestriction(test_fes, ir, nelem, indices,
+                                 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);
+      InitBasisAndRestriction(*mesh_fes, ir, nelem, indices,
+                              ceed, &mesh_basis, &mesh_restr);
 
       CeedInt trial_nqpts, test_nqpts;
       CeedBasisGetNumQuadraturePoints(trial_basis, &trial_nqpts);
@@ -177,7 +234,7 @@ public:
       MFEM_VERIFY(trial_nqpts == test_nqpts,
                   "Trial and test basis must have the same number of quadrature"
                   " points.");
-      nqpts = trial_nqpts;
+      CeedInt nqpts = trial_nqpts;
 
       const int qdatasize = op.qdatasize;
       InitStridedRestriction(*mesh_fes, nelem, nqpts, qdatasize,
@@ -221,8 +278,10 @@ public:
       CeedOperatorCreate(ceed, build_qfunc, NULL, NULL, &build_oper);
       if (GridCoefficient *gridCoeff = dynamic_cast<GridCoefficient*>(coeff))
       {
-         InitBasisAndRestriction(*gridCoeff->gf.FESpace(), irm, ceed,
-                                 &gridCoeff->basis, &gridCoeff->restr);
+         InitBasisAndRestriction(*gridCoeff->gf.FESpace(), ir,
+                                 nelem, indices, ceed,
+                                 &gridCoeff->basis,
+                                 &gridCoeff->restr);
          CeedOperatorSetField(build_oper, "coeff", gridCoeff->restr,
                               gridCoeff->basis, gridCoeff->coeffVector);
       }
@@ -231,7 +290,8 @@ public:
       {
          const int ncomp = quadCoeff->ncomp;
          CeedInt strides[3] = {ncomp, 1, ncomp*nqpts};
-         InitStridedRestriction(*mesh_fes, nelem, nqpts, ncomp, strides,
+         InitStridedRestriction(*mesh.GetNodalFESpace(),
+                                nelem, nqpts, ncomp, strides,
                                 &quadCoeff->restr);
          CeedOperatorSetField(build_oper, "coeff", quadCoeff->restr,
                               CEED_BASIS_COLLOCATED, quadCoeff->coeffVector);
@@ -254,22 +314,17 @@ public:
       switch (op.trial_op)
       {
          case EvalMode::None:
-            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
-                                  CEED_EVAL_NONE);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_NONE);
             break;
          case EvalMode::Interp:
-            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
-                                  CEED_EVAL_INTERP);
+            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_INTERP);
             break;
          case EvalMode::Grad:
-            CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
-                                  CEED_EVAL_GRAD);
+            CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim, CEED_EVAL_GRAD);
             break;
          case EvalMode::InterpAndGrad:
-            CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
-                                  CEED_EVAL_INTERP);
-            CeedQFunctionAddInput(apply_qfunc, "gu", trial_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
@@ -278,22 +333,17 @@ public:
       switch (op.test_op)
       {
          case EvalMode::None:
-            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
-                                   CEED_EVAL_NONE);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_NONE);
             break;
          case EvalMode::Interp:
-            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
-                                   CEED_EVAL_INTERP);
+            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_INTERP);
             break;
          case EvalMode::Grad:
-            CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
-                                   CEED_EVAL_GRAD);
+            CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim, CEED_EVAL_GRAD);
             break;
          case EvalMode::InterpAndGrad:
-            CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
-                                   CEED_EVAL_INTERP);
-            CeedQFunctionAddOutput(apply_qfunc, "gv", test_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);
@@ -308,18 +358,14 @@ public:
                                  CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
             break;
          case EvalMode::Interp:
-            CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
-                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "u", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
             break;
          case EvalMode::Grad:
-            CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
-                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gu", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
             break;
          case EvalMode::InterpAndGrad:
-            CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
-                                 CEED_VECTOR_ACTIVE);
-            CeedOperatorSetField(oper, "gu", trial_restr, trial_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
@@ -333,18 +379,14 @@ public:
                                  CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
             break;
          case EvalMode::Interp:
-            CeedOperatorSetField(oper, "v", test_restr, test_basis,
-                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "v", test_restr, test_basis, CEED_VECTOR_ACTIVE);
             break;
          case EvalMode::Grad:
-            CeedOperatorSetField(oper, "gv", test_restr, test_basis,
-                                 CEED_VECTOR_ACTIVE);
+            CeedOperatorSetField(oper, "gv", test_restr, test_basis, CEED_VECTOR_ACTIVE);
             break;
          case EvalMode::InterpAndGrad:
-            CeedOperatorSetField(oper, "v", test_restr, test_basis,
-                                 CEED_VECTOR_ACTIVE);
-            CeedOperatorSetField(oper, "gv", test_restr, test_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;
       }
 
@@ -402,6 +444,7 @@ protected:
    Coefficient *coeff;
    CeedQFunctionContext build_ctx;
 
+public:
    MFIntegrator()
       : Operator(),
         trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
@@ -410,23 +453,51 @@ protected:
         apply_qfunc(nullptr), node_coords(nullptr),
         qdata(nullptr), coeff(nullptr), build_ctx(nullptr) { }
 
-public:
-   /** This method assembles the MFIntegrator.
+   /** @brief 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.
 
-       @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`. */
+       @param[in] info is the structure describing the CeedOperator to assemble.
+       @param[in] fes is the finite element space.
+       @param[in] ir is the integration rule for the operator.
+       @param[in] Q is the coefficient from the `Integrator`. */
    template <typename CeedOperatorInfo, typename CoeffType>
    void Assemble(CeedOperatorInfo &info,
                  const mfem::FiniteElementSpace &fes,
-                 const mfem::IntegrationRule &irm,
+                 const mfem::IntegrationRule &ir,
                  CoeffType *Q)
    {
-      Assemble(info, fes, fes, irm, Q);
+      Assemble(info, fes, ir, fes.GetNE(), nullptr, Q);
+   }
+
+   /** @brief 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 for the elements given by the indices
+       @a indices.
+       The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
+       and contain a `Context` type relevant to the qFunctions.
+
+       @param[in] info is the structure describing the CeedOperator to assemble.
+       @param[in] fes is the finite element space.
+       @param[in] ir is the integration rule for the operator.
+       @param[in] nelem The number of elements.
+       @param[in] indices The indices of the elements of same type in the
+                          `FiniteElementSpace`. If `indices == nullptr`, assumes
+                          that the `FiniteElementSpace` is not mixed.
+       @param[in] Q is the coefficient from the `Integrator`. */
+   template <typename CeedOperatorInfo, typename CoeffType>
+   void Assemble(CeedOperatorInfo &info,
+                 const mfem::FiniteElementSpace &fes,
+                 const mfem::IntegrationRule &ir,
+                 int nelem,
+                 const int* indices,
+                 CoeffType *Q)
+   {
+      Assemble(info, fes, fes, ir, nelem, indices, Q);
    }
 
    /** This method assembles the MFIntegrator for mixed forms.
@@ -443,12 +514,40 @@ public:
    void Assemble(CeedOperatorInfo &info,
                  const mfem::FiniteElementSpace &trial_fes,
                  const mfem::FiniteElementSpace &test_fes,
-                 const mfem::IntegrationRule &irm,
+                 const mfem::IntegrationRule &ir,
+                 CoeffType *Q)
+   {
+      Assemble(info, trial_fes, test_fes, ir, trial_fes.GetNE(), nullptr, 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] nelem The number of elements,
+       @param[in] indices The indices of the elements of same type in the
+                          `FiniteElementSpace`. If `indices == nullptr`, assumes
+                          that the `FiniteElementSpace` is not mixed,
+       @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 &ir,
+                 int nelem,
+                 const int* indices,
                  CoeffType *Q)
    {
       Ceed ceed(internal::ceed);
       Mesh &mesh = *trial_fes.GetMesh();
-      InitCoefficient(Q, mesh, irm, coeff, info.ctx);
+      MFEM_VERIFY(!(!indices && mesh.GetNumGeometries(mesh.Dimension()) > 1),
+                  "Use ceed::MixedIntegrator on mixed meshes.");
+      InitCoefficient(Q, mesh, ir, nelem, indices, coeff, info.ctx);
       bool const_coeff = coeff->IsConstant();
       std::string apply_func = const_coeff ? info.apply_func_mf_const
                                : info.apply_func_mf_quad;
@@ -459,7 +558,7 @@ public:
                      info.trial_op,
                      info.test_op
                     };
-      CeedInt nqpts, nelem = mesh.GetNE();
+
       CeedInt dim = mesh.SpaceDimension();
       CeedInt trial_vdim = trial_fes.GetVDim();
       CeedInt test_vdim = test_fes.GetVDim();
@@ -467,22 +566,22 @@ public:
       mesh.EnsureNodes();
       if ( &trial_fes == &test_fes )
       {
-         InitBasisAndRestriction(trial_fes, irm, ceed,
+         InitBasisAndRestriction(trial_fes, ir, nelem, indices, ceed,
                                  &trial_basis, &trial_restr);
          test_basis = trial_basis;
          test_restr = trial_restr;
       }
       else
       {
-         InitBasisAndRestriction(trial_fes, irm, ceed,
+         InitBasisAndRestriction(trial_fes, ir, nelem, indices, ceed,
                                  &trial_basis, &trial_restr);
-         InitBasisAndRestriction(test_fes, irm, ceed,
+         InitBasisAndRestriction(test_fes, ir, nelem, indices, 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,
+      InitBasisAndRestriction(*mesh_fes, ir, nelem, indices, ceed, &mesh_basis,
                               &mesh_restr);
 
       CeedInt trial_nqpts, test_nqpts;
@@ -491,7 +590,7 @@ public:
       MFEM_VERIFY(trial_nqpts == test_nqpts,
                   "Trial and test basis must have the same number of quadrature"
                   " points.");
-      nqpts = trial_nqpts;
+      CeedInt nqpts = trial_nqpts;
 
       InitVector(*mesh.GetNodes(), node_coords);
 
@@ -572,8 +671,8 @@ public:
       // coefficient
       if (GridCoefficient *gridCoeff = dynamic_cast<GridCoefficient*>(coeff))
       {
-         InitBasisAndRestriction(*gridCoeff->gf.FESpace(), irm, ceed,
-                                 &gridCoeff->basis, &gridCoeff->restr);
+         InitBasisAndRestriction(*gridCoeff->gf.FESpace(), ir, nelem, indices,
+                                 ceed, &gridCoeff->basis, &gridCoeff->restr);
          CeedOperatorSetField(oper, "coeff", gridCoeff->restr,
                               gridCoeff->basis, gridCoeff->coeffVector);
       }

fem/ceed/interface/interface.hpp

@@ -22,6 +22,8 @@
 #include "coefficient.hpp"
 // PA or MF Operator using libCEED.
 #include "integrator.hpp"
+// PA Operator supporting mixed finite element spaces.
+#include "mixed_integrator.hpp"
 // Utility functions
 #include "util.hpp"
 // Wrapper to include <ceed.h>

fem/ceed/interface/mixed_integrator.hpp

@@ -0,0 +1,126 @@
+// 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_MIXED_INTEGRATOR
+#define MFEM_LIBCEED_MIXED_INTEGRATOR
+
+#include "ceed.hpp"
+#include "integrator.hpp"
+#include <unordered_map>
+
+namespace mfem
+{
+
+namespace ceed
+{
+
+/** @brief This class wraps a `ceed::PAIntegrator` or `ceed::MFIntegrator` to
+    support mixed finite element spaces. */
+template <typename CeedInteg>
+class MixedIntegrator : public ceed::Operator
+{
+#ifdef MFEM_USE_CEED
+   using ElementKey = std::pair<int, int>; //< Element::Type, Order >
+   struct key_hash
+   {
+      std::size_t operator()(const ElementKey& k) const
+      {
+         return k.first + 2 * k.second;
+      }
+   };
+   using ElementsMap = std::unordered_map<const ElementKey, int*, key_hash>;
+   std::vector<CeedInteg*> sub_ops;
+
+public:
+   template <typename Integrator, typename CeedOperatorInfo, typename CoeffType>
+   void Assemble(const Integrator &integ,
+                 CeedOperatorInfo &info,
+                 const mfem::FiniteElementSpace &fes,
+                 CoeffType *Q)
+   {
+      ElementsMap count;
+      ElementsMap element_indices;
+      ElementsMap offsets;
+
+      // Count the number of elements of each type
+      for (int i = 0; i < fes.GetNE(); i++)
+      {
+         ElementKey key(fes.GetElementType(i), fes.GetElementOrder(i));
+         auto value = count.find(key);
+         if (value == count.end())
+         {
+            count[key] = new int(1);
+         }
+         else
+         {
+            (*value->second)++;
+         }
+      }
+
+      // Initialization of the arrays
+      for ( const auto& value : count )
+      {
+         element_indices[value.first] = new int[*value.second];
+         offsets[value.first] = new int(0);
+      }
+
+      // Populates the indices arrays for each element type
+      for (int i = 0; i < fes.GetNE(); i++)
+      {
+         ElementKey key(fes.GetElementType(i), fes.GetElementOrder(i));
+         int &offset = *(offsets[key]);
+         int* indices_array = element_indices[key];
+         indices_array[offset] = i;
+         offset++;
+      }
+
+      // Create composite CeedOperator
+      CeedCompositeOperatorCreate(internal::ceed, &oper);
+
+      // Create each sub-CeedOperator
+      sub_ops.reserve(element_indices.size());
+      for (const auto& value : element_indices)
+      {
+         const int* indices = value.second;
+         const int first_index = indices[0];
+         const mfem::FiniteElement &el = *fes.GetFE(first_index);
+         auto &T = *fes.GetMesh()->GetElementTransformation(first_index);
+         MFEM_ASSERT(!integ.GetIntegrationRule(),
+                     "Mixed mesh integrators should not have an"
+                     " IntegrationRule.");
+         const IntegrationRule &ir = GetRule(integ, el, el, T);
+         auto sub_op = new CeedInteg();
+         int nelem = *count[value.first];
+         sub_op->Assemble(info, fes, ir, nelem, indices, Q);
+         sub_ops.push_back(sub_op);
+         CeedCompositeOperatorAddSub(oper, sub_op->GetCeedOperator());
+      }
+
+      const int ndofs = fes.GetVDim() * fes.GetNDofs();
+      CeedVectorCreate(internal::ceed, ndofs, &u);
+      CeedVectorCreate(internal::ceed, ndofs, &v);
+   }
+
+   virtual ~MixedIntegrator()
+   {
+      for (auto sub_op : sub_ops)
+      {
+         delete sub_op;
+      }
+   }
+#endif
+};
+
+} // namespace ceed
+
+} // namespace mfem
+
+#endif // MFEM_LIBCEED_MIXED_INTEGRATOR

fem/ceed/interface/restriction.cpp

@@ -20,8 +20,8 @@ namespace ceed
 
 #ifdef MFEM_USE_CEED
 
-static void InitNonTensorRestriction(const mfem::FiniteElementSpace &fes,
-                                     Ceed ceed, CeedElemRestriction *restr)
+static void InitNativeRestr(const mfem::FiniteElementSpace &fes,
+                            Ceed ceed, CeedElemRestriction *restr)
 {
    const mfem::FiniteElement *fe = fes.GetFE(0);
    const int P = fe->GetDof();
@@ -31,77 +31,173 @@ static void InitNonTensorRestriction(const mfem::FiniteElementSpace &fes,
    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 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++)
       {
-         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];
-         }
+         tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[dof_map[j]+el_offset];
       }
    }
+
    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)
+static void InitLexicoRestr(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 int stride = compstride == 1 ? fes.GetVDim() : 1;
+
+   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);
+}
+
+static void InitRestrictionImpl(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");
+   if ( tfe && tfe->GetDofMap().Size()>0 ) // Native ordering using dof_map
+   {
+      InitNativeRestr(fes, ceed, restr);
+   }
+   else  // Lexicographic ordering
+   {
+      InitLexicoRestr(fes, ceed, restr);
+   }
+}
+
+static void InitNativeRestrWithIndices(
+   const mfem::FiniteElementSpace &fes,
+   int nelem,
+   const int* indices,
+   Ceed ceed, CeedElemRestriction *restr)
+{
+   const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
+   const int P = fe->GetDof();
+   CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
+   mfem::Array<int> tp_el_dof(nelem*P);
+   const mfem::TensorBasisElement * tfe =
+      dynamic_cast<const mfem::TensorBasisElement *>(fe);
+   Array<int> dofs;
+   const int stride = compstride == 1 ? fes.GetVDim() : 1;
    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 < nelem; i++)
    {
-      for (int i = 0; i < fes.GetNE(); i++)
+      const int elem_index = indices[i];
+      fes.GetElementDofs(elem_index, dofs);
+      const int el_offset = P * i;
+      for (int j = 0; j < P; j++)
       {
-         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];
-         }
+         tp_el_dof[j + el_offset] = stride*dofs[dof_map[j]];
       }
    }
-   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(),
+
+   CeedElemRestrictionCreate(ceed, nelem, P, fes.GetVDim(),
                              compstride, (fes.GetVDim())*(fes.GetNDofs()),
                              CEED_MEM_HOST, CEED_COPY_VALUES,
                              tp_el_dof.GetData(), restr);
 }
 
+static void InitLexicoRestrWithIndices(
+   const mfem::FiniteElementSpace &fes,
+   int nelem,
+   const int* indices,
+   Ceed ceed, CeedElemRestriction *restr)
+{
+   const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
+   const int P = fe->GetDof();
+   CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
+   mfem::Array<int> tp_el_dof(nelem*P);
+   Array<int> dofs;
+   const int stride = compstride == 1 ? fes.GetVDim() : 1;
+
+   for (int i = 0; i < nelem; i++)
+   {
+      const int elem_index = indices[i];
+      fes.GetElementDofs(elem_index, dofs);
+      const int el_offset = P * i;
+      for (int j = 0; j < P; j++)
+      {
+         tp_el_dof[j + el_offset] = stride*dofs[j];
+      }
+   }
+
+   CeedElemRestrictionCreate(ceed, nelem, P, fes.GetVDim(),
+                             compstride, (fes.GetVDim())*(fes.GetNDofs()),
+                             CEED_MEM_HOST, CEED_COPY_VALUES,
+                             tp_el_dof.GetData(), restr);
+}
+
+static void InitRestrictionWithIndicesImpl(
+   const mfem::FiniteElementSpace &fes,
+   int nelem,
+   const int* indices,
+   Ceed ceed, CeedElemRestriction *restr)
+{
+   const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
+   const mfem::TensorBasisElement * tfe =
+      dynamic_cast<const mfem::TensorBasisElement *>(fe);
+   if ( tfe && tfe->GetDofMap().Size()>0 ) // Native ordering using dof_map
+   {
+      InitNativeRestrWithIndices(fes, nelem, indices, ceed, restr);
+   }
+   else  // Lexicographic ordering
+   {
+      InitLexicoRestrWithIndices(fes, nelem, indices, ceed, restr);
+   }
+}
+
+static void InitCoeffRestrictionWithIndicesImpl(
+   const mfem::FiniteElementSpace &fes,
+   int nelem,
+   const int* indices,
+   int nquads,
+   int ncomp,
+   Ceed ceed,
+   CeedElemRestriction *restr)
+{
+   mfem::Array<int> tp_el_dof(nelem*nquads);
+   const int stride_quad = ncomp;
+   const int stride_elem = ncomp*nquads;
+   // TODO generalize to support different #quads
+   for (int i = 0; i < nelem; i++)
+   {
+      const int elem_index = indices[i];
+      const int el_offset = elem_index * stride_elem;
+      for (int j = 0; j < nquads; j++)
+      {
+         tp_el_dof[j + nquads * i] = j * stride_quad + el_offset;
+      }
+   }
+   CeedElemRestrictionCreate(ceed, nelem, nquads, ncomp, 1,
+                             ncomp*fes.GetNE()*nquads,
+                             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,
@@ -139,14 +235,57 @@ void InitRestriction(const FiniteElementSpace &fes,
    // 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);
-      }
+      InitRestrictionImpl(fes, ceed, restr);
+      mfem::internal::ceed_restr_map[restr_key] = *restr;
+   }
+   else
+   {
+      *restr = restr_itr->second;
+   }
+}
+
+void InitRestrictionWithIndices(const FiniteElementSpace &fes,
+                                int nelem,
+                                const int* indices,
+                                Ceed ceed,
+                                CeedElemRestriction *restr)
+{
+   // Check for FES -> basis, restriction in hash tables
+   const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
+   const int P = fe->GetDof();
+   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())
+   {
+      InitRestrictionWithIndicesImpl(fes, nelem, indices, ceed, restr);
+      mfem::internal::ceed_restr_map[restr_key] = *restr;
+   }
+   else
+   {
+      *restr = restr_itr->second;
+   }
+}
+
+void InitCoeffRestrictionWithIndices(const FiniteElementSpace &fes,
+                                     int nelem,
+                                     const int* indices,
+                                     int nquads,
+                                     int ncomp,
+                                     Ceed ceed,
+                                     CeedElemRestriction *restr)
+{
+   // Check for FES -> basis, restriction in hash tables
+   RestrKey restr_key(&fes, nelem, nquads, ncomp, restr_type::Coeff);
+   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())
+   {
+      InitCoeffRestrictionWithIndicesImpl(fes, nelem, indices, nquads, ncomp,
+                                          ceed, restr);
       mfem::internal::ceed_restr_map[restr_key] = *restr;
    }
    else

fem/ceed/interface/restriction.hpp

@@ -21,37 +21,63 @@ 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 for non-mixed meshes.
 
-/** @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.
+    @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);
+/** @brief Initialize a CeedElemRestriction for mixed meshes.
+
+    @param[in] fes The finite element space.
+    @param[in] ceed The Ceed object.
+    @param[in] nelem The number of elements.
+    @param[in] indices The indices of the elements of same type in the
+                       `FiniteElementSpace`.
+    @param[out] restr The `CeedElemRestriction` to initialize. */
+void InitRestrictionWithIndices(const FiniteElementSpace &fes,
+                                int nelem,
+                                const int* indices,
+                                Ceed ceed,
+                                CeedElemRestriction *restr);
+
+/** @brief Initialize a strided CeedElemRestriction
+
+    @param[in] nelem is the number of elements.
+    @param[in] nqpts is the total number of quadrature points.
+    @param[in] qdatasize is the number of data per quadrature point.
+    @param[in] 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.
+    @param[out] restr The `CeedElemRestriction` to initialize. */
+void InitStridedRestriction(const mfem::FiniteElementSpace &fes,
+                            CeedInt nelem, CeedInt nqpts, CeedInt qdatasize,
+                            const CeedInt *strides,
+                            CeedElemRestriction *restr);
+
+/** @brief Initialize a CeedElemRestriction for a mfem::Coefficient on a mixed
+    mesh.
+
+    @param[in] fes The finite element space.
+    @param[in] nelem is the number of elements.
+    @param[in] indices The indices of the elements of same type in the
+                       `FiniteElementSpace`.
+    @param[in] nquads is the total number of quadrature points
+    @param[in] ncomp is the number of data per quadrature point
+    @param[in] ceed The Ceed object.
+    @param[out] restr The `CeedElemRestriction` to initialize. */
+void InitCoeffRestrictionWithIndices(const FiniteElementSpace &fes,
+                                     int nelem,
+                                     const int* indices,
+                                     int nquads,
+                                     int ncomp,
+                                     Ceed ceed,
+                                     CeedElemRestriction *restr);
 
 #endif
 

fem/ceed/interface/util.cpp

@@ -99,6 +99,34 @@ void InitBasisAndRestriction(const FiniteElementSpace &fes,
    InitRestriction(fes, ceed, restr);
 }
 
+void InitBasisAndRestrictionWithIndices(const FiniteElementSpace &fes,
+                                        const IntegrationRule &irm,
+                                        int nelem,
+                                        const int* indices,
+                                        Ceed ceed, CeedBasis *basis,
+                                        CeedElemRestriction *restr)
+{
+   InitBasisWithIndices(fes, irm, nelem, indices, ceed, basis);
+   InitRestrictionWithIndices(fes, nelem, indices, ceed, restr);
+}
+
+void InitBasisAndRestriction(const FiniteElementSpace &fes,
+                             const IntegrationRule &irm,
+                             int nelem,
+                             const int* indices,
+                             Ceed ceed, CeedBasis *basis,
+                             CeedElemRestriction *restr)
+{
+   if (indices)
+   {
+      InitBasisAndRestrictionWithIndices(fes,irm,nelem,indices,ceed,basis,restr);
+   }
+   else
+   {
+      InitBasisAndRestriction(fes,irm,ceed,basis,restr);
+   }
+}
+
 // Assumes a tensor-product operator with one active field
 int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
 {
@@ -158,6 +186,66 @@ int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
    return 0;
 }
 
+template <>
+const IntegrationRule & GetRule<MassIntegrator>(
+   const MassIntegrator &integ,
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &trans)
+{
+   return MassIntegrator::GetRule(trial_fe, test_fe, trans);
+}
+
+template <>
+const IntegrationRule & GetRule<VectorMassIntegrator>(
+   const VectorMassIntegrator &integ,
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &trans)
+{
+   return MassIntegrator::GetRule(trial_fe, test_fe, trans);
+}
+
+template <>
+const IntegrationRule & GetRule<ConvectionIntegrator>(
+   const ConvectionIntegrator &integ,
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &trans)
+{
+   return ConvectionIntegrator::GetRule(trial_fe, test_fe, trans);
+}
+
+template <>
+const IntegrationRule & GetRule<VectorConvectionNLFIntegrator>(
+   const VectorConvectionNLFIntegrator &integ,
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &trans)
+{
+   return VectorConvectionNLFIntegrator::GetRule(trial_fe, trans);
+}
+
+template <>
+const IntegrationRule & GetRule<DiffusionIntegrator>(
+   const DiffusionIntegrator &integ,
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &trans)
+{
+   return DiffusionIntegrator::GetRule(trial_fe, test_fe);
+}
+
+template <>
+const IntegrationRule & GetRule<VectorDiffusionIntegrator>(
+   const VectorDiffusionIntegrator &integ,
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &trans)
+{
+   return DiffusionIntegrator::GetRule(trial_fe, test_fe);
+}
+
 std::string ceed_path;
 
 const std::string &GetCeedPath()

fem/ceed/interface/util.hpp

@@ -26,7 +26,9 @@
 namespace mfem
 {
 
+class FiniteElement;
 class FiniteElementSpace;
+class ElementTransformation;
 class IntegrationRule;
 class Vector;
 
@@ -55,15 +57,51 @@ void RemoveBasisAndRestriction(const mfem::FiniteElementSpace *fes);
 /// Initialize a CeedVector from an mfem::Vector
 void InitVector(const mfem::Vector &v, CeedVector &cv);
 
-/** Initialize a CeedBasis and a CeedElemRestriction based on an
-    mfem::FiniteElementSpace @a fes, and an mfem::IntegrationRule @a ir. */
+/** @brief Initialize a CeedBasis and a CeedElemRestriction based on an
+    mfem::FiniteElementSpace @a fes, and an mfem::IntegrationRule @a ir.
+
+    @param[in] fes The finite element space.
+    @param[in] ir The integration rule.
+    @param[in] ceed The Ceed object.
+    @param[out] basis The `CeedBasis` to initialize.
+    @param[out] restr The `CeedElemRestriction` to initialize.
+
+    @warning Only for non-mixed finite element spaces. */
 void InitBasisAndRestriction(const mfem::FiniteElementSpace &fes,
                              const mfem::IntegrationRule &ir,
                              Ceed ceed, CeedBasis *basis,
                              CeedElemRestriction *restr);
 
+/** @brief Initialize a CeedBasis and a CeedElemRestriction based on an
+    mfem::FiniteElementSpace @a fes, and an mfem::IntegrationRule @a ir,
+    and a list of @a nelem elements of indices @a indices.
+
+    @param[in] fes The finite element space.
+    @param[in] ir The integration rule.
+    @param[in] nelem The number of elements.
+    @param[in] indices The indices of the elements of same type in the
+                       `FiniteElementSpace`. If `indices == nullptr`, assumes
+                       that the `FiniteElementSpace` is not mixed.
+    @param[in] ceed The Ceed object.
+    @param[out] basis The `CeedBasis` to initialize.
+    @param[out] restr The `CeedElemRestriction` to initialize. */
+void InitBasisAndRestriction(const FiniteElementSpace &fes,
+                             const IntegrationRule &ir,
+                             int nelem,
+                             const int* indices,
+                             Ceed ceed, CeedBasis *basis,
+                             CeedElemRestriction *restr);
+
 int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field);
 
+
+template <typename Integrator>
+const IntegrationRule & GetRule(
+   const Integrator &integ,
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans);
+
 /// Return the path to the libCEED q-function headers.
 const std::string &GetCeedPath();
 
@@ -87,7 +125,7 @@ struct BasisHash
 };
 using BasisMap = std::unordered_map<const BasisKey, CeedBasis, BasisHash>;
 
-enum restr_type {Standard, Strided};
+enum restr_type {Standard, Strided, Coeff};
 
 // Hash table for CeedElemRestriction
 using RestrKey =
@@ -117,6 +155,8 @@ namespace internal
 {
 
 #ifdef MFEM_USE_CEED
+/** @warning These maps have a tendency to create bugs when adding new "types"
+    of CeedBasis and CeedElemRestriction. */
 extern ceed::BasisMap ceed_basis_map;
 extern ceed::RestrMap ceed_restr_map;
 #endif

fem/ceed/solvers/algebraic.cpp

@@ -633,7 +633,7 @@ AlgebraicSpaceHierarchy::AlgebraicSpaceHierarchy(FiniteElementSpace &fes)
    current_order = order;
 
    Ceed ceed = internal::ceed;
-   InitTensorRestriction(fes, ceed, &fine_er);
+   InitRestriction(fes, ceed, &fine_er);
    CeedElemRestriction er = fine_er;
 
    int dim = fes.GetMesh()->Dimension();

fem/ceed/solvers/algebraic.hpp

@@ -134,7 +134,6 @@ public:
    }
    ~AlgebraicSpaceHierarchy()
    {
-      CeedElemRestrictionDestroy(&fine_er);
       for (int i=0; i<R_tr.Size(); ++i)
       {
          delete R_tr[i];

fem/datacollection.cpp

@@ -1078,7 +1078,7 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
          it->second->GetValues(i, RefG->RefPts, val, pmat);
          for (int j = 0; j < val.Size(); j++)
          {
-            WriteBinaryOrASCII(out, buf, val(j), "\n", pv_data_format);
+            WriteBinaryOrASCII(os, buf, val(j), "\n", pv_data_format);
          }
       }
    }
@@ -1094,7 +1094,7 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
          {
             for (int ii = 0; ii < vval.Height(); ii++)
             {
-               WriteBinaryOrASCII(out, buf, vval(ii,jj), " ", pv_data_format);
+               WriteBinaryOrASCII(os, buf, vval(ii,jj), " ", pv_data_format);
             }
             if (pv_data_format == VTKFormat::ASCII) { os << '\n'; }
          }

fem/fespace.hpp

@@ -980,14 +980,14 @@ public:
    void Save(std::ostream &out) const;
 };
 
+/// @brief Return true if the mesh contains only one topology and the elements are tensor elements.
 inline bool UsesTensorBasis(const FiniteElementSpace& fes)
 {
-   // TODO: mixed meshes: return true if there is at least one tensor-product
-   // Geometry in the global mesh and the FE collection returns a
-   // TensorBasisElement for that Geometry?
-
+   Mesh & mesh = *fes.GetMesh();
+   const bool mixed = mesh.GetNumGeometries(mesh.Dimension()) > 1;
    // Potential issue: empty local mesh --> no element 0.
-   return dynamic_cast<const mfem::TensorBasisElement *>(fes.GetFE(0))!=nullptr;
+   return !mixed &&
+          dynamic_cast<const mfem::TensorBasisElement *>(fes.GetFE(0))!=nullptr;
 }
 
 }

fem/gridfunc.cpp

@@ -4276,7 +4276,7 @@ void TensorProductLegendre(int dim,                // input
       poly1d.CalcLegendre(order, x3, poly_z);
    }
 
-   int basis_dimension = pow(order+1,dim);
+   int basis_dimension = static_cast<int>(pow(order+1,dim));
    poly.SetSize(basis_dimension);
    switch (dim)
    {
@@ -4458,7 +4458,7 @@ double LSZZErrorEstimator(BilinearFormIntegrator &blfi,  // input
       const int patch_order = max(ufes->GetElementOrder(el1),
                                   ufes->GetElementOrder(el2));
 
-      int num_basis_functions = pow(patch_order+1,dim);
+      int num_basis_functions = static_cast<int>(pow(patch_order+1,dim));
       int flux_order = 2*patch_order + 1;
       DenseMatrix A(num_basis_functions);
       Array<double> b(sdim * num_basis_functions);

fem/linearform.cpp

@@ -129,12 +129,7 @@ bool LinearForm::SupportsDevice()
    if (mesh_dim == 1 || mesh_dim != mesh.SpaceDimension()) { return false; }
 
    // tensor-product finite element space only
-   // with point values preserving scalar fields
-   for (int e = 0; e < fes->GetNE(); ++e)
-   {
-      const FiniteElement *fe = fes->GetFE(e);
-      if (!dynamic_cast<const TensorBasisElement*>(fe)) { return false; }
-   }
+   if (!UsesTensorBasis(*fes)) { return false; }
 
    return true;
 }

fem/lor/lor.cpp

@@ -467,7 +467,8 @@ void LORDiscretization::FormLORSpace()
    mesh = new Mesh(Mesh::MakeRefined(mesh_ho, refinements, ref_type));
 
    fec = fes_ho.FEColl()->Clone(GetLOROrder());
-   fes = new FiniteElementSpace(mesh, fec);
+   const int vdim = fes_ho.GetVDim();
+   fes = new FiniteElementSpace(mesh, fec, vdim);
    SetupProlongationAndRestriction();
 }
 
@@ -511,7 +512,8 @@ void ParLORDiscretization::FormLORSpace()
    mesh = pmesh;
 
    fec = pfes_ho.FEColl()->Clone(GetLOROrder());
-   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec);
+   const int vdim = fes_ho.GetVDim();
+   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec, vdim);
    fes = pfes;
    SetupProlongationAndRestriction();
 }

fem/lor/lor_ams.cpp

@@ -39,7 +39,7 @@ void BatchedLOR_AMS::Form2DEdgeToVertex_ND(Array<int> &edge2vert)
 {
    const int o = order;
    const int op1 = o + 1;
-   const int nedge = dim*o*pow(op1, dim-1);
+   const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
 
    edge2vert.SetSize(2*nedge);
    auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
@@ -73,7 +73,7 @@ void BatchedLOR_AMS::Form2DEdgeToVertex_RT(Array<int> &edge2vert)
 {
    const int o = order;
    const int op1 = o + 1;
-   const int nedge = dim*o*pow(op1, dim-1);
+   const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
 
    edge2vert.SetSize(2*nedge);
    auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
@@ -106,7 +106,7 @@ void BatchedLOR_AMS::Form3DEdgeToVertex(Array<int> &edge2vert)
 {
    const int o = order;
    const int op1 = o + 1;
-   const int nedge = dim*o*pow(op1, dim-1);
+   const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
 
    edge2vert.SetSize(2*nedge);
    auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
@@ -170,8 +170,8 @@ void BatchedLOR_AMS::FormGradientMatrix()
    MFEM_VERIFY(R_v != NULL && R_e != NULL, "");
 
    const int nel_ho = edge_fes.GetNE();
-   const int nedge_per_el = dim*order*pow(order + 1, dim - 1);
-   const int nvert_per_el = pow(order + 1, dim);
+   const int nedge_per_el = static_cast<int>(dim*order*pow(order + 1, dim - 1));
+   const int nvert_per_el = static_cast<int>(pow(order + 1, dim));
 
    const auto offsets_e = R_e->Offsets().Read();
    const auto indices_e = R_e->Indices().Read();
@@ -274,7 +274,7 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
 
    const int nel_ho = vert_fes.GetNE();
    const int ndp1 = order + 1;
-   const int ndof_per_el = pow(ndp1, dim);
+   const int ndof_per_el = static_cast<int>(pow(ndp1, dim));
    const int sdim = dim;
    const int ntdofs = R->Height();
 

fem/lor/lor_batched.cpp

@@ -80,7 +80,7 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
    const int nel_ho = mesh_ho.GetNE();
    const int order = fes_ho.GetMaxElementOrder();
    const int nd1d = order + 1;
-   const int ndof_per_el = pow(nd1d, dim);
+   const int ndof_per_el = static_cast<int>(pow(nd1d, dim));
 
    const GridFunction *nodal_gf = mesh_ho.GetNodes();
    const FiniteElementSpace *nodal_fes = nodal_gf->FESpace();
@@ -284,7 +284,7 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
          const bool plus = si_E >= 0;
          const int i_E = plus ? si_E : -1 - si_E;
          i_elts[e_i] = i_E/ndof_per_el;
-         const double i_Bi = i_E%ndof_per_el;
+         const int i_Bi = i_E % ndof_per_el;
          i_B[e_i] = plus ? i_Bi : -1 - i_Bi; // encode with sign
       }
       for (int j=0; j<nnz_per_row; ++j)
@@ -314,7 +314,7 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
                const bool plus = sj_E >= 0;
                const int j_E = plus ? sj_E : -1 - sj_E;
                j_elts[e_j] = j_E/ndof_per_el;
-               const double j_Bj = j_E%ndof_per_el;
+               const int j_Bj = j_E % ndof_per_el;
                j_B[e_j] = plus ? j_Bj : -1 - j_Bj; // encode with sign
             }
             const int min_e = GetMinElt(i_elts, i_ne, j_elts, j_ne);

fem/nonlinearform_ext.cpp

@@ -25,12 +25,16 @@ PANonlinearFormExtension::PANonlinearFormExtension(const NonlinearForm *nlf):
    NonlinearFormExtension(nlf),
    fes(*nlf->FESpace()),
    dnfi(*nlf->GetDNFI()),
-   elemR(fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC)),
+   elemR(nullptr),
    Grad(*this)
 {
-   // TODO: optimize for the case when 'elemR' is identity
-   xe.SetSize(elemR->Height(), Device::GetMemoryType());
-   ye.SetSize(elemR->Height(), Device::GetMemoryType());
+   if (!DeviceCanUseCeed())
+   {
+      elemR = fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      // TODO: optimize for the case when 'elemR' is identity
+      xe.SetSize(elemR->Height(), Device::GetMemoryType());
+      ye.SetSize(elemR->Height(), Device::GetMemoryType());
+   }
    ye.UseDevice(true);
 }
 
@@ -135,13 +139,16 @@ void PANonlinearFormExtension::Gradient::Update()
 MFNonlinearFormExtension::MFNonlinearFormExtension(const NonlinearForm *form):
    NonlinearFormExtension(form), fes(*form->FESpace())
 {
-   const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
-   elem_restrict_lex = fes.GetElementRestriction(ordering);
-   if (elem_restrict_lex) // replace with a check for not identity
+   if (!DeviceCanUseCeed())
    {
-      localX.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
-      localY.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
-      localY.UseDevice(true); // ensure 'localY = 0.0' is done on device
+      const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
+      elem_restrict_lex = fes.GetElementRestriction(ordering);
+      if (elem_restrict_lex) // replace with a check for not identity
+      {
+         localX.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
+         localY.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
+         localY.UseDevice(true); // ensure 'localY = 0.0' is done on device
+      }
    }
 }
 

fem/nonlininteg_vectorconvection.cpp

@@ -28,7 +28,16 @@ void VectorConvectionNLFIntegrator::AssemblePA(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::PAVectorConvectionNLFIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedPAVectorConvectionNLIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::PAVectorConvectionNLFIntegrator(fes, *ir, Q);
+      }
       return;
    }
    dim = mesh->Dimension();

fem/nonlininteg_vectorconvection_mf.cpp

@@ -28,7 +28,16 @@ void VectorConvectionNLFIntegrator::AssembleMF(const FiniteElementSpace &fes)
    if (DeviceCanUseCeed())
    {
       delete ceedOp;
-      ceedOp = new ceed::MFVectorConvectionNLFIntegrator(fes, *ir, Q);
+      const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+                         fes.IsVariableOrder();
+      if (mixed)
+      {
+         ceedOp = new ceed::MixedMFVectorConvectionNLIntegrator(*this, fes, Q);
+      }
+      else
+      {
+         ceedOp = new ceed::MFVectorConvectionNLFIntegrator(fes, *ir, Q);
+      }
       return;
    }
    MFEM_ABORT("Not yet implemented.");

fem/pfespace.cpp

@@ -194,12 +194,12 @@ void ParFiniteElementSpace::Construct()
 
 void ParFiniteElementSpace::PrintPartitionStats()
 {
-   long ltdofs = ltdof_size;
-   long min_ltdofs, max_ltdofs, sum_ltdofs;
+   long long ltdofs = ltdof_size;
+   long long min_ltdofs, max_ltdofs, sum_ltdofs;
 
-   MPI_Reduce(&ltdofs, &min_ltdofs, 1, MPI_LONG, MPI_MIN, 0, MyComm);
-   MPI_Reduce(&ltdofs, &max_ltdofs, 1, MPI_LONG, MPI_MAX, 0, MyComm);
-   MPI_Reduce(&ltdofs, &sum_ltdofs, 1, MPI_LONG, MPI_SUM, 0, MyComm);
+   MPI_Reduce(&ltdofs, &min_ltdofs, 1, MPI_LONG_LONG, MPI_MIN, 0, MyComm);
+   MPI_Reduce(&ltdofs, &max_ltdofs, 1, MPI_LONG_LONG, MPI_MAX, 0, MyComm);
+   MPI_Reduce(&ltdofs, &sum_ltdofs, 1, MPI_LONG_LONG, MPI_SUM, 0, MyComm);
 
    if (MyRank == 0)
    {
@@ -219,14 +219,14 @@ void ParFiniteElementSpace::PrintPartitionStats()
          for (int i = 1; i < NRanks; i++)
          {
             MPI_Status status;
-            MPI_Recv(&ltdofs, 1, MPI_LONG, i, 123, MyComm, &status);
+            MPI_Recv(&ltdofs, 1, MPI_LONG_LONG, i, 123, MyComm, &status);
             mfem::out << " " << ltdofs;
          }
          mfem::out << "\n";
       }
       else
       {
-         MPI_Send(&ltdofs, 1, MPI_LONG, 0, 123, MyComm);
+         MPI_Send(&ltdofs, 1, MPI_LONG_LONG, 0, 123, MyComm);
       }
    }
 }
@@ -959,6 +959,10 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
    SparseMatrix Pdiag;
    P->GetDiag(Pdiag);
    R = Transpose(Pdiag);
+
+   // The following call ensures that the action of the transpose of P is
+   // performed fast when HYPRE is built for GPUs.
+   P->EnsureMultTranspose();
 }
 
 HypreParMatrix *ParFiniteElementSpace::GetPartialConformingInterpolation()
@@ -2624,6 +2628,10 @@ int ParFiniteElementSpace
    {
       *P_ = MakeVDimHypreMatrix(pmatrix, ndofs, num_true_dofs,
                                 dof_offs, tdof_offs);
+
+      // The following call ensures that the action of the transpose of *P_ is
+      // performed fast when HYPRE is built for GPUs.
+      (*P_)->EnsureMultTranspose();
    }
 
    // clean up possible remaining messages in the queue to avoid receiving

fem/tbilininteg.hpp

@@ -361,8 +361,8 @@ struct TDiffusionKernel<2,2,complex_t>
        @param F Jt [M x Dim x SDim x NE] - Jacobian transposed, data member in F
        @param Q CoefficientEval<>::Type
        @param q CoefficientEval<>::Type::result_t
-       @param A [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm. matrices
-       @param A [M x Dim x Dim]       partially assembled Dim x Dim matrices
+       @param A either [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm.
+       matrices, or [M x Dim x Dim] partially assembled Dim x Dim matrices.
    */
    template <typename T_result_t, typename Q_t, typename q_t, typename asm_type>
    static inline MFEM_ALWAYS_INLINE

general/array.hpp

@@ -298,7 +298,7 @@ public:
    inline const T* end() const { return data + size; }
 
    /// Returns the number of bytes allocated for the array including any reserve.
-   long MemoryUsage() const { return Capacity() * sizeof(T); }
+   std::size_t MemoryUsage() const { return Capacity() * sizeof(T); }
 
    /// Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
    const T *Read(bool on_dev = true) const
@@ -509,7 +509,7 @@ public:
 
    void Swap(BlockArray<T> &other);
 
-   long MemoryUsage() const;
+   std::size_t MemoryUsage() const;
 
 protected:
    template <typename cA, typename cT>
@@ -1043,7 +1043,7 @@ void BlockArray<T>::Swap(BlockArray<T> &other)
 }
 
 template<typename T>
-long BlockArray<T>::MemoryUsage() const
+std::size_t BlockArray<T>::MemoryUsage() const
 {
    return (mask+1)*sizeof(T)*blocks.Size() + blocks.MemoryUsage();
 }

general/hash.hpp

@@ -284,7 +284,7 @@ public:
    void Reparent(int id, int new_p1, int new_p2, int new_p3, int new_p4 = -1);
 
    /// @brief Return total size of allocated memory (tables plus items), in bytes.
-   long MemoryUsage() const;
+   std::size_t MemoryUsage() const;
 
    /// @brief Write details of the memory usage to the mfem output stream.
    void PrintMemoryDetail() const;
@@ -875,7 +875,7 @@ void HashTable<T>::Reparent(int id,
 }
 
 template<typename T>
-long HashTable<T>::MemoryUsage() const
+std::size_t HashTable<T>::MemoryUsage() const
 {
    return (mask+1) * sizeof(int) + Base::MemoryUsage() + unused.MemoryUsage();
 }

general/table.cpp

@@ -398,7 +398,7 @@ void Table::Swap(Table & other)
    mfem::Swap(J, other.J);
 }
 
-long Table::MemoryUsage() const
+std::size_t Table::MemoryUsage() const
 {
    if (size < 0 || I == NULL) { return 0; }
    return (size+1 + I[size]) * sizeof(int);

general/table.hpp

@@ -191,7 +191,7 @@ public:
 
    void Clear();
 
-   long MemoryUsage() const;
+   std::size_t MemoryUsage() const;
 
    /// Destroys Table.
    ~Table();

linalg/densemat.hpp

@@ -379,7 +379,7 @@ public:
    /// Invert and print the numerical conditioning of the inversion.
    void TestInversion();
 
-   long MemoryUsage() const { return data.Capacity() * sizeof(double); }
+   std::size_t MemoryUsage() const { return data.Capacity() * sizeof(double); }
 
    /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
    const double *Read(bool on_dev = true) const
@@ -1004,7 +1004,7 @@ public:
    void Clear()
    { UseExternalData(NULL, 0, 0, 0); }
 
-   long MemoryUsage() const { return nk*Mk.MemoryUsage(); }
+   std::size_t MemoryUsage() const { return nk*Mk.MemoryUsage(); }
 
    /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
    const double *Read(bool on_dev = true) const

linalg/dtensor.hpp

@@ -28,8 +28,8 @@ public:
 #if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
       MFEM_ASSERT(first<sizes[N-1],"Trying to access out of boundary.");
 #endif
-      return first + sizes[N - 1] * TensorInd < N + 1, Dim, Args... >
-             ::result(sizes, args...);
+      return static_cast<int>(first + sizes[N - 1] * TensorInd < N + 1, Dim, Args... >
+                              ::result(sizes, args...));
    }
 };
 
@@ -44,7 +44,7 @@ public:
 #if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
       MFEM_ASSERT(first<sizes[Dim-1],"Trying to access out of boundary.");
 #endif
-      return first;
+      return static_cast<int>(first);
    }
 };
 

linalg/hypre.cpp

@@ -1711,6 +1711,16 @@ HypreParMatrix *HypreParMatrix::ExtractSubmatrix(const Array<int> &indices,
 }
 #endif
 
+void HypreParMatrix::EnsureMultTranspose() const
+{
+#if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
+    (MFEM_HYPRE_VERSION > 22500)
+#ifdef HYPRE_USING_GPU
+   hypre_ParCSRMatrixLocalTranspose(A);
+#endif
+#endif
+}
+
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                                double a, double b) const
 {
@@ -1833,6 +1843,14 @@ void HypreParMatrix::MultTranspose(double a, const Vector &x,
       }
    }
 
+#if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
+    (MFEM_HYPRE_VERSION > 22500)
+#ifdef HYPRE_USING_GPU
+   MFEM_VERIFY(A->diagT != NULL,
+               "Transpose action requires EnsureMultTranspose()");
+#endif
+#endif
+
    hypre_ParCSRMatrixMatvecT(a, A, *Y, b, *X);
 
    if (!yshallow) { y = *X; }  // Deep copy
@@ -1848,6 +1866,13 @@ HYPRE_Int HypreParMatrix::Mult(HYPRE_ParVector x, HYPRE_ParVector y,
 HYPRE_Int HypreParMatrix::MultTranspose(HypreParVector & x, HypreParVector & y,
                                         double a, double b) const
 {
+#if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
+    (MFEM_HYPRE_VERSION > 22500)
+#ifdef HYPRE_USING_GPU
+   MFEM_VERIFY(A->diagT != NULL,
+               "Transpose action requires EnsureMultTranspose()");
+#endif
+#endif
    x.HypreRead();
    (b == 0.0) ? y.HypreWrite() : y.HypreReadWrite();
    return hypre_ParCSRMatrixMatvecT(a, A, x, b, y);

linalg/hypre.hpp

@@ -658,6 +658,12 @@ public:
 
    virtual MemoryClass GetMemoryClass() const { return GetHypreMemoryClass(); }
 
+   /// Ensure the action of the transpose is performed fast.
+   /** When HYPRE is built for GPUs, this method will construct and store the
+       transposes of the 'diag' and 'offd' CSR matrices. When HYPRE is not built
+       for GPUs, this method is a no-op. */
+   void EnsureMultTranspose() const;
+
    /// Computes y = alpha * A * x + beta * y
    HYPRE_Int Mult(HypreParVector &x, HypreParVector &y,
                   double alpha = 1.0, double beta = 0.0) const;

linalg/hypre_parcsr.cpp

@@ -1560,7 +1560,9 @@ void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
 
    if (num_cols_offd)
    {
-#if MFEM_HYPRE_VERSION >= 21100
+      // Disable the use of offdT for now, until we implement
+      // hypre_CSRMatrixAbsMatvec on device.
+#if MFEM_HYPRE_VERSION >= 21100 && 0
       if (A->offdT)
       {
          // offdT is optional. Used only if it's present.
@@ -1575,7 +1577,9 @@ void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
 
    comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg, y_tmp, y_buf);
 
-#if MFEM_HYPRE_VERSION >= 21100
+   // Disable the use of diagT for now, until we implement
+   // hypre_CSRMatrixAbsMatvec on device.
+#if MFEM_HYPRE_VERSION >= 21100 && 0
    if (A->diagT)
    {
       // diagT is optional. Used only if it's present.
@@ -1697,7 +1701,9 @@ void hypre_ParCSRMatrixBooleanMatvecT(hypre_ParCSRMatrix *A,
 
    if (num_cols_offd)
    {
-#if MFEM_HYPRE_VERSION >= 21100
+      // Disable the use of offdT for now, until we implement
+      // hypre_CSRMatrixBooleanMatvec on device.
+#if MFEM_HYPRE_VERSION >= 21100 && 0
       if (A->offdT)
       {
          // offdT is optional. Used only if it's present.
@@ -1712,7 +1718,9 @@ void hypre_ParCSRMatrixBooleanMatvecT(hypre_ParCSRMatrix *A,
 
    comm_handle = hypre_ParCSRCommHandleCreate_bool(2, comm_pkg, y_tmp, y_buf);
 
-#if MFEM_HYPRE_VERSION >= 21100
+   // Disable the use of diagT for now, until we implement
+   // hypre_CSRMatrixBooleanMatvec on device.
+#if MFEM_HYPRE_VERSION >= 21100 && 0
    if (A->diagT)
    {
       // diagT is optional. Used only if it's present.

linalg/ode.cpp

@@ -188,7 +188,6 @@ void ExplicitRKSolver::Step(Vector &x, double &t, double &dt)
 
 ExplicitRKSolver::~ExplicitRKSolver()
 {
-   delete [] k;
 }
 
 const double RK6Solver::a[] =

linalg/ode.hpp

@@ -331,7 +331,6 @@ public:
    ~AdamsMoultonSolver()
    {
       if (RKsolver) { delete RKsolver; }
-      delete [] k;
    };
 };
 

linalg/operator.cpp

@@ -201,7 +201,7 @@ void Operator::PrintMatlab(std::ostream & os, int n, int m) const
       Mult(x, y);
       for (int j = 0; j < m; j++)
       {
-         if (y(j))
+         if (y(j) != 0)
          {
             os << j+1 << " " << i+1 << " " << y(j) << '\n';
          }

linalg/solvers.cpp

@@ -528,7 +528,7 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
 
    // Optimized preconditioned SLI with fixed number of iterations and given
    // initial guess
-   if (!rel_tol && iterative_mode && prec)
+   if (rel_tol == 0.0 && iterative_mode && prec)
    {
       for (i = 0; i < max_iter; i++)
       {
@@ -544,7 +544,7 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
 
    // Optimized preconditioned SLI with fixed number of iterations and zero
    // initial guess
-   if (!rel_tol && !iterative_mode && prec)
+   if (rel_tol == 0.0 && !iterative_mode && prec)
    {
       prec->Mult(b, x);     // x = B b (initial guess 0)
       for (i = 1; i < max_iter; i++)
@@ -817,7 +817,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
       if (print_options.iterations)
       {
          mfem::out << "   Iteration : " << setw(3) << i << "  (B r, r) = "
-                   << betanom << '\n';
+                   << betanom << std::endl;
       }
 
       Monitor(i, betanom, r, x);
@@ -861,7 +861,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
       }
       nom = betanom;
    }
-   if (print_options.first_and_last)
+   if (print_options.first_and_last && !print_options.iterations)
    {
       mfem::out << "   Iteration : " << setw(3) << final_iter << "  (B r, r) = "
                 << betanom << '\n';
@@ -1597,6 +1597,13 @@ void MINRESSolver::SetOperator(const Operator &op)
    {
       u1.SetSize(width);
    }
+
+   v0.UseDevice(true);
+   v1.UseDevice(true);
+   w0.UseDevice(true);
+   w1.UseDevice(true);
+   q.UseDevice(true);
+   u1.UseDevice(true);
 }
 
 void MINRESSolver::Mult(const Vector &b, Vector &x) const
@@ -1606,6 +1613,9 @@ void MINRESSolver::Mult(const Vector &b, Vector &x) const
    // by Henk A. van der Vorst, 2003.
    // Extended to support an SPD preconditioner.
 
+   b.UseDevice(true);
+   x.UseDevice(true);
+
    int it;
    double beta, eta, gamma0, gamma1, sigma0, sigma1;
    double alpha, delta, rho1, rho2, rho3, norm_goal;

linalg/symmat.hpp

@@ -99,7 +99,7 @@ public:
 
    DenseSymmetricMatrix &operator*=(double c);
 
-   long MemoryUsage() const { return data.Capacity() * sizeof(double); }
+   std::size_t MemoryUsage() const { return data.Capacity() * sizeof(double); }
 
    /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
    const double *Read(bool on_dev = true) const

mesh/mesh.hpp

@@ -954,10 +954,10 @@ public:
    virtual int GetNFbyType(FaceType type) const;
 
    /// Utility function: sum integers from all processors (Allreduce).
-   virtual long ReduceInt(int value) const { return value; }
+   virtual long long ReduceInt(int value) const { return value; }
 
    /// Return the total (global) number of elements.
-   long GetGlobalNE() const { return ReduceInt(NumOfElements); }
+   long long GetGlobalNE() const { return ReduceInt(NumOfElements); }
 
    /** @brief Return the mesh geometric factors corresponding to the given
        integration rule.
@@ -1531,7 +1531,7 @@ public:
    /** Replace the internal node GridFunction with a new GridFunction defined
        on the given FiniteElementSpace. The new node coordinates are projected
        (derived) from the current nodes/vertices. */
-   void SetNodalFESpace(FiniteElementSpace *nfes);
+   virtual void SetNodalFESpace(FiniteElementSpace *nfes);
    /** Replace the internal node GridFunction with the given GridFunction. The
        given GridFunction is updated with node coordinates projected (derived)
        from the current nodes/vertices. */

mesh/mesh_operators.cpp

@@ -59,10 +59,10 @@ ThresholdRefiner::ThresholdRefiner(ErrorEstimator &est)
    total_err_goal = 0.0;
    total_fraction = 0.5;
    local_err_goal = 0.0;
-   max_elements = std::numeric_limits<long>::max();
+   max_elements = std::numeric_limits<long long>::max();
 
    threshold = 0.0;
-   num_marked_elements = 0L;
+   num_marked_elements = 0LL;
    current_sequence = -1;
 
    non_conforming = -1;
@@ -84,11 +84,11 @@ double ThresholdRefiner::GetNorm(const Vector &local_err, Mesh &mesh) const
 int ThresholdRefiner::ApplyImpl(Mesh &mesh)
 {
    threshold = 0.0;
-   num_marked_elements = 0;
+   num_marked_elements = 0LL;
    marked_elements.SetSize(0);
    current_sequence = mesh.GetSequence();
 
-   const long num_elements = mesh.GetGlobalNE();
+   const long long num_elements = mesh.GetGlobalNE();
    if (num_elements >= max_elements) { return STOP; }
 
    const int NE = mesh.GetNE();
@@ -131,7 +131,7 @@ int ThresholdRefiner::ApplyImpl(Mesh &mesh)
    }
 
    num_marked_elements = mesh.ReduceInt(marked_elements.Size());
-   if (num_marked_elements == 0) { return STOP; }
+   if (num_marked_elements == 0LL) { return STOP; }
 
    mesh.GeneralRefinement(marked_elements, non_conforming, nc_limit);
    return CONTINUE + REFINED;
@@ -141,7 +141,7 @@ void ThresholdRefiner::Reset()
 {
    estimator.Reset();
    current_sequence = -1;
-   num_marked_elements = 0;
+   num_marked_elements = 0LL;
    // marked_elements.SetSize(0); // not necessary
 }
 

mesh/mesh_operators.hpp

@@ -182,10 +182,10 @@ protected:
    double total_err_goal;
    double total_fraction;
    double local_err_goal;
-   long   max_elements;
+   long long max_elements;
 
    double threshold;
-   long num_marked_elements;
+   long long num_marked_elements;
 
    Array<Refinement> marked_elements;
    long current_sequence;
@@ -230,7 +230,7 @@ public:
    /** @brief Set the maximum number of elements stopping criterion: stop when
        the input mesh has num_elements >= max_elem. The default value is
        LONG_MAX. */
-   void SetMaxElements(long max_elem) { max_elements = max_elem; }
+   void SetMaxElements(long long max_elem) { max_elements = max_elem; }
 
    /// Use nonconforming refinement, if possible (triangles, quads, hexes).
    void PreferNonconformingRefinement() { non_conforming = 1; }
@@ -248,7 +248,7 @@ public:
    }
 
    /// Get the number of marked elements in the last Apply() call.
-   long GetNumMarkedElements() const { return num_marked_elements; }
+   long long GetNumMarkedElements() const { return num_marked_elements; }
 
    /// Get the threshold used in the last Apply() call.
    double GetThreshold() const { return threshold; }
@@ -335,7 +335,7 @@ protected:
    int nc_limit = 1;
    int nonconforming = -1;
    int order;
-   long max_elements = std::numeric_limits<long>::max();
+   long long max_elements = std::numeric_limits<long long>::max();
    double threshold = 1.0e-2;
    double global_osc = NAN;
    Array<int> mesh_refinements;
@@ -378,7 +378,7 @@ public:
    /** @brief Set the maximum number of elements stopping criterion: stop when
        the input mesh has num_elements >= max_elem. The default value is
        LONG_MAX. */
-   void SetMaxElements(long max_elements_) { max_elements = max_elements_; }
+   void SetMaxElements(long long max_elements_) { max_elements = max_elements_; }
 
    /// Reset the function f
    void ResetCoefficient(Coefficient &coeff_)

mesh/ncmesh.cpp

@@ -5984,9 +5984,9 @@ void NCMesh::Trim()
    // maybe also of 'nodes' and 'faces'.
 }
 
-long NCMesh::NCList::MemoryUsage() const
+std::size_t NCMesh::NCList::MemoryUsage() const
 {
-   int pm_size = 0;
+   std::size_t pm_size = 0;
    for (int i = 0; i < Geometry::NumGeom; i++)
    {
       for (int j = 0; j < point_matrices[i].Size(); i++)
@@ -6002,9 +6002,9 @@ long NCMesh::NCList::MemoryUsage() const
           pm_size;
 }
 
-long CoarseFineTransformations::MemoryUsage() const
+std::size_t CoarseFineTransformations::MemoryUsage() const
 {
-   long mem = embeddings.MemoryUsage();
+   std::size_t mem = embeddings.MemoryUsage();
    for (int i = 0; i < Geometry::NumGeom; i++)
    {
       mem += point_matrices[i].MemoryUsage();
@@ -6012,7 +6012,7 @@ long CoarseFineTransformations::MemoryUsage() const
    return mem;
 }
 
-long NCMesh::MemoryUsage() const
+std::size_t NCMesh::MemoryUsage() const
 {
    return nodes.MemoryUsage() +
           faces.MemoryUsage() +

mesh/ncmesh.hpp

@@ -84,7 +84,7 @@ struct CoarseFineTransformations
 
    void Clear();
    bool IsInitialized() const;
-   long MemoryUsage() const;
+   std::size_t MemoryUsage() const;
 
    MFEM_DEPRECATED
    void GetCoarseToFineMap(const Mesh &fine_mesh, Table &coarse_to_fine) const
@@ -240,7 +240,7 @@ public:
       void Clear();
       bool Empty() const { return !conforming.Size() && !masters.Size(); }
       long TotalSize() const;
-      long MemoryUsage() const;
+      std::size_t MemoryUsage() const;
 
       const MeshId& LookUp(int index, int *type = NULL) const;
 
@@ -390,7 +390,7 @@ public:
    virtual void Trim();
 
    /// Return total number of bytes allocated.
-   long MemoryUsage() const;
+   std::size_t MemoryUsage() const;
 
    int PrintMemoryDetail() const;
 

mesh/pmesh.cpp

@@ -869,8 +869,9 @@ void ParMesh::ComputeGlobalElementOffset() const
 {
    if (glob_offset_sequence != sequence) // mesh has changed
    {
-      long local_elems = NumOfElements;
-      MPI_Scan(&local_elems, &glob_elem_offset, 1, MPI_LONG, MPI_SUM, MyComm);
+      long long local_elems = NumOfElements;
+      MPI_Scan(&local_elems, &glob_elem_offset, 1, MPI_LONG_LONG, MPI_SUM,
+               MyComm);
       glob_elem_offset -= local_elems;
 
       glob_offset_sequence = sequence; // don't recalculate until refinement etc.
@@ -1526,15 +1527,15 @@ void ParMesh::Finalize(bool refine, bool fix_orientation)
    FinalizeParTopo();
 }
 
-int ParMesh::GetLocalElementNum(long global_element_num) const
+int ParMesh::GetLocalElementNum(long long global_element_num) const
 {
    ComputeGlobalElementOffset();
-   long local = global_element_num - glob_elem_offset;
+   long long local = global_element_num - glob_elem_offset;
    if (local < 0 || local >= NumOfElements) { return -1; }
    return local;
 }
 
-long ParMesh::GetGlobalElementNum(int local_element_num) const
+long long ParMesh::GetGlobalElementNum(int local_element_num) const
 {
    ComputeGlobalElementOffset();
    return glob_elem_offset + local_element_num;
@@ -1990,6 +1991,25 @@ void ParMesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
    Nodes->MakeOwner(nfec);
 }
 
+void ParMesh::SetNodalFESpace(FiniteElementSpace *nfes)
+{
+   ParFiniteElementSpace *npfes = dynamic_cast<ParFiniteElementSpace*>(nfes);
+   if (npfes)
+   {
+      SetNodalFESpace(npfes);
+   }
+   else
+   {
+      Mesh::SetNodalFESpace(nfes);
+   }
+}
+
+void ParMesh::SetNodalFESpace(ParFiniteElementSpace *npfes)
+{
+   ParGridFunction *nodes = new ParGridFunction(npfes);
+   SetNodalGridFunction(nodes, true);
+}
+
 void ParMesh::EnsureParNodes()
 {
    if (Nodes && dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace()) == NULL)
@@ -3870,7 +3890,7 @@ bool ParMesh::NonconformingDerefinement(Array<double> &elem_error,
       if (error < threshold) { derefs.Append(i); }
    }
 
-   long glob_size = ReduceInt(derefs.Size());
+   long long glob_size = ReduceInt(derefs.Size());
    if (!glob_size) { return false; }
 
    // Destroy face-neighbor data only when actually de-refining.
@@ -5134,6 +5154,332 @@ void ParMesh::PrintAsOne(std::ostream &os) const
    }
 }
 
+void ParMesh::PrintAsSerial(std::ostream &os) const
+{
+   int save_rank = 0;
+   Mesh serialmesh = GetSerialMesh(save_rank);
+   if (MyRank == save_rank)
+   {
+      serialmesh.Printer(os);
+   }
+   MPI_Barrier(MyComm);
+}
+
+Mesh ParMesh::GetSerialMesh(int save_rank) const
+{
+   if (pncmesh || NURBSext)
+   {
+      MFEM_ABORT("Nonconforming meshes and NURBS meshes are not yet supported.");
+   }
+
+   // Define linear H1 space for vertex numbering
+   H1_FECollection fec_linear(1, Dim);
+   ParMesh *pm = const_cast<ParMesh *>(this);
+   ParFiniteElementSpace pfespace_linear(pm, &fec_linear);
+
+   long long ne_glob_l = GetGlobalNE(); // needs to be called by all ranks
+   MFEM_VERIFY(int(ne_glob_l) == ne_glob_l,
+               "overflow in the number of elements!");
+   int ne_glob = (save_rank == MyRank) ? int(ne_glob_l) : 0;
+
+   long long nvertices = pfespace_linear.GetTrueVSize();
+   long long nvertices_glob_l = 0;
+   MPI_Reduce(&nvertices, &nvertices_glob_l, 1, MPI_LONG_LONG, MPI_SUM,
+              save_rank, MyComm);
+   int nvertices_glob = int(nvertices_glob_l);
+   MFEM_VERIFY(nvertices_glob == nvertices_glob_l,
+               "overflow in the number of vertices!");
+
+   long long nbe = NumOfBdrElements;
+   long long nbe_glob_l = 0;
+   MPI_Reduce(&nbe, &nbe_glob_l, 1, MPI_LONG_LONG, MPI_SUM, save_rank, MyComm);
+   int nbe_glob = int(nbe_glob_l);
+   MFEM_VERIFY(nbe_glob == nbe_glob_l,
+               "overflow in the number of boundary elements!");
+
+   // On ranks other than save_rank, the *_glob variables are 0, so the serial
+   // mesh is empty.
+   Mesh serialmesh(Dim, nvertices_glob, ne_glob, nbe_glob, spaceDim);
+
+   int n_send_recv;
+   MPI_Status status;
+   Array<double> vert;
+   Array<int> ints, dofs;
+
+   // First set the connectivity of serial mesh using the True Dofs from
+   // the linear H1 space.
+   if (MyRank == save_rank)
+   {
+      for (int e = 0; e < NumOfElements; e++)
+      {
+         const int attr = elements[e]->GetAttribute();
+         const int geom_type = elements[e]->GetGeometryType();
+         pfespace_linear.GetElementDofs(e, dofs);
+         for (int j = 0; j < dofs.Size(); j++)
+         {
+            dofs[j] = pfespace_linear.GetGlobalTDofNumber(dofs[j]);
+         }
+         Element *elem = serialmesh.NewElement(geom_type);
+         elem->SetAttribute(attr);
+         elem->SetVertices(dofs);
+         serialmesh.AddElement(elem);
+      }
+
+      for (int p = 0; p < NRanks; p++)
+      {
+         if (p == save_rank) { continue; }
+         MPI_Recv(&n_send_recv, 1, MPI_INT, p, 444, MyComm, &status);
+         ints.SetSize(n_send_recv);
+         if (n_send_recv)
+         {
+            MPI_Recv(&ints[0], n_send_recv, MPI_INT, p, 445, MyComm, &status);
+         }
+         for (int i = 0; i < n_send_recv; )
+         {
+            int attr = ints[i++];
+            int geom_type = ints[i++];
+            Element *elem = serialmesh.NewElement(geom_type);
+            elem->SetAttribute(attr);
+            elem->SetVertices(&ints[i]); i += Geometry::NumVerts[geom_type];
+            serialmesh.AddElement(elem);
+         }
+      }
+   }
+   else
+   {
+      n_send_recv = 0;
+      for (int e = 0; e < NumOfElements; e++)
+      {
+         n_send_recv += 2 + elements[e]->GetNVertices();
+      }
+      MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 444, MyComm);
+      ints.Reserve(n_send_recv);
+      ints.SetSize(0);
+      for (int e = 0; e < NumOfElements; e++)
+      {
+         const int attr = elements[e]->GetAttribute();
+         const int geom_type = elements[e]->GetGeometryType();;
+         ints.Append(attr);
+         ints.Append(geom_type);
+         pfespace_linear.GetElementDofs(e, dofs);
+         for (int j = 0; j < dofs.Size(); j++)
+         {
+            ints.Append(pfespace_linear.GetGlobalTDofNumber(dofs[j]));
+         }
+      }
+      if (n_send_recv)
+      {
+         MPI_Send(&ints[0], n_send_recv, MPI_INT, save_rank, 445, MyComm);
+      }
+   }
+
+   // Write out boundary elements
+   if (MyRank == save_rank)
+   {
+      for (int e = 0; e < NumOfBdrElements; e++)
+      {
+         const int attr = boundary[e]->GetAttribute();
+         const int geom_type = boundary[e]->GetGeometryType();
+         pfespace_linear.GetBdrElementDofs(e, dofs);
+         for (int j = 0; j < dofs.Size(); j++)
+         {
+            dofs[j] = pfespace_linear.GetGlobalTDofNumber(dofs[j]);
+         }
+         Element *elem = serialmesh.NewElement(geom_type);
+         elem->SetAttribute(attr);
+         elem->SetVertices(dofs);
+         serialmesh.AddBdrElement(elem);
+      }
+
+      for (int p = 0; p < NRanks; p++)
+      {
+         if (p == save_rank) { continue; }
+         MPI_Recv(&n_send_recv, 1, MPI_INT, p, 446, MyComm, &status);
+         ints.SetSize(n_send_recv);
+         if (n_send_recv)
+         {
+            MPI_Recv(&ints[0], n_send_recv, MPI_INT, p, 447, MyComm, &status);
+         }
+         for (int i = 0; i < n_send_recv; )
+         {
+            int attr = ints[i++];
+            int geom_type = ints[i++];
+            Element *elem = serialmesh.NewElement(geom_type);
+            elem->SetAttribute(attr);
+            elem->SetVertices(&ints[i]); i += Geometry::NumVerts[geom_type];
+            serialmesh.AddBdrElement(elem);
+         }
+      }
+   } // MyRank == save_rank
+   else
+   {
+      n_send_recv = 0;
+      for (int e = 0; e < NumOfBdrElements; e++)
+      {
+         n_send_recv += 2 + GetBdrElement(e)->GetNVertices();
+      }
+      MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 446, MyComm);
+      ints.Reserve(n_send_recv);
+      ints.SetSize(0);
+      for (int e = 0; e < NumOfBdrElements; e++)
+      {
+         const int attr = boundary[e]->GetAttribute();
+         const int geom_type = boundary[e]->GetGeometryType();
+         ints.Append(attr);
+         ints.Append(geom_type);
+         pfespace_linear.GetBdrElementDofs(e, dofs);
+         for (int j = 0; j < dofs.Size(); j++)
+         {
+            ints.Append(pfespace_linear.GetGlobalTDofNumber(dofs[j]));
+         }
+      }
+      if (n_send_recv)
+      {
+         MPI_Send(&ints[0], n_send_recv, MPI_INT, save_rank, 447, MyComm);
+      }
+   } // MyRank != save_rank
+
+   if (MyRank == save_rank)
+   {
+      for (int v = 0; v < nvertices_glob; v++)
+      {
+         serialmesh.AddVertex(0.0); // all other coordinates are 0 by default
+      }
+      serialmesh.FinalizeTopology();
+   }
+
+   // From each processor, we send element-wise vertex/dof locations and
+   // overwrite the vertex/dof locations of the serial mesh.
+   if (MyRank == save_rank && Nodes)
+   {
+      FiniteElementSpace *fespace_serial = NULL;
+      // Duplicate the FE collection to make sure the serial mesh is completely
+      // independent of the parallel mesh:
+      auto fec_serial = FiniteElementCollection::New(
+                           GetNodalFESpace()->FEColl()->Name());
+      fespace_serial = new FiniteElementSpace(&serialmesh,
+                                              fec_serial,
+                                              spaceDim,
+                                              GetNodalFESpace()->GetOrdering());
+      serialmesh.SetNodalFESpace(fespace_serial);
+      serialmesh.GetNodes()->MakeOwner(fec_serial);
+      // The serial mesh owns its Nodes and they, in turn, own fec_serial and
+      // fespace_serial.
+   }
+
+   int elem_count = 0; // To keep track of element count in serial mesh
+   if (MyRank == save_rank)
+   {
+      Vector nodeloc;
+      Array<int> ints_serial;
+      for (int e = 0; e < NumOfElements; e++)
+      {
+         if (Nodes)
+         {
+            Nodes->GetElementDofValues(e, nodeloc);
+            serialmesh.GetNodalFESpace()->GetElementVDofs(elem_count++, dofs);
+            serialmesh.GetNodes()->SetSubVector(dofs, nodeloc);
+         }
+         else
+         {
+            GetElementVertices(e, ints);
+            serialmesh.GetElementVertices(elem_count++, ints_serial);
+            for (int i = 0; i < ints.Size(); i++)
+            {
+               const double *vdata = GetVertex(ints[i]);
+               double *vdata_serial = serialmesh.GetVertex(ints_serial[i]);
+               for (int d = 0; d < spaceDim; d++)
+               {
+                  vdata_serial[d] = vdata[d];
+               }
+            }
+         }
+      }
+
+      for (int p = 0; p < NRanks; p++)
+      {
+         if (p == save_rank) { continue; }
+         MPI_Recv(&n_send_recv, 1, MPI_INT, p, 448, MyComm, &status);
+         vert.SetSize(n_send_recv);
+         if (n_send_recv)
+         {
+            MPI_Recv(&vert[0], n_send_recv, MPI_DOUBLE, p, 449, MyComm, &status);
+         }
+         for (int i = 0; i < n_send_recv; )
+         {
+            if (Nodes)
+            {
+               serialmesh.GetNodalFESpace()->GetElementVDofs(elem_count++, dofs);
+               serialmesh.GetNodes()->SetSubVector(dofs, &vert[i]);
+               i += dofs.Size();
+            }
+            else
+            {
+               serialmesh.GetElementVertices(elem_count++, ints_serial);
+               for (int j = 0; j < ints_serial.Size(); j++)
+               {
+                  double *vdata_serial = serialmesh.GetVertex(ints_serial[j]);
+                  for (int d = 0; d < spaceDim; d++)
+                  {
+                     vdata_serial[d] = vert[i++];
+                  }
+               }
+            }
+         }
+      }
+   } // MyRank == save_rank
+   else
+   {
+      n_send_recv = 0;
+      Vector nodeloc;
+      for (int e = 0; e < NumOfElements; e++)
+      {
+         if (Nodes)
+         {
+            const FiniteElement *fe = Nodes->FESpace()->GetFE(e);
+            n_send_recv += spaceDim*fe->GetDof();
+         }
+         else
+         {
+            n_send_recv += elements[e]->GetNVertices()*spaceDim;
+         }
+      }
+      MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 448, MyComm);
+      vert.Reserve(n_send_recv);
+      vert.SetSize(0);
+      for (int e = 0; e < NumOfElements; e++)
+      {
+         if (Nodes)
+         {
+            Nodes->GetElementDofValues(e, nodeloc);
+            for (int j = 0; j < nodeloc.Size(); j++)
+            {
+               vert.Append(nodeloc(j));
+            }
+         }
+         else
+         {
+            GetElementVertices(e, ints);
+            for (int i = 0; i < ints.Size(); i++)
+            {
+               const double *vdata = GetVertex(ints[i]);
+               for (int d = 0; d < spaceDim; d++)
+               {
+                  vert.Append(vdata[d]);
+               }
+            }
+         }
+      }
+      if (n_send_recv)
+      {
+         MPI_Send(&vert[0], n_send_recv, MPI_DOUBLE, save_rank, 449, MyComm);
+      }
+   }
+
+   MPI_Barrier(MyComm);
+   return serialmesh;
+}
+
 void ParMesh::SaveAsOne(const char *fname, int precision) const
 {
    ofstream ofs;
@@ -5743,8 +6089,8 @@ void ParMesh::PrintInfo(std::ostream &os)
 
    // TODO: collect and print stats by geometry
 
-   long ldata[5]; // vert, edge, face, elem, neighbors;
-   long mindata[5], maxdata[5], sumdata[5];
+   long long ldata[5]; // vert, edge, face, elem, neighbors;
+   long long mindata[5], maxdata[5], sumdata[5];
 
    // count locally owned vertices, edges, and faces
    ldata[0] = GetNV();
@@ -5763,9 +6109,9 @@ void ParMesh::PrintInfo(std::ostream &os)
       }
    }
 
-   MPI_Reduce(ldata, mindata, 5, MPI_LONG, MPI_MIN, 0, MyComm);
-   MPI_Reduce(ldata, sumdata, 5, MPI_LONG, MPI_SUM, 0, MyComm);
-   MPI_Reduce(ldata, maxdata, 5, MPI_LONG, MPI_MAX, 0, MyComm);
+   MPI_Reduce(ldata, mindata, 5, MPI_LONG_LONG, MPI_MIN, 0, MyComm);
+   MPI_Reduce(ldata, sumdata, 5, MPI_LONG_LONG, MPI_SUM, 0, MyComm);
+   MPI_Reduce(ldata, maxdata, 5, MPI_LONG_LONG, MPI_MAX, 0, MyComm);
 
    if (MyRank == 0)
    {
@@ -5816,10 +6162,10 @@ void ParMesh::PrintInfo(std::ostream &os)
    }
 }
 
-long ParMesh::ReduceInt(int value) const
+long long ParMesh::ReduceInt(int value) const
 {
-   long local = value, global;
-   MPI_Allreduce(&local, &global, 1, MPI_LONG, MPI_SUM, MyComm);
+   long long local = value, global;
+   MPI_Allreduce(&local, &global, 1, MPI_LONG_LONG, MPI_SUM, MyComm);
    return global;
 }
 
@@ -6180,7 +6526,8 @@ void ParMesh::GetGlobalElementIndices(Array<HYPRE_BigInt> &gi) const
 {
    ComputeGlobalElementOffset();
 
-   const HYPRE_BigInt offset = glob_elem_offset; // Cast from long to HYPRE_BigInt
+   // Cast from long long to HYPRE_BigInt
+   const HYPRE_BigInt offset = glob_elem_offset;
 
    gi.SetSize(GetNE());
    for (int i=0; i<GetNE(); ++i)

mesh/pmesh.hpp

@@ -81,7 +81,8 @@ protected:
    IsoparametricTransformation FaceNbrTransformation;
 
    // glob_elem_offset + local element number defines a global element numbering
-   mutable long glob_elem_offset, glob_offset_sequence;
+   mutable long long glob_elem_offset;
+   mutable long glob_offset_sequence;
    void ComputeGlobalElementOffset() const;
 
    // Enable Print() to add the parallel interface as boundary (typically used
@@ -295,10 +296,10 @@ public:
 
    /** Map a global element number to a local element number. If the global
        element is not on this processor, return -1. */
-   int GetLocalElementNum(long global_element_num) const;
+   int GetLocalElementNum(long long global_element_num) const;
 
    /// Map a local element number to a global element number.
-   long GetGlobalElementNum(int local_element_num) const;
+   long long GetGlobalElementNum(int local_element_num) const;
 
    /** The following functions define global indices for all local vertices,
        edges, faces, or elements. The global indices have no meaning or
@@ -352,6 +353,12 @@ public:
    void SetCurvature(int order, bool discont = false, int space_dim = -1,
                      int ordering = 1) override;
 
+   /** Replace the internal node GridFunction with a new GridFunction defined
+       on the given FiniteElementSpace. The new node coordinates are projected
+       (derived) from the current nodes/vertices. */
+   void SetNodalFESpace(FiniteElementSpace *nfes) override;
+   void SetNodalFESpace(ParFiniteElementSpace *npfes);
+
    int GetNFaceNeighbors() const { return face_nbr_group.Size(); }
    int GetNFaceNeighborElements() const { return face_nbr_elements.Size(); }
    int GetFaceNbrGroup(int fn) const { return face_nbr_group[fn]; }
@@ -446,7 +453,7 @@ public:
    MFEM_DEPRECATED void ReorientTetMesh() override;
 
    /// Utility function: sum integers from all processors (Allreduce).
-   long ReduceInt(int value) const override;
+   long long ReduceInt(int value) const override;
 
    /** Load balance the mesh by equipartitioning the global space-filling
        sequence of elements. Works for nonconforming meshes only. */
@@ -491,6 +498,15 @@ public:
        attributes are set to the processor number.  */
    void PrintAsOne(std::ostream &out = mfem::out) const;
 
+   /** Write the mesh to the stream 'out' on Process 0 as a serial mesh. The
+       output mesh does not have any duplication of vertices/nodes at
+       processor boundaries. */
+   void PrintAsSerial(std::ostream &out = mfem::out) const;
+
+   /** Returns a Serial mesh on MPI rank @a save_rank that does not have any
+       duplication of vertices/nodes at processor boundaries. */
+   Mesh GetSerialMesh(int save_rank) const;
+
    /// Save the mesh as a single file (using ParMesh::PrintAsOne). The given
    /// @a precision is used for ASCII output.
    void SaveAsOne(const char *fname, int precision=16) const;

mesh/pncmesh.cpp

@@ -1330,8 +1330,8 @@ void ParNCMesh::LimitNCLevel(int max_nc_level)
       Array<Refinement> refinements;
       GetLimitRefinements(refinements, max_nc_level);
 
-      long size = refinements.Size(), glob_size;
-      MPI_Allreduce(&size, &glob_size, 1, MPI_LONG, MPI_SUM, MyComm);
+      long long size = refinements.Size(), glob_size;
+      MPI_Allreduce(&size, &glob_size, 1, MPI_LONG_LONG, MPI_SUM, MyComm);
 
       if (!glob_size) { break; }
 
@@ -2755,15 +2755,15 @@ void ParNCMesh::Trim()
    ClearAuxPM();
 }
 
-long ParNCMesh::RebalanceDofMessage::MemoryUsage() const
+std::size_t ParNCMesh::RebalanceDofMessage::MemoryUsage() const
 {
    return (elem_ids.capacity() + dofs.capacity()) * sizeof(int);
 }
 
 template<typename K, typename V>
-static long map_memory_usage(const std::map<K, V> &map)
+static std::size_t map_memory_usage(const std::map<K, V> &map)
 {
-   long result = 0;
+   std::size_t result = 0;
    for (typename std::map<K, V>::const_iterator
         it = map.begin(); it != map.end(); ++it)
    {
@@ -2773,9 +2773,9 @@ static long map_memory_usage(const std::map<K, V> &map)
    return result;
 }
 
-long ParNCMesh::GroupsMemoryUsage() const
+std::size_t ParNCMesh::GroupsMemoryUsage() const
 {
-   long groups_size = groups.capacity() * sizeof(CommGroup);
+   std::size_t groups_size = groups.capacity() * sizeof(CommGroup);
    for (unsigned i = 0; i < groups.size(); i++)
    {
       groups_size += groups[i].capacity() * sizeof(int);
@@ -2786,9 +2786,9 @@ long ParNCMesh::GroupsMemoryUsage() const
 }
 
 template<typename Type, int Size>
-static long arrays_memory_usage(const Array<Type> (&arrays)[Size])
+static std::size_t arrays_memory_usage(const Array<Type> (&arrays)[Size])
 {
-   long total = 0;
+   std::size_t total = 0;
    for (int i = 0; i < Size; i++)
    {
       total += arrays[i].MemoryUsage();
@@ -2796,9 +2796,9 @@ static long arrays_memory_usage(const Array<Type> (&arrays)[Size])
    return total;
 }
 
-long ParNCMesh::MemoryUsage(bool with_base) const
+std::size_t ParNCMesh::MemoryUsage(bool with_base) const
 {
-   long total_groups_owners = 0;
+   std::size_t total_groups_owners = 0;
    for (int i = 0; i < 3; i++)
    {
       total_groups_owners += entity_owner[i].MemoryUsage() +

mesh/pncmesh.hpp

@@ -240,7 +240,7 @@ public:
    virtual void Trim();
 
    /// Return total number of bytes allocated.
-   long MemoryUsage(bool with_base = true) const;
+   std::size_t MemoryUsage(bool with_base = true) const;
 
    int PrintMemoryDetail(bool with_base = true) const;
 
@@ -506,7 +506,7 @@ protected: // implementation
 
       void SetElements(const Array<int> &elems, NCMesh *ncmesh);
       void SetNCMesh(NCMesh* ncmesh) { eset.SetNCMesh(ncmesh); }
-      long MemoryUsage() const;
+      std::size_t MemoryUsage() const;
 
       typedef std::map<int, RebalanceDofMessage> Map;
 
@@ -540,7 +540,7 @@ protected: // implementation
    Array<DenseMatrix*> aux_pm_store;
    void ClearAuxPM();
 
-   long GroupsMemoryUsage() const;
+   std::size_t GroupsMemoryUsage() const;
 
    friend class NeighborRowMessage;
 };

mesh/vtk.cpp

@@ -106,9 +106,9 @@ int VTKGeometry::GetOrder(int vtk_geom, int npoints)
          case LAGRANGE_SEGMENT:
             return npoints - 1;
          case LAGRANGE_TRIANGLE:
-            return (std::sqrt(8*npoints + 1) - 3)/2;
+            return static_cast<int>(std::sqrt(8*npoints + 1) - 3)/2;
          case LAGRANGE_SQUARE:
-            return std::round(std::sqrt(npoints)) - 1;
+            return static_cast<int>(std::round(std::sqrt(npoints))) - 1;
          case LAGRANGE_TETRAHEDRON:
             switch (npoints)
             {
@@ -138,7 +138,7 @@ int VTKGeometry::GetOrder(int vtk_geom, int npoints)
                }
             }
          case LAGRANGE_CUBE:
-            return std::round(std::cbrt(npoints)) - 1;
+            return static_cast<int>(std::round(std::cbrt(npoints))) - 1;
          case LAGRANGE_PRISM:
          {
             const double n = npoints;
@@ -148,7 +148,7 @@ int VTKGeometry::GetOrder(int vtk_geom, int npoints)
             const double term =
                std::cbrt(third*sqrt(third)*sqrt((27.0*n - 2.0)*n) + n
                          - twentyseventh);
-            return std::round(term + ninth / term - 4*third);
+            return static_cast<int>(std::round(term + ninth / term - 4*third));
          }
          case LAGRANGE_PYRAMID:
             MFEM_ABORT("Lagrange pyramids not currently supported in VTK.");

miniapps/autodiff/par_example.cpp

@@ -526,11 +526,11 @@ int main(int argc, char *argv[])
       dacol->SetTime(pp);
       if (pp < 2.0)
       {
-         dacol->SetCycle(std::floor(pp));
+         dacol->SetCycle(static_cast<int>(std::floor(pp)));
       }
       else
       {
-         dacol->SetCycle(std::ceil(pp));
+         dacol->SetCycle(static_cast<int>(std::ceil(pp)));
       }
       dacol->Save();
    }

miniapps/autodiff/seq_example.cpp

@@ -461,11 +461,11 @@ int main(int argc, char *argv[])
       dacol->SetTime(pp);
       if (pp < 2.0)
       {
-         dacol->SetCycle(std::floor(pp));
+         dacol->SetCycle(static_cast<int>(std::floor(pp)));
       }
       else
       {
-         dacol->SetCycle(std::ceil(pp));
+         dacol->SetCycle(static_cast<int>(std::ceil(pp)));
       }
       dacol->Save();
    }

miniapps/electromagnetics/maxwell.cpp

@@ -479,7 +479,7 @@ void dipole_pulse(const Vector &x, double t, Vector &j)
       j = v;
    }
 
-   j *= a * (t - b) * exp(-0.5 * pow((t-b)/c, 2)) / (c * c);
+   j *= a * (t - b) * exp(-0.5 * pow((t-b)/c, 2.0)) / (c * c);
 }
 
 void
@@ -508,7 +508,7 @@ SnapTimeStep(double tmax, double dtmax, double & dt)
 {
    double dsteps = tmax/dtmax;
 
-   int nsteps = pow(10,(int)ceil(log10(dsteps)));
+   int nsteps = static_cast<int>(pow(10,(int)ceil(log10(dsteps))));
 
    for (int i=1; i<=5; i++)
    {

miniapps/electromagnetics/tesla.cpp

@@ -540,10 +540,10 @@ void halbach_array(const Vector &x, Vector &m)
    int ri = (int)ha_params_[7];
    int n  = (int)ha_params_[8];
 
-   int i = (int)n * (x[ai] - ha_params_[ai]) /
-           (ha_params_[ai+3] - ha_params_[ai]);
+   int i = static_cast<int>(n * (x[ai] - ha_params_[ai]) /
+                            (ha_params_[ai+3] - ha_params_[ai]));
 
-   m[(ri + 1 + (i % 2)) % 3] = pow(-1.0,i/2);
+   m[(ri + 1 + (i % 2)) % 3] = static_cast<int>(pow(-1.0,i/2));
 }
 
 // To produce a uniform magnetic flux the vector potential can be set

miniapps/electromagnetics/tesla_solver.cpp

@@ -414,7 +414,7 @@ TeslaSolver::GetErrorEstimates(Vector & errors)
    ParFiniteElementSpace flux_fes(pmesh_, &flux_fec);
 
    // Space for the smoothed (conforming) flux
-   double norm_p = 1;
+   int norm_p = 1;
    ND_FECollection smooth_flux_fec(order_, pmesh_->Dimension());
    ParFiniteElementSpace smooth_flux_fes(pmesh_, &smooth_flux_fec);
 

miniapps/electromagnetics/volta_solver.cpp

@@ -514,7 +514,7 @@ VoltaSolver::GetErrorEstimates(Vector & errors)
    ParFiniteElementSpace flux_fes(pmesh_, &flux_fec, pmesh_->SpaceDimension());
 
    // Space for the smoothed (conforming) flux
-   double norm_p = 1;
+   int norm_p = 1;
    RT_FECollection smooth_flux_fec(order_-1, pmesh_->Dimension());
    ParFiniteElementSpace smooth_flux_fes(pmesh_, &smooth_flux_fec);
 

miniapps/meshing/mesh-optimizer.cpp

@@ -864,7 +864,7 @@ int main(int argc, char *argv[])
       for (int i = 0; i < mesh->GetNE(); i++)
       {
          mat(i) = material_id(i, surf_fit_gf0);
-         mesh->SetAttribute(i, mat(i) + 1);
+         mesh->SetAttribute(i, static_cast<int>(mat(i) + 1));
       }
 
       GridFunctionCoefficient mat_coeff(&mat);