Modify ex26 with anisotropic diffusion and ILU smoothing

2022-02-15 12:07:23 -08:00
826 changed files with 7501 additions and 39718 deletions
@@ -5,7 +5,7 @@ channels:
 dependencies:
 - xeus-cling=0.13.0
 - xwidgets=0.26.0
-# NOTE: it is possible that these are not needed for the lab frontend
+# NOTE: it's possible these aren't needed for the lab frontend
 - widgetsnbextension=3.5.1
 - pip
 - pip:
@@ -1,48 +0,0 @@
-                    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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -166,7 +166,7 @@ jobs:

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

 examples/refined.mesh
 examples/displaced.mesh
@@ -78,7 +78,6 @@ 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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -20,8 +20,8 @@
    - |
      (
        date
-        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
-        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
+        # try to get an excusive 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 "Acquired lock on '$PWD/autotest.lock'"
+        echo "Aquired 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 acquire lock on '$PWD/autotest.lock' ..."
-        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
+        # try to get an excusive 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 "Acquired lock on '$PWD/autotest.lock'"
+        echo "Aquired 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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -45,8 +45,8 @@ setup_baseline:
    - |
      (
        date
-        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
-        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
+        # try to get an excusive 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 "Acquired lock on '$PWD/autotest.lock'"
+        echo "Aquired 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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -34,8 +34,8 @@ setup:
    - |
      (
        date
-        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
+        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
        # 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 "Acquired lock on '$PWD/mfem-data.lock'"
+        echo "Aquired 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 acquire lock on '$PWD/autotest.lock' ..."
-        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
+        # try to get an excusive 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 "Acquired lock on '$PWD/autotest.lock'"
+        echo "Aquired 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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -64,8 +64,8 @@ report_baseline:
    - |
      (
        date
-        echo "Waiting to acquire lock on '$PWD/autotest.lock' ..."
-        # try to get an exclusive lock on fd 9 (autotest.lock) repeating the try
+        echo "Waiting to aquire lock on '$PWD/autotest.lock' ..."
+        # try to get an excusive 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 "Acquired lock on '$PWD/autotest.lock'"
+        echo "Aquired lock on '$PWD/autotest.lock'"
        date
        # ----------------------
        cd ${AUTOTEST_ROOT}/autotest || \
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

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

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

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

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

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


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

-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.
+- Add hipSPARSE support for sparse mat-vec multiplications.

- 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.
+- 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 hipSPARSE support for sparse matrix-vector multiplications.
-
- More explicit and consistent formatting of the output of iterative solvers
+- More explicit and consistent formating of the output of iterative solvers
  with the new IterativeSolver::PrintLevel options. See linalg/solvers.hpp.

-Meshing improvements
--------------------
- New TMOP-based methods for hr-adaptivity, interface fitting, and tangential
-  relaxation of high-order meshes.
+- 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.

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

- Added a simpler interface to access mesh face information, see FaceInformation
-  and GetFaceInformation in the Mesh class.
+- Updated the hypre interface according to changes in hypre-2.22.1. The ADS
+  solver is now fully working on GPUs.

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

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

-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
+- Added new miniapps that use the ParELAG library, its hybrid smoothers, and the
  hierarchy of spaces created by the element-based AMG (AMGe) methodology in
  ParELAG to build multigrid solvers for H(curl) and H(div) forms. See the
  miniapps/parelag directory for more details.

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

- 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,
+- Remove the 'u' flag in the ar command, to update all files in the archive,
  avoiding file name collisions from different subdirectories.

- Various other simplifications, extensions, and bugfixes in the code.
+- 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`.


 Version 4.3, released on July 29, 2021
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -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.4.0)
+set(${PROJECT_NAME}_VERSION 4.3.1)

 # 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://docs.mfem.org/html/index.html"><img alt="Doxygen" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></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>
 </p>


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

 - The main mesh classes are:
-  + [`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)
+  + [`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)

 - The main finite element classes are:
-  + [`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)
+  + [`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)

 - The main linear algebra classes and sources are
-  + [`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)
+  + [`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)

 #### Parallel implementation

@@ -184,13 +183,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://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)
+  + [`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)

 #### GPU and general device support

@@ -199,10 +198,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://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
+  + [`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

 #### Utilities, building and documentation
 - The `general/` directory contains C++ classes that serve as utilities for
@@ -210,7 +209,8 @@ 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://docs.mfem.org.
+  that can either be built locally or browsed online at
+  https://mfem.github.io/doxygen/html/index.html.

 #### 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](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)
+- @v-dobrev (Veselin Dobrev)
+- @tzanio (Tzanio Kolev)
+- @pazner (Will Pazner)
+- @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 the `README.md` file and the `html` link in the `mfem/doxygen` repo.
+- [ ] Update and push documentation  to `mfem/doxygen`.
 - [ ] 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.24.0 release
+project webpage https://www.llnl.gov/casc/hypre. For example, the 2.20.0 release
 of hypre is available at

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

 The METIS dependency can be disabled but that is not generally recommended, see
 the option MFEM_USE_METIS.
@@ -60,6 +60,7 @@ 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
@@ -566,8 +567,7 @@ The specific libraries and their options are:
  Options: HYPRE_OPT, HYPRE_LIB.
  Versions: HYPRE >= 2.10.0b  (HYPRE built without CUDA)
            HYPRE >= 2.20.0   (HYPRE built with '--enable-mixedint')
-            HYPRE >= 2.22.1   (HYPRE built with CUDA)
-            HYPRE >= 2.23.0   (HYPRE built with HIP)
+            HYPRE >= 2.22.1   (HYPRE built with CUDA or 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 (optional), used with MFEM_USE_CODIPACK = YES
+- CoDiPack (optiobal), 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.10.
+  Versions: libCEED >= 0.8.

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

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

-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -32,7 +32,7 @@ groups_serial=(
 '"examples"
   "Examples:"
   "examples"
-   "ex{,1,2,3}[0-9].cpp"'
+   "ex{,1,2}[0-9].cpp"'
 #   "ex1.cpp"'
 '"sundials"
   "SUNDIALS examples:"
@@ -83,13 +83,11 @@ 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"
@@ -100,7 +98,7 @@ groups_parallel=(
 '"examples"
   "Examples:"
   "examples"
-   "ex{,1,2,3}[0-9]p.cpp"'
+   "ex{,1,2}[0-9]p.cpp"'
 #   "ex1p.cpp"'
 '"sundials"
   "SUNDIALS examples:"
@@ -186,7 +184,7 @@ groups_all=(
 '"examples"
   "Examples:"
   "examples"
-   "ex\"{,1,2,3}[0-9]\"{,p}.cpp"'
+   "ex\"{,1,2}[0-9]\"{,p}.cpp"'
 '"sundials"
   "SUNDIALS examples:"
   "examples/sundials"
@@ -240,8 +238,7 @@ 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"
@@ -261,13 +258,11 @@ 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-2022, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,43 +0,0 @@
-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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -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.4.0
+PROJECT_NUMBER         = v4.3.1

 # 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,11 +98,6 @@ 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
@@ -195,13 +190,12 @@ 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,7 +21,8 @@ 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://docs.mfem.org
+   https://mfem.github.io/doxygen/html/index.html
+

 Some building considerations:

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -16,7 +16,7 @@ DOXYGEN_CONF = CodeDocumentation.conf
 # doxygen uses: graphviz, latex
 html: $(DOXYGEN_CONF)
 	@# Generate the html documentation
-	@( cat $(DOXYGEN_CONF) ; echo "$(MFEM_DOXYGEN_FLAGS)" ) | doxygen -
+	@doxygen $(DOXYGEN_CONF)
 	@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,3 +29,4 @@ 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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -38,7 +38,6 @@ list(APPEND ALL_EXE_SRCS
  ex28.cpp
  ex29.cpp
  ex30.cpp
-  ex31.cpp
  )

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

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

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

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

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

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

   return 0;
 }
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,294 +0,0 @@
-//                 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());
-}
@@ -1,59 +0,0 @@
-# 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:
@@ -1,271 +0,0 @@
-//                 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());
-}
@@ -1,546 +0,0 @@
-//                   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;
-   }
-
-}
@@ -1,507 +0,0 @@
-//                   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.;
-}
@@ -1,59 +0,0 @@
-# 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:
@@ -1,649 +0,0 @@
-//                     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;
-}
@@ -1,698 +0,0 @@
-//                                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