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.")
@@ -428,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 "")
@@ -496,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/ex25.cpp

@@ -92,7 +92,7 @@ class PMLDiagMatrixCoefficient : public VectorCoefficient
 {
 private:
    CartesianPML * pml = nullptr;
-   void (*Function)(const Vector &, CartesianPML * , Vector &);
+   void (*Function)(const Vector &, CartesianPML *, Vector &);
 public:
    PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
                                               Vector &),

examples/ex25p.cpp

@@ -92,7 +92,7 @@ class PMLDiagMatrixCoefficient : public VectorCoefficient
 {
 private:
    CartesianPML * pml = nullptr;
-   void (*Function)(const Vector &, CartesianPML * , Vector &);
+   void (*Function)(const Vector &, CartesianPML *, Vector &);
 public:
    PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
                                               Vector &),

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

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

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/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_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/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

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

linalg/sundials.cpp

@@ -1677,7 +1677,7 @@ int KINSolver::GradientMult(N_Vector v, N_Vector Jv, N_Vector u,
 
 // Wrapper for evaluating linear systems J u = b
 int KINSolver::LinSysSetup(N_Vector u, N_Vector, SUNMatrix J,
-                           void *, N_Vector , N_Vector )
+                           void *, N_Vector, N_Vector )
 {
    const SundialsNVector mfem_u(u);
    KINSolver *self = static_cast<KINSolver*>(GET_CONTENT(J));

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,6 +2850,7 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          y = sfc[3*k + 1];
          z = sfc[3*k + 2];
 
+         // *INDENT-OFF*
          ind[0] = VTX(x  , y  , z  );
          ind[1] = VTX(x+1, y  , z  );
          ind[2] = VTX(x+1, y+1, z  );
@@ -2858,6 +2859,7 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          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*
 
          AddHex(ind, 1);
       }
@@ -2870,6 +2872,7 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          {
             for (x = 0; x < nx; x++)
             {
+               // *INDENT-OFF*
                ind[0] = VTX(x  , y  , z  );
                ind[1] = VTX(x+1, y  , z  );
                ind[2] = VTX(x+1, y+1, z  );
@@ -2878,6 +2881,7 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
                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++)
       {
+         // *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);
+         // *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++)
       {
+         // *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);
+         // *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++)
       {
+         // *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++)
       {
+         // *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);
+         // *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++)
       {
+         // *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);
+         // *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++)
       {
+         // *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/shifted/dist_solver.cpp

@@ -483,11 +483,11 @@ void ScreenedPoisson::AssembleElementVector(const FiniteElement &el,
       pval=shapef*elfun;
       if (fval>0.0)
       {
-         elvect.Add( -w , shapef);
+         elvect.Add( -w, shapef);
       }
       else if (fval<0.0)
       {
-         elvect.Add(  w , shapef);
+         elvect.Add(  w, shapef);
       }
    }
 }
@@ -523,7 +523,7 @@ void ScreenedPoisson::AssembleElementGrad(const FiniteElement &el,
       el.CalcPhysDShape(trans, B);
       el.CalcPhysShape(trans,shapef);
 
-      AddMult_a_VVt(w , shapef, elmat);
+      AddMult_a_VVt(w, shapef, elmat);
       AddMult_a_AAt(w * diffcoef, B, elmat);
    }
 }
@@ -683,11 +683,11 @@ void PUMPLaplacian::AssembleElementVector(const FiniteElement &el,
       // add the external load -1 if tval > 0.0; 1 if tval < 0.0;
       if (tval>0.0)
       {
-         elvect.Add( -w*fval , shapef);
+         elvect.Add( -w*fval, shapef);
       }
       else  if (tval<0.0)
       {
-         elvect.Add(  w*fval , shapef);
+         elvect.Add(  w*fval, shapef);
       }
    }
 }

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_doftrans.cpp

@@ -98,7 +98,7 @@ TEST_CASE("DoF Transformation Classes",
          double uAv = A.InnerProduct(v, u);
 
          REQUIRE(fabs(uAv -  At.InnerProduct(vt, u )) < tol * fabs(uAv));
-         REQUIRE(fabs(uAv -  tA.InnerProduct(v , ut)) < tol * fabs(uAv));
+         REQUIRE(fabs(uAv -  tA.InnerProduct(v, ut)) < tol * fabs(uAv));
          REQUIRE(fabs(uAv - tAt.InnerProduct(vt, ut)) < tol * fabs(uAv));
       }
       SECTION("Inner product of a primal vector and a dual vector")
@@ -119,7 +119,7 @@ TEST_CASE("DoF Transformation Classes",
          double fAv = A.InnerProduct(v, f);
 
          REQUIRE(fabs(fAv -  At.InnerProduct(vt, f )) < tol * fabs(fAv));
-         REQUIRE(fabs(fAv -  tA.InnerProduct(v , ft)) < tol * fabs(fAv));
+         REQUIRE(fabs(fAv -  tA.InnerProduct(v, ft)) < tol * fabs(fAv));
          REQUIRE(fabs(fAv - tAt.InnerProduct(vt, ft)) < tol * fabs(fAv));
       }
    }
@@ -185,9 +185,9 @@ TEST_CASE("DoF Transformation Functions",
 
       double fAv = A.InnerProduct(v, f);
 
-      REQUIRE(fabs(fAv - nAn.InnerProduct(v , f )) < tol * fabs(fAv));
+      REQUIRE(fabs(fAv - nAn.InnerProduct(v, f )) < tol * fabs(fAv));
       REQUIRE(fabs(fAv -  At.InnerProduct(vt, f )) < tol * fabs(fAv));
-      REQUIRE(fabs(fAv -  tA.InnerProduct(v , ft)) < tol * fabs(fAv));
+      REQUIRE(fabs(fAv -  tA.InnerProduct(v, ft)) < tol * fabs(fAv));
       REQUIRE(fabs(fAv - tAt.InnerProduct(vt, ft)) < tol * fabs(fAv));
    }
    SECTION("TransformDual")
@@ -217,9 +217,9 @@ TEST_CASE("DoF Transformation Functions",
 
       double uAv = A.InnerProduct(v, u);
 
-      REQUIRE(fabs(uAv - nAn.InnerProduct(v , u )) < tol * fabs(uAv));
+      REQUIRE(fabs(uAv - nAn.InnerProduct(v, u )) < tol * fabs(uAv));
       REQUIRE(fabs(uAv -  At.InnerProduct(vt, u )) < tol * fabs(uAv));
-      REQUIRE(fabs(uAv -  tA.InnerProduct(v , ut)) < tol * fabs(uAv));
+      REQUIRE(fabs(uAv -  tA.InnerProduct(v, ut)) < tol * fabs(uAv));
       REQUIRE(fabs(uAv - tAt.InnerProduct(vt, ut)) < tol * fabs(uAv));
    }
 }

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