testing dpg residual with cholesky factors

fix repo check
Merge branch 'master' into densemat-cholesky
2022-03-22 14:13:56 -07:00 · 2022-03-22 12:08:56 -07:00 · 2022-03-22 12:05:57 -07:00 · 2022-03-22 12:05:22 -07:00 · 2022-03-22 10:24:38 -07:00 · 2022-03-22 10:19:35 -07:00
825 changed files with 39711 additions and 7444 deletions
@@ -5,7 +5,7 @@ channels:
 dependencies:
 - xeus-cling=0.13.0
 - xwidgets=0.26.0
-# NOTE: it's possible these aren't needed for the lab frontend
+# NOTE: it is possible that these are not needed for the lab frontend
 - widgetsnbextension=3.5.1
 - pip
 - pip:
@@ -0,0 +1,48 @@
+                    Finite Element Discretization Library
+                                   __
+                       _ __ ___   / _|  ___  _ __ ___
+                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
+                      | | | | | ||  _||  __/| | | | | |
+                      |_| |_| |_||_|   \___||_| |_| |_|
+
+                               https://mfem.org
+
+
+This directory contains the GitHub CI scripts for MFEM.
+
+Note that some of these scripts use the shared MFEM GitHub Actions from the external mfem/github-actions repository:
+
+  https://github.com/mfem/github-actions
+
+For a particular action, e.g. `mfem/github-actions/build-mfem@v2.1`, the `v2.1` suffix denotes the branch in the above from which the action is taken.
+
+The current CI workflows are:
+
+### `repo-check.yml`
+
+Runs a number of static repository-level sanity checks.
+
+- `file-headers-check` checks copyright date, license, etc. using the `--copyright`, `--license` and `--release` options of the `config/githooks/pre-push` script.
+
+- `code-style` checks the code style using the `--style` option of the `config/githooks/pre-push` script.
+
+- `documentation` checks the documentation build using the `tests/scripts/documentation` script.
+
+- `branch-history` guards against accidental commits of large files using the `--history` option of the `config/githooks/pre-push` script.
+
+### `mfem-analysis.yml` (`build-analysis`)
+
+Checks if the code builds and satisfies minimal requirements.
+
+- `gitignore` builds hypre, METIS, and MFEM using `mfem/github-actions/build-hypre`, `mfem/github-actions/build-metis`, and `mfem/github-actions/build-mfem` and checks for correct `.gitignore` settings by running the `tests/scripts/gitignore` script.
+
+### `builds-and-tests.yml`
+
+Runs a matrix of builds and tests runs with different compilers, OS, mfem/hypre settings, etc. Also processes and upload Codecov reports.
+
+Uses the following GitHub Actions from https://github.com/mfem/github-actions:
+
+- `mfem/github-actions/build-hypre`
+- `mfem/github-actions/build-metis`
+- `mfem/github-actions/build-mfem`
+- `mfem/github-actions/upload-coverage`
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -166,7 +166,7 @@ jobs:

    # MFEM build and test
    - name: build
-      uses: mfem/github-actions/build-mfem@v2.0
+      uses: mfem/github-actions/build-mfem@v2.1
      with:
        os: ${{ matrix.os }}
        target: ${{ matrix.target }}
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -51,8 +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/ex3[0-9]
+examples/ex3[0-9]p

 examples/refined.mesh
 examples/displaced.mesh
@@ -78,6 +78,7 @@ examples/deformed.*
 examples/velocity.*
 examples/elastic_energy.*
 examples/mode_*
+examples/mode_deriv_*
 examples/ex5-p-*.bp
 examples/ex9-p-*.bp
 examples/ex12-p-*.bp
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -18,14 +18,14 @@ variables:
    - shell
    - corona
  rules:
-    # Don’t run corona jobs if...
+    # 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...
+    # Don't run autotest update if...
    - if: '$CI_JOB_NAME =~ /report/ && $AUTOTEST != "YES"'
      when: never
    # Report success on success status
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -20,8 +20,8 @@
    - |
      (
        date
-        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
-        # try to get an excusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
+        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
        # every 5 seconds; simply using no timeout, i.e. 'flock 9', causes the
        # command to hang indefinitely sometimes, so we use the timeout & retry
        # as a workaround; we may want to add a counter for the number of
@@ -29,7 +29,7 @@
        while ! flock -w 5 9; do
          true
        done
-        echo "Aquired lock on '$PWD/autotest.lock'"
+        echo "Acquired lock on '$PWD/autotest.lock'"
        date
        # Report SUCCESS while holding the file lock on 'autotest.lock'.
        # The next script uses the following environment variables:
@@ -55,8 +55,8 @@
    - |
      (
        date
-        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
-        # try to get an excusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
+        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
        # every 5 seconds; simply using no timeout, i.e. 'flock 9', causes the
        # command to hang indefinitely sometimes, so we use the timeout & retry
        # as a workaround; we may want to add a counter for the number of
@@ -64,7 +64,7 @@
        while ! flock -w 5 9; do
          true
        done
-        echo "Aquired lock on '$PWD/autotest.lock'"
+        echo "Acquired lock on '$PWD/autotest.lock'"
        date
        # Report FAILURE while holding the file lock on 'autotest.lock'.
        # The next script uses the following environment variables:
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -45,8 +45,8 @@ setup_baseline:
    - |
      (
        date
-        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
-        # try to get an excusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
+        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
        # every 5 seconds; simply using no timeout, i.e. 'flock 9', causes the
        # command to hang indefinitely sometimes, so we use the timeout & retry
        # as a workaround; we may want to add a counter for the number of
@@ -54,7 +54,7 @@ setup_baseline:
        while ! flock -w 5 9; do
          true
        done
-        echo "Aquired lock on '$PWD/autotest.lock'"
+        echo "Acquired lock on '$PWD/autotest.lock'"
        date
        # clone/update the autotest repo while holding the file lock on
        # 'autotest.lock'
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -34,8 +34,8 @@ setup:
    - |
      (
        date
-        echo "Waiting to aquire lock on '$PWD/mfem-data.lock' ..."
-        # try to get an excusive lock on fd 9 (mfem-data.lock) repeating the try
+        echo "Waiting to acquire lock on '$PWD/mfem-data.lock' ..."
+        # try to get an exclusive lock on fd 9 (mfem-data.lock) repeating the try
        # every 5 seconds; simply using no timeout, i.e. 'flock 9', causes the
        # command to hang indefinitely sometimes, so we use the timeout & retry
        # as a workaround; we may want to add a counter for the number of
@@ -43,7 +43,7 @@ setup:
        while ! flock -w 5 9; do
          true
        done
-        echo "Aquired lock on '$PWD/mfem-data.lock'"
+        echo "Acquired lock on '$PWD/mfem-data.lock'"
        date
        # clone/update the mfem/data repo while holding the file lock on
        # 'mfem-data.lock'
@@ -67,8 +67,8 @@ setup:
    - |
      (
        date
-        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
-        # try to get an excusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
+        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
        # every 5 seconds; simply using no timeout, i.e. 'flock 9', causes the
        # command to hang indefinitely sometimes, so we use the timeout & retry
        # as a workaround; we may want to add a counter for the number of
@@ -76,7 +76,7 @@ setup:
        while ! flock -w 5 9; do
          true
        done
-        echo "Aquired lock on '$PWD/autotest.lock'"
+        echo "Acquired lock on '$PWD/autotest.lock'"
        date
        # clone/update the autotest repo while holding the file lock on
        # 'autotest.lock'
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -64,8 +64,8 @@ report_baseline:
    - |
      (
        date
-        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
-        # try to get an excusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
+        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
        # every 5 seconds; simply using no timeout, i.e. 'flock 9', causes the
        # command to hang indefinitely sometimes, so we use the timeout & retry
        # as a workaround; we may want to add a counter for the number of
@@ -73,7 +73,7 @@ report_baseline:
        while ! flock -w 5 9; do
          true
        done
-        echo "Aquired lock on '$PWD/autotest.lock'"
+        echo "Acquired lock on '$PWD/autotest.lock'"
        date
        # ----------------------
        cd ${AUTOTEST_ROOT}/autotest || \
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -8,98 +8,145 @@
                               https://mfem.org


-Version 4.3.1 (development)
-===========================
+Version 4.4, released on March 21, 2022
+=======================================

- Add hipSPARSE support for sparse mat-vec multiplications.
+Linear and nonlinear solvers
+----------------------------
+- Added support for using the hypre library built with HIP support. Similar to
+  the existing hypre + CUDA support, most of the MFEM examples and miniapps work
+  transparently with hypre + HIP builds. This includes the BoomerAMG, AMS, and
+  ADS solvers.

- Added support for using the HYPRE library built with HIP support. Similar to
-  the HYPRE + CUDA support added earlier, most of the MFEM examples and miniapps
-  work transparently with HYPRE + HIP builds. This includes the BoomerAMG, AMS,
-  and ADS solvers.
+- Added a simple singleton class, Hypre, to automatically set hypre's global
+  parameters, particularly GPU-relevant options. Updated parallel example codes
+  and miniapps to call Hypre::Init() where appropriate.

- More explicit and consistent formating of the output of iterative solvers
+- Added hipSPARSE support for sparse matrix-vector multiplications.
+
+- More explicit and consistent formatting of the output of iterative solvers
  with the new IterativeSolver::PrintLevel options. See linalg/solvers.hpp.

- Added a miniapp for PDE-based extrapolation of finite element functions. See
-  miniapps/shifted/extrapolate.cpp.
-
- Added support for automatic differentiation. Users can select between native
-  implementation and external library implementation during configuration. One
-  parallel and two serial examples are implemented in the miniapps/autodiff/
-  directory.
-
- GridFunctionCoefficient (and the related vector, gradient, divergence, and
-  curl classes) now work properly with LORDiscretization and LORSolver.
-
- 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.
+Meshing improvements
+--------------------
+- New TMOP-based methods for hr-adaptivity, interface fitting, and tangential
+  relaxation of high-order meshes.

 - 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.
+- Added a simpler interface to access mesh face information, see FaceInformation
+  and GetFaceInformation in the Mesh class.

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

- 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 ParMesh adjacency set (adjset) creation support to the Conduit Mesh
+  Blueprint MFEM wrapper functions in ConduitDataCollection.

- Added new miniapps that use the ParELAG library, its hybrid smoothers, and the
+Discretization improvements
+---------------------------
+- Added general dof transformation to support high order Nedelec basis functions
+  on tetrahedral meshes without reordering. The ReorientTetMesh method of the
+  Mesh and ParMesh classes has been deprecated. See the new DofTransformation
+  class in fem/doftrans.hpp.
+
+- GPU-enabled partial (PA) and element (EA) assembly for discontinuous Galerkin
+  methods on nonconforming AMR meshes.
+
+- Support for arbitrary order Nedelec and Raviart-Thomas elements on wedges.
+
+- Added special Nedelec and Raviart-Thomas basis functions for modeling three
+  dimensional vector fields in 1D and 2D domains, see the new Example 31/31p.
+
+- GridFunctionCoefficient (and the related vector, gradient, divergence, and
+  curl classes) now work properly with LORDiscretization and LORSolver.
+
+- Added PA support for the action of MixedScalarCurlIntegrator in 2D and
+  MixedVectorGradientIntegrator in 2D and 3D, as well as their transposes.
+
+- Coefficient::SetTime now propagates the new time into internally stored
+  Coefficient objects.
+
+- Split the fem/fe.?pp files into separate files in the new fem/fe/ directory to
+  simplify and clarify the organization of FiniteElement classes.
+
+New and updated examples and miniapps
+-------------------------------------
+- Added two new miniapps with initial support for automatic differentiation (AD)
+  in the miniapps/autodiff/ directory. Users can select between external library
+  and native implementation during configuration. The support for AD will be
+  extended in future releases of MFEM.
+
+- Added Binder (mybinder.org) configuration files for C++ MFEM Jupyter Notebooks
+  with inline GLVis visualization in the new examples/jupyter/ directory with a
+  sample notebook based on Example 1. The implementation is based on xeus-cling,
+  see github.com/jupyter-xeus/xeus-cling and github.com/GLVis/xeus-glvis.
+
+- Added a new miniapp (Extrapolation) for PDE-based extrapolation of finite
+  element functions from known values in a set of elements to the rest of the
+  computational domain. See miniapps/shifted/extrapolate.cpp.
+
+- Added new miniapp that uses 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.

+- Added a new Example 30/30p demonstrating 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.
+
+- Added new Examples 31, 31p and 32p showing anisotropic definite Maxwell
+  serial/parallel solver and parallel eigensolver 1D, 2D, or 3D.
+
+- Updated the mesh-optimizer and pmesh-optimizer miniapps to demonstrate the
+  hr-adaptivity and interface fitting capability.
+
+- The HPC versions of ex1 and ex1p (in miniapps/performance) now support runtime
+  selection of either 2D or 3D meshes.
+
+Integrations, testing and documentation
+---------------------------------------
+- Doxygen documentation for all releases is now available at docs.mfem.org.
+
+- The following integrations have updated minimum version requirements:
+  * HIOP >= 0.4.6
+  * HYPRE >= 2.23.0 for HIP support
+  * libCEED >= 0.10
+  * PUMI >= 2.2.6
+  * RAJA >= 0.14.0
+  * Umpire >= 3.0.0
+  see INSTALL for more details.
+
+- Added new optional integrations with ParELAG and CoDiPack (version >= 1.9.3+).
+
+- Added initial support for Google Benchmark (version >= 1.5.6) in the
+  tests/benchmarks directory. It can be enabled with MFEM_USE_BENCHMARK=YES.
+
+- Switched from Artistic Style (astyle) version 2.05.1 to version 3.1 for code
+  formatting. See the "make style" target.
+
+Miscellaneous
+-------------
+- Added a simple singleton class, Mpi, as a replacement for MPI_Session. New
+  code should use Mpi::Init() and other Mpi methods instead of MPI_Session.
+
+- Added ParaView visualization of QuadratureFunction fields, through both
+  QuadratureFunction::SaveVTU and ParaViewDataCollection::RegisterQField.
+
 - Fixed several MinGW build issues on Windows.

- Remove the 'u' flag in the ar command, to update all files in the archive,
+- Added 'double' atomicAdd implementation for previous versions of CUDA.
+
+- HypreParVector and Vector now support C++ move semantics, and the copy
+  constructor for HypreParVector now copies the local vector data.
+
+- Removed 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.
-
- Added ParMesh Adjaceny Set (adjset) creation support to the Conduit Mesh
-  Blueprint MFEM wrapper functions in ConduitDataCollection.
-
- `HypreParVector` and `Vector` now support move semantics, and the copy
-  constructor for `HypreParVector` now copies the local vector data.
-
- The HPC versions of ex1 and ex1p (in miniapps/performance) now support
-  runtime selection of either 2D or 3D meshes.
-
- Added ParaView visualization of `QuadratureFunction` fields, through both
-  `QuadratureFunction::SaveVTU` and `ParaViewDataCollection::RegisterQField`.
+- Various other simplifications, extensions, and bugfixes in the code.


 Version 4.3, released on July 29, 2021
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -51,7 +51,7 @@ project(mfem NONE)
 # Current version of MFEM, see also `makefile`.
 #   mfem_VERSION = (string)
 #   MFEM_VERSION = (int)   [automatically derived from mfem_VERSION]
-set(${PROJECT_NAME}_VERSION 4.3.1)
+set(${PROJECT_NAME}_VERSION 4.4.0)

 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -8,7 +8,7 @@
 <a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amaster"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=master"></a>
 <a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amaster"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=master"></a>
 <a href="https://ci.appveyor.com/project/mfem/mfem"><img alt="Build Status" src="https://ci.appveyor.com/api/projects/status/19non9sqm6msi2wy?svg=true"></a>
-<a href="https://mfem.github.io/doxygen/html/index.html"><img alt="Doxygen" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
+<a href="https://docs.mfem.org/html/index.html"><img alt="Doxygen" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
 </p>


@@ -141,6 +141,7 @@ The MFEM source code has the following structure:
  │   ├── tools
  │   └── toys
  └── tests
+      ├── benchmarks
      ├── convergence
      ├── gitlab
      ├── mem_manager
@@ -156,24 +157,24 @@ classes implementing the finite element, mesh and linear algebra concepts
 respectively.

 - The main mesh classes are:
-  + [`Mesh`](https://mfem.github.io/doxygen/html/classmfem_1_1Mesh.html)
-  + [`NCMesh`](https://mfem.github.io/doxygen/html/classmfem_1_1NCMesh.html)
-  + [`Element`](https://mfem.github.io/doxygen/html/classmfem_1_1Element.html)
-  + [`ElementTransformation`](https://mfem.github.io/doxygen/html/classmfem_1_1ElementTransformation.html)
+  + [`Mesh`](https://docs.mfem.org/html/classmfem_1_1Mesh.html)
+  + [`NCMesh`](https://docs.mfem.org/html/classmfem_1_1NCMesh.html)
+  + [`Element`](https://docs.mfem.org/html/classmfem_1_1Element.html)
+  + [`ElementTransformation`](https://docs.mfem.org/html/classmfem_1_1ElementTransformation.html)

 - The main finite element classes are:
-  + [`FiniteElement`](https://mfem.github.io/doxygen/html/classmfem_1_1FiniteElement.html)
-  + [`FiniteElementCollection`](https://mfem.github.io/doxygen/html/classmfem_1_1FiniteElement.html)
-  + [`FiniteElementSpace`](https://mfem.github.io/doxygen/html/classmfem_1_1FiniteElementSpace.html)
-  + [`GridFunction`](https://mfem.github.io/doxygen/html/classmfem_1_1GridFunction.html)
-  + [`BilinearFormIntegrator`](https://mfem.github.io/doxygen/html/classmfem_1_1BilinearFormIntegrator.html) and [`LinearFormIntegrator`](https://mfem.github.io/doxygen/html/classmfem_1_1LinearFormIntegrator.html)
-  + [`LinearForm`](https://mfem.github.io/doxygen/html/classmfem_1_1LinearFormIntegrator.html), [`BilinearForm`](https://mfem.github.io/doxygen/html/classmfem_1_1BilinearForm.html) and [`MixedBilinearForm`](https://mfem.github.io/doxygen/html/classmfem_1_1MixedBilinearForm.html)
+  + [`FiniteElement`](https://docs.mfem.org/html/classmfem_1_1FiniteElement.html)
+  + [`FiniteElementCollection`](https://docs.mfem.org/html/classmfem_1_1FiniteElement.html)
+  + [`FiniteElementSpace`](https://docs.mfem.org/html/classmfem_1_1FiniteElementSpace.html)
+  + [`GridFunction`](https://docs.mfem.org/html/classmfem_1_1GridFunction.html)
+  + [`BilinearFormIntegrator`](https://docs.mfem.org/html/classmfem_1_1BilinearFormIntegrator.html) and [`LinearFormIntegrator`](https://docs.mfem.org/html/classmfem_1_1LinearFormIntegrator.html)
+  + [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearFormIntegrator.html), [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`MixedBilinearForm`](https://docs.mfem.org/html/classmfem_1_1MixedBilinearForm.html)

 - The main linear algebra classes and sources are
-  + [`Operator`](https://mfem.github.io/doxygen/html/classmfem_1_1Operator.html) and [`BilinearForm`](https://mfem.github.io/doxygen/html/classmfem_1_1BilinearForm.html)
-  + [`Vector`](https://mfem.github.io/doxygen/html/classmfem_1_1BilinearForm.html) and [`LinearForm`](https://mfem.github.io/doxygen/html/classmfem_1_1LinearForm.html)
-  + [`DenseMatrix`](https://mfem.github.io/doxygen/html/classmfem_1_1DenseMatrix.html) and [`SparseMatrix`](https://mfem.github.io/doxygen/html/classmfem_1_1SparseMatrix.html)
-  + Sparse [smoothers](https://mfem.github.io/doxygen/html/sparsesmoothers_8hpp.html) and linear [solvers](https://mfem.github.io/doxygen/html/solvers_8hpp.html)
+  + [`Operator`](https://docs.mfem.org/html/classmfem_1_1Operator.html) and [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html)
+  + [`Vector`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearForm.html)
+  + [`DenseMatrix`](https://docs.mfem.org/html/classmfem_1_1DenseMatrix.html) and [`SparseMatrix`](https://docs.mfem.org/html/classmfem_1_1SparseMatrix.html)
+  + Sparse [smoothers](https://docs.mfem.org/html/sparsesmoothers_8hpp.html) and linear [solvers](https://docs.mfem.org/html/solvers_8hpp.html)

 #### Parallel implementation

@@ -183,13 +184,13 @@ shared geometric entities between different tasks. The parallel source files
 have a `p` prefix, e.g. `pmesh.cpp` vs. the serial `mesh.cpp`.

 - The main parallel classes are
-  + [`ParMesh`](https://mfem.github.io/doxygen/html/solvers_8hpp.html)
-  + [`ParNCMesh`](https://mfem.github.io/doxygen/html/classmfem_1_1ParMesh.html)
-  + [`ParFiniteElementSpace`](https://mfem.github.io/doxygen/html/classmfem_1_1ParFiniteElementSpace.html)
-  + [`ParGridFunction`](https://mfem.github.io/doxygen/html/classmfem_1_1ParGridFunction.html)
-  + [`ParBilinearForm`](https://mfem.github.io/doxygen/html/classmfem_1_1ParBilinearForm.html) and [`ParLinearForm`](https://mfem.github.io/doxygen/html/classmfem_1_1ParLinearForm.html)
-  + [`HypreParMatrix`](https://mfem.github.io/doxygen/html/classmfem_1_1HypreParMatrix.html) and [`HypreParVector`](https://mfem.github.io/doxygen/html/classmfem_1_1HypreParVector.html)
-  + [`HypreSolver`](https://mfem.github.io/doxygen/html/classmfem_1_1HypreSolver.html) and other [hypre classes](https://mfem.github.io/doxygen/html/hypre_8hpp.html)
+  + [`ParMesh`](https://docs.mfem.org/html/solvers_8hpp.html)
+  + [`ParNCMesh`](https://docs.mfem.org/html/classmfem_1_1ParMesh.html)
+  + [`ParFiniteElementSpace`](https://docs.mfem.org/html/classmfem_1_1ParFiniteElementSpace.html)
+  + [`ParGridFunction`](https://docs.mfem.org/html/classmfem_1_1ParGridFunction.html)
+  + [`ParBilinearForm`](https://docs.mfem.org/html/classmfem_1_1ParBilinearForm.html) and [`ParLinearForm`](https://docs.mfem.org/html/classmfem_1_1ParLinearForm.html)
+  + [`HypreParMatrix`](https://docs.mfem.org/html/classmfem_1_1HypreParMatrix.html) and [`HypreParVector`](https://docs.mfem.org/html/classmfem_1_1HypreParVector.html)
+  + [`HypreSolver`](https://docs.mfem.org/html/classmfem_1_1HypreSolver.html) and other [hypre classes](https://docs.mfem.org/html/hypre_8hpp.html)

 #### GPU and general device support

@@ -198,10 +199,10 @@ backends (CUDA, OCCA, RAJA, OpenMP, etc.) and an internal lightweight
 device/host memory manager.

 - The main device-relevant classes and sources are:
-  + [`Device`](https://mfem.github.io/doxygen/html/device_8hpp.html)
-  + [`MemoryManager`](https://mfem.github.io/doxygen/html/mem_manager_8hpp.html)
-  + the [`MFEM_FORALL`](https://mfem.github.io/doxygen/html/forall_8hpp.html) macro
-  + the [`cuda.hpp`](https://mfem.github.io/doxygen/html/cuda_8hpp.html) and [`occa.hpp`](https://mfem.github.io/doxygen/html/occa_8hpp.html) files
+  + [`Device`](https://docs.mfem.org/html/device_8hpp.html)
+  + [`MemoryManager`](https://docs.mfem.org/html/mem_manager_8hpp.html)
+  + the [`MFEM_FORALL`](https://docs.mfem.org/html/forall_8hpp.html) macro
+  + the [`cuda.hpp`](https://docs.mfem.org/html/cuda_8hpp.html) and [`occa.hpp`](https://docs.mfem.org/html/occa_8hpp.html) files

 #### Utilities, building and documentation
 - The `general/` directory contains C++ classes that serve as utilities for
@@ -209,8 +210,7 @@ device/host memory manager.
 - The `config/` directory contains build-related files, both for the plain
  Makefile and the CMake build options.
 - The `doc/` directory contains configuration for the Doxygen code documentation
-  that can either be built locally or browsed online at
-  https://mfem.github.io/doxygen/html/index.html.
+  that can either be built locally or browsed online at https://docs.mfem.org.

 #### Examples and tests
 - `examples` and `miniapps` respectively gather simple and more fully-featured
@@ -430,10 +430,10 @@ The Pull Request (PR) approval process in MFEM is similar to the approval of pap

 The current list of MFEM editors is:

- @v-dobrev (Veselin Dobrev)
- @tzanio (Tzanio Kolev)
- @pazner (Will Pazner)
- @mlstowell (Mark Stowell)
+- [@v-dobrev](https://github.com/v-dobrev) (Veselin Dobrev)
+- [@tzanio](https://github.com/tzanio) (Tzanio Kolev)
+- [@pazner](https://github.com/pazner) (Will Pazner)
+- [@mlstowell](https://github.com/mlstowell) (Mark Stowell)

 **The responsibilities of the editors are:**

@@ -639,7 +639,7 @@ MFEM uses a `master`/`next`-branch workflow as described below:
  ```
 - [ ] Create the release tarball and push to `mfem/releases`.
 - [ ] Recreate the `next` branch as described in previous section.
- [ ] Update and push documentation  to `mfem/doxygen`.
+- [ ] Update and push documentation to `mfem/doxygen`. Update the `README.md` file and the `html` link in the `mfem/doxygen` repo.
 - [ ] Update URL shortlinks:
    - [ ] Create a shortlink at [http://bit.ly/](http://bit.ly/) for the release tarball, e.g. https://mfem.github.io/releases/mfem-3.1.tgz.
    - [ ] (LLNL only) Add and commit the new shortlink in the `links` and `links-mfem` files of the internal `mfem/downloads` repo.
@@ -19,10 +19,10 @@ requires an MPI C++ compiler, as well as the following external libraries:
  http://glaros.dtc.umn.edu/gkhome/metis/metis/overview

 The hypre dependency can be downloaded as a tarball from GitHub or from the
-project webpage https://www.llnl.gov/casc/hypre. For example, the 2.20.0 release
+project webpage https://www.llnl.gov/casc/hypre. For example, the 2.24.0 release
 of hypre is available at

-   https://github.com/hypre-space/hypre/archive/v2.20.0.tar.gz
+   https://github.com/hypre-space/hypre/archive/v2.24.0.tar.gz

 The METIS dependency can be disabled but that is not generally recommended, see
 the option MFEM_USE_METIS.
@@ -60,7 +60,6 @@ following package managers:
 - OpenHPC, http://openhpc.community
 - Conda-forge, https://conda-forge.org (pre-built binaries linked with
  OpenMPI/MPICH, hypre, and METIS)
- Homebrew/Science, https://github.com/Homebrew/homebrew-science (deprecated)

 We also recommend downloading and building the MFEM-based GLVis visualization
 tool which can be used to visualize the meshes and solution in MFEM's examples
@@ -567,7 +566,8 @@ 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.1   (HYPRE built with CUDA or HIP)
+            HYPRE >= 2.22.1   (HYPRE built with CUDA)
+            HYPRE >= 2.23.0   (HYPRE built with HIP)

 - METIS, used when MFEM_USE_METIS = YES. If using METIS 5, set
  MFEM_USE_METIS_5 = YES (default is to use METIS 4).
@@ -713,7 +713,7 @@ The specific libraries and their options are:
  Options: HIOP_OPT, HIOP_LIB.
  Versions: HIOP >= 0.4.6.

- CoDiPack (optiobal), used with MFEM_USE_CODIPACK = YES
+- CoDiPack (optional), used with MFEM_USE_CODIPACK = YES
  URL: https://www.scicomp.uni-kl.de/codi/
  Options: CODIPACK_OPT
  Versions: 1.9.3
@@ -751,7 +751,7 @@ The specific libraries and their options are:
  URL: https://github.com/CEED/libCEED
       https://ceed.exascaleproject.org/libceed
  Options: CEED_DIR, CEED_OPT, CEED_LIB.
-  Versions: libCEED >= 0.8.
+  Versions: libCEED >= 0.10.

 - RAJA (optional), used when MFEM_USE_RAJA = YES.
  Beginning with MFEM v4.3, only RAJA v0.14.0+ is supported.
@@ -1,6 +1,6 @@
 BSD 3-Clause License

-Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC
+Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC
 All rights reserved.

 Redistribution and use in source and binary forms, with or without
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -35,8 +35,8 @@ cd $(git rev-parse --show-toplevel)
 # copyright check
 copyright=true
 if [[ "${option}" == "--copyright" || "${option}" == "" ]]; then
-  if git grep -l "^\(#\|//\).*\(\-2020\|\ 2010,\)" > matches.txt; then
-    echo "Please update the following files to Copyright (c) 2010-2021:"
+  if git grep -l "^\(#\|//\).*\(\-2021\|\ 2010,\)" > matches.txt; then
+    echo "Please update the following files to Copyright (c) 2010-2022:"
    cat matches.txt
    copyright=false
  fi
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -32,7 +32,7 @@ groups_serial=(
 '"examples"
   "Examples:"
   "examples"
-   "ex{,1,2}[0-9].cpp"'
+   "ex{,1,2,3}[0-9].cpp"'
 #   "ex1.cpp"'
 '"sundials"
   "SUNDIALS examples:"
@@ -83,11 +83,13 @@ groups_serial=(
 '"tools"
   "Tools miniapps:"
   "miniapps/tools"
-   "convert-dc.cpp display-basis.cpp get-values.cpp load-dc.cpp lor-transfer.cpp"'
+   "convert-dc.cpp display-basis.cpp get-values.cpp load-dc.cpp
+    lor-transfer.cpp"'
 '"toys"
   "Toys miniapps:"
   "miniapps/toys"
-   "automata.cpp life.cpp lissajous.cpp mandel.cpp mondrian.cpp rubik.cpp snake.cpp"'
+   "automata.cpp life.cpp lissajous.cpp mandel.cpp mondrian.cpp rubik.cpp
+    snake.cpp"'
 '"convergence"
   "Convergence tests:"
   "tests/convergence"
@@ -98,7 +100,7 @@ groups_parallel=(
 '"examples"
   "Examples:"
   "examples"
-   "ex{,1,2}[0-9]p.cpp"'
+   "ex{,1,2,3}[0-9]p.cpp"'
 #   "ex1p.cpp"'
 '"sundials"
   "SUNDIALS examples:"
@@ -184,7 +186,7 @@ groups_all=(
 '"examples"
   "Examples:"
   "examples"
-   "ex\"{,1,2}[0-9]\"{,p}.cpp"'
+   "ex\"{,1,2,3}[0-9]\"{,p}.cpp"'
 '"sundials"
   "SUNDIALS examples:"
   "examples/sundials"
@@ -238,7 +240,8 @@ groups_all=(
 '"gslib"
   "GSLIB miniapps:"
   "miniapps/gslib"
-   "field-diff.cpp field-interp.cpp findpts.cpp schwarz_ex1.cpp pfindpts.cpp schwarz_ex1p.cpp"'
+   "field-diff.cpp field-interp.cpp findpts.cpp schwarz_ex1.cpp pfindpts.cpp
+    schwarz_ex1p.cpp"'
 '"navier"
   "Navier miniapps:"
   "miniapps/navier"
@@ -258,11 +261,13 @@ groups_all=(
 '"tools"
   "Tools miniapps:"
   "miniapps/tools"
-   "convert-dc.cpp display-basis.cpp get-values.cpp load-dc.cpp lor-transfer.cpp"'
+   "convert-dc.cpp display-basis.cpp get-values.cpp load-dc.cpp
+    lor-transfer.cpp"'
 '"toys"
   "Toys miniapps:"
   "miniapps/toys"
-   "automata.cpp life.cpp lissajous.cpp mandel.cpp mondrian.cpp rubik.cpp snake.cpp"'
+   "automata.cpp life.cpp lissajous.cpp mandel.cpp mondrian.cpp rubik.cpp
+    snake.cpp"'
 '"convergence"
   "Convergence tests:"
   "tests/convergence"
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -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
+#
+
+dimension
+2
+
+# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
+elements
+3
+
+1 3 0 1 2 3
+2 3 0 3 4 5
+3 3 0 5 6 1
+
+boundary
+6
+1 1 1 2
+2 1 2 3
+3 1 3 4
+4 1 4 5
+5 1 5 6
+6 1 6 1
+
+vertices
+7
+2
+0.0 0.0
+-0.50 -0.8660254037844386
+0.50 -0.8660254037844386
+1.0 0.0
+0.50 0.8660254037844386
+-0.50 0.8660254037844386
+-1.0 0.0
@@ -15,146 +15,146 @@ DataCollection/Fields/0/NumberOfVectorComponents: 2
 DataCollection/Fields/0/FieldDescriptorName: CoordsDescriptor
 DataCollection/Fields/0/Data/Size: 422
 DataCollection/Fields/0/Data/Type: FMS_DOUBLE
-DataCollection/Fields/0/Data/Values: [-0.016886, 1.000000, 0.309017, 
-1.309020, -0.809017, -0.500000, 
-0.809017, -1.618030, 0.309017, 
-0.500000, 1.309020, 0.519420, 
-1.154510, 0.809019, 0.147680, 
-0.095492, -0.654508, -0.415586, 
-1.213520, -1.213520, -0.392210, 
-0.654508, -0.095492, 0.139949, 
-0.809019, 1.154510, 0.660184, 
-0.264063, -0.800064, -0.231060, 
-0.663691, 0.183114, 0.317639, 
-0.543082, 0.598483, 0.345112, 
-0.478298, 0.027703, 0.095229, 
-0.012368, -0.092534, -0.334412, 
-0.313767, -0.140526, -0.293881, 
-0.534056, -0.660290, -0.537646, 
-0.655590, -0.121396, -0.274504, 
-0.346497, -0.296570, 0.004737, 
-0.098835, 0.069287, 0.082675, 
-0.318799, 0.467183, 0.564505, 
-0.595190, 0.846237, 0.671735, 
-1.051500, 1.103010, 0.964008, 
-0.821603, 1.257520, 1.206010, 
-1.142350, 0.975686, 0.781273, 
-0.717257, 0.475684, 0.642352, 
-0.268930, 0.211049, 0.174181, 
-0.039345, -0.147746, -0.177481, 
-0.365164, -0.230328, -0.551503, 
-0.603005, -0.497587, -0.389864, 
-0.757514, -0.706011, -0.675487, 
-0.528946, -0.943851, -1.078690, 
-1.087600, -0.955467, -1.483190, 
-1.348360, -1.483190, -1.348360, 
-1.085930, -0.938010, -0.943851, 
-1.078690, -0.681476, -0.540944, 
-0.757514, -0.706011, -0.540614, 
-0.367058, -0.551503, -0.603005, 
-0.365164, -0.230328, -0.138552, 
-0.206896, 0.174181, 0.039345, 
-0.268468, 0.222269, 0.475684, 
-0.642352, 0.759791, 0.719381, 
-1.142350, 0.975686, 1.257520, 
-1.206010, 0.972837, 0.836119, 
-1.051500, 1.103010, 0.214572, 
-0.407449, 0.288323, 0.449827, 
-0.086700, -0.027358, -0.200560, 
-0.166595, -0.271802, -0.418426, 
-0.426131, -0.551441, -0.096117, 
-0.206969, -0.027946, -0.184969, 
-0.211136, 0.260131, 0.407172, 
-0.430781, 0.718277, 0.885068, 
-0.753103, 0.957692, 0.866273, 
-1.024530, 0.934099, 1.093820, 
-0.348422, 0.524463, 0.404903, 
-0.587376, 0.054525, 0.146431, 
-0.078026, -0.007795, -0.329488, 
-0.302967, -0.488115, -0.439332, 
-0.498515, -0.453527, -0.633059, 
-0.570251, -0.655787, -0.791132, 
-0.802013, -0.956872, -1.094720, 
-1.208560, -1.207720, -1.339910, 
-0.693371, -0.795404, -0.803594, 
-0.945068, -0.464668, -0.631721, 
-0.458968, -0.547876, -0.332878, 
-0.485696, -0.286385, -0.408481, 
-0.053021, -0.058200, 0.110846, 
-0.002086, 0.381384, 0.416784, 
-0.551133, 0.613261, 0.872474, 
-0.901208, 1.038300, 1.084660, 
-0.737459, 0.751250, 0.890002, 
-0.915210, 0.010915, 0.000000, 
-0.951057, 0.951057, 0.587785, 
-1.538840, -0.587785, 0.000000, 
-0.951057, -1.538840, -0.951057, 
-0.015847, 0.475529, 0.951057, 
-0.492248, 1.244950, 1.063310, 
-0.274399, 0.293893, -0.293892, 
-0.296404, -1.063310, -1.244950, 
-0.453865, -0.951057, -0.475529, 
-0.466620, 0.792932, -0.013913, 
-0.748783, -0.497528, 0.021382, 
-0.017158, 0.172591, 0.330125, 
-0.458568, 0.457971, 0.137740, 
-0.299049, 0.588394, 0.667324, 
-0.432341, 0.634346, 0.117322, 
-0.193603, 0.211702, 0.098278, 
-0.199438, -0.077304, -0.082243, 
-0.216296, -0.458634, -0.592374, 
-0.563926, -0.680404, -0.135751, 
-0.302942, -0.469005, -0.453640, 
-0.182727, -0.314240, 0.024270, 
-0.021546, 0.158510, 0.317019, 
-0.485799, 0.492951, 0.792548, 
-0.634038, 0.951057, 0.951057, 
-0.777915, 0.613430, 0.951057, 
-0.951057, 0.793994, 0.635800, 
-1.049020, 1.146990, 1.084480, 
-0.924310, 1.440880, 1.342910, 
-1.380330, 1.221820, 0.948209, 
-0.856297, 0.746293, 0.904802, 
-0.476242, 0.393234, 0.489821, 
-0.391857, 0.194471, 0.075751, 
-0.097964, 0.195929, -0.097964, 
-0.195928, -0.173234, -0.078922, 
-0.489821, -0.391856, -0.467007, 
-0.397859, -0.746293, -0.904802, 
-0.945206, -0.849559, -1.380330, 
-1.221820, -1.440880, -1.342910, 
-1.100830, -0.923191, -1.049020, 
-1.146990, -0.774515, -0.621542, 
-0.951057, -0.951057, -0.803055, 
-0.635255, -0.951057, -0.951057, 
-0.792548, -0.634038, -0.454301, 
-0.479369, -0.158510, -0.317019, 
-0.149331, 0.178643, 0.295860, 
-0.306275, 0.246225, 0.406610, 
-0.367954, 0.498458, -0.014929, 
-0.084917, -0.094272, -0.000726, 
-0.245374, -0.340755, -0.435351, 
-0.490564, -0.176355, -0.294974, 
-0.148010, -0.328246, 0.136248, 
-0.182883, 0.328957, 0.309903, 
-0.646446, 0.622546, 0.800859, 
-0.801803, 0.643087, 0.654473, 
-0.796963, 0.816799, 0.748189, 
-0.898148, 0.819092, 0.970033, 
-1.049540, 1.162660, 1.115310, 
-1.266700, 0.539960, 0.728396, 
-0.656318, 0.783152, 0.302716, 
-0.403858, 0.219006, 0.309252, 
-0.009667, 0.107161, -0.118211, 
-0.011786, -0.284704, -0.179858, 
-0.411567, -0.296676, -0.570168, 
-0.635348, -0.726784, -0.793211, 
-1.016860, -1.112970, -1.170910, 
-1.281530, -0.737860, -0.851723, 
-0.883732, -0.995070, -0.653153, 
-0.769939, -0.631918, -0.797688, 
-0.616289, -0.806819, -0.638485, 
-0.790356, -0.136399, -0.322769, 
+DataCollection/Fields/0/Data/Values: [-0.016886, 1.000000, 0.309017,
+1.309020, -0.809017, -0.500000,
+-0.809017, -1.618030, 0.309017,
+-0.500000, 1.309020, 0.519420,
+1.154510, 0.809019, 0.147680,
+-0.095492, -0.654508, -0.415586,
+-1.213520, -1.213520, -0.392210,
+-0.654508, -0.095492, 0.139949,
+0.809019, 1.154510, 0.660184,
+-0.264063, -0.800064, -0.231060,
+0.663691, 0.183114, 0.317639,
+0.543082, 0.598483, 0.345112,
+0.478298, 0.027703, 0.095229,
+0.012368, -0.092534, -0.334412,
+-0.313767, -0.140526, -0.293881,
+-0.534056, -0.660290, -0.537646,
+-0.655590, -0.121396, -0.274504,
+-0.346497, -0.296570, 0.004737,
+-0.098835, 0.069287, 0.082675,
+0.318799, 0.467183, 0.564505,
+0.595190, 0.846237, 0.671735,
+1.051500, 1.103010, 0.964008,
+0.821603, 1.257520, 1.206010,
+1.142350, 0.975686, 0.781273,
+0.717257, 0.475684, 0.642352,
+0.268930, 0.211049, 0.174181,
+0.039345, -0.147746, -0.177481,
+-0.365164, -0.230328, -0.551503,
+-0.603005, -0.497587, -0.389864,
+-0.757514, -0.706011, -0.675487,
+-0.528946, -0.943851, -1.078690,
+-1.087600, -0.955467, -1.483190,
+-1.348360, -1.483190, -1.348360,
+-1.085930, -0.938010, -0.943851,
+-1.078690, -0.681476, -0.540944,
+-0.757514, -0.706011, -0.540614,
+-0.367058, -0.551503, -0.603005,
+-0.365164, -0.230328, -0.138552,
+-0.206896, 0.174181, 0.039345,
+0.268468, 0.222269, 0.475684,
+0.642352, 0.759791, 0.719381,
+1.142350, 0.975686, 1.257520,
+1.206010, 0.972837, 0.836119,
+1.051500, 1.103010, 0.214572,
+0.407449, 0.288323, 0.449827,
+-0.086700, -0.027358, -0.200560,
+-0.166595, -0.271802, -0.418426,
+-0.426131, -0.551441, -0.096117,
+-0.206969, -0.027946, -0.184969,
+0.211136, 0.260131, 0.407172,
+0.430781, 0.718277, 0.885068,
+0.753103, 0.957692, 0.866273,
+1.024530, 0.934099, 1.093820,
+0.348422, 0.524463, 0.404903,
+0.587376, 0.054525, 0.146431,
+-0.078026, -0.007795, -0.329488,
+-0.302967, -0.488115, -0.439332,
+-0.498515, -0.453527, -0.633059,
+-0.570251, -0.655787, -0.791132,
+-0.802013, -0.956872, -1.094720,
+-1.208560, -1.207720, -1.339910,
+-0.693371, -0.795404, -0.803594,
+-0.945068, -0.464668, -0.631721,
+-0.458968, -0.547876, -0.332878,
+-0.485696, -0.286385, -0.408481,
+0.053021, -0.058200, 0.110846,
+-0.002086, 0.381384, 0.416784,
+0.551133, 0.613261, 0.872474,
+0.901208, 1.038300, 1.084660,
+0.737459, 0.751250, 0.890002,
+0.915210, 0.010915, 0.000000,
+0.951057, 0.951057, 0.587785,
+1.538840, -0.587785, 0.000000,
+-0.951057, -1.538840, -0.951057,
+-0.015847, 0.475529, 0.951057,
+0.492248, 1.244950, 1.063310,
+0.274399, 0.293893, -0.293892,
+-0.296404, -1.063310, -1.244950,
+-0.453865, -0.951057, -0.475529,
+0.466620, 0.792932, -0.013913,
+-0.748783, -0.497528, 0.021382,
+-0.017158, 0.172591, 0.330125,
+0.458568, 0.457971, 0.137740,
+0.299049, 0.588394, 0.667324,
+0.432341, 0.634346, 0.117322,
+0.193603, 0.211702, 0.098278,
+-0.199438, -0.077304, -0.082243,
+-0.216296, -0.458634, -0.592374,
+-0.563926, -0.680404, -0.135751,
+-0.302942, -0.469005, -0.453640,
+-0.182727, -0.314240, 0.024270,
+0.021546, 0.158510, 0.317019,
+0.485799, 0.492951, 0.792548,
+0.634038, 0.951057, 0.951057,
+0.777915, 0.613430, 0.951057,
+0.951057, 0.793994, 0.635800,
+1.049020, 1.146990, 1.084480,
+0.924310, 1.440880, 1.342910,
+1.380330, 1.221820, 0.948209,
+0.856297, 0.746293, 0.904802,
+0.476242, 0.393234, 0.489821,
+0.391857, 0.194471, 0.075751,
+0.097964, 0.195929, -0.097964,
+-0.195928, -0.173234, -0.078922,
+-0.489821, -0.391856, -0.467007,
+-0.397859, -0.746293, -0.904802,
+-0.945206, -0.849559, -1.380330,
+-1.221820, -1.440880, -1.342910,
+-1.100830, -0.923191, -1.049020,
+-1.146990, -0.774515, -0.621542,
+-0.951057, -0.951057, -0.803055,
+-0.635255, -0.951057, -0.951057,
+-0.792548, -0.634038, -0.454301,
+-0.479369, -0.158510, -0.317019,
+0.149331, 0.178643, 0.295860,
+0.306275, 0.246225, 0.406610,
+0.367954, 0.498458, -0.014929,
+0.084917, -0.094272, -0.000726,
+-0.245374, -0.340755, -0.435351,
+-0.490564, -0.176355, -0.294974,
+-0.148010, -0.328246, 0.136248,
+0.182883, 0.328957, 0.309903,
+0.646446, 0.622546, 0.800859,
+0.801803, 0.643087, 0.654473,
+0.796963, 0.816799, 0.748189,
+0.898148, 0.819092, 0.970033,
+1.049540, 1.162660, 1.115310,
+1.266700, 0.539960, 0.728396,
+0.656318, 0.783152, 0.302716,
+0.403858, 0.219006, 0.309252,
+0.009667, 0.107161, -0.118211,
+0.011786, -0.284704, -0.179858,
+-0.411567, -0.296676, -0.570168,
+-0.635348, -0.726784, -0.793211,
+-1.016860, -1.112970, -1.170910,
+-1.281530, -0.737860, -0.851723,
+-0.883732, -0.995070, -0.653153,
+-0.769939, -0.631918, -0.797688,
+-0.616289, -0.806819, -0.638485,
+-0.790356, -0.136399, -0.322769,
 -0.165339, -0.309622]
 DataCollection/Mesh/PartitionInfo/Size: 2
 DataCollection/Mesh/PartitionInfo/Type: FMS_UINT64
@@ -170,69 +170,69 @@ DataCollection/Mesh/DomainNames/0/Domains/0/Entities/0/EntityType: FMS_EDGE
 DataCollection/Mesh/DomainNames/0/Domains/0/Entities/0/NumEntities: 50
 DataCollection/Mesh/DomainNames/0/Domains/0/Entities/0/Size: 100
 DataCollection/Mesh/DomainNames/0/Domains/0/Entities/0/Type: FMS_INT32
-DataCollection/Mesh/DomainNames/0/Domains/0/Entities/0/Values: [11, 0, 26, 
-11, 26, 14, 
-14, 0, 27, 
-14, 27, 17, 
-17, 0, 28, 
-17, 28, 20, 
-20, 0, 29, 
-20, 29, 23, 
-23, 0, 30, 
-23, 30, 11, 
-11, 1, 12, 
-1, 26, 12, 
-12, 3, 13, 
-3, 26, 13, 
-13, 2, 14, 
-2, 15, 2, 
-27, 15, 15, 
-5, 16, 5, 
-27, 16, 16, 
-4, 17, 4, 
-18, 4, 28, 
-18, 18, 7, 
-19, 7, 28, 
-19, 19, 6, 
-20, 6, 21, 
-6, 29, 21, 
-21, 9, 22, 
-9, 29, 22, 
-22, 8, 23, 
-8, 24, 8, 
-30, 24, 24, 
-10, 25, 10, 
+DataCollection/Mesh/DomainNames/0/Domains/0/Entities/0/Values: [11, 0, 26,
+11, 26, 14,
+14, 0, 27,
+14, 27, 17,
+17, 0, 28,
+17, 28, 20,
+20, 0, 29,
+20, 29, 23,
+23, 0, 30,
+23, 30, 11,
+11, 1, 12,
+1, 26, 12,
+12, 3, 13,
+3, 26, 13,
+13, 2, 14,
+2, 15, 2,
+27, 15, 15,
+5, 16, 5,
+27, 16, 16,
+4, 17, 4,
+18, 4, 28,
+18, 18, 7,
+19, 7, 28,
+19, 19, 6,
+20, 6, 21,
+6, 29, 21,
+21, 9, 22,
+9, 29, 22,
+22, 8, 23,
+8, 24, 8,
+30, 24, 24,
+10, 25, 10,
 30, 25, 25, 1]
 DataCollection/Mesh/DomainNames/0/Domains/0/Entities/1/EntityType: FMS_QUADRILATERAL
 DataCollection/Mesh/DomainNames/0/Domains/0/Entities/1/NumEntities: 20
 DataCollection/Mesh/DomainNames/0/Domains/0/Entities/1/Size: 80
 DataCollection/Mesh/DomainNames/0/Domains/0/Entities/1/Type: FMS_INT32
-DataCollection/Mesh/DomainNames/0/Domains/0/Entities/1/Values: [0, 1, 2, 
-3, 3, 4, 
-5, 6, 6, 
-7, 8, 9, 
-9, 10, 11, 
-12, 12, 13, 
-14, 0, 15, 
-16, 17, 1, 
-17, 18, 19, 
-20, 2, 20, 
-21, 22, 22, 
-23, 24, 4, 
-24, 25, 26, 
-27, 5, 27, 
-28, 29, 29, 
-30, 31, 7, 
-31, 32, 33, 
-34, 8, 34, 
-35, 36, 36, 
-37, 38, 10, 
-38, 39, 40, 
-41, 11, 41, 
-42, 43, 43, 
-44, 45, 13, 
-45, 46, 47, 
-48, 14, 48, 
+DataCollection/Mesh/DomainNames/0/Domains/0/Entities/1/Values: [0, 1, 2,
+3, 3, 4,
+5, 6, 6,
+7, 8, 9,
+9, 10, 11,
+12, 12, 13,
+14, 0, 15,
+16, 17, 1,
+17, 18, 19,
+20, 2, 20,
+21, 22, 22,
+23, 24, 4,
+24, 25, 26,
+27, 5, 27,
+28, 29, 29,
+30, 31, 7,
+31, 32, 33,
+34, 8, 34,
+35, 36, 36,
+37, 38, 10,
+38, 39, 40,
+41, 11, 41,
+42, 43, 43,
+44, 45, 13,
+45, 46, 47,
+48, 14, 48,
 49, 15]
 DataCollection/Mesh/Components/0/Name: volume
 DataCollection/Mesh/Components/0/Dimension: 2
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -38,7 +38,7 @@ PROJECT_NAME           = "MFEM"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.

-PROJECT_NUMBER         = v4.3.1
+PROJECT_NUMBER         = v4.4.0

 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
@@ -98,6 +98,11 @@ namespace mfem {
 * - <a class="el" href="ex28p_8cpp_source.html">Example 28p</a>: parallel sliding contact in elasticity
 * - <a class="el" href="ex29_8cpp_source.html">Example 29</a>: Laplace solve on a 3D-embedded surface
 * - <a class="el" href="ex29p_8cpp_source.html">Example 29p</a>: parallel Laplace solve on a 3D-embedded surface
+ * - <a class="el" href="ex30_8cpp_source.html">Example 30</a>: mesh preprocessing to resolve problem data
+ * - <a class="el" href="ex30p_8cpp_source.html">Example 30p</a>: parallel mesh preprocessing to resolve problem data
+ * - <a class="el" href="ex31_8cpp_source.html">Example 31</a>: Nedelec H(curl) FEM for the definite anisotropic Maxwell problem
+ * - <a class="el" href="ex31p_8cpp_source.html">Example 31p</a>: parallel Nedelec H(curl) FEM for the definite anisotropic Maxwell problem
+ * - <a class="el" href="ex32p_8cpp_source.html">Example 32p</a>: parallel anisotropic Maxwell eigensolver
 *
 * <H4>AmgX Examples</H4>
 * - Variants of Examples
@@ -190,12 +195,13 @@ namespace mfem {
 * - <a class="el" href="field-interp_8cpp_source.html">Field Interp</a>: transfer a grid functions between meshes
 * - <a class="el" href="distance_8cpp_source.html">Distance</a>: finite element distance function solver
 * - <a class="el" href="diffusion_8cpp_source.html">Shifted Diffusion</a>: shifted boundary diffusion solver
+ * - <a class="el" href="extrapolate_8cpp_source.html">Extrapolation</a>: PDE-based extrapolation of finite element functions
 * - <a class="el" href="distance_8cpp_source.html">Block Solvers</a>: comparison of saddle point system solvers
 * - <a class="el" href="parheat_8cpp_source.html">Optimization gradients</a>: Gradients of PDE-constrained function
+ * - <a class="el" href="par__example_8cpp_source.html">Parallel AD</a>: Parallel p-Laplacian example
+ * - <a class="el" href="seq__example_8cpp_source.html">Serial AD</a>: Serial p-Laplacian example
 * - <a class="el" href="miniapps_2performance_2ex1_8cpp_source.html">HPC Example 1</a>: high-performance nodal H1 FEM for the Laplace problem
 * - <a class="el" href="miniapps_2performance_2ex1p_8cpp_source.html">HPC Example 1p</a>: high-performance parallel nodal H1 FEM for the Laplace problem
- * - <a class="el" href="par__example_8cpp_source.html">Parallel pLaplacian example using AD</a>: Parallel pLaplacian example
- * - <a class="el" href="seq__example_8cpp_source.html">Serial pLaplacian example using AD</a>: Serial pLaplacian example
 *
 * See also the <a class="el" href="https://mfem.org/examples/">examples documentation</a> online.
 */
@@ -21,8 +21,7 @@ create the file CodeDocumentation.html, which can be viewed in any web browser.

 Alternatively, the (pre-build) documentation can be browsed directly online at

-   https://mfem.github.io/doxygen/html/index.html
-
+   https://docs.mfem.org

 Some building considerations:

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -16,7 +16,7 @@ DOXYGEN_CONF = CodeDocumentation.conf
 # doxygen uses: graphviz, latex
 html: $(DOXYGEN_CONF)
 	@# Generate the html documentation
-	@doxygen $(DOXYGEN_CONF)
+	@( cat $(DOXYGEN_CONF) ; echo "$(MFEM_DOXYGEN_FLAGS)" ) | doxygen -
 	@echo "<meta http-equiv=\"REFRESH\" content=\"0;URL=CodeDocumentation/html/index.html\">" > CodeDocumentation.html
 	@cat warnings.log 1>&2
 	@# Generate the log of undocumented methods
@@ -29,4 +29,3 @@ clean:
 $(DOXYGEN_CONF): $(MFEM_DIR)/doc/$(DOXYGEN_CONF).in
 	@sed -e 's%@MFEM_SOURCE_DIR@%$(MFEM_DIR)%g' $(<) \
 	  > $(DOXYGEN_CONF)
-
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -38,6 +38,7 @@ list(APPEND ALL_EXE_SRCS
  ex28.cpp
  ex29.cpp
  ex30.cpp
+  ex31.cpp
  )

 if (MFEM_USE_MPI)
@@ -72,6 +73,8 @@ if (MFEM_USE_MPI)
    ex28p.cpp
    ex29p.cpp
    ex30p.cpp
+    ex31p.cpp
+    ex32p.cpp
    )
 endif()

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -37,11 +37,11 @@ using namespace mfem;

 int main(int argc, char *argv[])
 {
-   // 1. 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. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();

   // 2. Parse command-line options.
   const char *mesh_file = "../../data/star.mesh";
@@ -87,7 +87,6 @@ int main(int argc, char *argv[])
      {
         args.PrintUsage(cout);
      }
-      MPI_Finalize();
      return 1;
   }
   if (myid == 0)
@@ -320,7 +319,6 @@ int main(int argc, char *argv[])
   {
      delete fec;
   }
-   MPI_Finalize();

   return 0;
 }
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -58,11 +58,11 @@ using namespace mfem;

 int main(int argc, char *argv[])
 {
-   // 1. 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. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
   // Define Caliper ConfigManager
   cali::ConfigManager mgr;
   // Caliper instrumentation
@@ -101,7 +101,6 @@ int main(int argc, char *argv[])
      {
         args.PrintUsage(cout);
      }
-      MPI_Finalize();
      return 1;
   }
   if (myid == 0)
@@ -293,7 +292,6 @@ int main(int argc, char *argv[])
   }
   // Flush output before MPI_finalize
   mgr.flush();
-   MPI_Finalize();

   return 0;
 }
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -0,0 +1,294 @@
+//                 MFEM FOSLS acoustics Example
+//
+// Compile with: make fosls
+//
+// Definite/Indefinite Helmholtz
+
+//     - Δ p ± ω^2 p = f̃ ,   in Ω
+//                 p = p_0, on ∂Ω
+
+// First Order System
+
+//   ∇ p - ω u = 0, in Ω
+// - ∇⋅u ± ω p = f, in Ω
+//           p = p_0, in ∂Ω
+// where f:=f̃/ω 
+
+// FOSLS:
+//       minimize  1/2(||∇p - ω u||^2 + ||-∇⋅u ± ω p - f||^2)
+
+// (p,u) ∈ H^1(Ω) × H(div,Ω)
+// -------------------------------------------------------------------
+// |   |            p           |            u             |   RHS    | 
+// -------------------------------------------------------------------
+// | q |  (∇ p,∇ q) + ω^2(p,q)  | ∓ ω (∇⋅u,q) - ω (u, ∇ q) | ± ω(f,q) |
+// |   |                        |                          |          |
+// | v | ∓ ω (p,∇⋅v) - ω (∇ p,v)|  (∇⋅u,∇⋅v) + ω^2 (u,v)   | -(f,∇⋅v) |
+
+// where (q,v) ∈ H^1(Ω) × H(div,Ω)
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// #define DEFINITE
+
+double p_exact(const Vector &x);
+void u_exact(const Vector &x, Vector & u);
+double rhs_func(const Vector &x);
+void gradp_exact(const Vector &x, Vector &gradu);
+double divu_exact(const Vector &x);
+double d2_exact(const Vector &x);
+
+int dim;
+double omega;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../../data/inline-quad.mesh";
+   int order = 1;
+   bool visualization = true;
+   double rnum=1.0;
+   int sr = 1;
+
+   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(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&rnum, "-rnum", "--number_of_wavelenths",
+                  "Number of wavelengths");      
+   args.AddOption(&sr, "-sr", "--serial_ref",
+                  "Number of serial refinements.");                               
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   omega = 2.0 * M_PI * rnum;
+
+
+   Mesh mesh(mesh_file, 1, 1);
+   dim = mesh.Dimension();
+
+   for (int i = 0; i < sr; i++ )
+   {
+      mesh.UniformRefinement();
+   }
+
+   FiniteElementCollection *H1fec = new H1_FECollection(order, dim);
+   FiniteElementCollection *RTfec = new RT_FECollection(order-1, dim);
+   FiniteElementSpace * H1fes = new FiniteElementSpace(&mesh, H1fec);
+   FiniteElementSpace * RTfes = new FiniteElementSpace(&mesh, RTfec);
+
+   Array<FiniteElementSpace *> fespaces(2);
+   fespaces[0] = H1fes;
+   fespaces[1] = RTfes;
+
+   Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
+   if (mesh.bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespaces[0]->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   BlockBilinearForm a(fespaces);
+   a.SetDiagonalPolicy(mfem::Operator::DIAG_KEEP);
+   cout << "H1 fespace = " << H1fes->GetTrueVSize() << endl;
+   cout << "RT fespace = " << RTfes->GetTrueVSize() << endl;
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient negone(-1.0);
+   ConstantCoefficient omeg(omega);
+   ConstantCoefficient negomeg(-omega);
+   ConstantCoefficient omeg2(omega*omega);
+
+
+   Array2D<BilinearFormIntegrator * > blfi(2,2);
+
+   // blfi(0,0) = (∇ p,∇ q) + ω^2(p,q)
+   SumIntegrator * integ00 = new SumIntegrator();
+   integ00->AddIntegrator(new DiffusionIntegrator(one));
+   integ00->AddIntegrator(new MassIntegrator(omeg2));
+   blfi(0,0) = integ00;
+
+   // blfi(0,1) = ∓ ω (∇⋅u,q) - ω (u, ∇ q)
+   SumIntegrator * integ01 = new SumIntegrator();
+#ifdef DEFINITE
+   // -ω (∇⋅u,q)
+   integ01->AddIntegrator(new MixedScalarDivergenceIntegrator(negomeg));
+#else   
+   // ω (∇⋅u,q)
+   integ01->AddIntegrator(new MixedScalarDivergenceIntegrator(omeg));
+#endif
+   // - ω (u, ∇ q)
+   integ01->AddIntegrator(new MixedVectorWeakDivergenceIntegrator(omeg));
+   blfi(0,1) = integ01;
+
+   // blfi(1,0) = ∓ ω (p,∇⋅v) - ω (∇ p,v)
+   SumIntegrator * integ10 = new SumIntegrator();
+#ifdef DEFINITE
+   // - ω (p,∇⋅v)
+   integ10->AddIntegrator(new MixedScalarWeakGradientIntegrator(omeg));
+#else
+   // ω (p,∇⋅v)
+   integ10->AddIntegrator(new MixedScalarWeakGradientIntegrator(negomeg));
+#endif
+   // - ω (∇ p,v)
+   integ10->AddIntegrator(new MixedVectorGradientIntegrator(negomeg));
+   blfi(1,0) = integ10;
+
+   // blfi(1,1) = (∇⋅u,∇⋅v) + ω^2 (u,v)
+   SumIntegrator * integ11 = new SumIntegrator();
+   integ11->AddIntegrator(new DivDivIntegrator(one));
+   integ11->AddIntegrator(new VectorFEMassIntegrator(omeg2));
+   blfi(1,1) = integ11;
+
+
+   BlockLinearForm b(fespaces);
+   Array<LinearFormIntegrator * > lfi(2);
+   // ± ω (f,q)
+   FunctionCoefficient f_rhs(rhs_func);
+#ifdef DEFINITE  
+   ProductCoefficient w_f(omeg,f_rhs);
+#else   
+   ProductCoefficient w_f(negomeg,f_rhs);
+#endif   
+   // lfi[0] = new DomainLFIntegrator(w_f);
+   lfi[0] = new DomainLFIntegrator(w_f);
+
+
+   //  -(f,∇⋅v)
+   ProductCoefficient neg_f(negone,f_rhs);
+   // lfi[1] = new VectorFEDomainLFDivIntegrator(f_rhs);
+   lfi[1] = new VectorFEDomainLFDivIntegrator(neg_f);
+
+   TestBlockBilinearFormIntegrator * integ = new TestBlockBilinearFormIntegrator();
+   integ->SetIntegrators(blfi);
+   a.AddDomainIntegrator(integ);
+   a.Assemble();
+
+   TestBlockLinearFormIntegrator * lininteg = new TestBlockLinearFormIntegrator();
+   lininteg->SetIntegrators(lfi);
+   b.AddDomainIntegrator(lininteg);
+   b.Assemble();
+
+   int size = 0;
+   for (int i = 0; i<fespaces.Size(); i++)
+   {
+      size += fespaces[i]->GetVSize();
+   }
+
+   Vector x(size);
+   x = 0.0;
+   FunctionCoefficient p_ex(p_exact);
+   VectorFunctionCoefficient gradp_ex(dim,gradp_exact);
+   VectorFunctionCoefficient u_ex(dim,u_exact);
+   FunctionCoefficient divu_ex(divu_exact);
+   GridFunction p_gf, u_gf;   
+   GridFunction pex_gf(H1fes);   
+
+   p_gf.MakeRef(H1fes,x,0);
+   // p_gf.ProjectBdrCoefficient(p_ex,ess_bdr);
+   p_gf.ProjectCoefficient(p_ex);
+   pex_gf.ProjectCoefficient(p_ex);
+
+   u_gf.MakeRef(RTfes,x,H1fes->GetVSize());
+   u_gf = 0.;
+
+
+   OperatorPtr A;
+   Vector X,B;
+   a.FormLinearSystem(ess_tdof_list,x,b,A,X,B);
+
+   GSSmoother M((SparseMatrix&)(*A));
+   CGSolver cg;
+   cg.SetRelTol(1e-10);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   cg.SetPreconditioner(M);
+   cg.SetOperator(*A);
+   cg.Mult(B, X);
+
+   a.RecoverFEMSolution(X,b,x);
+
+   p_gf.MakeRef(H1fes,x,0);
+   u_gf.MakeRef(RTfes,x,H1fes->GetVSize());
+
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream solu_sock(vishost, visport);
+      solu_sock.precision(8);
+      solu_sock << "solution\n" << mesh << p_gf <<
+                "window_title 'Numerical p' "
+                << flush;
+      // socketstream sols_sock(vishost, visport);
+      // sols_sock.precision(8);
+      // sols_sock << "solution\n" << mesh << u_gf <<
+               //  "window_title 'Numerical sigma' "
+               //  << flush;
+
+      socketstream solex_sock(vishost, visport);
+      solex_sock.precision(8);
+      solex_sock << "solution\n" << mesh << pex_gf <<
+                    "window_title 'Exact p' "
+                    << flush;          
+   }
+
+   return 0;
+}
+
+
+double rhs_func(const Vector &x)
+{
+   double p = p_exact(x);
+   double divu = divu_exact(x);
+   // f = - ∇⋅u ± ω p, 
+#ifdef DEFINITE   
+   return -divu + omega * p;
+#else
+   return -divu - omega * p;
+#endif   
+}
+
+double p_exact(const Vector &x)
+{
+   return sin(omega*x.Sum());
+}
+
+void gradp_exact(const Vector &x, Vector &grad)
+{
+   grad.SetSize(x.Size());
+   grad = omega * cos(omega * x.Sum());
+}
+
+void u_exact(const Vector &x, Vector & u)
+{
+   gradp_exact(x,u);
+   u *= 1./omega;
+}
+
+double divu_exact(const Vector &x)
+{
+   return d2_exact(x)/omega;
+}
+
+double d2_exact(const Vector &x)
+{
+   return -dim * omega * omega * sin(omega*x.Sum());
+}
@@ -0,0 +1,59 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../../..
+MFEM_BUILD_DIR ?= ../../..
+SRC = $(if $(MFEM_DIR:../../..=),$(MFEM_DIR)/examples/dpg_tests/acoustics,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = fosls uw_dpg  strong_dpg
+PAR_EXAMPLES = uw_dpgp
+
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES) 
+
+MFEM_TESTS = EXAMPLES
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel example)
+%-test-seq: %
+	@$(call mfem-test,$<,, Serial example)
+
+clean: clean-build clean-exec 
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
@@ -0,0 +1,271 @@
+//                 MFEM DPG_strong acoustics Example
+//
+// Compile with: make strong_dpg
+//
+// Definite/Indefinite Helmholtz
+
+//     - Δ p ± ω^2 p = f̃ ,   in Ω
+//                 p = p_0, on ∂Ω
+
+// First Order System
+
+//   ∇ p - ω u = 0, in Ω
+// - ∇⋅u ± ω p = f, in Ω
+//           p = p_0, in ∂Ω
+// where f:=f̃/ω 
+
+// Strong DPG formulation
+// (p,u) ∈ H^1(Ω) × H(div,Ω)
+// 
+//   (∇ p, v) - ω (u,v) = 0,     in Ω, ∀ v ∈ (L^2)^dim
+//  -(∇⋅u, q) ± ω (p,q) = (f,q), in Ω, ∀ q ∈  L^2 
+//                    p = p_0, in ∂Ω
+// 
+// ------------------------------------
+// |   |      p    |     u    |  RHS  | 
+// ------------------------------------
+// | q | ± ω (p,q) | -(∇⋅u,q) | (f,q) |
+// |   |           |          |       |
+// | v | (∇ p, v)  | -ω (u,v) |       |
+
+// where (q,v) ∈ L^2  × (L^2)^dim
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// #define DEFINITE
+
+double p_exact(const Vector &x);
+void u_exact(const Vector &x, Vector & u);
+double rhs_func(const Vector &x);
+void gradp_exact(const Vector &x, Vector &gradu);
+double divu_exact(const Vector &x);
+double d2_exact(const Vector &x);
+
+int dim;
+double omega;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../../data/inline-quad.mesh";
+   int order = 1;
+   int delta_order = 1;
+   bool visualization = true;
+   double rnum=1.0;
+   int ref = 1;
+
+   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(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&rnum, "-rnum", "--number_of_wavelenths",
+                  "Number of wavelengths");      
+   args.AddOption(&delta_order, "-do", "--delta_order",
+                  "Order enrichment for DPG test space.");                  
+   args.AddOption(&ref, "-ref", "--serial_ref",
+                  "Number of serial refinements.");                               
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   omega = 2.0 * M_PI * rnum;
+
+
+   Mesh mesh(mesh_file, 1, 1);
+   dim = mesh.Dimension();
+
+
+   for (int i = 0; i < ref; i++ )
+   {
+      mesh.UniformRefinement();
+   }
+
+   // Define spaces
+   // H1 space for p
+   FiniteElementCollection *p_fec = new H1_FECollection(order, dim);
+   FiniteElementSpace * p_fes = new FiniteElementSpace(&mesh, p_fec);
+
+   // H(div) for u
+   FiniteElementCollection *u_fec = new RT_FECollection(order-1, dim);
+   FiniteElementSpace * u_fes = new FiniteElementSpace(&mesh, u_fec);
+
+   // testspace fe collections
+   int test_order = order+delta_order;
+   FiniteElementCollection * q_fec = new L2_FECollection(test_order-1, dim);
+   FiniteElementCollection * v_fec = new L2_FECollection(test_order-1, dim);
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient negone(-1.0);
+   ConstantCoefficient omeg(omega);
+   ConstantCoefficient negomeg(-omega);
+
+   // Normal equation weak formulation
+   Array<FiniteElementSpace * > trial_fes; 
+   Array<FiniteElementCollection * > test_fec; 
+
+   trial_fes.Append(p_fes);
+   trial_fes.Append(u_fes);
+   test_fec.Append(q_fec);
+   test_fec.Append(v_fec);
+
+   NormalEquations * a = new NormalEquations(trial_fes,test_fec);
+   a->SetTestFECollVdim(1,dim);
+   
+   a->StoreMatrices(true);
+
+//  ± ω (p, q)
+#ifdef DEFINITE
+//  ω (p, q)
+   a->AddTrialIntegrator(new MassIntegrator(omeg),0,0);
+#else
+// -ω (p, q)
+   a->AddTrialIntegrator(new MassIntegrator(negomeg),0,0);
+#endif
+
+// -(∇⋅u, q) 
+   a->AddTrialIntegrator(new MixedScalarDivergenceIntegrator(negone),1,0);
+
+// -ω (u,v)
+   a->AddTrialIntegrator(new VectorFEMassIntegrator(negomeg),1,1);
+
+// (∇ p, v)
+   a->AddTrialIntegrator(new GradientIntegrator(one),0,1);
+
+// (v,δv)
+   a->AddTestIntegrator(new VectorMassIntegrator(one),1,1);
+    
+// (q,δq)
+   a->AddTestIntegrator(new MassIntegrator(one),0,0);
+
+   FunctionCoefficient f_rhs(rhs_func);
+   a->AddDomainLFIntegrator(new DomainLFIntegrator(f_rhs),0);
+
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   if (mesh.bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      p_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   FunctionCoefficient p_ex(p_exact);
+   VectorFunctionCoefficient gradp_ex(dim,gradp_exact);
+   VectorFunctionCoefficient u_ex(dim,u_exact);
+   FunctionCoefficient divu_ex(divu_exact);
+   GridFunction p_gf, u_gf;   
+   GridFunction pex_gf(p_fes);   
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = p_fes->GetVSize();
+   offsets[2] = u_fes->GetVSize();
+   offsets.PartialSum();
+   BlockVector x(offsets);
+   x = 0.0;
+
+   p_gf.MakeRef(p_fes,x.GetBlock(0));
+   p_gf.ProjectBdrCoefficient(p_ex,ess_bdr);
+
+   u_gf.MakeRef(u_fes,x.GetBlock(1));
+
+   a->Assemble();
+
+   OperatorPtr Ah;
+   Vector X,B;
+   a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
+
+   BlockMatrix * A = Ah.As<BlockMatrix>();
+
+   BlockDiagonalPreconditioner * M = new BlockDiagonalPreconditioner(A->RowOffsets());
+   M->owns_blocks = 1;
+   for (int i=0; i<A->NumRowBlocks(); i++)
+   {
+      M->SetDiagonalBlock(i,new UMFPackSolver(A->GetBlock(i,i)));
+   }
+
+   CGSolver cg;
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(3);
+   cg.SetPreconditioner(*M);
+   cg.SetOperator(*A);
+   cg.Mult(B, X);
+   delete M;
+
+   a->RecoverFEMSolution(X,x);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream p_out;
+      socketstream u_out;
+      p_out.open(vishost, visport);
+      u_out.open(vishost, visport);
+      p_out.precision(8);
+      p_out << "solution\n" << mesh << p_gf <<
+               "window_title 'Numerical p' "
+               << flush;
+
+      u_out.precision(8);
+      u_out << "solution\n" << mesh << u_gf <<
+            "window_title 'Numerical flux' "
+            << flush;
+   }
+
+   return 0;
+}
+
+
+double rhs_func(const Vector &x)
+{
+   double p = p_exact(x);
+   double divu = divu_exact(x);
+   // f = - ∇⋅u ± ω p, 
+#ifdef DEFINITE   
+   return -divu + omega * p;
+#else
+   return -divu - omega * p;
+#endif   
+}
+
+double p_exact(const Vector &x)
+{
+   return sin(omega*x.Sum());
+}
+
+void gradp_exact(const Vector &x, Vector &grad)
+{
+   grad.SetSize(x.Size());
+   grad = omega * cos(omega * x.Sum());
+}
+
+void u_exact(const Vector &x, Vector & u)
+{
+   gradp_exact(x,u);
+   u *= 1./omega;
+}
+
+double divu_exact(const Vector &x)
+{
+   return d2_exact(x)/omega;
+}
+
+double d2_exact(const Vector &x)
+{
+   return -dim * omega * omega * sin(omega*x.Sum());
+}
@@ -0,0 +1,546 @@
+//                   MFEM Ultraweak DPG acoustics example
+//
+// Compile with: make uw_dpg
+//
+// ./uw_dpg -m ../../../data/inline-quad.mesh -rnum 40 -theta 0.7 -prob 1 -graph-norm -ref 40 -o 3
+
+//     - Δ p ± ω^2 p = f̃ ,   in Ω
+//                 p = p_0, on ∂Ω
+
+// First Order System
+
+//   ∇ p - ω u = 0, in Ω
+// - ∇⋅u ± ω p = f, in Ω
+//           p = p_0, in ∂Ω
+// where f:=f̃/ω 
+
+// UW-DPG:
+// 
+// p ∈ L^2(Ω), u ∈ (L^2(Ω))^dim 
+// p̂ ∈ H^1/2(Ω), û ∈ H^-1/2(Ω)  
+// -(p,  ∇⋅v) - ω (u , v) + < p̂, v⋅n> = 0,      ∀ v ∈ H(div,Ω)      
+//  (u , ∇ q) ± ω (p , q) + < û, q >  = (f,q)   ∀ q ∈ H^1(Ω)
+//                                  p̂ = p_0     on ∂Ω 
+
+// Note: 
+// p̂ := p on Γ_h (skeleton)
+// û := -u on Γ_h  
+
+// -------------------------------------------------------------
+// |   |     p     |     u     |    p̂      |    û    |  RHS    |
+// -------------------------------------------------------------
+// | v | -(p, ∇⋅v) | - ω (u,v) | < p̂, v⋅n> |         |         |
+// |   |           |           |           |         |         |
+// | q | ± ω (p,q) | (u , ∇ q) |           | < û,q > |  (f,q)  |  
+
+// where (q,v) ∈  H^1(Ω) × H(div,Ω) 
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// #define DEFINITE
+void acoustics_solution(const Vector & X, double & p, Vector & dp, double & d2p);
+double p_exact(const Vector &x);
+void u_exact(const Vector &x, Vector & u);
+double rhs_func(const Vector &x);
+double divu_exact(const Vector &x);
+double hatp_exact(const Vector & X);
+void hatu_exact(const Vector & X, Vector & hatu);
+
+int dim;
+double omega;
+
+enum prob_type
+{
+   plane_wave,
+   gaussian_beam  
+};
+
+prob_type prob;
+
+int main(int argc, char *argv[])
+{
+   const char *mesh_file = "../../../data/inline-quad.mesh";
+   int order = 1;
+   int delta_order = 1;
+   bool visualization = true;
+   double rnum=1.0;
+   int ref = 1;
+   double theta = 0.0;
+   bool adjoint_graph_norm = false;
+   int iprob = 0;
+
+   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(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&rnum, "-rnum", "--number_of_wavelenths",
+                  "Number of wavelengths");      
+   args.AddOption(&delta_order, "-do", "--delta_order",
+                  "Order enrichment for DPG test space.");     
+   args.AddOption(&theta, "-theta", "--theta",
+                  "Theta parameter for AMR");    
+   args.AddOption(&iprob, "-prob", "--problem", "Problem case"
+                  " 0: plane wave, 1: Gaussian beam");                                    
+   args.AddOption(&adjoint_graph_norm, "-graph-norm", "--adjoint-graph-norm",
+                  "-no-graph-norm", "--no-adjoint-graph-norm",
+                  "Enable or disable Adjoint Graph Norm on the test space");                                
+   args.AddOption(&ref, "-ref", "--serial_ref",
+                  "Number of serial refinements.");                               
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   if (iprob > 1) { iprob = 0; }
+   prob = (prob_type)iprob;
+
+
+   omega = 2.0 * M_PI * rnum;
+
+
+   Mesh mesh(mesh_file, 1, 1);
+   dim = mesh.Dimension();
+
+
+   // Define spaces
+   // L2 space for p
+   FiniteElementCollection *p_fec = new L2_FECollection(order-1,dim);
+   FiniteElementSpace *p_fes = new FiniteElementSpace(&mesh,p_fec);
+
+   // Vector L2 space for u 
+   FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
+   FiniteElementSpace *u_fes = new FiniteElementSpace(&mesh,u_fec, dim); 
+
+   // H^1/2 space for p̂  
+   FiniteElementCollection * hatp_fec = new H1_Trace_FECollection(order,dim);
+   FiniteElementSpace *hatp_fes = new FiniteElementSpace(&mesh,hatp_fec);
+
+   // H^-1/2 space for û  
+   FiniteElementCollection * hatu_fec = new RT_Trace_FECollection(order-1,dim);   
+   FiniteElementSpace *hatu_fes = new FiniteElementSpace(&mesh,hatu_fec);
+
+   // testspace fe collections
+   int test_order = order+delta_order;
+   FiniteElementCollection * q_fec = new H1_FECollection(test_order, dim);
+   FiniteElementCollection * v_fec = new RT_FECollection(test_order-1, dim);
+
+
+   // Coefficients
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(0.0);
+   Vector vec0(dim); vec0 = 0.;
+   VectorConstantCoefficient vzero(vec0);
+   ConstantCoefficient negone(-1.0);
+   ConstantCoefficient omeg(omega);
+   ConstantCoefficient omeg2(omega*omega);
+   ConstantCoefficient negomeg(-omega);
+
+   // Normal equation weak formulation
+   Array<FiniteElementSpace * > trial_fes; 
+   Array<FiniteElementCollection * > test_fec; 
+
+   trial_fes.Append(p_fes);
+   trial_fes.Append(u_fes);
+   trial_fes.Append(hatp_fes);
+   trial_fes.Append(hatu_fes);
+
+   test_fec.Append(q_fec);
+   test_fec.Append(v_fec);
+
+   NormalEquations * a = new NormalEquations(trial_fes,test_fec);
+   a->StoreMatrices(true);
+
+
+   // ± ω (p,q)
+#ifdef DEFINITE
+   a->AddTrialIntegrator(new MixedScalarMassIntegrator(omeg),0,0);
+#else
+   a->AddTrialIntegrator(new MixedScalarMassIntegrator(negomeg),0,0);
+#endif   
+
+// (u , ∇ q)
+   a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),1,0);
+
+// -(p, ∇⋅v)
+   a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(one),0,1);
+
+// - ω (u,v)
+   a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(negomeg)),1,1);
+
+// < p̂, v⋅n>
+   a->AddTrialIntegrator(new NormalTraceIntegrator,2,1);
+
+// < û,q >
+   a->AddTrialIntegrator(new TraceIntegrator,3,0);
+
+
+// test integrators 
+
+   //space-induced norm for H(div) × H1
+   // (∇q,∇δq)
+   a->AddTestIntegrator(new DiffusionIntegrator(one),0,0);
+   // (q,δq)
+   a->AddTestIntegrator(new MassIntegrator(one),0,0);
+   // (∇⋅v,∇⋅δv)
+   a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+   // (v,δv)
+   a->AddTestIntegrator(new VectorFEMassIntegrator(one),1,1);
+
+   // additional integrators for the adjoint graph norm
+   if (adjoint_graph_norm)
+   {   
+      // -ω (∇q,δv)
+      a->AddTestIntegrator(new MixedVectorGradientIntegrator(negomeg),0,1);
+      // -ω (v,δq)
+      a->AddTestIntegrator(new MixedVectorWeakDivergenceIntegrator(omeg),1,0);
+      // ω^2 (v,δv)
+      a->AddTestIntegrator(new VectorFEMassIntegrator(omeg2),1,1);
+
+#ifdef DEFINITE
+      // - ω (∇⋅v,δq)   
+      a->AddTestIntegrator(new VectorFEDivergenceIntegrator(negomeg),1,0);
+      // - ω (q,∇⋅v)   
+      a->AddTestIntegrator(new MixedScalarWeakGradientIntegrator(omeg),0,1);
+#else
+      // ω (∇⋅v,δq)   
+      a->AddTestIntegrator(new VectorFEDivergenceIntegrator(omeg),1,0);
+      // ω (q,∇⋅v)   
+      a->AddTestIntegrator(new MixedScalarWeakGradientIntegrator(negomeg),0,1);
+#endif
+      // ω^2 (q,δq)
+      a->AddTestIntegrator(new MassIntegrator(omeg2),0,0);
+   }
+
+   // RHS
+   FunctionCoefficient f_rhs(rhs_func);
+   a->AddDomainLFIntegrator(new DomainLFIntegrator(f_rhs),0);
+
+
+   FunctionCoefficient hatpex(hatp_exact);
+   FunctionCoefficient pex(p_exact);
+   VectorFunctionCoefficient uex(dim,u_exact);
+   Array<int> elements_to_refine;
+   GridFunction hatp_gf;
+
+
+   socketstream p_out;
+   // socketstream u_out;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      p_out.open(vishost, visport);
+      // u_out.open(vishost, visport);
+   }
+
+   double res0 = 0.;
+   double err0 = 0.;
+   int dof0;
+   mfem::out << " Refinement |" 
+             << "    Dofs    |" 
+             << "  L2 Error  |" 
+             << " Relative % |" 
+             << "  Rate  |" 
+             << "  Residual  |" 
+             << "  Rate  |" << endl;
+   mfem::out << " --------------------"      
+             <<  "-------------------"    
+             <<  "-------------------"    
+             <<  "-------------------" << endl;   
+
+
+   for (int i = 0; i<ref; i++)
+   {
+      a->Assemble();
+
+      Array<int> ess_tdof_list;
+      Array<int> ess_bdr;
+      if (mesh.bdr_attributes.Size())
+      {
+         ess_bdr.SetSize(mesh.bdr_attributes.Max());
+         ess_bdr = 1;
+         hatp_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+
+   // shift the ess_tdofs
+      for (int i = 0; i < ess_tdof_list.Size(); i++)
+      {
+         ess_tdof_list[i] += p_fes->GetTrueVSize() + u_fes->GetTrueVSize();
+      }
+
+      Array<int> offsets(5);
+      offsets[0] = 0;
+      offsets[1] = p_fes->GetVSize();
+      offsets[2] = u_fes->GetVSize();
+      offsets[3] = hatp_fes->GetVSize();
+      offsets[4] = hatu_fes->GetVSize();
+      offsets.PartialSum();
+      BlockVector x(offsets);
+      x = 0.0;
+      hatp_gf.MakeRef(hatp_fes,x.GetBlock(2));
+      hatp_gf.ProjectBdrCoefficient(hatpex,ess_bdr);
+
+      OperatorPtr Ah;
+      Vector X,B;
+      a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
+
+      BlockMatrix * A = Ah.As<BlockMatrix>();
+
+      BlockDiagonalPreconditioner * M = new BlockDiagonalPreconditioner(A->RowOffsets());
+      M->owns_blocks = 1;
+      for (int i=0; i<A->NumRowBlocks(); i++)
+      {
+         M->SetDiagonalBlock(i,new UMFPackSolver(A->GetBlock(i,i)));
+      }
+
+      CGSolver cg;
+      cg.SetRelTol(1e-8);
+      cg.SetMaxIter(20000);
+      cg.SetPrintLevel(3);
+      cg.SetPreconditioner(*M);
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+      delete M;
+
+      a->RecoverFEMSolution(X,x);
+      Vector & residuals = a->ComputeResidual(x);
+
+      double residual = residuals.Norml2();
+
+      elements_to_refine.SetSize(0);
+      double max_resid = residuals.Max();
+      for (int iel = 0; iel<mesh.GetNE(); iel++)
+      {
+         if (residuals[iel] > theta * max_resid)
+         {
+            elements_to_refine.Append(iel);
+         }
+      }
+      
+      GridFunction p_gf;
+      p_gf.MakeRef(p_fes,x.GetBlock(0));
+
+      GridFunction u_gf;
+      u_gf.MakeRef(u_fes,x.GetBlock(1));
+
+      GridFunction pex_gf(p_fes);
+      GridFunction uex_gf(u_fes);
+      pex_gf.ProjectCoefficient(pex);
+      uex_gf.ProjectCoefficient(uex);
+
+
+      // Error
+      int dofs = X.Size();
+      double p_err = p_gf.ComputeL2Error(pex);
+      double p_norm = uex_gf.ComputeL2Error(zero);
+      double u_err = u_gf.ComputeL2Error(uex);
+      double u_norm = u_gf.ComputeL2Error(vzero);
+
+      double L2Error = sqrt(p_err*p_err + u_err*u_err);
+      double L2norm = sqrt(p_norm * p_norm + u_norm * u_norm);
+
+      double rel_error = L2Error/L2norm;
+
+      double rate_err = (i) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0;
+      double rate_res = (i) ? dim*log(res0/residual)/log((double)dof0/dofs) : 0.0;
+
+      err0 = L2Error;
+      res0 = residual;
+      dof0 = dofs;
+      mfem::out << std::right << std::setw(11) << i << " | " 
+                << std::setw(10) <<  dof0 << " | " 
+                << std::setprecision(3) 
+                << std::setw(10) << std::scientific <<  err0 << " | " 
+                << std::setprecision(3) 
+                << std::setw(10) << std::fixed <<  rel_error * 100. << " | " 
+                << std::setprecision(2) 
+                << std::setw(6) << std::fixed << rate_err << " | " 
+                << std::setprecision(3) 
+                << std::setw(10) << std::scientific <<  res0 << " | " 
+                << std::setprecision(2) 
+                << std::setw(6) << std::fixed << rate_res << " | " 
+                << std::resetiosflags(std::ios::showbase)
+                << std::endl;
+
+      if (visualization)
+      {
+         p_out.precision(8);
+         p_out << "solution\n" << mesh << p_gf <<
+                  "window_title 'Numerical presure' "
+                  << flush;
+
+         // u_out.precision(8);
+         // u_out << "solution\n" << mesh << u_gf <<
+         //       "window_title 'Numerical velocity' "
+         //       << flush;
+      }
+
+      if (i == ref)
+         break;
+
+      mesh.GeneralRefinement(elements_to_refine,1,1);
+      for (int i =0; i<trial_fes.Size(); i++)
+      {
+         trial_fes[i]->Update(false);
+      }
+      a->Update();
+   }
+
+   delete a;
+   delete q_fec;
+   delete v_fec;
+   delete hatp_fes;
+   delete hatp_fec;
+   delete hatu_fes;
+   delete hatu_fec;
+   delete u_fec;
+   delete p_fec;
+   delete u_fes;
+   delete p_fes;
+
+   return 0;
+}
+
+double rhs_func(const Vector &x)
+{
+   double p = p_exact(x);
+   double divu = divu_exact(x);
+   // f = - ∇⋅u ± ω p, 
+#ifdef DEFINITE   
+   return -divu + omega * p;
+#else
+   return -divu - omega * p;
+#endif   
+}
+
+double p_exact(const Vector &x)
+{
+   double p, d2p;
+   Vector dp;
+   acoustics_solution(x,p,dp,d2p);
+   return p;
+}
+
+void u_exact(const Vector &x, Vector & u)
+{
+   double p, d2p;
+   acoustics_solution(x,p,u,d2p);
+   u *= 1./omega;
+}
+
+double divu_exact(const Vector &x)
+{
+   double p, d2p;
+   Vector dp;
+   acoustics_solution(x,p,dp,d2p);
+   return d2p/omega;
+}
+
+double hatp_exact(const Vector & X)
+{
+   return p_exact(X);
+}
+
+void hatu_exact(const Vector & X, Vector & hatu)
+{
+   u_exact(X,hatu);
+   hatu *= -1.;
+}
+
+void acoustics_solution(const Vector & X, double & p, Vector & dp, double & d2p)
+{
+   dp.SetSize(X.Size());
+   switch (prob)
+   {
+   case plane_wave:
+   {
+      p = sin(omega*X.Sum());
+      dp = omega * cos(omega * X.Sum());
+      d2p = -dim * omega * omega * sin(omega*X.Sum());
+   }
+      break;
+   default:
+   {
+      double rk = omega;
+      double alpha = 45 * M_PI/180.;
+      double sina = sin(alpha); 
+      double cosa = cos(alpha);
+      // shift the origin
+      double xprim=X(0) + 0.1; 
+      double yprim=X(1) + 0.1;
+
+      double  x = xprim*sina - yprim*cosa;
+      double  y = xprim*cosa + yprim*sina;
+      double  dxdxprim = sina, dxdyprim = -cosa;
+      double  dydxprim = cosa, dydyprim =  sina;
+      //wavelength
+      double rl = 2.*M_PI/rk;
+
+      // beam waist radius
+      double w0 = 0.05;
+
+      // function w
+      double fact = rl/M_PI/(w0*w0);
+      double aux = 1. + (fact*y)*(fact*y);
+
+      double w = w0*sqrt(aux);
+      double dwdy = w0*fact*fact*y/sqrt(aux);
+      double d2wdydy = w0*fact*fact*(1. - (fact*y)*(fact*y)/aux)/sqrt(aux);
+
+      double phi0 = atan(fact*y);
+      double dphi0dy = cos(phi0)*cos(phi0)*fact;
+      double d2phi0dydy = -2.*cos(phi0)*sin(phi0)*fact*dphi0dy;
+
+      double r = y + 1./y/(fact*fact);
+      double drdy = 1. - 1./(y*y)/(fact*fact);
+      double d2rdydy = 2./(y*y*y)/(fact*fact);
+
+      // pressure
+      complex<double> zi = complex<double>(0., 1.);
+      complex<double> ze = - x*x/(w*w) - zi*rk*y - zi * M_PI * x * x/rl/r + zi*phi0/2.;
+
+      complex<double> zdedx = -2.*x/(w*w) - 2.*zi*M_PI*x/rl/r;
+      complex<double> zdedy = 2.*x*x/(w*w*w)*dwdy - zi*rk + zi*M_PI*x*x/rl/(r*r)*drdy + zi*dphi0dy/2.;
+      complex<double> zd2edxdx = -2./(w*w) - 2.*zi*M_PI/rl/r;
+      complex<double> zd2edxdy = 4.*x/(w*w*w)*dwdy + 2.*zi*M_PI*x/rl/(r*r)*drdy;
+      complex<double> zd2edydx = zd2edxdy;
+      complex<double> zd2edydy = -6.*x*x/(w*w*w*w)*dwdy*dwdy + 2.*x*x/(w*w*w)*d2wdydy - 2.*zi*M_PI*x*x/rl/(r*r*r)*drdy*drdy
+                              + zi*M_PI*x*x/rl/(r*r)*d2rdydy + zi/2.*d2phi0dydy;
+
+      double pf = pow(2.0/M_PI/(w*w),0.25);
+      double dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
+      double d2pfdydy = -1./M_PI*pow(2./M_PI,-0.75)*(-1.5*pow(w,-2.5)
+                        *dwdy*dwdy + pow(w,-1.5)*d2wdydy);
+
+
+      complex<double> zp = pf*exp(ze);
+      complex<double> zdpdx = zp*zdedx;
+      complex<double> zdpdy = dpfdy*exp(ze)+zp*zdedy;
+      complex<double> zd2pdxdx = zdpdx*zdedx + zp*zd2edxdx;
+      complex<double> zd2pdxdy = zdpdy*zdedx + zp*zd2edxdy;
+      complex<double> zd2pdydx = dpfdy*exp(ze)*zdedx + zdpdx*zdedy + zp*zd2edydx;
+      complex<double> zd2pdydy = d2pfdydy*exp(ze) + dpfdy*exp(ze)*zdedy + zdpdy*zdedy + zp*zd2edydy;
+
+      p = zp.real();
+      dp[0] = (zdpdx*dxdxprim + zdpdy*dydxprim).real();
+      dp[1] = (zdpdx*dxdyprim + zdpdy*dydyprim).real();
+
+      d2p = (  (zd2pdxdx*dxdxprim + zd2pdydx*dydxprim)*dxdxprim + (zd2pdxdy*dxdxprim + zd2pdydy*dydxprim)*dydxprim
+            + (zd2pdxdx*dxdyprim + zd2pdydx*dydyprim)*dxdyprim + (zd2pdxdy*dxdyprim + zd2pdydy*dydyprim)*dydyprim ).real();
+   }
+   break;
+   }
+
+}
@@ -0,0 +1,507 @@
+//                   MFEM Ultraweak DPG MPI acoustics (Helmholtz) example
+//
+// Compile with: make uw_dpgp
+//
+//     - Δ p ± ω^2 p = f̃ ,   in Ω
+//                 p = p_0, on ∂Ω
+//
+// First Order System
+
+//   ∇ p - ω u = 0, in Ω
+// - ∇⋅u ± ω p = f, in Ω
+//           p = p_0, in ∂Ω
+// where f:=f̃/ω 
+// 
+// UW-DPG:
+// 
+// p ∈ L^2(Ω), u ∈ (L^2(Ω))^dim 
+// p̂ ∈ H^1/2(Ω), û ∈ H^-1/2(Ω)  
+// -(p,  ∇⋅v) - ω (u , v) + < p̂, v⋅n> = 0,      ∀ v ∈ H(div,Ω)      
+//  (u , ∇ q) ± ω (p , q) + < û, q >  = (f,q)   ∀ q ∈ H^1(Ω)
+//                                  p̂ = p_0     on ∂Ω 
+
+// Note: 
+// p̂ := p on Γ_h (skeleton)
+// û := -u on Γ_h  
+
+// -------------------------------------------------------------
+// |   |     p     |     u     |    p̂      |    û    |  RHS    |
+// -------------------------------------------------------------
+// | v | -(p, ∇⋅v) | - ω (u,v) | < p̂, v⋅n> |         |         |
+// |   |           |           |           |         |         |
+// | q | ± ω (p,q) | (u , ∇ q) |           | < û,q > |  (f,q)  |  
+
+// where (q,v) ∈  H^1(Ω) × H(div,Ω) 
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// #define DEFINITE
+
+double p_exact(const Vector &x);
+void u_exact(const Vector &x, Vector & u);
+double rhs_func(const Vector &x);
+void gradp_exact(const Vector &x, Vector &gradu);
+double divu_exact(const Vector &x);
+double d2_exact(const Vector &x);
+double hatp_exact(const Vector & X);
+void hatu_exact(const Vector & X, Vector & hatu);
+
+int dim;
+double omega;
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi;
+   int num_procs = mpi.WorldSize();
+   int myid = mpi.WorldRank();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../../data/inline-quad.mesh";
+   int order = 1;
+   int delta_order = 1;
+   bool visualization = true;
+   double rnum=1.0;
+   int ref = 1;
+   double theta = 0.0;
+   bool adjoint_graph_norm = false;
+
+   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(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&rnum, "-rnum", "--number_of_wavelenths",
+                  "Number of wavelengths");      
+   args.AddOption(&delta_order, "-do", "--delta_order",
+                  "Order enrichment for DPG test space.");     
+   args.AddOption(&theta, "-theta", "--theta",
+                  "Theta parameter for AMR");                    
+   args.AddOption(&adjoint_graph_norm, "-graph-norm", "--adjoint-graph-norm",
+                  "-no-graph-norm", "--no-adjoint-graph-norm",
+                  "Enable or disable Adjoint Graph Norm on the test space");                                
+   args.AddOption(&ref, "-ref", "--serial_ref",
+                  "Number of serial refinements.");  
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   omega = 2.0 * M_PI * rnum;
+
+   Mesh mesh(mesh_file, 1, 1);
+   dim = mesh.Dimension();
+
+   mesh.EnsureNCMesh();
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // Define spaces
+   // L2 space for p
+   FiniteElementCollection *p_fec = new L2_FECollection(order-1,dim);
+   ParFiniteElementSpace *p_fes = new ParFiniteElementSpace(&pmesh,p_fec);
+
+   // Vector L2 space for u 
+   FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
+   ParFiniteElementSpace *u_fes = new ParFiniteElementSpace(&pmesh,u_fec, dim); 
+
+   // H^1/2 space for p̂  
+   FiniteElementCollection * hatp_fec = new H1_Trace_FECollection(order,dim);
+   ParFiniteElementSpace *hatp_fes = new ParFiniteElementSpace(&pmesh,hatp_fec);
+
+   // H^-1/2 space for û  
+   FiniteElementCollection * hatu_fec = new RT_Trace_FECollection(order-1,dim);   
+   ParFiniteElementSpace *hatu_fes = new ParFiniteElementSpace(&pmesh,hatu_fec);
+
+   // testspace fe collections
+   int test_order = order+delta_order;
+   FiniteElementCollection * q_fec = new H1_FECollection(test_order, dim);
+   FiniteElementCollection * v_fec = new RT_FECollection(test_order-1, dim);
+
+
+   Array<ParFiniteElementSpace * > trial_fes; 
+   trial_fes.Append(p_fes);
+   trial_fes.Append(u_fes);
+   trial_fes.Append(hatp_fes);
+   trial_fes.Append(hatu_fes);
+
+   Array<FiniteElementCollection * > test_fec; 
+   test_fec.Append(q_fec);
+   test_fec.Append(v_fec);
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(0.0);
+   Vector vec0(dim); vec0 = 0.;
+   VectorConstantCoefficient vzero(vec0);
+   ConstantCoefficient negone(-1.0);
+   ConstantCoefficient omeg(omega);
+   ConstantCoefficient omeg2(omega*omega);
+   ConstantCoefficient negomeg(-omega);
+
+   ParNormalEquations * a = new ParNormalEquations(trial_fes,test_fec);
+   a->StoreMatrices(true);
+
+
+   // Integrators
+
+   // ± ω (p,q)
+#ifdef DEFINITE
+   a->AddTrialIntegrator(new MixedScalarMassIntegrator(omeg),0,0);
+#else
+   a->AddTrialIntegrator(new MixedScalarMassIntegrator(negomeg),0,0);
+#endif   
+
+// (u , ∇ q)
+   a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),1,0);
+
+// -(p, ∇⋅v)
+   a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(one),0,1);
+
+// - ω (u,v)
+   a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(negomeg)),1,1);
+
+// < p̂, v⋅n>
+   a->AddTrialIntegrator(new NormalTraceIntegrator,2,1);
+
+// < û,q >
+   a->AddTrialIntegrator(new TraceIntegrator,3,0);
+
+
+   // test integrators 
+
+   //space-induced norm for H(div) × H1
+   // (∇q,∇δq)
+   a->AddTestIntegrator(new DiffusionIntegrator(one),0,0);
+   // (q,δq)
+   a->AddTestIntegrator(new MassIntegrator(one),0,0);
+   // (∇⋅v,∇⋅δv)
+   a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+   // (v,δv)
+   a->AddTestIntegrator(new VectorFEMassIntegrator(one),1,1);
+
+   // additional integrators for the adjoint graph norm
+   if (adjoint_graph_norm)
+   {   
+      // -ω (∇q,δv)
+      a->AddTestIntegrator(new MixedVectorGradientIntegrator(negomeg),0,1);
+      // -ω (v,δq)
+      a->AddTestIntegrator(new MixedVectorWeakDivergenceIntegrator(omeg),1,0);
+      // ω^2 (v,δv)
+      a->AddTestIntegrator(new VectorFEMassIntegrator(omeg2),1,1);
+
+#ifdef DEFINITE
+      // - ω (∇⋅v,δq)   
+      a->AddTestIntegrator(new VectorFEDivergenceIntegrator(negomeg),1,0);
+      // - ω (q,∇⋅v)   
+      a->AddTestIntegrator(new MixedScalarWeakGradientIntegrator(omeg),0,1);
+#else
+      // ω (∇⋅v,δq)   
+      a->AddTestIntegrator(new VectorFEDivergenceIntegrator(omeg),1,0);
+      // ω (q,∇⋅v)   
+      a->AddTestIntegrator(new MixedScalarWeakGradientIntegrator(negomeg),0,1);
+#endif
+      // ω^2 (q,δq)
+      a->AddTestIntegrator(new MassIntegrator(omeg2),0,0);
+   }   
+
+   // RHS
+   FunctionCoefficient f_rhs(rhs_func);
+   a->AddDomainLFIntegrator(new DomainLFIntegrator(f_rhs),0);
+
+
+   FunctionCoefficient hatpex(hatp_exact);
+   FunctionCoefficient pex(p_exact);
+   VectorFunctionCoefficient uex(dim,u_exact);
+   Array<int> elements_to_refine;
+   ParGridFunction hatp_gf;
+
+
+
+
+   socketstream p_out;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      p_out.open(vishost, visport);
+   }
+   double res0 = 0.;
+   double err0 = 0.;
+   int dof0;
+   if (myid == 0)
+   {
+      mfem::out << " Refinement |" 
+               << "    Dofs    |" 
+               << "  L2 Error  |" 
+               << " Relative % |" 
+               << "  Rate  |" 
+               << "  Residual  |" 
+               << "  Rate  |" << endl;
+      mfem::out << " --------------------"      
+               <<  "-------------------"    
+               <<  "-------------------"    
+               <<  "-------------------" << endl;   
+   }
+
+
+   for (int i = 0; i<ref; i++)
+   {
+      a->Assemble();
+
+      Array<int> ess_tdof_list;
+      Array<int> ess_bdr;
+      if (pmesh.bdr_attributes.Size())
+      {
+         ess_bdr.SetSize(pmesh.bdr_attributes.Max());
+         ess_bdr = 1;
+         hatp_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+
+      // shift the ess_tdofs
+      for (int i = 0; i < ess_tdof_list.Size(); i++)
+      {
+         ess_tdof_list[i] += p_fes->GetTrueVSize() + u_fes->GetTrueVSize();
+      }
+
+      Array<int> offsets(5);
+      offsets[0] = 0;
+      offsets[1] = p_fes->GetVSize();
+      offsets[2] = u_fes->GetVSize();
+      offsets[3] = hatp_fes->GetVSize();
+      offsets[4] = hatu_fes->GetVSize();
+      offsets.PartialSum();
+      BlockVector x(offsets);
+      x = 0.0;
+      hatp_gf.MakeRef(hatp_fes,x.GetBlock(2));
+      hatp_gf.ProjectBdrCoefficient(hatpex,ess_bdr);
+
+      Vector X,B;
+      OperatorPtr Ah;
+      a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
+
+      BlockOperator * A = Ah.As<BlockOperator>();
+
+      BlockDiagonalPreconditioner * M = new BlockDiagonalPreconditioner(A->RowOffsets());
+      M->owns_blocks = 1;
+
+      HypreBoomerAMG * amg0 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(0,0));
+      HypreBoomerAMG * amg1 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(1,1));
+      HypreBoomerAMG * amg2 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(2,2));
+      amg0->SetPrintLevel(0);
+      amg1->SetPrintLevel(0);
+      amg2->SetPrintLevel(0);
+      amg0->SetRelaxType(16);
+      amg1->SetRelaxType(16);
+      amg2->SetRelaxType(16);
+
+      M->SetDiagonalBlock(0,amg0);
+      M->SetDiagonalBlock(1,amg1);
+      M->SetDiagonalBlock(2,amg2);
+      // for (int i = 0; i < 3; i++)
+      // {
+      //    MUMPSSolver * mumps = new MUMPSSolver;
+      //    mumps->SetOperator(A->GetBlock(i,i));
+      //    M->SetDiagonalBlock(i,mumps);
+      // }
+
+      HypreSolver * prec;
+      if (dim == 2)
+      {
+         prec = new HypreAMS((HypreParMatrix &)A->GetBlock(3,3), hatu_fes);
+      }
+      else
+      {
+         prec = new HypreADS((HypreParMatrix &)A->GetBlock(3,3), hatu_fes);
+      }
+      M->SetDiagonalBlock(3,prec);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-10);
+      cg.SetMaxIter(20000);
+      cg.SetPrintLevel(-1);
+      cg.SetPreconditioner(*M);
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+      delete M;
+
+      a->RecoverFEMSolution(X,x);
+
+      Vector & residuals = a->ComputeResidual(x);
+
+      double residual = residuals.Norml2();
+      
+      double maxresidual = residuals.Max(); 
+      double globalresidual = residual * residual; 
+
+      MPI_Allreduce(MPI_IN_PLACE,&maxresidual,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
+      MPI_Allreduce(MPI_IN_PLACE,&globalresidual,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+
+      globalresidual = sqrt(globalresidual);
+
+
+      elements_to_refine.SetSize(0);
+      for (int iel = 0; iel<pmesh.GetNE(); iel++)
+      {
+         if (residuals[iel] > theta * maxresidual)
+         {
+            elements_to_refine.Append(iel);
+         }
+      }
+
+
+      ParGridFunction p_gf;
+      p_gf.MakeRef(p_fes,x.GetBlock(0));
+
+      ParGridFunction u_gf;
+      u_gf.MakeRef(u_fes,x.GetBlock(1));
+
+
+      ParGridFunction pex_gf(p_fes);
+      ParGridFunction uex_gf(u_fes);
+      pex_gf.ProjectCoefficient(pex);
+      uex_gf.ProjectCoefficient(uex);
+
+      int dofs = p_fes->GlobalTrueVSize()
+               + u_fes->GlobalTrueVSize()
+               + hatp_fes->GlobalTrueVSize()
+               + hatu_fes->GlobalTrueVSize();
+
+      double p_err = p_gf.ComputeL2Error(pex);
+      double p_norm = pex_gf.ComputeL2Error(zero);
+      double u_err = u_gf.ComputeL2Error(uex);
+      double u_norm = uex_gf.ComputeL2Error(vzero);
+
+      double L2Error = sqrt(p_err*p_err + u_err*u_err);
+      double L2norm = sqrt(p_norm * p_norm + u_norm * u_norm);
+
+      double rel_error = L2Error/L2norm;
+
+      double rate_err = (i) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0;
+      double rate_res = (i) ? dim*log(res0/globalresidual)/log((double)dof0/dofs) : 0.0;
+
+      err0 = L2Error;
+      res0 = globalresidual;
+      dof0 = dofs;
+
+         std::ios oldState(nullptr);
+      if (myid == 0)
+      {
+         mfem::out << std::right << std::setw(11) << i << " | " 
+                   << std::setw(10) <<  dof0 << " | " 
+                   << std::setprecision(3) 
+                   << std::setw(10) << std::scientific <<  err0 << " | " 
+                   << std::setprecision(3) 
+                   << std::setw(10) << std::fixed <<  rel_error * 100. << " | " 
+                   << std::setprecision(2) 
+                   << std::setw(6) << std::fixed << rate_err << " | " 
+                   << std::setprecision(3) 
+                   << std::setw(10) << std::scientific <<  res0 << " | " 
+                   << std::setprecision(2) 
+                   << std::setw(6) << std::fixed << rate_res << " | " 
+                   << std::setprecision(5) 
+                   << std::scientific 
+                   << std::endl;
+      }   
+
+      if (visualization)
+      {
+         p_out << "parallel " << num_procs << " " << myid << "\n";
+         p_out.precision(8);
+         p_out << "solution\n" << pmesh << p_gf <<
+                  "window_title 'Numerical pressure' "
+                  << flush;
+      }
+
+      if (i == ref)
+         break;
+
+      pmesh.GeneralRefinement(elements_to_refine,1,1);
+      for (int i =0; i<trial_fes.Size(); i++)
+      {
+         trial_fes[i]->Update(false);
+      }
+      a->Update();
+
+   }
+
+   delete a;
+   delete q_fec;
+   delete v_fec;
+   delete hatp_fes;
+   delete hatp_fec;
+   delete hatu_fes;
+   delete hatu_fec;
+   delete u_fec;
+   delete p_fec;
+   delete u_fes;
+   delete p_fes;
+
+   return 0;
+}
+
+double rhs_func(const Vector &x)
+{
+   double p = p_exact(x);
+   double divu = divu_exact(x);
+   // f = - ∇⋅u ± ω p, 
+#ifdef DEFINITE   
+   return -divu + omega * p;
+#else
+   return -divu - omega * p;
+#endif   
+}
+
+double p_exact(const Vector &x)
+{
+   return sin(omega*x.Sum());
+}
+
+void gradp_exact(const Vector &x, Vector &grad)
+{
+   grad.SetSize(x.Size());
+   grad = omega * cos(omega * x.Sum());
+}
+
+void u_exact(const Vector &x, Vector & u)
+{
+   gradp_exact(x,u);
+   u *= 1./omega;
+}
+
+double divu_exact(const Vector &x)
+{
+   return d2_exact(x)/omega;
+}
+
+double d2_exact(const Vector &x)
+{
+   return -dim * omega * omega * sin(omega*x.Sum());
+}
+
+double hatp_exact(const Vector & X)
+{
+   return p_exact(X);
+}
+
+void hatu_exact(const Vector & X, Vector & hatu)
+{
+   u_exact(X,hatu);
+   hatu *= -1.;
+}
@@ -0,0 +1,59 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../../..
+MFEM_BUILD_DIR ?= ../../..
+SRC = $(if $(MFEM_DIR:../../..=),$(MFEM_DIR)/examples/dpg_tests/convection-diffusion,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = uw_dpg
+PAR_EXAMPLES = uw_dpgp
+
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES) 
+
+MFEM_TESTS = EXAMPLES
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel example)
+%-test-seq: %
+	@$(call mfem-test,$<,, Serial example)
+
+clean: clean-build clean-exec 
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
@@ -0,0 +1,649 @@
+//                     MFEM Ultraweak DPG example
+//
+// Compile with: make uw_dpg
+//
+// sample runs 
+// ./uw_dpg  -m ../../../data/inline-quad.mesh -o 3 -ref 10 -test-norm 2 -do 1 -prob 1 -eps 1e-4
+//     - εΔu + ∇⋅(βu) = f,   in Ω
+//                  u = u_0, on ∂Ω
+
+// First Order System
+
+//     - ∇⋅σ + ∇⋅(βu) = f,   in Ω
+//        1/ε σ - ∇u  = 0,   in Ω
+//                  u = u_0, on ∂Ω
+
+// UW-DPG:
+// 
+// u ∈ L^2(Ω), σ ∈ (L^2(Ω))^dim 
+// û ∈ H^1/2, σ̂ ∈ H^-1/2  
+// -(βu , ∇v)  + (σ , ∇v)     + < f̂ ,  v  > = (f,v),   ∀ v ∈ H^1(Ω)      
+//   (u , ∇⋅τ) + 1/ε (σ , τ)  + < û , τ⋅n > = 0,       ∀ τ ∈ H(div,Ω)
+//                                        û = u_0  on ∂Ω 
+
+// Note: 
+// f̂ := βu - σ
+// û := -u
+
+// -------------------------------------------------------------
+// |   |     u     |     σ     |   û       |     f̂    |  RHS    |
+// -------------------------------------------------------------
+// | v |-(βu , ∇v) | (σ , ∇v)  |           | < f̂ ,v > |  (f,v)  |
+// |   |           |           |           |          |         |
+// | τ | (u ,∇⋅τ)  | 1/ε(σ , τ)|  <û,τ⋅n>  |          |    0    |  
+
+// where (v,τ) ∈  H^1(Ω_h) × H(div,Ω_h)  
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+enum prob_type
+{
+   polynomial,
+   EJ,   
+   general  
+};
+
+enum test_norm_type
+{
+   standard,   
+   adjoint_graph,
+   robust  
+};
+
+prob_type prob;
+test_norm_type test_norm;
+Vector beta;
+double epsilon;
+// Function returns the solution u, and gradient du and the Laplacian d2u
+void solution(const Vector & x, double & u, Vector & du, double & d2u);
+double exact_u(const Vector & X);
+void exact_sigma(const Vector & X, Vector & sigma);
+double exact_hatu(const Vector & X);
+void exact_hatf(const Vector & X, Vector & hatf);
+double f_exact(const Vector & X);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../../data/inline-quad.mesh";
+   int order = 1;
+   int delta_order = 1;
+   int ref = 1;
+   bool visualization = true;
+   int iprob = 0;
+   int itest_norm = 0;
+   double theta = 0.7;
+   epsilon = 1e0;
+
+   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(&delta_order, "-do", "--delta_order",
+                  "Order enrichment for DPG test space.");
+   args.AddOption(&epsilon, "-eps", "--epsilon",
+                  "Epsilon coefficient");               
+   args.AddOption(&ref, "-ref", "--num_refinements",
+                  "Number of uniform refinements");        
+   args.AddOption(&theta, "-theta", "--theta",
+                  "Theta parameter for AMR");                             
+   args.AddOption(&iprob, "-prob", "--problem", "Problem case"
+                  " 0: polynomial, 1: EJ ,2: General");           
+   args.AddOption(&itest_norm, "-test-norm", "--test-norm", "Choice of test norm"
+                  " 0: Standard, 1: Adjoint Graph, 2: Robust");       
+   args.AddOption(&beta, "-beta", "--beta",
+                  "Vector Coefficient beta");                                              
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   if (iprob > 2) { iprob = 2; }
+   prob = (prob_type)iprob;
+   test_norm = (test_norm_type)itest_norm;
+
+   if (prob == prob_type::EJ)
+   {
+      mesh_file = "../../../data/inline-quad.mesh";
+   }
+
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   if (beta.Size() == 0)
+   {
+      beta.SetSize(dim);
+      beta[0] = 1.;
+      beta[1] = 0.;
+   }
+
+   // Define spaces
+   // L2 space for u
+   FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
+   FiniteElementSpace *u_fes = new FiniteElementSpace(&mesh,u_fec);
+
+   // Vector L2 space for σ 
+   FiniteElementCollection *sigma_fec = new L2_FECollection(order-1,dim);
+   FiniteElementSpace *sigma_fes = new FiniteElementSpace(&mesh,sigma_fec, dim); 
+
+   // H^1/2 space for û 
+   FiniteElementCollection * hatu_fec = new H1_Trace_FECollection(order,dim);
+   FiniteElementSpace *hatu_fes = new FiniteElementSpace(&mesh,hatu_fec);
+
+   // H^-1/2 space for σ̂ 
+   FiniteElementCollection * hatf_fec = new RT_Trace_FECollection(order-1,dim);   
+   FiniteElementSpace *hatf_fes = new FiniteElementSpace(&mesh,hatf_fec);
+
+   // testspace fe collections
+   int test_order = order+delta_order;
+   FiniteElementCollection * v_fec = new H1_FECollection(test_order, dim);
+   FiniteElementCollection * tau_fec = new RT_FECollection(test_order-1, dim);
+
+   // Coefficients
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient negone(-1.0);
+   ConstantCoefficient eps(epsilon);
+   ConstantCoefficient eps1(1./epsilon);
+   ConstantCoefficient negeps1(-1./epsilon);
+   ConstantCoefficient eps2(1/(epsilon*epsilon));
+
+   ConstantCoefficient negeps(-epsilon);
+   VectorConstantCoefficient betacoeff(beta);
+   Vector negbeta = beta;
+   negbeta.Neg();
+
+   ConstantCoefficient zero(0.0);
+   Vector vec0(dim); vec0 = 0.;
+   VectorConstantCoefficient vzero(vec0);
+
+
+   DenseMatrix bbt(beta.Size());
+   MultVVt(beta, bbt); 
+   MatrixConstantCoefficient bbtcoeff(bbt);
+
+
+   VectorConstantCoefficient negbetacoeff(negbeta);
+   // Normal equation weak formulation
+   Array<FiniteElementSpace * > trial_fes; 
+   Array<FiniteElementCollection * > test_fec; 
+
+   trial_fes.Append(u_fes);
+   trial_fes.Append(sigma_fes);
+   trial_fes.Append(hatu_fes);
+   trial_fes.Append(hatf_fes);
+   test_fec.Append(v_fec);
+   test_fec.Append(tau_fec);
+
+
+   FiniteElementCollection *coeff_fec = new L2_FECollection(0,dim);
+   FiniteElementSpace *coeff_fes = new FiniteElementSpace(&mesh,coeff_fec);
+   GridFunction c1_gf, c2_gf;
+   GridFunctionCoefficient c1_coeff(&c1_gf);
+   GridFunctionCoefficient c2_coeff(&c2_gf);
+
+
+   NormalEquations * a = new NormalEquations(trial_fes,test_fec);
+   a->StoreMatrices(true);
+
+   //-(βu , ∇v)
+   a->AddTrialIntegrator(new MixedScalarWeakDivergenceIntegrator(betacoeff),0,0);
+   
+   // (σ,∇ v)
+   a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),1,0);
+
+   // (u ,∇⋅τ)
+   a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(negone),0,1);
+
+   // 1/ε (σ,τ)
+   a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(eps1)),1,1);
+
+   //  <û,τ⋅n>
+   a->AddTrialIntegrator(new NormalTraceIntegrator,2,1);
+
+   // <f̂ ,v> 
+   a->AddTrialIntegrator(new TraceIntegrator,3,0);
+
+
+   switch (test_norm)
+   {
+      case standard:
+      {   
+         // (∇v,∇δv)
+         mfem::out << "\n Test norm: Standard" << endl;
+         a->AddTestIntegrator(new DiffusionIntegrator(one),0,0);
+         // (v,δv)
+         a->AddTestIntegrator(new MassIntegrator(one),0,0);
+         // (∇⋅τ,∇⋅δτ)
+         a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+         // (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(one),1,1);
+      }
+      break;
+      case adjoint_graph:
+      {
+         mfem::out << "\n Test norm: Adjoint Graph" << endl;
+         // (∇v,∇δv)
+         a->AddTestIntegrator(new DiffusionIntegrator(one),0,0);
+         // (β⋅∇v, β⋅∇δv)   
+         a->AddTestIntegrator(new DiffusionIntegrator(bbtcoeff), 0,0);
+         // (v,δv)
+         a->AddTestIntegrator(new MassIntegrator(one),0,0);
+         // (∇⋅τ,∇⋅δτ)
+         a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+         // (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(one),1,1);
+         // 1/ε^2 (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(eps2),1,1);
+         // 1/ε (∇v, δτ)
+         a->AddTestIntegrator(new MixedVectorGradientIntegrator(eps1),0,1);
+         // - (β ⋅ ∇v,∇⋅δτ)  
+         a->AddTestIntegrator(new MixedGradDivIntegrator(betacoeff),0,1);
+         // 1/ε (τ,∇δv) 
+         a->AddTestIntegrator(new MixedVectorWeakDivergenceIntegrator(negeps1),1,0);
+         // -(β ∇⋅τ ,∇⋅δv)  
+         a->AddTestIntegrator(new MixedDivGradIntegrator(betacoeff),1,0);
+      }
+      break;
+      default:
+      {
+         mfem::out << "\n Test norm: Robust" << endl;
+         c1_gf.SetSpace(coeff_fes);
+         c2_gf.SetSpace(coeff_fes);
+         Array<int> dofs;
+         for (int i =0; i < mesh.GetNE(); i++)
+         {
+            double volume = mesh.GetElementVolume(i);
+            double c1 = min(epsilon/volume, 1.);
+            double c2 = min(1./epsilon, 1./volume);
+            // double c2 = 1.;
+            coeff_fes->GetElementDofs(i,dofs);
+            c1_gf.SetSubVector(dofs,c1);
+            c2_gf.SetSubVector(dofs,c2);
+         }
+         // c1 (v,δv)
+         a->AddTestIntegrator(new MassIntegrator(c1_coeff),0,0);
+         // ε (∇v,∇δv)
+         a->AddTestIntegrator(new DiffusionIntegrator(eps),0,0);
+         // (β⋅∇v, β⋅∇δv)   
+         a->AddTestIntegrator(new DiffusionIntegrator(bbtcoeff), 0,0);
+         // c2 (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(c2_coeff),1,1);
+         // (∇⋅τ,∇⋅δτ)
+         a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+      }
+      break;
+   }
+
+
+   FunctionCoefficient f(f_exact);
+   // if (prob != prob_type::EJ)
+   // {
+      a->AddDomainLFIntegrator(new DomainLFIntegrator(f),0);
+   // }
+
+   FunctionCoefficient hatuex(exact_hatu);
+   VectorFunctionCoefficient hatfex(dim,exact_hatf);
+   Array<int> elements_to_refine;
+   FunctionCoefficient uex(exact_u);
+   VectorFunctionCoefficient sigmaex(dim,exact_sigma);
+   GridFunction hatu_gf;
+   GridFunction hatf_gf;
+
+   // socketstream uex_out;
+   socketstream u_out;
+   // socketstream sigma_out;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      u_out.open(vishost, visport);
+      // uex_out.open(vishost, visport);
+      // sigma_out.open(vishost, visport);
+   }
+
+   double res0 = 0.;
+   double err0 = 0.;
+   int dof0;
+   mfem::out << " Refinement |" 
+             << "    Dofs    |" 
+             << "  L2 Error  |" 
+             << " Relative % |" 
+             << "  Rate  |" 
+             << "  Residual  |" 
+             << "  Rate  |" << endl;
+   mfem::out << " --------------------"      
+             <<  "-------------------"    
+             <<  "-------------------"    
+             <<  "-------------------" << endl;   
+
+
+   for (int i = 0; i<=ref; i++)
+   {
+      a->Assemble();
+
+      Array<int> ess_tdof_list_uhat;
+      Array<int> ess_tdof_list_fhat;
+      Array<int> ess_bdr_uhat;
+      Array<int> ess_bdr_fhat;
+      if (mesh.bdr_attributes.Size())
+      {
+         ess_bdr_uhat.SetSize(mesh.bdr_attributes.Max());
+         ess_bdr_fhat.SetSize(mesh.bdr_attributes.Max());
+         // ess_bdr_uhat = 1;
+         // ess_bdr_fhat = 0;
+         ess_bdr_uhat = 0;
+         ess_bdr_fhat = 1;
+         ess_bdr_uhat[1] = 1;
+         ess_bdr_fhat[1] = 0;
+         hatu_fes->GetEssentialTrueDofs(ess_bdr_uhat, ess_tdof_list_uhat);
+         hatf_fes->GetEssentialTrueDofs(ess_bdr_fhat, ess_tdof_list_fhat);
+      }
+
+      // shift the ess_tdofs
+      int n = ess_tdof_list_uhat.Size();
+      int m = ess_tdof_list_fhat.Size();
+      Array<int> ess_tdof_list(n+m);
+      for (int i = 0; i < n; i++)
+      {
+         ess_tdof_list[i] = ess_tdof_list_uhat[i] 
+                          + u_fes->GetTrueVSize() 
+                          + sigma_fes->GetTrueVSize();
+      }
+      for (int i = 0; i < m; i++)
+      {
+         ess_tdof_list[i+n] = ess_tdof_list_fhat[i] 
+                            + u_fes->GetTrueVSize() 
+                            + sigma_fes->GetTrueVSize()
+                            + hatu_fes->GetTrueVSize();
+      }
+
+      Array<int> offsets(5);
+      offsets[0] = 0;
+      offsets[1] = u_fes->GetVSize();
+      offsets[2] = sigma_fes->GetVSize();
+      offsets[3] = hatu_fes->GetVSize();
+      offsets[4] = hatf_fes->GetVSize();
+      offsets.PartialSum();
+      BlockVector x(offsets);
+      x = 0.0;
+      hatu_gf.MakeRef(hatu_fes,x.GetBlock(2));
+
+      hatf_gf.MakeRef(hatf_fes,x.GetBlock(3));
+
+      hatu_gf.ProjectBdrCoefficient(hatuex,ess_bdr_uhat);
+      hatf_gf.ProjectBdrCoefficientNormal(hatfex,ess_bdr_fhat);
+
+      OperatorPtr Ah;
+      Vector X,B;
+      a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
+
+      BlockMatrix * A = Ah.As<BlockMatrix>();
+
+      BlockDiagonalPreconditioner * M = new BlockDiagonalPreconditioner(A->RowOffsets());
+      M->owns_blocks = 1;
+      for (int i=0; i<A->NumRowBlocks(); i++)
+      {
+         M->SetDiagonalBlock(i,new UMFPackSolver(A->GetBlock(i,i)));
+      }
+
+      CGSolver cg;
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(200000);
+      cg.SetPrintLevel(0);
+      cg.SetPreconditioner(*M);
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+      delete M;
+
+      a->RecoverFEMSolution(X,x);
+      Vector & residuals = a->ComputeResidual(x);
+
+      double residual = residuals.Norml2();
+
+      elements_to_refine.SetSize(0);
+      double max_resid = residuals.Max();
+      for (int iel = 0; iel<mesh.GetNE(); iel++)
+      {
+         if (residuals[iel] > theta * max_resid)
+         {
+            elements_to_refine.Append(iel);
+         }
+      }   
+
+      GridFunction uex_gf(u_fes);
+      uex_gf.ProjectCoefficient(uex);
+
+      GridFunction sigmaex_gf(sigma_fes);
+      sigmaex_gf.ProjectCoefficient(sigmaex);
+
+      GridFunction u_gf;
+      u_gf.MakeRef(u_fes,x.GetBlock(0));
+
+      GridFunction sigma_gf;
+      sigma_gf.MakeRef(sigma_fes,x.GetBlock(1));
+
+      int dofs = X.Size();
+      double u_err = u_gf.ComputeL2Error(uex);
+      double u_norm = uex_gf.ComputeL2Error(zero);
+      // mfem::out << "u_err = " << u_err << endl;
+      double sigma_err = sigma_gf.ComputeL2Error(sigmaex);
+      double sigma_norm = sigmaex_gf.ComputeL2Error(vzero);
+      // mfem::out << "sigma_err = " << sigma_err << endl;
+      double L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
+      double L2norm = sqrt(u_norm * u_norm + sigma_norm * sigma_norm);
+
+      double rel_error = L2Error/L2norm;
+
+      double rate_err = (i) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0;
+      double rate_res = (i) ? dim*log(res0/residual)/log((double)dof0/dofs) : 0.0;
+
+      err0 = L2Error;
+      res0 = residual;
+      dof0 = dofs;
+      mfem::out << std::right << std::setw(11) << i << " | " 
+                << std::setw(10) <<  dof0 << " | " 
+                << std::setprecision(3) 
+                << std::setw(10) << std::scientific <<  err0 << " | " 
+                << std::setprecision(3) 
+                << std::setw(10) << std::fixed <<  rel_error * 100. << " | " 
+                << std::setprecision(2) 
+                << std::setw(6) << std::fixed << rate_err << " | " 
+                << std::setprecision(3) 
+                << std::setw(10) << std::scientific <<  res0 << " | " 
+                << std::setprecision(2) 
+                << std::setw(6) << std::fixed << rate_res << " | " 
+                << std::resetiosflags(std::ios::showbase)
+                << std::endl;
+
+
+      if (visualization)
+      {
+         // uex_out.precision(8);
+         // uex_out << "solution\n" << mesh << uex_gf <<
+         //          "window_title 'Exact u' "
+         //          << flush;
+         u_out.precision(8);
+         u_out << "solution\n" << mesh << u_gf <<
+                  "window_title 'Numerical u' "
+                  << flush;
+         // sigma_out.precision(8);
+         // sigma_out << "solution\n" << mesh << sigma_gf <<
+         //       "window_title 'Numerical flux' "
+         //       << flush;
+      }
+
+      if (i == ref)
+         break;
+
+      mesh.GeneralRefinement(elements_to_refine,1,1);
+      for (int i =0; i<trial_fes.Size(); i++)
+      {
+         trial_fes[i]->Update(false);
+      }
+      a->Update();
+
+      if (test_norm == test_norm_type::robust)
+      {
+         coeff_fes->Update();
+         c1_gf.Update();
+         c2_gf.Update();
+         Array<int> dofs;
+         for (int i = 0; i < mesh.GetNE(); i++)
+         {
+            double volume = mesh.GetElementVolume(i);
+            double c1 = min(epsilon/volume, 1.);
+            double c2 = min(1./epsilon, 1./volume);
+            // double c2 = 1.;
+            coeff_fes->GetElementDofs(i,dofs);
+            c1_gf.SetSubVector(dofs,c1);
+            c2_gf.SetSubVector(dofs,c2);
+         }
+      }
+
+   }
+
+   delete a;
+   delete tau_fec;
+   delete v_fec;
+   delete hatf_fes;
+   delete hatf_fec;
+   delete hatu_fes;
+   delete hatu_fec;
+   delete sigma_fec;
+   delete u_fec;
+   delete u_fes;
+
+   return 0;
+}
+
+
+void solution(const Vector & X, double & u, Vector & du, double & d2u)
+{
+   double x = X[0];
+   double y = X[1];
+   double z = 0.;
+   if (X.Size() == 3) z = X[2];
+   du.SetSize(X.Size());
+   du = 0.;
+   d2u = 0.;
+
+   switch(prob)
+   {
+      case polynomial:
+      {
+         int n=2;
+         int m=2;
+         u = pow(x,n)*pow(y,m);
+         du[0] = n * pow(x,n-1) * pow(y,m);
+         du[1] = m * pow(x,n) * pow(y,m-1);
+         d2u = n * (n-1) * pow(x,n-2) * pow(y,m)
+             + m * (m-1) * pow(x,n) * pow(y,m-2);
+      }
+      break;
+      case EJ:
+      {
+         double alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
+         double r1 = (1. + alpha) / (2.*epsilon);
+         double r2 = (1. - alpha) / (2.*epsilon);
+         double denom = exp(-r2) - exp(-r1);
+
+
+         double g1 = exp(r2*(x-1.));
+         double g1_x = r2*g1;
+         double g1_xx = r2*g1_x;
+         double g2 = exp(r1*(x-1.));
+         double g2_x = r1*g2;
+         double g2_xx = r1*g2_x;
+         double g = g1-g2;
+         double g_x = g1_x - g2_x;
+         double g_xx = g1_xx - g2_xx;
+
+
+         u = g * cos(M_PI * y)/denom;
+         double u_x = g_x * cos(M_PI * y)/denom;
+         double u_xx = g_xx * cos(M_PI * y)/denom;
+         double u_y = -M_PI * g * sin(M_PI*y)/denom;
+         double u_yy = -M_PI * M_PI * u;
+         du[0] = u_x;
+         du[1] = u_y;
+         d2u = u_xx + u_yy;
+
+      }
+      break;
+      default:
+      {
+         double alpha = M_PI * (x + y + z);
+         u = sin(alpha);
+         du.SetSize(X.Size());
+         for (int i = 0; i<du.Size(); i++)
+         {
+            du[i] = M_PI * cos(alpha);
+         }
+         d2u = - M_PI*M_PI * u * du.Size();
+      }
+      break;
+   }
+}
+
+
+double exact_u(const Vector & X)
+{
+   double u, d2u;
+   Vector du;
+   solution(X,u,du,d2u);
+   return u;
+}
+
+void exact_sigma(const Vector & X, Vector & sigma)
+{
+   double u, d2u;
+   Vector du;
+   solution(X,u,du,d2u);
+   // σ = ε ∇ u
+   sigma = du;
+   sigma *= epsilon;
+}
+
+double exact_hatu(const Vector & X)
+{
+   return -exact_u(X);
+}
+
+void exact_hatf(const Vector & X, Vector & hatf)
+{
+   Vector sigma;
+   exact_sigma(X,sigma);
+   double u = exact_u(X);
+   hatf.SetSize(X.Size());
+   for (int i = 0; i<hatf.Size(); i++)
+   {
+      hatf[i] = beta[i] * u - sigma[i];
+   }
+}
+
+double f_exact(const Vector & X)
+{
+   // f = - εΔu + ∇⋅(βu)
+   double u, d2u;
+   Vector du;
+   solution(X,u,du,d2u);
+
+   double s = 0;
+   for (int i = 0; i<du.Size(); i++)
+   {
+      s += beta[i] * du[i];
+   }
+   return -epsilon * d2u + s;
+}
@@ -0,0 +1,698 @@
+//                                MFEM Ultraweak DPG example
+//
+// Compile with: make  uw_dpgp
+//
+// sample runs 
+// mpirun -np 6 ./uw_dpgp  -m ../../../data/inline-quad.mesh -o 3 -ref 10 -test-norm 2 -do 1 -prob 1 -eps 1e-4
+//     - εΔu + ∇⋅(βu) = f,   in Ω
+//                  u = u_0, on ∂Ω
+
+// First Order System
+
+//     - ∇⋅σ + ∇⋅(βu) = f,   in Ω
+//        1/ε σ - ∇u  = 0,   in Ω
+//                  u = u_0, on ∂Ω
+
+// UW-DPG:
+// 
+// u ∈ L^2(Ω), σ ∈ (L^2(Ω))^dim 
+// û ∈ H^1/2, f̂ ∈ H^-1/2  
+// -(βu , ∇v)  + (σ , ∇v)     + < f̂ ,  v  > = (f,v),   ∀ v ∈ H^1(Ω)      
+//   (u , ∇⋅τ) + 1/ε (σ , τ)  + < û , τ⋅n > = 0,       ∀ τ ∈ H(div,Ω)
+//                                        û = u_0  on ∂Ω 
+
+// Note: 
+// f̂ := βu - σ
+// û := -u
+
+// -------------------------------------------------------------
+// |   |     u     |     σ     |   û       |     f̂    |  RHS    |
+// -------------------------------------------------------------
+// | v |-(βu , ∇v) | (σ , ∇v)  |           | < f̂ ,v > |  (f,v)  |
+// |   |           |           |           |          |         |
+// | τ | (u ,∇⋅τ)  | 1/ε(σ , τ)|  <û,τ⋅n>  |          |    0    |  
+
+// where (v,τ) ∈  H^1(Ω_h) × H(div,Ω_h)  
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+enum prob_type
+{
+   polynomial,
+   EJ,   
+   general  
+};
+
+enum test_norm_type
+{
+   standard,   
+   adjoint_graph,
+   robust  
+};
+
+prob_type prob;
+test_norm_type test_norm;
+Vector beta;
+double epsilon;
+// Function returns the solution u, and gradient du and the Laplacian d2u
+void solution(const Vector & x, double & u, Vector & du, double & d2u);
+double exact_u(const Vector & X);
+void exact_sigma(const Vector & X, Vector & sigma);
+double exact_hatu(const Vector & X);
+void exact_hatf(const Vector & X, Vector & hatf);
+double f_exact(const Vector & X);
+
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi;
+   int num_procs = mpi.WorldSize();
+   int myid = mpi.WorldRank();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../../data/inline-quad.mesh";
+   int order = 1;
+   int delta_order = 1;
+   int ref = 1;
+   bool visualization = true;
+   int iprob = 0;
+   int itest_norm = 0;
+   double theta = 0.7;
+   bool static_cond = false;
+   epsilon = 1e0;
+
+   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(&delta_order, "-do", "--delta_order",
+                  "Order enrichment for DPG test space.");
+   args.AddOption(&epsilon, "-eps", "--epsilon",
+                  "Epsilon coefficient");               
+   args.AddOption(&ref, "-ref", "--num_refinements",
+                  "Number of uniform refinements");     
+   args.AddOption(&theta, "-theta", "--theta",
+                  "Theta parameter for AMR");                           
+   args.AddOption(&iprob, "-prob", "--problem", "Problem case"
+                  " 0: lshape, 1: General");                         
+   args.AddOption(&itest_norm, "-test-norm", "--test-norm", "Choice of test norm"
+                  " 0: Standard, 1: Adjoint Graph, 2: Robust");     
+   args.AddOption(&beta, "-beta", "--beta",
+                  "Vector Coefficient beta");          
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");                            
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   if (iprob > 2) { iprob = 2; }
+   prob = (prob_type)iprob;
+
+   test_norm = (test_norm_type)itest_norm;
+
+   if (prob == prob_type::EJ)
+   {
+      mesh_file = "../../../data/inline-quad.mesh";
+   }
+
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   if (beta.Size() == 0)
+   {
+      beta.SetSize(dim);
+      beta[0] = 1.;
+      beta[1] = 0.;
+   }
+
+   mesh.EnsureNCMesh();
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // Define spaces
+   // L2 space for u
+   FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
+   ParFiniteElementSpace *u_fes = new ParFiniteElementSpace(&pmesh,u_fec);
+
+   // Vector L2 space for σ 
+   FiniteElementCollection *sigma_fec = new L2_FECollection(order-1,dim);
+   ParFiniteElementSpace *sigma_fes = new ParFiniteElementSpace(&pmesh,sigma_fec, dim); 
+
+   // H^1/2 space for û 
+   FiniteElementCollection * hatu_fec = new H1_Trace_FECollection(order,dim);
+   ParFiniteElementSpace *hatu_fes = new ParFiniteElementSpace(&pmesh,hatu_fec);
+
+   // H^-1/2 space for σ̂ 
+   FiniteElementCollection * hatf_fec = new RT_Trace_FECollection(order-1,dim);   
+   ParFiniteElementSpace *hatf_fes = new ParFiniteElementSpace(&pmesh,hatf_fec);
+
+   // testspace fe collections
+   int test_order = order+delta_order;
+   FiniteElementCollection * v_fec = new H1_FECollection(test_order, dim);
+   FiniteElementCollection * tau_fec = new RT_FECollection(test_order-1, dim);
+
+   // Coefficients
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient negone(-1.0);
+   ConstantCoefficient eps(epsilon);
+   ConstantCoefficient eps1(1./epsilon);
+   ConstantCoefficient negeps1(-1./epsilon);
+   ConstantCoefficient eps2(1/(epsilon*epsilon));
+
+   ConstantCoefficient negeps(-epsilon);
+   VectorConstantCoefficient betacoeff(beta);
+   Vector negbeta = beta;
+   negbeta.Neg();
+
+   DenseMatrix bbt(beta.Size());
+   MultVVt(beta, bbt); 
+   MatrixConstantCoefficient bbtcoeff(bbt);
+
+
+   VectorConstantCoefficient negbetacoeff(negbeta);
+   // Normal equation weak formulation
+   Array<ParFiniteElementSpace * > trial_fes; 
+   Array<FiniteElementCollection * > test_fec; 
+
+   trial_fes.Append(u_fes);
+   trial_fes.Append(sigma_fes);
+   trial_fes.Append(hatu_fes);
+   trial_fes.Append(hatf_fes);
+   test_fec.Append(v_fec);
+   test_fec.Append(tau_fec);
+
+   ParNormalEquations * a = new ParNormalEquations(trial_fes,test_fec);
+   a->StoreMatrices(true);
+
+   //-(βu , ∇v)
+   a->AddTrialIntegrator(new MixedScalarWeakDivergenceIntegrator(betacoeff),0,0);
+   
+   // (σ,∇ v)
+   a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),1,0);
+
+   // (u ,∇⋅τ)
+   a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(negone),0,1);
+
+   // 1/ε (σ,τ)
+   a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(eps1)),1,1);
+
+   //  <û,τ⋅n>
+   a->AddTrialIntegrator(new NormalTraceIntegrator,2,1);
+
+   // <f̂ ,v> 
+   a->AddTrialIntegrator(new TraceIntegrator,3,0);
+
+
+   FiniteElementCollection *coeff_fec = new L2_FECollection(0,dim);
+   ParFiniteElementSpace *coeff_fes = new ParFiniteElementSpace(&pmesh,coeff_fec);
+   ParGridFunction c1_gf, c2_gf;
+   GridFunctionCoefficient c1_coeff(&c1_gf);
+   GridFunctionCoefficient c2_coeff(&c2_gf);
+
+   switch (test_norm)
+   {
+      case standard:
+      {   
+         if (myid == 0)
+         {
+            mfem::out << "\n Test norm: Standard" << endl;
+         }
+         // (∇v,∇δv)
+         a->AddTestIntegrator(new DiffusionIntegrator(one),0,0);
+         // (v,δv)
+         a->AddTestIntegrator(new MassIntegrator(one),0,0);
+         // (∇⋅τ,∇⋅δτ)
+         a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+         // (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(one),1,1);
+      }
+      break;
+      case adjoint_graph:
+      {
+         if (myid == 0)
+         {
+            mfem::out << "\n Test norm: Adjoint Graph" << endl;
+         }
+         // (∇v,∇δv)
+         a->AddTestIntegrator(new DiffusionIntegrator(one),0,0);
+         // (β⋅∇v, β⋅∇δv)   
+         a->AddTestIntegrator(new DiffusionIntegrator(bbtcoeff), 0,0);
+         // (v,δv)
+         a->AddTestIntegrator(new MassIntegrator(one),0,0);
+         // (∇⋅τ,∇⋅δτ)
+         a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+         // (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(one),1,1);
+         // 1/ε^2 (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(eps2),1,1);
+         // 1/ε (∇v, δτ)
+         a->AddTestIntegrator(new MixedVectorGradientIntegrator(eps1),0,1);
+         // - (β ⋅ ∇v,∇⋅δτ)  
+         a->AddTestIntegrator(new MixedGradDivIntegrator(betacoeff),0,1);
+         // 1/ε (τ,∇δv) 
+         a->AddTestIntegrator(new MixedVectorWeakDivergenceIntegrator(negeps1),1,0);
+         // -(β ∇⋅τ ,∇⋅δv)  
+         a->AddTestIntegrator(new MixedDivGradIntegrator(betacoeff),1,0);
+      }
+      break;
+      default:
+      {
+         if (myid == 0)
+         {
+            mfem::out << "\n Test norm: Robust" << endl;
+         }
+         c1_gf.SetSpace(coeff_fes);
+         c2_gf.SetSpace(coeff_fes);
+         Array<int> dofs;
+         for (int i =0; i < pmesh.GetNE(); i++)
+         {
+            double volume = pmesh.GetElementVolume(i);
+            double c1 = min(epsilon/volume, 1.);
+            double c2 = min(1./epsilon, 1./volume);
+            coeff_fes->GetElementDofs(i,dofs);
+            c1_gf.SetSubVector(dofs,c1);
+            c2_gf.SetSubVector(dofs,c2);
+         }
+         // c1 (v,δv)
+         a->AddTestIntegrator(new MassIntegrator(c1_coeff),0,0);
+         // ε (∇v,∇δv)
+         a->AddTestIntegrator(new DiffusionIntegrator(eps),0,0);
+         // (β⋅∇v, β⋅∇δv)   
+         a->AddTestIntegrator(new DiffusionIntegrator(bbtcoeff), 0,0);
+         // c2 (τ,δτ)
+         a->AddTestIntegrator(new VectorFEMassIntegrator(c2_coeff),1,1);
+         // (∇⋅τ,∇⋅δτ)
+         a->AddTestIntegrator(new DivDivIntegrator(one),1,1);
+      }
+      break;
+   }
+
+
+   FunctionCoefficient f(f_exact);
+   // if (prob != prob_type::EJ)
+   // {
+   a->AddDomainLFIntegrator(new DomainLFIntegrator(f),0);
+   // }
+
+   FunctionCoefficient hatuex(exact_hatu);
+   VectorFunctionCoefficient hatfex(dim,exact_hatf);
+   Array<int> elements_to_refine;
+   FunctionCoefficient uex(exact_u);
+   VectorFunctionCoefficient sigmaex(dim,exact_sigma);
+
+   ParGridFunction hatu_gf;
+   ParGridFunction hatf_gf;
+
+   // socketstream uex_out;
+   socketstream u_out;
+   // socketstream sigma_out;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      u_out.open(vishost, visport);
+      // uex_out.open(vishost, visport);
+      // sigma_out.open(vishost, visport);
+   }
+
+   double res0 = 0.;
+   double err0 = 0.;
+   int dof0;
+   if (myid == 0)
+   {
+      mfem::out << " Refinement |" 
+               << "    Dofs    |" 
+               << "  L2 Error  |" 
+               << "  Rate  |" 
+               << "  Residual  |" 
+               << "  Rate  |" << endl;
+      mfem::out << " --------------------"      
+               <<  "-------------------"    
+               <<  "-------------------"    
+               <<  "-------------------" << endl;   
+   }
+
+
+   for (int i = 0; i<=ref; i++)
+   {
+      if (static_cond) { a->EnableStaticCondensation(); }
+      a->Assemble();
+
+      Array<int> ess_tdof_list_uhat;
+      Array<int> ess_tdof_list_fhat;
+      Array<int> ess_bdr_uhat;
+      Array<int> ess_bdr_fhat;
+      if (pmesh.bdr_attributes.Size())
+      {
+         ess_bdr_uhat.SetSize(pmesh.bdr_attributes.Max());
+         ess_bdr_fhat.SetSize(pmesh.bdr_attributes.Max());
+         // ess_bdr_uhat = 1;
+         // ess_bdr_fhat = 0;
+         ess_bdr_uhat = 0;
+         ess_bdr_fhat = 1;
+         ess_bdr_uhat[1] = 1;
+         ess_bdr_fhat[1] = 0;
+         hatu_fes->GetEssentialTrueDofs(ess_bdr_uhat, ess_tdof_list_uhat);
+         hatf_fes->GetEssentialTrueDofs(ess_bdr_fhat, ess_tdof_list_fhat);
+      }
+
+      // shift the ess_tdofs
+      int n = ess_tdof_list_uhat.Size();
+      int m = ess_tdof_list_fhat.Size();
+      Array<int> ess_tdof_list(n+m);
+      for (int i = 0; i < n; i++)
+      {
+         ess_tdof_list[i] = ess_tdof_list_uhat[i] 
+                          + u_fes->GetTrueVSize() 
+                          + sigma_fes->GetTrueVSize();
+      }
+      for (int i = 0; i < m; i++)
+      {
+         ess_tdof_list[i+n] = ess_tdof_list_fhat[i] 
+                            + u_fes->GetTrueVSize() 
+                            + sigma_fes->GetTrueVSize()
+                            + hatu_fes->GetTrueVSize();
+      }
+
+      Array<int> offsets(5);
+      offsets[0] = 0;
+      offsets[1] = u_fes->GetVSize();
+      offsets[2] = sigma_fes->GetVSize();
+      offsets[3] = hatu_fes->GetVSize();
+      offsets[4] = hatf_fes->GetVSize();
+      offsets.PartialSum();
+      BlockVector x(offsets);
+      x = 0.0;
+      hatu_gf.MakeRef(hatu_fes,x.GetBlock(2));
+      hatu_gf.ProjectBdrCoefficient(hatuex,ess_bdr_uhat);
+
+      hatf_gf.MakeRef(hatf_fes,x.GetBlock(3));
+      hatf_gf.ProjectBdrCoefficientNormal(hatfex,ess_bdr_fhat);
+
+      OperatorPtr Ah;
+      Vector X,B;
+      a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
+
+      BlockOperator * A = Ah.As<BlockOperator>();
+
+
+      BlockDiagonalPreconditioner * M = new BlockDiagonalPreconditioner(A->RowOffsets());
+      M->owns_blocks = 1;
+      int skip = 0;
+      if (!static_cond)
+      {
+         HypreBoomerAMG * amg0 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(0,0));
+         HypreBoomerAMG * amg1 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(1,1));
+         amg0->SetPrintLevel(0);
+         amg1->SetPrintLevel(0);
+         M->SetDiagonalBlock(0,amg0);
+         M->SetDiagonalBlock(1,amg1);
+         skip = 2;
+      }
+      HypreBoomerAMG * amg2 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(skip,skip));
+      amg2->SetPrintLevel(0);
+      M->SetDiagonalBlock(skip,amg2);
+
+      HypreSolver * prec;
+      if (dim == 2)
+      {
+         prec = new HypreAMS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatf_fes);
+      }
+      else
+      {
+         prec = new HypreADS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatf_fes);
+      }
+      M->SetDiagonalBlock(skip+1,prec);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-6);
+      cg.SetMaxIter(200000);
+      cg.SetPrintLevel(0);
+      cg.SetPreconditioner(*M);
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+      delete M;
+
+      a->RecoverFEMSolution(X,x);
+      Vector & residuals = a->ComputeResidual(x);
+
+      double residual = residuals.Norml2();
+      double maxresidual = residuals.Max(); 
+
+      double gresidual = residual * residual; 
+
+      MPI_Allreduce(MPI_IN_PLACE,&maxresidual,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
+      MPI_Allreduce(MPI_IN_PLACE,&gresidual,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+
+      gresidual = sqrt(gresidual);
+
+      elements_to_refine.SetSize(0);
+      for (int iel = 0; iel<pmesh.GetNE(); iel++)
+      {
+         if (residuals[iel] > theta * maxresidual)
+         {
+            elements_to_refine.Append(iel);
+         }
+      }
+
+      ParGridFunction u_gf;
+      u_gf.MakeRef(u_fes,x.GetBlock(0));
+
+      ParGridFunction sigma_gf;
+      sigma_gf.MakeRef(sigma_fes,x.GetBlock(1));
+
+      int dofs = u_fes->GlobalTrueVSize()
+               + sigma_fes->GlobalTrueVSize()
+               + hatu_fes->GlobalTrueVSize()
+               + hatf_fes->GlobalTrueVSize();
+
+      double u_err = u_gf.ComputeL2Error(uex);
+      double sigma_err = sigma_gf.ComputeL2Error(sigmaex);
+      double L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
+
+      double rate_err = (i) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0;
+      double rate_res = (i) ? dim*log(res0/gresidual)/log((double)dof0/dofs) : 0.0;
+
+      err0 = L2Error;
+      res0 = gresidual;
+      dof0 = dofs;
+
+      if (myid == 0)
+      {
+         mfem::out << std::right << std::setw(11) << i << " | " 
+                  << std::setw(10) <<  dof0 << " | " 
+                  << std::setprecision(3) 
+                  << std::setw(10) << std::scientific <<  err0 << " | " 
+                  << std::setprecision(2) 
+                  << std::setw(6) << std::fixed << rate_err << " | " 
+                  << std::setprecision(3) 
+                  << std::setw(10) << std::scientific <<  res0 << " | " 
+                  << std::setprecision(2) 
+                  << std::setw(6) << std::fixed << rate_res << " | " 
+                  << std::resetiosflags(std::ios::showbase)
+                  << std::endl;
+      }
+
+
+      if (visualization)
+      {
+         // uex_out.precision(8);
+         // uex_out << "parallel " << num_procs << " " << myid << "\n";
+         // uex_out << "solution\n" << pmesh << uex_gf <<
+         //          "window_title 'Exact u' "
+         //          << flush;
+
+         u_out << "parallel " << num_procs << " " << myid << "\n";
+         u_out.precision(8);
+         u_out << "solution\n" << pmesh << u_gf <<
+                  "window_title 'Numerical u' "
+                  << flush;
+
+         // sigma_out << "parallel " << num_procs << " " << myid << "\n";
+         // sigma_out.precision(8);
+         // sigma_out << "solution\n" << pmesh << sigma_gf <<
+         //       "window_title 'Numerical flux' "
+         //       << flush;
+      }
+
+      if (i == ref)
+         break;
+
+      pmesh.GeneralRefinement(elements_to_refine,1,1);
+      for (int i =0; i<trial_fes.Size(); i++)
+      {
+         trial_fes[i]->Update(false);
+      }
+      a->Update();
+
+      if (test_norm == test_norm_type::robust)
+      {
+         coeff_fes->Update();
+         c1_gf.Update();
+         c2_gf.Update();
+         Array<int> edofs;
+         for (int i = 0; i < pmesh.GetNE(); i++)
+         {
+            double volume = pmesh.GetElementVolume(i);
+            double c1 = min(epsilon/volume, 1.);
+            double c2 = min(1./epsilon, 1./volume);
+            coeff_fes->GetElementDofs(i,edofs);
+            c1_gf.SetSubVector(edofs,c1);
+            c2_gf.SetSubVector(edofs,c2);
+         }
+      }
+   }
+
+
+   delete a;
+   delete tau_fec;
+   delete v_fec;
+   delete hatf_fes;
+   delete hatf_fec;
+   delete hatu_fes;
+   delete hatu_fec;
+   delete sigma_fec;
+   delete u_fec;
+   delete u_fes;
+
+   return 0;
+}
+
+void solution(const Vector & X, double & u, Vector & du, double & d2u)
+{
+   double x = X[0];
+   double y = X[1];
+   double z = 0.;
+   if (X.Size() == 3) z = X[2];
+   du.SetSize(X.Size());
+   du = 0.;
+   d2u = 0.;
+
+   switch(prob)
+   {
+      case polynomial:
+      {
+         int n=2;
+         int m=2;
+         u = pow(x,n)*pow(y,m);
+         du[0] = n * pow(x,n-1) * pow(y,m);
+         du[1] = m * pow(x,n) * pow(y,m-1);
+         d2u = n * (n-1) * pow(x,n-2) * pow(y,m)
+             + m * (m-1) * pow(x,n) * pow(y,m-2);
+      }
+      break;
+      case EJ:
+      {
+         double alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
+         double r1 = (1. + alpha) / (2.*epsilon);
+         double r2 = (1. - alpha) / (2.*epsilon);
+         double denom = exp(-r2) - exp(-r1);
+
+
+         double g1 = exp(r2*(x-1.));
+         double g1_x = r2*g1;
+         double g1_xx = r2*g1_x;
+         double g2 = exp(r1*(x-1.));
+         double g2_x = r1*g2;
+         double g2_xx = r1*g2_x;
+         double g = g1-g2;
+         double g_x = g1_x - g2_x;
+         double g_xx = g1_xx - g2_xx;
+
+
+         u = g * cos(M_PI * y)/denom;
+         double u_x = g_x * cos(M_PI * y)/denom;
+         double u_xx = g_xx * cos(M_PI * y)/denom;
+         double u_y = -M_PI * g * sin(M_PI*y)/denom;
+         double u_yy = -M_PI * M_PI * u;
+         du[0] = u_x;
+         du[1] = u_y;
+         d2u = u_xx + u_yy;
+
+      }
+      break;
+      default:
+      {
+         double alpha = M_PI * (x + y + z);
+         u = sin(alpha);
+         du.SetSize(X.Size());
+         for (int i = 0; i<du.Size(); i++)
+         {
+            du[i] = M_PI * cos(alpha);
+         }
+         d2u = - M_PI*M_PI * u * du.Size();
+      }
+      break;
+   }
+}
+
+
+double exact_u(const Vector & X)
+{
+   double u, d2u;
+   Vector du;
+   solution(X,u,du,d2u);
+   return u;
+}
+
+void exact_sigma(const Vector & X, Vector & sigma)
+{
+   double u, d2u;
+   Vector du;
+   solution(X,u,du,d2u);
+   // σ = ε ∇ u
+   sigma = du;
+   sigma *= epsilon;
+}
+
+double exact_hatu(const Vector & X)
+{
+   return -exact_u(X);
+}
+
+void exact_hatf(const Vector & X, Vector & hatf)
+{
+   Vector sigma;
+   exact_sigma(X,sigma);
+   double u = exact_u(X);
+   hatf.SetSize(X.Size());
+   for (int i = 0; i<hatf.Size(); i++)
+   {
+      hatf[i] = beta[i] * u - sigma[i];
+   }
+}
+
+double f_exact(const Vector & X)
+{
+   // f = - εΔu + ∇⋅(βu)
+   double u, d2u;
+   Vector du;
+   solution(X,u,du,d2u);
+
+   double s = 0;
+   for (int i = 0; i<du.Size(); i++)
+   {
+      s += beta[i] * du[i];
+   }
+   return -epsilon * d2u + s;
+}
--- a/Show More
+++ b/Show More