Merge remote-tracking branch 'origin/master' into ierus-dev

cmake: install examples when enabled

.binder/environment.yml

@@ -0,0 +1,12 @@
+# extends https://github.com/jupyterhub/repo2docker/blob/main/repo2docker/buildpacks/conda/environment.yml
+# see https://mybinder.readthedocs.io/en/latest/using/config_files.html#environment-yml-install-a-conda-environment
+channels:
+- conda-forge
+dependencies:
+- xeus-cling=0.13.0
+- xwidgets=0.26.0
+# NOTE: it's possible these aren't needed for the lab frontend
+- widgetsnbextension=3.5.1
+- pip
+- pip:
+  - glvis==0.3.2

.binder/postBuild

@@ -0,0 +1,26 @@
+#!/bin/bash
+
+set -e
+
+# cling is installed here (in bin) and will look in {dir}/include and {dir}/lib
+# without extra intervention (jk it doesn't look in {dir}/lib unless something
+# has been #included from {dir}/include first...)
+install_dir=/srv/conda/envs/notebook
+mkdir -p $install_dir
+
+# build and install mfem, which is the directory we start in
+make serial SHARED=YES -j8
+make install PREFIX=$install_dir
+
+# install xeus-glvis
+git clone https://github.com/GLVis/xeus-glvis.git
+pushd xeus-glvis
+make install prefix=$install_dir
+popd
+
+# install jupyter-lab extension
+jupyter labextension install @jupyter-widgets/jupyterlab-manager --no-build
+jupyter labextension install glvis-jupyter
+
+# fixup kernelspec, we could probably do this from sh but ¯\_(ツ)_/¯
+python .binder/update_kernel_env.py

.binder/update_kernel_env.py

@@ -0,0 +1,14 @@
+# Update the LD_LIBRARY_PATH of the C++14 kernel so it can find mfem without
+# extra pragma cling statements
+
+import json
+
+kernelspec = "/srv/conda/envs/notebook/share/jupyter/kernels/xcpp14/kernel.json"
+
+with open(kernelspec, "r") as f:
+    obj = json.load(f)
+
+obj["env"] = {"LD_LIBRARY_PATH": "/srv/conda/envs/notebook/lib"}
+
+with open(kernelspec, "w") as f:
+    json.dump(obj, f)

.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: |
@@ -105,6 +105,9 @@ jobs:
 
     - name: branch-history
       run: |
-        git fetch origin master:master
+        # We override origin to make sure we point to the main repo.
+        # This is to have consistent test results on PRs from forks.
+        git remote remove origin
+        git remote add origin https://github.com/mfem/mfem.git
         git checkout -b gh-actions-branch-history
         ./config/githooks/pre-push --history

.gitignore

@@ -51,6 +51,8 @@ examples/ex1[04-9]
 examples/ex1[0-9]p
 examples/ex2[0-9]
 examples/ex2[0-9]p
+examples/ex30
+examples/ex30p
 
 examples/refined.mesh
 examples/displaced.mesh
@@ -288,10 +290,15 @@ miniapps/solvers/ParaView
 miniapps/solvers/mesh.*
 miniapps/solvers/sol.*
 
+miniapps/parelag/MultilevelHcurlHdivSolver
+miniapps/parelag/*.mesh
+
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests
 tests/unit/punit_tests
+tests/unit/cunit_tests
+tests/unit/pcunit_tests
 tests/unit/sedov_tests_*
 tests/unit/psedov_tests_*
 tests/unit/tmop_pa_tests_*
@@ -299,6 +306,12 @@ tests/unit/ptmop_pa_tests_*
 tests/unit/ceed_tests
 tests/unit/debug_device_tests
 
+# Benchmark binaries
+tests/benchmarks/bench_ceed
+tests/benchmarks/bench_tmop
+tests/benchmarks/bench_vector
+tests/benchmarks/bench_virtuals
+
 # Test script output
 tests/scripts/*.err
 tests/scripts/*.out
@@ -315,3 +328,6 @@ build-*/*
 # PETSc automated build
 petsc-build/*
 pkg.gitcommit
+
+# Jupyter Notebook Checkpoints
+.ipynb_checkpoints

.gitlab-ci.yml

@@ -13,242 +13,52 @@
 # at Lawrence Livermore National Laboratory (LLNL). This entire pipeline is
 # LLNL-specific!
 
-# We define the following GitLab pipeline variables:
-#
-# BUILD_ROOT:
-# The path to the shared resources 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
-# also means that the BUILD_ROOT directory will never be cleaned.
-# TODO: add a clean-up mechanism
-#
-# REBASELINE:
-# Defines the default choice for updating the saved baseline results. By default
-# the baseline can only be updated from the master branch. This variable offers
-# the option to manually ask for rebaselining from another branch if necessary.
-#
-# 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
-# are sure to retrieve it.
-#
-# TPLS_REPO & TESTS_REPO:
-# Git repositories used in the pipeline
-#
-# ARTIFACTS_DIR:
-# Directory used to place artifacts.
-
-variables:
-  BUILD_ROOT: ${CI_BUILDS_DIR}/MFEM/${CI_PROJECT_NAME}_${CI_COMMIT_REF_SLUG}_${CI_PIPELINE_ID}
-  AUTOTEST_ROOT: ${CI_BUILDS_DIR}/MFEM
-  REBASELINE: "NO"
-  AUTOTEST: "NO"
-  ALLOC_NAME: ${CI_PROJECT_NAME}_ci_${CI_PIPELINE_ID}
-  TPLS_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/tpls.git
-  TESTS_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/tests.git
-  AUTOTEST_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/autotest.git
-  MFEM_DATA_REPO: https://github.com/mfem/data.git
-  ARTIFACTS_DIR: artifacts
 
 # The pipeline is divided into stages. Usually, jobs in a given stage wait for
 # the preceding stages to complete before to start. However, we sometimes use
 # 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_build_and_test
-  - q_release_resources
-  - l_build_and_test
-  - c_build_and_test
-  - setup_baseline
-  - baseline_check
-  - baseline_to_autotest
-  - baseline_publish
+  - sub-pipelines
 
-# setup clones the mfem/data repo in ${BUILD_ROOT}. The build_and_test script
-# then symlinks the repo to the parent directory of the MFEM source directory.
-# Unit tests that depend on the mfem/data repo will then detect that this
-# directory is present and be enabled.
-setup:
-  tags:
-    - shell
-    - quartz
-  stage: setup
+variables:
+  CUSTOM_CI_BUILDS_DIR: "/usr/workspace/mfem/gitlab-runner"
+
+# Trigger subpipelines:
+quartz-build-and-test:
+  stage: sub-pipelines
   variables:
-    GIT_STRATEGY: none
-  script:
-    - mkdir -p ${BUILD_ROOT} && cd ${BUILD_ROOT}
-    - if [ ! -d data ]; then git clone ${MFEM_DATA_REPO}; fi
+    _AUTOTEST: $AUTOTEST
+  trigger:
+    include: .gitlab/quartz-build-and-test.yml
+    strategy: depend
 
-# The setup_baseline job in setup stage_baseline doesn't rely on MFEM git repo.
-# It prepares a pipeline-wide working directory downloading/updating external
-# repos. TODO: updating tests and tpls is not necessary anymore since pipelines
-# are now using unique directories so repo are never shared with another
-# pipeline. This is not memory efficient (we keep a lot of data), hence this
-# reminder.
-# Note: This job can start immediately.
-setup_baseline:
-  tags:
-    - shell
-    - quartz
-  stage: setup_baseline
+quartz-baseline:
+  stage: sub-pipelines
   variables:
-    GIT_STRATEGY: none
-  script:
-    - mkdir -p ${BUILD_ROOT} && cd ${BUILD_ROOT}
-    - if [ ! -d "tpls" ]; then git clone ${TPLS_REPO}; fi
-    - if [ ! -d "tests" ]; then git clone ${TESTS_REPO}; fi
-    - cd tpls && git pull && cd ..
-    - cd tests && git pull && cd ..
-    - cd ${AUTOTEST_ROOT}
-    - if [ ! -d "autotest" ]; then git clone ${AUTOTEST_REPO}; fi
-    - cd autotest && git pull && cd ..
-  needs: []
+    _AUTOTEST: $AUTOTEST
+  trigger:
+    include: .gitlab/quartz-baseline.yml
+    strategy: depend
 
-.build_toss_3_x86_64_ib_script:
-  script:
-    - export THREADS=12
-    - 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
-
-.build_toss_3_x86_64_ib_corona_script:
-  script:
-    - srun -p mi60 -t 15 -N 1 tests/gitlab/build_and_test
-
-# 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
-# speed-up the allocation.  However this would not be scalable to multiple
-# builds.
-.build_blueos_3_ppc64le_ib_script:
-  script:
-    - lalloc 1 -W 30 -q pdebug tests/gitlab/build_and_test
-
-# Shared script for baseline and sample-run-baseline, the value of BASELINE_TEST
-# differentiates between the two tests.
-.baseline_script: &baseline_script |
-  # locals
-  _glob_err=${BASELINE_TEST}.err
-  _base_diff=${BASELINE_TEST}-${SYS_TYPE}.diff
-  _base_patch=${BASELINE_TEST}-${SYS_TYPE}.patch
-  _base_out=${BASELINE_TEST}-${SYS_TYPE}.out
-  # prepare
-  cd ${BUILD_ROOT}
-  ln -snf ${CI_PROJECT_DIR} mfem
-  cd tests
-  mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
-  # run
-  srun --nodes=1 -p pdebug ../runtest ../../mfem "${BASELINE_TEST} ${ADDITIONAL_DIR}"
-  # post
-  mkdir ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}
-  if [[ -s ${_glob_err} ]]
-  then
-    echo "ERROR during ${BASELINE_TEST} execution";
-    echo "Here is the ${_glob_err} file content";
-    cat ${_glob_err}
-    cp ${_glob_err} ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/${_glob_err}
-    exit 1;
-  elif [[ ! -f ${_base_patch} && ! -f ${_base_out} ]]
-  then
-    echo "Something went WRONG in ${BASELINE_TEST}:";
-    echo "Either ${_base_patch} or ${_base_out} should exists";
-    exit 1;
-  elif [[ -f ${_base_patch} ]]
-  then
-    echo "${BASELINE_TEST}: Differences found, patch generated"
-    cp ${_base_patch} ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/${_base_patch}
-  elif [[ -f ${_base_out} ]]
-  then
-    echo "${BASELINE_TEST}: Differences found, replacement file generated"
-    cp ${_base_out} ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/${_base_out}
-  fi
-  # _base_diff won't even exist if there is no difference.
-  if [[ -f ${_base_diff} ]]
-  then
-    echo "${BASELINE_TEST}: Relevant differences (filtered diff) ..."
-    cat ${_base_diff}
-    cp ${_base_diff} ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/${_base_diff}
-    # We create a .err file, because that's how we signal that there was a diff.
-    cp ${_base_diff} ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/gitlab-${BASELINE_TEST}-${SYS_TYPE}.err
-  fi
-  if [[ ! -s ${_base_diff} ]]
-  then
-    echo "${BASELINE_TEST}: PASSED"
-    true
-  else
-    echo "${BASELINE_TEST}: FAILED"
-    false
-  fi
-
-# Actual templates for baseline checks
-.baselinecheck_mfem:
-  stage: baseline_check
+lassen-build-and-test:
+  stage: sub-pipelines
   variables:
-    BASELINE_TEST: baseline
-    ADDITIONAL_DIR: ${BUILD_ROOT}/tpls
-  script:
-    - *baseline_script
-  artifacts:
-    when: always
-    paths:
-      - ${ARTIFACTS_DIR}
-  allow_failure: true
+    _AUTOTEST: $AUTOTEST
+  trigger:
+    include: .gitlab/lassen-build-and-test.yml
+    strategy: depend
 
-.samplebaselinecheck_mfem:
-  stage: baseline_check
+corona-build-and-test:
+  stage: sub-pipelines
   variables:
-    BASELINE_TEST: sample-runs-baseline
-    ADDITIONAL_DIR: ""
-  script:
-    - *baseline_script
-  timeout: 4h
-  artifacts:
-    when: always
-    paths:
-      - ${ARTIFACTS_DIR}
-  allow_failure: true
-
-# This job can only be manually triggered on a pipeline for master branch, or if
-# the pipeline was triggered with REBASELINE="YES"
-.rebaseline_mfem:
-  stage: baseline_publish
-  rules:
-    - if: '$CI_COMMIT_BRANCH == "master" || $REBASELINE == "YES"'
-      when: manual
-  script:
-    - export PATCH_FILE=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/baseline-${SYS_TYPE}.patch
-    - export FULL_FILE=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/baseline-${SYS_TYPE}.out
-    - export DIFF_FILE=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/baseline-${SYS_TYPE}.diff
-    - cd ${BUILD_ROOT}/tests
-    - |
-      if [[ ! -f "${DIFF_FILE}" ]]
-      then
-          echo "Nothing to be done: no relevant change in baseline"
-          exit 0
-      elif [[ -f "${PATCH_FILE}" ]]
-      then
-          patch "./baseline-${SYS_TYPE}.saved" < "${PATCH_FILE}"
-      elif [[ -f "${FULL_FILE}t" ]]
-      then
-          cp "${FULL_FILE}" "./baseline-${SYS_TYPE}.saved"
-      else
-          echo "File missing: expected ${PATCH_FILE} or ${FULL_FILE}"
-          exit 1
-      fi
-    - git add baseline-${SYS_TYPE}.saved
-    - git commit -m "${SYS_TYPE} rebaselined in GitLab pipeline ${CI_PIPELINE_ID}"
-    - git push origin master
-
-# The list on jobs is defined in machine-specific files.
-include:
-  - local: .gitlab/quartz.yml
-  - local: .gitlab/lassen.yml
+    _AUTOTEST: $AUTOTEST
+  trigger:
+    include: .gitlab/corona-build-and-test.yml
+    strategy: depend

.gitlab/README.md

@@ -0,0 +1,94 @@
+                    Finite Element Discretization Library
+                                   __
+                       _ __ ___   / _|  ___  _ __ ___
+                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
+                      | | | | | ||  _||  __/| | | | | |
+                      |_| |_| |_||_|   \___||_| |_| |_|
+
+                               https://mfem.org
+
+
+This directory contains most of the GitLab CI configuration. MFEM runs both PR
+and nightly testing on GitLab.
+
+# Structure
+
+## Top level
+
+The root configuration file is `.gitlab-ci.yml` at the root of MFEM repo.
+This file only defines one stage, in which we trigger several
+sub-pipelines.
+
+We use sub-pipelines to isolate the test for one combination of `machine`
+and `test type`.
+
+Machines typically include:
+
+* Quartz: Intel bi-socket x86
+* Lassen: Power9 + Nvidia GPU
+* Corona: AMD GPU
+
+Test types include:
+
+* Build and test: Spack driven build of dependencies, mfem build, mfem
+  test
+* Baseline: Script driven build of dependencies, thorough testing
+
+⚠️ The sub-pipeline design allows to add a new machine or a new test type without
+altering the scheduling, execution and displaying of the others.
+
+## Sub-pipelines
+
+Each file is this directory is the root configuration file for one
+sub-pipeline. The naming reflects the corresponding couple (`machine`,
+`test_type`).
+
+Those files define the *stages* and the *jobs* for the sub-pipeline. They
+also contain any configuration that cannot be shared. For the most part
+though, the configuration is shared and is placed in `.gitlab/configs`.
+
+We try to keep scripts out of the CI config and share them among similar
+jobs. They are gathered in `.gitlab/scripts`.
+
+## Scripts
+
+Scripts specific to the CI only are in `.gitlab/scripts`. It is best practice
+to keep scripts outside the CI configuration (no bash scripts embedded in a
+yaml file) because it helps with readability, maintenance and also with
+transition to another CI system.
+
+⚠️ Most of the scripts there are driven by environment variables and do not have a
+usage function. This should be improved.
+
+
+# More testing
+
+## Adding a new target to a build_and_test pipeline
+
+`build_and_test` pipelines rely on Spack to install dependencies. Spack is
+driven by Uberenv which helps freezing Spack configuration: the goal being to
+point to specific commit in Spack and isolate its configuration so that it is
+not influenced by the user environment. More documentation about this can be
+found in `tests/gitlab`.
+
+In the end, the MFEM target for which to build the dependencies is expressed
+with a spack spec of MFEM, within the limits permitted by the MFEM spack
+package.
+
+In any build-and-test sub-pipeline a job basically consists in defining the
+spack spec to use. Adding a job on quartz for example resumes to:
+
+```yaml
+<job_name>:
+  variables:
+    SPEC: "<spack_spec>"
+  extends: .build_and_test_on_quartz
+```
+
+The remaining and non trivial work is to make sure this spec is working. To
+test a spec before adding it, or reproduce a CI configuration, please refer to
+`tests/gitlab/reproduce-ci-jobs-interactively.md`.
+
+⚠️ It is assumed that the spack spec applies to `mfem@develop`. That's why in the
+CI all the specs start with the compiler or the variants to apply to mfem. The
+mechanism still works with a full spec.

.gitlab/configs/common.yml

@@ -0,0 +1,44 @@
+# 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.
+
+# We define the following GitLab pipeline variables:
+variables:
+
+# 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
+# also means that the BUILD_ROOT directory will never be cleaned.
+# TODO: add a clean-up mechanism
+  BUILD_ROOT: ${CI_BUILDS_DIR}/MFEM_${MACHINE_NAME}/${CI_PROJECT_NAME}_${CI_COMMIT_REF_SLUG}_${CI_PIPELINE_ID}
+
+# On LLNL's quartz, there is only one allocation shared among jobs in order to
+# save time and resource. This allocation has to be uniquely named so that we
+# are sure to retrieve it.
+  ALLOC_NAME: ${CI_PROJECT_NAME}_ci_${CI_PIPELINE_ID}
+
+# Defines the default choice for updating the saved baseline results. By default
+# the baseline can only be updated from the master branch. This variable offers
+# the option to manually ask for rebaselining from another branch if necessary.
+  _REBASELINE: "NO"
+  _AUTOTEST: "NO"
+
+# Git repositories used in the pipeline
+  TPLS_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/tpls.git
+  TESTS_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/tests.git
+  AUTOTEST_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/autotest.git
+  MFEM_DATA_REPO: https://github.com/mfem/data.git
+
+# Directory used to place artifacts.
+  ARTIFACTS_DIR: artifacts
+  SLURM_OVERLAP: 1
+
+

.gitlab/configs/corona-config.yml

@@ -0,0 +1,54 @@
+# 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
+variables:
+  MACHINE_NAME: corona
+
+.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 =~ /report/ && $_AUTOTEST != "YES"'
+      when: never
+    # Report success on success status
+    - if: '$CI_JOB_NAME =~ /report_job_success/ && $_AUTOTEST == "YES"'
+      when: on_success
+    # Report failure on failure status
+    - if: '$CI_JOB_NAME =~ /report_job_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
+
+# Spack helped builds
+# Generic corona build job, extending build script
+.build_and_test_on_corona:
+  extends: [.on_corona]
+  stage: build_and_test
+  script:
+    - export THREADS=12
+    - echo ${ALLOC_NAME}
+    - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
+    - echo ${JOBID}
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -t 15 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --build-root "${BUILD_ROOT}" --data
+

.gitlab/configs/lassen-config.yml

@@ -0,0 +1,43 @@
+# 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 pipelines configurations for the Lassen machine at LLNL
+variables:
+  MACHINE_NAME: lassen
+
+.on_lassen:
+  tags:
+    - shell
+    - lassen
+  rules:
+    - if: '$CI_COMMIT_BRANCH =~ /_lnone/ || $ON_LASSEN == "OFF"' #run except if ...
+      when: never
+    # Don't run autotest update if...
+    - if: '$CI_JOB_NAME =~ /report/ && $_AUTOTEST != "YES"'
+      when: never
+    # Report success on success status
+    - if: '$CI_JOB_NAME =~ /report_job_success/ && $_AUTOTEST == "YES"'
+      when: on_success
+    # Report failure on failure status
+    - if: '$CI_JOB_NAME =~ /report_job_failure/ && $_AUTOTEST == "YES"'
+      when: on_failure
+    - when: on_success
+
+# 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 speed-up the allocation.  However this would not be scalable to
+# multiple builds.
+.build_and_test_on_lassen:
+  extends: [.on_lassen]
+  stage: build_and_test
+  script:
+    - lalloc 1 -W 30 -q pdebug tests/gitlab/build_and_test --spec "${SPEC}" --build-root "${BUILD_ROOT}" --data
+  needs: [setup]

.gitlab/configs/quartz-config.yml

@@ -0,0 +1,50 @@
+# 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 pipelines configurations for the Quartz machine at LLNL
+variables:
+  MACHINE_NAME: quartz
+
+.on_quartz:
+  tags:
+    - shell
+    - quartz
+  rules:
+    # Don't run quartz jobs if...
+    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_QUARTZ == "OFF"'
+      when: never
+    # Don't run autotest update if...
+    - if: '$CI_JOB_NAME =~ /report/ && $_AUTOTEST != "YES"'
+      when: never
+    # Report success on success status
+    - if: '$CI_JOB_NAME =~ /report_job_success/ && $_AUTOTEST == "YES"'
+      when: on_success
+    # Report failure on failure status
+    - if: '$CI_JOB_NAME =~ /report_job_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
+
+# Spack helped builds
+# Generic quartz build job, extending build script
+.build_and_test_on_quartz:
+  extends: [.on_quartz]
+  stage: build_and_test
+  script:
+    - export THREADS=12
+    - 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 --spec "${SPEC}" --build-root "${BUILD_ROOT}" --data
+

.gitlab/configs/setup-baseline.yml

@@ -0,0 +1,42 @@
+# 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.
+
+# TPLS_DIR is used in .gitlab/scripts/baseline to provide the tpls location
+# when call the runtest script in MFEM test repo.
+# Note: the value must be consistent with what setup_baseline does.
+variables:
+  TPLS_DIR: ${BUILD_ROOT}/tpls
+  AUTOTEST_ROOT: ${CI_BUILDS_DIR}/MFEM_${MACHINE_NAME}_baseline
+
+# The setup_baseline job doesn't rely on MFEM git repo. It prepares a
+# pipeline-wide working directory downloading/updating external repos.
+# TODO:
+# updating tests and tpls is not necessary anymore since pipelines are now
+# using unique directories so repo are never shared with another pipeline. This
+# is not memory efficient (we keep a lot of data), hence this reminder.
+setup_baseline:
+  tags:
+    - shell
+    - quartz
+  stage: setup
+  variables:
+    GIT_STRATEGY: none
+  script:
+    - echo "BUILD_ROOT ${BUILD_ROOT}"
+    - mkdir -p ${BUILD_ROOT} && cd ${BUILD_ROOT}
+    - if [ ! -d "tpls" ]; then git clone ${TPLS_REPO}; fi
+    - if [ ! -d "tests" ]; then git clone ${TESTS_REPO}; fi
+    - cd tpls && git pull && cd ..
+    - cd tests && git pull origin && cd ..
+    - echo "AUTOTEST_ROOT ${AUTOTEST_ROOT}"
+    - mkdir -p ${AUTOTEST_ROOT} && cd ${AUTOTEST_ROOT}
+    - if [ ! -d "autotest" ]; then git clone ${AUTOTEST_REPO}; fi
+    - cd autotest && git pull && cd ..

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

@@ -0,0 +1,34 @@
+# 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.
+
+variables:
+  AUTOTEST_ROOT: ${CI_BUILDS_DIR}/MFEM_${MACHINE_NAME}_build_and_test
+
+# setup clones the mfem/data repo in ${BUILD_ROOT}. The build_and_test script
+# then symlinks the repo to the parent directory of the MFEM source directory.
+# Unit tests that depend on the mfem/data repo will then detect that this
+# directory is present and be enabled.
+setup:
+  tags:
+    - shell
+    - quartz
+  stage: setup
+  variables:
+    GIT_STRATEGY: none
+  script:
+    - echo "BUILD_ROOT ${BUILD_ROOT}"
+    - mkdir -p ${BUILD_ROOT} && cd ${BUILD_ROOT}
+    - if [ ! -d data ]; then git clone ${MFEM_DATA_REPO}; fi
+    - echo "AUTOTEST_ROOT ${AUTOTEST_ROOT}"
+    - mkdir -p ${AUTOTEST_ROOT} && cd ${AUTOTEST_ROOT}
+    - if [ ! -d "autotest" ]; then git clone ${AUTOTEST_REPO}; fi
+    - cd autotest && git pull && cd ..
+

.gitlab/corona-build-and-test.yml

@@ -0,0 +1,63 @@
+# 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.
+
+stages:
+  - setup
+  - allocate_resource
+  - build_and_test
+  - release_resource_and_report
+
+# Slurm shared allocation
+allocate_resource:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_corona
+  stage: allocate_resource
+  script:
+    - salloc --exclusive --nodes=1 --partition=mi60 --time=30 --no-shell --job-name=${ALLOC_NAME}
+  timeout: 6h
+  needs: [setup]
+
+# Build and test jobs, simply provide a spec
+rocm_gcc_8.3.1:
+  variables:
+    SPEC: "@develop%gcc@8.3.1+rocm amdgpu_target=gfx906"
+  extends: .build_and_test_on_corona
+  needs: [allocate_resource]
+
+# Release slurm allocation
+release_resource:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_corona
+  stage: release_resource_and_report
+  script:
+    - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
+    - ([[ -n "${JOBID}" ]] && scancel ${JOBID})
+  needs: [rocm_gcc_8.3.1]
+
+# Jobs report
+report_job_success:
+  extends: .on_corona
+  stage: release_resource_and_report
+  script:
+    - .gitlab/scripts/report_build_and_test_success
+
+report_job_failure:
+  extends: .on_corona
+  stage: release_resource_and_report
+  script:
+    - .gitlab/scripts/report_build_and_test_failure
+
+include:
+  - local: .gitlab/configs/common.yml
+  - local: .gitlab/configs/corona-config.yml
+  - local: .gitlab/configs/setup-build-and-test.yml

.gitlab/lassen-build-and-test.yml

@@ -9,26 +9,30 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-# GitLab pipelines configurations for the Lassen machine at LLNL
-
-.on_lassen:
-  tags:
-    - shell
-    - lassen
-  rules:
-    - if: '$CI_COMMIT_BRANCH =~ /_lnone/ || $ON_LASSEN == "OFF"' #run except if ...
-      when: never
-    - when: on_success
-
-# Spack helped builds
-# Generic lassen build job, extending build script
-# Note: Lassen jobs can start as soon as the setup job is complete.
-.build_and_test_on_lassen:
-  extends: [.build_blueos_3_ppc64le_ib_script, .on_lassen]
-  stage: l_build_and_test
-  needs: [setup]
+stages:
+  - setup
+  - build_and_test
+  - report
 
 opt_mpi_cuda_xl_16_1_1_8:
   variables:
-    SPEC: "%xl@16.1.1.8 +mpi +cuda cuda_arch=sm_70"
+    SPEC: "%xl@16.1.1.8 +mpi +cuda cuda_arch=70"
   extends: .build_and_test_on_lassen
+
+# Jobs report
+report_job_success:
+  extends: .on_lassen
+  stage: report
+  script:
+    - .gitlab/scripts/report_build_and_test_success
+
+report_job_failure:
+  extends: .on_lassen
+  stage: report
+  script:
+    - .gitlab/scripts/report_build_and_test_failure
+
+include:
+  - local: .gitlab/configs/common.yml
+  - local: .gitlab/configs/lassen-config.yml
+  - local: .gitlab/configs/setup-build-and-test.yml

.gitlab/quartz-baseline.yml

@@ -0,0 +1,64 @@
+# 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.
+
+variables:
+  BASELINE_TEST: baseline
+
+stages:
+  - setup
+  - baseline_check
+  - baseline_report
+  - baseline_publish
+
+baselinecheck_mfem_intel_quartz:
+  extends: [.on_quartz]
+  stage: baseline_check
+  script:
+    - .gitlab/scripts/baseline
+  artifacts:
+    when: always
+    paths:
+      - ${ARTIFACTS_DIR}
+  allow_failure: true
+
+report_baseline:
+  extends: [.on_quartz]
+  stage: baseline_report
+  script:
+    - cd ${AUTOTEST_ROOT}/autotest && git pull
+    - mkdir -p ${MACHINE_NAME}
+    - rundir="${MACHINE_NAME}/$(date +%Y-%m-%d)-gitlab-${BASELINE_TEST}-${CI_COMMIT_REF_SLUG}"
+    - rundir=$(${CI_PROJECT_DIR}/.gitlab/scripts/safe_create_rundir ${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).
+    - |
+      if [[ -f ${rundir}/*.err ]]
+      then
+        echo "See the pipeline here -> $CI_PIPELINE_URL" >> ${rundir}/*.err
+        cp ${rundir}/*.err ${rundir}/autotest-email.html
+      fi
+    - git add ${rundir}
+    - git commit -am "GitLab CI log for ${BASELINE_TEST} on ${MACHINE_NAME} with intel ($(date +%Y-%m-%d))"
+    - git push origin master
+
+baselinepublish_mfem_quartz:
+  extends: [.on_quartz]
+  stage: baseline_publish
+  rules:
+    - if: '$CI_COMMIT_BRANCH == "master" || $REBASELINE == "YES"'
+      when: manual
+  script:
+    - .gitlab/scripts/rebaseline
+
+include:
+  - local: .gitlab/configs/common.yml
+  - local: .gitlab/configs/quartz-config.yml
+  - local: .gitlab/configs/setup-baseline.yml

.gitlab/quartz-build-and-test.yml

@@ -0,0 +1,95 @@
+# 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.
+
+stages:
+  - setup
+  - allocate_resource
+  - build_and_test
+  - release_resource_and_report
+
+# Allocate
+allocate_resource:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_quartz
+  stage: allocate_resource
+  script:
+    - salloc --exclusive --nodes=1 --partition=pdebug --time=30 --no-shell --job-name=${ALLOC_NAME}
+  timeout: 6h
+
+# GitLab jobs for the Quartz machine at LLNL
+debug_ser_gcc_4_9_3:
+  variables:
+    SPEC: "%gcc@4.9.3 +debug~mpi"
+  extends: .build_and_test_on_quartz
+
+debug_ser_gcc_6_1_0:
+  variables:
+    SPEC: "%gcc@6.1.0 +debug~mpi"
+  extends: .build_and_test_on_quartz
+
+debug_par_gcc_6_1_0:
+  variables:
+    SPEC: "%gcc@6.1.0 +debug+mpi"
+  extends: .build_and_test_on_quartz
+
+opt_ser_gcc_6_1_0:
+  variables:
+    SPEC: "%gcc@6.1.0 ~mpi"
+  extends: .build_and_test_on_quartz
+
+opt_par_gcc_6_1_0:
+  variables:
+    SPEC: "%gcc@6.1.0"
+  extends: .build_and_test_on_quartz
+
+opt_par_gcc_6_1_0_sundials:
+  variables:
+    SPEC: "%gcc@6.1.0 +sundials"
+  extends: .build_and_test_on_quartz
+
+opt_par_gcc_6_1_0_petsc:
+  variables:
+    SPEC: "%gcc@6.1.0 +petsc ^petsc+mumps"
+  extends: .build_and_test_on_quartz
+
+opt_par_gcc_6_1_0_pumi:
+  variables:
+    SPEC: "%gcc@6.1.0 +pumi"
+  extends: .build_and_test_on_quartz
+
+# Release
+release_resource:
+  variables:
+    GIT_STRATEGY: none
+  extends: .on_quartz
+  stage: release_resource_and_report
+  script:
+    - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
+    - ([[ -n "${JOBID}" ]] && scancel ${JOBID})
+
+# Jobs report
+report_job_success:
+  extends: .on_quartz
+  stage: release_resource_and_report
+  script:
+    - .gitlab/scripts/report_build_and_test_success
+
+report_job_failure:
+  extends: .on_quartz
+  stage: release_resource_and_report
+  script:
+    - .gitlab/scripts/report_build_and_test_failure
+
+include:
+  - local: .gitlab/configs/common.yml
+  - local: .gitlab/configs/quartz-config.yml
+  - local: .gitlab/configs/setup-build-and-test.yml

.gitlab/quartz.yml

@@ -1,184 +0,0 @@
-# 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 pipelines configurations for the Quartz machine at LLNL
-
-.on_quartz:
-  tags:
-    - shell
-    - quartz
-  rules:
-    # Don't run quartz jobs if...
-    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_QUARTZ == "OFF"'
-      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 =~ /q_report/ && $AUTOTEST != "YES"'
-      when: never
-    # Report success on success status
-    - if: '$CI_JOB_NAME =~ /q_report_success/ && $AUTOTEST == "YES"'
-      when: on_success
-    # 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/'
-      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
-q_allocate_resources:
-  variables:
-    GIT_STRATEGY: none
-  extends: .on_quartz
-  stage: q_allocate_resources
-  script:
-    - salloc --exclusive --nodes=1 --partition=pdebug --time=30 --no-shell --job-name=${ALLOC_NAME}
-  timeout: 6h
-
-# Release
-q_release_resources:
-  variables:
-    GIT_STRATEGY: none
-  extends: .on_quartz
-  stage: q_release_resources
-  script:
-    - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
-    - ([[ -n "${JOBID}" ]] && scancel ${JOBID})
-
-# Release
-q_report_success:
-  variables:
-    GIT_STRATEGY: none
-  extends: .on_quartz
-  stage: q_release_resources
-  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}"
-    - *safe_create_rundir
-    - echo "The Quartz 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 quartz with intel ($(date +%Y-%m-%d))"
-    - git push origin master
-
-q_report_failure:
-  variables:
-    GIT_STRATEGY: none
-  extends: .on_quartz
-  stage: q_release_resources
-  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}"
-    - *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
-    - cp ${rundir}/gitlab.err ${rundir}/autotest-email.html
-    - git add ${rundir}
-    - git commit -am "Gitlab CI log for baseline on quartz with intel ($(date +%Y-%m-%d))"
-    - git push origin master
-
-# Spack helped builds
-# Generic quartz build job, extending build script
-.build_and_test_on_quartz:
-  extends: [.build_toss_3_x86_64_ib_script, .on_quartz]
-  stage: q_build_and_test
-
-# Build MFEM
-debug_ser_gcc_4_9_3:
-  variables:
-    SPEC: "%gcc@4.9.3 +debug~mpi"
-  extends: .build_and_test_on_quartz
-
-debug_ser_gcc_6_1_0:
-  variables:
-    SPEC: "%gcc@6.1.0 +debug~mpi"
-  extends: .build_and_test_on_quartz
-
-debug_par_gcc_6_1_0:
-  variables:
-    SPEC: "%gcc@6.1.0 +debug+mpi"
-  extends: .build_and_test_on_quartz
-
-opt_ser_gcc_6_1_0:
-  variables:
-    SPEC: "%gcc@6.1.0 ~mpi"
-  extends: .build_and_test_on_quartz
-
-opt_par_gcc_6_1_0:
-  variables:
-    SPEC: "%gcc@6.1.0"
-  extends: .build_and_test_on_quartz
-
-opt_par_gcc_6_1_0_sundials:
-  variables:
-    SPEC: "%gcc@6.1.0 +sundials"
-  extends: .build_and_test_on_quartz
-
-opt_par_gcc_6_1_0_petsc:
-  variables:
-    SPEC: "%gcc@6.1.0 +petsc ^petsc+mumps"
-  extends: .build_and_test_on_quartz
-
-opt_par_gcc_6_1_0_pumi:
-  variables:
-    SPEC: "%gcc@6.1.0 +pumi"
-  extends: .build_and_test_on_quartz
-
-# Baseline jobs form an independent set of jobs. We use `needs:[]` to specify
-# that "setup-baseline" can start immediately. Then, we have to use needs for
-# each one of the baseline jobs, otherwise they will wait for the rest of the
-# pipeline.
-
-baselinecheck_mfem_intel_quartz:
-  extends: [.baselinecheck_mfem, .on_quartz]
-  needs: [setup_baseline]
-
-update_autotest:
-  extends: [.on_quartz]
-  needs: [baselinecheck_mfem_intel_quartz]
-  stage: baseline_to_autotest
-  script:
-    - cd ${AUTOTEST_ROOT}/autotest && git pull
-    - rundir="quartz/$(date +%Y-%m-%d)-github-${CI_COMMIT_REF_SLUG}"
-    - *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).
-    - |
-      if [[ -f ${rundir}/*.err ]]
-      then
-        echo "See the pipeline here -> $CI_PIPELINE_URL" >> ${rundir}/*.err
-        cp ${rundir}/*.err ${rundir}/autotest-email.html
-      fi
-    - git add ${rundir}
-    - git commit -am "Gitlab CI log for baseline on quartz with intel ($(date +%Y-%m-%d))"
-    - git push origin master
-
-baselinepublish_mfem_quartz:
-  extends: [.on_quartz, .rebaseline_mfem]
-  needs: [baselinecheck_mfem_intel_quartz]

.gitlab/scripts/baseline

@@ -0,0 +1,80 @@
+#!/bin/bash
+
+# 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.
+
+# locals
+glob_err=${BASELINE_TEST}.err
+base=${BASELINE_TEST}-${SYS_TYPE}
+base_diff=${base}.diff
+base_patch=${base}.patch
+base_out=${base}.out
+artifacts_path=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}
+
+# prepare
+cd ${BUILD_ROOT}
+ln -snf ${CI_PROJECT_DIR} mfem
+cd tests
+[[ -d _${BASELINE_TEST} ]] && rm -rf _${BASELINE_TEST}
+mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
+
+# run
+if [[ "${MACHINE_NAME}" == "quartz" || "${MACHINE_NAME}" == "ruby" ]]; then
+  srun --nodes=1 -p pdebug ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+elif [[ ${MACHINE_NAME} == "corona" ]]; then
+  srun --nodes=1 -t 60 -p mi60 ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+elif [[ ${MACHINE_NAME} == "lassen" ]]; then
+  lalloc 1 -q pdebug ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+fi
+
+# post
+mkdir ${artifacts_path}
+
+if [[ -s ${glob_err} ]]
+then
+  echo "ERROR during ${BASELINE_TEST} execution";
+  echo "Here is the ${glob_err} file content";
+  cat ${glob_err}
+  cp ${glob_err} ${artifacts_path}/${glob_err}
+  exit 1;
+elif [[ ! -f ${base_patch} && ! -f ${base_out} ]]
+then
+  echo "Something went WRONG in ${BASELINE_TEST}:";
+  echo "Either ${base_patch} or ${base_out} should exists";
+  exit 1;
+elif [[ -f ${base_patch} ]]
+then
+  echo "${BASELINE_TEST}: Differences found, patch generated"
+  cp ${base_patch} ${artifacts_path}/${base_patch}
+elif [[ -f ${base_out} ]]
+then
+  echo "${BASELINE_TEST}: Differences found, replacement file generated"
+  cp ${base_out} ${artifacts_path}/${base_out}
+fi
+
+# base_diff won't even exist if there is no difference.
+if [[ -f ${base_diff} ]]
+then
+  echo "${BASELINE_TEST}: Relevant differences (filtered diff) ..."
+  cat ${base_diff}
+  cp ${base_diff} ${artifacts_path}/${base_diff}
+  # We create a .err file, because that's how we signal that there was a diff.
+  cp ${base_diff} ${artifacts_path}/gitlab-${BASELINE_TEST}-${MACHINE_NAME}.err
+fi
+
+if [[ ! -s ${base_diff} ]]
+then
+  echo "${BASELINE_TEST}: PASSED"
+  true
+else
+  echo "${BASELINE_TEST}: FAILED"
+  false
+fi

.gitlab/scripts/rebaseline

@@ -0,0 +1,49 @@
+#!/bin/bash
+
+# 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.
+
+# There will be collision between corona and quartz baselines.
+# Once the corresponding files have been generated, we can switch to machine
+# specific ref.
+ARTIFACT_PATH=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/baseline-${SYS_TYPE}
+#ARTIFACT_PATH=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/baseline-${SYS_TYPE}-${MACHINE_NAME}
+
+PATCH_FILE=${ARTIFACT_PATH}.patch
+FULL_FILE=${ARTIFACT_PATH}.out
+DIFF_FILE=${ARTIFACT_PATH}.diff
+
+# There will be collision between corona and quartz baselines.
+# Once the corresponding files have been generated, we can switch to machine
+# specific ref.
+SAVED_NAME=baseline-${SYS_TYPE}.saved
+#SAVED_NAME=baseline-${SYS_TYPE}-${MACHINE_NAME}.saved
+
+cd ${BUILD_ROOT}/tests
+
+if [[ ! -f "${DIFF_FILE}" ]]
+then
+    echo "Nothing to be done: no relevant change in baseline"
+    exit 0
+elif [[ -f "${PATCH_FILE}" ]]
+then
+    patch "${SAVED_NAME}" < "${PATCH_FILE}"
+elif [[ -f "${FULL_FILE}" ]]
+then
+    cp "${FULL_FILE}" "${SAVED_NAME}"
+else
+    echo "File missing: expected ${PATCH_FILE} or ${FULL_FILE}"
+    exit 1
+fi
+
+git add "${SAVED_NAME}"
+git commit -m "${SYS_TYPE} (${MACHINE_NAME}) rebaselined in GitLab pipeline ${CI_PIPELINE_ID}"
+git push origin master

.gitlab/scripts/report_build_and_test_failure

@@ -0,0 +1,32 @@
+#!/bin/bash
+
+# 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.
+
+echo "Runs if there was at least one failure on ${MACHINE_NAME}"
+
+cd ${AUTOTEST_ROOT}/autotest && git pull
+mkdir -p ${MACHINE_NAME}
+
+rundir="${MACHINE_NAME}/$(date +%Y-%m-%d)-gitlab-ci-${CI_COMMIT_REF_SLUG}"
+rundir=$(${CI_PROJECT_DIR}/.gitlab/scripts/safe_create_rundir $rundir)
+
+echo "There was an error while running CI on ${MACHINE_NAME}" > ${rundir}/gitlab.err
+echo "See the pipeline here -> $CI_PIPELINE_URL" >> ${rundir}/gitlab.err
+
+msg="GitLab CI log for build-and-test on ${MACHINE_NAME} ($(date +%Y-%m-%d))"
+
+cp ${rundir}/gitlab.err ${rundir}/autotest-email.html
+
+git pull
+git add ${rundir}
+git commit -am "${msg}"
+git push origin master

.gitlab/scripts/report_build_and_test_success

@@ -0,0 +1,30 @@
+#!/bin/bash
+
+# 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.
+
+echo "Can only run if all the ${MACHINE_NAME} jobs passed"
+
+cd ${AUTOTEST_ROOT}/autotest && git pull
+mkdir -p ${MACHINE_NAME}
+
+rundir="${MACHINE_NAME}/$(date +%Y-%m-%d)-gitlab-ci-${CI_COMMIT_REF_SLUG}"
+rundir=$(${CI_PROJECT_DIR}/.gitlab/scripts/safe_create_rundir $rundir)
+
+echo "The ${MACHINE_NAME} jobs were successful" > ${rundir}/gitlab.out
+echo "See the pipeline here -> $CI_PIPELINE_URL" >> ${rundir}/gitlab.err
+
+msg="GitLab CI log for build-and-test on ${MACHINE_NAME} ($(date +%Y-%m-%d))"
+
+git pull
+git add ${rundir}
+git commit -am "${msg}"
+git push origin master

.gitlab/scripts/safe_create_rundir

@@ -0,0 +1,42 @@
+#!/bin/bash
+
+# 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.
+
+# This script takes a seed for a directory name and appends it with a counter
+# incremented until it can create a new directory with it.
+
+# Usage:
+#
+# Expects 1 argument: a string that is use as a seed for the directory name.
+#
+# > rundir="desired_name"
+# > rundir=$(./safe_create_rundir $rundir)
+
+set -o errexit
+set -o nounset
+
+rundir=${1:-""}
+if [[ -z ${rundir} ]]; then
+  >&2 echo "The script expects a string as argument for directory creation."
+  exit 1
+fi
+
+if ! mkdir ${rundir}; then
+  n=1
+  while ! mkdir ${rundir}_${n}
+  do
+    n=$((n+1))
+  done
+  rundir=${rundir}_${n}
+fi
+
+echo $rundir

CHANGELOG

@@ -10,6 +10,62 @@
 
 Version 4.3.1 (development)
 ===========================
+- Added support for mesh preprocessing to resolve fine scale problem data
+  before simulation. This feature uses adaptive mesh refinement to control the
+  associated data oscillation error. See the new Example 30/30p.
+
+- Switched from Artistic Style (astyle) version 2.05.1 to version 3.1 for code
+  formatting. See the "make style" target.
+
+- Split the fem/fe.?pp files into separate files in the new fem/fe/ directory
+  to simplify and clarify the organization of FiniteElement classes.
+
+- Added support for hr-adaptivity using TMOP-based error estimator.
+
+- Coefficient::SetTime now propagates the new time into internally stored
+  Coefficient objects.
+
+- Added initial support for google-benchmarks in the tests/benchmarks directory.
+  It can be enabled with MFEM_USE_BENCHMARK=YES.
+
+- Added Binder (mybinder.org) configuration files for C++ MFEM Jupyter Notebooks
+  with inline GLVis visualization as well as a new examples/jupyter/ directory
+  with a sample notebook based on Example 1. Implementation based on xeus-cling,
+  github.com/jupyter-xeus/xeus-cling + xeus-glvis, github.com/GLVis/xeus-glvis.
+
+- Added 'double' atomicAdd implementation for previous versions of CUDA.
+
+- Adding lowest order Nedelec and Raviart-Thomas basis functions on wedge
+  shaped elements.
+
+- Added initial support for meshes with pyramidal elements, including several
+  pyramidal meshes in the data/ directory and support for the lowest order H1,
+  Nedelec, Raviart-Thomas, and L2 basis functions on pyramids.
+
+- Updated the hypre interface according to changes in hypre-2.22.1. The ADS
+  solver is now fully working on GPUs.
+
+- Tetrahedral meshes no longer need to be reordered to support high order
+  Nedelec basis functions. This will allow future support for Nedelec basis
+  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.
+
+- Added new miniapps that use the ParELAG library, its hybrid smoothers, and the
+  hierarchy of spaces created by the element-based AMG (AMGe) methodology in
+  ParELAG to build multigrid solvers for H(curl) and H(div) forms. See the
+  miniapps/parelag directory for more details.
+
+- Fixed several MinGW build issues on Windows.
+
+- Remove the 'u' flag in the ar command, to update all files in the archive,
+  avoiding file name collisions from different subdirectories.
+
+- Added initial TMOP-based capabilities for surface fitting and tangential
+  relaxation in the mesh-optimizer and pmesh-optimizer miniapps.
 
 
 Version 4.3, released on July 29, 2021

CITATION.cff

@@ -0,0 +1,66 @@
+cff-version: 1.2.0
+message: "If you use MFEM, please cite it as follows."
+authors:
+- family-names: "MFEM Team"
+title: "MFEM: Modular Finite Element Methods [Software]"
+doi: 10.11578/dc.20171025.1248
+url: "https://mfem.org"
+preferred-citation:
+  type: article
+  authors:
+  - family-names: "Anderson"
+    given-names: "Robert"
+    orcid: "https://orcid.org/0000-0002-3508-9944"
+  - family-names: "Andrej"
+    given-names: "Julian"
+    orcid: "https://orcid.org/0000-0001-7661-4840"
+  - family-names: "Barker"
+    given-names: "Andrew"
+    orcid: "https://orcid.org/0000-0003-3572-911X"
+  - family-names: "Bramwell"
+    given-names: "Jamie"
+  - family-names: "Camier"
+    given-names: "Jean-Sylvain"
+    orcid: "https://orcid.org/0000-0003-2421-1999"
+  - family-names: "Cerveny"
+    given-names: "Jakub"
+    orcid: "https://orcid.org/0000-0003-4231-2531"
+  - family-names: "Dobrev"
+    given-names: "Veselin"
+    orcid: "https://orcid.org/0000-0003-1793-5622"
+  - family-names: "Dudouit"
+    given-names: "Yohann"
+    orcid: "https://orcid.org/0000-0001-5831-561X"
+  - family-names: "Fisher"
+    given-names: "Aaron"
+  - family-names: "Kolev"
+    given-names: "Tzanio"
+    orcid: "https://orcid.org/0000-0002-2810-3090"
+  - family-names: "Pazner"
+    given-names: "Will"
+    orcid: "https://orcid.org/0000-0003-4885-2934"
+  - family-names: "Stowell"
+    given-names: "Mark"
+    orcid: "https://orcid.org/0000-0002-5389-7435"
+  - family-names: "Tomov"
+    given-names: "Vladimir"
+    orcid: "https://orcid.org/0000-0002-1846-6816"
+  - family-names: "Akkerman"
+    given-names: "Ido"
+    orcid: "https://orcid.org/0000-0002-5937-0300"
+  - family-names: "Dahm"
+    given-names: "Johann"
+    orcid: "https://orcid.org/0000-0001-9657-3564"
+  - family-names: "Medina"
+    given-names: "David"
+  - family-names: "Zampini"
+    given-names: "Stefano"
+    orcid: "https://orcid.org/0000-0002-0435-0433"
+  doi: "10.1016/j.camwa.2020.06.009"
+  journal: "Computers \\& Mathematics with Applications"
+  month: 1
+  start: 42 # First page number
+  end: 74 # Last page number
+  title: "MFEM: A Modular Finite Element Methods Library"
+  volume: 81
+  year: 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.")
@@ -328,11 +333,11 @@ if (MFEM_USE_AMGX)
 endif()
 
 if (MFEM_USE_CONDUIT)
-    find_package(Conduit REQUIRED conduit relay blueprint )
+    find_package(Conduit REQUIRED conduit relay blueprint)
 endif()
 
 if (MFEM_USE_FMS)
-    find_package(FMS REQUIRED fms )
+    find_package(FMS REQUIRED fms)
 endif()
 
 # Axom/Sidre
@@ -377,6 +382,11 @@ if (MFEM_USE_UMPIRE)
    find_package(UMPIRE REQUIRED)
 endif()
 
+# GOOGLE-BENCHMARK
+if (MFEM_USE_BENCHMARK)
+   find_package(Benchmark REQUIRED)
+endif()
+
 # Caliper
 if (MFEM_USE_CALIPER)
   find_package(Caliper REQUIRED)
@@ -405,6 +415,11 @@ if (MFEM_USE_MKL_CPARDISO)
   endif()
 endif()
 
+# PARELAG
+if (MFEM_USE_PARELAG)
+   find_package(PARELAG REQUIRED)
+endif()
+
 # MFEM_TIMER_TYPE
 if (NOT DEFINED MFEM_TIMER_TYPE)
   if (APPLE)
@@ -428,10 +443,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 BENCHMARK PARELAG MPI_CXX)
 
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
@@ -496,13 +511,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}")
@@ -541,15 +552,21 @@ endif()
 set(MFEM_CUSTOM_TARGET_PREFIX CACHE STRING "")
 
 #-------------------------------------------------------------------------------
-# Examples, miniapps, and testing
+# Examples, miniapps, benchmarks and testing
 #-------------------------------------------------------------------------------
 
-# Enable testing if required
+# Enable testing and benchmarks if required
 if (MFEM_ENABLE_TESTING)
   enable_testing()
   set(MFEM_ALL_TESTS_TARGET_NAME tests)
   add_mfem_target(${MFEM_ALL_TESTS_TARGET_NAME} OFF)
   add_subdirectory(tests EXCLUDE_FROM_ALL)
+
+  # Create a target for all benchmarks and, optionally, enable it.
+  set(MFEM_ALL_BENCHMARKS_TARGET_NAME benchmarks)
+  add_mfem_target(${MFEM_ALL_BENCHMARKS_TARGET_NAME}
+                  ${MFEM_ENABLE_GOOGLE_BENCHMARKS})
+  add_subdirectory(tests/benchmarks EXCLUDE_FROM_ALL)
 endif()
 
 # Define a target that all examples and miniapps will depend on.
@@ -559,7 +576,11 @@ add_custom_target(${MFEM_EXEC_PREREQUISITES_TARGET_NAME})
 # Create a target for all examples and, optionally, enable it.
 set(MFEM_ALL_EXAMPLES_TARGET_NAME examples)
 add_mfem_target(${MFEM_ALL_EXAMPLES_TARGET_NAME} ${MFEM_ENABLE_EXAMPLES})
-add_subdirectory(examples EXCLUDE_FROM_ALL)
+if (MFEM_ENABLE_EXAMPLES)
+  add_subdirectory(examples) #install examples if enabled
+else()
+  add_subdirectory(examples EXCLUDE_FROM_ALL)
+endif()
 
 # Create a target for all miniapps and, optionally, enable it.
 set(MFEM_ALL_MINIAPPS_TARGET_NAME miniapps)
@@ -569,6 +590,7 @@ add_subdirectory(miniapps EXCLUDE_FROM_ALL)
 # Target to build all executables, i.e. everything.
 add_custom_target(exec)
 add_dependencies(exec
+  ${MFEM_ALL_BENCHMARKS_TARGET_NAME}
   ${MFEM_ALL_EXAMPLES_TARGET_NAME}
   ${MFEM_ALL_MINIAPPS_TARGET_NAME}
   ${MFEM_ALL_TESTS_TARGET_NAME})

CONTRIBUTING.md

@@ -42,6 +42,7 @@ back to them before issuing pull requests:
   - [New Feature Development](#new-feature-development)
   - [Developer Guidelines](#developer-guidelines)
   - [Pull Requests](#pull-requests)
+  - [MFEM PR Rules](#mfem-pr-rules)
   - [Pull Request Checklist](#pull-request-checklist)
   - [Master/Next Workflow](#masternext-workflow)
   - [Releases](#releases)
@@ -67,8 +68,9 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
   with regards to documentation and code styling.
 - Please do not commit large/binary files to the central repository (use a fork
   instead).
-- Pull requests  should be issued toward `mfem:master`. Make sure
-  to check the items off the [Pull Request Checklist](#pull-request-checklist).
+- Pull requests should be issued toward `mfem:master`. Make sure
+  to check the items off the [Pull Request Checklist](#pull-request-checklist) and
+  follow the [MFEM PR Rules](#mfem-pr-rules).
 - When your contribution is fully working and ready to be reviewed, add
   the `ready-for-review` label.
 - PRs are treated similarly to journal submission with an "editor" assigning two
@@ -105,12 +107,14 @@ The MFEM source code has the following structure:
   │   ├── caliper
   │   ├── ginkgo
   │   ├── hiop
+  │   ├── jupyter
   │   ├── petsc
   │   ├── pumi
   │   ├── sundials
   |   └── superlu
   ├── fem
   │   ├── ceed
+  │   ├── fe
   │   ├── qinterp
   │   └── tmop
   ├── general
@@ -126,6 +130,7 @@ The MFEM source code has the following structure:
   │   ├── mtop
   │   ├── navier
   │   ├── nurbs
+  │   ├── parelag
   │   ├── performance
   │   ├── shifted
   │   ├── solvers
@@ -323,15 +328,22 @@ Before you can start, you need a GitHub account, here are a few suggestions:
     change the code by default.
 
 - Code specifics
-  - All significant new classes, methods and functions have Doxygen-style
-    documentation in source comments.
+  - All new public, protected, and private classes, methods, data members, and
+    functions have Doxygen-style documentation in source comments.
+  - In addition to arguments and functionality, documentation should include the
+    current limitations of the code, any background information that is
+    implicitly assumed in the implementation, and the ownership and lifetime
+    of data.
   - 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
     conflicted files in the commit message.
+  - All significant new features and changes should be documented in CHANGELOG.
+  - New examples and miniapps should have documentation on the MFEM webpage.
+
 
 ### Pull Requests
 
@@ -397,6 +409,83 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 - If triggered, track the status of the LLNL GitLab tests. If failing, ask
   one of the _LLNL developers_ for details.
 
+
+### MFEM PR Rules
+
+The Pull Request (PR) approval process in MFEM is similar to the approval of papers in a peer-reviewed journal. In particular:
+
+1. There is an MFEM board of "editors" that evaluates new PRs and assigns "reviewers" for each PR.
+
+2. The assigned reviewers are responsible to carefully review and test the proposed PR.
+
+3. A PR can be (manually) merged in the *next* branch only if 2 of the assigned reviewers have approved it and it has passed internal testing. This merge can be performed by any of the assigned reviewers or by any of the editors.
+
+4. A PR can be merged in the *master* branch only if it has been tested successfully for a week in *next* and an editor has (optionally) taken a final look. This merge can be performed only by one of the editors.
+
+#### Responsibilities of Editors
+
+The current list of MFEM editors is:
+
+- @v-dobrev (Veselin Dobrev)
+- @tzanio (Tzanio Kolev)
+- @pazner (Will Pazner)
+- @mlstowell (Mark Stowell)
+
+**The responsibilities of the editors are:**
+
+1. To assign appropriate milestone and labels for new PRs, e.g. *bugfix*, *minor*, *api-change*, *high-impact*, etc.
+
+2. To assign at least 2 reviewers for new PRs. An editor can also be a reviewer. The editor, reviewers, and author should be listed as "Assignees" on the GitHub PR page. After assignment, the `in-review` label should be added.
+
+3. To complete the initial PR evaluation and assignments in a timely manner: 1 week from submission.
+
+4. To assist reviewers when they need help with their reviews (but also to stay out of the way when they don't).
+
+5. To remind the reviewers about timely completion of their review.
+
+6. To take a final look and complete the PR merge in *master*. The final look step is optional and shouldn't take more than 3 days.
+
+7. The assignment of bugfixes should be expedited proportional to their importance, e.g. in some cases the editor can assign much shorter review window.
+
+#### Responsibilities of Reviewers
+
+Everyone on the MFEM team can be asked to serve as a reviewer on a PR in their area of expertise.
+
+**The responsibilities of the reviewers are:**
+
+1. To let the editors know if the proposed assignment is not a good match for them.
+
+2. To communicate with the PR author, provide feedback and work with them to resolve issues.
+
+3. To ensure the quality of the PR by making sure that the code adheres to the [Developer Guidelines](#developer-guidelines), e.g. all methods, data members, and functions have documentation, including data ownership and lifetime, new examples/miniapps have a corresponding PR in mfem/web, major features have `CHANGELOG` entries, etc.
+
+3. To seek help from the editors in case of difficulties.
+
+4. To complete the review in a timely manner: 3 weeks from assignment.
+
+5. To test the PR thoroughly before merging in *next*. The PR author is also encouraged to perform testing and inform the reviewers about the results.
+
+6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *master*.
+
+7. The review of bugfixes should be expedited proportional to their importance. The review window can be much less than three weeks in such cases.
+
+#### Responsibilities of Authors
+
+Authors should clearly indicate when a PR is ready for review (before that the PR should be marked as `Draft` or `[WIP]`).
+
+**The responsibilities of the authors are:**
+
+1. To follow the instructions and PR checklist in the `CONTRIBUTING.md` document in the MFEM repository.
+
+2. To respond to reviewer feedback in a timely manner.
+
+3. Authors are encouraged to perform testing and inform the reviewers about the results.
+
+4. Authors can use the "Reviewers" section of the GitHub PR page to suggest reviewers, but the "Assignees" section will show who the editor has assigned to do the reviews.
+
+5. To indicate when the PR is ready for review by adding the `ready-for-review` label.
+
+
 ### Pull Request Checklist
 
 Before a PR can be merged, it should satisfy the following:
@@ -450,7 +539,9 @@ Before a PR can be merged, it should satisfy the following:
      - [ ] The miniapps go at the end of the page, and are usually listed only under a specific "Application (PDE)" category.
      - [ ] Add a short description of the miniapp in the "Extensive Examples" section of `features.md`.
 - [ ] New capability:
-   - [ ] All significant new classes, methods and functions have Doxygen-style documentation in source comments.
+   - [ ] All new public, protected, and private classes, methods, data members, and functions have full Doxygen-style documentation in source comments. Documentation should include descriptions of member data, function arguments and return values, template parameters, and prerequisites for calling new functions.
+   - [ ] Pointer arguments and return values must specify whether ownership is being transferred or lent with the call.
+   - [ ] Any new functions should include descriptions of their intended use e.g. for internal use only, user-facing, etc., along with references to example code whenever possible/appropriate.
    - [ ] Consider adding new sample runs in existing examples to highlight the new capability.
    - [ ] Consider saving cool simulation pictures with the new capability in the Confluence gallery (LLNL only) or submitting them, via pull request, to the gallery section of the `mfem/web` repo.
    - [ ] If this is a major new feature, consider mentioning it in the short summary inside `README` *(rare)*.
@@ -461,6 +552,7 @@ Before a PR can be merged, it should satisfy the following:
 - [ ] (LLNL only) After merging:
    - [ ] Update internal tests to include the new features.
 
+
 ### Master/Next Workflow
 
 MFEM uses a `master`/`next`-branch workflow as described below:
@@ -552,8 +644,10 @@ MFEM uses a `master`/`next`-branch workflow as described below:
     - Update version and shortlinks in `src/index.md` and `src/download.md`.
     -  Use [cloc-1.62.pl](http://cloc.sourceforge.net/) and `ls -lh` to estimate the SLOC and the tarball size in `src/download.md`.
 
+
 ## LLNL Workflow
 
+
 ### Mirroring on Bitbucket
 
 - The GitHub `master` and `next` branches are mirrored to the LLNL institutional
@@ -573,16 +667,17 @@ MFEM uses a `master`/`next`-branch workflow as described below:
   - `mfem:gh-next` -- Bleeding-edge development version, may be broken, use at
      your own risk.
 
+
 ### Mirroring on GitLab
 
 - MFEM repository is also mirrored on the LLNL GitLab instance, in a
   semi-automated manner.
 
 - This instance is meant to complete CI testing with tests on Livermore
-  Computing systems. Gitlab pipeline status is reported in the corresponding
+  Computing systems. GitLab pipeline status is reported in the corresponding
   GitHub pull request.
 
-- In Gitlab pipelines, TPLs (dependencies) are built using Spack, driven by Uberenv.
+- In GitLab pipelines, TPLs (dependencies) are built using Spack, driven by Uberenv.
 
 - No change to the MFEM repo can be made on this instance.
 
@@ -595,6 +690,7 @@ In addition, developers can set local git hooks to run some quick checks on
 commit or push, see the [README](config/githooks/README.md) in the `config/githooks`
 directory.
 
+
 ### Linux and Mac smoke tests
 We use GitHub Actions to drive the default tests on the `master` and `next`
 branches. See the `.github/workflows` files and the logs at
@@ -606,6 +702,7 @@ constraint on jobs. Two virtual machines are configured - Mac (OS X) and Linux.
 - Tests on the `master` branch are triggered whenever a PR is issued on this branch.
 - Tests on the `next` branch are currently scheduled to run each night.
 
+
 ### Windows smoke test
 We use Appveyor to test building with the MS Visual C++ compiler in a Windows
 environment, as well as to test the CMake build. See the `.appveyor` file and the
@@ -615,6 +712,7 @@ build logs at
 CMake is used to generate the MSVC Project files and drive the build.  A release
 and debug build is performed with a simple run of `ex1` to verify the executable.
 
+
 ### Tests at LLNL
 
 - We mirror the `master` and `next` branches internally (to `gh-master` and

INSTALL

@@ -459,6 +459,10 @@ MFEM_USE_UMPIRE = YES/NO
    discovery, provision, and management of memory on machines with multiple
    memory devices like NUMA and GPUs.
 
+MFEM_USE_BENCHMARK = YES/NO
+   Enables support for Google Benchmark, a library to support the benchmarking
+   of functions, in the tests/benchmarks directory.
+
 MFEM_USE_HIOP = YES/NO
    Enable the usage of HiOp (https://github.com/LLNL/hiop) in MFEM. HiOp is an
    HPC solver for nonlinear optimization problems.
@@ -523,6 +527,11 @@ MFEM_USE_FMS = YES/NO
    convetion routines between FMS's FmsDataCollection structure and MFEM's
    DataCollection class, see the header file fem/fmsconvert.hpp.
 
+MFEM_USE_PARELAG = YES/NO
+   Enables the miniapps that use the ParELAG library. MFEM does not currently
+   use ParELAG. In fact, ParELAG is dependent on MFEM. Therefore, this option
+   currently only concerns the miniapps.
+
 MFEM_BUILD_TAG = (any value)
    An optional tag to characterize the build. Exported to config/config.mk.
    Can be used to identify the MFEM build from other makefiles.
@@ -549,7 +558,7 @@ The specific libraries and their options are:
   Options: HYPRE_OPT, HYPRE_LIB.
   Versions: HYPRE >= 2.10.0b  (HYPRE built without CUDA)
             HYPRE >= 2.20.0   (HYPRE built with '--enable-mixedint')
-            HYPRE >= 2.22.0   (HYPRE built with CUDA)
+            HYPRE >= 2.22.1   (HYPRE built with CUDA)
 
 - METIS, used when MFEM_USE_METIS = YES. If using METIS 5, set
   MFEM_USE_METIS_5 = YES (default is to use METIS 4).
@@ -692,7 +701,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
@@ -730,10 +739,10 @@ The specific libraries and their options are:
   Versions: libCEED >= 0.8.
 
 - RAJA (optional), used when MFEM_USE_RAJA = YES.
-  Beginning with MFEM v4.3, only RAJA v0.13.0+ is supported.
+  Beginning with MFEM v4.3, only RAJA v0.14.0+ is supported.
   URL: https://github.com/LLNL/RAJA
   Options: RAJA_DIR, RAJA_OPT, RAJA_LIB.
-  Versions: RAJA >= 0.13.0.
+  Versions: RAJA >= 0.14.0.
 
 - Caliper (optional), used when MFEM_USE_CALIPER = YES.
   URL: https://github.com/LLNL/Caliper
@@ -744,7 +753,12 @@ The specific libraries and their options are:
   Umpire requires camp when the Umpire version is >= 3.0.0.
   URL: https://github.com/LLNL/Umpire
   Options: UMPIRE_DIR, UMPIRE_OPT, UMPIRE_LIB.
-  Versions: Umpire >= 2.0.0.
+  Versions: Umpire >= 3.0.0.
+
+- Benchmark, used when MFEM_USE_BENCHMARK = YES.
+  URL: https://github.com/google/benchmark
+  Options: BENCHMARK_DIR, BENCHMARK_LIB.
+  Versions: Benchmark >= 1.5.6.
 
 - MPFR (optional), used when MFEM_USE_MPFR = YES.
   URL: http://mpfr.org, it depends on the GMP library: https://gmplib.org
@@ -766,6 +780,10 @@ The specific libraries and their options are:
   Options: FMS_OPT, FMS_LIB.
   Versions: FMS >= 0.2.
 
+- ParELAG, used when MFEM_USE_PARELAG = YES.
+  URL: https://github.com/LLNL/parelag
+  Options: PARELAG_DIR, PARELAG_OPT, PARELAG_LIB.
+
 Building with CMake
 ===================
 The MFEM build system consists of two steps: configuration and compilation.
@@ -898,6 +916,8 @@ MFEM_USE_UMPIRE
 MFEM_USE_SIDRE
 MFEM_USE_CALIPER
 MFEM_USE_FMS
+MFEM_USE_BENCHMARK
+MFEM_USE_PARELAG
 
 The following options are CMake specific:
 
@@ -953,6 +973,8 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - AXOM - Used when MFEM_USE_SIDRE is enabled
  - CALIPER
  - FMS
+ - BENCHMARK
+ - ParELAG
 
 The following built-in CMake packages are also used:
 

config/XSDKDefaults.cmake

@@ -283,3 +283,11 @@ ENDIF()
 IF (DEFINED TPL_ENABLE_UMPIRE)
   SET(MFEM_USE_UMPIRE ${TPL_ENABLE_UMPIRE} CACHE BOOL "Enable Umpire" FORCE)
 ENDIF()
+
+IF (DEFINED TPL_ENABLE_BENCHMARK)
+  SET(MFEM_USE_BENCHMARK ${TPL_ENABLE_BENCHMARK} CACHE BOOL "Enable Google-Benchmark" FORCE)
+ENDIF()
+
+IF (DEFINED TPL_ENABLE_PARELAG)
+  SET(MFEM_USE_PARELAG ${TPL_ENABLE_PARELAG} CACHE BOOL "Enable ParELAG" FORCE)
+ENDIF()

config/cmake/MFEMConfig.cmake.in

@@ -55,6 +55,8 @@ set(MFEM_USE_UMPIRE @MFEM_USE_UMPIRE@)
 set(MFEM_USE_SIMD @MFEM_USE_SIMD@)
 set(MFEM_USE_ADIOS2 @MFEM_USE_ADIOS2@)
 set(MFEM_USE_CALIPER @MFEM_USE_CALIPER@)
+set(MFEM_USE_BENCHMARK @MFEM_USE_BENCHMARK@)
+set(MFEM_USE_PARELAG @MFEM_USE_PARELAG@)
 
 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
 set(MFEM_CXX_FLAGS "@CMAKE_CXX_FLAGS@")

config/cmake/config.hpp.in

@@ -175,4 +175,7 @@
 // Enable interface to the MKL CPardiso library.
 #cmakedefine MFEM_USE_MKL_CPARDISO
 
+// Enable MFEM functionality based on the Google Benchmark library.
+#cmakedefine MFEM_USE_BENCHMARK
+
 #endif // MFEM_CONFIG_HEADER

config/cmake/modules/FindBenchmark.cmake

@@ -0,0 +1,22 @@
+# 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.
+
+# Defines the following variables:
+#   - BENCHMARK_FOUND
+#   - BENCHMARK_LIBRARIES
+#   - BENCHMARK_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(Benchmark BENCHMARK BENCHMARK_DIR
+        "include" "benchmark/benchmark.h"
+        "lib" "benchmark"
+        "Paths to headers required by Google Benchmark."
+        "Libraries required by Google Benchmark.")

config/cmake/modules/FindPARELAG.cmake

@@ -0,0 +1,19 @@
+# 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.
+
+# Defines the following variables:
+#   - PARELAG_FOUND
+#   - PARELAG_LIBRARIES
+#   - PARELAG_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(PARELAG PARELAG PARELAG_DIR "" "" "" ""
+  "Paths to headers required by ParELAG." "Libraries required by ParELAG.")

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -100,6 +100,8 @@ macro(add_mfem_examples EXE_SRCS)
 
     string(REPLACE ".cpp" "" EXE_NAME "${EXE_PREFIX}${SRC_FILENAME}")
     mfem_add_executable(${EXE_NAME} ${SRC_FILE})
+    install(TARGETS ${EXE_NAME}
+            RUNTIME DESTINATION examples)
     add_dependencies(${MFEM_ALL_EXAMPLES_TARGET_NAME} ${EXE_NAME})
     if (EXE_NEEDED_BY)
       add_dependencies(${EXE_NEEDED_BY} ${EXE_NAME})
@@ -764,7 +766,7 @@ function(mfem_export_mk_files)
       MFEM_USE_GNUTLS MFEM_USE_GSLIB MFEM_USE_NETCDF MFEM_USE_PETSC
       MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT MFEM_USE_PUMI
       MFEM_USE_CUDA MFEM_USE_OCCA MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD
-      MFEM_USE_ADIOS2)
+      MFEM_USE_ADIOS2 MFEM_USE_BENCHMARK MFEM_USE_PARELAG)
   foreach(var ${CONFIG_MK_BOOL_VARS})
     if (${var})
       set(${var} YES)

config/config.hpp.in

@@ -180,4 +180,7 @@
 // Enable interface to the MKL CPardiso library.
 // #define MFEM_USE_MKL_CPARDISO
 
+// Enable functionality based on the Google Benchmark library.
+// #define MFEM_USE_BENCHMARK
+
 #endif // MFEM_CONFIG_HEADER

config/config.mk.in

@@ -58,6 +58,8 @@ MFEM_USE_UMPIRE        = @MFEM_USE_UMPIRE@
 MFEM_USE_SIMD          = @MFEM_USE_SIMD@
 MFEM_USE_ADIOS2        = @MFEM_USE_ADIOS2@
 MFEM_USE_MKL_CPARDISO  = @MFEM_USE_MKL_CPARDISO@
+MFEM_USE_BENCHMARK     = @MFEM_USE_BENCHMARK@
+MFEM_USE_PARELAG       = @MFEM_USE_PARELAG@
 
 # Compiler, compile options, and link options
 MFEM_CXX       = @MFEM_CXX@

config/defaults.cmake

@@ -58,6 +58,8 @@ option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" OFF)
 option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)
 option(MFEM_USE_CALIPER "Enable Caliper support" OFF)
 option(MFEM_USE_MKL_CPARDISO "Enable MKL CPardiso" OFF)
+option(MFEM_USE_BENCHMARK "Enable Google Benchmark" OFF)
+option(MFEM_USE_PARELAG "Enable ParELAG" OFF)
 
 # Optional overrides for autodetected MPIEXEC and MPIEXEC_NUMPROC_FLAG
 # set(MFEM_MPIEXEC "mpirun" CACHE STRING "Command for running MPI tests")
@@ -74,6 +76,7 @@ set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 option(MFEM_ENABLE_TESTING "Enable the ctest framework for testing" ON)
 option(MFEM_ENABLE_EXAMPLES "Build all of the examples" OFF)
 option(MFEM_ENABLE_MINIAPPS "Build all of the miniapps" OFF)
+option(MFEM_ENABLE_GOOGLE_BENCHMARKS "Build all of the Google benchmarks" OFF)
 
 # Setting CXX/MPICXX on the command line or in user.cmake will overwrite the
 # autodetected C++ compiler.
@@ -224,6 +227,16 @@ set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
 set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
 set(UMPIRE_DIR "${MFEM_DIR}/../umpire" CACHE PATH "Path to Umpire")
 set(CALIPER_DIR "${MFEM_DIR}/../caliper" CACHE PATH "Path to Caliper")
+set(BENCHMARK_DIR "${MFEM_DIR}/../google-benchmark" CACHE PATH
+    "Path to Google Benchmark")
+
+# Provide paths, since ParELAG is dependent on MFEM and MFEM needs to be
+# compiled (or at least cmake needs to succeed) before compiling ParELAG.
+set(PARELAG_DIR "${MFEM_DIR}/../parelag" CACHE PATH "Path to ParELAG")
+set(PARELAG_INCLUDE_DIRS "${PARELAG_DIR}/src;${PARELAG_DIR}/build/src" CACHE
+    STRING "Path to ParELAG headers.")
+set(PARELAG_LIBRARIES "${PARELAG_DIR}/build/src/libParELAG.a" CACHE STRING
+    "The ParELAG library.")
 
 set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
 set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")

config/defaults.mk

@@ -61,7 +61,7 @@ HIP_XLINKER   = -Wl,
 
 ifneq ($(NOTMAC),)
    AR      = ar
-   ARFLAGS = cruv
+   ARFLAGS = crv
    RANLIB  = ranlib
    PICFLAG = $(XCOMPILER)-fPIC
    SO_EXT  = so
@@ -73,7 +73,7 @@ ifneq ($(NOTMAC),)
 else
    # Silence "has no symbols" warnings on Mac OS X
    AR      = ar
-   ARFLAGS = Scruv
+   ARFLAGS = Scrv
    RANLIB  = ranlib -no_warning_for_no_symbols
    PICFLAG = $(XCOMPILER)-fPIC
    SO_EXT  = dylib
@@ -151,6 +151,8 @@ MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
 MFEM_USE_MKL_CPARDISO  = NO
+MFEM_USE_BENCHMARK     = NO
+MFEM_USE_PARELAG       = NO
 
 # MPI library compile and link flags
 # These settings are used only when building MFEM with MPI + HIP
@@ -429,6 +431,11 @@ CALIPER_DIR = @MFEM_DIR@/../caliper
 CALIPER_OPT = -I$(CALIPER_DIR)/include
 CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib64 -lcaliper
 
+# BENCHMARK library configuration
+BENCHMARK_DIR = @MFEM_DIR@/../google-benchmark
+BENCHMARK_OPT = -I$(BENCHMARK_DIR)/include
+BENCHMARK_LIB = -L$(BENCHMARK_DIR)/lib -lbenchmark -lpthread
+
 # libCEED library configuration
 CEED_DIR ?= @MFEM_DIR@/../libCEED
 CEED_OPT = -I$(CEED_DIR)/include
@@ -459,6 +466,11 @@ MKL_CPARDISO_LIB = $(XLINKER)-rpath,$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR)\
    -L$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR) -l$(MKL_MPI_WRAPPER)\
    -lmkl_intel_lp64 -lmkl_sequential -lmkl_core
 
+# PARELAG library configuration
+PARELAG_DIR = @MFEM_DIR@/../parelag
+PARELAG_OPT = -I$(PARELAG_DIR)/src -I$(PARELAG_DIR)/build/src
+PARELAG_LIB = -L$(PARELAG_DIR)/build/src -lParELAG
+
 # If YES, enable some informational messages
 VERBOSE = NO
 

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
 

data/inline-pyramid.mesh

@@ -0,0 +1,9 @@
+MFEM INLINE mesh v1.0
+
+type = pyramid
+nx = 4
+ny = 4
+nz = 4
+sx = 1.0
+sy = 1.0
+sz = 1.0

data/octahedron.mesh

@@ -0,0 +1,43 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+2
+1 7 4 3 2 1 0
+1 7 1 2 3 4 5
+
+boundary
+8
+1 2 0 2 1
+2 2 0 3 2
+3 2 0 4 3
+4 2 0 1 4
+5 2 1 2 5
+6 2 2 3 5
+7 2 3 4 5
+8 2 4 1 5
+
+vertices
+6
+3
+0 0 -1
+1 0 0
+0 1 0
+-1 0 0
+0 -1 0
+0 0 1

data/ref-pyramid.mesh

@@ -0,0 +1,38 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+1
+1 7 0 1 2 3 4
+
+boundary
+5
+1 3 3 2 1 0
+2 2 0 1 4
+3 2 1 2 4
+4 2 2 3 4
+5 2 3 0 4
+
+vertices
+5
+3
+0 0 0
+1 0 0
+1 1 0
+0 1 0
+0 0 1

doc/CodeDocumentation.conf.in

@@ -765,6 +765,7 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                          @MFEM_SOURCE_DIR@/linalg \
                          @MFEM_SOURCE_DIR@/mesh \
                          @MFEM_SOURCE_DIR@/fem \
+                         @MFEM_SOURCE_DIR@/fem/fe \
                          @MFEM_SOURCE_DIR@/examples \
                          @MFEM_SOURCE_DIR@/examples/caliper \
                          @MFEM_SOURCE_DIR@/examples/amgx \
@@ -786,7 +787,8 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                          @MFEM_SOURCE_DIR@/miniapps/shifted \
                          @MFEM_SOURCE_DIR@/miniapps/solvers \
                          @MFEM_SOURCE_DIR@/miniapps/tools \
-                         @MFEM_SOURCE_DIR@/miniapps/toys
+                         @MFEM_SOURCE_DIR@/miniapps/toys \
+                         @MFEM_SOURCE_DIR@/miniapps/parelag
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses

examples/CMakeLists.txt

@@ -37,6 +37,7 @@ list(APPEND ALL_EXE_SRCS
   ex27.cpp
   ex28.cpp
   ex29.cpp
+  ex30.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -70,6 +71,7 @@ if (MFEM_USE_MPI)
     ex27p.cpp
     ex28p.cpp
     ex29p.cpp
+    ex30p.cpp
     )
 endif()
 

examples/caliper/ex1.cpp

@@ -206,9 +206,9 @@ int main(int argc, char *argv[])
    cout << "Size of linear system: " << A->Height() << endl;
 
    // 11. Solve the linear system A X = B.
-   MFEM_PERF_BEGIN("Solve A X=B");
    if (!pa)
    {
+      MFEM_PERF_SCOPE("Solve A X=B (FA)");
 #ifndef MFEM_USE_SUITESPARSE
       // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
       GSSmoother M((SparseMatrix&)(*A));
@@ -223,6 +223,7 @@ int main(int argc, char *argv[])
    }
    else // Jacobi preconditioning in partial assembly mode
    {
+      MFEM_PERF_SCOPE("Solve A X=B (PA)");
       if (UsesTensorBasis(fespace))
       {
          OperatorJacobiSmoother M(a, ess_tdof_list);
@@ -233,7 +234,6 @@ int main(int argc, char *argv[])
          CG(*A, B, X, 1, 400, 1e-12, 0.0);
       }
    }
-   MFEM_PERF_END("Solve A X=B");
    // 12. Recover the solution as a finite element grid function.
    a.RecoverFEMSolution(X, b, x);
 

examples/caliper/ex1p.cpp

@@ -231,28 +231,29 @@ int main(int argc, char *argv[])
    // 13. Solve the linear system A X = B.
    //     * With full assembly, use the BoomerAMG preconditioner from hypre.
    //     * With partial assembly, use Jacobi smoothing, for now.
-   MFEM_PERF_BEGIN("Solve A X = B");
-   Solver *prec = NULL;
-   if (pa)
    {
-      if (UsesTensorBasis(fespace))
+      MFEM_PERF_SCOPE("Solve A X=B");
+      Solver *prec = NULL;
+      if (pa)
       {
-         prec = new OperatorJacobiSmoother(a, ess_tdof_list);
+         if (UsesTensorBasis(fespace))
+         {
+            prec = new OperatorJacobiSmoother(a, ess_tdof_list);
+         }
       }
+      else
+      {
+         prec = new HypreBoomerAMG;
+      }
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(2000);
+      cg.SetPrintLevel(1);
+      if (prec) { cg.SetPreconditioner(*prec); }
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+      delete prec;
    }
-   else
-   {
-      prec = new HypreBoomerAMG;
-   }
-   CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
-   cg.SetPrintLevel(1);
-   if (prec) { cg.SetPreconditioner(*prec); }
-   cg.SetOperator(*A);
-   cg.Mult(B, X);
-   delete prec;
-   MFEM_PERF_END("Solve A X = B");
    // 14. Recover the parallel grid function corresponding to X. This is the
    //     local finite element solution on each processor.
    a.RecoverFEMSolution(X, b, x);

examples/ex1.cpp

@@ -9,6 +9,7 @@
 //               ex1 -m ../data/fichera.mesh
 //               ex1 -m ../data/fichera-mixed.mesh
 //               ex1 -m ../data/toroid-wedge.mesh
+//               ex1 -m ../data/octahedron.mesh -o 1
 //               ex1 -m ../data/periodic-annulus-sector.msh
 //               ex1 -m ../data/periodic-torus-sector.msh
 //               ex1 -m ../data/square-disc-p2.vtk -o 2

examples/ex13p.cpp

@@ -118,7 +118,6 @@ int main(int argc, char *argv[])
    {
       pmesh->UniformRefinement();
    }
-   pmesh->ReorientTetMesh();
 
    // 7. Define a parallel finite element space on the parallel mesh. Here we
    //    use the Nedelec finite elements of the specified order.

examples/ex14.cpp

@@ -135,8 +135,8 @@ int main(int argc, char *argv[])
    a->AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
    if (eta > 0)
    {
-      a->AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
-      a->AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+      a->AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
+      a->AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
    }
    a->Assemble();
    a->Finalize();

examples/ex14p.cpp

@@ -199,8 +199,8 @@ int main(int argc, char *argv[])
    a->AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
    if (eta > 0)
    {
-      a->AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
-      a->AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+      a->AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
+      a->AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
    }
    a->Assemble();
    a->Finalize();
@@ -221,7 +221,7 @@ int main(int argc, char *argv[])
    {
       HyprePCG pcg(*A);
       pcg.SetTol(1e-12);
-      pcg.SetMaxIter(200);
+      pcg.SetMaxIter(500);
       pcg.SetPrintLevel(2);
       pcg.SetPreconditioner(*amg);
       pcg.Mult(*B, *X);
@@ -232,7 +232,7 @@ int main(int argc, char *argv[])
       GMRESSolver gmres(MPI_COMM_WORLD);
       gmres.SetAbsTol(0.0);
       gmres.SetRelTol(1e-12);
-      gmres.SetMaxIter(200);
+      gmres.SetMaxIter(500);
       gmres.SetKDim(10);
       gmres.SetPrintLevel(1);
       gmres.SetOperator(*A);

examples/ex16.cpp

@@ -24,7 +24,10 @@
 //               class ConductionOperator defining C(u)), as well as their
 //               implicit time integration. Note that implementing the method
 //               ConductionOperator::ImplicitSolve is the only requirement for
-//               high-order implicit (SDIRK) time integration.
+//               high-order implicit (SDIRK) time integration. In this example,
+//               the diffusion operator is linearized by evaluating with the
+//               lagged solution from the previous timestep, so there is only
+//               a linear solve.
 //
 //               We recommend viewing examples 2, 9 and 10 before viewing this
 //               example.
@@ -326,8 +329,8 @@ ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 {
    // Compute:
-   //    du_dt = M^{-1}*-K(u)
-   // for du_dt
+   //    du_dt = M^{-1}*-Ku
+   // for du_dt, where K is linearized by using u from the previous timestep
    Kmat.Mult(u, z);
    z.Neg(); // z = -z
    M_solver.Mult(z, du_dt);
@@ -338,7 +341,7 @@ void ConductionOperator::ImplicitSolve(const double dt,
 {
    // Solve the equation:
    //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
+   // for du_dt, where K is linearized by using u from the previous timestep
    if (!T)
    {
       T = Add(1.0, Mmat, dt, Kmat);

examples/ex16p.cpp

@@ -24,8 +24,11 @@
 //               class ConductionOperator defining C(u)), as well as their
 //               implicit time integration. Note that implementing the method
 //               ConductionOperator::ImplicitSolve is the only requirement for
-//               high-order implicit (SDIRK) time integration. Optional saving
-//               with ADIOS2 (adios2.readthedocs.io) is also illustrated.
+//               high-order implicit (SDIRK) time integration. In this example,
+//               the diffusion operator is linearized by evaluating with the
+//               lagged solution from the previous timestep, so there is only
+//               a linear solve. Optional saving with ADIOS2
+//               (adios2.readthedocs.io) is also illustrated.
 //
 //               We recommend viewing examples 2, 9 and 10 before viewing this
 //               example.
@@ -420,8 +423,8 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 {
    // Compute:
-   //    du_dt = M^{-1}*-K(u)
-   // for du_dt
+   //    du_dt = M^{-1}*-Ku
+   // for du_dt, where K is linearized by using u from the previous timestep
    Kmat.Mult(u, z);
    z.Neg(); // z = -z
    M_solver.Mult(z, du_dt);
@@ -432,7 +435,7 @@ void ConductionOperator::ImplicitSolve(const double dt,
 {
    // Solve the equation:
    //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
+   // for du_dt, where K is linearized by using u from the previous timestep
    if (!T)
    {
       T = Add(1.0, Mmat, dt, Kmat);

examples/ex1p.cpp

@@ -9,6 +9,7 @@
 //               mpirun -np 4 ex1p -m ../data/fichera.mesh
 //               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
 //               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
+//               mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
 //               mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
 //               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2

examples/ex22.cpp

@@ -13,6 +13,8 @@
 //               ex22 -m ../data/inline-hex.mesh -o 2 -p 1
 //               ex22 -m ../data/inline-hex.mesh -o 2 -p 2
 //               ex22 -m ../data/inline-hex.mesh -o 2 -p 2 -pa
+//               ex22 -m ../data/inline-wedge.mesh -o 1
+//               ex22 -m ../data/inline-pyramid.mesh -o 1
 //               ex22 -m ../data/star.mesh -r 1 -o 2 -sigma 10.0
 //
 // Device sample runs:

examples/ex22p.cpp

@@ -13,6 +13,8 @@
 //               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2 -p 1
 //               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2 -p 2
 //               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 1 -p 2 -pa
+//               mpirun -np 4 ex22p -m ../data/inline-wedge.mesh -o 1
+//               mpirun -np 4 ex22p -m ../data/inline-pyramid.mesh -o 1
 //               mpirun -np 4 ex22p -m ../data/star.mesh -o 2 -sigma 10.0
 //
 // Device sample runs:

examples/ex24.cpp

@@ -113,7 +113,6 @@ int main(int argc, char *argv[])
          mesh->UniformRefinement();
       }
    }
-   mesh->ReorientTetMesh();
 
    // 5. Define a finite element space on the mesh. Here we use Nedelec or
    //    Raviart-Thomas finite elements of the specified order.

examples/ex24p.cpp

@@ -141,7 +141,6 @@ int main(int argc, char *argv[])
          pmesh->UniformRefinement();
       }
    }
-   pmesh->ReorientTetMesh();
 
    // 7. Define a parallel finite element space on the parallel mesh. Here we
    //    use Nedelec or Raviart-Thomas finite elements of the specified order.

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 &),
@@ -277,10 +277,8 @@ int main(int argc, char *argv[])
       mesh->UniformRefinement();
    }
 
-   // 6. Reorient mesh in case of a tet mesh
-   mesh->ReorientTetMesh();
-
-   // Set element attributes in order to distinguish elements in the PML region
+   // 6. Set element attributes in order to distinguish elements in the
+   //    PML region
    pml->SetAttributes(mesh);
 
    // 7. Define a finite element space on the mesh. Here we use the Nedelec

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 &),
@@ -316,9 +316,6 @@ int main(int argc, char *argv[])
       }
    }
 
-   // 7a. Reorient mesh in case of a tet mesh
-   pmesh->ReorientTetMesh();
-
    // 8. Set element attributes in order to distinguish elements in the PML
    pml->SetAttributes(pmesh);
 

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;

examples/ex3.cpp

@@ -16,6 +16,8 @@
 //               ex3 -m ../data/beam-hex-nurbs.mesh
 //               ex3 -m ../data/amr-hex.mesh
 //               ex3 -m ../data/fichera-amr.mesh
+//               ex3 -m ../data/ref-prism.mesh -o 1
+//               ex3 -m ../data/octahedron.mesh -o 1
 //               ex3 -m ../data/star-surf.mesh -o 1
 //               ex3 -m ../data/mobius-strip.mesh -f 0.1
 //               ex3 -m ../data/klein-bottle.mesh -f 0.1
@@ -113,7 +115,6 @@ int main(int argc, char *argv[])
          mesh->UniformRefinement();
       }
    }
-   mesh->ReorientTetMesh();
 
    // 5. Define a finite element space on the mesh. Here we use the Nedelec
    //    finite elements of the specified order.

examples/ex30.cpp

@@ -0,0 +1,195 @@
+//                            MFEM Example 30
+//
+// Compile with: make ex30
+//
+// Sample runs:  ex30 -m ../data/square-disc.mesh -o 1
+//               ex30 -m ../data/square-disc.mesh -o 2
+//               ex30 -m ../data/square-disc.mesh -o 2 -me 1e3
+//               ex30 -m ../data/square-disc-nurbs.mesh -o 2
+//               ex30 -m ../data/star.mesh -o 2 -eo 4
+//               ex30 -m ../data/fichera.mesh -o 2 -me 1e4
+//               ex30 -m ../data/disc-nurbs.mesh -o 2
+//               ex30 -m ../data/ball-nurbs.mesh -o 2 -eo 3 -e 1e-2 -me 1e4
+//               ex30 -m ../data/star-surf.mesh -o 2
+//               ex30 -m ../data/square-disc-surf.mesh -o 2
+//               ex30 -m ../data/amr-quad.mesh -l 2
+//
+// Description:  This is an example of adaptive mesh refinement preprocessing
+//               which lowers the data oscillation [1] to a user-defined
+//               relative threshold. There is no PDE being solved.
+//
+//               MFEM's capability to work with both conforming and
+//               nonconforming meshes is demonstrated in example 6. In some
+//               problems, the material data or loading data is not sufficiently
+//               resolved on the initial mesh. This missing fine scale data
+//               reduces the accuracy of the solution as well as the accuracy
+//               of some local error estimators. By preprocessing the mesh
+//               before solving the PDE, many issues can be avoided.
+//
+//               [1] Morin, P., Nochetto, R. H., & Siebert, K. G. (2000).
+//                   Data oscillation and convergence of adaptive FEM. SIAM
+//                   Journal on Numerical Analysis, 38(2), 466-488.
+//
+//               [2] Mitchell, W. F. (2013). A collection of 2D elliptic
+//                   problems for testing adaptive grid refinement algorithms.
+//                   Applied mathematics and computation, 220, 350-364.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Piecewise-affine function which is sometimes mesh-conforming
+double affine_function(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   if (x < 0.0)
+   {
+      return 1.0 + x + y;
+   }
+   else
+   {
+      return 1.0;
+   }
+}
+
+// Piecewise-constant function which is never mesh-conforming
+double jump_function(const Vector &p)
+{
+   if (p.Normlp(2.0) > 0.4 && p.Normlp(2.0) < 0.6) { return 1.0; }
+   return 5.0;
+}
+
+// Singular function derived from the Laplacian of the "steep wavefront"
+// problem in [2].
+double singular_function(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   double alpha = 1000.0;
+   double xc = 0.75, yc = 0.5;
+   double r0 = 0.7;
+   double r = sqrt(pow(x - xc,2.0) + pow(y - yc,2.0));
+   double num = - ( alpha - pow(alpha,3) * (pow(r,2) - pow(r0,2)) );
+   double denom = pow(r * ( pow(alpha,2) * pow(r0,2) + pow(alpha,2) * pow(r,2) \
+                            - 2 * pow(alpha,2) * r0 * r + 1.0 ),2);
+   denom = max(denom,1e-8);
+   return num / denom;
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int nc_limit = 1;
+   int max_elems = 1e5;
+   double double_max_elems = double(max_elems);
+   bool visualization = true;
+   double osc_threshold = 1e-3;
+   int enriched_order = 5;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&nc_limit, "-l", "--nc-limit",
+                  "Maximum level of hanging nodes.");
+   args.AddOption(&double_max_elems, "-me", "--max-elems",
+                  "Stop after reaching this many elements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&osc_threshold, "-e", "--error",
+                  "relative data oscillation threshold.");
+   args.AddOption(&enriched_order, "-eo", "--enriched_order",
+                  "Enriched quadrature order.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   max_elems = int(double_max_elems);
+   Mesh mesh(mesh_file, 1, 1);
+
+   // 2. Since a NURBS mesh can currently only be refined uniformly, we need to
+   //    convert it to a piecewise-polynomial curved mesh. First we refine the
+   //    NURBS mesh a bit more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext)
+   {
+      for (int i = 0; i < 2; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.SetCurvature(2);
+   }
+
+   // 3. Define functions and refiner.
+   FunctionCoefficient affine_coeff(affine_function);
+   FunctionCoefficient jump_coeff(jump_function);
+   FunctionCoefficient singular_coeff(singular_function);
+   CoefficientRefiner coeffrefiner(affine_coeff, order);
+
+   // 4. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost, visport);
+   }
+
+   // 5. Define custom integration rule (optional).
+   const IntegrationRule *irs[Geometry::NumGeom];
+   int order_quad = 2*order + enriched_order;
+   for (int i = 0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+
+   // 6. Apply custom refiner settings.
+   coeffrefiner.SetIntRule(irs);
+   coeffrefiner.SetMaxElements(max_elems);
+   coeffrefiner.SetThreshold(osc_threshold);
+   coeffrefiner.SetNCLimit(nc_limit);
+   coeffrefiner.PrintWarnings();
+
+   // 7. Preprocess mesh to control osc (piecewise-affine function).
+   //    This is mostly just a verification check. The oscillation should
+   //    be zero if the function is mesh-conforming and order > 0.
+   coeffrefiner.PreprocessMesh(mesh);
+
+   mfem::out  << "\n";
+   mfem::out << "Function 0 (affine) \n";
+   mfem::out << "Number of Elements " << mesh.GetNE() << "\n";
+   mfem::out << "Osc error " << coeffrefiner.GetOsc() << "\n";
+
+   // 8. Preprocess mesh to control osc (jump function).
+   coeffrefiner.ResetCoefficient(jump_coeff);
+   coeffrefiner.PreprocessMesh(mesh);
+
+   mfem::out  << "\n";
+   mfem::out << "Function 1 (discontinuous) \n";
+   mfem::out << "Number of Elements " << mesh.GetNE() << "\n";
+   mfem::out << "Osc error " << coeffrefiner.GetOsc() << "\n";
+
+   // 9. Preprocess mesh to control osc (singular function).
+   coeffrefiner.ResetCoefficient(singular_coeff);
+   coeffrefiner.PreprocessMesh(mesh);
+
+   mfem::out  << "\n";
+   mfem::out << "Function 2 (singular) \n";
+   mfem::out << "Number of Elements " << mesh.GetNE() << "\n";
+   mfem::out << "Osc error " << coeffrefiner.GetOsc() << "\n";
+
+   sol_sock.precision(8);
+   sol_sock << "mesh\n" << mesh << flush;
+
+   return 0;
+}

examples/ex30p.cpp

@@ -0,0 +1,241 @@
+//                        MFEM Example 30 - Parallel Version
+//
+// Compile with: make ex30p
+//
+// Sample runs:  mpirun -np 4 ex30p -m ../data/square-disc.mesh -o 1
+//               mpirun -np 4 ex30p -m ../data/square-disc.mesh -o 2
+//               mpirun -np 4 ex30p -m ../data/square-disc.mesh -o 2 -me 1e3
+//               mpirun -np 4 ex30p -m ../data/square-disc-nurbs.mesh -o 2
+//               mpirun -np 4 ex30p -m ../data/star.mesh -o 2 -eo 4
+//               mpirun -np 4 oscp -m ../data/fichera.mesh -o 2 -me 1e4
+//               mpirun -np 4 ex30p -m ../data/disc-nurbs.mesh -o 2
+//               mpirun -np 4 ex30p -m ../data/ball-nurbs.mesh -o 2 -eo 3 -e 1e-2
+//               mpirun -np 4 ex30p -m ../data/star-surf.mesh -o 2
+//               mpirun -np 4 ex30p -m ../data/square-disc-surf.mesh -o 2
+//               mpirun -np 4 ex30p -m ../data/amr-quad.mesh -l 2
+//
+// Description:  This is an example of adaptive mesh refinement preprocessing
+//               which lowers the data oscillation [1] to a user-defined
+//               relative threshold. There is no PDE being solved.
+//
+//               MFEM's capability to work with both conforming and
+//               nonconforming meshes is demonstrated in example 6. In some
+//               problems, the material data or loading data is not sufficiently
+//               resolved on the initial mesh. This missing fine scale data
+//               reduces the accuracy of the solution as well as the accuracy
+//               of some local error estimators. By preprocessing the mesh
+//               before solving the PDE, many issues can be avoided.
+//
+//               [1] Morin, P., Nochetto, R. H., & Siebert, K. G. (2000).
+//                   Data oscillation and convergence of adaptive FEM. SIAM
+//                   Journal on Numerical Analysis, 38(2), 466-488.
+//
+//               [2] Mitchell, W. F. (2013). A collection of 2D elliptic
+//                   problems for testing adaptive grid refinement algorithms.
+//                   Applied mathematics and computation, 220, 350-364.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Piecewise-affine function which is sometimes mesh-conforming
+double affine_function(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   if (x < 0.0)
+   {
+      return 1.0 + x + y;
+   }
+   else
+   {
+      return 1.0;
+   }
+}
+
+// Piecewise-constant function which is never mesh-conforming
+double jump_function(const Vector &p)
+{
+   if (p.Normlp(2.0) > 0.4 && p.Normlp(2.0) < 0.6) { return 1.0; }
+   return 5.0;
+}
+
+// Singular function derived from the Laplacian of the "steep wavefront"
+// problem in [2].
+double singular_function(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   double alpha = 1000.0;
+   double xc = 0.75, yc = 0.5;
+   double r0 = 0.7;
+   double r = sqrt(pow(x - xc,2.0) + pow(y - yc,2.0));
+   double num = - ( alpha - pow(alpha,3) * (pow(r,2) - pow(r0,2)) );
+   double denom = pow(r * ( pow(alpha,2) * pow(r0,2) + pow(alpha,2) * pow(r,2) \
+                            - 2 * pow(alpha,2) * r0 * r + 1.0 ),2);
+   denom = max(denom,1e-8);
+   return num / denom;
+}
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int nc_limit = 1;
+   int max_elems = 1e5;
+   double double_max_elems = double(max_elems);
+   bool visualization = true;
+   bool nc_simplices = true;
+   double osc_threshold = 1e-3;
+   int enriched_order = 5;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&nc_limit, "-l", "--nc-limit",
+                  "Maximum level of hanging nodes.");
+   args.AddOption(&double_max_elems, "-me", "--max-elems",
+                  "Stop after reaching this many elements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&osc_threshold, "-e", "--error",
+                  "relative data oscillation threshold.");
+   args.AddOption(&enriched_order, "-eo", "--enriched_order",
+                  "Enriched quadrature order.");
+   args.AddOption(&nc_simplices, "-ns", "--nonconforming-simplices",
+                  "-cs", "--conforming-simplices",
+                  "For simplicial meshes, enable/disable nonconforming"
+                  " refinement");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   max_elems = int(double_max_elems);
+   Mesh mesh(mesh_file, 1, 1);
+
+   // 2. Since a NURBS mesh can currently only be refined uniformly, we need to
+   //    convert it to a piecewise-polynomial curved mesh. First we refine the
+   //    NURBS mesh a bit more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext)
+   {
+      for (int i = 0; i < 2; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.SetCurvature(2);
+   }
+
+   // 3. Make sure the mesh is in the non-conforming mode to enable local
+   //    refinement of quadrilaterals/hexahedra. Simplices can be refined
+   //    either in conforming or in non-conforming mode. The conforming
+   //    mode however does not support dynamic partitioning.
+   mesh.EnsureNCMesh(nc_simplices);
+
+   // 4. Define a parallel mesh by partitioning the serial mesh.
+   //    Once the parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // 5. Define functions and refiner.
+   FunctionCoefficient affine_coeff(affine_function);
+   FunctionCoefficient jump_coeff(jump_function);
+   FunctionCoefficient singular_coeff(singular_function);
+   CoefficientRefiner coeffrefiner(affine_coeff,order);
+
+   // 6. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost, visport);
+   }
+
+   // 7. Define custom integration rule (optional).
+   const IntegrationRule *irs[Geometry::NumGeom];
+   int order_quad = 2*order + enriched_order;
+   for (int i=0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+
+   // 8. Apply custom refiner settings.
+   coeffrefiner.SetIntRule(irs);
+   coeffrefiner.SetMaxElements(max_elems);
+   coeffrefiner.SetThreshold(osc_threshold);
+   coeffrefiner.SetNCLimit(nc_limit);
+   coeffrefiner.PrintWarnings();
+
+   // 9. Preprocess mesh to control osc (piecewise-affine function).
+   //    This is mostly just a verification check. The oscillation should
+   //    be zero if the function is mesh-conforming and order > 0.
+   coeffrefiner.PreprocessMesh(pmesh);
+
+   int globalNE = pmesh.GetGlobalNE();
+   double osc = coeffrefiner.GetOsc();
+   if (myid == 0)
+   {
+      mfem::out  << "\n";
+      mfem::out << "Function 0 (affine) \n";
+      mfem::out << "Number of Elements " << globalNE << "\n";
+      mfem::out << "Osc error " << osc << "\n";
+   }
+
+   // 10. Preprocess mesh to control osc (jump function).
+   coeffrefiner.ResetCoefficient(jump_coeff);
+   coeffrefiner.PreprocessMesh(pmesh);
+
+   globalNE = pmesh.GetGlobalNE();
+   osc = coeffrefiner.GetOsc();
+   if (myid == 0)
+   {
+      mfem::out  << "\n";
+      mfem::out << "Function 1 (discontinuous) \n";
+      mfem::out << "Number of Elements " << globalNE << "\n";
+      mfem::out << "Osc error " << osc << "\n";
+   }
+
+   // 11. Preprocess mesh to control osc (singular function).
+   coeffrefiner.ResetCoefficient(singular_coeff);
+   coeffrefiner.PreprocessMesh(pmesh);
+
+   globalNE = pmesh.GetGlobalNE();
+   osc = coeffrefiner.GetOsc();
+   if (myid == 0)
+   {
+      mfem::out  << "\n";
+      mfem::out << "Function 2 (singular) \n";
+      mfem::out << "Number of Elements " << globalNE << "\n";
+      mfem::out << "Osc error " << osc << "\n";
+   }
+
+   sol_sock.precision(8);
+   sol_sock << "parallel " << num_procs << " " << myid << "\n";
+   sol_sock << "mesh\n" << pmesh << flush;
+
+   MPI_Finalize();
+   return 0;
+}

examples/ex3p.cpp

@@ -16,6 +16,8 @@
 //               mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
 //               mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
 //               mpirun -np 4 ex3p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/ref-prism.mesh -o 1
+//               mpirun -np 4 ex3p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
 //               mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
 //               mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1
@@ -139,9 +141,7 @@ int main(int argc, char *argv[])
 
    // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
    //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
-   //    meshes need to be reoriented before we can define high-order Nedelec
-   //    spaces on them.
+   //    parallel mesh is defined, the serial mesh can be deleted.
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
    {
@@ -151,7 +151,6 @@ int main(int argc, char *argv[])
          pmesh->UniformRefinement();
       }
    }
-   pmesh->ReorientTetMesh();
 
    // 7. Define a parallel finite element space on the parallel mesh. Here we
    //    use the Nedelec finite elements of the specified order.

examples/ex4.cpp

@@ -19,6 +19,8 @@
 //               ex4 -m ../data/amr-hex.mesh
 //               ex4 -m ../data/amr-hex.mesh -o 2 -hb
 //               ex4 -m ../data/fichera-amr.mesh -o 2 -sc
+//               ex4 -m ../data/ref-prism.mesh -o 1
+//               ex4 -m ../data/octahedron.mesh -o 1
 //               ex4 -m ../data/star-surf.mesh -o 1
 //
 // Device sample runs:

examples/ex4p.cpp

@@ -19,6 +19,8 @@
 //               mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
+//               mpirun -np 4 ex4p -m ../data/ref-prism.mesh -o 1
+//               mpirun -np 4 ex4p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
 //
 // Device sample runs:
@@ -135,9 +137,7 @@ int main(int argc, char *argv[])
 
    // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
    //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
-   //    meshes need to be reoriented before we can define high-order Nedelec
-   //    spaces on them (this is needed in the ADS solver below).
+   //    parallel mesh is defined, the serial mesh can be deleted.
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
    {
@@ -147,7 +147,6 @@ int main(int argc, char *argv[])
          pmesh->UniformRefinement();
       }
    }
-   pmesh->ReorientTetMesh();
 
    // 7. Define a parallel finite element space on the parallel mesh. Here we
    //    use the Raviart-Thomas finite elements of the specified order.

examples/ex8p.cpp

@@ -106,7 +106,6 @@ int main(int argc, char *argv[])
          pmesh->UniformRefinement();
       }
    }
-   pmesh->ReorientTetMesh();
 
    // 6. Define the trial, interfacial (trace) and test DPG spaces:
    //    - The trial space, x0_space, contains the non-interfacial unknowns and

examples/jupyter/README.md

@@ -0,0 +1,31 @@
+# Jupyter Notebooks using xeus-cling
+
+[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mfem/mfem/master?filepath=examples%2Fjupyter%2Fex.ipynb)
+
+[xeus-cling](https://github.com/jupyter-xeus/xeus-cling) is a C++ Jupyter Kernel based on [cling](https://github.com/root-project/cling),
+which can be used to create interactive C++ MFEM and GLVis notebooks.
+
+Click on the `binder` button above for an interactive example.
+
+## Installing Locally
+
+In order to run notebooks locally you will need `xeus-cling` along with `mfem` and `xglvis`. We recommend you use
+[miniconda](https://docs.conda.io/en/latest/miniconda.html) or, if you already have it installed,
+[conda](https://docs.conda.io/projects/conda/en/latest/).
+
+1. Follow the install steps on https://github.com/jupyter-xeus/xeus-cling to install the C++ kernels
+2. Build and install a _shared_ version of mfem
+    * for example: `make serial SHARED=YES`
+3. Install [pyglvis](https://github.com/glvis/pyglvis)
+    * for the widget frontend
+4. Get [xeus-glvis](https://github.com/glvis/xeus-glvis) and `cp` the header to `{PREFIX}/glvis/xglvis.hpp`
+    * (this could be improved)
+
+## Running Locally
+
+Once you've installed Jupyter, the C++ Kernel, mfem, and glvis start the notebook server (`jupyter-notebook`)
+and open an existing example or a new `C++ 1x` kernel.
+
+You will _always_ need to `#pragma cling load("mfem")` and you may need to point the `cling` runtime at your
+mfem and/or glvis installs, do this with the
+`#pragma cling` [statements](https://xeus-cling.readthedocs.io/en/latest/build_options.html#using-third-party-libraries).

examples/jupyter/ex.ipynb

@@ -0,0 +1,155 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "owned-extraction",
+   "metadata": {},
+   "source": [
+    "## Load the MFEM library\n",
+    "\n",
+    "Any non-default libraries must be loaded before you can `#include` files that use them. For more info see the [xeus-cling help](https://xeus-cling.readthedocs.io/en/latest/build_options.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "waiting-portrait",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#pragma cling load(\"mfem\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "foreign-recycling",
+   "metadata": {},
+   "source": [
+    "## MFEM Example 1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "public-white",
+   "metadata": {},
+   "source": [
+    "This is the simplest MFEM example and a good starting point for new users. The example demonstrates the use of MFEM to define and solve an $H^1$ finite element discretization of the Laplace problem\n",
+    "\n",
+    "$$\n",
+    "-\\Delta u = 1\n",
+    "$$\n",
+    "\n",
+    "with homogeneous Dirichlet boundary conditions $u=0$.\n",
+    "\n",
+    "The example illustrates the use of the basic MFEM classes for defining the mesh, finite element space, as well as linear and bilinear forms corresponding to the left-hand side and right-hand side of the discrete linear system.\n",
+    "\n",
+    "Compare with MFEM's [ex1.cpp](https://github.com/mfem/mfem/blob/master/examples/ex1.cpp) and PyMFEM's [ex1.py](https://github.com/mfem/PyMFEM/blob/master/examples/ex1.py)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "protective-darkness",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#include <fstream>\n",
+    "#include <iostream>\n",
+    "#include <sstream>\n",
+    "\n",
+    "#include <mfem.hpp>\n",
+    "#include <glvis/xglvis.hpp>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "falling-monkey",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using namespace std;\n",
+    "using namespace mfem;\n",
+    "\n",
+    "Mesh mesh = Mesh::MakeCartesian2D(5, 5, Element::TRIANGLE);\n",
+    "mesh.UniformRefinement();\n",
+    "\n",
+    "H1_FECollection fec(2, mesh.Dimension());\n",
+    "\n",
+    "FiniteElementSpace fespace(&mesh, &fec);\n",
+    "cout << \"Number of finite element unknowns: \" << fespace.GetTrueVSize() << endl;\n",
+    "\n",
+    "Array<int> ess_tdof_list;\n",
+    "if (mesh.bdr_attributes.Size())\n",
+    "{\n",
+    "  Array<int> ess_bdr(mesh.bdr_attributes.Max());\n",
+    "  ess_bdr = 1;\n",
+    "  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);\n",
+    "}\n",
+    "\n",
+    "LinearForm b(&fespace);\n",
+    "ConstantCoefficient one(1.0);\n",
+    "b.AddDomainIntegrator(new DomainLFIntegrator(one));\n",
+    "b.Assemble();\n",
+    "\n",
+    "GridFunction x(&fespace);\n",
+    "x = 0.0;\n",
+    "\n",
+    "BilinearForm a(&fespace);\n",
+    "a.AddDomainIntegrator(new DiffusionIntegrator(one));\n",
+    "a.Assemble();\n",
+    "\n",
+    "OperatorPtr A;\n",
+    "Vector B, X;\n",
+    "a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);\n",
+    "\n",
+    "cout << \"Size of linear system: \" << A->Height() << endl;\n",
+    "\n",
+    "GSSmoother M((SparseMatrix&)(*A));\n",
+    "PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);\n",
+    "a.RecoverFEMSolution(X, b, x);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "hawaiian-republican",
+   "metadata": {},
+   "source": [
+    "## GLVis Visualization\n",
+    "\n",
+    "For now we save the computational mesh and finite element solution in a string and pass that to the glvis widget, see https://github.com/glvis/xeus-glvis for the widget backend and https://github.com/GLVis/pyglvis/tree/master/js for the widget frontend."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ordinary-equation",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "std::stringstream ss;\n",
+    "ss << \"solution\\n\" << mesh << x << flush;\n",
+    "\n",
+    "auto glv = glvis::glvis();\n",
+    "glv.plot(ss.str() + \"keys Rjml\"); // the `+ \"keys ....\"' is optional\n",
+    "glv"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "C++14",
+   "language": "C++14",
+   "name": "xcpp14"
+  },
+  "language_info": {
+   "codemirror_mode": "text/x-c++src",
+   "file_extension": ".cpp",
+   "mimetype": "text/x-c++src",
+   "name": "c++",
+   "version": "14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

examples/makefile

@@ -22,10 +22,10 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
-  ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29
+  ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
   ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
-  ex25p ex26p ex27p ex28p ex29p
+  ex25p ex26p ex27p ex28p ex29p ex30p
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p \
   ex24p ex25p ex26p

examples/petsc/ex3p.cpp

@@ -121,9 +121,7 @@ int main(int argc, char *argv[])
 
    // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
    //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
-   //    meshes need to be reoriented before we can define high-order Nedelec
-   //    spaces on them.
+   //    parallel mesh is defined, the serial mesh can be deleted.
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
    {
@@ -133,7 +131,6 @@ int main(int argc, char *argv[])
          pmesh->UniformRefinement();
       }
    }
-   pmesh->ReorientTetMesh();
 
    // 6. Define a parallel finite element space on the parallel mesh. Here we
    //    use the Nedelec finite elements of the specified order.

examples/petsc/ex4p.cpp

@@ -122,9 +122,7 @@ int main(int argc, char *argv[])
 
    // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
    //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
-   //    meshes need to be reoriented before we can define high-order Nedelec
-   //    spaces on them (this is needed in the ADS solver below).
+   //    parallel mesh is defined, the serial mesh can be deleted.
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
    {
@@ -134,7 +132,6 @@ int main(int argc, char *argv[])
          pmesh->UniformRefinement();
       }
    }
-   pmesh->ReorientTetMesh();
 
    // 6. Define a parallel finite element space on the parallel mesh. Here we
    //    use the Raviart-Thomas finite elements of the specified order.

fem/CMakeLists.txt

@@ -39,9 +39,19 @@ set(SRCS
   complex_fem.cpp
   convergence.cpp
   datacollection.cpp
+  doftrans.cpp
   eltrans.cpp
   estimators.cpp
   fe.cpp
+  fe/fe_base.cpp
+  fe/fe_fixed_order.cpp
+  fe/fe_h1.cpp
+  fe/fe_l2.cpp
+  fe/fe_nd.cpp
+  fe/fe_nurbs.cpp
+  fe/fe_pos.cpp
+  fe/fe_rt.cpp
+  fe/fe_ser.cpp
   fe_coll.cpp
   fespace.cpp
   geom.cpp
@@ -105,6 +115,7 @@ set(SRCS
   tmop/tmop_pa_w3.cpp
   tmop/tmop_pa_w3_c0.cpp
   tmop_tools.cpp
+  tmop_amr.cpp
   gslib.cpp
   transfer.cpp
   lor.cpp
@@ -118,9 +129,19 @@ set(HDRS
   complex_fem.hpp
   convergence.hpp
   datacollection.hpp
+  doftrans.hpp
   eltrans.hpp
   estimators.hpp
   fe.hpp
+  fe/fe_base.hpp
+  fe/fe_fixed_order.hpp
+  fe/fe_h1.hpp
+  fe/fe_l2.hpp
+  fe/fe_nd.hpp
+  fe/fe_nurbs.hpp
+  fe/fe_pos.hpp
+  fe/fe_rt.hpp
+  fe/fe_ser.hpp
   fe_coll.hpp
   fem.hpp
   fespace.hpp
@@ -164,6 +185,7 @@ set(HDRS
   tmop.hpp
   tmop/tmop_pa.hpp
   tmop_tools.hpp
+  tmop_amr.hpp
   gslib.hpp
   transfer.hpp
   lor.hpp

fem/bilinearform.cpp

@@ -391,6 +391,7 @@ void BilinearForm::Assemble(int skip_zeros)
    }
 
    ElementTransformation *eltrans;
+   DofTransformation * doftrans;
    Mesh *mesh = fes -> GetMesh();
    DenseMatrix elmat, *elmat_p;
 
@@ -424,7 +425,7 @@ void BilinearForm::Assemble(int skip_zeros)
       for (int i = 0; i < fes -> GetNE(); i++)
       {
          int elem_attr = fes->GetMesh()->GetAttribute(i);
-         fes->GetElementVDofs(i, vdofs);
+         doftrans = fes->GetElementVDofs(i, vdofs);
          if (element_matrices)
          {
             elmat_p = &(*element_matrices)(i);
@@ -458,6 +459,11 @@ void BilinearForm::Assemble(int skip_zeros)
             {
                elmat_p = &elmat;
             }
+            if (doftrans)
+            {
+               doftrans->TransformDual(elmat);
+            }
+            elmat_p = &elmat;
          }
          if (static_cond)
          {
@@ -503,7 +509,7 @@ void BilinearForm::Assemble(int skip_zeros)
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
 
          const FiniteElement &be = *fes->GetBE(i);
-         fes -> GetBdrElementVDofs (i, vdofs);
+         doftrans = fes -> GetBdrElementVDofs (i, vdofs);
          eltrans = fes -> GetBdrElementTransformation (i);
          int k = 0;
          for (; k < boundary_integs.Size(); k++)
@@ -523,17 +529,22 @@ void BilinearForm::Assemble(int skip_zeros)
             boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
             elmat += elemmat;
          }
+         if (doftrans)
+         {
+            doftrans->TransformDual(elmat);
+         }
+         elmat_p = &elmat;
          if (!static_cond)
          {
-            mat->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
+            mat->AddSubMatrix(vdofs, vdofs, *elmat_p, skip_zeros);
             if (hybridization)
             {
-               hybridization->AssembleBdrMatrix(i, elmat);
+               hybridization->AssembleBdrMatrix(i, *elmat_p);
             }
          }
          else
          {
-            static_cond->AssembleBdrMatrix(i, elmat);
+            static_cond->AssembleBdrMatrix(i, *elmat_p);
          }
       }
    }
@@ -958,6 +969,7 @@ void BilinearForm::EliminateVDofs(const Array<int> &vdofs,
                                   const Vector &sol, Vector &rhs,
                                   DiagonalPolicy dpolicy)
 {
+   vdofs.HostRead();
    for (int i = 0; i < vdofs.Size(); i++)
    {
       int vdof = vdofs[i];
@@ -1318,9 +1330,10 @@ void MixedBilinearForm::Assemble (int skip_zeros)
       return;
    }
 
-   Array<int> tr_vdofs, te_vdofs;
    ElementTransformation *eltrans;
-   DenseMatrix elemmat;
+   DofTransformation * dom_dof_trans;
+   DofTransformation * ran_dof_trans;
+   DenseMatrix elmat;
 
    Mesh *mesh = test_fes -> GetMesh();
 
@@ -1333,16 +1346,24 @@ void MixedBilinearForm::Assemble (int skip_zeros)
    {
       for (int i = 0; i < test_fes -> GetNE(); i++)
       {
-         trial_fes -> GetElementVDofs (i, tr_vdofs);
-         test_fes  -> GetElementVDofs (i, te_vdofs);
+         dom_dof_trans = trial_fes -> GetElementVDofs (i, trial_vdofs);
+         ran_dof_trans = test_fes  -> GetElementVDofs (i, test_vdofs);
          eltrans = test_fes -> GetElementTransformation (i);
+
+         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
+         elmat = 0.0;
          for (int k = 0; k < domain_integs.Size(); k++)
          {
             domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
                                                         *test_fes  -> GetFE(i),
                                                         *eltrans, elemmat);
-            mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
+            elmat += elemmat;
          }
+         if (ran_dof_trans || dom_dof_trans)
+         {
+            TransformDual(ran_dof_trans, dom_dof_trans, elmat);
+         }
+         mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
       }
    }
 
@@ -1374,9 +1395,12 @@ void MixedBilinearForm::Assemble (int skip_zeros)
          const int bdr_attr = mesh->GetBdrAttribute(i);
          if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
 
-         trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
-         test_fes  -> GetBdrElementVDofs (i, te_vdofs);
+         dom_dof_trans = trial_fes -> GetBdrElementVDofs (i, trial_vdofs);
+         ran_dof_trans = test_fes  -> GetBdrElementVDofs (i, test_vdofs);
          eltrans = test_fes -> GetBdrElementTransformation (i);
+
+         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
+         elmat = 0.0;
          for (int k = 0; k < boundary_integs.Size(); k++)
          {
             if (boundary_integs_marker[k] &&
@@ -1385,29 +1409,34 @@ void MixedBilinearForm::Assemble (int skip_zeros)
             boundary_integs[k]->AssembleElementMatrix2 (*trial_fes -> GetBE(i),
                                                         *test_fes  -> GetBE(i),
                                                         *eltrans, elemmat);
-            mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
+            elmat += elemmat;
          }
+         if (ran_dof_trans || dom_dof_trans)
+         {
+            TransformDual(ran_dof_trans, dom_dof_trans, elmat);
+         }
+         mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
       }
    }
 
    if (trace_face_integs.Size())
    {
       FaceElementTransformations *ftr;
-      Array<int> te_vdofs2;
+      Array<int> test_vdofs2;
       const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
 
       int nfaces = mesh->GetNumFaces();
       for (int i = 0; i < nfaces; i++)
       {
          ftr = mesh->GetFaceElementTransformations(i);
-         trial_fes->GetFaceVDofs(i, tr_vdofs);
-         test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
+         trial_fes->GetFaceVDofs(i, trial_vdofs);
+         test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
          trial_face_fe = trial_fes->GetFaceElement(i);
          test_fe1 = test_fes->GetFE(ftr->Elem1No);
          if (ftr->Elem2No >= 0)
          {
-            test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2);
-            te_vdofs.Append(te_vdofs2);
+            test_fes->GetElementVDofs(ftr->Elem2No, test_vdofs2);
+            test_vdofs.Append(test_vdofs2);
             test_fe2 = test_fes->GetFE(ftr->Elem2No);
          }
          else
@@ -1421,7 +1450,7 @@ void MixedBilinearForm::Assemble (int skip_zeros)
          {
             trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1,
                                                      *test_fe2, *ftr, elemmat);
-            mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros);
+            mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
          }
       }
    }
@@ -1461,8 +1490,8 @@ void MixedBilinearForm::Assemble (int skip_zeros)
          ftr = mesh->GetBdrFaceTransformations(i);
          if (ftr)
          {
-            trial_fes->GetFaceVDofs(ftr->ElementNo, tr_vdofs);
-            test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs);
+            trial_fes->GetFaceVDofs(ftr->ElementNo, trial_vdofs);
+            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
             trial_face_fe = trial_fes->GetFaceElement(ftr->ElementNo);
             test_fe1 = test_fes->GetFE(ftr->Elem1No);
             // The test_fe2 object is really a dummy and not used on the
@@ -1479,7 +1508,7 @@ void MixedBilinearForm::Assemble (int skip_zeros)
                                                                  *test_fe1,
                                                                  *test_fe2,
                                                                  *ftr, elemmat);
-               mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros);
+               mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
             }
          }
       }
@@ -1841,6 +1870,8 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
 
    Array<int> dom_vdofs, ran_vdofs;
    ElementTransformation *T;
+   DofTransformation * dom_dof_trans;
+   DofTransformation * ran_dof_trans;
    const FiniteElement *dom_fe, *ran_fe;
    DenseMatrix totelmat, elmat;
 
@@ -1853,8 +1884,8 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
    {
       for (int i = 0; i < test_fes->GetNE(); i++)
       {
-         trial_fes->GetElementVDofs(i, dom_vdofs);
-         test_fes->GetElementVDofs(i, ran_vdofs);
+         dom_dof_trans = trial_fes->GetElementVDofs(i, dom_vdofs);
+         ran_dof_trans = test_fes->GetElementVDofs(i, ran_vdofs);
          T = test_fes->GetElementTransformation(i);
          dom_fe = trial_fes->GetFE(i);
          ran_fe = test_fes->GetFE(i);
@@ -1867,6 +1898,10 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
                                                      elmat);
             totelmat += elmat;
          }
+         if (ran_dof_trans || dom_dof_trans)
+         {
+            TransformPrimal(ran_dof_trans, dom_dof_trans, totelmat);
+         }
          mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
       }
    }

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);
@@ -984,6 +985,8 @@ void DiffusionIntegrator::ComputeElementFlux
                   "Unexpected height for MatrixCoefficient");
    }
 
+   MFEM_VERIFY(!SMQ, "SymmetricMatrixCoefficient not supported here");
+
 #ifdef MFEM_THREAD_SAFE
    DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim);
    DenseMatrix M(MQ ? spaceDim : 0);
@@ -996,7 +999,7 @@ void DiffusionIntegrator::ComputeElementFlux
 #endif
    vec.SetSize(dim);
    vecdxt.SetSize(spaceDim);
-   pointflux.SetSize(MQ ? spaceDim : 0);
+   pointflux.SetSize(MQ || VQ ? spaceDim : 0);
 
    const IntegrationRule &ir = fluxelem.GetNodes();
    fnd = ir.GetNPoints();
@@ -1012,36 +1015,45 @@ void DiffusionIntegrator::ComputeElementFlux
       CalcInverse(Trans.Jacobian(), invdfdx);
       invdfdx.MultTranspose(vec, vecdxt);
 
-      if (!MQ && !VQ)
+      if (with_coef)
       {
-         if (Q && with_coef)
+         if (!MQ && !VQ)
          {
-            vecdxt *= Q->Eval(Trans,ip);
+            if (Q)
+            {
+               vecdxt *= Q->Eval(Trans,ip);
+            }
+            for (j = 0; j < spaceDim; j++)
+            {
+               flux(fnd*j+i) = vecdxt(j);
+            }
          }
-         for (j = 0; j < spaceDim; j++)
+         else
          {
-            flux(fnd*j+i) = vecdxt(j);
+            if (MQ)
+            {
+               MQ->Eval(M, Trans, ip);
+               M.Mult(vecdxt, pointflux);
+            }
+            else
+            {
+               VQ->Eval(D, Trans, ip);
+               for (int j=0; j<spaceDim; ++j)
+               {
+                  pointflux[j] = D[j] * vecdxt[j];
+               }
+            }
+            for (j = 0; j < spaceDim; j++)
+            {
+               flux(fnd*j+i) = pointflux(j);
+            }
          }
       }
       else
       {
-         if (MQ)
-         {
-            MQ->Eval(M, Trans, ip);
-            M.Mult(vecdxt, pointflux);
-         }
-         else
-         {
-            VQ->Eval(D, Trans, ip);
-            for (int j=0; j<spaceDim; ++j)
-            {
-               pointflux[j] = D[j] * vecdxt[j];
-            }
-
-         }
          for (j = 0; j < spaceDim; j++)
          {
-            flux(fnd*j+i) = pointflux(j);
+            flux(fnd*j+i) = vecdxt(j);
          }
       }
    }
@@ -1057,8 +1069,13 @@ double DiffusionIntegrator::ComputeFluxEnergy
 
 #ifdef MFEM_THREAD_SAFE
    DenseMatrix M;
+   Vector D(VQ ? VQ->GetVDim() : 0);
+#else
+   D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
 
+   MFEM_VERIFY(!SMQ, "SymmetricMatrixCoefficient not supported here");
+
    shape.SetSize(nd);
    pointflux.SetSize(spaceDim);
    if (d_energy) { vec.SetSize(spaceDim); }
@@ -1087,17 +1104,23 @@ double DiffusionIntegrator::ComputeFluxEnergy
       Trans.SetIntPoint(&ip);
       double w = Trans.Weight() * ip.weight;
 
-      if (!MQ)
+      if (MQ)
+      {
+         MQ->Eval(M, Trans, ip);
+         energy += w * M.InnerProduct(pointflux, pointflux);
+      }
+      else if (VQ)
+      {
+         VQ->Eval(D, Trans, ip);
+         D *= pointflux;
+         energy += w * (D * pointflux);
+      }
+      else
       {
          double e = (pointflux * pointflux);
          if (Q) { e *= Q->Eval(Trans, ip); }
          energy += w * e;
       }
-      else
-      {
-         MQ->Eval(M, Trans, ip);
-         energy += w * M.InnerProduct(pointflux, pointflux);
-      }
 
       if (d_energy)
       {
@@ -1107,7 +1130,7 @@ double DiffusionIntegrator::ComputeFluxEnergy
          {
             (*d_energy)[k] += w * vec[k] * vec[k];
          }
-         // TODO: Q, MQ
+         // TODO: Q, VQ, MQ
       }
    }
 

fem/bilininteg.hpp

@@ -1982,24 +1982,32 @@ private:
 
 public:
    /// Construct a diffusion integrator with coefficient Q = 1
-   DiffusionIntegrator()
-      : Q(NULL), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(const IntegrationRule *ir = nullptr)
+      : BilinearFormIntegrator(ir),
+        Q(NULL), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
 
    /// Construct a diffusion integrator with a scalar coefficient q
-   DiffusionIntegrator(Coefficient &q)
-      : Q(&q), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(Coefficient &q, const IntegrationRule *ir = nullptr)
+      : BilinearFormIntegrator(ir),
+        Q(&q), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
 
    /// Construct a diffusion integrator with a vector coefficient q
-   DiffusionIntegrator(VectorCoefficient &q)
-      : Q(NULL), VQ(&q), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(VectorCoefficient &q,
+                       const IntegrationRule *ir = nullptr)
+      : BilinearFormIntegrator(ir),
+        Q(NULL), VQ(&q), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { }
 
    /// Construct a diffusion integrator with a matrix coefficient q
-   DiffusionIntegrator(MatrixCoefficient &q)
-      : Q(NULL), VQ(NULL), MQ(&q), SMQ(NULL), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(MatrixCoefficient &q,
+                       const IntegrationRule *ir = nullptr)
+      : BilinearFormIntegrator(ir),
+        Q(NULL), VQ(NULL), MQ(&q), SMQ(NULL), maps(NULL), geom(NULL) { }
 
    /// Construct a diffusion integrator with a symmetric matrix coefficient q
-   DiffusionIntegrator(SymmetricMatrixCoefficient &q)
-      : Q(NULL), VQ(NULL), MQ(NULL), SMQ(&q), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(SymmetricMatrixCoefficient &q,
+                       const IntegrationRule *ir = nullptr)
+      : BilinearFormIntegrator(ir),
+        Q(NULL), VQ(NULL), MQ(NULL), SMQ(&q), maps(NULL), geom(NULL) { }
 
    /** Given a particular Finite Element computes the element stiffness matrix
        elmat. */
@@ -2668,6 +2676,9 @@ public:
    VectorDiffusionIntegrator(Coefficient &q)
       : Q(&q) { }
 
+   VectorDiffusionIntegrator(Coefficient &q, const IntegrationRule *ir)
+      : BilinearFormIntegrator(ir), Q(&q) { }
+
    /** \brief Integrator with scalar coefficient for caller-specified vector
        dimension.
 
@@ -2927,10 +2938,11 @@ public:
 
     sum_e eta (r_e([u]), r_e([v]))
 
-    where r_e is the lifting operator defined on each edge e. The parameter eta
-    can be chosen to be one to obtain a stable discretization. The constructor
-    for this integrator requires the finite element space because the lifting
-    operator depends on the element-wise inverse mass matrix.
+    where r_e is the lifting operator defined on each edge e (potentially
+    weighted by a coefficient Q). The parameter eta can be chosen to be one to
+    obtain a stable discretization. The constructor for this integrator requires
+    the finite element space because the lifting operator depends on the
+    element-wise inverse mass matrix.
 
     BR2 stands for the second method of Bassi and Rebay:
 
@@ -2953,14 +2965,28 @@ protected:
    Array<int> ipiv;
    Array<int> ipiv_offsets, Minv_offsets;
 
+   Coefficient *Q;
+
    Vector shape1, shape2;
 
    DenseMatrix R11, R12, R21, R22;
    DenseMatrix MinvR11, MinvR12, MinvR21, MinvR22;
    DenseMatrix Re, MinvRe;
 
+   /// Precomputes the inverses (LU factorizations) of the local mass matrices.
+   /** @a fes must be a DG space, so the mass matrix is block diagonal, and its
+       inverse can be computed locally. This is required for the computation of
+       the lifting operators @a r_e.
+   */
+   void PrecomputeMassInverse(class FiniteElementSpace &fes);
+
 public:
-   DGDiffusionBR2Integrator(class FiniteElementSpace *fes, double e = 1.0);
+   DGDiffusionBR2Integrator(class FiniteElementSpace &fes, double e = 1.0);
+   DGDiffusionBR2Integrator(class FiniteElementSpace &fes, Coefficient &Q_,
+                            double e = 1.0);
+   MFEM_DEPRECATED DGDiffusionBR2Integrator(class FiniteElementSpace *fes,
+                                            double e = 1.0);
+
    using BilinearFormIntegrator::AssembleFaceMatrix;
    virtual void AssembleFaceMatrix(const FiniteElement &el1,
                                    const FiniteElement &el2,

fem/bilininteg_br2.cpp

@@ -16,20 +16,39 @@
 namespace mfem
 {
 
-DGDiffusionBR2Integrator::DGDiffusionBR2Integrator(FiniteElementSpace *fes,
-                                                   double e) : eta(e)
+DGDiffusionBR2Integrator::DGDiffusionBR2Integrator(
+   FiniteElementSpace &fes, double e) : eta(e), Q(NULL)
 {
+   PrecomputeMassInverse(fes);
+}
+
+DGDiffusionBR2Integrator::DGDiffusionBR2Integrator(
+   FiniteElementSpace &fes, Coefficient &Q_, double e) : eta(e), Q(&Q_)
+{
+   PrecomputeMassInverse(fes);
+}
+
+DGDiffusionBR2Integrator::DGDiffusionBR2Integrator(
+   FiniteElementSpace *fes, double e) : eta(e), Q(NULL)
+{
+   PrecomputeMassInverse(*fes);
+}
+
+void DGDiffusionBR2Integrator::PrecomputeMassInverse(FiniteElementSpace &fes)
+{
+   MFEM_VERIFY(fes.IsDGSpace(),
+               "The BR2 integrator is only defined for DG spaces.");
    // Precompute local mass matrix inverses needed for the lifting operators
    // First compute offsets and total size needed (e.g. for mixed meshes or
    // p-refinement)
-   int nel = fes->GetNE();
+   int nel = fes.GetNE();
    Minv_offsets.SetSize(nel+1);
    ipiv_offsets.SetSize(nel+1);
    ipiv_offsets[0] = 0;
    Minv_offsets[0] = 0;
    for (int i=0; i<nel; ++i)
    {
-      int dof = fes->GetFE(i)->GetDof();
+      int dof = fes.GetFE(i)->GetDof();
       ipiv_offsets[i+1] = ipiv_offsets[i] + dof;
       Minv_offsets[i+1] = Minv_offsets[i] + dof*dof;
    }
@@ -37,7 +56,7 @@ DGDiffusionBR2Integrator::DGDiffusionBR2Integrator(FiniteElementSpace *fes,
 #ifdef MFEM_USE_MPI
    // When running in parallel, we also need to compute the local mass matrices
    // of face neighbor elements
-   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace *>(fes);
+   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace *>(&fes);
    if (pfes != NULL)
    {
       ParMesh *pmesh = pfes->GetParMesh();
@@ -64,15 +83,15 @@ DGDiffusionBR2Integrator::DGDiffusionBR2Integrator(FiniteElementSpace *fes,
    {
       const FiniteElement *fe = NULL;
       ElementTransformation *tr = NULL;
-      if (i < fes->GetNE())
+      if (i < fes.GetNE())
       {
-         fe = fes->GetFE(i);
-         tr = fes->GetElementTransformation(i);
+         fe = fes.GetFE(i);
+         tr = fes.GetElementTransformation(i);
       }
       else
       {
 #ifdef MFEM_USE_MPI
-         int inbr = i - fes->GetNE();
+         int inbr = i - fes.GetNE();
          fe = pfes->GetFaceNbrFE(inbr);
          tr = pfes->GetParMesh()->GetFaceNbrElementTransformation(inbr);
 #endif
@@ -151,21 +170,24 @@ void DGDiffusionBR2Integrator::AssembleFaceMatrix(
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);
-      IntegrationPoint eip1, eip2;
+      Trans.SetAllIntPoints(&ip);
 
-      Trans.Loc1.Transform(ip, eip1);
+      const IntegrationPoint &eip1 = Trans.Elem1->GetIntPoint();
       el1.CalcShape(eip1, shape1);
+      double q = Q ? Q->Eval(*Trans.Elem1, eip1) : 1.0;
       if (ndof2)
       {
-         Trans.Loc2.Transform(ip, eip2);
+         const IntegrationPoint &eip2 = Trans.Elem2->GetIntPoint();
          el2.CalcShape(eip2, shape2);
+         // Set coefficient value q to the average of the values on either side
+         if (Q) { q = 0.5*(q + Q->Eval(*Trans.Elem2, eip2)); }
       }
-
-      double w = factor*sqrt(eta)*ip.weight*Trans.Face->Weight();
-      if (ndof2)
-      {
-         w /= 2;
-      }
+      // Take sqrt here because
+      //    eta (r_e([u]), r_e([v])) = (sqrt(eta) r_e([u]), sqrt(eta) r_e([v]))
+      double w = sqrt((factor + 1)*eta*q)*ip.weight*Trans.Face->Weight();
+      // r_e is defined by, (r_e([u]), tau) = <[u], {tau}>, so we pick up a
+      // factor of 0.5 on interior faces from the average term.
+      if (ndof2) { w *= 0.5; }
 
       for (int i = 0; i < ndof1; i++)
       {

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
@@ -903,9 +903,11 @@ static void PADiffusionAssembleDiagonal(const int dim,
    {
       switch ((D1D << 4 ) | Q1D)
       {
+         case 0x22: return SmemPADiffusionDiagonal3D<2,2>(NE,symm,B,G,D,Y);
          case 0x23: return SmemPADiffusionDiagonal3D<2,3>(NE,symm,B,G,D,Y);
          case 0x34: return SmemPADiffusionDiagonal3D<3,4>(NE,symm,B,G,D,Y);
          case 0x45: return SmemPADiffusionDiagonal3D<4,5>(NE,symm,B,G,D,Y);
+         case 0x46: return SmemPADiffusionDiagonal3D<4,6>(NE,symm,B,G,D,Y);
          case 0x56: return SmemPADiffusionDiagonal3D<5,6>(NE,symm,B,G,D,Y);
          case 0x67: return SmemPADiffusionDiagonal3D<6,7>(NE,symm,B,G,D,Y);
          case 0x78: return SmemPADiffusionDiagonal3D<7,8>(NE,symm,B,G,D,Y);
@@ -1877,6 +1879,7 @@ static void PADiffusionApply(const int dim,
    {
       switch (ID)
       {
+         case 0x22: return SmemPADiffusionApply3D<2,2>(NE,symm,B,G,D,X,Y);
          case 0x23: return SmemPADiffusionApply3D<2,3>(NE,symm,B,G,D,X,Y);
          case 0x34: return SmemPADiffusionApply3D<3,4>(NE,symm,B,G,D,X,Y);
          case 0x45: return SmemPADiffusionApply3D<4,5>(NE,symm,B,G,D,X,Y);

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/bilininteg_mass_pa.cpp

@@ -1203,8 +1203,10 @@ static void PAMassApply(const int dim,
    {
       switch (id)
       {
+         case 0x22: return SmemPAMassApply3D<2,2>(NE,B,Bt,D,X,Y);
          case 0x23: return SmemPAMassApply3D<2,3>(NE,B,Bt,D,X,Y);
          case 0x24: return SmemPAMassApply3D<2,4>(NE,B,Bt,D,X,Y);
+         case 0x26: return SmemPAMassApply3D<2,6>(NE,B,Bt,D,X,Y);
          case 0x34: return SmemPAMassApply3D<3,4>(NE,B,Bt,D,X,Y);
          case 0x35: return SmemPAMassApply3D<3,5>(NE,B,Bt,D,X,Y);
          case 0x36: return SmemPAMassApply3D<3,6>(NE,B,Bt,D,X,Y);

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/ceed/util.cpp

@@ -186,14 +186,18 @@ static void InitTensorBasis(const mfem::FiniteElementSpace &fes,
    const int ndofs = maps.ndof;
    const int nqpts = maps.nqpt;
    mfem::Vector qX(nqpts), qW(nqpts);
-   const mfem::IntegrationRule &ir1d =
-      IntRules.Get(Geometry::SEGMENT, ir.GetOrder());
+   // The x-coordinates of the first `nqpts` points of the integration rule are
+   // the points of the corresponding 1D rule. We also scale the weights
+   // accordingly.
+   double w_sum = 0.0;
    for (int i = 0; i < nqpts; i++)
    {
-      const mfem::IntegrationPoint &ip = ir1d.IntPoint(i);
+      const mfem::IntegrationPoint &ip = ir.IntPoint(i);
       qX(i) = ip.x;
       qW(i) = ip.weight;
+      w_sum += ip.weight;
    }
+   qW *= 1.0/w_sum;
    CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes.GetVDim(), ndofs,
                            nqpts, maps.Bt.GetData(),
                            maps.Gt.GetData(), qX.GetData(),

fem/coefficient.cpp

@@ -52,6 +52,13 @@ double GridFunctionCoefficient::Eval (ElementTransformation &T,
    return GridF -> GetValue (T, ip, Component);
 }
 
+void TransformedCoefficient::SetTime(double t)
+{
+   if (Q1) { Q1->SetTime(t); }
+   if (Q2) { Q2->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
 double TransformedCoefficient::Eval(ElementTransformation &T,
                                     const IntegrationPoint &ip)
 {
@@ -66,6 +73,12 @@ double TransformedCoefficient::Eval(ElementTransformation &T,
    }
 }
 
+void DeltaCoefficient::SetTime(double t)
+{
+   if (weight) { weight->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
 void DeltaCoefficient::SetDeltaCenter(const Vector& vcenter)
 {
    MFEM_VERIFY(vcenter.Size() <= 3,
@@ -87,6 +100,12 @@ double DeltaCoefficient::EvalDelta(ElementTransformation &T,
    return weight ? weight->Eval(T, ip, GetTime())*w : w;
 }
 
+void RestrictedCoefficient::SetTime(double t)
+{
+   if (c) { c->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
 void VectorCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
                              const IntegrationRule &ir)
 {
@@ -134,6 +153,15 @@ VectorArrayCoefficient::VectorArrayCoefficient (int dim)
    }
 }
 
+void VectorArrayCoefficient::SetTime(double t)
+{
+   for (int i = 0; i < vdim; i++)
+   {
+      if (Coeff[i]) { Coeff[i]->SetTime(t); }
+   }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void VectorArrayCoefficient::Set(int i, Coefficient *c, bool own)
 {
    if (ownCoeff[i]) { delete Coeff[i]; }
@@ -247,6 +275,12 @@ double DivergenceGridFunctionCoefficient::Eval(ElementTransformation &T,
    return GridFunc->GetDivergence(T);
 }
 
+void VectorDeltaCoefficient::SetTime(double t)
+{
+   d.SetTime(t);
+   this->VectorCoefficient::SetTime(t);
+}
+
 void VectorDeltaCoefficient::SetDirection(const Vector &d_)
 {
    dir = d_;
@@ -261,6 +295,12 @@ void VectorDeltaCoefficient::EvalDelta(
    V *= d.EvalDelta(T, ip);
 }
 
+void VectorRestrictedCoefficient::SetTime(double t)
+{
+   if (c) { c->SetTime(t); }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void VectorRestrictedCoefficient::Eval(Vector &V, ElementTransformation &T,
                                        const IntegrationPoint &ip)
 {
@@ -291,6 +331,12 @@ void VectorRestrictedCoefficient::Eval(
    }
 }
 
+void MatrixFunctionCoefficient::SetTime(double t)
+{
+   if (Q) { Q->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void MatrixFunctionCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                      const IntegrationPoint &ip)
 {
@@ -371,6 +417,12 @@ void MatrixFunctionCoefficient::EvalSymmetric(Vector &K,
    }
 }
 
+void SymmetricMatrixFunctionCoefficient::SetTime(double t)
+{
+   if (Q) { Q->SetTime(t); }
+   this->SymmetricMatrixCoefficient::SetTime(t);
+}
+
 void SymmetricMatrixFunctionCoefficient::Eval(DenseSymmetricMatrix &K,
                                               ElementTransformation &T,
                                               const IntegrationPoint &ip)
@@ -413,6 +465,15 @@ MatrixArrayCoefficient::MatrixArrayCoefficient (int dim)
    }
 }
 
+void MatrixArrayCoefficient::SetTime(double t)
+{
+   for (int i=0; i < height*width; i++)
+   {
+      if (Coeff[i]) { Coeff[i]->SetTime(t); }
+   }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void MatrixArrayCoefficient::Set(int i, int j, Coefficient * c, bool own)
 {
    if (ownCoeff[i*width+j]) { delete Coeff[i*width+j]; }
@@ -431,6 +492,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++)
@@ -440,6 +502,12 @@ void MatrixArrayCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
    }
 }
 
+void MatrixRestrictedCoefficient::SetTime(double t)
+{
+   if (c) { c->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void MatrixRestrictedCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                        const IntegrationPoint &ip)
 {
@@ -455,6 +523,33 @@ void MatrixRestrictedCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
    }
 }
 
+void SumCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+void ProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+void RatioCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+void PowerCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
 InnerProductCoefficient::InnerProductCoefficient(VectorCoefficient &A,
                                                  VectorCoefficient &B)
    : a(&A), b(&B)
@@ -464,6 +559,13 @@ InnerProductCoefficient::InnerProductCoefficient(VectorCoefficient &A,
                "Arguments have incompatible dimensions.");
 }
 
+void InnerProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
 double InnerProductCoefficient::Eval(ElementTransformation &T,
                                      const IntegrationPoint &ip)
 {
@@ -481,6 +583,13 @@ VectorRotProductCoefficient::VectorRotProductCoefficient(VectorCoefficient &A,
                "Arguments must have dimension equal to two.");
 }
 
+void VectorRotProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
 double VectorRotProductCoefficient::Eval(ElementTransformation &T,
                                          const IntegrationPoint &ip)
 {
@@ -497,6 +606,12 @@ DeterminantCoefficient::DeterminantCoefficient(MatrixCoefficient &A)
                "Argument must be a square matrix.");
 }
 
+void DeterminantCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
 double DeterminantCoefficient::Eval(ElementTransformation &T,
                                     const IntegrationPoint &ip)
 {
@@ -545,6 +660,15 @@ VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &A_,
                "Arguments must have the same dimension.");
 }
 
+void VectorSumCoefficient::SetTime(double t)
+{
+   if (ACoef) { ACoef->SetTime(t); }
+   if (BCoef) { BCoef->SetTime(t); }
+   if (alphaCoef) { alphaCoef->SetTime(t); }
+   if (betaCoef) { betaCoef->SetTime(t); }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void VectorSumCoefficient::Eval(Vector &V, ElementTransformation &T,
                                 const IntegrationPoint &ip)
 {
@@ -568,6 +692,13 @@ ScalarVectorProductCoefficient::ScalarVectorProductCoefficient(
    : VectorCoefficient(B.GetVDim()), aConst(0.0), a(&A), b(&B)
 {}
 
+void ScalarVectorProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void ScalarVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
                                           const IntegrationPoint &ip)
 {
@@ -581,6 +712,12 @@ NormalizedVectorCoefficient::NormalizedVectorCoefficient(VectorCoefficient &A,
    : VectorCoefficient(A.GetVDim()), a(&A), tol(tol_)
 {}
 
+void NormalizedVectorCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void NormalizedVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
                                        const IntegrationPoint &ip)
 {
@@ -599,6 +736,13 @@ VectorCrossProductCoefficient::VectorCrossProductCoefficient(
                "Arguments must have dimension equal to three.");
 }
 
+void VectorCrossProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void VectorCrossProductCoefficient::Eval(Vector &V, ElementTransformation &T,
                                          const IntegrationPoint &ip)
 {
@@ -620,6 +764,13 @@ MatrixVectorProductCoefficient::MatrixVectorProductCoefficient(
                "Arguments have incompatible dimensions.");
 }
 
+void MatrixVectorProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->VectorCoefficient::SetTime(t);
+}
+
 void MatrixVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
                                           const IntegrationPoint &ip)
 {
@@ -649,6 +800,13 @@ MatrixSumCoefficient::MatrixSumCoefficient(MatrixCoefficient &A,
                "Arguments must have the same dimensions.");
 }
 
+void MatrixSumCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void MatrixSumCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
                                 const IntegrationPoint &ip)
 {
@@ -670,6 +828,13 @@ ScalarMatrixProductCoefficient::ScalarMatrixProductCoefficient(
    : MatrixCoefficient(B.GetHeight(), B.GetWidth()), aConst(0.0), a(&A), b(&B)
 {}
 
+void ScalarMatrixProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void ScalarMatrixProductCoefficient::Eval(DenseMatrix &M,
                                           ElementTransformation &T,
                                           const IntegrationPoint &ip)
@@ -683,6 +848,12 @@ TransposeMatrixCoefficient::TransposeMatrixCoefficient(MatrixCoefficient &A)
    : MatrixCoefficient(A.GetWidth(), A.GetHeight()), a(&A)
 {}
 
+void TransposeMatrixCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void TransposeMatrixCoefficient::Eval(DenseMatrix &M,
                                       ElementTransformation &T,
                                       const IntegrationPoint &ip)
@@ -699,6 +870,12 @@ InverseMatrixCoefficient::InverseMatrixCoefficient(MatrixCoefficient &A)
                "Argument must be a square matrix.");
 }
 
+void InverseMatrixCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void InverseMatrixCoefficient::Eval(DenseMatrix &M,
                                     ElementTransformation &T,
                                     const IntegrationPoint &ip)
@@ -713,6 +890,13 @@ OuterProductCoefficient::OuterProductCoefficient(VectorCoefficient &A,
      va(A.GetVDim()), vb(B.GetVDim())
 {}
 
+void OuterProductCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (b) { b->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void OuterProductCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
                                    const IntegrationPoint &ip)
 {
@@ -739,6 +923,13 @@ CrossCrossCoefficient::CrossCrossCoefficient(Coefficient &A,
      vk(K.GetVDim())
 {}
 
+void CrossCrossCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   if (k) { k->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
 void CrossCrossCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
                                  const IntegrationPoint &ip)
 {

fem/coefficient.hpp

@@ -45,7 +45,7 @@ public:
    Coefficient() { time = 0.; }
 
    /// Set the time for time dependent coefficients
-   void SetTime(double t) { time = t; }
+   virtual void SetTime(double t) { time = t; }
 
    /// Get the time for time dependent coefficients
    double GetTime() { return time; }
@@ -217,6 +217,9 @@ public:
                            double (*F)(double,double))
       : Q1(q1), Q2(q2), Transform2(F) { Transform1 = 0; }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Evaluate the coefficient at @a ip.
    virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
 };
@@ -269,6 +272,9 @@ public:
       weight = NULL; sdim = 3; tdf = NULL;
    }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Set the center location of the delta function.
    void SetDeltaCenter(const Vector& center);
 
@@ -333,6 +339,9 @@ public:
    RestrictedCoefficient(Coefficient &c_, Array<int> &attr)
    { c = &c_; attr.Copy(active_attr); }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Evaluate the coefficient at @a ip.
    virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
    { return active_attr[T.Attribute-1] ? c->Eval(T, ip, GetTime()) : 0.0; }
@@ -350,7 +359,7 @@ public:
    VectorCoefficient(int vd) { vdim = vd; time = 0.; }
 
    /// Set the time for time dependent coefficients
-   void SetTime(double t) { time = t; }
+   virtual void SetTime(double t) { time = t; }
 
    /// Get the time for time dependent coefficients
    double GetTime() { return time; }
@@ -456,6 +465,9 @@ public:
        still need to be added with Set(). */
    explicit VectorArrayCoefficient(int dim);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Returns i'th coefficient.
    Coefficient* GetCoeff(int i) { return Coeff[i]; }
 
@@ -632,6 +644,9 @@ public:
                           double s)
       : VectorCoefficient(dir_.Size()), dir(dir_), d(x,y,z,s) { }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Replace the associated DeltaCoefficient with a new DeltaCoefficient.
    /** The new DeltaCoefficient cannot have a specified weight Coefficient, i.e.
        DeltaCoefficient::Weight() should return NULL. */
@@ -677,6 +692,9 @@ public:
       : VectorCoefficient(vc.GetVDim())
    { c = &vc; attr.Copy(active_attr); }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Evaluate the vector coefficient at @a ip.
    virtual void Eval(Vector &V, ElementTransformation &T,
                      const IntegrationPoint &ip);
@@ -708,7 +726,7 @@ public:
       height(h), width(w), time(0.), symmetric(symm) { }
 
    /// Set the time for time dependent coefficients
-   void SetTime(double t) { time = t; }
+   virtual void SetTime(double t) { time = t; }
 
    /// Get the time for time dependent coefficients
    double GetTime() { return time; }
@@ -817,6 +835,9 @@ public:
       : MatrixCoefficient(dim), TDFunction(std::move(TDF)), Q(q)
    { }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Evaluate the matrix coefficient at @a ip.
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
                      const IntegrationPoint &ip);
@@ -844,6 +865,9 @@ public:
        actual coefficients still need to be added with Set(). */
    explicit MatrixArrayCoefficient (int dim);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Get the coefficient located at (i,j) in the matrix.
    Coefficient* GetCoeff (int i, int j) { return Coeff[i*width+j]; }
 
@@ -881,6 +905,9 @@ public:
       : MatrixCoefficient(mc.GetHeight(), mc.GetWidth())
    { c = &mc; attr.Copy(active_attr); }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Evaluate the matrix coefficient at @a ip.
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
                      const IntegrationPoint &ip);
@@ -911,6 +938,9 @@ public:
                   double alpha_ = 1.0, double beta_ = 1.0)
       : aConst(0.0), a(&A), b(&B), alpha(alpha_), beta(beta_) { }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first term in the linear combination as a constant
    void SetAConst(double A) { a = NULL; aConst = A; }
    /// Return the first term in the linear combination
@@ -959,7 +989,7 @@ public:
    { dim = dimension; time = 0.; }
 
    /// Set the time for time dependent coefficients
-   void SetTime(double t) { time = t; }
+   virtual void SetTime(double t) { time = t; }
 
    /// Get the time for time dependent coefficients
    double GetTime() { return time; }
@@ -1037,6 +1067,9 @@ public:
       : SymmetricMatrixCoefficient(dim), TDFunction(std::move(TDF)), Q(q)
    { }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Evaluate the matrix coefficient at @a ip.
    virtual void Eval(DenseSymmetricMatrix &K, ElementTransformation &T,
                      const IntegrationPoint &ip);
@@ -1063,6 +1096,9 @@ public:
    ProductCoefficient(Coefficient &A, Coefficient &B)
       : aConst(0.0), a(&A), b(&B) { }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first term in the product as a constant
    void SetAConst(double A) { a = NULL; aConst = A; }
    /// Return the first term in the product
@@ -1108,6 +1144,9 @@ public:
    RatioCoefficient(Coefficient &A, double B)
       : aConst(0.0), bConst(B), a(&A), b(NULL) { }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the numerator in the ratio as a constant
    void SetAConst(double A) { a = NULL; aConst = A; }
    /// Return the numerator of the ratio
@@ -1151,6 +1190,9 @@ public:
    PowerCoefficient(Coefficient &A, double p_)
       : a(&A), p(p_) { }
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the base coefficient
    void SetACoef(Coefficient &A) { a = &A; }
    /// Return the base coefficient
@@ -1181,6 +1223,9 @@ public:
    /// Construct with the two vector coefficients.  Result is \f$ A \cdot B \f$.
    InnerProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first vector in the inner product
    void SetACoef(VectorCoefficient &A) { a = &A; }
    /// Return the first vector coefficient in the inner product
@@ -1210,6 +1255,9 @@ public:
    /// Constructor with two vector coefficients.  Result is \f$ A_x B_y - A_y * B_x; \f$.
    VectorRotProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first vector in the product
    void SetACoef(VectorCoefficient &A) { a = &A; }
    /// Return the first vector of the product
@@ -1237,6 +1285,9 @@ public:
    /// Construct with the matrix.
    DeterminantCoefficient(MatrixCoefficient &A);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the matrix coefficient
    void SetACoef(MatrixCoefficient &A) { a = &A; }
    /// Return the matrix coefficient
@@ -1280,6 +1331,9 @@ public:
    VectorSumCoefficient(VectorCoefficient &A_, VectorCoefficient &B_,
                         Coefficient &alpha_, Coefficient &beta_);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first vector coefficient
    void SetACoef(VectorCoefficient &A) { ACoef = &A; }
    /// Return the first vector coefficient
@@ -1341,6 +1395,9 @@ public:
    /// Constructor with two coefficients.  Result is A * B.
    ScalarVectorProductCoefficient(Coefficient &A, VectorCoefficient &B);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the scalar factor as a constant
    void SetAConst(double A) { a = NULL; aConst = A; }
    /// Return the scalar factor
@@ -1379,6 +1436,9 @@ public:
    */
    NormalizedVectorCoefficient(VectorCoefficient &A, double tol = 1e-6);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the vector coefficient
    void SetACoef(VectorCoefficient &A) { a = &A; }
    /// Return the vector coefficient
@@ -1404,6 +1464,9 @@ public:
    /// Construct with the two coefficients.  Result is A x B.
    VectorCrossProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first term in the product
    void SetACoef(VectorCoefficient &A) { a = &A; }
    /// Return the first term in the product
@@ -1435,6 +1498,9 @@ public:
    /// Constructor with two coefficients.  Result is A*B.
    MatrixVectorProductCoefficient(MatrixCoefficient &A, VectorCoefficient &B);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the matrix coefficient
    void SetACoef(MatrixCoefficient &A) { a = &A; }
    /// Return the matrix coefficient
@@ -1487,6 +1553,9 @@ public:
    MatrixSumCoefficient(MatrixCoefficient &A, MatrixCoefficient &B,
                         double alpha_ = 1.0, double beta_ = 1.0);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first matrix coefficient
    void SetACoef(MatrixCoefficient &A) { a = &A; }
    /// Return the first matrix coefficient
@@ -1528,6 +1597,9 @@ public:
    /// Constructor with two coefficients.  Result is A*B.
    ScalarMatrixProductCoefficient(Coefficient &A, MatrixCoefficient &B);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the scalar factor as a constant
    void SetAConst(double A) { a = NULL; aConst = A; }
    /// Return the scalar factor
@@ -1558,6 +1630,9 @@ public:
    /// Construct with the matrix coefficient.  Result is \f$ A^T \f$.
    TransposeMatrixCoefficient(MatrixCoefficient &A);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the matrix coefficient
    void SetACoef(MatrixCoefficient &A) { a = &A; }
    /// Return the matrix coefficient
@@ -1578,6 +1653,9 @@ public:
    /// Construct with the matrix coefficient.  Result is \f$ A^{-1} \f$.
    InverseMatrixCoefficient(MatrixCoefficient &A);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the matrix coefficient
    void SetACoef(MatrixCoefficient &A) { a = &A; }
    /// Return the matrix coefficient
@@ -1602,6 +1680,9 @@ public:
    /// Construct with two vector coefficients.  Result is \f$ A B^T \f$.
    OuterProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the first vector in the outer product
    void SetACoef(VectorCoefficient &A) { a = &A; }
    /// Return the first vector coefficient in the outer product
@@ -1637,6 +1718,9 @@ public:
    CrossCrossCoefficient(double A, VectorCoefficient &K);
    CrossCrossCoefficient(Coefficient &A, VectorCoefficient &K);
 
+   /// Set the time for internally stored coefficients
+   void SetTime(double t);
+
    /// Reset the scalar factor as a constant
    void SetAConst(double A) { a = NULL; aConst = A; }
    /// Return the scalar factor

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/doftrans.cpp

@@ -0,0 +1,358 @@
+// 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.
+
+#include "fem.hpp"
+
+namespace mfem
+{
+
+void DofTransformation::TransformPrimal(Vector &v) const
+{
+   TransformPrimal(v.GetData());
+}
+
+void DofTransformation::TransformPrimalCols(DenseMatrix &V) const
+{
+   for (int c=0; c<V.Width(); c++)
+   {
+      TransformPrimal(V.GetColumn(c));
+   }
+}
+
+void DofTransformation::TransformDual(Vector &v) const
+{
+   TransformDual(v.GetData());
+}
+
+void DofTransformation::TransformDual(DenseMatrix &V) const
+{
+   TransformDualCols(V);
+   TransformDualRows(V);
+}
+
+void DofTransformation::TransformDualRows(DenseMatrix &V) const
+{
+   Vector row;
+   for (int r=0; r<V.Height(); r++)
+   {
+      V.GetRow(r, row);
+      TransformDual(row);
+      V.SetRow(r, row);
+   }
+}
+
+void DofTransformation::TransformDualCols(DenseMatrix &V) const
+{
+   for (int c=0; c<V.Width(); c++)
+   {
+      TransformDual(V.GetColumn(c));
+   }
+}
+
+void DofTransformation::InvTransformPrimal(Vector &v) const
+{
+   InvTransformPrimal(v.GetData());
+}
+
+void TransformPrimal(const DofTransformation *ran_dof_trans,
+                     const DofTransformation *dom_dof_trans,
+                     DenseMatrix &elmat)
+{
+   if (ran_dof_trans && dom_dof_trans)
+   {
+      ran_dof_trans->TransformPrimalCols(elmat);
+      dom_dof_trans->TransformDualRows(elmat);
+   }
+   else if (ran_dof_trans)
+   {
+      ran_dof_trans->TransformPrimalCols(elmat);
+   }
+   else if (dom_dof_trans)
+   {
+      dom_dof_trans->TransformDualRows(elmat);
+   }
+   else
+   {
+      // If both transformations are NULL this function should not be called
+   }
+}
+
+void TransformDual(const DofTransformation *ran_dof_trans,
+                   const DofTransformation *dom_dof_trans,
+                   DenseMatrix &elmat)
+{
+   if (ran_dof_trans && dom_dof_trans)
+   {
+      ran_dof_trans->TransformDualCols(elmat);
+      dom_dof_trans->TransformDualRows(elmat);
+   }
+   else if (ran_dof_trans)
+   {
+      ran_dof_trans->TransformDualCols(elmat);
+   }
+   else if (dom_dof_trans)
+   {
+      dom_dof_trans->TransformDualRows(elmat);
+   }
+   else
+   {
+      // If both transformations are NULL this function should not be called
+   }
+}
+
+void VDofTransformation::TransformPrimal(double *v) const
+{
+   int size = doftrans_->Size();
+
+   if ((Ordering::Type)ordering_ == Ordering::byNODES || vdim_ == 1)
+   {
+      for (int i=0; i<vdim_; i++)
+      {
+         doftrans_->TransformPrimal(&v[i*size]);
+      }
+   }
+   else
+   {
+      Vector vec(size);
+      for (int i=0; i<vdim_; i++)
+      {
+         for (int j=0; j<size; j++)
+         {
+            vec(j) = v[j*vdim_+i];
+         }
+         doftrans_->TransformPrimal(vec);
+         for (int j=0; j<size; j++)
+         {
+            v[j*vdim_+i] = vec(j);
+         }
+      }
+   }
+}
+
+void VDofTransformation::InvTransformPrimal(double *v) const
+{
+   int size = doftrans_->Height();
+
+   if ((Ordering::Type)ordering_ == Ordering::byNODES)
+   {
+      for (int i=0; i<vdim_; i++)
+      {
+         doftrans_->InvTransformPrimal(&v[i*size]);
+      }
+   }
+   else
+   {
+      Vector vec(size);
+      for (int i=0; i<vdim_; i++)
+      {
+         for (int j=0; j<size; j++)
+         {
+            vec(j) = v[j*vdim_+i];
+         }
+         doftrans_->InvTransformPrimal(vec);
+         for (int j=0; j<size; j++)
+         {
+            v[j*vdim_+i] = vec(j);
+         }
+      }
+   }
+}
+
+void VDofTransformation::TransformDual(double *v) const
+{
+   int size = doftrans_->Size();
+
+   if ((Ordering::Type)ordering_ == Ordering::byNODES)
+   {
+      for (int i=0; i<vdim_; i++)
+      {
+         doftrans_->TransformDual(&v[i*size]);
+      }
+   }
+   else
+   {
+      Vector vec(size);
+      for (int i=0; i<vdim_; i++)
+      {
+         for (int j=0; j<size; j++)
+         {
+            vec(j) = v[j*vdim_+i];
+         }
+         doftrans_->TransformDual(vec);
+         for (int j=0; j<size; j++)
+         {
+            v[j*vdim_+i] = vec(j);
+         }
+      }
+   }
+}
+
+const double ND_DofTransformation::T_data[24] =
+{
+   1.0,  0.0,  0.0,  1.0,
+   -1.0, -1.0,  0.0,  1.0,
+   0.0,  1.0, -1.0, -1.0,
+   1.0,  0.0, -1.0, -1.0,
+   -1.0, -1.0,  1.0,  0.0,
+   0.0,  1.0,  1.0,  0.0
+};
+
+const DenseTensor ND_DofTransformation
+::T(const_cast<double*>(ND_DofTransformation::T_data), 2, 2, 6);
+
+const double ND_DofTransformation::TInv_data[24] =
+{
+   1.0,  0.0,  0.0,  1.0,
+   -1.0, -1.0,  0.0,  1.0,
+   -1.0, -1.0,  1.0,  0.0,
+   1.0,  0.0, -1.0, -1.0,
+   0.0,  1.0, -1.0, -1.0,
+   0.0,  1.0,  1.0,  0.0
+};
+
+const DenseTensor ND_DofTransformation
+::TInv(const_cast<double*>(TInv_data), 2, 2, 6);
+
+ND_DofTransformation::ND_DofTransformation(int size, int p)
+   : DofTransformation(size),
+     order(p)
+{
+}
+
+ND_TriDofTransformation::ND_TriDofTransformation(int p)
+   : ND_DofTransformation(p*(p + 2), p)
+{
+}
+
+void ND_TriDofTransformation::TransformPrimal(double *v) const
+{
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face
+
+   double data[2];
+   Vector v2(data, 2);
+
+   // Transform face DoFs
+   for (int f=0; f<1; f++)
+   {
+      for (int i=0; i<nfdofs/2; i++)
+      {
+         v2 = &v[3*nedofs + f*nfdofs + 2*i];
+         T(Fo[f]).Mult(v2, &v[3*nedofs + f*nfdofs + 2*i]);
+      }
+   }
+}
+
+void
+ND_TriDofTransformation::InvTransformPrimal(double *v) const
+{
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face
+
+   double data[2];
+   Vector v2(data, 2);
+
+   // Transform face DoFs
+   for (int f=0; f<1; f++)
+   {
+      for (int i=0; i<nfdofs/2; i++)
+      {
+         v2 = &v[3*nedofs + f*nfdofs + 2*i];
+         TInv(Fo[f]).Mult(v2, &v[3*nedofs + f*nfdofs + 2*i]);
+      }
+   }
+}
+
+void
+ND_TriDofTransformation::TransformDual(double *v) const
+{
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face
+
+   double data[2];
+   Vector v2(data, 2);
+
+   // Transform face DoFs
+   for (int f=0; f<1; f++)
+   {
+      for (int i=0; i<nfdofs/2; i++)
+      {
+         v2 = &v[3*nedofs + f*nfdofs + 2*i];
+         TInv(Fo[f]).MultTranspose(v2, &v[3*nedofs + f*nfdofs + 2*i]);
+      }
+   }
+}
+
+ND_TetDofTransformation::ND_TetDofTransformation(int p)
+   : ND_DofTransformation(p*(p + 2)*(p + 3)/2, p)
+{
+}
+
+void ND_TetDofTransformation::TransformPrimal(double *v) const
+{
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face
+
+   double data[2];
+   Vector v2(data, 2);
+
+   // Transform face DoFs
+   for (int f=0; f<4; f++)
+   {
+      for (int i=0; i<nfdofs/2; i++)
+      {
+         v2 = &v[6*nedofs + f*nfdofs + 2*i];
+         T(Fo[f]).Mult(v2, &v[6*nedofs + f*nfdofs + 2*i]);
+      }
+   }
+}
+
+void
+ND_TetDofTransformation::InvTransformPrimal(double *v) const
+{
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face
+
+   double data[2];
+   Vector v2(data, 2);
+
+   // Transform face DoFs
+   for (int f=0; f<4; f++)
+   {
+      for (int i=0; i<nfdofs/2; i++)
+      {
+         v2 = &v[6*nedofs + f*nfdofs + 2*i];
+         TInv(Fo[f]).Mult(v2, &v[6*nedofs + f*nfdofs + 2*i]);
+      }
+   }
+}
+
+void
+ND_TetDofTransformation::TransformDual(double *v) const
+{
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face
+
+   double data[2];
+   Vector v2(data, 2);
+
+   // Transform face DoFs
+   for (int f=0; f<4; f++)
+   {
+      for (int i=0; i<nfdofs/2; i++)
+      {
+         v2 = &v[6*nedofs + f*nfdofs + 2*i];
+         TInv(Fo[f]).MultTranspose(v2, &v[6*nedofs + f*nfdofs + 2*i]);
+      }
+   }
+}
+
+} // namespace mfem

fem/doftrans.hpp

@@ -0,0 +1,277 @@
+// 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.
+
+#ifndef MFEM_DOFTRANSFORM
+#define MFEM_DOFTRANSFORM
+
+#include "../config/config.hpp"
+#include "../linalg/linalg.hpp"
+#include "intrules.hpp"
+#include "fe.hpp"
+
+namespace mfem
+{
+
+/** The DofTransformation class is an abstract base class for a family of
+    transformations that map local degrees of freedom (DoFs), contained within
+    individual elements, to global degrees of freedom, stored within
+    GridFunction objects. These transformations are necessary to ensure that
+    basis functions in neighboring elements align correctly. Closely related but
+    complementary transformations are required for the entries stored in
+    LinearForm and BilinearForm objects. The DofTransformation class is designed
+    to apply the action of both of these types of DoF transformations.
+
+    Let the "primal transformation" be given by the operator T. This means that
+    given a local element vector v the data that must be placed into a
+    GridFunction object is v_t = T * v.
+
+    We also need the inverse of the primal transformation T^{-1} so that we can
+    recover the local element vector from data read out of a GridFunction
+    e.g. v = T^{-1} * v_t.
+
+    We need to preserve the action of our linear forms applied to primal
+    vectors. In other words, if f is the local vector computed by a linear
+    form then f * v = f_t * v_t (where "*" represents an inner product of
+    vectors). This requires that f_t = T^{-T} * f i.e. the "dual transform" is
+    given by the transpose of the inverse of the primal transformation.
+
+    For bilinear forms we require that v^T * A * v = v_t^T * A_t * v_t. This
+    implies that A_t = T^{-T} * A * T^{-1}. This can be accomplished by
+    performing dual transformations of the rows and columns of the matrix A.
+
+    For discrete linear operators the range must be modified with the primal
+    transformation rather than the dual transformation because the result is a
+    primal vector rather than a dual vector. This leads to the transformation
+    D_t = T * D * T^{-1}. This can be accomplished by using a primal
+    transformation on the columns of D and a dual transformation on its rows.
+*/
+class DofTransformation
+{
+protected:
+   int size_;
+
+   Array<int> Fo;
+
+   DofTransformation(int size)
+      : size_(size) {}
+
+public:
+
+   inline int Size() const { return size_; }
+   inline int Height() const { return size_; }
+   inline int NumRows() const { return size_; }
+   inline int Width() const { return size_; }
+   inline int NumCols() const { return size_; }
+
+   /** @brief Configure the transformation using face orientations for the
+       current element. */
+   /// The face_orientation array can be obtained from Mesh::GetElementFaces.
+   inline void SetFaceOrientations(const Array<int> & face_orientation)
+   { Fo = face_orientation; }
+
+   inline const Array<int> & GetFaceOrientations() const { return Fo; }
+
+   /** Transform local DoFs to align with the global DoFs. For example, this
+       transformation can be used to map the local vector computed by
+       FiniteElement::Project() to the transformed vector stored within a
+       GridFunction object. */
+   virtual void TransformPrimal(double *v) const = 0;
+   virtual void TransformPrimal(Vector &v) const;
+
+   /// Transform groups of DoFs stored as dense matrices
+   virtual void TransformPrimalCols(DenseMatrix &V) const;
+
+   /** Inverse transform local DoFs. Used to transform DoFs from a global vector
+       back to their element-local form. For example, this must be used to
+       transform the vector obtained using GridFunction::GetSubVector before it
+       can be used to compute a local interpolation.
+   */
+   virtual void InvTransformPrimal(double *v) const = 0;
+   virtual void InvTransformPrimal(Vector &v) const;
+
+   /** Transform dual DoFs as computed by a LinearFormIntegrator before summing
+       into a LinearForm object. */
+   virtual void TransformDual(double *v) const = 0;
+   virtual void TransformDual(Vector &v) const;
+
+   /** Transform a matrix of dual DoFs entries as computed by a
+       BilinearFormIntegrator before summing into a BilinearForm object. */
+   virtual void TransformDual(DenseMatrix &V) const;
+
+   /// Transform groups of dual DoFs stored as dense matrices
+   virtual void TransformDualRows(DenseMatrix &V) const;
+   virtual void TransformDualCols(DenseMatrix &V) const;
+
+   virtual ~DofTransformation() {}
+};
+
+/** Transform a matrix of DoFs entries from different finite element spaces as
+    computed by a DiscreteInterpolator before copying into a
+    DiscreteLinearOperator.
+*/
+void TransformPrimal(const DofTransformation *ran_dof_trans,
+                     const DofTransformation *dom_dof_trans,
+                     DenseMatrix &elmat);
+
+/** Transform a matrix of dual DoFs entries from different finite element spaces
+    as computed by a BilinearFormIntegrator before summing into a
+    MixedBilinearForm object.
+*/
+void TransformDual(const DofTransformation *ran_dof_trans,
+                   const DofTransformation *dom_dof_trans,
+                   DenseMatrix &elmat);
+
+/** The VDofTransformation class implements a nested transformation where an
+    arbitrary DofTransformation is replicated with a vdim >= 1.
+*/
+class VDofTransformation : public DofTransformation
+{
+private:
+   int vdim_;
+   int ordering_;
+   DofTransformation * doftrans_;
+
+public:
+   /** @brief Default constructor which requires that SetDofTransformation be
+       called before use. */
+   VDofTransformation(int vdim = 1, int ordering = 0)
+      : DofTransformation(0),
+        vdim_(vdim), ordering_(ordering),
+        doftrans_(NULL) {}
+
+   /// Constructor with a known DofTransformation
+   VDofTransformation(DofTransformation & doftrans, int vdim = 1,
+                      int ordering = 0)
+      : DofTransformation(vdim * doftrans.Size()),
+        vdim_(vdim), ordering_(ordering),
+        doftrans_(&doftrans) {}
+
+   /// Set or change the vdim parameter
+   inline void SetVDim(int vdim)
+   {
+      vdim_ = vdim;
+      if (doftrans_)
+      {
+         size_ = vdim_ * doftrans_->Size();
+      }
+   }
+
+   /// Return the current vdim value
+   inline int GetVDim() const { return vdim_; }
+
+   /// Set or change the nested DofTransformation object
+   inline void SetDofTransformation(DofTransformation & doftrans)
+   {
+      size_ = vdim_ * doftrans.Size();
+      doftrans_ = &doftrans;
+   }
+
+   /// Return the nested DofTransformation object
+   inline DofTransformation * GetDofTransformation() const { return doftrans_; }
+
+   inline void SetFaceOrientation(const Array<int> & face_orientation)
+   { Fo = face_orientation; doftrans_->SetFaceOrientations(face_orientation); }
+
+   using DofTransformation::TransformPrimal;
+   using DofTransformation::InvTransformPrimal;
+   using DofTransformation::TransformDual;
+
+   void TransformPrimal(double *v) const;
+   void InvTransformPrimal(double *v) const;
+   void TransformDual(double *v) const;
+};
+
+/** Abstract base class for high-order Nedelec spaces on elements with
+    triangular faces.
+
+    The Nedelec DoFs on the interior of triangular faces come in pairs which
+    share an interpolation point but have different vector directions. These
+    directions depend on the orientation of the face and can therefore differ in
+    neighboring elements. The mapping required to transform these DoFs can be
+    implemented as series of 2x2 linear transformations. The raw data for these
+    linear transformations is stored in the T_data and TInv_data arrays and can
+    be accessed as DenseMatrices using the GetFaceTransform() and
+    GetFaceInverseTransform() methods.
+*/
+class ND_DofTransformation : public DofTransformation
+{
+protected:
+   static const double T_data[24];
+   static const double TInv_data[24];
+   static const DenseTensor T, TInv;
+   int order;
+
+   ND_DofTransformation(int size, int order);
+
+public:
+   // Return the 2x2 transformation operator for the given face orientation
+   static const DenseMatrix & GetFaceTransform(int ori) { return T(ori); }
+
+   // Return the 2x2 inverse transformation operator
+   static const DenseMatrix & GetFaceInverseTransform(int ori)
+   { return TInv(ori); }
+};
+
+/// DoF transformation implementation for the Nedelec basis on triangles
+class ND_TriDofTransformation : public ND_DofTransformation
+{
+public:
+   ND_TriDofTransformation(int order);
+
+   using DofTransformation::TransformPrimal;
+   using DofTransformation::InvTransformPrimal;
+   using DofTransformation::TransformDual;
+
+   void TransformPrimal(double *v) const;
+
+   void InvTransformPrimal(double *v) const;
+
+   void TransformDual(double *v) const;
+};
+
+/// DoF transformation implementation for the Nedelec basis on tetrahedra
+class ND_TetDofTransformation : public ND_DofTransformation
+{
+public:
+   ND_TetDofTransformation(int order);
+
+   using DofTransformation::TransformPrimal;
+   using DofTransformation::InvTransformPrimal;
+   using DofTransformation::TransformDual;
+
+   void TransformPrimal(double *v) const;
+
+   void InvTransformPrimal(double *v) const;
+
+   void TransformDual(double *v) const;
+};
+
+/// DoF transformation implementation for the Nedelec basis on wedge elements
+/** TODO: (Under development) */
+class ND_WedgeDofTransformation : public ND_DofTransformation
+{
+public:
+   ND_WedgeDofTransformation(int order);
+
+   using DofTransformation::TransformPrimal;
+   using DofTransformation::InvTransformPrimal;
+   using DofTransformation::TransformDual;
+
+   void TransformPrimal(double *v) const;
+
+   void InvTransformPrimal(double *v) const;
+
+   void TransformDual(double *v) const;
+};
+
+} // namespace mfem
+
+#endif // MFEM_DOFTRANSFORM

fem/eltrans.cpp

@@ -380,6 +380,7 @@ void IsoparametricTransformation::SetIdentityTransformation(
       case Geometry::TETRAHEDRON : FElem = &TetrahedronFE; break;
       case Geometry::CUBE :        FElem = &HexahedronFE; break;
       case Geometry::PRISM :       FElem = &WedgeFE; break;
+      case Geometry::PYRAMID :     FElem = &PyramidFE; break;
       default:
          MFEM_ABORT("unknown Geometry::Type!");
    }

fem/estimators.cpp

@@ -329,7 +329,7 @@ void KellyErrorEstimator::ComputeEstimates()
          error_estimates(e) = sqrt(factor * error_estimates(e));
       }
 
-      total_error = error_estimates.Sum();
+      total_error = error_estimates.Norml2();
       delete flux;
       return;
    }
@@ -452,9 +452,10 @@ void KellyErrorEstimator::ComputeEstimates()
    auto pfes = dynamic_cast<ParFiniteElementSpace*>(xfes);
    MFEM_VERIFY(pfes, "xfes is not a ParFiniteElementSpace pointer");
 
-   double process_local_error = error_estimates.Sum();
+   double process_local_error = pow(error_estimates.Norml2(),2.0);
    MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
                  MPI_SUM, pfes->GetComm());
+   total_error = sqrt(total_error);
 #endif // MFEM_USE_MPI
 }
 

fem/fe/fe_h1.hpp

@@ -0,0 +1,153 @@
+// 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.
+
+#ifndef MFEM_FE_H1
+#define MFEM_FE_H1
+
+#include "fe_base.hpp"
+
+namespace mfem
+{
+
+/// Arbitrary order H1 elements in 1D
+class H1_SegmentElement : public NodalTensorFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, dshape_x, d2shape_x;
+#endif
+
+public:
+   /// Construct the H1_SegmentElement of order @a p and BasisType @a btype
+   H1_SegmentElement(const int p, const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &Hessian) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+};
+
+
+/// Arbitrary order H1 elements in 2D on a square
+class H1_QuadrilateralElement : public NodalTensorFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, dshape_x, dshape_y, d2shape_x, d2shape_y;
+#endif
+
+public:
+   /// Construct the H1_QuadrilateralElement of order @a p and BasisType @a btype
+   H1_QuadrilateralElement(const int p,
+                           const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &Hessian) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+};
+
+
+/// Arbitrary order H1 elements in 3D on a cube
+class H1_HexahedronElement : public NodalTensorFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z, dshape_x, dshape_y, dshape_z,
+           d2shape_x, d2shape_y, d2shape_z;
+#endif
+
+public:
+   /// Construct the H1_HexahedronElement of order @a p and BasisType @a btype
+   H1_HexahedronElement(const int p, const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &Hessian) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+};
+
+
+/// Arbitrary order H1 elements in 2D on a triangle
+class H1_TriangleElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
+   mutable Vector ddshape_x, ddshape_y, ddshape_l;
+   mutable DenseMatrix du, ddu;
+#endif
+   DenseMatrixInverse Ti;
+
+public:
+   /// Construct the H1_TriangleElement of order @a p and BasisType @a btype
+   H1_TriangleElement(const int p, const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &ddshape) const;
+};
+
+
+/// Arbitrary order H1 elements in 3D  on a tetrahedron
+class H1_TetrahedronElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z, shape_l;
+   mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
+   mutable Vector ddshape_x, ddshape_y, ddshape_z, ddshape_l;
+   mutable DenseMatrix du, ddu;
+#endif
+   DenseMatrixInverse Ti;
+
+public:
+   /// Construct the H1_TetrahedronElement of order @a p and BasisType @a btype
+   H1_TetrahedronElement(const int p,
+                         const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &ddshape) const;
+};
+
+
+
+/// Arbitrary order H1 elements in 3D on a wedge
+class H1_WedgeElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector t_shape, s_shape;
+   mutable DenseMatrix t_dshape, s_dshape;
+#endif
+   Array<int> t_dof, s_dof;
+
+   H1_TriangleElement TriangleFE;
+   H1_SegmentElement  SegmentFE;
+
+public:
+   /// Construct the H1_WedgeElement of order @a p and BasisType @a btype
+   H1_WedgeElement(const int p,
+                   const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+};
+
+} // namespace mfem
+
+#endif

fem/fe/fe_l2.cpp

@@ -0,0 +1,694 @@
+// 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.
+
+// L2 Finite Element classes
+
+#include "fe_l2.hpp"
+#include "fe_h1.hpp"
+
+namespace mfem
+{
+
+using namespace std;
+
+L2_SegmentElement::L2_SegmentElement(const int p, const int btype)
+   : NodalTensorFiniteElement(1, p, VerifyOpen(btype), L2_DOF_MAP)
+{
+   const double *op = poly1d.OpenPoints(p, btype);
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(p + 1);
+   dshape_x.SetDataAndSize(NULL, p + 1);
+#endif
+
+   for (int i = 0; i <= p; i++)
+   {
+      Nodes.IntPoint(i).x = op[i];
+   }
+}
+
+void L2_SegmentElement::CalcShape(const IntegrationPoint &ip,
+                                  Vector &shape) const
+{
+   basis1d.Eval(ip.x, shape);
+}
+
+void L2_SegmentElement::CalcDShape(const IntegrationPoint &ip,
+                                   DenseMatrix &dshape) const
+{
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(dof), dshape_x(dshape.Data(), dof);
+#else
+   dshape_x.SetData(dshape.Data());
+#endif
+   basis1d.Eval(ip.x, shape_x, dshape_x);
+}
+
+void L2_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
+{
+   const int p = order;
+   const double *op = poly1d.OpenPoints(p, b_type);
+
+   switch (vertex)
+   {
+      case 0:
+         for (int i = 0; i <= p; i++)
+         {
+            dofs(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
+         }
+         break;
+
+      case 1:
+         for (int i = 0; i <= p; i++)
+         {
+            dofs(i) = poly1d.CalcDelta(p,op[i]);
+         }
+         break;
+   }
+}
+
+
+L2_QuadrilateralElement::L2_QuadrilateralElement(const int p, const int btype)
+   : NodalTensorFiniteElement(2, p, VerifyOpen(btype), L2_DOF_MAP)
+{
+   const double *op = poly1d.OpenPoints(p, b_type);
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(p + 1);
+   shape_y.SetSize(p + 1);
+   dshape_x.SetSize(p + 1);
+   dshape_y.SetSize(p + 1);
+#endif
+
+   for (int o = 0, j = 0; j <= p; j++)
+      for (int i = 0; i <= p; i++)
+      {
+         Nodes.IntPoint(o++).Set2(op[i], op[j]);
+      }
+}
+
+void L2_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
+                                        Vector &shape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1), shape_y(p+1);
+#endif
+
+   basis1d.Eval(ip.x, shape_x);
+   basis1d.Eval(ip.y, shape_y);
+
+   for (int o = 0, j = 0; j <= p; j++)
+      for (int i = 0; i <= p; i++)
+      {
+         shape(o++) = shape_x(i)*shape_y(j);
+      }
+}
+
+void L2_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
+                                         DenseMatrix &dshape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
+#endif
+
+   basis1d.Eval(ip.x, shape_x, dshape_x);
+   basis1d.Eval(ip.y, shape_y, dshape_y);
+
+   for (int o = 0, j = 0; j <= p; j++)
+      for (int i = 0; i <= p; i++)
+      {
+         dshape(o,0) = dshape_x(i)* shape_y(j);
+         dshape(o,1) =  shape_x(i)*dshape_y(j);  o++;
+      }
+}
+
+void L2_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
+{
+   const int p = order;
+   const double *op = poly1d.OpenPoints(p, b_type);
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1), shape_y(p+1);
+#endif
+
+   for (int i = 0; i <= p; i++)
+   {
+      shape_x(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
+      shape_y(i) = poly1d.CalcDelta(p,op[i]);
+   }
+
+   switch (vertex)
+   {
+      case 0:
+         for (int o = 0, j = 0; j <= p; j++)
+            for (int i = 0; i <= p; i++)
+            {
+               dofs[o++] = shape_x(i)*shape_x(j);
+            }
+         break;
+      case 1:
+         for (int o = 0, j = 0; j <= p; j++)
+            for (int i = 0; i <= p; i++)
+            {
+               dofs[o++] = shape_y(i)*shape_x(j);
+            }
+         break;
+      case 2:
+         for (int o = 0, j = 0; j <= p; j++)
+            for (int i = 0; i <= p; i++)
+            {
+               dofs[o++] = shape_y(i)*shape_y(j);
+            }
+         break;
+      case 3:
+         for (int o = 0, j = 0; j <= p; j++)
+            for (int i = 0; i <= p; i++)
+            {
+               dofs[o++] = shape_x(i)*shape_y(j);
+            }
+         break;
+   }
+}
+
+
+L2_HexahedronElement::L2_HexahedronElement(const int p, const int btype)
+   : NodalTensorFiniteElement(3, p, VerifyOpen(btype), L2_DOF_MAP)
+{
+   const double *op = poly1d.OpenPoints(p, btype);
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(p + 1);
+   shape_y.SetSize(p + 1);
+   shape_z.SetSize(p + 1);
+   dshape_x.SetSize(p + 1);
+   dshape_y.SetSize(p + 1);
+   dshape_z.SetSize(p + 1);
+#endif
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j <= p; j++)
+         for (int i = 0; i <= p; i++)
+         {
+            Nodes.IntPoint(o++).Set3(op[i], op[j], op[k]);
+         }
+}
+
+void L2_HexahedronElement::CalcShape(const IntegrationPoint &ip,
+                                     Vector &shape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
+#endif
+
+   basis1d.Eval(ip.x, shape_x);
+   basis1d.Eval(ip.y, shape_y);
+   basis1d.Eval(ip.z, shape_z);
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j <= p; j++)
+         for (int i = 0; i <= p; i++)
+         {
+            shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
+         }
+}
+
+void L2_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
+                                      DenseMatrix &dshape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1),  shape_y(p+1),  shape_z(p+1);
+   Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
+#endif
+
+   basis1d.Eval(ip.x, shape_x, dshape_x);
+   basis1d.Eval(ip.y, shape_y, dshape_y);
+   basis1d.Eval(ip.z, shape_z, dshape_z);
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j <= p; j++)
+         for (int i = 0; i <= p; i++)
+         {
+            dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
+            dshape(o,1) =  shape_x(i)*dshape_y(j)* shape_z(k);
+            dshape(o,2) =  shape_x(i)* shape_y(j)*dshape_z(k);  o++;
+         }
+}
+
+void L2_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
+{
+   const int p = order;
+   const double *op = poly1d.OpenPoints(p, b_type);
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p+1), shape_y(p+1);
+#endif
+
+   for (int i = 0; i <= p; i++)
+   {
+      shape_x(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
+      shape_y(i) = poly1d.CalcDelta(p,op[i]);
+   }
+
+   switch (vertex)
+   {
+      case 0:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
+               }
+         break;
+      case 1:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
+               }
+         break;
+      case 2:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
+               }
+         break;
+      case 3:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
+               }
+         break;
+      case 4:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
+               }
+         break;
+      case 5:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
+               }
+         break;
+      case 6:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
+               }
+         break;
+      case 7:
+         for (int o = 0, k = 0; k <= p; k++)
+            for (int j = 0; j <= p; j++)
+               for (int i = 0; i <= p; i++)
+               {
+                  dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
+               }
+         break;
+   }
+}
+
+
+L2_TriangleElement::L2_TriangleElement(const int p, const int btype)
+   : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
+                        FunctionSpace::Pk)
+{
+   const double *op = poly1d.OpenPoints(p, VerifyOpen(btype));
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(p + 1);
+   shape_y.SetSize(p + 1);
+   shape_l.SetSize(p + 1);
+   dshape_x.SetSize(p + 1);
+   dshape_y.SetSize(p + 1);
+   dshape_l.SetSize(p + 1);
+   u.SetSize(dof);
+   du.SetSize(dof, dim);
+#else
+   Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
+#endif
+
+   for (int o = 0, j = 0; j <= p; j++)
+      for (int i = 0; i + j <= p; i++)
+      {
+         double w = op[i] + op[j] + op[p-i-j];
+         Nodes.IntPoint(o++).Set2(op[i]/w, op[j]/w);
+      }
+
+   DenseMatrix T(dof);
+   for (int k = 0; k < dof; k++)
+   {
+      IntegrationPoint &ip = Nodes.IntPoint(k);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
+
+      for (int o = 0, j = 0; j <= p; j++)
+         for (int i = 0; i + j <= p; i++)
+         {
+            T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
+         }
+   }
+
+   Ti.Factor(T);
+   // mfem::out << "L2_TriangleElement(" << p << ") : "; Ti.TestInversion();
+}
+
+void L2_TriangleElement::CalcShape(const IntegrationPoint &ip,
+                                   Vector &shape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(dof);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
+
+   for (int o = 0, j = 0; j <= p; j++)
+      for (int i = 0; i + j <= p; i++)
+      {
+         u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
+      }
+
+   Ti.Mult(u, shape);
+}
+
+void L2_TriangleElement::CalcDShape(const IntegrationPoint &ip,
+                                    DenseMatrix &dshape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector  shape_x(p + 1),  shape_y(p + 1),  shape_l(p + 1);
+   Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
+   DenseMatrix du(dof, dim);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
+
+   for (int o = 0, j = 0; j <= p; j++)
+      for (int i = 0; i + j <= p; i++)
+      {
+         int k = p - i - j;
+         du(o,0) = ((dshape_x(i)* shape_l(k)) -
+                    ( shape_x(i)*dshape_l(k)))*shape_y(j);
+         du(o,1) = ((dshape_y(j)* shape_l(k)) -
+                    ( shape_y(j)*dshape_l(k)))*shape_x(i);
+         o++;
+      }
+
+   Ti.Mult(du, dshape);
+}
+
+void L2_TriangleElement::ProjectDelta(int vertex, Vector &dofs) const
+{
+   switch (vertex)
+   {
+      case 0:
+         for (int i = 0; i < dof; i++)
+         {
+            const IntegrationPoint &ip = Nodes.IntPoint(i);
+            dofs[i] = pow(1.0 - ip.x - ip.y, order);
+         }
+         break;
+      case 1:
+         for (int i = 0; i < dof; i++)
+         {
+            const IntegrationPoint &ip = Nodes.IntPoint(i);
+            dofs[i] = pow(ip.x, order);
+         }
+         break;
+      case 2:
+         for (int i = 0; i < dof; i++)
+         {
+            const IntegrationPoint &ip = Nodes.IntPoint(i);
+            dofs[i] = pow(ip.y, order);
+         }
+         break;
+   }
+}
+
+
+L2_TetrahedronElement::L2_TetrahedronElement(const int p, const int btype)
+   : NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
+                        p, FunctionSpace::Pk)
+{
+   const double *op = poly1d.OpenPoints(p, VerifyNodal(VerifyOpen(btype)));
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(p + 1);
+   shape_y.SetSize(p + 1);
+   shape_z.SetSize(p + 1);
+   shape_l.SetSize(p + 1);
+   dshape_x.SetSize(p + 1);
+   dshape_y.SetSize(p + 1);
+   dshape_z.SetSize(p + 1);
+   dshape_l.SetSize(p + 1);
+   u.SetSize(dof);
+   du.SetSize(dof, dim);
+#else
+   Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
+#endif
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j + k <= p; j++)
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            double w = op[i] + op[j] + op[k] + op[p-i-j-k];
+            Nodes.IntPoint(o++).Set3(op[i]/w, op[j]/w, op[k]/w);
+         }
+
+   DenseMatrix T(dof);
+   for (int m = 0; m < dof; m++)
+   {
+      IntegrationPoint &ip = Nodes.IntPoint(m);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, ip.z, shape_z);
+      poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
+
+      for (int o = 0, k = 0; k <= p; k++)
+         for (int j = 0; j + k <= p; j++)
+            for (int i = 0; i + j + k <= p; i++)
+            {
+               T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
+            }
+   }
+
+   Ti.Factor(T);
+   // mfem::out << "L2_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
+}
+
+void L2_TetrahedronElement::CalcShape(const IntegrationPoint &ip,
+                                      Vector &shape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
+   Vector u(dof);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j + k <= p; j++)
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
+         }
+
+   Ti.Mult(u, shape);
+}
+
+void L2_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
+                                       DenseMatrix &dshape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector  shape_x(p + 1),  shape_y(p + 1),  shape_z(p + 1),  shape_l(p + 1);
+   Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
+   DenseMatrix du(dof, dim);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
+
+   for (int o = 0, k = 0; k <= p; k++)
+      for (int j = 0; j + k <= p; j++)
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            int l = p - i - j - k;
+            du(o,0) = ((dshape_x(i)* shape_l(l)) -
+                       ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
+            du(o,1) = ((dshape_y(j)* shape_l(l)) -
+                       ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
+            du(o,2) = ((dshape_z(k)* shape_l(l)) -
+                       ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
+            o++;
+         }
+
+   Ti.Mult(du, dshape);
+}
+
+void L2_TetrahedronElement::ProjectDelta(int vertex, Vector &dofs) const
+{
+   switch (vertex)
+   {
+      case 0:
+         for (int i = 0; i < dof; i++)
+         {
+            const IntegrationPoint &ip = Nodes.IntPoint(i);
+            dofs[i] = pow(1.0 - ip.x - ip.y - ip.z, order);
+         }
+         break;
+      case 1:
+         for (int i = 0; i < dof; i++)
+         {
+            const IntegrationPoint &ip = Nodes.IntPoint(i);
+            dofs[i] = pow(ip.x, order);
+         }
+         break;
+      case 2:
+         for (int i = 0; i < dof; i++)
+         {
+            const IntegrationPoint &ip = Nodes.IntPoint(i);
+            dofs[i] = pow(ip.y, order);
+         }
+         break;
+      case 3:
+         for (int i = 0; i < dof; i++)
+         {
+            const IntegrationPoint &ip = Nodes.IntPoint(i);
+            dofs[i] = pow(ip.z, order);
+         }
+         break;
+   }
+}
+
+
+L2_WedgeElement::L2_WedgeElement(const int p, const int btype)
+   : NodalFiniteElement(3, Geometry::PRISM, ((p + 1)*(p + 1)*(p + 2))/2,
+                        p, FunctionSpace::Qk),
+     TriangleFE(p, btype),
+     SegmentFE(p, btype)
+{
+#ifndef MFEM_THREAD_SAFE
+   t_shape.SetSize(TriangleFE.GetDof());
+   s_shape.SetSize(SegmentFE.GetDof());
+   t_dshape.SetSize(TriangleFE.GetDof(), 2);
+   s_dshape.SetSize(SegmentFE.GetDof(), 1);
+#endif
+
+   t_dof.SetSize(dof);
+   s_dof.SetSize(dof);
+
+   // Interior DoFs
+   int m=0;
+   for (int k=0; k<=p; k++)
+   {
+      int l=0;
+      for (int j=0; j<=p; j++)
+      {
+         for (int i=0; i<=j; i++)
+         {
+            t_dof[m] = l;
+            s_dof[m] = k;
+            l++; m++;
+         }
+      }
+   }
+
+   // Define Nodes
+   const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
+   const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
+   for (int i=0; i<dof; i++)
+   {
+      Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
+      Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
+      Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
+   }
+}
+
+void L2_WedgeElement::CalcShape(const IntegrationPoint &ip,
+                                Vector &shape) const
+{
+#ifdef MFEM_THREAD_SAFE
+   Vector t_shape(TriangleFE.GetDof());
+   Vector s_shape(SegmentFE.GetDof());
+#endif
+
+   IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
+
+   TriangleFE.CalcShape(ip, t_shape);
+   SegmentFE.CalcShape(ipz, s_shape);
+
+   for (int i=0; i<dof; i++)
+   {
+      shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
+   }
+}
+
+void L2_WedgeElement::CalcDShape(const IntegrationPoint &ip,
+                                 DenseMatrix &dshape) const
+{
+#ifdef MFEM_THREAD_SAFE
+   Vector      t_shape(TriangleFE.GetDof());
+   DenseMatrix t_dshape(TriangleFE.GetDof(), 2);
+   Vector      s_shape(SegmentFE.GetDof());
+   DenseMatrix s_dshape(SegmentFE.GetDof(), 1);
+#endif
+
+   IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
+
+   TriangleFE.CalcShape(ip, t_shape);
+   TriangleFE.CalcDShape(ip, t_dshape);
+   SegmentFE.CalcShape(ipz, s_shape);
+   SegmentFE.CalcDShape(ipz, s_dshape);
+
+   for (int i=0; i<dof; i++)
+   {
+      dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
+      dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
+      dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
+   }
+}
+
+}

fem/fe/fe_l2.hpp

@@ -0,0 +1,151 @@
+// 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.
+
+#ifndef MFEM_FE_L2
+#define MFEM_FE_L2
+
+#include "fe_base.hpp"
+
+namespace mfem
+{
+
+/// Arbitrary order L2 elements in 1D on a segment
+class L2_SegmentElement : public NodalTensorFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, dshape_x;
+#endif
+
+public:
+   /// Construct the L2_SegmentElement of order @a p and BasisType @a btype
+   L2_SegmentElement(const int p, const int btype = BasisType::GaussLegendre);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+};
+
+
+/// Arbitrary order L2 elements in 2D on a square
+class L2_QuadrilateralElement : public NodalTensorFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, dshape_x, dshape_y;
+#endif
+
+public:
+   /// Construct the L2_QuadrilateralElement of order @a p and BasisType @a btype
+   L2_QuadrilateralElement(const int p,
+                           const int btype = BasisType::GaussLegendre);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   virtual void ProjectCurl(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &curl) const
+   { ProjectCurl_2D(fe, Trans, curl); }
+};
+
+
+/// Arbitrary order L2 elements in 3D on a cube
+class L2_HexahedronElement : public NodalTensorFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z, dshape_x, dshape_y, dshape_z;
+#endif
+
+public:
+   /// Construct the L2_HexahedronElement of order @a p and BasisType @a btype
+   L2_HexahedronElement(const int p,
+                        const int btype = BasisType::GaussLegendre);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+};
+
+
+/// Arbitrary order L2 elements in 2D on a triangle
+class L2_TriangleElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
+   mutable DenseMatrix du;
+#endif
+   DenseMatrixInverse Ti;
+
+public:
+   /// Construct the L2_TriangleElement of order @a p and BasisType @a btype
+   L2_TriangleElement(const int p,
+                      const int btype = BasisType::GaussLegendre);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   virtual void ProjectCurl(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &curl) const
+   { ProjectCurl_2D(fe, Trans, curl); }
+};
+
+
+/// Arbitrary order L2 elements in 3D on a tetrahedron
+class L2_TetrahedronElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z, shape_l;
+   mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
+   mutable DenseMatrix du;
+#endif
+   DenseMatrixInverse Ti;
+
+public:
+   /// Construct the L2_TetrahedronElement of order @a p and BasisType @a btype
+   L2_TetrahedronElement(const int p,
+                         const int btype = BasisType::GaussLegendre);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const;
+};
+
+
+/// Arbitrary order L2 elements in 3D on a wedge
+class L2_WedgeElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector t_shape, s_shape;
+   mutable DenseMatrix t_dshape, s_dshape;
+#endif
+   Array<int> t_dof, s_dof;
+
+   L2_TriangleElement TriangleFE;
+   L2_SegmentElement  SegmentFE;
+
+public:
+   /// Construct the L2_WedgeElement of order @a p and BasisType @a btype
+   L2_WedgeElement(const int p,
+                   const int btype = BasisType::GaussLegendre);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+};
+
+} // namespace mfem
+
+#endif