Adding first draft of MixedSesquilinearForm class

Fixing HiOp int types

.github/workflows/repo-check.yml

@@ -63,7 +63,7 @@ jobs:
         exit 1
 
   code-style:
-    runs-on: ubuntu-16.04 # needed for astyle 2.05.1
+    runs-on: ubuntu-18.04
 
     steps:
     - name: checkout mfem
@@ -71,7 +71,7 @@ jobs:
 
     - name: get astyle
       run: |
-        sudo apt-get install astyle=2.05.1-0ubuntu1
+        sudo apt-get install astyle=3.1-1ubuntu2
 
     - name: style check
       run: |

.gitlab-ci.yml

@@ -16,7 +16,7 @@
 # We define the following GitLab pipeline variables:
 #
 # BUILD_ROOT:
-# The path to the shared resources between all jobs. For example, external
+# The path to the shared resource between all jobs. For example, external
 # repositories like 'tests' and 'tpls' are cloned here. Also, 'tpls' is built
 # once for all targets, so that build happen here. The BUILD_ROOT is unique to
 # the pipeline, preventing any form of concurrency with other pipelines. This
@@ -30,7 +30,7 @@
 #
 # MFEM_ALLOC_NAME:
 # On LLNL's quartz, there is only one allocation shared among jobs in order to
-# save time and resources. This allocation has to be uniquely named so that we
+# save time and resource. This allocation has to be uniquely named so that we
 # are sure to retrieve it.
 #
 # TPLS_REPO & TESTS_REPO:
@@ -57,18 +57,20 @@ variables:
 # the "needs" keyword and express the DAG of jobs for more efficiency.
 # - We use setup and setup_baseline phases to download content outside of mfem
 #   directory.
-# - Allocate/Release is where quartz resources are allocated/released once for all.
+# - Allocate/Release is where quartz resource are allocated/released once for all.
 # - Build and Test is where we build and MFEM for multiple toolchains.
 # - Baseline_checks gathers baseline-type test suites execution
 # - Baseline_publish, only available on master, allows to update baseline
 #   results
 stages:
   - setup
-  - q_allocate_resources
+  - q_allocate_resource
   - q_build_and_test
-  - q_release_resources
+  - q_release_resource
   - l_build_and_test
+  - c_allocate_resource
   - c_build_and_test
+  - c_release_resource
   - setup_baseline
   - baseline_check
   - baseline_to_autotest
@@ -120,11 +122,11 @@ setup_baseline:
     - echo ${ALLOC_NAME}
     - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
     - echo ${JOBID}
-    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -t 30 -N 1 tests/gitlab/build_and_test
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -t 30 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --build-root "${BUILD_ROOT}" --data
 
 .build_toss_3_x86_64_ib_corona_script:
   script:
-    - srun -p mi60 -t 15 -N 1 tests/gitlab/build_and_test
+    - srun -p mi60 -t 15 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --build-root "${BUILD_ROOT}" --data
 
 # Lassen uses a different job scheduler (spectrum lsf) that does not allow
 # pre-allocation the same way slurm does.  We use pdebug queue on lassen to
@@ -132,7 +134,7 @@ setup_baseline:
 # builds.
 .build_blueos_3_ppc64le_ib_script:
   script:
-    - lalloc 1 -W 30 -q pdebug tests/gitlab/build_and_test
+    - lalloc 1 -W 30 -q pdebug tests/gitlab/build_and_test --spec "${SPEC}" --build-root "${BUILD_ROOT}" --data
 
 # Shared script for baseline and sample-run-baseline, the value of BASELINE_TEST
 # differentiates between the two tests.
@@ -253,3 +255,4 @@ setup_baseline:
 include:
   - local: .gitlab/quartz.yml
   - local: .gitlab/lassen.yml
+  - local: .gitlab/corona.yml

.gitlab/corona.yml

@@ -0,0 +1,127 @@
+# Copyright (c) 2010-2021, 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.
+
+# GitLab pipeline configuration for the Corona machine at LLNL
+
+.on_corona:
+  tags:
+    - shell
+    - corona
+  rules:
+    # Don’t run corona jobs if...
+    # Note: This makes corona an "opt-in" machine. To activate builds on corona
+    # for a given GitLab clone of MFEM, go to Setting/CI-CD/variables, and set
+    # "ON_CORONA" to "ON". An LC account on for corona is required to trigger a
+    # pipeline there.
+    - if: '$CI_COMMIT_BRANCH =~ /_cnone/ || $ON_CORONA != "ON"'
+      when: never
+    # Don’t run autotest update if...
+    - if: '$CI_JOB_NAME =~ /update_autotest/ && $AUTOTEST != "YES"'
+      when: never
+    # Don’t run autotest update if...
+    - if: '$CI_JOB_NAME =~ /c_report/ && $AUTOTEST != "YES"'
+      when: never
+    # Report success on success status
+    - if: '$CI_JOB_NAME =~ /c_report_success/ && $AUTOTEST == "YES"'
+      when: on_success
+    # Report failure on failure status
+    - if: '$CI_JOB_NAME =~ /c_report_failure/ && $AUTOTEST == "YES"'
+      when: on_failure
+    # Always release resource
+    - if: '$CI_JOB_NAME =~ /release_resource/'
+      when: always
+    # Default is to run if previous stage succeeded
+    - when: on_success
+
+# This is a yaml anchor, it can be used to avoid duplication like here.
+# The code below will simply be pasted wherever the anchor is placed.
+.safe_create_rundir: &safe_create_rundir |
+  if ! mkdir ${rundir}; then
+    n=1
+    while ! mkdir ${rundir}_${n}
+    do
+      n=$((n+1))
+    done
+    rundir=${rundir}_${n}
+  fi
+
+# Allocate
+c_allocate_resource:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_corona
+  stage: c_allocate_resource
+  script:
+    - salloc --exclusive --nodes=1 --partition=mi60 --time=30 --no-shell --job-name=${ALLOC_NAME}
+  timeout: 6h
+  needs: [setup]
+
+# Release
+c_release_resource:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_corona
+  stage: c_release_resource
+  script:
+    - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
+    - ([[ -n "${JOBID}" ]] && scancel ${JOBID})
+  needs: [rocm_gcc_8.3.1]
+
+# Release
+c_report_success:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_corona
+  stage: c_release_resource
+  script:
+    - echo "Can only run if all the corona jobs passed"
+    - cd ${AUTOTEST_ROOT}/autotest && git pull
+    - rundir="corona/$(date +%Y-%m-%d)-gitlab-ci-${CI_COMMIT_REF_SLUG}"
+    - mkdir -p corona
+    - *safe_create_rundir
+    - echo "The Corona jobs were successful" > ${rundir}/gitlab.out
+    - echo "See the pipeline here -> $CI_PIPELINE_URL" >> ${rundir}/gitlab.err
+    - git add ${rundir}
+    - git commit -am "Gitlab CI log for baseline on corona with intel ($(date +%Y-%m-%d))"
+    - git push origin master
+  needs: [rocm_gcc_8.3.1]
+
+c_report_failure:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_corona
+  stage: c_release_resource
+  script:
+    - echo "Runs if there was at least one failure on corona"
+    - cd ${AUTOTEST_ROOT}/autotest && git pull
+    - rundir="corona/$(date +%Y-%m-%d)-gitlab-ci-${CI_COMMIT_REF_SLUG}"
+    - mkdir -p corona
+    - *safe_create_rundir
+    - echo "There was an error while running CI on Corona" > ${rundir}/gitlab.err
+    - cp ${rundir}/gitlab.err ${rundir}/autotest-email.html
+    - git add ${rundir}
+    - git commit -am "Gitlab CI log for baseline on corona with intel ($(date +%Y-%m-%d))"
+    - git push origin master
+  needs: [rocm_gcc_8.3.1]
+
+# Spack helped builds
+# Generic corona build job, extending build script
+.build_and_test_on_corona:
+  extends: [.build_toss_3_x86_64_ib_corona_script, .on_corona]
+  stage: c_build_and_test
+
+# Build MFEM
+rocm_gcc_8.3.1:
+  variables:
+    SPEC: "@develop%gcc@8.3.1+rocm amdgpu_target=gfx906"
+  extends: .build_and_test_on_corona
+  needs: [c_allocate_resource]
+

.gitlab/quartz.yml

@@ -31,8 +31,8 @@
     # Report failure on failure status
     - if: '$CI_JOB_NAME =~ /q_report_failure/ && $AUTOTEST == "YES"'
       when: on_failure
-    # Always release resources
-    - if: '$CI_JOB_NAME =~ /release_resources/'
+    # Always release resource
+    - if: '$CI_JOB_NAME =~ /release_resource/'
       when: always
     # Default is to run if previous stage succeeded
     - when: on_success
@@ -50,21 +50,21 @@
   fi
 
 # Allocate
-q_allocate_resources:
+q_allocate_resource:
   variables:
     GIT_STRATEGY: none
   extends: .on_quartz
-  stage: q_allocate_resources
+  stage: q_allocate_resource
   script:
     - salloc --exclusive --nodes=1 --partition=pdebug --time=30 --no-shell --job-name=${ALLOC_NAME}
   timeout: 6h
 
 # Release
-q_release_resources:
+q_release_resource:
   variables:
     GIT_STRATEGY: none
   extends: .on_quartz
-  stage: q_release_resources
+  stage: q_release_resource
   script:
     - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
     - ([[ -n "${JOBID}" ]] && scancel ${JOBID})
@@ -74,11 +74,12 @@ q_report_success:
   variables:
     GIT_STRATEGY: none
   extends: .on_quartz
-  stage: q_release_resources
+  stage: q_release_resource
   script:
     - echo "Can only run if all the quartz jobs passed"
     - cd ${AUTOTEST_ROOT}/autotest && git pull
-    - rundir="gitlab/$(date +%Y-%m-%d)-github-${CI_COMMIT_REF_SLUG}"
+    - rundir="quartz/$(date +%Y-%m-%d)-gitlab-ci-${CI_COMMIT_REF_SLUG}"
+    - mkdir -p quartz
     - *safe_create_rundir
     - echo "The Quartz jobs were successful" > ${rundir}/gitlab.out
     - echo "See the pipeline here -> $CI_PIPELINE_URL" >> ${rundir}/gitlab.err
@@ -90,11 +91,12 @@ q_report_failure:
   variables:
     GIT_STRATEGY: none
   extends: .on_quartz
-  stage: q_release_resources
+  stage: q_release_resource
   script:
     - echo "Runs if there was at least one failure on quartz"
     - cd ${AUTOTEST_ROOT}/autotest && git pull
-    - rundir="gitlab/$(date +%Y-%m-%d)-github-${CI_COMMIT_REF_SLUG}"
+    - rundir="quartz/$(date +%Y-%m-%d)-gitlab-ci-${CI_COMMIT_REF_SLUG}"
+    - mkdir -p quartz
     - *safe_create_rundir
     - echo "There was an error while running CI on Quartz" > ${rundir}/gitlab.err
     - echo "See the pipeline here -> $CI_PIPELINE_URL" >> ${rundir}/gitlab.err
@@ -165,7 +167,8 @@ update_autotest:
   stage: baseline_to_autotest
   script:
     - cd ${AUTOTEST_ROOT}/autotest && git pull
-    - rundir="quartz/$(date +%Y-%m-%d)-github-${CI_COMMIT_REF_SLUG}"
+    - rundir="quartz/$(date +%Y-%m-%d)-gitlab-baseline-${CI_COMMIT_REF_SLUG}"
+    - mkdir -p quartz
     - *safe_create_rundir
     - cp ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/* ${rundir}
     # We create an autotest-email.html file, because that's how we signal that there was a diff (temporary).

CHANGELOG

@@ -10,8 +10,13 @@
 
 Version 4.3.1 (development)
 ===========================
+- Switched from Artistic Style (astyle) version 2.05.1 to version 3.1 for code
+  formatting. See the "make style" target.
+
 - Added support for hr-adaptivity using TMOP-based error estimator.
 
+- Added 'double' atomicAdd implementation for previous versions of CUDA.
+
 - Adding lowest order Nedelec and Raviart-Thomas basis functions on wedge
   shaped elements.
 
@@ -27,6 +32,12 @@ Version 4.3.1 (development)
   functions on wedges and pyramids which are not amenable to reordering. The
   ReorientTetMesh method of the Mesh and ParMesh classes has been deprecated.
 
+- 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.
+
+- Fixed several MinGW build issues on Windows.
+
 
 Version 4.3, released on July 29, 2021
 ======================================

CMakeLists.txt

@@ -90,6 +90,11 @@ include("${CMAKE_CURRENT_SOURCE_DIR}/config/XSDKDefaults.cmake")
 
 # Enable languages.
 enable_language(CXX)
+if (MINGW)
+   # MinGW GCC does not expose the functions jn/_jn, yn/_yn (used in Example
+   # 25/25p) unless we use '-std=gnu++11':
+   set(CMAKE_CXX_EXTENSIONS ON)
+endif()
 if (MFEM_USE_CUDA)
    if (MFEM_USE_HIP)
       message(FATAL_ERROR " *** MFEM_USE_HIP cannot be combined with MFEM_USE_CUDA.")
@@ -175,10 +180,6 @@ else()
   set(MFEM_DEBUG OFF)
 endif()
 
-if (WIN32)
-	add_definitions(-D_USE_MATH_DEFINES)
-endif()
-
 # MPI -> hypre; PETSc (optional)
 if (MFEM_USE_MPI)
   find_package(MPI REQUIRED)
@@ -432,10 +433,10 @@ endif()
 #    With newer versions of SuiteSparse which include METIS header using 64-bit
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
-set(MFEM_TPLS MPI_CXX OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS PETSC
+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)
+    CUSPARSE MKL_CPARDISO AMGX CALIPER MPI_CXX)
 
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
@@ -500,13 +501,9 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX} CACHE PATH
 # Declaring the library
 mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
-if (CMAKE_VERSION VERSION_GREATER 2.8.11)
-  target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES})
-else()
-  target_link_libraries(mfem ${TPL_LIBRARIES})
-endif()
+target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES})
 if (MINGW)
-  target_link_libraries(mfem ws2_32)
+  target_link_libraries(mfem PRIVATE ws2_32)
 endif()
 set_target_properties(mfem PROPERTIES VERSION "${mfem_VERSION}")
 set_target_properties(mfem PROPERTIES SOVERSION "${mfem_VERSION}")

CONTRIBUTING.md

@@ -327,7 +327,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
     documentation in source comments.
   - Consistent code styling is enforced with `make style` in the top-level
     directory. This requires [Artistic Style](http://astyle.sourceforge.net) (we
-    specifically use version 2.05.1). See also the file `config/mfem.astylerc`.
+    specifically use version 3.1). See also the file `config/mfem.astylerc`.
   - Use `mfem::out` and `mfem::err` instead of `std::cout` and `std::cerr` in
     internal library code. (You can use `std` in examples and miniapps.)
   - When manually resolving conflicts during a merge, make sure to mention the

INSTALL

@@ -692,7 +692,7 @@ The specific libraries and their options are:
 - HiOp (optional), used when MFEM_USE_HIOP = YES.
   URL: https://github.com/LLNL/hiop
   Options: HIOP_OPT, HIOP_LIB.
-  Versions: HIOP >= 0.4.
+  Versions: HIOP >= 0.4.6.
 
 - GSLIB (optional), used when MFEM_USE_GSLIB = YES. The gslib library must be
   built prior to the MFEM build, as follows: download gslib-1.0.7, untar it at

config/githooks/pre-push

@@ -87,12 +87,12 @@ fi
 ## style check
 #if [[ "${option}" == "--style" || "${option}" == "" ]]; then
 if [[ "${option}" == "--style" ]]; then
-  if which astyle && [[ "$(astyle --version)" == "Artistic Style Version 2.05.1" ]]; then
+  if which astyle && [[ "$(astyle --version)" == "Artistic Style Version 3.1" ]]; then
     cd tests/scripts
     if ! ./runtest code-style; then code=1; fi
     cd -
   else
-    echo "Warning: astyle not found or version is not 2.05.1"
+    echo "Warning: astyle not found or version is not 3.1"
   fi
 fi
 

examples/ex27.cpp

@@ -295,17 +295,7 @@ int main(int argc, char *argv[])
    //     element solution.
    a.RecoverFEMSolution(X, b, u);
 
-   // 13. Build a mass matrix to help solve for n.Grad(u) where 'n' is a surface
-   //     normal.
-   BilinearForm m(&fespace);
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-
-   ess_tdof_list.SetSize(0);
-   OperatorPtr M;
-   m.FormSystemMatrix(ess_tdof_list, M);
-
-   // 14. Compute the various boundary integrals.
+   // 13. Compute the various boundary integrals.
    mfem::out << endl
              << "Verifying boundary conditions" << endl
              << "=============================" << endl;
@@ -361,7 +351,7 @@ int main(int argc, char *argv[])
                 << " error " << err << endl;
    }
 
-   // 15. Save the refined mesh and the solution. This output can be viewed
+   // 14. Save the refined mesh and the solution. This output can be viewed
    //     later using GLVis: "glvis -m refined.mesh -g sol.gf".
    {
       ofstream mesh_ofs("refined.mesh");
@@ -372,7 +362,7 @@ int main(int argc, char *argv[])
       u.Save(sol_ofs);
    }
 
-   // 16. Send the solution by socket to a GLVis server.
+   // 15. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       string title_str = h1 ? "H1" : "DG";
@@ -385,7 +375,7 @@ int main(int argc, char *argv[])
                << " keys 'mmc'" << flush;
    }
 
-   // 17. Free the used memory.
+   // 16. Free the used memory.
    delete fec;
    delete mesh;
 

examples/ex27p.cpp

@@ -314,17 +314,7 @@ int main(int argc, char *argv[])
    //     local finite element solution on each processor.
    a.RecoverFEMSolution(X, b, u);
 
-   // 14. Build a mass matrix to help solve for n.Grad(u) where 'n' is a surface
-   //     normal.
-   ParBilinearForm m(&fespace);
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-
-   ess_tdof_list.SetSize(0);
-   OperatorPtr M;
-   m.FormSystemMatrix(ess_tdof_list, M);
-
-   // 15. Compute the various boundary integrals.
+   // 14. Compute the various boundary integrals.
    mfem::out << endl
              << "Verifying boundary conditions" << endl
              << "=============================" << endl;
@@ -380,7 +370,7 @@ int main(int argc, char *argv[])
                 << " error " << err << endl;
    }
 
-   // 16. Save the refined mesh and the solution in parallel. This output can be
+   // 15. Save the refined mesh and the solution in parallel. This output can be
    //     viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
    {
       ostringstream mesh_name, sol_name;
@@ -396,7 +386,7 @@ int main(int argc, char *argv[])
       u.Save(sol_ofs);
    }
 
-   // 17. Send the solution by socket to a GLVis server.
+   // 16. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       string title_str = h1 ? "H1" : "DG";
@@ -411,7 +401,7 @@ int main(int argc, char *argv[])
                << " keys 'mmc'" << flush;
    }
 
-   // 18. Free the used memory.
+   // 17. Free the used memory.
    delete fec;
 
    return 0;

fem/bilininteg.cpp

@@ -747,6 +747,7 @@ void DiffusionIntegrator::AssembleElementMatrix
 #ifdef MFEM_THREAD_SAFE
    DenseMatrix dshape(nd, dim), dshapedxt(nd, spaceDim);
    DenseMatrix dshapedxt_m(nd, MQ ? spaceDim : 0);
+   DenseMatrix M(MQ ? spaceDim : 0);
    Vector D(VQ ? VQ->GetVDim() : 0);
 #else
    dshape.SetSize(nd, dim);

fem/bilininteg_convection_ea.cpp

@@ -17,14 +17,14 @@ namespace mfem
 {
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EAConvectionAssemble1D(const int NE,
-                            const Array<double> &b,
-                            const Array<double> &g,
-                            const Vector &padata,
-                            Vector &eadata,
-                            const bool add,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void EAConvectionAssemble1D(const int NE,
+                                   const Array<double> &b,
+                                   const Array<double> &g,
+                                   const Vector &padata,
+                                   Vector &eadata,
+                                   const bool add,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -69,14 +69,14 @@ void EAConvectionAssemble1D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EAConvectionAssemble2D(const int NE,
-                            const Array<double> &b,
-                            const Array<double> &g,
-                            const Vector &padata,
-                            Vector &eadata,
-                            const bool add,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void EAConvectionAssemble2D(const int NE,
+                                   const Array<double> &b,
+                                   const Array<double> &g,
+                                   const Vector &padata,
+                                   Vector &eadata,
+                                   const bool add,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -146,14 +146,14 @@ void EAConvectionAssemble2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EAConvectionAssemble3D(const int NE,
-                            const Array<double> &b,
-                            const Array<double> &g,
-                            const Vector &padata,
-                            Vector &eadata,
-                            const bool add,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void EAConvectionAssemble3D(const int NE,
+                                   const Array<double> &b,
+                                   const Array<double> &g,
+                                   const Vector &padata,
+                                   Vector &eadata,
+                                   const bool add,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_convection_pa.cpp

@@ -21,13 +21,13 @@ namespace mfem
 // PA Convection Integrator
 
 // PA Convection Assemble 2D kernel
-void PAConvectionSetup2D(const int NQ,
-                         const int NE,
-                         const Array<double> &w,
-                         const Vector &j,
-                         const Vector &vel,
-                         const double alpha,
-                         Vector &op)
+static void PAConvectionSetup2D(const int NQ,
+                                const int NE,
+                                const Array<double> &w,
+                                const Vector &j,
+                                const Vector &vel,
+                                const double alpha,
+                                Vector &op)
 {
    constexpr int DIM = 2;
 
@@ -60,13 +60,13 @@ void PAConvectionSetup2D(const int NQ,
 }
 
 // PA Convection Assemble 3D kernel
-void PAConvectionSetup3D(const int NQ,
-                         const int NE,
-                         const Array<double> &w,
-                         const Vector &j,
-                         const Vector &vel,
-                         const double alpha,
-                         Vector &op)
+static void PAConvectionSetup3D(const int NQ,
+                                const int NE,
+                                const Array<double> &w,
+                                const Vector &j,
+                                const Vector &vel,
+                                const double alpha,
+                                Vector &op)
 {
    constexpr int DIM = 3;
    constexpr int SDIM = DIM;
@@ -135,7 +135,7 @@ static void PAConvectionSetup(const int dim,
 }
 
 // PA Convection Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
+template<int T_D1D = 0, int T_Q1D = 0> static
 void PAConvectionApply2D(const int ne,
                          const Array<double> &b,
                          const Array<double> &g,
@@ -254,7 +254,7 @@ void PAConvectionApply2D(const int ne,
 }
 
 // Optimized PA Convection Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
+template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
 void SmemPAConvectionApply2D(const int ne,
                              const Array<double> &b,
                              const Array<double> &g,
@@ -382,7 +382,7 @@ void SmemPAConvectionApply2D(const int ne,
 }
 
 // PA Convection Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
+template<int T_D1D = 0, int T_Q1D = 0> static
 void PAConvectionApply3D(const int ne,
                          const Array<double> &b,
                          const Array<double> &g,
@@ -563,7 +563,7 @@ void PAConvectionApply3D(const int ne,
 }
 
 // Optimized PA Convection Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
+template<int T_D1D = 0, int T_Q1D = 0> static
 void SmemPAConvectionApply3D(const int ne,
                              const Array<double> &b,
                              const Array<double> &g,

fem/bilininteg_dgtrace_ea.cpp

@@ -16,12 +16,12 @@
 namespace mfem
 {
 
-void EADGTraceAssemble1DInt(const int NF,
-                            const Array<double> &basis,
-                            const Vector &padata,
-                            Vector &eadata_int,
-                            Vector &eadata_ext,
-                            const bool add)
+static void EADGTraceAssemble1DInt(const int NF,
+                                   const Array<double> &basis,
+                                   const Vector &padata,
+                                   Vector &eadata_int,
+                                   Vector &eadata_ext,
+                                   const bool add)
 {
    auto D = Reshape(padata.Read(), 2, 2, NF);
    auto A_int = Reshape(eadata_int.ReadWrite(), 2, NF);
@@ -50,11 +50,11 @@ void EADGTraceAssemble1DInt(const int NF,
    });
 }
 
-void EADGTraceAssemble1DBdr(const int NF,
-                            const Array<double> &basis,
-                            const Vector &padata,
-                            Vector &eadata_bdr,
-                            const bool add)
+static void EADGTraceAssemble1DBdr(const int NF,
+                                   const Array<double> &basis,
+                                   const Vector &padata,
+                                   Vector &eadata_bdr,
+                                   const bool add)
 {
    auto D = Reshape(padata.Read(), 2, 2, NF);
    auto A_bdr = Reshape(eadata_bdr.ReadWrite(), NF);
@@ -72,14 +72,14 @@ void EADGTraceAssemble1DBdr(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EADGTraceAssemble2DInt(const int NF,
-                            const Array<double> &basis,
-                            const Vector &padata,
-                            Vector &eadata_int,
-                            Vector &eadata_ext,
-                            const bool add,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void EADGTraceAssemble2DInt(const int NF,
+                                   const Array<double> &basis,
+                                   const Vector &padata,
+                                   Vector &eadata_int,
+                                   Vector &eadata_ext,
+                                   const bool add,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -128,13 +128,13 @@ void EADGTraceAssemble2DInt(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EADGTraceAssemble2DBdr(const int NF,
-                            const Array<double> &basis,
-                            const Vector &padata,
-                            Vector &eadata_bdr,
-                            const bool add,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void EADGTraceAssemble2DBdr(const int NF,
+                                   const Array<double> &basis,
+                                   const Vector &padata,
+                                   Vector &eadata_bdr,
+                                   const bool add,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -170,14 +170,14 @@ void EADGTraceAssemble2DBdr(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EADGTraceAssemble3DInt(const int NF,
-                            const Array<double> &basis,
-                            const Vector &padata,
-                            Vector &eadata_int,
-                            Vector &eadata_ext,
-                            const bool add,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void EADGTraceAssemble3DInt(const int NF,
+                                   const Array<double> &basis,
+                                   const Vector &padata,
+                                   Vector &eadata_int,
+                                   Vector &eadata_ext,
+                                   const bool add,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -268,13 +268,13 @@ void EADGTraceAssemble3DInt(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EADGTraceAssemble3DBdr(const int NF,
-                            const Array<double> &basis,
-                            const Vector &padata,
-                            Vector &eadata_bdr,
-                            const bool add,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void EADGTraceAssemble3DBdr(const int NF,
+                                   const Array<double> &basis,
+                                   const Vector &padata,
+                                   Vector &eadata_bdr,
+                                   const bool add,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_dgtrace_pa.cpp

@@ -19,16 +19,16 @@ using namespace std;
 namespace mfem
 {
 // PA DG Trace Integrator
-void PADGTraceSetup2D(const int Q1D,
-                      const int NF,
-                      const Array<double> &w,
-                      const Vector &det,
-                      const Vector &nor,
-                      const Vector &rho,
-                      const Vector &vel,
-                      const double alpha,
-                      const double beta,
-                      Vector &op)
+static void PADGTraceSetup2D(const int Q1D,
+                             const int NF,
+                             const Array<double> &w,
+                             const Vector &det,
+                             const Vector &nor,
+                             const Vector &rho,
+                             const Vector &vel,
+                             const double alpha,
+                             const double beta,
+                             Vector &op)
 {
    const int VDIM = 2;
 
@@ -61,16 +61,16 @@ void PADGTraceSetup2D(const int Q1D,
    });
 }
 
-void PADGTraceSetup3D(const int Q1D,
-                      const int NF,
-                      const Array<double> &w,
-                      const Vector &det,
-                      const Vector &nor,
-                      const Vector &rho,
-                      const Vector &vel,
-                      const double alpha,
-                      const double beta,
-                      Vector &op)
+static void PADGTraceSetup3D(const int Q1D,
+                             const int NF,
+                             const Array<double> &w,
+                             const Vector &det,
+                             const Vector &nor,
+                             const Vector &rho,
+                             const Vector &vel,
+                             const double alpha,
+                             const double beta,
+                             Vector &op)
 {
    const int VDIM = 3;
 
@@ -301,7 +301,7 @@ void DGTraceIntegrator::AssemblePABoundaryFaces(const FiniteElementSpace& fes)
 }
 
 // PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0>
+template<int T_D1D = 0, int T_Q1D = 0> static
 void PADGTraceApply2D(const int NF,
                       const Array<double> &b,
                       const Array<double> &bt,
@@ -392,7 +392,7 @@ void PADGTraceApply2D(const int NF,
 }
 
 // PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0>
+template<int T_D1D = 0, int T_Q1D = 0> static
 void PADGTraceApply3D(const int NF,
                       const Array<double> &b,
                       const Array<double> &bt,
@@ -537,7 +537,7 @@ void PADGTraceApply3D(const int NF,
 }
 
 // Optimized PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
+template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
 void SmemPADGTraceApply3D(const int NF,
                           const Array<double> &b,
                           const Array<double> &bt,
@@ -701,7 +701,7 @@ static void PADGTraceApply(const int dim,
 }
 
 // PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0>
+template<int T_D1D = 0, int T_Q1D = 0> static
 void PADGTraceApplyTranspose2D(const int NF,
                                const Array<double> &b,
                                const Array<double> &bt,
@@ -797,7 +797,7 @@ void PADGTraceApplyTranspose2D(const int NF,
 }
 
 // PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0>
+template<int T_D1D = 0, int T_Q1D = 0> static
 void PADGTraceApplyTranspose3D(const int NF,
                                const Array<double> &b,
                                const Array<double> &bt,
@@ -953,7 +953,7 @@ void PADGTraceApplyTranspose3D(const int NF,
 }
 
 // Optimized PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
+template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
 void SmemPADGTraceApplyTranspose3D(const int NF,
                                    const Array<double> &b,
                                    const Array<double> &bt,

fem/bilininteg_diffusion_ea.cpp

@@ -17,14 +17,14 @@ namespace mfem
 {
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EADiffusionAssemble1D(const int NE,
-                           const Array<double> &b,
-                           const Array<double> &g,
-                           const Vector &padata,
-                           Vector &eadata,
-                           const bool add,
-                           const int d1d = 0,
-                           const int q1d = 0)
+static void EADiffusionAssemble1D(const int NE,
+                                  const Array<double> &b,
+                                  const Array<double> &g,
+                                  const Vector &padata,
+                                  Vector &eadata,
+                                  const bool add,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -68,14 +68,14 @@ void EADiffusionAssemble1D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EADiffusionAssemble2D(const int NE,
-                           const Array<double> &b,
-                           const Array<double> &g,
-                           const Vector &padata,
-                           Vector &eadata,
-                           const bool add,
-                           const int d1d = 0,
-                           const int q1d = 0)
+static void EADiffusionAssemble2D(const int NE,
+                                  const Array<double> &b,
+                                  const Array<double> &g,
+                                  const Vector &padata,
+                                  Vector &eadata,
+                                  const bool add,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -145,14 +145,14 @@ void EADiffusionAssemble2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EADiffusionAssemble3D(const int NE,
-                           const Array<double> &b,
-                           const Array<double> &g,
-                           const Vector &padata,
-                           Vector &eadata,
-                           const bool add,
-                           const int d1d = 0,
-                           const int q1d = 0)
+static void EADiffusionAssemble3D(const int NE,
+                                  const Array<double> &b,
+                                  const Array<double> &g,
+                                  const Vector &padata,
+                                  Vector &eadata,
+                                  const bool add,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_diffusion_pa.cpp

@@ -125,7 +125,7 @@ void PADiffusionSetup2D<2>(const int Q1D,
                D(qx,qy,0,e) = w_detJ * ( J22*R11 - J12*R21); // 1,1
                D(qx,qy,1,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1
                D(qx,qy,2,e) = w_detJ * (symmetric ? (-J21*R12 + J11*R22) :
-               (J22*R12 - J12*R22)); // 2,2 or 1,2
+                                        (J22*R12 - J12*R22)); // 2,2 or 1,2
                if (!symmetric)
                {
                   D(qx,qy,3,e) = w_detJ * (-J21*R12 + J11*R22); // 2,2
@@ -496,14 +496,14 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PADiffusionDiagonal2D(const int NE,
-                           const bool symmetric,
-                           const Array<double> &b,
-                           const Array<double> &g,
-                           const Vector &d,
-                           Vector &y,
-                           const int d1d = 0,
-                           const int q1d = 0)
+static void PADiffusionDiagonal2D(const int NE,
+                                  const bool symmetric,
+                                  const Array<double> &b,
+                                  const Array<double> &g,
+                                  const Vector &d,
+                                  Vector &y,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -562,14 +562,14 @@ void PADiffusionDiagonal2D(const int NE,
 
 // Shared memory PA Diffusion Diagonal 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-void SmemPADiffusionDiagonal2D(const int NE,
-                               const bool symmetric,
-                               const Array<double> &b_,
-                               const Array<double> &g_,
-                               const Vector &d_,
-                               Vector &y_,
-                               const int d1d = 0,
-                               const int q1d = 0)
+static void SmemPADiffusionDiagonal2D(const int NE,
+                                      const bool symmetric,
+                                      const Array<double> &b_,
+                                      const Array<double> &g_,
+                                      const Vector &d_,
+                                      Vector &y_,
+                                      const int d1d = 0,
+                                      const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -656,14 +656,14 @@ void SmemPADiffusionDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PADiffusionDiagonal3D(const int NE,
-                           const bool symmetric,
-                           const Array<double> &b,
-                           const Array<double> &g,
-                           const Vector &d,
-                           Vector &y,
-                           const int d1d = 0,
-                           const int q1d = 0)
+static void PADiffusionDiagonal3D(const int NE,
+                                  const bool symmetric,
+                                  const Array<double> &b,
+                                  const Array<double> &g,
+                                  const Vector &d,
+                                  Vector &y,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    constexpr int DIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -757,14 +757,14 @@ void PADiffusionDiagonal3D(const int NE,
 
 // Shared memory PA Diffusion Diagonal 3D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-void SmemPADiffusionDiagonal3D(const int NE,
-                               const bool symmetric,
-                               const Array<double> &b_,
-                               const Array<double> &g_,
-                               const Vector &d_,
-                               Vector &y_,
-                               const int d1d = 0,
-                               const int q1d = 0)
+static void SmemPADiffusionDiagonal3D(const int NE,
+                                      const bool symmetric,
+                                      const Array<double> &b_,
+                                      const Array<double> &g_,
+                                      const Vector &d_,
+                                      Vector &y_,
+                                      const int d1d = 0,
+                                      const int q1d = 0)
 {
    constexpr int DIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -1034,17 +1034,17 @@ static void OccaPADiffusionApply3D(const int D1D,
 
 // PA Diffusion Apply 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-void PADiffusionApply2D(const int NE,
-                        const bool symmetric,
-                        const Array<double> &b_,
-                        const Array<double> &g_,
-                        const Array<double> &bt_,
-                        const Array<double> &gt_,
-                        const Vector &d_,
-                        const Vector &x_,
-                        Vector &y_,
-                        const int d1d = 0,
-                        const int q1d = 0)
+static void PADiffusionApply2D(const int NE,
+                               const bool symmetric,
+                               const Array<double> &b_,
+                               const Array<double> &g_,
+                               const Array<double> &bt_,
+                               const Array<double> &gt_,
+                               const Vector &d_,
+                               const Vector &x_,
+                               Vector &y_,
+                               const int d1d = 0,
+                               const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1156,15 +1156,15 @@ void PADiffusionApply2D(const int NE,
 
 // Shared memory PA Diffusion Apply 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-void SmemPADiffusionApply2D(const int NE,
-                            const bool symmetric,
-                            const Array<double> &b_,
-                            const Array<double> &g_,
-                            const Vector &d_,
-                            const Vector &x_,
-                            Vector &y_,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void SmemPADiffusionApply2D(const int NE,
+                                   const bool symmetric,
+                                   const Array<double> &b_,
+                                   const Array<double> &g_,
+                                   const Vector &d_,
+                                   const Vector &x_,
+                                   Vector &y_,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1314,16 +1314,16 @@ void SmemPADiffusionApply2D(const int NE,
 
 // PA Diffusion Apply 3D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-void PADiffusionApply3D(const int NE,
-                        const bool symmetric,
-                        const Array<double> &b,
-                        const Array<double> &g,
-                        const Array<double> &bt,
-                        const Array<double> &gt,
-                        const Vector &d_,
-                        const Vector &x_,
-                        Vector &y_,
-                        int d1d = 0, int q1d = 0)
+static void PADiffusionApply3D(const int NE,
+                               const bool symmetric,
+                               const Array<double> &b,
+                               const Array<double> &g,
+                               const Array<double> &bt,
+                               const Array<double> &gt,
+                               const Vector &d_,
+                               const Vector &x_,
+                               Vector &y_,
+                               int d1d = 0, int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1533,15 +1533,15 @@ static MFEM_HOST_DEVICE inline double sign(const int q, const int d)
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void SmemPADiffusionApply3D(const int NE,
-                            const bool symmetric,
-                            const Array<double> &b_,
-                            const Array<double> &g_,
-                            const Vector &d_,
-                            const Vector &x_,
-                            Vector &y_,
-                            const int d1d = 0,
-                            const int q1d = 0)
+static void SmemPADiffusionApply3D(const int NE,
+                                   const bool symmetric,
+                                   const Array<double> &b_,
+                                   const Array<double> &g_,
+                                   const Vector &d_,
+                                   const Vector &x_,
+                                   Vector &y_,
+                                   const int d1d = 0,
+                                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_divergence.cpp

@@ -21,12 +21,12 @@ namespace mfem
 // PA Divergence Integrator
 
 // PA Divergence Assemble 2D kernel
-void PADivergenceSetup2D(const int Q1D,
-                         const int NE,
-                         const Array<double> &w,
-                         const Vector &j,
-                         const double COEFF,
-                         Vector &op)
+static void PADivergenceSetup2D(const int Q1D,
+                                const int NE,
+                                const Array<double> &w,
+                                const Vector &j,
+                                const double COEFF,
+                                Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -51,12 +51,12 @@ void PADivergenceSetup2D(const int Q1D,
 }
 
 // PA Divergence Assemble 3D kernel
-void PADivergenceSetup3D(const int Q1D,
-                         const int NE,
-                         const Array<double> &w,
-                         const Vector &j,
-                         const double COEFF,
-                         Vector &op)
+static void PADivergenceSetup3D(const int Q1D,
+                                const int NE,
+                                const Array<double> &w,
+                                const Vector &j,
+                                const double COEFF,
+                                Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -160,16 +160,16 @@ void VectorDivergenceIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
 // PA Divergence Apply 2D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-void PADivergenceApply2D(const int NE,
-                         const Array<double> &b,
-                         const Array<double> &g,
-                         const Array<double> &bt,
-                         const Vector &op_,
-                         const Vector &x_,
-                         Vector &y_,
-                         const int tr_d1d = 0,
-                         const int te_d1d = 0,
-                         const int q1d = 0)
+static void PADivergenceApply2D(const int NE,
+                                const Array<double> &b,
+                                const Array<double> &g,
+                                const Array<double> &bt,
+                                const Vector &op_,
+                                const Vector &x_,
+                                Vector &y_,
+                                const int tr_d1d = 0,
+                                const int te_d1d = 0,
+                                const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -281,16 +281,16 @@ void PADivergenceApply2D(const int NE,
 // Shared memory PA Divergence Apply 2D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0,
          const int T_NBZ = 0>
-void SmemPADivergenceApply2D(const int NE,
-                             const Array<double> &b_,
-                             const Array<double> &g_,
-                             const Array<double> &bt_,
-                             const Vector &op_,
-                             const Vector &x_,
-                             Vector &y_,
-                             const int tr_d1d = 0,
-                             const int te_d1d = 0,
-                             const int q1d = 0)
+static void SmemPADivergenceApply2D(const int NE,
+                                    const Array<double> &b_,
+                                    const Array<double> &g_,
+                                    const Array<double> &bt_,
+                                    const Vector &op_,
+                                    const Vector &x_,
+                                    Vector &y_,
+                                    const int tr_d1d = 0,
+                                    const int te_d1d = 0,
+                                    const int q1d = 0)
 {
    // TODO
    MFEM_ASSERT(false, "SHARED MEM NOT PROGRAMMED YET");
@@ -298,16 +298,16 @@ void SmemPADivergenceApply2D(const int NE,
 
 // PA Divergence Apply 2D kernel transpose
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-void PADivergenceApplyTranspose2D(const int NE,
-                                  const Array<double> &bt,
-                                  const Array<double> &gt,
-                                  const Array<double> &b,
-                                  const Vector &op_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int tr_d1d = 0,
-                                  const int te_d1d = 0,
-                                  const int q1d = 0)
+static void PADivergenceApplyTranspose2D(const int NE,
+                                         const Array<double> &bt,
+                                         const Array<double> &gt,
+                                         const Array<double> &b,
+                                         const Vector &op_,
+                                         const Vector &x_,
+                                         Vector &y_,
+                                         const int tr_d1d = 0,
+                                         const int te_d1d = 0,
+                                         const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -414,16 +414,16 @@ void PADivergenceApplyTranspose2D(const int NE,
 
 // PA Vector Divergence Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-void PADivergenceApply3D(const int NE,
-                         const Array<double> &b,
-                         const Array<double> &g,
-                         const Array<double> &bt,
-                         const Vector &op_,
-                         const Vector &x_,
-                         Vector &y_,
-                         int tr_d1d = 0,
-                         int te_d1d = 0,
-                         int q1d = 0)
+static void PADivergenceApply3D(const int NE,
+                                const Array<double> &b,
+                                const Array<double> &g,
+                                const Array<double> &bt,
+                                const Vector &op_,
+                                const Vector &x_,
+                                Vector &y_,
+                                int tr_d1d = 0,
+                                int te_d1d = 0,
+                                int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -597,16 +597,16 @@ void PADivergenceApply3D(const int NE,
 
 // PA Vector Divergence Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-void PADivergenceApplyTranspose3D(const int NE,
-                                  const Array<double> &bt,
-                                  const Array<double> &gt,
-                                  const Array<double> &b,
-                                  const Vector &op_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  int tr_d1d = 0,
-                                  int te_d1d = 0,
-                                  int q1d = 0)
+static void PADivergenceApplyTranspose3D(const int NE,
+                                         const Array<double> &bt,
+                                         const Array<double> &gt,
+                                         const Array<double> &b,
+                                         const Vector &op_,
+                                         const Vector &x_,
+                                         Vector &y_,
+                                         int tr_d1d = 0,
+                                         int te_d1d = 0,
+                                         int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -775,16 +775,16 @@ void PADivergenceApplyTranspose3D(const int NE,
 
 // Shared memory PA Vector Divergence Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-void SmemPADivergenceApply3D(const int NE,
-                             const Array<double> &b_,
-                             const Array<double> &g_,
-                             const Array<double> &bt_,
-                             const Vector &q_,
-                             const Vector &x_,
-                             Vector &y_,
-                             const int tr_d1d = 0,
-                             const int te_d1d = 0,
-                             const int q1d = 0)
+static void SmemPADivergenceApply3D(const int NE,
+                                    const Array<double> &b_,
+                                    const Array<double> &g_,
+                                    const Array<double> &bt_,
+                                    const Vector &q_,
+                                    const Vector &x_,
+                                    Vector &y_,
+                                    const int tr_d1d = 0,
+                                    const int te_d1d = 0,
+                                    const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;

fem/bilininteg_gradient.cpp

@@ -70,12 +70,12 @@ namespace mfem
    the \b MFEM_SHARED keyword for local arrays. */
 
 // PA Gradient Assemble 2D kernel
-void PAGradientSetup2D(const int Q1D,
-                       const int NE,
-                       const Array<double> &w,
-                       const Vector &j,
-                       const Vector &c,
-                       Vector &op)
+static void PAGradientSetup2D(const int Q1D,
+                              const int NE,
+                              const Array<double> &w,
+                              const Vector &j,
+                              const Vector &c,
+                              Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -105,12 +105,12 @@ void PAGradientSetup2D(const int Q1D,
 }
 
 // PA Gradient Assemble 3D kernel
-void PAGradientSetup3D(const int Q1D,
-                       const int NE,
-                       const Array<double> &w,
-                       const Vector &j,
-                       const Vector &c,
-                       Vector &op)
+static void PAGradientSetup3D(const int Q1D,
+                              const int NE,
+                              const Array<double> &w,
+                              const Vector &j,
+                              const Vector &c,
+                              Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -254,16 +254,16 @@ void GradientIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
 // PA Gradient Apply 2D kernel
 template<int T_TR_D1D = 0, int T_TE_D1D = 0, int T_Q1D = 0>
-void PAGradientApply2D(const int NE,
-                       const Array<double> &b,
-                       const Array<double> &g,
-                       const Array<double> &bt,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_,
-                       const int tr_d1d = 0,
-                       const int te_d1d = 0,
-                       const int q1d = 0)
+static void PAGradientApply2D(const int NE,
+                              const Array<double> &b,
+                              const Array<double> &g,
+                              const Array<double> &bt,
+                              const Vector &op_,
+                              const Vector &x_,
+                              Vector &y_,
+                              const int tr_d1d = 0,
+                              const int te_d1d = 0,
+                              const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -384,16 +384,16 @@ static void PAGradientApplyTranspose2D(const int NE,
 
 // PA Gradient Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-void PAGradientApply3D(const int NE,
-                       const Array<double> &b,
-                       const Array<double> &g,
-                       const Array<double> &bt,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_,
-                       int tr_d1d = 0,
-                       int te_d1d = 0,
-                       int q1d = 0)
+static void PAGradientApply3D(const int NE,
+                              const Array<double> &b,
+                              const Array<double> &g,
+                              const Array<double> &bt,
+                              const Vector &op_,
+                              const Vector &x_,
+                              Vector &y_,
+                              int tr_d1d = 0,
+                              int te_d1d = 0,
+                              int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -579,16 +579,16 @@ static void PAGradientApplyTranspose3D(const int NE,
 
 // Shared memory PA Gradient Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-void SmemPAGradientApply3D(const int NE,
-                           const Array<double> &b_,
-                           const Array<double> &g_,
-                           const Array<double> &bt_,
-                           const Vector &d_,
-                           const Vector &x_,
-                           Vector &y_,
-                           const int tr_d1d = 0,
-                           const int te_d1d = 0,
-                           const int q1d = 0)
+static void SmemPAGradientApply3D(const int NE,
+                                  const Array<double> &b_,
+                                  const Array<double> &g_,
+                                  const Array<double> &bt_,
+                                  const Vector &d_,
+                                  const Vector &x_,
+                                  Vector &y_,
+                                  const int tr_d1d = 0,
+                                  const int te_d1d = 0,
+                                  const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;

fem/bilininteg_hcurl.cpp

@@ -186,7 +186,7 @@ void PAHcurlMassAssembleDiagonal2D(const int D1D,
                   const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
 
                   mass[qx] += wy * wy * ((c == 0) ? op(qx,qy,0,e) :
-                  op(qx,qy,symmetric ? 2 : 3, e));
+                                         op(qx,qy,symmetric ? 2 : 3, e));
                }
             }
 
@@ -237,7 +237,7 @@ void PAHcurlMassAssembleDiagonal3D(const int D1D,
          const int D1Dx = (c == 0) ? D1D - 1 : D1D;
 
          const int opc = (c == 0) ? 0 : ((c == 1) ? (symmetric ? 3 : 4) :
-         (symmetric ? 5 : 8));
+                                         (symmetric ? 5 : 8));
 
          double mass[MAX_Q1D];
 
@@ -791,12 +791,12 @@ void SmemPAHcurlMassApply3D(const int D1D,
 }
 
 // PA H(curl) curl-curl assemble 2D kernel
-void PACurlCurlSetup2D(const int Q1D,
-                       const int NE,
-                       const Array<double> &w,
-                       const Vector &j,
-                       Vector &coeff,
-                       Vector &op)
+static void PACurlCurlSetup2D(const int Q1D,
+                              const int NE,
+                              const Array<double> &w,
+                              const Vector &j,
+                              Vector &coeff,
+                              Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -818,13 +818,13 @@ void PACurlCurlSetup2D(const int Q1D,
 }
 
 // PA H(curl) curl-curl assemble 3D kernel
-void PACurlCurlSetup3D(const int Q1D,
-                       const int coeffDim,
-                       const int NE,
-                       const Array<double> &w,
-                       const Vector &j,
-                       Vector &coeff,
-                       Vector &op)
+static void PACurlCurlSetup3D(const int Q1D,
+                              const int coeffDim,
+                              const int NE,
+                              const Array<double> &w,
+                              const Vector &j,
+                              Vector &coeff,
+                              Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    const bool symmetric = (coeffDim != 9);
@@ -1045,16 +1045,16 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
    }
 }
 
-void PACurlCurlApply2D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const Array<double> &bo,
-                       const Array<double> &bot,
-                       const Array<double> &gc,
-                       const Array<double> &gct,
-                       const Vector &pa_data,
-                       const Vector &x,
-                       Vector &y)
+static void PACurlCurlApply2D(const int D1D,
+                              const int Q1D,
+                              const int NE,
+                              const Array<double> &bo,
+                              const Array<double> &bot,
+                              const Array<double> &gc,
+                              const Array<double> &gct,
+                              const Vector &pa_data,
+                              const Vector &x,
+                              Vector &y)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HCURL_MAX_D1D;
@@ -1166,19 +1166,19 @@ void PACurlCurlApply2D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void PACurlCurlApply3D(const int D1D,
-                       const int Q1D,
-                       const bool symmetric,
-                       const int NE,
-                       const Array<double> &bo,
-                       const Array<double> &bc,
-                       const Array<double> &bot,
-                       const Array<double> &bct,
-                       const Array<double> &gc,
-                       const Array<double> &gct,
-                       const Vector &pa_data,
-                       const Vector &x,
-                       Vector &y)
+static void PACurlCurlApply3D(const int D1D,
+                              const int Q1D,
+                              const bool symmetric,
+                              const int NE,
+                              const Array<double> &bo,
+                              const Array<double> &bc,
+                              const Array<double> &bot,
+                              const Array<double> &bct,
+                              const Array<double> &gc,
+                              const Array<double> &gct,
+                              const Vector &pa_data,
+                              const Vector &x,
+                              Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -1677,19 +1677,19 @@ void PACurlCurlApply3D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void SmemPACurlCurlApply3D(const int D1D,
-                           const int Q1D,
-                           const bool symmetric,
-                           const int NE,
-                           const Array<double> &bo,
-                           const Array<double> &bc,
-                           const Array<double> &bot,
-                           const Array<double> &bct,
-                           const Array<double> &gc,
-                           const Array<double> &gct,
-                           const Vector &pa_data,
-                           const Vector &x,
-                           Vector &y)
+static void SmemPACurlCurlApply3D(const int D1D,
+                                  const int Q1D,
+                                  const bool symmetric,
+                                  const int NE,
+                                  const Array<double> &bo,
+                                  const Array<double> &bc,
+                                  const Array<double> &bot,
+                                  const Array<double> &bct,
+                                  const Array<double> &gc,
+                                  const Array<double> &gct,
+                                  const Vector &pa_data,
+                                  const Vector &x,
+                                  Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -2032,13 +2032,13 @@ void CurlCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
-void PACurlCurlAssembleDiagonal2D(const int D1D,
-                                  const int Q1D,
-                                  const int NE,
-                                  const Array<double> &bo,
-                                  const Array<double> &gc,
-                                  const Vector &pa_data,
-                                  Vector &diag)
+static void PACurlCurlAssembleDiagonal2D(const int D1D,
+                                         const int Q1D,
+                                         const int NE,
+                                         const Array<double> &bo,
+                                         const Array<double> &gc,
+                                         const Vector &pa_data,
+                                         Vector &diag)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
@@ -2087,16 +2087,16 @@ void PACurlCurlAssembleDiagonal2D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void PACurlCurlAssembleDiagonal3D(const int D1D,
-                                  const int Q1D,
-                                  const bool symmetric,
-                                  const int NE,
-                                  const Array<double> &bo,
-                                  const Array<double> &bc,
-                                  const Array<double> &go,
-                                  const Array<double> &gc,
-                                  const Vector &pa_data,
-                                  Vector &diag)
+static void PACurlCurlAssembleDiagonal3D(const int D1D,
+                                         const int Q1D,
+                                         const bool symmetric,
+                                         const int NE,
+                                         const Array<double> &bo,
+                                         const Array<double> &bc,
+                                         const Array<double> &go,
+                                         const Array<double> &gc,
+                                         const Vector &pa_data,
+                                         Vector &diag)
 {
    constexpr static int VDIM = 3;
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
@@ -2273,16 +2273,16 @@ void PACurlCurlAssembleDiagonal3D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void SmemPACurlCurlAssembleDiagonal3D(const int D1D,
-                                      const int Q1D,
-                                      const bool symmetric,
-                                      const int NE,
-                                      const Array<double> &bo,
-                                      const Array<double> &bc,
-                                      const Array<double> &go,
-                                      const Array<double> &gc,
-                                      const Vector &pa_data,
-                                      Vector &diag)
+static void SmemPACurlCurlAssembleDiagonal3D(const int D1D,
+                                             const int Q1D,
+                                             const bool symmetric,
+                                             const int NE,
+                                             const Array<double> &bo,
+                                             const Array<double> &bc,
+                                             const Array<double> &go,
+                                             const Array<double> &gc,
+                                             const Vector &pa_data,
+                                             Vector &diag)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -2955,18 +2955,18 @@ void MixedVectorCurlIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 // Apply to x corresponding to DOF's in H(curl) (trial), whose curl is
 // integrated against H(curl) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void PAHcurlL2Apply3D(const int D1D,
-                      const int Q1D,
-                      const int coeffDim,
-                      const int NE,
-                      const Array<double> &bo,
-                      const Array<double> &bc,
-                      const Array<double> &bot,
-                      const Array<double> &bct,
-                      const Array<double> &gc,
-                      const Vector &pa_data,
-                      const Vector &x,
-                      Vector &y)
+static void PAHcurlL2Apply3D(const int D1D,
+                             const int Q1D,
+                             const int coeffDim,
+                             const int NE,
+                             const Array<double> &bo,
+                             const Array<double> &bc,
+                             const Array<double> &bot,
+                             const Array<double> &bct,
+                             const Array<double> &gc,
+                             const Vector &pa_data,
+                             const Vector &x,
+                             Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -3297,16 +3297,16 @@ void PAHcurlL2Apply3D(const int D1D,
 // Apply to x corresponding to DOF's in H(curl) (trial), whose curl is
 // integrated against H(curl) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void SmemPAHcurlL2Apply3D(const int D1D,
-                          const int Q1D,
-                          const int coeffDim,
-                          const int NE,
-                          const Array<double> &bo,
-                          const Array<double> &bc,
-                          const Array<double> &gc,
-                          const Vector &pa_data,
-                          const Vector &x,
-                          Vector &y)
+static void SmemPAHcurlL2Apply3D(const int D1D,
+                                 const int Q1D,
+                                 const int coeffDim,
+                                 const int NE,
+                                 const Array<double> &bo,
+                                 const Array<double> &bc,
+                                 const Array<double> &gc,
+                                 const Vector &pa_data,
+                                 const Vector &x,
+                                 Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -3585,18 +3585,18 @@ void SmemPAHcurlL2Apply3D(const int D1D,
 // Apply to x corresponding to DOF's in H(curl) (trial), whose curl is
 // integrated against H(div) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void PAHcurlHdivApply3D(const int D1D,
-                        const int D1Dtest,
-                        const int Q1D,
-                        const int NE,
-                        const Array<double> &bo,
-                        const Array<double> &bc,
-                        const Array<double> &bot,
-                        const Array<double> &bct,
-                        const Array<double> &gc,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y)
+static void PAHcurlHdivApply3D(const int D1D,
+                               const int D1Dtest,
+                               const int Q1D,
+                               const int NE,
+                               const Array<double> &bo,
+                               const Array<double> &bc,
+                               const Array<double> &bot,
+                               const Array<double> &bct,
+                               const Array<double> &gc,
+                               const Vector &pa_data,
+                               const Vector &x,
+                               Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -4071,18 +4071,18 @@ void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
 // Apply to x corresponding to DOF's in H(curl) (trial), integrated against curl
 // of H(curl) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void PAHcurlL2Apply3DTranspose(const int D1D,
-                               const int Q1D,
-                               const int coeffDim,
-                               const int NE,
-                               const Array<double> &bo,
-                               const Array<double> &bc,
-                               const Array<double> &bot,
-                               const Array<double> &bct,
-                               const Array<double> &gct,
-                               const Vector &pa_data,
-                               const Vector &x,
-                               Vector &y)
+static void PAHcurlL2Apply3DTranspose(const int D1D,
+                                      const int Q1D,
+                                      const int coeffDim,
+                                      const int NE,
+                                      const Array<double> &bo,
+                                      const Array<double> &bc,
+                                      const Array<double> &bot,
+                                      const Array<double> &bct,
+                                      const Array<double> &gct,
+                                      const Vector &pa_data,
+                                      const Vector &x,
+                                      Vector &y)
 {
    // See PAHcurlL2Apply3D for comments.
 
@@ -4413,16 +4413,16 @@ void PAHcurlL2Apply3DTranspose(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-void SmemPAHcurlL2Apply3DTranspose(const int D1D,
-                                   const int Q1D,
-                                   const int coeffDim,
-                                   const int NE,
-                                   const Array<double> &bo,
-                                   const Array<double> &bc,
-                                   const Array<double> &gc,
-                                   const Vector &pa_data,
-                                   const Vector &x,
-                                   Vector &y)
+static void SmemPAHcurlL2Apply3DTranspose(const int D1D,
+                                          const int Q1D,
+                                          const int coeffDim,
+                                          const int NE,
+                                          const Array<double> &bo,
+                                          const Array<double> &bc,
+                                          const Array<double> &gc,
+                                          const Vector &pa_data,
+                                          const Vector &x,
+                                          Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -4675,13 +4675,13 @@ void MixedVectorWeakCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
 // Apply to x corresponding to DOFs in H^1 (domain) the (topological) gradient
 // to get a dof in H(curl) (range). You can think of the range as the "test" space
 // and the domain as the "trial" space, but there's no integration.
-void PAHcurlApplyGradient2D(const int c_dofs1D,
-                            const int o_dofs1D,
-                            const int NE,
-                            const Array<double> &B_,
-                            const Array<double> &G_,
-                            const Vector &x_,
-                            Vector &y_)
+static void PAHcurlApplyGradient2D(const int c_dofs1D,
+                                   const int o_dofs1D,
+                                   const int NE,
+                                   const Array<double> &B_,
+                                   const Array<double> &G_,
+                                   const Vector &x_,
+                                   Vector &y_)
 {
    auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
@@ -4753,12 +4753,12 @@ void PAHcurlApplyGradient2D(const int c_dofs1D,
 }
 
 // Specialization of PAHcurlApplyGradient2D to the case where B is identity
-void PAHcurlApplyGradient2DBId(const int c_dofs1D,
-                               const int o_dofs1D,
-                               const int NE,
-                               const Array<double> &G_,
-                               const Vector &x_,
-                               Vector &y_)
+static void PAHcurlApplyGradient2DBId(const int c_dofs1D,
+                                      const int o_dofs1D,
+                                      const int NE,
+                                      const Array<double> &G_,
+                                      const Vector &x_,
+                                      Vector &y_)
 {
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
 
@@ -4822,7 +4822,7 @@ void PAHcurlApplyGradient2DBId(const int c_dofs1D,
    });
 }
 
-void PAHcurlApplyGradientTranspose2D(
+static void PAHcurlApplyGradientTranspose2D(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &B_, const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -4898,7 +4898,7 @@ void PAHcurlApplyGradientTranspose2D(
 
 // Specialization of PAHcurlApplyGradientTranspose2D to the case where
 // B is identity
-void PAHcurlApplyGradientTranspose2DBId(
+static void PAHcurlApplyGradientTranspose2DBId(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -4965,13 +4965,13 @@ void PAHcurlApplyGradientTranspose2DBId(
    });
 }
 
-void PAHcurlApplyGradient3D(const int c_dofs1D,
-                            const int o_dofs1D,
-                            const int NE,
-                            const Array<double> &B_,
-                            const Array<double> &G_,
-                            const Vector &x_,
-                            Vector &y_)
+static void PAHcurlApplyGradient3D(const int c_dofs1D,
+                                   const int o_dofs1D,
+                                   const int NE,
+                                   const Array<double> &B_,
+                                   const Array<double> &G_,
+                                   const Vector &x_,
+                                   Vector &y_)
 {
    auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
@@ -5154,12 +5154,12 @@ void PAHcurlApplyGradient3D(const int c_dofs1D,
 }
 
 // Specialization of PAHcurlApplyGradient3D to the case where
-void PAHcurlApplyGradient3DBId(const int c_dofs1D,
-                               const int o_dofs1D,
-                               const int NE,
-                               const Array<double> &G_,
-                               const Vector &x_,
-                               Vector &y_)
+static void PAHcurlApplyGradient3DBId(const int c_dofs1D,
+                                      const int o_dofs1D,
+                                      const int NE,
+                                      const Array<double> &G_,
+                                      const Vector &x_,
+                                      Vector &y_)
 {
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
 
@@ -5322,7 +5322,7 @@ void PAHcurlApplyGradient3DBId(const int c_dofs1D,
    });
 }
 
-void PAHcurlApplyGradientTranspose3D(
+static void PAHcurlApplyGradientTranspose3D(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &B_, const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -5507,7 +5507,7 @@ void PAHcurlApplyGradientTranspose3D(
 }
 
 // Specialization of PAHcurlApplyGradientTranspose3D to the case where
-void PAHcurlApplyGradientTranspose3DBId(
+static void PAHcurlApplyGradientTranspose3DBId(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -5789,14 +5789,14 @@ void GradientInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
    }
 }
 
-void PAHcurlVecH1IdentityApply3D(const int c_dofs1D,
-                                 const int o_dofs1D,
-                                 const int NE,
-                                 const Array<double> &Bclosed,
-                                 const Array<double> &Bopen,
-                                 const Vector &pa_data,
-                                 const Vector &x_,
-                                 Vector &y_)
+static void PAHcurlVecH1IdentityApply3D(const int c_dofs1D,
+                                        const int o_dofs1D,
+                                        const int NE,
+                                        const Array<double> &Bclosed,
+                                        const Array<double> &Bopen,
+                                        const Vector &pa_data,
+                                        const Vector &x_,
+                                        Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
@@ -6002,14 +6002,14 @@ void PAHcurlVecH1IdentityApply3D(const int c_dofs1D,
    });
 }
 
-void PAHcurlVecH1IdentityApplyTranspose3D(const int c_dofs1D,
-                                          const int o_dofs1D,
-                                          const int NE,
-                                          const Array<double> &Bclosed,
-                                          const Array<double> &Bopen,
-                                          const Vector &pa_data,
-                                          const Vector &x_,
-                                          Vector &y_)
+static void PAHcurlVecH1IdentityApplyTranspose3D(const int c_dofs1D,
+                                                 const int o_dofs1D,
+                                                 const int NE,
+                                                 const Array<double> &Bclosed,
+                                                 const Array<double> &Bopen,
+                                                 const Vector &pa_data,
+                                                 const Vector &x_,
+                                                 Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
@@ -6228,14 +6228,14 @@ void PAHcurlVecH1IdentityApplyTranspose3D(const int c_dofs1D,
    });
 }
 
-void PAHcurlVecH1IdentityApply2D(const int c_dofs1D,
-                                 const int o_dofs1D,
-                                 const int NE,
-                                 const Array<double> &Bclosed,
-                                 const Array<double> &Bopen,
-                                 const Vector &pa_data,
-                                 const Vector &x_,
-                                 Vector &y_)
+static void PAHcurlVecH1IdentityApply2D(const int c_dofs1D,
+                                        const int o_dofs1D,
+                                        const int NE,
+                                        const Array<double> &Bclosed,
+                                        const Array<double> &Bopen,
+                                        const Vector &pa_data,
+                                        const Vector &x_,
+                                        Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
@@ -6327,14 +6327,14 @@ void PAHcurlVecH1IdentityApply2D(const int c_dofs1D,
    });
 }
 
-void PAHcurlVecH1IdentityApplyTranspose2D(const int c_dofs1D,
-                                          const int o_dofs1D,
-                                          const int NE,
-                                          const Array<double> &Bclosed,
-                                          const Array<double> &Bopen,
-                                          const Vector &pa_data,
-                                          const Vector &x_,
-                                          Vector &y_)
+static void PAHcurlVecH1IdentityApplyTranspose2D(const int c_dofs1D,
+                                                 const int o_dofs1D,
+                                                 const int NE,
+                                                 const Array<double> &Bclosed,
+                                                 const Array<double> &Bopen,
+                                                 const Vector &pa_data,
+                                                 const Vector &x_,
+                                                 Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);

fem/bilininteg_hdiv.cpp

@@ -539,12 +539,12 @@ void PAHdivMassApply3D(const int D1D,
 
 // PA H(div) div-div assemble 2D kernel
 // NOTE: this is identical to PACurlCurlSetup3D
-void PADivDivSetup2D(const int Q1D,
-                     const int NE,
-                     const Array<double> &w,
-                     const Vector &j,
-                     Vector &coeff_,
-                     Vector &op)
+static void PADivDivSetup2D(const int Q1D,
+                            const int NE,
+                            const Array<double> &w,
+                            const Vector &j,
+                            Vector &coeff_,
+                            Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -565,12 +565,12 @@ void PADivDivSetup2D(const int Q1D,
    });
 }
 
-void PADivDivSetup3D(const int Q1D,
-                     const int NE,
-                     const Array<double> &w,
-                     const Vector &j,
-                     Vector &coeff_,
-                     Vector &op)
+static void PADivDivSetup3D(const int Q1D,
+                            const int NE,
+                            const Array<double> &w,
+                            const Vector &j,
+                            Vector &coeff_,
+                            Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -599,16 +599,16 @@ void PADivDivSetup3D(const int Q1D,
    });
 }
 
-void PADivDivApply2D(const int D1D,
-                     const int Q1D,
-                     const int NE,
-                     const Array<double> &Bo_,
-                     const Array<double> &Gc_,
-                     const Array<double> &Bot_,
-                     const Array<double> &Gct_,
-                     const Vector &op_,
-                     const Vector &x_,
-                     Vector &y_)
+static void PADivDivApply2D(const int D1D,
+                            const int Q1D,
+                            const int NE,
+                            const Array<double> &Bo_,
+                            const Array<double> &Gc_,
+                            const Array<double> &Bot_,
+                            const Array<double> &Gct_,
+                            const Vector &op_,
+                            const Vector &x_,
+                            Vector &y_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HDIV_MAX_D1D;
@@ -718,16 +718,16 @@ void PADivDivApply2D(const int D1D,
    }); // end of element loop
 }
 
-void PADivDivApply3D(const int D1D,
-                     const int Q1D,
-                     const int NE,
-                     const Array<double> &Bo_,
-                     const Array<double> &Gc_,
-                     const Array<double> &Bot_,
-                     const Array<double> &Gct_,
-                     const Vector &op_,
-                     const Vector &x_,
-                     Vector &y_)
+static void PADivDivApply3D(const int D1D,
+                            const int Q1D,
+                            const int NE,
+                            const Array<double> &Bo_,
+                            const Array<double> &Gc_,
+                            const Array<double> &Bot_,
+                            const Array<double> &Gct_,
+                            const Vector &op_,
+                            const Vector &x_,
+                            Vector &y_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -967,13 +967,13 @@ void DivDivIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
-void PADivDivAssembleDiagonal2D(const int D1D,
-                                const int Q1D,
-                                const int NE,
-                                const Array<double> &Bo_,
-                                const Array<double> &Gc_,
-                                const Vector &op_,
-                                Vector &diag_)
+static void PADivDivAssembleDiagonal2D(const int D1D,
+                                       const int Q1D,
+                                       const int NE,
+                                       const Array<double> &Bo_,
+                                       const Array<double> &Gc_,
+                                       const Vector &op_,
+                                       Vector &diag_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_Q1D = HDIV_MAX_Q1D;
@@ -1023,13 +1023,13 @@ void PADivDivAssembleDiagonal2D(const int D1D,
    });
 }
 
-void PADivDivAssembleDiagonal3D(const int D1D,
-                                const int Q1D,
-                                const int NE,
-                                const Array<double> &Bo_,
-                                const Array<double> &Gc_,
-                                const Vector &op_,
-                                Vector &diag_)
+static void PADivDivAssembleDiagonal3D(const int D1D,
+                                       const int Q1D,
+                                       const int NE,
+                                       const Array<double> &Bo_,
+                                       const Array<double> &Gc_,
+                                       const Vector &op_,
+                                       Vector &diag_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1104,11 +1104,11 @@ void DivDivIntegrator::AssembleDiagonalPA(Vector& diag)
 }
 
 // PA H(div)-L2 (div u, p) assemble 2D kernel
-void PADivL2Setup2D(const int Q1D,
-                    const int NE,
-                    const Array<double> &w,
-                    Vector &coeff_,
-                    Vector &op)
+static void PADivL2Setup2D(const int Q1D,
+                           const int NE,
+                           const Array<double> &w,
+                           Vector &coeff_,
+                           Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -1123,11 +1123,11 @@ void PADivL2Setup2D(const int Q1D,
    });
 }
 
-void PADivL2Setup3D(const int Q1D,
-                    const int NE,
-                    const Array<double> &w,
-                    Vector &coeff_,
-                    Vector &op)
+static void PADivL2Setup3D(const int Q1D,
+                           const int NE,
+                           const Array<double> &w,
+                           Vector &coeff_,
+                           Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -1225,16 +1225,16 @@ VectorFEDivergenceIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
 // Apply to x corresponding to DOF's in H(div) (trial), whose divergence is
 // integrated against L_2 test functions corresponding to y.
-void PAHdivL2Apply3D(const int D1D,
-                     const int Q1D,
-                     const int L2D1D,
-                     const int NE,
-                     const Array<double> &Bo_,
-                     const Array<double> &Gc_,
-                     const Array<double> &L2Bot_,
-                     const Vector &op_,
-                     const Vector &x_,
-                     Vector &y_)
+static void PAHdivL2Apply3D(const int D1D,
+                            const int Q1D,
+                            const int L2D1D,
+                            const int NE,
+                            const Array<double> &Bo_,
+                            const Array<double> &Gc_,
+                            const Array<double> &L2Bot_,
+                            const Vector &op_,
+                            const Vector &x_,
+                            Vector &y_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1388,16 +1388,16 @@ void PAHdivL2Apply3D(const int D1D,
 
 // Apply to x corresponding to DOF's in H(div) (trial), whose divergence is
 // integrated against L_2 test functions corresponding to y.
-void PAHdivL2Apply2D(const int D1D,
-                     const int Q1D,
-                     const int L2D1D,
-                     const int NE,
-                     const Array<double> &Bo_,
-                     const Array<double> &Gc_,
-                     const Array<double> &L2Bot_,
-                     const Vector &op_,
-                     const Vector &x_,
-                     Vector &y_)
+static void PAHdivL2Apply2D(const int D1D,
+                            const int Q1D,
+                            const int L2D1D,
+                            const int NE,
+                            const Array<double> &Bo_,
+                            const Array<double> &Gc_,
+                            const Array<double> &L2Bot_,
+                            const Vector &op_,
+                            const Vector &x_,
+                            Vector &y_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HDIV_MAX_D1D;
@@ -1494,16 +1494,16 @@ void PAHdivL2Apply2D(const int D1D,
    }); // end of element loop
 }
 
-void PAHdivL2ApplyTranspose3D(const int D1D,
-                              const int Q1D,
-                              const int L2D1D,
-                              const int NE,
-                              const Array<double> &L2Bo_,
-                              const Array<double> &Gct_,
-                              const Array<double> &Bot_,
-                              const Vector &op_,
-                              const Vector &x_,
-                              Vector &y_)
+static void PAHdivL2ApplyTranspose3D(const int D1D,
+                                     const int Q1D,
+                                     const int L2D1D,
+                                     const int NE,
+                                     const Array<double> &L2Bo_,
+                                     const Array<double> &Gct_,
+                                     const Array<double> &Bot_,
+                                     const Vector &op_,
+                                     const Vector &x_,
+                                     Vector &y_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1656,16 +1656,16 @@ void PAHdivL2ApplyTranspose3D(const int D1D,
    }); // end of element loop
 }
 
-void PAHdivL2ApplyTranspose2D(const int D1D,
-                              const int Q1D,
-                              const int L2D1D,
-                              const int NE,
-                              const Array<double> &L2Bo_,
-                              const Array<double> &Gct_,
-                              const Array<double> &Bot_,
-                              const Vector &op_,
-                              const Vector &x_,
-                              Vector &y_)
+static void PAHdivL2ApplyTranspose2D(const int D1D,
+                                     const int Q1D,
+                                     const int L2D1D,
+                                     const int NE,
+                                     const Array<double> &L2Bo_,
+                                     const Array<double> &Gct_,
+                                     const Array<double> &Bot_,
+                                     const Vector &op_,
+                                     const Vector &x_,
+                                     Vector &y_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HDIV_MAX_D1D;
@@ -1791,16 +1791,16 @@ void VectorFEDivergenceIntegrator::AddMultTransposePA(const Vector &x,
    }
 }
 
-void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
-                                      const int Q1D,
-                                      const int L2D1D,
-                                      const int NE,
-                                      const Array<double> &L2Bo_,
-                                      const Array<double> &Gct_,
-                                      const Array<double> &Bot_,
-                                      const Vector &op_,
-                                      const Vector &D_,
-                                      Vector &diag_)
+static void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
+                                             const int Q1D,
+                                             const int L2D1D,
+                                             const int NE,
+                                             const Array<double> &L2Bo_,
+                                             const Array<double> &Gct_,
+                                             const Array<double> &Bot_,
+                                             const Vector &op_,
+                                             const Vector &D_,
+                                             Vector &diag_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1916,16 +1916,16 @@ void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
    }); // end of element loop
 }
 
-void PAHdivL2AssembleDiagonal_ADAt_2D(const int D1D,
-                                      const int Q1D,
-                                      const int L2D1D,
-                                      const int NE,
-                                      const Array<double> &L2Bo_,
-                                      const Array<double> &Gct_,
-                                      const Array<double> &Bot_,
-                                      const Vector &op_,
-                                      const Vector &D_,
-                                      Vector &diag_)
+static void PAHdivL2AssembleDiagonal_ADAt_2D(const int D1D,
+                                             const int Q1D,
+                                             const int L2D1D,
+                                             const int NE,
+                                             const Array<double> &L2Bo_,
+                                             const Array<double> &Gct_,
+                                             const Array<double> &Bot_,
+                                             const Vector &op_,
+                                             const Vector &D_,
+                                             Vector &diag_)
 {
    constexpr static int VDIM = 2;
 

fem/bilininteg_mass_ea.cpp

@@ -17,13 +17,13 @@ namespace mfem
 {
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EAMassAssemble1D(const int NE,
-                      const Array<double> &basis,
-                      const Vector &padata,
-                      Vector &eadata,
-                      const bool add,
-                      const int d1d = 0,
-                      const int q1d = 0)
+static void EAMassAssemble1D(const int NE,
+                             const Array<double> &basis,
+                             const Vector &padata,
+                             Vector &eadata,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -67,13 +67,13 @@ void EAMassAssemble1D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EAMassAssemble2D(const int NE,
-                      const Array<double> &basis,
-                      const Vector &padata,
-                      Vector &eadata,
-                      const bool add,
-                      const int d1d = 0,
-                      const int q1d = 0)
+static void EAMassAssemble2D(const int NE,
+                             const Array<double> &basis,
+                             const Vector &padata,
+                             Vector &eadata,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -139,13 +139,13 @@ void EAMassAssemble2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void EAMassAssemble3D(const int NE,
-                      const Array<double> &basis,
-                      const Vector &padata,
-                      Vector &eadata,
-                      const bool add,
-                      const int d1d = 0,
-                      const int q1d = 0)
+static void EAMassAssemble3D(const int NE,
+                             const Array<double> &basis,
+                             const Vector &padata,
+                             Vector &eadata,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_mass_pa.cpp

@@ -155,12 +155,12 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAMassAssembleDiagonal2D(const int NE,
-                              const Array<double> &b,
-                              const Vector &d,
-                              Vector &y,
-                              const int d1d = 0,
-                              const int q1d = 0)
+static void PAMassAssembleDiagonal2D(const int NE,
+                                     const Array<double> &b,
+                                     const Vector &d,
+                                     Vector &y,
+                                     const int d1d = 0,
+                                     const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -201,12 +201,12 @@ void PAMassAssembleDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-void SmemPAMassAssembleDiagonal2D(const int NE,
-                                  const Array<double> &b_,
-                                  const Vector &d_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+static void SmemPAMassAssembleDiagonal2D(const int NE,
+                                         const Array<double> &b_,
+                                         const Vector &d_,
+                                         Vector &y_,
+                                         const int d1d = 0,
+                                         const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -267,12 +267,12 @@ void SmemPAMassAssembleDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAMassAssembleDiagonal3D(const int NE,
-                              const Array<double> &b,
-                              const Vector &d,
-                              Vector &y,
-                              const int d1d = 0,
-                              const int q1d = 0)
+static void PAMassAssembleDiagonal3D(const int NE,
+                                     const Array<double> &b,
+                                     const Vector &d,
+                                     Vector &y,
+                                     const int d1d = 0,
+                                     const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -336,12 +336,12 @@ void PAMassAssembleDiagonal3D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void SmemPAMassAssembleDiagonal3D(const int NE,
-                                  const Array<double> &b_,
-                                  const Vector &d_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+static void SmemPAMassAssembleDiagonal3D(const int NE,
+                                         const Array<double> &b_,
+                                         const Vector &d_,
+                                         Vector &y_,
+                                         const int d1d = 0,
+                                         const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -569,14 +569,14 @@ static void OccaPAMassApply3D(const int D1D,
 #endif // MFEM_USE_OCCA
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAMassApply2D(const int NE,
-                   const Array<double> &b_,
-                   const Array<double> &bt_,
-                   const Vector &d_,
-                   const Vector &x_,
-                   Vector &y_,
-                   const int d1d = 0,
-                   const int q1d = 0)
+static void PAMassApply2D(const int NE,
+                          const Array<double> &b_,
+                          const Array<double> &bt_,
+                          const Vector &d_,
+                          const Vector &x_,
+                          Vector &y_,
+                          const int d1d = 0,
+                          const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -661,14 +661,14 @@ void PAMassApply2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-void SmemPAMassApply2D(const int NE,
-                       const Array<double> &b_,
-                       const Array<double> &bt_,
-                       const Vector &d_,
-                       const Vector &x_,
-                       Vector &y_,
-                       const int d1d = 0,
-                       const int q1d = 0)
+static void SmemPAMassApply2D(const int NE,
+                              const Array<double> &b_,
+                              const Array<double> &bt_,
+                              const Vector &d_,
+                              const Vector &x_,
+                              Vector &y_,
+                              const int d1d = 0,
+                              const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -784,14 +784,14 @@ void SmemPAMassApply2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAMassApply3D(const int NE,
-                   const Array<double> &b_,
-                   const Array<double> &bt_,
-                   const Vector &d_,
-                   const Vector &x_,
-                   Vector &y_,
-                   const int d1d = 0,
-                   const int q1d = 0)
+static void PAMassApply3D(const int NE,
+                          const Array<double> &b_,
+                          const Array<double> &bt_,
+                          const Vector &d_,
+                          const Vector &x_,
+                          Vector &y_,
+                          const int d1d = 0,
+                          const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -925,14 +925,14 @@ void PAMassApply3D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void SmemPAMassApply3D(const int NE,
-                       const Array<double> &b_,
-                       const Array<double> &bt_,
-                       const Vector &d_,
-                       const Vector &x_,
-                       Vector &y_,
-                       const int d1d = 0,
-                       const int q1d = 0)
+static void SmemPAMassApply3D(const int NE,
+                              const Array<double> &b_,
+                              const Array<double> &bt_,
+                              const Vector &d_,
+                              const Vector &x_,
+                              Vector &y_,
+                              const int d1d = 0,
+                              const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;

fem/bilininteg_vecdiffusion.cpp

@@ -22,12 +22,12 @@ namespace mfem
 // PA Vector Diffusion Integrator
 
 // PA Diffusion Assemble 2D kernel
-void PAVectorDiffusionSetup2D(const int Q1D,
-                              const int NE,
-                              const Array<double> &w,
-                              const Vector &j,
-                              const Vector &c,
-                              Vector &op)
+static void PAVectorDiffusionSetup2D(const int Q1D,
+                                     const int NE,
+                                     const Array<double> &w,
+                                     const Vector &j,
+                                     const Vector &c,
+                                     Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -59,12 +59,12 @@ void PAVectorDiffusionSetup2D(const int Q1D,
 }
 
 // PA Diffusion Assemble 3D kernel
-void PAVectorDiffusionSetup3D(const int Q1D,
-                              const int NE,
-                              const Array<double> &w,
-                              const Vector &j,
-                              const Vector &c,
-                              Vector &op)
+static void PAVectorDiffusionSetup3D(const int Q1D,
+                                     const int NE,
+                                     const Array<double> &w,
+                                     const Vector &j,
+                                     const Vector &c,
+                                     Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -251,7 +251,7 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 }
 
 // PA Diffusion Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_VDIM = 0>
+template<int T_D1D = 0, int T_Q1D = 0, int T_VDIM = 0> static
 void PAVectorDiffusionApply2D(const int NE,
                               const Array<double> &b,
                               const Array<double> &g,
@@ -374,7 +374,7 @@ void PAVectorDiffusionApply2D(const int NE,
 
 // PA Diffusion Apply 3D kernel
 template<const int T_D1D = 0,
-         const int T_Q1D = 0>
+         const int T_Q1D = 0> static
 void PAVectorDiffusionApply3D(const int NE,
                               const Array<double> &b,
                               const Array<double> &g,
@@ -606,13 +606,13 @@ void VectorDiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAVectorDiffusionDiagonal2D(const int NE,
-                                 const Array<double> &b,
-                                 const Array<double> &g,
-                                 const Vector &d,
-                                 Vector &y,
-                                 const int d1d = 0,
-                                 const int q1d = 0)
+static void PAVectorDiffusionDiagonal2D(const int NE,
+                                        const Array<double> &b,
+                                        const Array<double> &g,
+                                        const Vector &d,
+                                        Vector &y,
+                                        const int d1d = 0,
+                                        const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -673,13 +673,13 @@ void PAVectorDiffusionDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAVectorDiffusionDiagonal3D(const int NE,
-                                 const Array<double> &b,
-                                 const Array<double> &g,
-                                 const Vector &d,
-                                 Vector &y,
-                                 const int d1d = 0,
-                                 const int q1d = 0)
+static void PAVectorDiffusionDiagonal3D(const int NE,
+                                        const Array<double> &b,
+                                        const Array<double> &g,
+                                        const Vector &d,
+                                        Vector &y,
+                                        const int d1d = 0,
+                                        const int q1d = 0)
 {
    constexpr int DIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;

fem/bilininteg_vecmass.cpp

@@ -104,14 +104,14 @@ void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
-void PAVectorMassApply2D(const int NE,
-                         const Array<double> &B_,
-                         const Array<double> &Bt_,
-                         const Vector &op_,
-                         const Vector &x_,
-                         Vector &y_,
-                         const int d1d = 0,
-                         const int q1d = 0)
+static void PAVectorMassApply2D(const int NE,
+                                const Array<double> &B_,
+                                const Array<double> &Bt_,
+                                const Vector &op_,
+                                const Vector &x_,
+                                Vector &y_,
+                                const int d1d = 0,
+                                const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -201,14 +201,14 @@ void PAVectorMassApply2D(const int NE,
 
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
-void PAVectorMassApply3D(const int NE,
-                         const Array<double> &B_,
-                         const Array<double> &Bt_,
-                         const Vector &op_,
-                         const Vector &x_,
-                         Vector &y_,
-                         const int d1d = 0,
-                         const int q1d = 0)
+static void PAVectorMassApply3D(const int NE,
+                                const Array<double> &B_,
+                                const Array<double> &Bt_,
+                                const Vector &op_,
+                                const Vector &x_,
+                                Vector &y_,
+                                const int d1d = 0,
+                                const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -379,13 +379,13 @@ void VectorMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 }
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
-void PAVectorMassAssembleDiagonal2D(const int NE,
-                                    const Array<double> &B_,
-                                    const Array<double> &Bt_,
-                                    const Vector &op_,
-                                    Vector &diag_,
-                                    const int d1d = 0,
-                                    const int q1d = 0)
+static void PAVectorMassAssembleDiagonal2D(const int NE,
+                                           const Array<double> &B_,
+                                           const Array<double> &Bt_,
+                                           const Vector &op_,
+                                           Vector &diag_,
+                                           const int d1d = 0,
+                                           const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -431,13 +431,13 @@ void PAVectorMassAssembleDiagonal2D(const int NE,
 }
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
-void PAVectorMassAssembleDiagonal3D(const int NE,
-                                    const Array<double> &B_,
-                                    const Array<double> &Bt_,
-                                    const Vector &op_,
-                                    Vector &diag_,
-                                    const int d1d = 0,
-                                    const int q1d = 0)
+static void PAVectorMassAssembleDiagonal3D(const int NE,
+                                           const Array<double> &B_,
+                                           const Array<double> &Bt_,
+                                           const Vector &op_,
+                                           Vector &diag_,
+                                           const int d1d = 0,
+                                           const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/ceed/solvers-atpmg.cpp

@@ -178,10 +178,10 @@ int CeedATPMGElemRestriction(int order,
                   {
                      left_in_edof = i*P1d + 0;
                      right_in_edof = i*P1d + (P1d - 1);
-                     left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]
-                                                + e*in_layout[2]] + rounding_guard;
-                     right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]
-                                                 + e*in_layout[2]] + rounding_guard;
+                     left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]+e*in_layout[2]]
+                                    + rounding_guard;
+                     right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]+e*in_layout[2]]
+                                     + rounding_guard;
                      coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
                      coarse_j = (left_in_ldof < right_in_ldof) ?
                                 coarse_1d_edof(j, P1d, coarse_P1d) : reverse_coarse_1d_edof(j, P1d, coarse_P1d);
@@ -190,10 +190,10 @@ int CeedATPMGElemRestriction(int order,
                   {
                      left_in_edof = 0*P1d + j;
                      right_in_edof = (P1d - 1)*P1d + j;
-                     left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]
-                                                + e*in_layout[2]] + rounding_guard;
-                     right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]
-                                                 + e*in_layout[2]] + rounding_guard;
+                     left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]+e*in_layout[2]]
+                                    + rounding_guard;
+                     right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]+e*in_layout[2]]
+                                     + rounding_guard;
                      coarse_i = (left_in_ldof < right_in_ldof) ?
                                 coarse_1d_edof(i, P1d, coarse_P1d) : reverse_coarse_1d_edof(i, P1d, coarse_P1d);
                      coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
@@ -234,8 +234,8 @@ int CeedATPMGElemRestriction(int order,
                   // Determine topology; is this edof on the outside of the element
                   // in the i, j, or k direction?
                   int in_edof = i*P1d*P1d + j*P1d + k;
-                  int in_ldof = in_elem_dof[in_edof*in_layout[0] +
-                                            e*in_layout[2]] + rounding_guard;
+                  int in_ldof = in_elem_dof[in_edof*in_layout[0]+e*in_layout[2]]
+                                + rounding_guard;
                   bool i_edge = (i == 0 || i == P1d - 1);
                   bool j_edge = (j == 0 || j == P1d - 1);
                   bool k_edge = (k == 0 || k == P1d - 1);
@@ -265,10 +265,10 @@ int CeedATPMGElemRestriction(int order,
                      {
                         left_in_edof = 0*P1d*P1d + j*P1d + k;
                         right_in_edof = (P1d - 1)*P1d*P1d + j*P1d + k;
-                        left_in_ldof = in_elem_dof[left_in_edof*in_layout[0] +
-                                                   e*in_layout[2]] + rounding_guard;
-                        right_in_ldof = in_elem_dof[right_in_edof*in_layout[0] +
-                                                    e*in_layout[2]] + rounding_guard;
+                        left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]+e*in_layout[2]]
+                                       + rounding_guard;
+                        right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]+e*in_layout[2]]
+                                        + rounding_guard;
                         coarse_i = (left_in_ldof < right_in_ldof) ?
                                    coarse_1d_edof(i, P1d, coarse_P1d) : reverse_coarse_1d_edof(i, P1d, coarse_P1d);
                         coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
@@ -278,10 +278,10 @@ int CeedATPMGElemRestriction(int order,
                      {
                         left_in_edof = i*P1d*P1d + 0*P1d + k;
                         right_in_edof = i*P1d*P1d + (P1d - 1)*P1d + k;
-                        left_in_ldof = in_elem_dof[left_in_edof*in_layout[0] +
-                                                   e*in_layout[2]] + rounding_guard;
-                        right_in_ldof = in_elem_dof[right_in_edof*in_layout[0] +
-                                                    e*in_layout[2]] + rounding_guard;
+                        left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]+e*in_layout[2]]
+                                       + rounding_guard;
+                        right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]+e*in_layout[2]]
+                                        + rounding_guard;
                         coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
                         coarse_j = (left_in_ldof < right_in_ldof) ?
                                    coarse_1d_edof(j, P1d, coarse_P1d) : reverse_coarse_1d_edof(j, P1d, coarse_P1d);
@@ -296,10 +296,10 @@ int CeedATPMGElemRestriction(int order,
                         }
                         left_in_edof = i*P1d*P1d + j*P1d + 0;
                         right_in_edof = i*P1d*P1d + j*P1d + (P1d - 1);
-                        left_in_ldof = in_elem_dof[left_in_edof*in_layout[0] +
-                                                   e*in_layout[2]] + rounding_guard;
-                        right_in_ldof = in_elem_dof[right_in_edof*in_layout[0] +
-                                                    e*in_layout[2]] + rounding_guard;
+                        left_in_ldof = in_elem_dof[left_in_edof*in_layout[0]+e*in_layout[2]]
+                                       + rounding_guard;
+                        right_in_ldof = in_elem_dof[right_in_edof*in_layout[0]+e*in_layout[2]]
+                                        + rounding_guard;
                         coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
                         coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
                         coarse_k = (left_in_ldof < right_in_ldof) ?
@@ -323,14 +323,14 @@ int CeedATPMGElemRestriction(int order,
                         bottom_right_edof = i*P1d*P1d + 0*P1d + (P1d - 1);
                         top_right_edof = i*P1d*P1d + (P1d - 1)*P1d + (P1d - 1);
                         top_left_edof = i*P1d*P1d + (P1d - 1)*P1d + 0;
-                        bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0] +
-                                                       e*in_layout[2]] + rounding_guard;
-                        bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0] +
-                                                        e*in_layout[2]] + rounding_guard;
-                        top_right_ldof = in_elem_dof[top_right_edof*in_layout[0] +
-                                                     e*in_layout[2]] + rounding_guard;
-                        top_left_ldof = in_elem_dof[top_left_edof*in_layout[0] +
-                                                    e*in_layout[2]] + rounding_guard;
+                        bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0]+e*in_layout[2]]
+                                           + rounding_guard;
+                        bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0]+e*in_layout[2]]
+                                            + rounding_guard;
+                        top_right_ldof = in_elem_dof[top_right_edof*in_layout[0]+e*in_layout[2]]
+                                         + rounding_guard;
+                        top_left_ldof = in_elem_dof[top_left_edof*in_layout[0]+e*in_layout[2]]
+                                        + rounding_guard;
                         int m = min4(bottom_left_ldof, bottom_right_ldof, top_right_ldof,
                                      top_left_ldof);
                         coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
@@ -361,14 +361,14 @@ int CeedATPMGElemRestriction(int order,
                         bottom_right_edof = 0*P1d*P1d + j*P1d + (P1d - 1);
                         top_right_edof = (P1d - 1)*P1d*P1d + j*P1d + (P1d - 1);
                         top_left_edof = (P1d - 1)*P1d*P1d + j*P1d + 0;
-                        bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0] +
-                                                       e*in_layout[2]] + rounding_guard;
-                        bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0] +
-                                                        e*in_layout[2]] + rounding_guard;
-                        top_right_ldof = in_elem_dof[top_right_edof*in_layout[0] +
-                                                     e*in_layout[2]] + rounding_guard;
-                        top_left_ldof = in_elem_dof[top_left_edof*in_layout[0] +
-                                                    e*in_layout[2]] + rounding_guard;
+                        bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0]+e*in_layout[2]]
+                                           + rounding_guard;
+                        bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0]+e*in_layout[2]]
+                                            + rounding_guard;
+                        top_right_ldof = in_elem_dof[top_right_edof*in_layout[0]+e*in_layout[2]]
+                                         + rounding_guard;
+                        top_left_ldof = in_elem_dof[top_left_edof*in_layout[0]+e*in_layout[2]]
+                                        + rounding_guard;
                         int m = min4(bottom_left_ldof, bottom_right_ldof, top_right_ldof,
                                      top_left_ldof);
                         coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
@@ -404,14 +404,14 @@ int CeedATPMGElemRestriction(int order,
                         bottom_right_edof = 0*P1d*P1d + (P1d - 1)*P1d + k;
                         top_right_edof = (P1d - 1)*P1d*P1d + (P1d - 1)*P1d + k;
                         top_left_edof = (P1d - 1)*P1d*P1d + 0*P1d + k;
-                        bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0] +
-                                                       e*in_layout[2]] + rounding_guard;
-                        bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0] +
-                                                        e*in_layout[2]] + rounding_guard;
-                        top_right_ldof = in_elem_dof[top_right_edof*in_layout[0] +
-                                                     e*in_layout[2]] + rounding_guard;
-                        top_left_ldof = in_elem_dof[top_left_edof*in_layout[0] +
-                                                    e*in_layout[2]] + rounding_guard;
+                        bottom_left_ldof = in_elem_dof[bottom_left_edof*in_layout[0]+e*in_layout[2]]
+                                           + rounding_guard;
+                        bottom_right_ldof = in_elem_dof[bottom_right_edof*in_layout[0]+e*in_layout[2]]
+                                            + rounding_guard;
+                        top_right_ldof = in_elem_dof[top_right_edof*in_layout[0]+e*in_layout[2]]
+                                         + rounding_guard;
+                        top_left_ldof = in_elem_dof[top_left_edof*in_layout[0]+e*in_layout[2]]
+                                        + rounding_guard;
                         int m = min4(bottom_left_ldof, bottom_right_ldof,
                                      top_right_ldof, top_left_ldof);
                         coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);

fem/coefficient.cpp

@@ -431,6 +431,7 @@ MatrixArrayCoefficient::~MatrixArrayCoefficient ()
 void MatrixArrayCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                   const IntegrationPoint &ip)
 {
+   K.SetSize(height, width);
    for (int i = 0; i < height; i++)
    {
       for (int j = 0; j < width; j++)

fem/complex_fem.cpp

@@ -635,6 +635,320 @@ SesquilinearForm::Update(FiniteElementSpace *nfes)
 }
 
 
+bool MixedSesquilinearForm::RealInteg()
+{
+   int nint = blfr->GetTFBFI()->Size() + blfr->GetDBFI()->Size() +
+              blfr->GetBBFI()->Size() + blfr->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+bool MixedSesquilinearForm::ImagInteg()
+{
+   int nint = blfi->GetTFBFI()->Size() + blfi->GetDBFI()->Size() +
+              blfi->GetBBFI()->Size() + blfi->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+MixedSesquilinearForm::MixedSesquilinearForm(FiniteElementSpace *tr_f,
+                                             FiniteElementSpace *te_f,
+                                             ComplexOperator::Convention
+                                             convention)
+   : conv(convention),
+     blfr(new MixedBilinearForm(tr_f, te_f)),
+     blfi(new MixedBilinearForm(tr_f, te_f))
+{}
+
+MixedSesquilinearForm::MixedSesquilinearForm(FiniteElementSpace *tr_f,
+                                             FiniteElementSpace *te_f,
+                                             MixedBilinearForm *bfr,
+                                             MixedBilinearForm *bfi,
+                                             ComplexOperator::Convention
+                                             convention)
+   : conv(convention),
+     blfr(new MixedBilinearForm(tr_f, te_f, bfr)),
+     blfi(new MixedBilinearForm(tr_f, te_f, bfi))
+{}
+
+MixedSesquilinearForm::~MixedSesquilinearForm()
+{
+   delete blfr;
+   delete blfi;
+}
+
+void
+MixedSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                                           BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { blfr->AddDomainIntegrator(bfi_real); }
+   if (bfi_imag) { blfi->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+MixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                             BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { blfr->AddBoundaryIntegrator(bfi_real); }
+   if (bfi_imag) { blfi->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+MixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                             BilinearFormIntegrator *bfi_imag,
+                                             Array<int> & bdr_marker)
+{
+   if (bfi_real) { blfr->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { blfi->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+MixedSesquilinearForm::AddTraceFaceIntegrator(BilinearFormIntegrator
+                                              *bfi_real,
+                                              BilinearFormIntegrator
+                                              *bfi_imag)
+{
+   if (bfi_real) { blfr->AddTraceFaceIntegrator(bfi_real); }
+   if (bfi_imag) { blfi->AddTraceFaceIntegrator(bfi_imag); }
+}
+
+void
+MixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                 *bfi_real,
+                                                 BilinearFormIntegrator
+                                                 *bfi_imag)
+{
+   if (bfi_real) { blfr->AddBdrTraceFaceIntegrator(bfi_real); }
+   if (bfi_imag) { blfi->AddBdrTraceFaceIntegrator(bfi_imag); }
+}
+
+void
+MixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                 *bfi_real,
+                                                 BilinearFormIntegrator
+                                                 *bfi_imag,
+                                                 Array<int> &bdr_marker)
+{
+   if (bfi_real) { blfr->AddBdrTraceFaceIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { blfi->AddBdrTraceFaceIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+MixedSesquilinearForm::Assemble(int skip_zeros)
+{
+   blfr->Assemble(skip_zeros);
+   blfi->Assemble(skip_zeros);
+}
+
+void
+MixedSesquilinearForm::Finalize(int skip_zeros)
+{
+   blfr->Finalize(skip_zeros);
+   blfi->Finalize(skip_zeros);
+}
+
+ComplexSparseMatrix *
+MixedSesquilinearForm::AssembleComplexSparseMatrix()
+{
+   return new ComplexSparseMatrix(&blfr->SpMat(),
+                                  &blfi->SpMat(),
+                                  false, false, conv);
+}
+
+void
+MixedSesquilinearForm::FormRectangularLinearSystem(const Array<int>
+                                                   &trial_tdof_list,
+                                                   const Array<int>
+                                                   &test_tdof_list,
+                                                   Vector &x, Vector &b,
+                                                   OperatorHandle &A,
+                                                   Vector &X, Vector &B)
+{
+   FiniteElementSpace *tr_fes = blfr->TrialFESpace();
+   FiniteElementSpace *te_fes = blfr->TestFESpace();
+   const int tr_vsize = tr_fes->GetVSize();
+   const int te_vsize = te_fes->GetVSize();
+
+   // Allocate temporary vector
+   Vector b_0;
+   b_0.UseDevice(true);
+   b_0.SetSize(tr_vsize);
+   b_0 = 0.0;
+
+   // Extract the real and imaginary parts of the input vectors
+   MFEM_ASSERT(x.Size() == 2 * te_vsize,
+               "Input GridFunction of incorrect size!");
+   x.Read();
+   Vector x_r; x_r.MakeRef(x, 0, te_vsize);
+   Vector x_i; x_i.MakeRef(x, te_vsize, te_vsize);
+
+   MFEM_ASSERT(b.Size() == 2 * tr_vsize, "Input LinearForm of incorrect size!");
+   b.Read();
+   Vector b_r; b_r.MakeRef(b, 0, tr_vsize);
+   Vector b_i; b_i.MakeRef(b, tr_vsize, tr_vsize);
+
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC) { b_i *= -1.0; }
+
+   const int tr_tvsize = tr_fes->GetTrueVSize();
+   const int te_tvsize = te_fes->GetTrueVSize();
+   OperatorHandle A_r, A_i;
+
+   X.UseDevice(true);
+   X.SetSize(2 * te_tvsize);
+   X = 0.0;
+
+   B.UseDevice(true);
+   B.SetSize(2 * tr_tvsize);
+   B = 0.0;
+
+   Vector X_r; X_r.MakeRef(X, 0, te_tvsize);
+   Vector X_i; X_i.MakeRef(X, te_tvsize, te_tvsize);
+   Vector B_r; B_r.MakeRef(B, 0, tr_tvsize);
+   Vector B_i; B_i.MakeRef(B, tr_tvsize, tr_tvsize);
+
+   Vector X_0, B_0;
+
+   if (RealInteg())
+   {
+      b_0 = b_r;
+      blfr->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
+                                        x_r, b_0, A_r, X_0, B_0);
+      X_r = X_0; B_r = B_0;
+
+      b_0 = b_i;
+      blfr->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
+                                        x_i, b_0, A_r, X_0, B_0);
+      X_i = X_0; B_i = B_0;
+
+      if (ImagInteg())
+      {
+         b_0 = 0.0;
+         blfi->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
+                                           x_i, b_0, A_i, X_0, B_0);
+         B_r -= B_0;
+
+         b_0 = 0.0;
+         blfi->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
+                                           x_r, b_0, A_i, X_0, B_0);
+         B_i += B_0;
+      }
+   }
+   else if (ImagInteg())
+   {
+      b_0 = b_i;
+      blfi->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
+                                        x_r, b_0, A_i, X_0, B_0);
+      X_r = X_0; B_i = B_0;
+
+      b_0 = b_r; b_0 *= -1.0;
+      blfi->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
+                                        x_i, b_0, A_i, X_0, B_0);
+      X_i = X_0; B_r = B_0; B_r *= -1.0;
+   }
+   else
+   {
+      MFEM_ABORT("Real and Imaginary part of the Sesquilinear form are empty");
+   }
+
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      B_i *= -1.0;
+      b_i *= -1.0;
+   }
+
+   x_r.SyncAliasMemory(x);
+   x_i.SyncAliasMemory(x);
+   b_r.SyncAliasMemory(b);
+   b_i.SyncAliasMemory(b);
+
+   X_r.SyncAliasMemory(X);
+   X_i.SyncAliasMemory(X);
+   B_r.SyncAliasMemory(B);
+   B_i.SyncAliasMemory(B);
+
+   // A = A_r + i A_i
+   A.Clear();
+   if ( A_r.Type() == Operator::MFEM_SPARSEMAT ||
+        A_i.Type() == Operator::MFEM_SPARSEMAT )
+   {
+      ComplexSparseMatrix * A_sp =
+         new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
+                                 A_i.As<SparseMatrix>(),
+                                 A_r.OwnsOperator(),
+                                 A_i.OwnsOperator(),
+                                 conv);
+      A.Reset<ComplexSparseMatrix>(A_sp, true);
+   }
+   else
+   {
+      ComplexOperator * A_op =
+         new ComplexOperator(A_r.Ptr(),
+                             A_i.Ptr(),
+                             A_r.OwnsOperator(),
+                             A_i.OwnsOperator(),
+                             conv);
+      A.Reset<ComplexOperator>(A_op, true);
+   }
+   A_r.SetOperatorOwner(false);
+   A_i.SetOperatorOwner(false);
+}
+
+void
+MixedSesquilinearForm::FormRectangularSystemMatrix(const Array<int>
+                                                   &trial_tdof_list,
+                                                   const Array<int>
+                                                   &test_tdof_list,
+                                                   OperatorHandle &A)
+
+{
+   OperatorHandle A_r, A_i;
+   if (RealInteg())
+   {
+      blfr->FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A_r);
+   }
+   if (ImagInteg())
+   {
+      blfi->FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A_i);
+   }
+   if (!RealInteg() && !ImagInteg())
+   {
+      MFEM_ABORT("Both Real and Imaginary part of the mixed Sesquilinear form "
+                 "are empty");
+   }
+
+   // A = A_r + i A_i
+   A.Clear();
+   if ( A_r.Type() == Operator::MFEM_SPARSEMAT ||
+        A_i.Type() == Operator::MFEM_SPARSEMAT )
+   {
+      ComplexSparseMatrix * A_sp =
+         new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
+                                 A_i.As<SparseMatrix>(),
+                                 A_r.OwnsOperator(),
+                                 A_i.OwnsOperator(),
+                                 conv);
+      A.Reset<ComplexSparseMatrix>(A_sp, true);
+   }
+   else
+   {
+      ComplexOperator * A_op =
+         new ComplexOperator(A_r.Ptr(),
+                             A_i.Ptr(),
+                             A_r.OwnsOperator(),
+                             A_i.OwnsOperator(),
+                             conv);
+      A.Reset<ComplexOperator>(A_op, true);
+   }
+   A_r.SetOperatorOwner(false);
+   A_i.SetOperatorOwner(false);
+}
+
+void
+MixedSesquilinearForm::Update()
+{
+   if ( blfr ) { blfr->Update(); }
+   if ( blfi ) { blfi->Update(); }
+}
+
+
 #ifdef MFEM_USE_MPI
 
 ParComplexGridFunction::ParComplexGridFunction(ParFiniteElementSpace *pfes)

fem/complex_fem.hpp

@@ -324,6 +324,150 @@ public:
    virtual ~SesquilinearForm();
 };
 
+/** Class for mixed sesquilinear form
+
+    A mixed sesquilinear form is a generalization of a mixed bilinear form to
+    complex-valued fields. Sesquilinear forms are linear in the second argument
+    but the first argument involves a complex conjugate in the sense that:
+
+                a(alpha u, beta v) = conj(alpha) beta a(u, v)
+
+    The 'Mixed' keyword indicates that the domain and range spaces may differ
+    leading to rectangular operators.
+
+    The @a convention argument in the class's constructor is documented in the
+    mfem::ComplexOperator class found in linalg/complex_operator.hpp.
+
+    When supplying integrators to the MixedSesquilinearForm either the real or
+    imaginary integrator can be NULL. This indicates that the corresponding
+    portion of the complex-valued material coefficient is equal to zero.
+*/
+class MixedSesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv;
+
+   MixedBilinearForm *blfr;
+   MixedBilinearForm *blfi;
+
+   /* These methods check if the real/imag parts of the sesquilinear form are
+      not empty */
+   bool RealInteg();
+   bool ImagInteg();
+
+public:
+   MixedSesquilinearForm(FiniteElementSpace *tr_fes,
+                         FiniteElementSpace *te_fes,
+                         ComplexOperator::Convention
+                         convention = ComplexOperator::HERMITIAN);
+   /** @brief Create a MixedSesquilinearForm given a trial space,
+       FiniteElementSpace @a tr_fes, and a test space,
+       FiniteElementSpace @a te_fes, using the same integrators as the
+       MixedBilinearForms @a bfr and @a bfi .
+
+       The pointers @a tr_fes and @a te_fes are not owned by the newly
+       constructed object.
+
+       The integrators are copied as pointers and they are not owned by the
+       newly constructed MixedSesquilinearForm. */
+   MixedSesquilinearForm(FiniteElementSpace *tr_fes,
+                         FiniteElementSpace *te_fes,
+                         MixedBilinearForm *bfr, MixedBilinearForm *bfi,
+                         ComplexOperator::Convention
+                         convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv; }
+   void SetConvention(const ComplexOperator::Convention &
+                      convention) { conv = convention; }
+
+   /// Set the desired assembly level.
+   /** Valid choices are:
+
+       - AssemblyLevel::LEGACY (default)
+       - AssemblyLevel::FULL
+       - AssemblyLevel::PARTIAL
+       - AssemblyLevel::ELEMENT
+       - AssemblyLevel::NONE
+
+       This method must be called before assembly. */
+   void SetAssemblyLevel(AssemblyLevel assembly_level)
+   {
+      blfr->SetAssemblyLevel(assembly_level);
+      blfi->SetAssemblyLevel(assembly_level);
+   }
+
+   MixedBilinearForm & real() { return *blfr; }
+   MixedBilinearForm & imag() { return *blfi; }
+   const MixedBilinearForm & real() const { return *blfr; }
+   const MixedBilinearForm & imag() const { return *blfi; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                            BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag,
+                              Array<int> &bdr_marker);
+
+   /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
+   void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                               BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                                  BilinearFormIntegrator *bfi_imag);
+
+   /** @brief Adds new boundary Face Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The array @a bdr_marker is stored internally as a pointer to the given
+       Array<int> object. */
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                                  BilinearFormIntegrator *bfi_imag,
+                                  Array<int> &bdr_marker);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   /// Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexSparseMatrix *AssembleComplexSparseMatrix();
+
+   /// Return the trial FE space associated with the BilinearForm.
+   FiniteElementSpace *TrialFESpace() { return blfr->TrialFESpace(); }
+   /// Read-only access to the associated trial FiniteElementSpace.
+   const FiniteElementSpace *TrialFESpace() const
+   { return blfr->TrialFESpace(); }
+
+   /// Return the test FE space associated with the BilinearForm.
+   FiniteElementSpace *TestFESpace() { return blfr->TestFESpace(); }
+   /// Read-only access to the associated test FiniteElementSpace.
+   const FiniteElementSpace *TestFESpace() const { return blfr->TestFESpace(); }
+
+   void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
+                                    const Array<int> &test_tdof_list,
+                                    Vector &x, Vector &b,
+                                    OperatorHandle &A, Vector &X, Vector &B);
+
+   void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
+                                    const Array<int> &test_tdof_list,
+                                    OperatorHandle &A);
+
+   virtual void Update();
+
+   virtual ~MixedSesquilinearForm();
+};
+
 #ifdef MFEM_USE_MPI
 
 /// Class for parallel complex-valued grid function - real + imaginary part

fem/datacollection.cpp

@@ -482,7 +482,7 @@ void VisItDataCollection::SaveRootFile()
    std::string root_name = prefix_path + name + "_" +
                            to_padded_string(cycle, pad_digits_cycle) +
                            ".mfem_root";
-   std::ofstream root_file(root_name.c_str());
+   std::ofstream root_file(root_name);
    root_file << GetVisItRootString();
    if (!root_file)
    {
@@ -548,7 +548,7 @@ void VisItDataCollection::Load(int cycle_)
 
 void VisItDataCollection::LoadVisItRootFile(const std::string& root_name)
 {
-   std::ifstream root_file(root_name.c_str());
+   std::ifstream root_file(root_name);
    std::stringstream buffer;
    buffer << root_file.rdbuf();
    if (!buffer)
@@ -853,6 +853,7 @@ void ParaViewDataCollection::Save()
       std::string dpath=GenerateCollectionPath();
       std::string pvdname=dpath+"/"+GeneratePVDFileName();
 
+      bool write_header = true;
       std::ifstream pvd_in;
       if (restart_mode && (pvd_in.open(pvdname,std::ios::binary),pvd_in.good()))
       {
@@ -879,20 +880,34 @@ void ParaViewDataCollection::Save()
                pos_end = pvd_in.tellg();
             }
          }
+         // Since pvd_in is opened in binary mode, count will store the number
+         // of bytes from the beginning of the file until the desired insertion
+         // point (in text mode on Windows this is not the case).
          size_t count = pos_end - pos_begin;
-         std::vector<char> buf(count);
-         pvd_in.clear();
-         pvd_in.seekg(pos_begin);
-         pvd_in.read(buf.data(), count);
-         pvd_in.close();
-         pvd_stream.open(pvdname.c_str(),std::ios::out);
-         pvd_stream.write(buf.data(), count);
+         if (count != 0)
+         {
+            write_header = false;
+            std::vector<char> buf(count);
+            // Read the contents of the PVD file, from the beginning to the
+            // insertion point.
+            pvd_in.clear();
+            pvd_in.seekg(pos_begin);
+            pvd_in.read(buf.data(), count);
+            pvd_in.close();
+            // Open the PVD file in truncate mode to delete the previous
+            // contents. Open in binary mode to write the data buffer without
+            // converting \r\n to \r\r\n on Windows.
+            pvd_stream.open(pvdname,std::ios::out|std::ios::trunc|std::ios::binary);
+            pvd_stream.write(buf.data(), count);
+            // Close and reopen the file in text mode, appending to the end.
+            pvd_stream.close();
+            pvd_stream.open(pvdname,std::ios::in|std::ios::out|std::ios::ate);
+         }
       }
-      else
+      if (write_header)
       {
-         // initialize new pvd file
-         pvd_stream.open(pvdname.c_str(),std::ios::out);
-         // initialize the file
+         // Initialize new pvd file.
+         pvd_stream.open(pvdname,std::ios::out|std::ios::trunc);
          pvd_stream << "<?xml version=\"1.0\"?>\n";
          pvd_stream << "<VTKFile type=\"Collection\" version=\"0.1\"";
          pvd_stream << " byte_order=\"" << VTKByteOrder() << "\">\n";
@@ -904,7 +919,7 @@ void ParaViewDataCollection::Save()
    {
       std::string fname = GenerateCollectionPath()+"/"+GenerateVTUPath()+"/"
                           +GenerateVTUFileName();
-      std::fstream out(fname.c_str(), std::ios::out);
+      std::fstream out(fname, std::ios::out);
       out.precision(precision);
       SaveDataVTU(out,levels_of_detail);
       out.close();
@@ -915,7 +930,7 @@ void ParaViewDataCollection::Save()
    {
       std::string fname = GenerateCollectionPath()+"/"+GeneratePVTUPath()+"/"
                           +GeneratePVTUFileName();
-      std::fstream out(fname.c_str(), std::ios::out);
+      std::fstream out(fname, std::ios::out);
 
       out << "<?xml version=\"1.0\"?>\n";
       out << "<VTKFile type=\"PUnstructuredGrid\"";
@@ -973,6 +988,7 @@ void ParaViewDataCollection::Save()
       pvd_stream << "<DataSet timestep=\"" << GetTime();  // GetCycle();
       pvd_stream << "\" group=\"\" part=\"" << 0 << "\" file=\"";
       pvd_stream << fname << "\"/>\n";
+      pvd_stream.flush();
       std::fstream::pos_type pos = pvd_stream.tellp();
       pvd_stream << "</Collection>\n";
       pvd_stream << "</VTKFile>" << std::endl;

fem/fe.cpp

@@ -4272,6 +4272,7 @@ void RT1QuadFiniteElement::Project (
 
 const double RT2TriangleFiniteElement::M[15][15] =
 {
+   // *INDENT-OFF*
    {
       0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
       0, 24.442740046346700787, -16.647580015448900262, -12.,
@@ -4325,6 +4326,7 @@ const double RT2TriangleFiniteElement::M[15][15] =
    { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
    { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
    { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
+   // *INDENT-ON*
 };
 
 RT2TriangleFiniteElement::RT2TriangleFiniteElement()
@@ -12134,6 +12136,7 @@ void RT_QuadrilateralElement::CalcVShape(const IntegrationPoint &ip,
 
 #ifdef MFEM_THREAD_SAFE
    Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
+   Vector dshape_cx, dshape_cy;
 #endif
 
    if (obasis1d.IsIntegratedType())
@@ -12459,6 +12462,7 @@ void RT_HexahedronElement::CalcVShape(const IntegrationPoint &ip,
 #ifdef MFEM_THREAD_SAFE
    Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
    Vector shape_cz(pp1 + 1), shape_oz(pp1);
+   Vector dshape_cx, dshape_cy, dshape_cz;
 #endif
 
    if (obasis1d.IsIntegratedType())
@@ -13300,6 +13304,7 @@ void ND_HexahedronElement::CalcVShape(const IntegrationPoint &ip,
 #ifdef MFEM_THREAD_SAFE
    Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
    Vector shape_cz(p + 1), shape_oz(p);
+   Vector dshape_cx, dshape_cy, dshape_cz;
 #endif
 
    if (obasis1d.IsIntegratedType())
@@ -13726,6 +13731,7 @@ void ND_QuadrilateralElement::CalcVShape(const IntegrationPoint &ip,
 
 #ifdef MFEM_THREAD_SAFE
    Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
+   Vector dshape_cx, dshape_cy;
 #endif
 
    if (obasis1d.IsIntegratedType())

fem/gridfunc.cpp

@@ -502,25 +502,24 @@ const
    {
       doftrans->InvTransformPrimal(loc_data);
    }
-   for (int k = 0; k < n; k++)
-      if (FElem->GetMapType() == FiniteElement::VALUE)
+   if (FElem->GetMapType() == FiniteElement::VALUE)
+   {
+      for (int k = 0; k < n; k++)
       {
-         for (int k = 0; k < n; k++)
-         {
-            FElem->CalcShape(ir.IntPoint(k), DofVal);
-            vals(k) = DofVal * loc_data;
-         }
+         FElem->CalcShape(ir.IntPoint(k), DofVal);
+         vals(k) = DofVal * loc_data;
       }
-      else
+   }
+   else
+   {
+      ElementTransformation *Tr = fes->GetElementTransformation(i);
+      for (int k = 0; k < n; k++)
       {
-         ElementTransformation *Tr = fes->GetElementTransformation(i);
-         for (int k = 0; k < n; k++)
-         {
-            Tr->SetIntPoint(&ir.IntPoint(k));
-            FElem->CalcPhysShape(*Tr, DofVal);
-            vals(k) = DofVal * loc_data;
-         }
+         Tr->SetIntPoint(&ir.IntPoint(k));
+         FElem->CalcPhysShape(*Tr, DofVal);
+         vals(k) = DofVal * loc_data;
       }
+   }
 }
 
 void GridFunction::GetValues(int i, const IntegrationRule &ir, Vector &vals,

fem/intrules.cpp

@@ -955,7 +955,7 @@ const IntegrationRule &IntegrationRules::Get(int GeomType, int Order)
             IntegrationRule *ir = GenerateIntegrationRule(GeomType, Order);
             int RealOrder = Order;
             while (RealOrder+1 < ir_array->Size() &&
-            /*  */ (*ir_array)[RealOrder+1] == ir)
+                   (*ir_array)[RealOrder+1] == ir)
             {
                RealOrder++;
             }

fem/nonlininteg_vectorconvection.cpp

@@ -116,15 +116,15 @@ void VectorConvectionNLFIntegrator::AssemblePA(const FiniteElementSpace &fes)
 
 // PA Convection NL 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAConvectionNLApply2D(const int NE,
-                           const Array<double> &b,
-                           const Array<double> &g,
-                           const Array<double> &bt,
-                           const Vector &q_,
-                           const Vector &x_,
-                           Vector &y_,
-                           const int d1d = 0,
-                           const int q1d = 0)
+static void PAConvectionNLApply2D(const int NE,
+                                  const Array<double> &b,
+                                  const Array<double> &g,
+                                  const Array<double> &bt,
+                                  const Vector &q_,
+                                  const Vector &x_,
+                                  Vector &y_,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -252,15 +252,15 @@ void PAConvectionNLApply2D(const int NE,
 
 // PA Convection NL 3D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-void PAConvectionNLApply3D(const int NE,
-                           const Array<double> &b,
-                           const Array<double> &g,
-                           const Array<double> &bt,
-                           const Vector &q_,
-                           const Vector &x_,
-                           Vector &y_,
-                           const int d1d = 0,
-                           const int q1d = 0)
+static void PAConvectionNLApply3D(const int NE,
+                                  const Array<double> &b,
+                                  const Array<double> &g,
+                                  const Array<double> &bt,
+                                  const Vector &q_,
+                                  const Vector &x_,
+                                  Vector &y_,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    constexpr int VDIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -558,14 +558,14 @@ void PAConvectionNLApply3D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0, int T_MAX_D1D =0, int T_MAX_Q1D =0>
-void SmemPAConvectionNLApply3D(const int NE,
-                               const Array<double> &b_,
-                               const Array<double> &g_,
-                               const Vector &d_,
-                               const Vector &x_,
-                               Vector &y_,
-                               const int d1d = 0,
-                               const int q1d = 0)
+static void SmemPAConvectionNLApply3D(const int NE,
+                                      const Array<double> &b_,
+                                      const Array<double> &g_,
+                                      const Vector &d_,
+                                      const Vector &x_,
+                                      Vector &y_,
+                                      const int d1d = 0,
+                                      const int q1d = 0)
 {
    constexpr int VDIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;

fem/qinterp/det.cpp

@@ -27,14 +27,14 @@ namespace quadrature_interpolator
 {
 
 template<int T_D1D = 0, int T_Q1D = 0, int MAX_D1D = 0, int MAX_Q1D = 0>
-void Det2D(const int NE,
-           const double *b,
-           const double *g,
-           const double *x,
-           double *y,
-           const int vdim = 1,
-           const int d1d = 0,
-           const int q1d = 0)
+static void Det2D(const int NE,
+                  const double *b,
+                  const double *g,
+                  const double *x,
+                  double *y,
+                  const int vdim = 1,
+                  const int d1d = 0,
+                  const int q1d = 0)
 {
    constexpr int DIM = 2;
    static constexpr int NBZ = 1;
@@ -79,15 +79,15 @@ void Det2D(const int NE,
 
 template<int T_D1D = 0, int T_Q1D = 0, int MAX_D1D = 0, int MAX_Q1D = 0,
          bool SMEM = true>
-void Det3D(const int NE,
-           const double *b,
-           const double *g,
-           const double *x,
-           double *y,
-           const int vdim = 1,
-           const int d1d = 0,
-           const int q1d = 0,
-           Vector *d_buff = nullptr) // used only with SMEM = false
+static void Det3D(const int NE,
+                  const double *b,
+                  const double *g,
+                  const double *x,
+                  double *y,
+                  const int vdim = 1,
+                  const int d1d = 0,
+                  const int q1d = 0,
+                  Vector *d_buff = nullptr) // used only with SMEM = false
 {
    constexpr int DIM = 3;
    static constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;

fem/qinterp/eval.hpp

@@ -31,13 +31,13 @@ namespace quadrature_interpolator
 template<QVectorLayout Q_LAYOUT,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int T_NBZ = 1, int MAX_D1D = 0, int MAX_Q1D = 0>
-void Values2D(const int NE,
-              const double *b_,
-              const double *x_,
-              double *y_,
-              const int vdim = 0,
-              const int d1d = 0,
-              const int q1d = 0)
+static void Values2D(const int NE,
+                     const double *b_,
+                     const double *x_,
+                     double *y_,
+                     const int vdim = 0,
+                     const int d1d = 0,
+                     const int q1d = 0)
 {
    static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
 
@@ -95,13 +95,13 @@ void Values2D(const int NE,
 template<QVectorLayout Q_LAYOUT,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int MAX_D1D = 0, int MAX_Q1D = 0>
-void Values3D(const int NE,
-              const double *b_,
-              const double *x_,
-              double *y_,
-              const int vdim = 0,
-              const int d1d = 0,
-              const int q1d = 0)
+static void Values3D(const int NE,
+                     const double *b_,
+                     const double *x_,
+                     double *y_,
+                     const int vdim = 0,
+                     const int d1d = 0,
+                     const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/qinterp/grad.hpp

@@ -31,15 +31,15 @@ namespace quadrature_interpolator
 template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int T_NBZ = 1, int MAX_D1D = 0, int MAX_Q1D = 0>
-void Derivatives2D(const int NE,
-                   const double *b_,
-                   const double *g_,
-                   const double *j_,
-                   const double *x_,
-                   double *y_,
-                   const int vdim = 0,
-                   const int d1d = 0,
-                   const int q1d = 0)
+static void Derivatives2D(const int NE,
+                          const double *b_,
+                          const double *g_,
+                          const double *j_,
+                          const double *x_,
+                          double *y_,
+                          const int vdim = 0,
+                          const int d1d = 0,
+                          const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -139,15 +139,15 @@ void Derivatives2D(const int NE,
 template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int MAX_D1D = 0, int MAX_Q1D = 0>
-void Derivatives3D(const int NE,
-                   const double *b_,
-                   const double *g_,
-                   const double *j_,
-                   const double *x_,
-                   double *y_,
-                   const int vdim = 0,
-                   const int d1d = 0,
-                   const int q1d = 0)
+static void Derivatives3D(const int NE,
+                          const double *b_,
+                          const double *g_,
+                          const double *j_,
+                          const double *x_,
+                          double *y_,
+                          const int vdim = 0,
+                          const int d1d = 0,
+                          const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/quadinterpolator.cpp

@@ -61,16 +61,16 @@ namespace quadrature_interpolator
 // * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
 // * assumes 'maps.mode == FULL'.
 template<const int T_VDIM, const int T_ND, const int T_NQ>
-void Eval2D(const int NE,
-            const int vdim,
-            const QVectorLayout q_layout,
-            const GeometricFactors *geom,
-            const DofToQuad &maps,
-            const Vector &e_vec,
-            Vector &q_val,
-            Vector &q_der,
-            Vector &q_det,
-            const int eval_flags)
+static void Eval2D(const int NE,
+                   const int vdim,
+                   const QVectorLayout q_layout,
+                   const GeometricFactors *geom,
+                   const DofToQuad &maps,
+                   const Vector &e_vec,
+                   Vector &q_val,
+                   Vector &q_der,
+                   Vector &q_det,
+                   const int eval_flags)
 {
    using QI = QuadratureInterpolator;
 
@@ -209,16 +209,16 @@ void Eval2D(const int NE,
 // * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
 // * assumes 'maps.mode == FULL'.
 template<const int T_VDIM, const int T_ND, const int T_NQ>
-void Eval3D(const int NE,
-            const int vdim,
-            const QVectorLayout q_layout,
-            const GeometricFactors *geom,
-            const DofToQuad &maps,
-            const Vector &e_vec,
-            Vector &q_val,
-            Vector &q_der,
-            Vector &q_det,
-            const int eval_flags)
+static void Eval3D(const int NE,
+                   const int vdim,
+                   const QVectorLayout q_layout,
+                   const GeometricFactors *geom,
+                   const DofToQuad &maps,
+                   const Vector &e_vec,
+                   Vector &q_val,
+                   Vector &q_der,
+                   Vector &q_det,
+                   const int eval_flags)
 {
    using QI = QuadratureInterpolator;
 

fem/restriction.cpp

@@ -167,8 +167,8 @@ void ElementRestriction::MultTranspose(const Vector& x, Vector& y) const
          for (int j = offset; j < nextOffset; ++j)
          {
             const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
-            dofValue += (d_indices[j] >= 0) ? d_x(idx_j % nd, c,
-            idx_j / nd) : -d_x(idx_j % nd, c, idx_j / nd);
+            dofValue += ((d_indices[j] >= 0) ? d_x(idx_j % nd, c, idx_j / nd) :
+                         -d_x(idx_j % nd, c, idx_j / nd));
          }
          d_y(t?c:i,t?i:c) = dofValue;
       }

fem/tmop.cpp

@@ -1306,7 +1306,7 @@ namespace internal
 // MFEM_FORALL-based copy kernel -- used by protected methods below.
 // Needed as a workaround for the nvcc restriction that methods with MFEM_FORALL
 // in them must to be public.
-inline void device_copy(double *d_dest, const double *d_src, int size)
+static inline void device_copy(double *d_dest, const double *d_src, int size)
 {
    MFEM_FORALL(i, size, d_dest[i] = d_src[i];);
 }

fem/tmop_amr.cpp

@@ -234,7 +234,7 @@ void TMOPRefinerEstimator::SetTriIntRules()
 
    // Reftype = 0 // original element
    const int Nvert = 3, NEsplit = 1;
-   Mesh meshsplit(2, Nvert, NEsplit, 0,2);
+   Mesh meshsplit(2, Nvert, NEsplit, 0, 2);
    const double tri_v[3][2] =
    {
       {0, 0}, {1, 0}, {0, 1}

general/backends.hpp

@@ -58,6 +58,23 @@
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=0; i<N; i++)
 #endif
 
+// 'double' atomicAdd implementation for previous versions of CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600
+MFEM_DEVICE double atomicAdd(double *add, double val)
+{
+   unsigned long long int *ptr = (unsigned long long int *) add;
+   unsigned long long int old = *ptr, reg;
+   do
+   {
+      reg = old;
+      old = atomicCAS(ptr, reg,
+                      __double_as_longlong(val + __longlong_as_double(reg)));
+   }
+   while (reg != old);
+   return __longlong_as_double(old);
+}
+#endif
+
 template <typename T>
 MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
 {
@@ -66,6 +83,9 @@ MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
    return atomicAdd(&add,val);
 #else
    T old = add;
+#ifdef MFEM_USE_OPENMP
+   #pragma omp atomic
+#endif
    add += val;
    return old;
 #endif

general/error.hpp

@@ -164,13 +164,7 @@ __device__ void abort_msg(T & msg)
 #endif
 
 // Abort inside a device kernel
-#if defined(__CUDA_ARCH__) && defined(_WIN32)
-#define MFEM_ABORT_KERNEL(msg) \
-   {                           \
-      printf(msg);             \
-      __debugbreak();          \
-   }
-#elif defined(__CUDA_ARCH__)
+#if defined(__CUDA_ARCH__)
 #define MFEM_ABORT_KERNEL(msg) \
    {                           \
       printf(msg);             \

general/isockstream.cpp

@@ -23,7 +23,9 @@
 #include <unistd.h>
 #else
 #include <winsock.h>
+#ifdef _MSC_VER
 typedef int ssize_t;
+#endif
 typedef int socklen_t;
 #define close closesocket
 // Link with ws2_32.lib
@@ -79,13 +81,7 @@ int isockstream::establish()
    int on=1;
    setsockopt(port, SOL_SOCKET, SO_REUSEADDR, (char *)(&on), sizeof(on));
 
-   if (bind(
-#ifdef _WIN32
-          (SOCKET)port
-#else
-          port
-#endif
-          ,(const sockaddr*)&sa,(socklen_t)sizeof(struct sockaddr_in)) < 0)
+   if (bind(port,(const sockaddr*)&sa,(socklen_t)sizeof(struct sockaddr_in)) < 0)
    {
       mfem::err << "isockstream::establish(): bind() failed!" << endl;
       close(port);

general/socketstream.cpp

@@ -28,7 +28,9 @@
 #define closesocket (::close)
 #else
 #include <winsock.h>
+#ifdef _MSC_VER
 typedef int ssize_t;
+#endif
 // Link with ws2_32.lib
 #pragma comment(lib, "ws2_32.lib")
 #endif

linalg/hiop.cpp

@@ -22,7 +22,7 @@ using namespace hiop;
 namespace mfem
 {
 
-bool HiopOptimizationProblem::get_prob_sizes(long long &n, long long &m)
+bool HiopOptimizationProblem::get_prob_sizes(size_type &n, size_type &m)
 {
    n = ntdofs_glob;
    m = problem.GetNumConstraints();
@@ -30,7 +30,7 @@ bool HiopOptimizationProblem::get_prob_sizes(long long &n, long long &m)
    return true;
 }
 
-bool HiopOptimizationProblem::get_starting_point(const long long &n, double *x0)
+bool HiopOptimizationProblem::get_starting_point(const size_type &n, double *x0)
 {
    MFEM_ASSERT(x_start != NULL && ntdofs_loc == x_start->Size(),
                "Starting point is not set properly.");
@@ -40,7 +40,7 @@ bool HiopOptimizationProblem::get_starting_point(const long long &n, double *x0)
    return true;
 }
 
-bool HiopOptimizationProblem::get_vars_info(const long long &n,
+bool HiopOptimizationProblem::get_vars_info(const size_type &n,
                                             double *xlow, double *xupp,
                                             NonlinearityType *type)
 {
@@ -55,7 +55,7 @@ bool HiopOptimizationProblem::get_vars_info(const long long &n,
    return true;
 }
 
-bool HiopOptimizationProblem::get_cons_info(const long long &m,
+bool HiopOptimizationProblem::get_cons_info(const size_type &m,
                                             double *clow, double *cupp,
                                             NonlinearityType *type)
 {
@@ -79,7 +79,7 @@ bool HiopOptimizationProblem::get_cons_info(const long long &m,
    return true;
 }
 
-bool HiopOptimizationProblem::eval_f(const long long &n, const double *x,
+bool HiopOptimizationProblem::eval_f(const size_type &n, const double *x,
                                      bool new_x, double &obj_value)
 {
    MFEM_ASSERT(n == ntdofs_glob, "Global input mismatch.");
@@ -93,7 +93,7 @@ bool HiopOptimizationProblem::eval_f(const long long &n, const double *x,
    return true;
 }
 
-bool HiopOptimizationProblem::eval_grad_f(const long long &n, const double *x,
+bool HiopOptimizationProblem::eval_grad_f(const size_type &n, const double *x,
                                           bool new_x, double *gradf)
 {
    MFEM_ASSERT(n == ntdofs_glob, "Global input mismatch.");
@@ -108,9 +108,9 @@ bool HiopOptimizationProblem::eval_grad_f(const long long &n, const double *x,
    return true;
 }
 
-bool HiopOptimizationProblem::eval_cons(const long long &n, const long long &m,
-                                        const long long &num_cons,
-                                        const long long *idx_cons,
+bool HiopOptimizationProblem::eval_cons(const size_type &n, const size_type &m,
+                                        const size_type &num_cons,
+                                        const index_type *idx_cons,
                                         const double *x, bool new_x,
                                         double *cons)
 {
@@ -134,10 +134,10 @@ bool HiopOptimizationProblem::eval_cons(const long long &n, const long long &m,
    return true;
 }
 
-bool HiopOptimizationProblem::eval_Jac_cons(const long long &n,
-                                            const long long &m,
-                                            const long long &num_cons,
-                                            const long long *idx_cons,
+bool HiopOptimizationProblem::eval_Jac_cons(const size_type &n,
+                                            const size_type &m,
+                                            const size_type &num_cons,
+                                            const index_type *idx_cons,
                                             const double *x, bool new_x,
                                             double *Jac)
 {
@@ -165,16 +165,16 @@ bool HiopOptimizationProblem::eval_Jac_cons(const long long &n,
    return true;
 }
 
-bool HiopOptimizationProblem::get_vecdistrib_info(long long global_n,
-                                                  long long *cols)
+bool HiopOptimizationProblem::get_vecdistrib_info(size_type global_n,
+                                                  index_type *cols)
 {
 #ifdef MFEM_USE_MPI
    int nranks;
    MPI_Comm_size(comm, &nranks);
 
-   long long *sizes = new long long[nranks];
-   MPI_Allgather(&ntdofs_loc, 1, MPI_LONG_LONG_INT, sizes, 1,
-                 MPI_LONG_LONG_INT, comm);
+   size_type *sizes = new size_type[nranks];
+   MPI_Allgather(&ntdofs_loc, 1, MPI_HIOP_SIZE_TYPE, sizes, 1,
+                 MPI_HIOP_SIZE_TYPE, comm);
    cols[0] = 0;
    for (int r = 1; r <= nranks; r++)
    {
@@ -302,6 +302,8 @@ void HiopNlpOptimizer::Mult(const Vector &xt, Vector &x) const
    hiopInstance.options->SetNumericValue("fixed_var_tolerance", 1e-20);
    hiopInstance.options->SetNumericValue("fixed_var_perturb", 1e-9);
 
+   hiopInstance.options->SetNumericValue("mu0", 1e-1);
+
    // 0: no output; 3: not too much
    hiopInstance.options->SetIntegerValue("verbosity_level", print_level);
 

linalg/hiop.hpp

@@ -41,8 +41,8 @@ private:
    const OptimizationProblem &problem;
 
    // Local and global number of variables and constraints.
-   const long long ntdofs_loc, m_total;
-   long long ntdofs_glob;
+   const hiop::size_type ntdofs_loc, m_total;
+   hiop::size_type ntdofs_glob;
 
    // Initial guess.
    const Vector *x_start;
@@ -78,8 +78,7 @@ public:
         constr_vals(m_total), constr_grads(m_total, ntdofs_loc),
         constr_info_is_current(false)
    {
-      MPI_Allreduce(&ntdofs_loc, &ntdofs_glob, 1, MPI_LONG_LONG_INT,
-                    MPI_SUM, comm);
+      MPI_Allreduce(&ntdofs_loc, &ntdofs_glob, 1, MPI_HIOP_SIZE_TYPE, MPI_SUM, comm);
    }
 #endif
 
@@ -87,27 +86,27 @@ public:
 
    /** Extraction of problem dimensions:
     *  n is the number of variables, m is the number of constraints. */
-   virtual bool get_prob_sizes(long long int& n, long long int& m);
+   virtual bool get_prob_sizes(hiop::size_type& n, hiop::size_type& m);
 
    /** Provide an primal starting point. This point is subject to adjustments
     *  internally in HiOp. */
-   virtual bool get_starting_point(const long long &n, double *x0);
+   virtual bool get_starting_point(const hiop::size_type &n, double *x0);
 
-   virtual bool get_vars_info(const long long& n, double *xlow, double* xupp,
+   virtual bool get_vars_info(const hiop::size_type &n, double *xlow, double* xupp,
                               NonlinearityType* type);
 
    /** bounds on the constraints
     *  (clow<=-1e20 means no lower bound, cupp>=1e20 means no upper bound) */
-   virtual bool get_cons_info(const long long &m, double *clow, double *cupp,
+   virtual bool get_cons_info(const hiop::size_type &m, double *clow, double *cupp,
                               NonlinearityType* type);
 
    /** Objective function evaluation.
     *  Each rank returns the global objective value. */
-   virtual bool eval_f(const long long& n, const double *x, bool new_x,
+   virtual bool eval_f(const hiop::size_type &n, const double *x, bool new_x,
                        double& obj_value);
 
    /** Gradient of the objective function (local chunk). */
-   virtual bool eval_grad_f(const long long &n, const double *x, bool new_x,
+   virtual bool eval_grad_f(const hiop::size_type &n, const double *x, bool new_x,
                             double *gradf);
 
    /** Evaluates a subset of the constraints cons(x). The subset is of size
@@ -132,8 +131,9 @@ public:
     *  When MPI enabled, every rank populates cons, since the constraints are
     *  not distributed.
     */
-   virtual bool eval_cons(const long long &n, const long long &m,
-                          const long long &num_cons, const long long *idx_cons,
+   virtual bool eval_cons(const hiop::size_type &n, const hiop::size_type &m,
+                          const hiop::size_type &num_cons,
+                          const hiop::index_type *idx_cons,
                           const double *x, bool new_x, double *cons);
 
    /** Evaluates the Jacobian of the subset of constraints indicated by
@@ -147,9 +147,9 @@ public:
     *  When MPI enabled, each rank computes only the local columns of the
     *  Jacobian, that is the partials with respect to local variables.
     */
-   virtual bool eval_Jac_cons(const long long &n, const long long &m,
-                              const long long &num_cons,
-                              const long long *idx_cons,
+   virtual bool eval_Jac_cons(const hiop::size_type &n, const hiop::size_type &m,
+                              const hiop::size_type &num_cons,
+                              const hiop::index_type *idx_cons,
                               const double *x, bool new_x, double *Jac);
 
    /** Specifies column partitioning for distributed memory vectors.
@@ -158,7 +158,8 @@ public:
     *  Example: for a vector x of 6 entries (globally) on 3 ranks, the uniform
     *  column partitioning is cols=[0,2,4,6].
     */
-   virtual bool get_vecdistrib_info(long long global_n, long long *cols);
+   virtual bool get_vecdistrib_info(hiop::size_type global_n,
+                                    hiop::index_type *cols);
 
 #ifdef MFEM_USE_MPI
    virtual bool get_MPI_comm(MPI_Comm &comm_out)

linalg/kernels.hpp

@@ -378,7 +378,7 @@ void Swap(T &a, T &b)
    b = tmp;
 }
 
-constexpr double Epsilon = std::numeric_limits<double>::epsilon();
+const double Epsilon = std::numeric_limits<double>::epsilon();
 
 /// Utility function used in CalcSingularvalue<3>.
 MFEM_HOST_DEVICE static inline

linalg/sparsemat.cpp

@@ -2449,52 +2449,67 @@ void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
    });
 }
 
+template <bool useFabs>
+static void JacobiDispatch(const Vector &b, const Vector &x0, Vector &x1,
+                           const Memory<int> &I, const Memory<int> &J,
+                           const Memory<double> &A, const int height,
+                           const double sc)
+{
+   const bool useDevice = b.UseDevice() || x0.UseDevice() || x1.UseDevice();
+
+   const auto bp  = b.Read(useDevice);
+   const auto x0p = x0.Read(useDevice);
+   auto       x1p = x1.Write(useDevice);
+
+   const auto Ip = Read(I, height+1, useDevice);
+   const auto Jp = Read(J, J.Capacity(), useDevice);
+   const auto Ap = Read(A, J.Capacity(), useDevice);
+
+   MFEM_FORALL_SWITCH(useDevice, i, height,
+   {
+      double resi = bp[i], norm = 0.0;
+      for (int j = Ip[i]; j < Ip[i+1]; j++)
+      {
+         resi -= Ap[j] * x0p[Jp[j]];
+         if (useFabs)
+         {
+            norm += fabs(Ap[j]);
+         }
+         else
+         {
+            norm += Ap[j];
+         }
+      }
+      if (norm > 0.0)
+      {
+         x1p[i] = x0p[i] + sc * resi / norm;
+      }
+      else
+      {
+         if (useFabs)
+         {
+            MFEM_ABORT_KERNEL("L1 norm of row is zero.");
+         }
+         else
+         {
+            MFEM_ABORT_KERNEL("sum of row is zero.");
+         }
+      }
+   });
+}
+
 void SparseMatrix::Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
                            double sc) const
 {
    MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
-
-   for (int i = 0; i < height; i++)
-   {
-      double resi = b(i), norm = 0.0;
-      for (int j = I[i]; j < I[i+1]; j++)
-      {
-         resi -= A[j] * x0(J[j]);
-         norm += fabs(A[j]);
-      }
-      if (norm > 0.0)
-      {
-         x1(i) = x0(i) + sc * resi / norm;
-      }
-      else
-      {
-         MFEM_ABORT("L1 norm of row " << i << " is zero.");
-      }
-   }
+   JacobiDispatch<true>(b,x0,x1,I,J,A,height,sc);
 }
 
 void SparseMatrix::Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
                            double sc) const
 {
    MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
-
-   for (int i = 0; i < height; i++)
-   {
-      double resi = b(i), sum = 0.0;
-      for (int j = I[i]; j < I[i+1]; j++)
-      {
-         resi -= A[j] * x0(J[j]);
-         sum  += A[j];
-      }
-      if (sum > 0.0)
-      {
-         x1(i) = x0(i) + sc * resi / sum;
-      }
-      else
-      {
-         MFEM_ABORT("sum of row " << i << " is zero.");
-      }
-   }
+   JacobiDispatch<false>(b,x0,x1,I,J,A,height,sc);
 }
 
 void SparseMatrix::AddSubMatrix(const Array<int> &rows, const Array<int> &cols,

makefile

@@ -703,7 +703,7 @@ status info:
 
 ASTYLE_BIN = astyle
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
-ASTYLE_VER = "Artistic Style Version 2.05.1"
+ASTYLE_VER = "Artistic Style Version 3.1"
 FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
 FORMAT_FILES += tests/unit/*.cpp
 UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed

mesh/mesh.cpp

@@ -2850,14 +2850,16 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          y = sfc[3*k + 1];
          z = sfc[3*k + 2];
 
-         ind[0] = VTX(x, y, z  );
-         ind[1] = VTX(x+1, y, z  );
+         // *INDENT-OFF*
+         ind[0] = VTX(x  , y  , z  );
+         ind[1] = VTX(x+1, y  , z  );
          ind[2] = VTX(x+1, y+1, z  );
-         ind[3] = VTX(x, y+1, z  );
-         ind[4] = VTX(x, y, z+1);
-         ind[5] = VTX(x+1, y, z+1);
+         ind[3] = VTX(x  , y+1, z  );
+         ind[4] = VTX(x  , y  , z+1);
+         ind[5] = VTX(x+1, y  , z+1);
          ind[6] = VTX(x+1, y+1, z+1);
-         ind[7] = VTX(x, y+1, z+1);
+         ind[7] = VTX(x  , y+1, z+1);
+         // *INDENT-ON*
 
          AddHex(ind, 1);
       }
@@ -2870,14 +2872,16 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          {
             for (x = 0; x < nx; x++)
             {
-               ind[0] = VTX(x, y, z  );
-               ind[1] = VTX(x+1, y, z  );
+               // *INDENT-OFF*
+               ind[0] = VTX(x  , y  , z  );
+               ind[1] = VTX(x+1, y  , z  );
                ind[2] = VTX(x+1, y+1, z  );
-               ind[3] = VTX(x, y+1, z  );
-               ind[4] = VTX(x, y, z+1);
-               ind[5] = VTX(x+1, y, z+1);
+               ind[3] = VTX(x  , y+1, z  );
+               ind[4] = VTX(x  , y  , z+1);
+               ind[5] = VTX(x+1, y  , z+1);
                ind[6] = VTX(x+1, y+1, z+1);
                ind[7] = VTX(  x, y+1, z+1);
+               // *INDENT-ON*
                if (type == Element::TETRAHEDRON)
                {
                   AddHexAsTets(ind, 1);
@@ -2888,7 +2892,7 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
                }
                else if (type == Element::PYRAMID)
                {
-                  ind[8] = VTXP(  x, y, z);
+                  ind[8] = VTXP(x, y, z);
                   AddHexAsPyramids(ind, 1);
                }
                else
@@ -2906,10 +2910,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (x = 0; x < nx; x++)
       {
-         ind[0] = VTX(x, y, 0);
-         ind[1] = VTX(x, y+1, 0);
+         // *INDENT-OFF*
+         ind[0] = VTX(x  , y  , 0);
+         ind[1] = VTX(x  , y+1, 0);
          ind[2] = VTX(x+1, y+1, 0);
-         ind[3] = VTX(x+1, y, 0);
+         ind[3] = VTX(x+1, y  , 0);
+         // *INDENT-ON*
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 1);
@@ -2929,10 +2935,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (x = 0; x < nx; x++)
       {
-         ind[0] = VTX(x, y, nz);
-         ind[1] = VTX(x+1, y, nz);
+         // *INDENT-OFF*
+         ind[0] = VTX(x  , y  , nz);
+         ind[1] = VTX(x+1, y  , nz);
          ind[2] = VTX(x+1, y+1, nz);
-         ind[3] = VTX(x, y+1, nz);
+         ind[3] = VTX(x  , y+1, nz);
+         // *INDENT-ON*
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 6);
@@ -2952,10 +2960,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (y = 0; y < ny; y++)
       {
-         ind[0] = VTX(0, y, z  );
-         ind[1] = VTX(0, y, z+1);
-         ind[2] = VTX(0, y+1, z+1);
-         ind[3] = VTX(0, y+1, z  );
+         // *INDENT-OFF*
+         ind[0] = VTX(0  , y  , z  );
+         ind[1] = VTX(0  , y  , z+1);
+         ind[2] = VTX(0  , y+1, z+1);
+         ind[3] = VTX(0  , y+1, z  );
+         // *INDENT-ON*
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 5);
@@ -2971,10 +2981,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (y = 0; y < ny; y++)
       {
-         ind[0] = VTX(nx, y, z  );
+         // *INDENT-OFF*
+         ind[0] = VTX(nx, y  , z  );
          ind[1] = VTX(nx, y+1, z  );
          ind[2] = VTX(nx, y+1, z+1);
-         ind[3] = VTX(nx, y, z+1);
+         ind[3] = VTX(nx, y  , z+1);
+         // *INDENT-ON*
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 3);
@@ -2990,10 +3002,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (z = 0; z < nz; z++)
       {
-         ind[0] = VTX(x, 0, z  );
+         // *INDENT-OFF*
+         ind[0] = VTX(x  , 0, z  );
          ind[1] = VTX(x+1, 0, z  );
          ind[2] = VTX(x+1, 0, z+1);
-         ind[3] = VTX(x, 0, z+1);
+         ind[3] = VTX(x  , 0, z+1);
+         // *INDENT-ON*
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 2);
@@ -3009,10 +3023,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (z = 0; z < nz; z++)
       {
-         ind[0] = VTX(x, ny, z  );
-         ind[1] = VTX(x, ny, z+1);
+         // *INDENT-OFF*
+         ind[0] = VTX(x  , ny, z  );
+         ind[1] = VTX(x  , ny, z+1);
          ind[2] = VTX(x+1, ny, z+1);
          ind[3] = VTX(x+1, ny, z  );
+         // *INDENT-ON*
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 4);
@@ -6348,9 +6364,12 @@ void Mesh::GenerateNCFaceInfo()
       const NCMesh::Master &master = list.masters[i];
       if (master.index >= nfaces) { continue; }
 
-      faces_info[master.index].NCFace = nc_faces_info.Size();
+      FaceInfo &master_fi = faces_info[master.index];
+      master_fi.NCFace = nc_faces_info.Size();
       nc_faces_info.Append(NCFaceInfo(false, master.local, NULL));
       // NOTE: one of the unused members stores local face no. to be used below
+      MFEM_ASSERT(master_fi.Elem2No == -1, "internal error");
+      MFEM_ASSERT(master_fi.Elem2Inf == -1, "internal error");
    }
 
    // add records for slave faces

mesh/mesh.hpp

@@ -107,18 +107,23 @@ protected:
    // element (slave or master, in the non-conforming case) that generated the
    // face.
    // Classification of a local (non-ghost) face based on its FaceInfo:
-   // - Elem2No >= 0 --> local internal face; can be either:
+   // - Elem2No >= 0 --> local interior face; can be either:
    //    - NCFace == -1 --> conforming face, or
-   //    - NCFace >= 0 --> non-conforming slave face.
+   //    - NCFace >= 0 --> non-conforming slave face; Elem2No is the index of
+   //      the master volume element; Elem2Inf%64 is 0, see the note in
+   //      Mesh::GenerateNCFaceInfo().
    // - Elem2No < 0 --> local "boundary" face; can be one of:
    //    - NCFace == -1 --> conforming face; can be either:
    //       - Elem2Inf < 0 --> true boundary face (no element on side 2)
    //       - Elem2Inf >= 0 --> shared face where element 2 is a face-neighbor
    //         element with index -1-Elem2No. This state is initialized by
    //         ParMesh::ExchangeFaceNbrData().
-   //    - NCFace >= 0 --> non-conforming master face. Elem2No is -1 or, in the
-   //      case of a shared face, -1-Elem2No is the index of one of the adjacent
-   //      (the last one?) slave ghost elements. Elem2Inf is -1.
+   //    - NCFace >= 0 --> non-conforming face; can be one of:
+   //       - Elem2Inf < 0 --> master non-conforming face, interior or shared;
+   //         In this case, Elem2No is -1; see GenerateNCFaceInfo().
+   //       - Elem2Inf >= 0 --> shared slave non-conforming face where element 2
+   //         is the master face-neighbor element with index -1-Elem2No; see
+   //         ParNCMesh::GetFaceNeighbors().
    //
    // A ghost face is a non-conforming face that is generated by a non-local,
    // i.e. ghost, element. A ghost face has index i in faces_info such that
@@ -130,7 +135,7 @@ protected:
    //   master side element, i.e. side 1 IS NOT the side that generated the
    //   face. Elem2No is < 0 and -1-Elem2No is the index of the ghost
    //   face-neighbor element that generated this slave ghost face. In this
-   //   case, Elem2Inf >= 0.
+   //   case, Elem2Inf >= 0 and NCFace >= 0.
    // Relevant methods: GenerateFaces(), GenerateNCFaceInfo(),
    //                   ParNCMesh::GetFaceNeighbors(),
    //                   ParMesh::ExchangeFaceNbrData()

mesh/mesh_readers.cpp

@@ -1895,6 +1895,9 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
          ho_wdg[2] = wdg18; ho_wdg[3] = wdg40;
          ho_pyr[2] = pyr14; ho_pyr[3] = pyr30;
 
+         bool has_nonpositive_phys_domain = false;
+         bool has_positive_phys_domain = false;
+
          if (binary)
          {
             int n_elem_part = 0; // partial sum of elements that are read
@@ -1945,17 +1948,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                      vert_indices[vi] = it->second;
                   }
 
-                  // non-positive attributes are not allowed in MFEM
+                  // Non-positive attributes are not allowed in MFEM. However,
+                  // by default, Gmsh sets the physical domain of all elements
+                  // to zero. In the case that all elements have physical domain
+                  // zero, we will given them attribute 1. If only some elements
+                  // have physical domain zero, we will throw an error.
                   if (phys_domain <= 0)
                   {
-                     MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
-                                "By default Gmsh sets element tags (attributes)"
-                                " to '0' but MFEM requires that they be"
-                                " positive integers.\n"
-                                "Use \"Physical Curve\", \"Physical Surface\","
-                                " or \"Physical Volume\" to set tags/attributes"
-                                " for all curves, surfaces, or volumes in your"
-                                " Gmsh geometry to values which are >= 1.");
+                     has_nonpositive_phys_domain = true;
+                     phys_domain = 1;
+                  }
+                  else
+                  {
+                     has_positive_phys_domain = true;
                   }
 
                   // initialize the mesh element
@@ -2172,17 +2177,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                   vert_indices[vi] = it->second;
                }
 
-               // non-positive attributes are not allowed in MFEM
+               // Non-positive attributes are not allowed in MFEM. However,
+               // by default, Gmsh sets the physical domain of all elements
+               // to zero. In the case that all elements have physical domain
+               // zero, we will given them attribute 1. If only some elements
+               // have physical domain zero, we will throw an error.
                if (phys_domain <= 0)
                {
-                  MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
-                             "By default Gmsh sets element tags (attributes)"
-                             " to '0' but MFEM requires that they be"
-                             " positive integers.\n"
-                             "Use \"Physical Curve\", \"Physical Surface\","
-                             " or \"Physical Volume\" to set tags/attributes"
-                             " for all curves, surfaces, or volumes in your"
-                             " Gmsh geometry to values which are >= 1.");
+                  has_nonpositive_phys_domain = true;
+                  phys_domain = 1;
+               }
+               else
+               {
+                  has_positive_phys_domain = true;
                }
 
                // initialize the mesh element
@@ -2367,6 +2374,24 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
             } // el (all elements)
          } // if ASCII
 
+         if (has_positive_phys_domain && has_nonpositive_phys_domain)
+         {
+            MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
+                       "By default Gmsh sets element tags (attributes)"
+                       " to '0' but MFEM requires that they be"
+                       " positive integers.\n"
+                       "Use \"Physical Curve\", \"Physical Surface\","
+                       " or \"Physical Volume\" to set tags/attributes"
+                       " for all curves, surfaces, or volumes in your"
+                       " Gmsh geometry to values which are >= 1.");
+         }
+         else if (has_nonpositive_phys_domain)
+         {
+            mfem::out << "\nGmsh reader: all element attributes were zero.\n"
+                      << "MFEM only supports positive element attributes.\n"
+                      << "Setting element attributes to 1.\n\n";
+         }
+
          if (!elements_3D.empty())
          {
             Dim = 3;

mesh/ncmesh.cpp

@@ -612,7 +612,7 @@ int NCMesh::NewSegment(int n0, int n1, int attr, int vattr1, int vattr2)
    // get (degenerate) faces and assign face attributes
    int v0 = el.node[0], v1 = el.node[1];
    faces.Get(v0, v0, v0, v0)->attribute = vattr1;
-   faces.Get(v1, v1, v1,v1)->attribute = vattr2;
+   faces.Get(v1, v1, v1, v1)->attribute = vattr2;
 
    return new_id;
 }

miniapps/meshing/minimal-surface.cpp

@@ -1159,8 +1159,8 @@ static double u0(const Vector &x) { return sin(3.0 * PI * (x[1] + x[0])); }
 
 enum {NORM, AREA};
 
-double qf(const int order, const int ker, Mesh &m,
-          FiniteElementSpace &fes, GridFunction &u)
+static double qf(const int order, const int ker, Mesh &m,
+                 FiniteElementSpace &fes, GridFunction &u)
 {
    const Geometry::Type type = m.GetElementBaseGeometry(0);
    const IntegrationRule &ir(IntRules.Get(type, order));

tests/convergence/prates.cpp

@@ -117,7 +117,7 @@ int main(int argc, char *argv[])
       MPI_Finalize();
       return 1;
    }
-   if (prob >3 || prob <0) { prob = 0; } // default problem = H1
+   if (prob >3 || prob <0) prob = 0; // default problem = H1
    if (prob == 3)
    {
       if (kappa < 0)
@@ -221,7 +221,7 @@ int main(int argc, char *argv[])
          gradu = new VectorFunctionCoefficient(dim,gradu_exact);
          b.AddDomainIntegrator(new DomainLFIntegrator(*f));
          b.AddBdrFaceIntegrator(
-            new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
+      new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
          a.AddDomainIntegrator(new DiffusionIntegrator(one));
          a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
          a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
@@ -291,8 +291,8 @@ int main(int argc, char *argv[])
 
       x = *X;
       JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
-                     : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
-                     : JumpScaling::CONSTANT);
+                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                        : JumpScaling::CONSTANT);
       switch (prob)
       {
          case 0: rates.AddH1GridFunction(&x,scalar_u,gradu); break;
@@ -305,7 +305,7 @@ int main(int argc, char *argv[])
       delete B;
       delete A;
 
-      if (l==pr) { break; }
+      if (l==pr) break;
 
       pmesh->UniformRefinement();
       fespace->Update();

tests/convergence/rates.cpp

@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
       args.PrintUsage(cout);
       return 1;
    }
-   if (prob >3 || prob <0) { prob = 0; } // default problem = H1
+   if (prob >3 || prob <0) prob = 0; // default problem = H1
    if (prob == 3)
    {
       if (kappa < 0)
@@ -194,7 +194,7 @@ int main(int argc, char *argv[])
          gradu = new VectorFunctionCoefficient(dim,gradu_exact);
          b.AddDomainIntegrator(new DomainLFIntegrator(*f));
          b.AddBdrFaceIntegrator(
-            new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
+      new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
          a.AddDomainIntegrator(new DiffusionIntegrator(one));
          a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
          a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
@@ -224,8 +224,8 @@ int main(int argc, char *argv[])
       }
 
       JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
-                     : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
-                     : JumpScaling::CONSTANT);
+                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                        : JumpScaling::CONSTANT);
 
       switch (prob)
       {
@@ -235,7 +235,7 @@ int main(int argc, char *argv[])
          case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one,js); break;
       }
 
-      if (l==sr) { break; }
+      if (l==sr) break;
 
       mesh->UniformRefinement();
       fespace->Update();

tests/gitlab/README.md

@@ -15,8 +15,9 @@ This directory contains utility scripts related to Gitlab testing at LLNL.
    While designed to be used in CI context, this script can also be used
    standalone on LLNL's LC in order to reproduce a similar build. The script
    uses environment variables for configuration (a place for improvement), such
-   as "BUILD_ROOT", "HOST_CONFIG", "SPEC", etc. Some are mandatory, while others
-   have default values.
+   as "BUILD_ROOT", "HOST_CONFIG", "SPEC", etc. Some are mandatory, while
+   others have default values.
+   Please refer to tests/gitlab/reproduce-ci-jobs-interactively.md for details.
 
 * `get_mfem_uberenv` sets uberenv up for use with MFEM, notably to install TPLs
   with Spack.
@@ -27,8 +28,8 @@ This directory contains utility scripts related to Gitlab testing at LLNL.
 
 * `generate_spack_upstream` can be used to generate a Spack upstream instance.
 
-  This script addresses a much less common use case, where the TPLs for any MFEM
-  target appearing in CI are built in Spack instance using Uberenv. Configuring
-  permissions accordingly allows this instance to be reused to prevent multiple
-  installations. Such an upstream instance can be found on LC in
-  `/usr/workspace/mfem/spack-upstream`.
+  This script addresses a much less common use case, where the TPLs for any
+  MFEM target appearing in CI are built in Spack instance using Uberenv.
+  Configuring permissions accordingly allows this instance to be reused to
+  prevent multiple installations. Such an upstream instance can be found on LC
+  in `/usr/workspace/mfem/spack-upstream`.

tests/gitlab/build_and_test

@@ -14,21 +14,103 @@
 set -o errexit
 set -o nounset
 
-option=${1:-""}
+script_name=$(basename $0)
+
+function usage()
+{
+  echo "MFEM CI build and test driver"
+  echo ""
+  echo "Syntax:"
+  echo "> ${script_name} --spec \"spack spec\" [--deps-only] [--data]"
+  echo "                 [--build-root /path/to/build/resource]"
+  echo ""
+  echo "> ${script_name} --build-only [--data]"
+  echo "                 [--build-root /path/to/build/resource]"
+  echo ""
+  echo "> ${script_name} --test-only [--data]"
+  echo "                 [--build-root /path/to/build/resource]"
+  echo ""
+  echo "Options:"
+  echo " --spec"
+  echo "    The most common usage is to pass a spack spec to the script."
+  echo "    This spec may or may not specify \"mfem\" as the package, since it"
+  echo "    is the default."
+  echo ""
+  echo " --deps-only"
+  echo "    Only install MFEM dependencies. Will clone mfem-uberenv to drive a"
+  echo "    local instance on spack with which we install the dependencies."
+  echo "    This local spack instance is tuned to be completely independent"
+  echo "    from the user environment."
+  echo ""
+  echo " --build-only"
+  echo "    Only build MFEM, based on the existing configuration."
+  echo ""
+  echo " --test-only"
+  echo "    Only test MFEM, based on the existing configuration."
+  echo ""
+  echo " --data"
+  echo "    Make running tests with \"--data\" a requirement: will fail in the"
+  echo "    data directory is not present in the parent of the mfem root directory."
+  echo "    Note: default behavior is to run data tests if data dir is present."
+  echo ""
+  echo " --build-root=/path/to/build/resource"
+  echo "    The script will use this directory to find the external resource"
+  echo "    needed, e.g. the data directory. Defaults to the parent location"
+  echo "    of the MFEM clone."
+  echo ""
+}
+
 hostname="$(hostname)"
 project_dir="$(pwd)"
 
-build_root=${BUILD_ROOT:-""}
-hostconfig=${HOST_CONFIG:-""}
-spec=${SPEC:-""}
+mode=""
+build_root=""
+spec=""
+with_data=false
 
 sys_type=${SYS_TYPE:-""}
 py_env_path=${PYTHON_ENVIRONMENT_PATH:-""}
+ci_context=${CI:-""}
 
 threads=${THREADS:-""}
 
+# Options
+while [[ $# -gt 0 ]]
+do
+    key="$1"
+
+    case $key in
+        --spec)
+            spec="$2"
+            shift # past argument
+            shift # past value
+            ;;
+        --deps-only|--build-only|--test-only)
+            mode="$key"
+            shift # past argument
+            ;;
+        --data)
+            with_data=true
+            shift # past argument
+            ;;
+        --build-root)
+            build_root="$2"
+            shift # past argument
+            shift # past value
+            ;;
+        -h|--help)
+            usage
+            exit 0;
+            ;;
+        *)  # unknown option
+            echo "ERROR: option $key is unknown"
+            exit 1;
+            ;;
+    esac
+done
+
 # Dependencies
-if [[ "${option}" != "--build-only" && "${option}" != "--test-only" ]]
+if [[ "${mode}" != "--build-only" && "${mode}" != "--test-only" ]]
 then
     echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
     echo "~~~~~ Building Dependencies"
@@ -36,19 +118,21 @@ then
 
     if [[ -z ${spec} ]]
     then
-        echo "SPEC is undefined, aborting..."
+        echo "--spec is mandatory to build dependencies, aborting..."
         exit 1
     fi
 
     upstream_opt=""
     if [[ ${CI_COMMIT_BRANCH:-"master"} != "master" && ${sys_type} != "" ]]
     then
-        upstream_opt="--upstream=/usr/workspace/mfem/spack-upstream/uberenv_libs"
+        upstream_opt="--upstream=/usr/workspace/mfem/mfem-dependencies/install"
     fi
 
     prefix_opt=""
 
-    if [[ -d "/dev/shm" ]]
+    # Do not run on /dev/shm if only installing dependencies: they will be lost
+    # otherwise.
+    if [[ -d "/dev/shm" && "${mode}" != "--deps-only" ]]
     then
         prefix="/dev/shm/${hostname}/${CI_PIPELINE_ID:-"NONE"}_${RANDOM}"
         mkdir -p ${prefix}
@@ -59,66 +143,69 @@ then
     echo "Fetching uberenv."
     tests/gitlab/get_mfem_uberenv || ( echo "Error fetching Uberenv" && exit 1 );
 
+    echo "Removing existing configuration"
+    make distclean
+
     echo "Installing dependencies."
     python tests/uberenv/uberenv.py --spec="${spec}" "${upstream_opt}" "${prefix_opt}"
 
+    # Make sure that a configuration was generated by spack (part 1).
+    cp config/config.mk config/spack-config.mk
+    cp config/_config.hpp config/spack_config.hpp
 fi
 
-# Host config file
-if [[ -z ${hostconfig} ]]
+# Configuration
+if [[ "${mode}" != "--deps-only" ]]
 then
-    # If no host config file was provided, we assume it was generated.
-    # This means we are looking of a unique one in project dir.
-    hostconfigs=( $( ls "${project_dir}/"hc-*.mk ) )
-    if [[ ${#hostconfigs[@]} == 1 ]]
+# Host config file
+# Make sure that a configuration was generated by spack (part 2).
+# Without this running the rest of the script would generate a new config.
+    if [[ -f "config/spack-config.mk" && -f "config/spack_config.hpp" ]]
     then
-        hostconfig_path=${hostconfigs[0]}
-        echo "Found host config file: ${hostconfig_path}"
-    elif [[ ${#hostconfigs[@]} == 0 ]]
-    then
-        echo "No result for: ${project_dir}/hc-*.mk"
-        echo "Spack generated host-config not found."
-        exit 1
+        cp config/spack-config.mk config/config.mk
+        cp config/spack_config.hpp config/_config.hpp
     else
-        echo "More than one result for: ${project_dir}/hc-*.mk"
-        echo "${hostconfigs[@]}"
-        echo "Please specify one with HOST_CONFIG variable"
+        echo "No result for at least one of"
+        echo "    ${project_dir}/config/spack-config.mk"
+        echo "    ${project_dir}/config/spack_config.hpp"
+        echo "ERROR: Spack generated configuration not found."
         exit 1
     fi
-else
-    # Using provided host-config file.
-    hostconfig_path="${project_dir}/host-configs/${hostconfig}"
+
+    # Build and Data Directories
+    if [[ -z ${build_root} ]]
+    then
+        # By default, build_root is the project parent dir.
+        build_root="${project_dir}/.."
+    else
+        # build_root is specified, so we need to link its content into the
+        # project parent dir.
+        ln -sf ${build_root}/data ${project_dir}/../
+    fi
+
+    if [[ "$with_data" == "true" && ! -d ${build_root}/data ]]
+    then
+        echo "ERROR: ${build_root}/data not found while asking for --data".
+    fi
 fi
 
-# Build Directory
-if [[ -z ${build_root} ]]
-then
-    build_root=$(pwd)
-fi
-
-build_dir="${build_root}/build_${hostconfig//.mk/}"
-
-# Build
-if [[ "${option}" != "--deps-only" && "${option}" != "--test-only" ]]
+# Build (also build when asked for testing, to make sure we test what we want)
+if [[ "${mode}" != "--deps-only" ]]
 then
     echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
-    echo "~ Host-config: ${hostconfig_path}"
-    echo "~ Build Dir:   ${build_dir}"
     echo "~ Project Dir: ${project_dir}"
+    echo "~ Build Root: ${build_root}"
     echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
 
     echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
     echo "~~~~~ Building MFEM"
     echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
 
-    cp ${hostconfig_path} ${project_dir}/config/
-    ln -sf ${build_root}/data ${project_dir}/../
-
     make all -j ${threads}
 fi
 
 # Test
-if [[ "${option}" != "--build-only" ]]
+if [[ "${mode}" != "--deps-only" && "${mode}" != "--build-only" ]]
 then
     echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
     echo "~~~~~ Testing MFEM"

tests/gitlab/generate_spack_upstream

@@ -21,11 +21,14 @@ cd tests/uberenv
 git apply upstream-permission.patch
 cd -
 
+# Get the hostname
+hostid=${HOSTNAME//[[:digit:]]/}
+
 # call uberenv for all specs in CI
-git grep -e "^[^#]" .gitlab | grep "SPEC" \
-                            | cut -d' ' -f6- \
-                            | sed 's/"//g' \
-                            | while read -r line; do
+git grep -e "^[^#]" .gitlab/$hostid* | grep "SPEC" \
+                                     | cut -d' ' -f6- \
+                                     | sed 's/"//g' \
+                                     | while read -r line; do
   python ./tests/uberenv/uberenv.py --spec="$line"
 done
 

tests/gitlab/get_mfem_uberenv

@@ -18,7 +18,7 @@ set -o errexit
 set -o nounset
 
 uberenv_url="https://github.com/mfem/mfem-uberenv.git"
-uberenv_ref="222e3743377d8c903e044f80ae44624e1668feb4"
+uberenv_ref="a268844932472ff9bec9caaee449ffecd4bdf2f3"
 
 [[ ! -d tests/uberenv ]] && git clone ${uberenv_url} tests/uberenv
 cd tests/uberenv

tests/gitlab/reproduce-ci-jobs-interactively.md

@@ -0,0 +1,121 @@
+# HowTo: Reproduce CI jobs interactively.
+
+We rely on Spack, driven by uberenv, to build MFEM dependencies automatically
+in CI. Then we use a script that will be retrieve the configuration files
+generated and use it to build MFEM before running the tests.
+
+This process is straightforward to reproduce. However the steps are not easy
+to extract from the CI configuration, hence this article.
+
+**WARNING**
+
+This will only work on the same machines as used for CI. To add a new machine
+you will need to add a corresponding configuration in mfem-uberenv repo.
+
+**NOTE**
+
+The `build_and_test` script controlling the build of MFEM and its dependencies
+has two modes:
+- The CI mode will build the deps, then MFEM and then run tests for the
+  specified spec.
+- The Interactive mode, presented here, where we first build the dependencies
+  with the `--deps-only` option and then use the spack generated configuration
+  files to build MFEM.
+
+Those modes are essentially the same, but we emphasize building the
+dependencies as a first isolated step because it makes it clear what is
+happening and how to use this workflow.
+
+## Prerequisite: Retrieve Uberenv
+
+```bash
+tests/gitlab/get_mfem_uberenv
+```
+
+We have a script to automatically download and install an MFEM-tuned uberenv
+instance.
+
+Uberenv will be placed in `tests/uberenv`, along with spack configuration files
+for each machine we currently support, and possibly some Spack packages we
+patch.
+
+
+## Install Dependencies
+
+### Option #1: Using the CI script
+
+```bash
+./tests/gitlab/build_and_test --deps-only --spec "%gcc@6.1.0 +sundials"
+```
+
+The CI script has three steps that can be run individually with the options
+`--deps-only`, `--build-only`, `--test-only`.
+
+We ask to build only the dependencies, and we need to provide a spack spec
+through the `--spec` option.
+
+Virtually any spec can be provided, but you should check which compilers are
+defined in the spack configuration
+(`tests/uberenv/spack-configs/<sys_type>/compilers.yaml).
+
+As a result, dependencies will be installed under `uberenv_libs`.  The
+configuration files `spack-config.mk` and `spack_config.hpp` will be generated
+in `config`. They are just copies of the `config.mk` and `\_config.hpp`
+generated by the MFEM spack package. This copy is meant to make sure in CI that
+those files were generated, otherwise `make all` would just regenerate them.
+
+**NOTE**
+
+The `build_and_test` script behaves slightly differently between CI context and
+elsewhere (depending on environment variable $CI). In CI, and if launched on
+quartz, ruby or corona, the script will build and install dependencies in
+`/dev/shm` for better performance. However, this is only valid if we don’t want
+the installation to persist. Installation will happen locally to the uberenv
+directory if not in CI context.
+
+### Option #2: Calling uberenv directly
+
+```bash
+python ./tests/uberenv/uberenv.py --spec="%gcc@6.1.0 +sundials"
+```
+
+This is essentially the command the CI script runs in the end.
+
+**NOTE**
+
+When using this command, the configuration files will not be copied (saved),
+but still ready to use.
+
+## Build and test MFEM
+
+### Option #1: Without using scripts
+
+```bash
+\# optional: the spack config is saved in copies (in case you reconfigure MFEM
+\# in the meantime).
+cp config/spack-config.mk config/config.mk
+cp config/spack_config.hpp config/_config.hpp
+```
+
+```bash
+make all -j 8
+make test
+```
+
+The key here is that the configuration files were generated during the spack
+run: they will point to the dependencies location, and apply any option
+selected with MFEM spec variants in Spack.
+
+### Option #2: Using the CI script
+
+```bash
+./tests/gitlab/build_and_test --build-only
+```
+
+We could use `--test-only` option, which would also build MFEM to make
+sure to use the spack generated configuration file.
+
+**NOTE**
+
+Using the CI script requires `config/spack-config.mk` and
+`config/spack_config.hpp` to exist.

tests/scripts/code-style

@@ -22,7 +22,7 @@ case $1 in
           -h|-help    prints this usage information and exits
 
    This script checks if the code in mfem_dir conforms to the MFEM C++ style
-   (currently based on astyle-2.05.1). Note that the comparison is done with
+   (currently based on astyle-3.1). Note that the comparison is done with
    "git diff", so all differences between mfem_dir and the repository will be
    reported. The local changes can be erased (lost!) with "git checkout .".
 
@@ -74,7 +74,7 @@ if [ -s "$test_name.err" ]; then
 
 To correct this error, run 'make style' and make sure that:
 
-1) You are using Artistic Style 2.05.1 ("astyle --version").
+1) You are using Artistic Style 3.1 ("astyle --version").
 2) There is no use of std::cout or std::cerr inside the library code (examples
    and miniapps are OK). Use mfem::out and mfem::err instead.
 3) Run "make style" to apply astyle formatting and verify your code.

tests/unit/fem/test_datacollection.cpp

@@ -306,7 +306,7 @@ TEST_CASE("ParaView restart mode", "[ParaView]")
    VerifyDataset(dataset, 1.0);
    dataset = dataset->NextSiblingElement();
    VerifyDataset(dataset, 1.5);
-   REQUIRE(dataset->NextSiblingElement() == NULL);
+   REQUIRE(dataset->NextSiblingElement() == nullptr);
 
    // Clean up
    for (int c=0; c<=2; ++c)

tests/unit/fem/test_inversetransform.cpp

@@ -60,7 +60,7 @@ TEST_CASE("InverseElementTransformation",
    Mesh mesh( meshStr );
 
    REQUIRE( mesh.GetNE() == 1 );
-   REQUIRE( mesh.GetNodes() != NULL );
+   REQUIRE( mesh.GetNodes() != nullptr );
 
    // Optionally, dump mesh to disk
    bool dumpMesh = false;

tests/unit/miniapps/test_sedov.cpp

@@ -1461,12 +1461,12 @@ void QKernel(const int nzones,
             MFEM_FOREACH_THREAD(qy,y,Q1D)
             {
                QBody<dim>(nzones, z, nqp, qx + qy * Q1D,
-               gamma, use_viscosity, h0, h1order, cfl, infinity,
-               Jinv,stress,sgrad_v,eig_val_data,eig_vec_data,
-               compr_dir,Jpi,ph_dir,stressJiT,
-               d_weights, d_Jacobians, d_rho0DetJ0w,
-               d_e_quads, d_grad_v_ext, d_Jac0inv,
-               d_dt_est, d_stressJinvT);
+                          gamma, use_viscosity, h0, h1order, cfl, infinity,
+                          Jinv,stress,sgrad_v,eig_val_data,eig_vec_data,
+                          compr_dir,Jpi,ph_dir,stressJiT,
+                          d_weights, d_Jacobians, d_rho0DetJ0w,
+                          d_e_quads, d_grad_v_ext, d_Jac0inv,
+                          d_dt_est, d_stressJinvT);
             }
          }
          MFEM_SYNC_THREAD;
@@ -1494,12 +1494,12 @@ void QKernel(const int nzones,
                MFEM_FOREACH_THREAD(qz,z,Q1D)
                {
                   QBody<dim>(nzones, z, nqp, qx + Q1D * (qy + qz * Q1D),
-                  gamma, use_viscosity, h0, h1order, cfl, infinity,
-                  Jinv,stress,sgrad_v,eig_val_data,eig_vec_data,
-                  compr_dir,Jpi,ph_dir,stressJiT,
-                  d_weights, d_Jacobians, d_rho0DetJ0w,
-                  d_e_quads, d_grad_v_ext, d_Jac0inv,
-                  d_dt_est, d_stressJinvT);
+                             gamma, use_viscosity, h0, h1order, cfl, infinity,
+                             Jinv,stress,sgrad_v,eig_val_data,eig_vec_data,
+                             compr_dir,Jpi,ph_dir,stressJiT,
+                             d_weights, d_Jacobians, d_rho0DetJ0w,
+                             d_e_quads, d_grad_v_ext, d_Jac0inv,
+                             d_dt_est, d_stressJinvT);
                }
             }
          }