Merge branch 'master' into main-dev

In config/githooks/pre-push, fetch 'main' instead of 'master'
Renamed the default branch from master to main
2025-10-16 07:06:45 -07:00 · 2022-12-29 17:14:55 -08:00 · 2022-12-29 12:03:23 -08:00
26 changed files with 346 additions and 1418 deletions
@@ -12,7 +12,7 @@ coverage:
        threshold: 0%
        base: auto
        branches:
-          - master
+          - main
        if_ci_failed: error
        informational: true
        only_pulls: true
@@ -22,7 +22,7 @@ coverage:
        threshold: 1%  # allows variations around the target
        base: auto
        branches:
-          - master
+          - main
        if_ci_failed: error
        only_pulls: true

@@ -40,6 +40,8 @@ Discretization improvements
  provided that neighboring hexahedra are not refined in conflicting directions.
  A new ParMesh method is added to check for such conflicts, before refinement.

+- Renamed the default GitHub branch from "master" to "main".
+
 Meshing improvements
 --------------------

@@ -652,7 +652,7 @@ if (MFEM_USE_CUDA)
 endif()

 add_subdirectory(config)
-set(MASTER_HEADERS
+set(MAIN_HEADERS
  ${PROJECT_SOURCE_DIR}/mfem.hpp
  ${PROJECT_SOURCE_DIR}/mfem-performance.hpp)

@@ -689,7 +689,7 @@ set(MFEM_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})

 # Declaring the library
-mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
+mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MAIN_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
 target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
 if (TPL_TARGETS)
@@ -888,12 +888,12 @@ install(TARGETS ${PROJECT_NAME}
  LIBRARY DESTINATION ${INSTALL_LIB_DIR}
  ARCHIVE DESTINATION ${INSTALL_LIB_DIR})

-# Install the master headers
+# Install the main headers
 foreach(Header mfem.hpp mfem-performance.hpp)
  install(FILES ${PROJECT_BINARY_DIR}/InstallHeaders/${Header}
    DESTINATION ${INSTALL_INCLUDE_DIR})
 endforeach()
-install(FILES ${MASTER_HEADERS} DESTINATION ${INSTALL_INCLUDE_DIR}/mfem)
+install(FILES ${MAIN_HEADERS} DESTINATION ${INSTALL_INCLUDE_DIR}/mfem)

 # Install the headers (except common miniapp which is installed from its subdir)
 install(DIRECTORY ${MFEM_SOURCE_DIRS}
@@ -3,10 +3,10 @@
 </p>

 <p align="center">
-<a href="https://github.com/mfem/mfem/blob/master/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
-<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amaster"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=master"></a>
-<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://github.com/mfem/mfem/blob/main/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amain"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=main"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amain"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=main"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amain"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=main"></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>
 </p>
@@ -27,7 +27,7 @@ in the MFEM community, you agree to abide by its rules.
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull
-request (PR) toward the `mfem:master` branch to propose your contribution. If
+request (PR) toward the `mfem:main` branch to propose your contribution. If
 you are planning significant code changes or have questions, you may want to
 open an [issue](https://github.com/mfem/mfem/issues) before issuing a PR. In
 addition to technical contributions, we are also interested in your results and
@@ -47,7 +47,7 @@ back to them before issuing pull requests:
  - [Pull Requests](#pull-requests)
  - [MFEM PR Rules](#mfem-pr-rules)
  - [Pull Request Checklist](#pull-request-checklist)
-  - [Master/Next Workflow](#masternext-workflow)
+  - [Main/Next Workflow](#mainnext-workflow)
  - [Releases](#releases)
  - [Release Checklist](#release-checklist)
 - [LLNL Workflow](#llnl-workflow)
@@ -66,12 +66,12 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
 ## Quick Summary

 - We encourage you to [join the MFEM organization](#mfem-organization) and create
-  development branches off `mfem:master`.
+  development branches off `mfem:main`.
 - Please follow the [developer guidelines](#developer-guidelines), in particular
  with regards to documentation and code styling.
 - Please do not commit large/binary files to the central repository (use a fork
  instead).
- Pull requests should be issued toward `mfem:master`. Make sure
+- Pull requests should be issued toward `mfem:main`. Make sure
  to check the items off the [Pull Request Checklist](#pull-request-checklist) and
  follow the [MFEM PR Rules](#mfem-pr-rules).
 - When your contribution is fully working and ready to be reviewed, add
@@ -81,8 +81,8 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
 - The reviewers have 3 weeks to evaluate the PR and work with the author to
  fix issues and implement improvements.
 - During review there should be no force pushes/rewriting history in the branch.
- After approval, MFEM developers merge the PR manually in the [mfem:next branch](#masternext-workflow).
- After a week of testing in `mfem:next`, the original PR is merged in `mfem:master`.
+- After approval, MFEM developers merge the PR manually in the [mfem:next branch](#mainnext-workflow).
+- After a week of testing in `mfem:next`, the original PR is merged in `mfem:main`.
 - We use [milestones](https://github.com/mfem/mfem/milestones) to coordinate the
  work on different PRs toward a release.
 - Don't hesitate to [contact us](#contact-information) if you have any questions.
@@ -289,7 +289,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  author, is willing to work on it and be its champion.

 - The author creates a branch for the new feature (with suffix `-dev`), off
-  the `master` branch, or another existing feature branch, for example:
+  the `main` branch, or another existing feature branch, for example:

  ```
  # Clone assuming you have setup your ssh keys on GitHub:
@@ -298,8 +298,8 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  # Alternatively, clone using the "https" protocol:
  git clone https://github.com/mfem/mfem.git

-  # Create a new feature branch starting from "master":
-  git checkout master
+  # Create a new feature branch starting from "main":
+  git checkout main
  git pull
  git checkout -b feature-dev

@@ -375,7 +375,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 ### Pull Requests

 - When your branch is ready for other developers to review / comment on
-  the code, create a pull request towards `mfem:master`.
+  the code, create a pull request towards `mfem:main`.

 - Pull request typically have titles like:

@@ -411,7 +411,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  assigned, the PR is considered under review. To help with the review process
  there should be no force pushes/rewriting history in the branch.

- After approval, the PR is [tested](#masternext-workflow) for a week with
+- After approval, the PR is [tested](#mainnext-workflow) for a week with
  other approved PRs in the `mfem:next` branch.

 - Consider manually running the tests in `tests/scripts` before merging in
@@ -447,7 +447,7 @@ The Pull Request (PR) approval process in MFEM is similar to the approval of pap

 3. A PR can be (manually) merged in the *next* branch only if 2 of the assigned reviewers have approved it and it has passed internal testing. This merge can be performed by any of the assigned reviewers or by any of the editors.

-4. A PR can be merged in the *master* branch only if it has been tested successfully for a week in *next* and an editor has (optionally) taken a final look. This merge can be performed only by one of the editors.
+4. A PR can be merged in the *main* branch only if it has been tested successfully for a week in *next* and an editor has (optionally) taken a final look. This merge can be performed only by one of the editors.

 #### Responsibilities of Editors

@@ -470,7 +470,7 @@ The current list of MFEM editors is:

 5. To remind the reviewers about timely completion of their review.

-6. To take a final look and complete the PR merge in *master*. The final look step is optional and shouldn't take more than 3 days.
+6. To take a final look and complete the PR merge in *main*. The final look step is optional and shouldn't take more than 3 days.

 7. The assignment of bugfixes should be expedited proportional to their importance, e.g. in some cases the editor can assign much shorter review window.

@@ -492,7 +492,7 @@ Everyone on the MFEM team can be asked to serve as a reviewer on a PR in their a

 5. To test the PR thoroughly before merging in *next*. The PR author is also encouraged to perform testing and inform the reviewers about the results.

-6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *master*.
+6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *main*.

 7. The review of bugfixes should be expedited proportional to their importance. The review window can be much less than three weeks in such cases.

@@ -582,11 +582,11 @@ Before a PR can be merged, it should satisfy the following:
   - [ ] Update internal tests to include the new features.


-### Master/Next Workflow
+### Main/Next Workflow

-MFEM uses a `master`/`next`-branch workflow as described below:
+MFEM uses a `main`/`next`-branch workflow as described below:

- The `master` branch should always be of release quality and changes should not
+- The `main` branch should always be of release quality and changes should not
  be merged until they have been fully tested. This branch is protected, and
  changes can only be made through pull requests.

@@ -613,20 +613,20 @@ MFEM uses a `master`/`next`-branch workflow as described below:

 - After a week of testing in `next` (excluding bugfixes), both on GitHub, as
  well as [internally](#tests-at-llnl) at LLNL, the original PR is merged into
-  `master` (provided there are no issues).
+  `main` (provided there are no issues).

 - After the merge, the feature branch is deleted (unless it is a long-term
  project with periodic PRs).

 - The `next` branch is used just for integrated testing of all PRs approved for
-  merging into `master` to verify that each works individually and that all of
+  merging into `main` to verify that each works individually and that all of
  them work as a group. This branch can be discarded at any time, though we
  typically do that only at the end of a [release cycle](#releases).


 ### Releases

- Releases are just tags in the `master` branch, e.g. https://github.com/mfem/mfem/releases/tag/v3.3.2,
+- Releases are just tags in the `main` branch, e.g. https://github.com/mfem/mfem/releases/tag/v3.3.2,
  and have a version that ends in an even "patch" number, e.g. `v3.2.2` or
  `v3.4` (by convention `v3.4` is the same as `v3.4.0`.)  Between releases, the
  version ends in an odd "patch" number, e.g. `v3.3.3`.
@@ -679,19 +679,19 @@ MFEM uses a `master`/`next`-branch workflow as described below:

 ### Mirroring on Bitbucket

- The GitHub `master` and `next` branches are mirrored to the LLNL institutional
-  Bitbucket repository as `gh-master` and `gh-next`.
+- The GitHub `main` and `next` branches are mirrored to the LLNL institutional
+  Bitbucket repository as `gh-main` and `gh-next`.

- `gh-master` is merged into LLNL's internal `master` through pull requests; write
-  permissions to `master` are restricted to ensure this is the only way in which it
+- `gh-main` is merged into LLNL's internal `main` through pull requests; write
+  permissions to `main` are restricted to ensure this is the only way in which it
  gets updated.

 - We never push directly from LLNL to GitHub.

 - Versions of the code on LLNL's internal server, from most to least stable:
  - MFEM official release on mfem.org -- Most stable, tested in many apps.
-  - `mfem:master` -- Recent development version, guaranteed to work.
-  - `mfem:gh-master` -- Stable development version, passed testing, you can use
+  - `mfem:main` -- Recent development version, guaranteed to work.
+  - `mfem:gh-main` -- Stable development version, passed testing, you can use
     it to build your code between releases.
  - `mfem:gh-next` -- Bleeding-edge development version, may be broken, use at
     your own risk.
@@ -721,14 +721,14 @@ directory.


 ### Linux and Mac smoke tests
-We use GitHub Actions to drive the default tests on the `master` and `next`
+We use GitHub Actions to drive the default tests on the `main` and `next`
 branches. See the `.github/workflows` files and the logs at
 [https://github.com/mfem/mfem/actions](https://github.com/mfem/mfem/actions).

 Testing using GitHub Actions should be kept lightweight, as there is a time
 constraint on jobs. Two virtual machines are configured - Mac (OS X) and Linux.

- Tests on the `master` branch are triggered whenever a PR is issued on this branch.
+- Tests on the `main` branch are triggered whenever a PR is issued on this branch.
 - Tests on the `next` branch are currently scheduled to run each night.


@@ -744,7 +744,7 @@ and debug build is performed with a simple run of `ex1` to verify the executable

 ### Tests at LLNL

- We mirror the `master` and `next` branches internally (to `gh-master` and
+- We mirror the `main` and `next` branches internally (to `gh-main` and
  `gh-next`) and run longer nightly tests via cron. On the weekends, a more
  extensive test is run which extracts and executes all the different sample
  runs from each example and most miniapps.
@@ -33,8 +33,8 @@ RUN mkdir -p /opt/mfem-env \
 RUN cd /opt/mfem-env && \
    . /opt/spack/share/spack/setup-env.sh && \
    spack env activate . && \
-    spack develop --path /code mfem@master+examples+miniapps && \
-    spack add mfem@master+examples+miniapps && \
+    spack develop --path /code mfem@main+examples+miniapps && \
+    spack add mfem@main+examples+miniapps # && \
    spack install

 # ensure mfem always on various paths
@@ -1,11 +1,11 @@
 spack:
-  specs: [mfem@master+examples+miniapps]
-  view: 
+  specs: [mfem@main+examples+miniapps]
+  view:
    mfem:
-      root: /opt/mfem-view 
+      root: /opt/mfem-view
      link_type: copy
  concretization: together
  develop:
    mfem:
      path: /code
-      spec: mfem@master+examples+miniapps
+      spec: mfem@main+examples+miniapps
@@ -115,7 +115,7 @@ fi

 # branch-history
 if [[ "${option}" == "--history" || "${option}" == "" ]]; then
-  git fetch origin master:master
+  git fetch origin main:main
  cd tests/scripts
  if ! ./runtest branch-history; then code=1; fi
  cd -
@@ -1,6 +1,6 @@
 # Jupyter Notebooks using xeus-cling

-[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mfem/mfem/master?filepath=examples%2Fjupyter%2Fex.ipynb)
+[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mfem/mfem/main?filepath=examples%2Fjupyter%2Fex.ipynb)

 [xeus-cling](https://github.com/jupyter-xeus/xeus-cling) is a C++ Jupyter Kernel based on [cling](https://github.com/root-project/cling),
 which can be used to create interactive C++ MFEM and GLVis notebooks.
@@ -28,8 +28,10 @@
 //
 //               The example demonstrates the use of nonlinear operators (the
 //               class ConductionOperator defining C(u)), as well as their
-//               implicit time integration. By default, this example uses the
-//               SUNDIALS ODE solvers from CVODE and ARKODE.
+//               implicit time integration. Note that implementing the method
+//               ConductionOperator::ImplicitSolve is the only requirement for
+//               high-order implicit (SDIRK) time integration. By default, this
+//               example uses the SUNDIALS ODE solvers from CVODE and ARKODE.
 //
 //               We recommend viewing examples 2, 9 and 10 before viewing this
 //               example.
@@ -49,16 +51,15 @@ using namespace mfem;
 *  and K(u) is the diffusion operator with diffusivity depending on u:
 *  (\kappa + \alpha u).
 *
- *  Class ConductionOperator represents the above ODE operator as a
- *  TimeDependentOperator for use with native MFEM integrators and CVODE
- *  integrators, i.e., F(u, k, t) = G(u, t) with F(u, du/dt, t) = du/dt and
- *  G(u, t) = -K(u) u
+ *  Class ConductionOperatorOperator represents the above ODE operator in the
+ *  general form F(u, k, t) = G(u, t) where
 *
- *  Class ConductionOperator represents the above ODE operator as an
- *  ARKStepODE for use with ARKODE integrators, i.e., either M du/dt = -K(u) u
- *  (mass form) or du/dt = -inv(M) K(u) u (MFEM form)
+ *    1. F(u, du/dt, t) = du/dt            (ODE is expressed in EXPLICIT form)
+ *       G(u, t) = - inv(M) K(u) u
+ *    2. F(u, du/dt, t) = M du/dt          (ODE is expressed in IMPLICIT form)
+ *       G(u, t) = - K(u) u
 */
-class ConductionOperator : public TimeDependentOperator, public ARKStepODE
+class ConductionOperator : public TimeDependentOperator
 {
   FiniteElementSpace &fespace;
   Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
@@ -80,90 +81,50 @@ class ConductionOperator : public TimeDependentOperator, public ARKStepODE

   mutable Vector z; // auxiliary vector

-   const bool use_mass_form;
-
 public:

   ConductionOperator(FiniteElementSpace &f, const real_t alpha,
                      const real_t kappa, const Vector &u,
-                      const bool use_mass_form);
+                      const Type &ode_expression_type);

   // Compute K(u_n) for use as an approximation in - K(u) u
   void SetConductionTensor(const Vector &u);

-   // ********* methods for MFEM native time integrators *********
+   /** Compute G(u, t) as defined in the IMPLICIT expression form of the ODE
+    operator, i.e., @a v = - K(u_n) @a u. Note that K(u_n) is an
+      approximation to K(u). */
+   void ExplicitMult(const Vector &u, Vector &v) const override;

-   /** Solve for k in F(u, k, t) = G(u, t), i.e., @a k = - inv(M) K(u_n) @a u.
+   /** Solve for k in F(u, k, t) = G(u, t) for either EXPLICIT or IMPLICIT
+       expression forms of the ODE operator, i.e., @a k = - inv(M) K(u_n) @a u.
       Note that K(u_n) is an approximation to K(u). */
   void Mult(const Vector &u, Vector &k) const override;

-   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t), i.e.,
-       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u .
-       Note that K(u_n) is an approximation to K(u). */
+   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t) for either EXPLICIT
+       or IMPLICIT expression forms of the ODE operator, i.e.,
+       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u . Note that K(u_n) is an
+       approximation to K(u). */
   void ImplicitSolve(const real_t gam, const Vector &u, Vector &k) override;

-      // ********* methods for ARKODE time integrators *********
-
-   // TODO: add comments
-   int ARKSize() const override;
-
-   // TODO: add comments
-   bool ARKInMassForm() const override;
-
-   // TODO: add comments
-   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const override;
-
-   // TODO: add comments
-   int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &fu,
-                        int jok, int *jcur, real_t gam) override;
-
-   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
-       with the residual @a r providing either
-        1. @a r = G - F = inv(M) f(u) - k       (MFEM form)
-        1. @a r = G - F = f(u) - M k            (mass form)
-       */
-   int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
-
-   int ARKMassSetup(const real_t t) override;
-
-   int ARKMassSolve(const Vector &b, Vector &x, real_t tol) override;
-
-   int ARKMassMult(const Vector &x, Vector &v) override;
-
-   // ********* methods for CVODE time integrators *********
-   // note these methods merely call the corresponding ARKStepODE methods until
-   // the CVODESolver is refactored to use specialized interface like ARKStepODE
-
-   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F, i.e.,
+   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F for
+       either EXPLICIT or IMPLICIT expression forms of the ODE operator, i.e.,
       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
-                        real_t gam) override
-   {
-      return ARKImplicitSetup(u, 0.0, fu, jok, jcur, gam); // the ODE is autonomous
-   }
+                        real_t gam) override;

   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
-       with the residual @a r providing @a r = G - F = inv(M) f(u) - k. */
-   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override
-   {
-      return ARKImplicitSolve(r, dk, tol);
-   }
+       with the residual @a r providing either
+        1. @a r = G - F = inv(M) f(u) - k       (EXPLICIT expression form)
+        1. @a r = G - F = f(u) - M k            (IMPLICIT expression form)
+       */
+   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;

-   int SUNMassSetup() override
-   {
-      return ARKMassSetup(0.0); // the ODE is autonomous
-   }
+   int SUNMassSetup() override;

-   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override
-   {
-      return ARKMassSolve(b, x, tol);
-   }
+   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override;

-   int SUNMassMult(const Vector &x, Vector &v) override
-   {
-      return ARKMassMult(x, v);
-   }
+   int SUNMassMult(const Vector &x, Vector &v) override;
 };

 real_t InitialTemperature(const Vector &x)
@@ -284,7 +245,16 @@ int main(int argc, char *argv[])
   u_gf.GetTrueDofs(u);

   // 6. Initialize the conduction ODE operator and the visualization.
-   ConductionOperator oper(fespace, alpha, kappa, u, use_mass_solver);
+   ConductionOperator::Type ode_expression_type;
+   if (use_mass_solver)
+   {
+      ode_expression_type = ConductionOperator::Type::IMPLICIT;
+   }
+   else
+   {
+      ode_expression_type = ConductionOperator::Type::EXPLICIT;
+   }
+   ConductionOperator oper(fespace, alpha, kappa, u, ode_expression_type);

   u_gf.SetFromTrueDofs(u);
   {
@@ -382,7 +352,7 @@ int main(int argc, char *argv[])
         }
         std::unique_ptr<ARKStepSolver> arkode(
            new ARKStepSolver(arkode_solver_type));
-         arkode->Init(&oper);
+         arkode->Init(oper);
         arkode->SetSStolerances(reltol, abstol);
         arkode->SetMaxStep(dt);
         if (ode_solver_type == 11 || ode_solver_type == 14)
@@ -475,10 +445,9 @@ int main(int argc, char *argv[])
 ConductionOperator::ConductionOperator(FiniteElementSpace &fes,
                                       const real_t alpha, const real_t kappa,
                                       const Vector &u,
-                                       const bool use_mass_form)
-   : TimeDependentOperator(fes.GetTrueVSize(), 0.0),
-     fespace(fes), M(&fespace), alpha(alpha), kappa(kappa), z(height),
-     use_mass_form(use_mass_form)
+                                       const Type &ode_expression_type)
+   : TimeDependentOperator(fes.GetTrueVSize(), 0.0, ode_expression_type),
+     fespace(fes), M(&fespace), alpha(alpha), kappa(kappa), z(height)
 {
   // specify a relative tolerance for all solves with MFEM integrators
   const real_t rel_tol = 1e-8;
@@ -505,16 +474,6 @@ ConductionOperator::ConductionOperator(FiniteElementSpace &fes,
   SetConductionTensor(u);
 }

-int ConductionOperator::ARKSize() const
-{
-   return z.Size();
-}
-
-bool ConductionOperator::ARKInMassForm() const
-{
-   return use_mass_form;
-}
-
 void ConductionOperator::SetConductionTensor(const Vector &u)
 {
   // Compute K(u_n).
@@ -532,27 +491,17 @@ void ConductionOperator::SetConductionTensor(const Vector &u)
   K->FormSystemMatrix(ess_tdof_list, Kmat);
 }

-void ConductionOperator::ARKEvaluateRHS(const Vector &u, const real_t t,
-                                        Vector &result) const
+void ConductionOperator::ExplicitMult(const Vector &u, Vector &v) const
 {
-   if (use_mass_form) // compute -K(u_n) u.
-   {
-      Kmat.Mult(u, result);
-      result.Neg();
-   }
-   else // compute -inv(M) K(u_n) u
-   {
-      Kmat.Mult(u, z);
-      z.Neg();
-      M_solver.Mult(z, result);
-   }
+   // Compute - K(u_n) u.
+   Kmat.Mult(u, v);
+   v.Neg();
 }

 void ConductionOperator::Mult(const Vector &u, Vector &k) const
 {
   // Compute - inv(M) K(u_n) u.
-   Kmat.Mult(u, z);
-   z.Neg();
+   ExplicitMult(u, z);
   M_solver.Mult(z, k);
 }

@@ -560,16 +509,14 @@ void ConductionOperator::ImplicitSolve(const real_t gam, const Vector &u,
                                       Vector &k)
 {
   // Solve for k in M k = - K(u_n) [u + gam*k].
-   Kmat.Mult(u, z);
-   z.Neg();
+   ExplicitMult(u, z);
   T = std::unique_ptr<SparseMatrix>(Add(1.0, Mmat, gam, Kmat));
   T_solver.SetOperator(*T);
   T_solver.Mult(z, k);
 }

-int ConductionOperator::ARKImplicitSetup(const Vector &u, const real_t t,
-                                         const Vector &fu, int jok, int *jcur,
-                                         real_t gam)
+int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
+                                         int jok, int *jcur, real_t gam)
 {
   // Compute T = M + gamma K(u_n).
   T = std::unique_ptr<SparseMatrix>(Add(1.0, Mmat, gam, Kmat));
@@ -578,23 +525,23 @@ int ConductionOperator::ARKImplicitSetup(const Vector &u, const real_t t,
   return SUN_SUCCESS;
 }

-int ConductionOperator::ARKImplicitSolve(const Vector &r, Vector &dk,
+int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
                                         real_t tol)
 {
   // Solve the system [M + gamma K(u_n)] dk = - K(u_n) u - M k.
   // What value r is providing depends on the ODE expression form:
-   //   MFEM form: r = -inv(M) K(u_n) u - k
-   //   mass form: r = -K(u_n) u - M k
+   //   EXPLICIT form: r = -inv(M) K(u_n) u - k
+   //   IMPLICIT form: r = -K(u_n) u - M k
   T_solver.SetRelTol(tol);
-   if (use_mass_form)
-   {
-      T_solver.Mult(r, dk);
-   }
-   else
+   if (isExplicit())
   {
      Mmat.Mult(r, z);
      T_solver.Mult(z, dk);
   }
+   else
+   {
+      T_solver.Mult(r, dk);
+   }
   if (T_solver.GetConverged())
   {
      return SUN_SUCCESS;
@@ -605,13 +552,13 @@ int ConductionOperator::ARKImplicitSolve(const Vector &r, Vector &dk,
   }
 }

-int ConductionOperator::ARKMassSetup(const real_t t)
+int ConductionOperator::SUNMassSetup()
 {
   // Do nothing b/c mass solver was setup in constructor.
   return SUN_SUCCESS;
 }

-int ConductionOperator::ARKMassSolve(const Vector &b, Vector &x, real_t tol)
+int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
 {
   // Solve the system M x = b.
   M_solver.SetRelTol(tol);
@@ -626,7 +573,7 @@ int ConductionOperator::ARKMassSolve(const Vector &b, Vector &x, real_t tol)
   }
 }

-int ConductionOperator::ARKMassMult(const Vector &x, Vector &v)
+int ConductionOperator::SUNMassMult(const Vector &x, Vector &v)
 {
   // Compute M x.
   Mmat.Mult(x, v);
@@ -29,8 +29,10 @@
 //
 //               The example demonstrates the use of nonlinear operators (the
 //               class ConductionOperator defining C(u)), as well as their
-//               implicit time integration. By default, this example uses the
-//               SUNDIALS ODE solvers from CVODE and ARKODE.
+//               implicit time integration. Note that implementing the method
+//               ConductionOperator::ImplicitSolve is the only requirement for
+//               high-order implicit (SDIRK) time integration. By default, this
+//               example uses the SUNDIALS ODE solvers from CVODE and ARKODE.
 //
 //               We recommend viewing examples 2, 9 and 10 before viewing this
 //               example.
@@ -50,16 +52,15 @@ using namespace mfem;
 *  and K(u) is the diffusion operator with diffusivity depending on u:
 *  (\kappa + \alpha u).
 *
- *  Class ConductionOperator represents the above ODE operator as a
- *  TimeDependentOperator for use with native MFEM integrators and CVODE
- *  integrators, i.e., F(u, k, t) = G(u, t) with F(u, du/dt, t) = du/dt and
- *  G(u, t) = -K(u) u
+ *  Class ConductionOperatorOperator represents the above ODE operator in the
+ *  general form F(u, k, t) = G(u, t) where either
 *
- *  Class ConductionOperator represents the above ODE operator as an
- *  ARKStepODE for use with ARKODE integrators, i.e., either M du/dt = -K(u) u
- *  (mass form) or du/dt = -inv(M) K(u) u (MFEM form)
+ *    1. F(u, du/dt, t) = du/dt            (ODE is expressed in EXPLICIT form)
+ *       G(u, t) = - inv(M) K(u) u
+ *    2. F(u, du/dt, t) = M du/dt          (ODE is expressed in IMPLICIT form)
+ *       G(u, t) = - K(u) u
 */
-class ConductionOperator : public TimeDependentOperator, public ARKStepODE
+class ConductionOperator : public TimeDependentOperator
 {
   ParFiniteElementSpace &fespace;
   Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
@@ -81,90 +82,50 @@ class ConductionOperator : public TimeDependentOperator, public ARKStepODE

   mutable Vector z; // auxiliary vector

-   const bool use_mass_form;
-
 public:

   ConductionOperator(ParFiniteElementSpace &f, const real_t alpha,
                      const real_t kappa, const Vector &u,
-                      const bool use_mass_form);
+                      const Type &ode_expression_type);

   // Compute K(u_n) for use as an approximation in - K(u) u
   void SetConductionTensor(const Vector &u);

-   // ********* methods for MFEM native time integrators *********
+   /** Compute G(u, t) as defined in the IMPLICIT expression form of the ODE
+       operator, i.e., @a v = - K(u_n) @a u. Note that K(u_n) is an
+       approximation to K(u). */
+   void ExplicitMult(const Vector &u, Vector &v) const override;

-   /** Solve for k in F(u, k, t) = G(u, t), i.e., @a k = - inv(M) K(u_n) @a u.
+   /** Solve for k in F(u, k, t) = G(u, t) for either EXPLICIT or IMPLICIT
+       expression forms of the ODE operator, i.e., @a k = - inv(M) K(u_n) @a u.
       Note that K(u_n) is an approximation to K(u). */
   void Mult(const Vector &u, Vector &k) const override;

-   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t), i.e.,
-       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u .
-       Note that K(u_n) is an approximation to K(u). */
+   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t) for either EXPLICIT
+       or IMPLICIT expression forms of the ODE operator, i.e.,
+       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u . Note that K(u_n) is an
+       approximation to K(u). */
   void ImplicitSolve(const real_t gam, const Vector &u, Vector &k) override;

-   // ********* methods for ARKODE time integrators *********
-
-   // TODO: add comments
-   int ARKSize() const override;
-
-   // TODO: add comments
-   bool ARKInMassForm() const override;
-
-   // TODO: add comments
-   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const override;
-
-   // TODO: add comments
-   int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &fu,
-                        int jok, int *jcur, real_t gam) override;
-
-   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
-       with the residual @a r providing either
-        1. @a r = G - F = inv(M) f(u) - k       (MFEM form)
-        1. @a r = G - F = f(u) - M k            (mass form)
-       */
-   int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
-
-   int ARKMassSetup(const real_t t) override;
-
-   int ARKMassSolve(const Vector &b, Vector &x, real_t tol) override;
-
-   int ARKMassMult(const Vector &x, Vector &v) override;
-
-   // ********* methods for CVODE time integrators *********
-   // note these methods merely call the corresponding ARKStepODE methods until
-   // the CVODESolver is refactored to use specialized interface like ARKStepODE
-
-   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F, i.e.,
+   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F for
+       either EXPLICIT or IMPLICIT expression forms of the ODE operator, i.e.,
       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
-                        real_t gam) override
-   {
-      return ARKImplicitSetup(u, 0.0, fu, jok, jcur, gam); // the ODE is autonomous
-   }
+                        real_t gam) override;

   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
-       with the residual @a r providing @a r = G - F = inv(M) f(u) - k. */
-   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override
-   {
-      return ARKImplicitSolve(r, dk, tol);
-   }
+       with the residual @a r providing either
+        1. @a r = G - F = inv(M) f(u) - k       (EXPLICIT expression form)
+        1. @a r = G - F = f(u) - M k            (IMPLICIT expression form)
+       */
+   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;

-   int SUNMassSetup() override
-   {
-      return ARKMassSetup(0.0); // the ODE is autonomous
-   }
+   int SUNMassSetup() override;

-   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override
-   {
-      return ARKMassSolve(b, x, tol);
-   }
+   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override;

-   int SUNMassMult(const Vector &x, Vector &v) override
-   {
-      return ARKMassMult(x, v);
-   }
+   int SUNMassMult(const Vector &x, Vector &v) override;
 };

 real_t InitialTemperature(const Vector &x)
@@ -312,7 +273,16 @@ int main(int argc, char *argv[])
   u_gf.GetTrueDofs(u);

   // 8. Initialize the conduction ODE operator and the visualization.
-   ConductionOperator oper(fespace, alpha, kappa, u, use_mass_solver);
+   ConductionOperator::Type ode_expression_type;
+   if (use_mass_solver)
+   {
+      ode_expression_type = ConductionOperator::Type::IMPLICIT;
+   }
+   else
+   {
+      ode_expression_type = ConductionOperator::Type::EXPLICIT;
+   }
+   ConductionOperator oper(fespace, alpha, kappa, u, ode_expression_type);

   u_gf.SetFromTrueDofs(u);
   {
@@ -424,7 +394,7 @@ int main(int argc, char *argv[])
         }
         std::unique_ptr<ARKStepSolver> arkode(
            new ARKStepSolver(MPI_COMM_WORLD, arkode_solver_type));
-         arkode->Init(&oper);
+         arkode->Init(oper);
         arkode->SetSStolerances(reltol, abstol);
         arkode->SetMaxStep(dt);
         if (ode_solver_type == 11 || ode_solver_type == 14)
@@ -527,11 +497,10 @@ int main(int argc, char *argv[])
 ConductionOperator::ConductionOperator(ParFiniteElementSpace &fes,
                                       const real_t alpha, const real_t kappa,
                                       const Vector &u,
-                                       const bool use_mass_form)
-   : TimeDependentOperator(fes.GetTrueVSize(), 0.0),
+                                       const Type &ode_expression_type)
+   : TimeDependentOperator(fes.GetTrueVSize(), 0.0, ode_expression_type),
     fespace(fes), M(&fespace), alpha(alpha), kappa(kappa),
-     M_solver(fes.GetComm()), T_solver(fes.GetComm()), z(height),
-     use_mass_form(use_mass_form)
+     M_solver(fes.GetComm()), T_solver(fes.GetComm()), z(height)
 {
   // specify a relative tolerance for all solves with MFEM integrators
   const real_t rel_tol = 1e-8;
@@ -559,16 +528,6 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &fes,
   SetConductionTensor(u);
 }

-int ConductionOperator::ARKSize() const
-{
-   return z.Size();
-}
-
-bool ConductionOperator::ARKInMassForm() const
-{
-   return use_mass_form;
-}
-
 void ConductionOperator::SetConductionTensor(const Vector &u)
 {
   // Compute K(u_n).
@@ -586,27 +545,17 @@ void ConductionOperator::SetConductionTensor(const Vector &u)
   K->FormSystemMatrix(ess_tdof_list, Kmat);
 }

-void ConductionOperator::ARKEvaluateRHS(const Vector &u, const real_t t,
-                                        Vector &result) const
+void ConductionOperator::ExplicitMult(const Vector &u, Vector &v) const
 {
-   if (use_mass_form) // compute -K(u_n) u.
-   {
-      Kmat.Mult(u, result);
-      result.Neg();
-   }
-   else // compute -inv(M) K(u_n) u
-   {
-      Kmat.Mult(u, z);
-      z.Neg();
-      M_solver.Mult(z, result);
-   }
+   // Compute - K(u_n) u.
+   Kmat.Mult(u, v);
+   v.Neg();
 }

 void ConductionOperator::Mult(const Vector &u, Vector &k) const
 {
   // Compute - inv(M) K(u_n) u.
-   Kmat.Mult(u, z);
-   z.Neg();
+   ExplicitMult(u, z);
   M_solver.Mult(z, k);
 }

@@ -614,16 +563,14 @@ void ConductionOperator::ImplicitSolve(const real_t gam, const Vector &u,
                                       Vector &k)
 {
   // Solve for k in M k = - K(u_n) [u + gam*k].
-   Kmat.Mult(u, z);
-   z.Neg();
+   ExplicitMult(u, z);
   T = std::unique_ptr<HypreParMatrix>(Add(1.0, Mmat, gam, Kmat));
   T_solver.SetOperator(*T);
   T_solver.Mult(z, k);
 }

-int ConductionOperator::ARKImplicitSetup(const Vector &u, const real_t t,
-                                         const Vector &fu, int jok, int *jcur,
-                                         real_t gam)
+int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
+                                         int jok, int *jcur, real_t gam)
 {
   // Compute T = M + gamma K(u_n).
   T = std::unique_ptr<HypreParMatrix>(Add(1.0, Mmat, gam, Kmat));
@@ -632,23 +579,23 @@ int ConductionOperator::ARKImplicitSetup(const Vector &u, const real_t t,
   return SUN_SUCCESS;
 }

-int ConductionOperator::ARKImplicitSolve(const Vector &r, Vector &dk,
+int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
                                         real_t tol)
 {
   // Solve the system [M + gamma K(u_n)] dk = - K(u_n) u - M k.
   // What value r is providing depends on the ODE expression form:
-   //   MFEM form: r = -inv(M) K(u_n) u - k
-   //   mass form: r = -K(u_n) u - M k
+   //   EXPLICIT form: r = -inv(M) K(u_n) u - k
+   //   IMPLICIT form: r = -K(u_n) u - M k
   T_solver.SetRelTol(tol);
-   if (use_mass_form)
-   {
-      T_solver.Mult(r, dk);
-   }
-   else
+   if (isExplicit())
   {
      Mmat.Mult(r, z);
      T_solver.Mult(z, dk);
   }
+   else
+   {
+      T_solver.Mult(r, dk);
+   }
   if (T_solver.GetConverged())
   {
      return SUN_SUCCESS;
@@ -659,13 +606,13 @@ int ConductionOperator::ARKImplicitSolve(const Vector &r, Vector &dk,
   }
 }

-int ConductionOperator::ARKMassSetup(const real_t t)
+int ConductionOperator::SUNMassSetup()
 {
   // Do nothing b/c mass solver was setup in constructor.
   return SUN_SUCCESS;
 }

-int ConductionOperator::ARKMassSolve(const Vector &b, Vector &x, real_t tol)
+int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
 {
   // Solve the system M x = b.
   M_solver.SetRelTol(tol);
@@ -680,7 +627,7 @@ int ConductionOperator::ARKMassSolve(const Vector &b, Vector &x, real_t tol)
   }
 }

-int ConductionOperator::ARKMassMult(const Vector &x, Vector &v)
+int ConductionOperator::SUNMassMult(const Vector &x, Vector &v)
 {
   // Compute M x.
   Mmat.Mult(x, v);
@@ -119,7 +119,7 @@ public:
    and advection matrices, and b describes the flow on the boundary. This can
    be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
    used to evaluate the right-hand side. */
-class FE_Evolution : public TimeDependentOperator, public ARKStepODE
+class FE_Evolution : public TimeDependentOperator
 {
 private:
   BilinearForm &M, &K;
@@ -133,14 +133,9 @@ private:
 public:
   FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);

-   // TimeDependentOperator methods for MFEM native and CVODE time integrators
   virtual void Mult(const Vector &x, Vector &y) const;
   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);

-   // ARKStepODE methods for ARKODE time integrators
-   int ARKSize() const override;
-   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector& result) const override;
-
   virtual ~FE_Evolution();
 };

@@ -409,14 +404,14 @@ int main(int argc, char *argv[])
         ode_solver = cvode; break;
      case 8:
         arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
-         arkode->Init(&adv);
+         arkode->Init(adv);
         arkode->SetSStolerances(reltol, abstol);
         arkode->SetMaxStep(dt);
         arkode->SetOrder(4);
         ode_solver = arkode; break;
      case 9:
         arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
-         arkode->Init(&adv);
+         arkode->Init(adv);
         arkode->SetSStolerances(reltol, abstol);
         arkode->SetMaxStep(dt);
         arkode->SetERKTableNum(ARKODE_FEHLBERG_13_7_8);
@@ -525,19 +520,6 @@ void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
   dg_solver->Mult(z, k);
 }

-int FE_Evolution::ARKSize() const
-{
-   return z.Size();
-}
-
-void FE_Evolution::ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const
-{
-   // y = M^{-1} (K x + b)
-   K.Mult(u, z);
-   z += b;
-   M_solver.Mult(z, result);
-}
-
 FE_Evolution::~FE_Evolution()
 {
   delete M_prec;
@@ -206,7 +206,7 @@ public:
    and advection matrices, and b describes the flow on the boundary. This can
    be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
    used to evaluate the right-hand side. */
-class FE_Evolution : public TimeDependentOperator, public ARKStepODE
+class FE_Evolution : public TimeDependentOperator
 {
 private:
   OperatorHandle M, K;
@@ -221,14 +221,9 @@ public:
   FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_, const Vector &b_,
                PrecType prec_type);

-   // TimeDependentOperator methods for MFEM native and CVODE time integrators
   virtual void Mult(const Vector &x, Vector &y) const;
   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);

-   // ARKStepODE methods for ARKODE time integrators
-   int ARKSize() const override;
-   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector& result) const override;
-
   virtual ~FE_Evolution();
 };

@@ -580,7 +575,7 @@ int main(int argc, char *argv[])
      case 8:
      case 9:
         arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::EXPLICIT);
-         arkode->Init(&adv);
+         arkode->Init(adv);
         arkode->SetSStolerances(reltol, abstol);
         arkode->SetMaxStep(dt);
         if (ode_solver_type == 9)
@@ -748,19 +743,6 @@ void FE_Evolution::Mult(const Vector &x, Vector &y) const
   M_solver.Mult(z, y);
 }

-int FE_Evolution::ARKSize() const
-{
-   return z.Size();
-}
-
-void FE_Evolution::ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const
-{
-   // y = M^{-1} (K x + b)
-   K->Mult(u, z);
-   z += b;
-   M_solver.Mult(z, result);
-}
-
 FE_Evolution::~FE_Evolution()
 {
   delete M_prec;
@@ -436,7 +436,7 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
   // In parallel, the result is in 'py' which is an alias for 'aux2'.
 }

-Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
+Operator &NonlinearForm::GetGradient(const Vector &x) const
 {
   if (ext)
   {
@@ -644,8 +644,6 @@ Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
      }
   }

-   if (!finalize) { return *Grad; }
-
   if (!Grad->Finalized())
   {
      Grad->Finalize(skip_zeros);
@@ -1205,14 +1203,7 @@ const BlockVector &BlockNonlinearForm::Prolongate(const BlockVector &bx) const
      aux1.Update(block_offsets);
      for (int s = 0; s < fes.Size(); s++)
      {
-         if (P[s])
-         {
-            P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
-         }
-         else
-         {
-            aux1.GetBlock(s) = bx.GetBlock(s);
-         }
+         P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
      }
      return aux1;
   }
@@ -1241,16 +1232,11 @@ void BlockNonlinearForm::Mult(const Vector &x, Vector &y) const
      {
         cP[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
      }
-      else if (needs_prolongation)
-      {
-         by.GetBlock(s) = pby.GetBlock(s);
-      }
      by.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
   }
 }

-void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
-                                                bool finalize) const
+void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
 {
   const int skip_zeros = 0;
   Array<Array<int> *> vdofs(fes.Size());
@@ -1504,7 +1490,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
      }
   }

-   if (finalize && !Grads(0,0)->Finalized())
+   if (!Grads(0,0)->Finalized())
   {
      for (int i=0; i<fes.Size(); ++i)
      {
@@ -1543,23 +1529,7 @@ Operator &BlockNonlinearForm::GetGradient(const Vector &x) const
         for (int s2 = 0; s2 < fes.Size(); ++s2)
         {
            delete cGrads(s1, s2);
-            if (cP[s1] && cP[s2])
-            {
-               cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
-            }
-            else if (cP[s1])
-            {
-               cGrads(s1, s2) = TransposeMult(*cP[s1], *Grads(s1, s2));
-            }
-            else if (cP[s2])
-            {
-               cGrads(s1, s2) = mfem::Mult(*Grads(s1, s2), *cP[s2]);
-            }
-            else
-            {
-               cGrads(s1, s2) = NULL;
-               continue;
-            }
+            cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
            mGrads(s1, s2) = cGrads(s1, s2);
         }
      }
@@ -217,12 +217,7 @@ public:
       In general, @a x may have non-homogeneous essential boundary values.

       The state @a x must be a true-dof vector. */
-   Operator &GetGradient(const Vector &x) const override { return GetGradient(x, true); }
-
-   /** @brief Compute the gradient Operator of the NonlinearForm corresponding
-       to the state @a x with optional finalization and elimintaion. */
-   /** @see GetGradient(const Vector &) */
-   Operator &GetGradient(const Vector &x, bool finalize) const;
+   Operator &GetGradient(const Vector &x) const override;

   /// Update the NonlinearForm to propagate updates of the associated FE space.
   /** After calling this method, the essential boundary conditions need to be
@@ -313,7 +308,7 @@ protected:
   void MultBlocked(const BlockVector &bx, BlockVector &by) const;

   /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx, bool finalize = true) const;
+   void ComputeGradientBlocked(const BlockVector &bx) const;

 public:
   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
@@ -151,15 +151,6 @@ void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
   }
 }

-HypreParMatrix *ParBilinearForm::ParallelAssembleInternalMatrix()
-{
-   if (p_mat.Ptr() == NULL)
-   {
-      ParallelAssemble(p_mat, mat);
-   }
-   return p_mat.As<HypreParMatrix>();
-}
-
 void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
 {
   A.Clear();
@@ -342,15 +333,6 @@ void ParBilinearForm
   A.EliminateRowsCols(dof_list, X, B);
 }

-void ParBilinearForm::ParallelEliminateEssentialBC(
-   const Array<int> &bdr_attr_is_ess, const HypreParVector &X, HypreParVector &B)
-{
-   Array<int> dof_list;
-   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
-
-   p_mat.As<HypreParMatrix>()->EliminateRowsCols(dof_list, X, B);
-}
-
 HypreParMatrix *ParBilinearForm::
 ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                             HypreParMatrix &A) const
@@ -362,26 +344,6 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
   return A.EliminateRowsCols(dof_list);
 }

-void ParBilinearForm::ParallelEliminateEssentialBC(const Array<int>
-                                                   &bdr_attr_is_ess)
-{
-   Array<int> tdofs_list;
-   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, tdofs_list);
-
-   ParallelEliminateTDofs(tdofs_list);
-}
-
-void ParBilinearForm::ParallelEliminateTDofs(const Array<int> &tdofs_list)
-{
-   p_mat_e.EliminateRowsCols(p_mat, tdofs_list);
-}
-
-void ParBilinearForm::ParallelEliminateTDofsInRHS(
-   const Array<int> &tdofs_list, const Vector &x, Vector &b)
-{
-   p_mat.EliminateBC(p_mat_e, tdofs_list, x, b);
-}
-
 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const real_t a)
 const
 {
@@ -523,7 +485,7 @@ void ParBilinearForm::FormLinearSystem(
      HypreParVector true_X(pfes), true_B(pfes);
      P.MultTranspose(b, true_B);
      R.Mult(x, true_X);
-      ParallelEliminateTDofsInRHS(ess_tdof_list, true_X, true_B);
+      p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
      R.MultTranspose(true_B, b);
      hybridization->ReduceRHS(true_B, B);
      X.SetSize(B.Size());
@@ -536,11 +498,17 @@ void ParBilinearForm::FormLinearSystem(
      B.SetSize(X.Size());
      P.MultTranspose(b, B);
      R.Mult(x, X);
-      ParallelEliminateTDofsInRHS(ess_tdof_list, X, B);
+      p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
   }
 }

+void ParBilinearForm::EliminateVDofsInRHS(
+   const Array<int> &vdofs, const Vector &x, Vector &b)
+{
+   p_mat.EliminateBC(p_mat_e, vdofs, x, b);
+}
+
 void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
                                       OperatorHandle &A)
 {
@@ -585,7 +553,7 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
         mat = NULL;
         delete mat_e;
         mat_e = NULL;
-         ParallelEliminateTDofs(ess_tdof_list);
+         p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
      }
      if (hybridization)
      {
@@ -647,180 +615,36 @@ void ParBilinearForm::Update(FiniteElementSpace *nfes)
   p_mat_e.Clear();
 }

-void ParMixedBilinearForm::pAllocMat()
-{
-   const int trial_nbr_size = trial_pfes->GetFaceNbrVSize();
-   const int test_nbr_size = test_pfes->GetFaceNbrVSize();

-   if (keep_nbr_block)
-   {
-      mat = new SparseMatrix(height + test_nbr_size, width + trial_nbr_size);
-   }
-   else
-   {
-      mat = new SparseMatrix(height, width + trial_nbr_size);
-   }
+HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
+{
+   // construct the block-diagonal matrix A
+   HypreParMatrix *A =
+      new HypreParMatrix(trial_pfes->GetComm(),
+                         test_pfes->GlobalVSize(),
+                         trial_pfes->GlobalVSize(),
+                         test_pfes->GetDofOffsets(),
+                         trial_pfes->GetDofOffsets(),
+                         mat);
+
+   HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
+                             trial_pfes->Dof_TrueDof_Matrix());
+
+   delete A;
+
+   return rap;
 }

-void ParMixedBilinearForm::AssembleSharedFaces(int skip_zeros)
+void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
 {
-   ParMesh *pmesh = trial_pfes->GetParMesh();
-   FaceElementTransformations *T;
-   Array<int> tr_vdofs1, tr_vdofs2, tr_vdofs_all;
-   Array<int> te_vdofs1, te_vdofs2, te_vdofs_all;
-   DenseMatrix elemmat;
-
-   int nfaces = pmesh->GetNSharedFaces();
-   for (int i = 0; i < nfaces; i++)
-   {
-      T = pmesh->GetSharedFaceTransformations(i);
-      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
-      trial_pfes->GetElementVDofs(T->Elem1No, tr_vdofs1);
-      test_pfes->GetElementVDofs(T->Elem1No, te_vdofs1);
-      trial_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, tr_vdofs2);
-      test_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, te_vdofs2);
-
-      tr_vdofs1.Copy(tr_vdofs_all);
-      for (int j = 0; j < tr_vdofs2.Size(); j++)
-      {
-         if (tr_vdofs2[j] >= 0)
-         {
-            tr_vdofs2[j] += width;
-         }
-         else
-         {
-            tr_vdofs2[j] -= width;
-         }
-      }
-      tr_vdofs_all.Append(tr_vdofs2);
-
-      if (keep_nbr_block)
-      {
-         te_vdofs1.Copy(te_vdofs_all);
-         for (int j = 0; j < te_vdofs2.Size(); j++)
-         {
-            if (te_vdofs2[j] >= 0)
-            {
-               te_vdofs2[j] += height;
-            }
-            else
-            {
-               te_vdofs2[j] -= height;
-            }
-         }
-         te_vdofs_all.Append(te_vdofs2);
-      }
-
-      for (int k = 0; k < interior_face_integs.Size(); k++)
-      {
-         interior_face_integs[k]->
-         AssembleFaceMatrix(*trial_pfes->GetFE(T->Elem1No),
-                            *test_pfes->GetFE(T->Elem1No),
-                            *trial_pfes->GetFaceNbrFE(Elem2NbrNo),
-                            *test_pfes->GetFaceNbrFE(Elem2NbrNo),
-                            *T, elemmat);
-         if (keep_nbr_block)
-         {
-            mat->AddSubMatrix(te_vdofs_all, tr_vdofs_all, elemmat, skip_zeros);
-         }
-         else
-         {
-            mat->AddSubMatrix(te_vdofs1, tr_vdofs_all, elemmat, skip_zeros);
-         }
-      }
-   }
-}
-
-void ParMixedBilinearForm::Assemble(int skip_zeros)
-{
-   if (interior_face_integs.Size())
-   {
-      trial_pfes->ExchangeFaceNbrData();
-      test_pfes->ExchangeFaceNbrData();
-      if (!ext && mat == NULL)
-      {
-         pAllocMat();
-      }
-   }
-
-   MixedBilinearForm::Assemble(skip_zeros);
-
-   if (!ext && interior_face_integs.Size() > 0)
-   {
-      AssembleSharedFaces(skip_zeros);
-   }
-}
-
-HypreParMatrix *ParMixedBilinearForm::ParallelAssembleInternalMatrix()
-{
-   if (p_mat.Ptr() == NULL)
-   {
-      ParallelAssemble(p_mat, mat);
-   }
-   return p_mat.As<HypreParMatrix>();
-}
-
-HypreParMatrix *ParMixedBilinearForm::ParallelAssemble(SparseMatrix *m)
-{
-   OperatorHandle Mh(Operator::Hypre_ParCSR);
-   ParallelAssemble(Mh, m);
-   Mh.SetOperatorOwner(false);
-   return Mh.As<HypreParMatrix>();
-}
-
-void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A,
-                                            SparseMatrix *A_local)
-{
-   A.Clear();
-
-   if (A_local == NULL) { return; }
-   MFEM_VERIFY(A_local->Finalized(), "the local matrix must be finalized");
-
-   OperatorHandle dA(A.Type()), hdA;
-
-   if (interior_face_integs.Size() == 0)
-   {
-      // construct the rectangular block-diagonal matrix dA
-      dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
-                                  test_pfes->GlobalVSize(),
-                                  trial_pfes->GlobalVSize(),
-                                  test_pfes->GetDofOffsets(),
-                                  trial_pfes->GetDofOffsets(),
-                                  A_local);
-   }
-   else
-   {
-      // handle the case when 'a' contains off-diagonal
-      const int lvrows = test_pfes->GetVSize();
-      const int lvcols = trial_pfes->GetVSize();
-      const HYPRE_BigInt *face_nbr_glob_lcol = trial_pfes->GetFaceNbrGlobalDofMap();
-      const HYPRE_BigInt lcol_offset = trial_pfes->GetMyDofOffset();
-
-      Array<HYPRE_BigInt> glob_J(A_local->NumNonZeroElems());
-      const int *J = A_local->GetJ();
-      for (int i = 0; i < glob_J.Size(); i++)
-      {
-         if (J[i] < lvcols)
-         {
-            glob_J[i] = J[i] + lcol_offset;
-         }
-         else
-         {
-            glob_J[i] = face_nbr_glob_lcol[J[i] - lvcols];
-         }
-      }
-
-      // TODO - construct dA directly in the A format
-      hdA.Reset(
-         new HypreParMatrix(trial_pfes->GetComm(), lvrows, test_pfes->GlobalVSize(),
-                            trial_pfes->GlobalVSize(), A_local->GetI(), glob_J,
-                            A_local->GetData(), test_pfes->GetDofOffsets(),
-                            trial_pfes->GetDofOffsets()));
-      // - hdA owns the new HypreParMatrix
-      // - the above constructor copies all input arrays
-      glob_J.DeleteAll();
-      dA.ConvertFrom(hdA);
-   }
+   // construct the rectangular block-diagonal matrix dA
+   OperatorHandle dA(A.Type());
+   dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
+                               test_pfes->GlobalVSize(),
+                               trial_pfes->GlobalVSize(),
+                               test_pfes->GetDofOffsets(),
+                               trial_pfes->GetDofOffsets(),
+                               mat);

   OperatorHandle P_test(A.Type()), P_trial(A.Type());

@@ -846,44 +670,6 @@ void ParMixedBilinearForm::TrueAddMult(const Vector &x, Vector &y,
   test_pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Yaux, 1.0, y);
 }

-void ParMixedBilinearForm::ParallelEliminateTrialEssentialBC(
-   const Array<int> &bdr_attr_is_ess)
-{
-   Array<int> trial_tdof_list;
-   trial_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, trial_tdof_list);
-
-   ParallelEliminateTrialTDofs(trial_tdof_list);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTrialTDofs(
-   const Array<int> &trial_tdof_list)
-{
-   HypreParMatrix *temp = p_mat.As<HypreParMatrix>()->EliminateCols(
-                             trial_tdof_list);
-   p_mat_e.Reset(temp, true);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTrialTDofsInRHS(
-   const Array<int> &trial_tdof_list, const Vector &x, Vector &b)
-{
-   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, x, 1.0, b);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTestEssentialBC(
-   const Array<int> &bdr_attr_is_ess)
-{
-   Array<int> test_tdof_list;
-   test_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, test_tdof_list);
-
-   ParallelEliminateTestTDofs(test_tdof_list);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTestTDofs(
-   const Array<int> &test_tdof_list)
-{
-   p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
-}
-
 void ParMixedBilinearForm::FormRectangularSystemMatrix(
   const Array<int>
   &trial_tdof_list,
@@ -904,8 +690,10 @@ void ParMixedBilinearForm::FormRectangularSystemMatrix(
      mat = NULL;
      delete mat_e;
      mat_e = NULL;
-      ParallelEliminateTrialTDofs(trial_tdof_list);
-      ParallelEliminateTestTDofs(test_tdof_list);
+      HypreParMatrix *temp =
+         p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
+      p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
+      p_mat_e.Reset(temp, true);
   }

   A = p_mat;
@@ -935,7 +723,7 @@ void ParMixedBilinearForm::FormRectangularLinearSystem(
   test_P->MultTranspose(b, B);
   trial_R->Mult(x, X);

-   ParallelEliminateTrialTDofsInRHS(trial_tdof_list, X, B);
+   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
   B.SetSubVector(test_tdof_list, 0.0);
 }

@@ -73,7 +73,7 @@ public:
   /** When set to true and the ParBilinearForm has interior face integrators,
       the local SparseMatrix will include the rows (in addition to the columns)
       corresponding to face-neighbor dofs. The default behavior is to disregard
-       those rows. Must be called before the first Assemble() call. */
+       those rows. Must be called before the first Assemble call. */
   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }

   /** @brief Set the operator type id for the parallel matrix/operator when
@@ -101,14 +101,6 @@ public:
       diagonal for this case. */
   void AssembleDiagonal(Vector &diag) const override;

-   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
-   /** The returned matrix is the internal one, owned by the form. It is not
-       reassembled if it has been already constructed. If FormSystemMatrix()
-       has been called before, it is the system matrix with eliminated
-       essential DOFs, otherwise the parallel matrix is assembled here without
-       the elimination process. */
-   HypreParMatrix *ParallelAssembleInternalMatrix();
-
   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
   /** The returned matrix has to be deleted by the caller. */
   HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
@@ -154,13 +146,6 @@ public:
                                     const HypreParVector &X,
                                     HypreParVector &B) const;

-   /// Eliminate essential boundary DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. */
-   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
-                                     const HypreParVector &X,
-                                     HypreParVector &B);
-
   /// Eliminate essential boundary DOFs from a parallel assembled matrix @a A.
   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
       the essential part of the boundary. The eliminated part is stored in a
@@ -172,12 +157,6 @@ public:
   HypreParMatrix *ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                                                HypreParMatrix &A) const;

-   /// Eliminate essential boundary DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. This method relies on
-       ParallelEliminateTDofs(const Array<int> &), see it for details. */
-   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess);
-
   /// Eliminate essential true DOFs from a parallel assembled matrix @a A.
   /** Given a list of essential true dofs and the parallel assembled matrix
       @a A, eliminate the true dofs from the matrix, storing the eliminated
@@ -190,28 +169,6 @@ public:
                                          HypreParMatrix &A) const
   { return A.EliminateRowsCols(tdofs_list); }

-   /// Eliminate essential true DOFs from the parallel system matrix.
-   /** Given a list of essential true dofs, eliminate the true dofs from
-       the parallel assembled system matrix, storing the eliminated part
-       internally. This method works in conjunction with
-       ParallelEliminateTDofsInRHS() and allows elimination of boundary
-       conditions in multiple right-hand sides. */
-   void ParallelEliminateTDofs(const Array<int> &tdofs_list);
-
-   /** @brief Use the stored eliminated part of the parallel system matrix for
-       elimination of boundary conditions in the r.h.s. */
-   /** Given a list of essential true dofs, eliminate the true dofs from the
-       right-hand side @a b using the solution vector @a x and the previously
-       stored eliminated part of the parallel assembled system matrix produced
-       by ParallelEliminateTDofs(const Array<int> &). */
-   void ParallelEliminateTDofsInRHS(const Array<int> &tdofs, const Vector &x,
-                                    Vector &b);
-
-   /// @deprecated Use ParallelEliminateTDofsInRHS() instead.
-   MFEM_DEPRECATED void EliminateVDofsInRHS(const Array<int> &vdofs,
-                                            const Vector &x, Vector &b)
-   { ParallelEliminateTDofsInRHS(vdofs, x, b); }
-
   /** @brief Compute @a y += @a a (P^t A P) @a x, where @a x and @a y are
       vectors on the true dofs. */
   void TrueAddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
@@ -281,6 +238,8 @@ public:

   void Update(FiniteElementSpace *nfes = NULL) override;

+   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
+
   virtual ~ParBilinearForm() { }
 };

@@ -298,13 +257,6 @@ protected:
   /// Matrix and eliminated matrix
   OperatorHandle p_mat, p_mat_e;

-   bool keep_nbr_block;
-
-   // Allocate mat - called when (mat == NULL && fbfi.Size() > 0)
-   void pAllocMat();
-
-   void AssembleSharedFaces(int skip_zeros = 1);
-
 private:
   /// Copy construction is not supported; body is undefined.
   ParMixedBilinearForm(const ParMixedBilinearForm &);
@@ -324,7 +276,6 @@ public:
   {
      trial_pfes = trial_fes;
      test_pfes  = test_fes;
-      keep_nbr_block = false;
   }

   /** @brief Create a ParMixedBilinearForm on the given FiniteElementSpace%s
@@ -344,89 +295,15 @@ public:
   {
      trial_pfes = trial_fes;
      test_pfes  = test_fes;
-      keep_nbr_block = false;
   }

-   /** When set to true and the ParMixedBilinearForm has interior face
-       integrators, the local SparseMatrix will include the rows (in addition
-       to the columns) corresponding to face-neighbor dofs. The default
-       behavior is to disregard those rows. Must be called before the first
-       Assemble() call. */
-   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
-
-   /// Assemble the local matrix
-   void Assemble(int skip_zeros = 1);
-
   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
-   /** The returned matrix is the internal one, owned by the form. It is not
-       reassembled if it has been already constructed. If
-       FormRectangularSystemMatrix() has been called before, it is the system
-       matrix with eliminated essential DOFs, otherwise the parallel matrix is
-       assembled here without the elimination process. */
-   HypreParMatrix *ParallelAssembleInternalMatrix();
-
-   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
-   /** The returned matrix has to be deleted by the caller. */
-   HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
-
-   /** @brief Returns the eliminated matrix assembled on the true dofs, i.e.
-       P_test^t A_local P_trial. */
-   /** The returned matrix has to be deleted by the caller. */
-   HypreParMatrix *ParallelAssembleElim() { return ParallelAssemble(mat_e); }
-
-   /** @brief Return the matrix @a m assembled on the true dofs, i.e. P_test^t
-       A_local P_trial. */
-   /** The returned matrix has to be deleted by the caller. */
-   HypreParMatrix *ParallelAssemble(SparseMatrix *m);
+   HypreParMatrix *ParallelAssemble();

   /** @brief Returns the matrix assembled on the true dofs, i.e.
       @a A = P_test^t A_local P_trial, in the format (type id) specified by
       @a A. */
-   void ParallelAssemble(OperatorHandle &A) { ParallelAssemble(A, mat); }
-
-   /** Returns the eliminated matrix assembled on the true dofs, i.e.
-       @a A_elim = P^t A_elim_local P in the format (type id) specified by @a A.
-    */
-   void ParallelAssembleElim(OperatorHandle &A_elim)
-   { ParallelAssemble(A_elim, mat_e); }
-
-   /** Returns the matrix @a A_local assembled on the true dofs, i.e.
-       @a A = P_test^t A_local P_trial in the format (type id) specified by
-       @a A. */
-   void ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local);
-
-   /// Eliminate essential boundary trial DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. This method relies on
-       ParallelEliminateTrialTDofs(const Array<int> &), see it for details. */
-   void ParallelEliminateTrialEssentialBC(const Array<int> &bdr_attr_is_ess);
-
-   /// Eliminate essential trial true DOFs from the parallel system matrix.
-   /** Given a list of essential trial true dofs, eliminate the trial true dofs
-       from the parallel assembled system matrix, storing the eliminated part
-       internally. This method works in conjunction with
-       ParallelEliminateTrialTDofsInRHS() and allows elimination of boundary
-       conditions in multiple right-hand sides. */
-   void ParallelEliminateTrialTDofs(const Array<int> &trial_tdof_list);
-
-   /** @brief Use the stored eliminated part of the parallel system matrix for
-       elimination of boundary conditions in the r.h.s. */
-   /** Given a list of essential trial true dofs, eliminate the trial true dofs
-       from the right-hand side @a B using the solution vector @a X and the
-       previously stored eliminated part of the parallel assembled system
-       matrix produced by ParallelEliminateTrialTDofs(const Array<int> &). */
-   void ParallelEliminateTrialTDofsInRHS(const Array<int> &trial_tdof_list,
-                                         const Vector &X, Vector &B);
-
-   /// Eliminate essential boundary test DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. */
-   void ParallelEliminateTestEssentialBC(const Array<int> &bdr_attr_is_ess);
-
-   /// Eliminate essential test true DOFs from the parallel system matrix.
-   /** Given a list of essential test true dofs, eliminate the test true dofs
-       from the parallel assembled system matrix. */
-   void ParallelEliminateTestTDofs(const Array<int> &test_tdof_list);
+   void ParallelAssemble(OperatorHandle &A);

   using MixedBilinearForm::FormRectangularSystemMatrix;
   using MixedBilinearForm::FormRectangularLinearSystem;
@@ -105,59 +105,6 @@ const SparseMatrix &ParNonlinearForm::GetLocalGradient(const Vector &x) const
   return *Grad;
 }

-void ParNonlinearForm::GradientSharedFaces(const Vector &x,
-                                           int skip_zeros) const
-{
-   ParFiniteElementSpace *pfes = ParFESpace();
-   ParMesh *pmesh = pfes->GetParMesh();
-   FaceElementTransformations *T;
-   Array<int> vdofs1, vdofs2, vdofs_all;
-   DenseMatrix elemmat;
-   Vector el_x, nbr_x, face_x;
-   const Vector &px = Prolongate(x);
-
-   ParGridFunction pgf(pfes, const_cast<Vector&>(px), 0);
-   pgf.ExchangeFaceNbrData();
-
-   int nfaces = pmesh->GetNSharedFaces();
-   for (int i = 0; i < nfaces; i++)
-   {
-      T = pmesh->GetSharedFaceTransformations(i);
-      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
-
-      pfes->GetElementVDofs(T->Elem1No, vdofs1);
-      pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
-      face_x.SetSize(vdofs1.Size() + vdofs2.Size());
-
-      el_x.MakeRef(face_x, 0, vdofs1.Size());
-      pgf.GetSubVector(vdofs1, el_x);
-
-      nbr_x.MakeRef(face_x, vdofs1.Size(), vdofs2.Size());
-      pgf.FaceNbrData().GetSubVector(vdofs2, nbr_x);
-
-      vdofs1.Copy(vdofs_all);
-      for (int j = 0; j < vdofs2.Size(); j++)
-      {
-         if (vdofs2[j] >= 0)
-         {
-            vdofs2[j] += height;
-         }
-         else
-         {
-            vdofs2[j] -= height;
-         }
-      }
-      vdofs_all.Append(vdofs2);
-      for (int k = 0; k < fnfi.Size(); k++)
-      {
-         fnfi[k]->AssembleFaceGrad(*pfes->GetFE(T->Elem1No),
-                                   *pfes->GetFaceNbrFE(Elem2NbrNo),
-                                   *T, face_x, elemmat);
-         Grad->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
-      }
-   }
-}
-
 Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 {
   if (NonlinearForm::ext) { return NonlinearForm::GetGradient(x); }
@@ -165,61 +112,19 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
   ParFiniteElementSpace *pfes = ParFESpace();

   pGrad.Clear();
-   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type()), hdA;

-   if (fnfi.Size())
+   NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
+
+   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
+
+   if (fnfi.Size() == 0)
   {
-      const int skip_zeros = 0;
-
-      pfes->ExchangeFaceNbrData();
-      if (Grad == NULL)
-      {
-         int nbr_size = pfes->GetFaceNbrVSize();
-         Grad = new SparseMatrix(pfes->GetVSize(), pfes->GetVSize() + nbr_size);
-      }
-
-      NonlinearForm::GetGradient(x, false); // (re)assemble Grad, no b.c.
-
-      GradientSharedFaces(x, skip_zeros);
-
-      Grad->Finalize(skip_zeros);
-
-      // handle the case when 'a' contains off-diagonal
-      int lvsize = pfes->GetVSize();
-      const HYPRE_BigInt *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
-      HYPRE_BigInt ldof_offset = pfes->GetMyDofOffset();
-
-      Array<HYPRE_BigInt> glob_J(Grad->NumNonZeroElems());
-      int *J = Grad->GetJ();
-      for (int i = 0; i < glob_J.Size(); i++)
-      {
-         if (J[i] < lvsize)
-         {
-            glob_J[i] = J[i] + ldof_offset;
-         }
-         else
-         {
-            glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
-         }
-      }
-
-      // TODO - construct dA directly in the A format
-      hdA.Reset(
-         new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
-                            pfes->GlobalVSize(), Grad->GetI(), glob_J,
-                            Grad->GetData(), pfes->GetDofOffsets(),
-                            pfes->GetDofOffsets()));
-      // - hdA owns the new HypreParMatrix
-      // - the above constructor copies all input arrays
-      glob_J.DeleteAll();
-      dA.ConvertFrom(hdA);
+      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
+                             pfes->GetDofOffsets(), Grad);
   }
   else
   {
-      NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
-
-      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
-                             pfes->GetDofOffsets(), Grad);
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
   }

   // RAP the local gradient dA.
@@ -366,70 +271,7 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const

   if (fnfi.Size() > 0)
   {
-      // Terms over shared interior faces in parallel.
-      ParMesh *pmesh = ParFESpace(0)->GetParMesh();
-      FaceElementTransformations *tr;
-
-      Array<Array<int> *>vdofs(fes.Size());
-      Array<Array<int> *>vdofs2(fes.Size());
-      Array<Vector *> el_x(fes.Size());
-      Array<const Vector *> el_x_const(fes.Size());
-      Array<Vector *> el_y(fes.Size());
-      Array<const FiniteElement *> fe(fes.Size());
-      Array<const FiniteElement *> fe2(fes.Size());
-      Array<ParGridFunction *> pgfs(fes.Size());
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         el_x_const[s] = el_x[s] = new Vector();
-         el_y[s] = new Vector();
-         vdofs[s] = new Array<int>;
-         vdofs2[s] = new Array<int>;
-         pgfs[s] = new ParGridFunction(const_cast<ParFiniteElementSpace*>(ParFESpace(s)),
-                                       xs.GetBlock(s));
-         pgfs[s]->ExchangeFaceNbrData();
-      }
-
-      const int n_shared_faces = pmesh->GetNSharedFaces();
-      for (int i = 0; i < n_shared_faces; i++)
-      {
-         tr = pmesh->GetSharedFaceTransformations(i, true);
-         int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
-
-         for (int s=0; s<fes.Size(); ++s)
-         {
-            const ParFiniteElementSpace *pfes = ParFESpace(s);
-            fe[s] = pfes->GetFE(tr->Elem1No);
-            fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
-
-            pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
-            pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
-
-            el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
-            xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
-            pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
-                                                el_x[s]->GetData() + vdofs[s]->Size());
-         }
-
-         for (int k = 0; k < fnfi.Size(); ++k)
-         {
-            fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
-
-            for (int s=0; s<fes.Size(); ++s)
-            {
-               if (el_y[s]->Size() == 0) { continue; }
-               ys.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
-            }
-         }
-      }
-
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         delete pgfs[s];
-         delete vdofs2[s];
-         delete vdofs[s];
-         delete el_y[s];
-         delete el_x[s];
-      }
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
   }

   for (int s=0; s<fes.Size(); ++s)
@@ -486,106 +328,6 @@ void ParBlockNonlinearForm::SetGradientType(Operator::Type tid)
   }
 }

-void ParBlockNonlinearForm::GradientSharedFaces(const BlockVector &xs,
-                                                int skip_zeros) const
-{
-   // Terms over shared interior faces in parallel.
-   ParMesh *pmesh = ParFESpace(0)->GetParMesh();
-   FaceElementTransformations *tr;
-
-   Array<Array<int> *>vdofs(fes.Size());
-   Array<Array<int> *>vdofs2(fes.Size());
-   Array<Array<int> *>vdofs_all(fes.Size());
-   Array<Vector *> el_x(fes.Size());
-   Array<const Vector *> el_x_const(fes.Size());
-   Array2D<DenseMatrix *> elmats(fes.Size(), fes.Size());
-   Array<const FiniteElement *> fe(fes.Size());
-   Array<const FiniteElement *> fe2(fes.Size());
-   Array<ParGridFunction *> pgfs(fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      el_x_const[s1] = el_x[s1] = new Vector();
-      vdofs[s1] = new Array<int>;
-      vdofs2[s1] = new Array<int>;
-      vdofs_all[s1] = new Array<int>;
-      pgfs[s1] = new ParGridFunction(
-         const_cast<ParFiniteElementSpace*>(ParFESpace(s1)),
-         const_cast<Vector&>(xs.GetBlock(s1)));
-      pgfs[s1]->ExchangeFaceNbrData();
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         elmats(s1,s2) = new DenseMatrix();
-      }
-   }
-
-   const int n_shared_faces = pmesh->GetNSharedFaces();
-   for (int i = 0; i < n_shared_faces; i++)
-   {
-      tr = pmesh->GetSharedFaceTransformations(i, true);
-      int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
-
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         const ParFiniteElementSpace *pfes = ParFESpace(s);
-         fe[s] = pfes->GetFE(tr->Elem1No);
-         fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
-
-         pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
-         pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
-
-         el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
-         xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
-         pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
-                                             el_x[s]->GetData() + vdofs[s]->Size());
-
-         vdofs[s]->Copy(*vdofs_all[s]);
-
-         const int lvsize = pfes->GetVSize();
-         for (int j = 0; j < vdofs2[s]->Size(); j++)
-         {
-            if ((*vdofs2[s])[j] >= 0)
-            {
-               (*vdofs2[s])[j] += lvsize;
-            }
-            else
-            {
-               (*vdofs2[s])[j] -= lvsize;
-            }
-         }
-         vdofs_all[s]->Append(*(vdofs2[s]));
-      }
-
-      for (int k = 0; k < fnfi.Size(); ++k)
-      {
-         fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
-
-         for (int s1=0; s1<fes.Size(); ++s1)
-         {
-            for (int s2=0; s2<fes.Size(); ++s2)
-            {
-               if (elmats(s1,s2)->Height() == 0) { continue; }
-               Grads(s1,s2)->AddSubMatrix(*vdofs[s1], *vdofs_all[s2],
-                                          *elmats(s1,s2), skip_zeros);
-            }
-         }
-      }
-   }
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      delete pgfs[s1];
-      delete vdofs_all[s1];
-      delete vdofs2[s1];
-      delete vdofs[s1];
-      delete el_x[s1];
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete elmats(s1,s2);
-      }
-   }
-}
-
 BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
 {
   if (pBlockGrad == NULL)
@@ -605,155 +347,49 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
      }
   }

-   // xs_true is not modified, so const_cast is okay
-   xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
+   GetLocalGradient(x); // gradients are stored in 'Grads'

   if (fnfi.Size() > 0)
   {
-      const int skip_zeros = 0;
-
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         const_cast<ParFiniteElementSpace*>(pfes[s])->ExchangeFaceNbrData();
-      }
-
-      for (int s1=0; s1<fes.Size(); ++s1)
-      {
-         for (int s2=0; s2<fes.Size(); ++s2)
-         {
-            if (Grads(s1,s2) == NULL)
-            {
-               int nbr_size = pfes[s2]->GetFaceNbrVSize();
-               Grads(s1,s2) = new SparseMatrix(pfes[s1]->GetVSize(),
-                                               pfes[s2]->GetVSize() + nbr_size);
-            }
-         }
-      }
-
-      // (re)assemble Grad without b.c. into 'Grads'
-      BlockNonlinearForm::ComputeGradientBlocked(xs, false);
-
-      GradientSharedFaces(xs, skip_zeros);
-
-      // finalize the gradients
-      for (int s1=0; s1<fes.Size(); ++s1)
-         for (int s2=0; s2<fes.Size(); ++s2)
-         {
-            Grads(s1,s2)->Finalize(skip_zeros);
-         }
-
-      for (int s1=0; s1<fes.Size(); ++s1)
-      {
-         for (int s2=0; s2<fes.Size(); ++s2)
-         {
-            OperatorHandle hdA;
-            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                           Ph(phBlockGrad(s1,s2)->Type()),
-                           Rh(phBlockGrad(s1,s2)->Type());
-
-            // handle the case when 'a' contains off-diagonal
-            int lvsize = pfes[s2]->GetVSize();
-            const HYPRE_BigInt *face_nbr_glob_ldof =
-               const_cast<ParFiniteElementSpace*>(pfes[s2])->GetFaceNbrGlobalDofMap();
-            HYPRE_BigInt ldof_offset = pfes[s2]->GetMyDofOffset();
-
-            Array<HYPRE_BigInt> glob_J(Grads(s1,s2)->NumNonZeroElems());
-            int *J = Grads(s1,s2)->GetJ();
-            for (int i = 0; i < glob_J.Size(); i++)
-            {
-               if (J[i] < lvsize)
-               {
-                  glob_J[i] = J[i] + ldof_offset;
-               }
-               else
-               {
-                  glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
-               }
-            }
-
-            // TODO - construct dA directly in the A format
-            hdA.Reset(
-               new HypreParMatrix(pfes[s2]->GetComm(), pfes[s1]->GetVSize(),
-                                  pfes[s1]->GlobalVSize(), pfes[s2]->GlobalVSize(),
-                                  Grads(s1,s2)->GetI(), glob_J, Grads(s1,s2)->GetData(),
-                                  pfes[s1]->GetDofOffsets(), pfes[s2]->GetDofOffsets()));
-            // - hdA owns the new HypreParMatrix
-            // - the above constructor copies all input arrays
-            glob_J.DeleteAll();
-            dA.ConvertFrom(hdA);
-
-            if (s1 == s2)
-            {
-               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-               OperatorHandle Ae;
-               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-            }
-            else
-            {
-               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-            }
-
-            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
-         }
-      }
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
   }
-   else
+
+   for (int s1=0; s1<fes.Size(); ++s1)
   {
-      // (re)assemble Grad without b.c. into 'Grads'
-      BlockNonlinearForm::ComputeGradientBlocked(xs);
-
-      for (int s1=0; s1<fes.Size(); ++s1)
+      for (int s2=0; s2<fes.Size(); ++s2)
      {
-         for (int s2=0; s2<fes.Size(); ++s2)
+         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                        Ph(phBlockGrad(s1,s2)->Type()),
+                        Rh(phBlockGrad(s1,s2)->Type());
+
+         if (s1 == s2)
         {
-            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                           Ph(phBlockGrad(s1,s2)->Type()),
-                           Rh(phBlockGrad(s1,s2)->Type());
+            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
+                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
+            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);

-            if (s1 == s2)
-            {
-               dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                      pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-               OperatorHandle Ae;
-               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-            }
-            else
-            {
-               dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                           pfes[s1]->GlobalVSize(),
-                                           pfes[s2]->GlobalVSize(),
-                                           pfes[s1]->GetDofOffsets(),
-                                           pfes[s2]->GetDofOffsets(),
-                                           Grads(s1,s2));
-               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-            }
-
-            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+            OperatorHandle Ae;
+            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
         }
+         else
+         {
+            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
+                                        pfes[s1]->GlobalVSize(),
+                                        pfes[s2]->GlobalVSize(),
+                                        pfes[s1]->GetDofOffsets(),
+                                        pfes[s2]->GetDofOffsets(),
+                                        Grads(s1,s2));
+            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+         }
+
+         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
      }
   }

@@ -29,8 +29,6 @@ protected:
   mutable ParGridFunction X, Y;
   mutable OperatorHandle pGrad;

-   void GradientSharedFaces(const Vector &x, int skip_zeros = 1) const;
-
 public:
   ParNonlinearForm(ParFiniteElementSpace *pf);

@@ -83,8 +81,6 @@ protected:
   mutable Array2D<OperatorHandle *> phBlockGrad;
   mutable BlockOperator *pBlockGrad;

-   void GradientSharedFaces(const BlockVector &xs, int skip_zeros) const;
-
 public:
   /// Computes the energy of the system
   real_t GetEnergy(const Vector &x) const override;
@@ -561,8 +561,7 @@ void CopyMemory(Memory<T> &src, Memory<T> &dst, MemoryClass dst_mc,
    this function. In particular, @a dst should be empty or deleted before
    calling this function. */
 template <typename SrcT, typename DstT>
-void CopyConvertMemory(const Memory<SrcT> &src, MemoryClass dst_mc,
-                       Memory<DstT> &dst)
+void CopyConvertMemory(Memory<SrcT> &src, MemoryClass dst_mc, Memory<DstT> &dst)
 {
   auto capacity = src.Capacity();
   dst.New(capacity, GetMemoryType(dst_mc));
@@ -843,8 +842,8 @@ static int GetPartitioningArraySize(MPI_Comm comm)
 ///
 /// Both @a row and @a col are partitioning arrays, whose length is returned by
 /// GetPartitioningArraySize(), see @ref hypre_partitioning_descr.
-static bool RowAndColStartsAreEqual(MPI_Comm comm, const HYPRE_BigInt *rows,
-                                    const HYPRE_BigInt *cols)
+static bool RowAndColStartsAreEqual(MPI_Comm comm, HYPRE_BigInt *rows,
+                                    HYPRE_BigInt *cols)
 {
   const int part_size = GetPartitioningArraySize(comm);
   bool are_equal = true;
@@ -1132,7 +1131,7 @@ HypreParMatrix::HypreParMatrix(
 HypreParMatrix::HypreParMatrix(MPI_Comm comm,
                               HYPRE_BigInt *row_starts,
                               HYPRE_BigInt *col_starts,
-                               const SparseMatrix *sm_a)
+                               SparseMatrix *sm_a)
 {
   MFEM_ASSERT(sm_a != NULL, "invalid input");
   MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
@@ -1146,7 +1145,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,

   hypre_CSRMatrixSetDataOwner(csr_a,0);
   MemoryIJData mem_a;
-   CopyCSR(const_cast<SparseMatrix*>(sm_a), mem_a, csr_a, false);
+   CopyCSR(sm_a, mem_a, csr_a, false);
   hypre_CSRMatrixSetRownnz(csr_a);

   // NOTE: this call creates a matrix on host even when device support is
@@ -1308,11 +1307,10 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows,
                               HYPRE_BigInt glob_nrows,
                               HYPRE_BigInt glob_ncols,
-                               const int *I,
-                               const HYPRE_BigInt *J,
-                               const real_t *data,
-                               const HYPRE_BigInt *rows,
-                               const HYPRE_BigInt *cols)
+                               int *I, HYPRE_BigInt *J,
+                               real_t *data,
+                               HYPRE_BigInt *rows,
+                               HYPRE_BigInt *cols)
 {
   Init();

@@ -565,7 +565,7 @@ public:
       partitioning arrays @a row_starts and @a col_starts. */
   HypreParMatrix(MPI_Comm comm, HYPRE_BigInt *row_starts,
                  HYPRE_BigInt *col_starts,
-                  const SparseMatrix *a); // constructor with 4 arguments, v2
+                  SparseMatrix *a); // constructor with 4 arguments, v2

   /// Creates boolean block-diagonal rectangular parallel matrix.
   /** The new HypreParMatrix does not take ownership of any of the input
@@ -594,9 +594,9 @@ public:
       arrays (so they can be deleted). See @ref hypre_partitioning_descr "here"
       for a description of the partitioning arrays @a rows and @a cols. */
   HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_BigInt glob_nrows,
-                  HYPRE_BigInt glob_ncols, const int *I, const HYPRE_BigInt *J,
-                  const real_t *data, const HYPRE_BigInt *rows,
-                  const HYPRE_BigInt *cols); // constructor with 9 arguments
+                  HYPRE_BigInt glob_ncols, int *I, HYPRE_BigInt *J,
+                  real_t *data, HYPRE_BigInt *rows,
+                  HYPRE_BigInt *cols); // constructor with 9 arguments

   /** @brief Copy constructor for a ParCSR matrix which creates a deep copy of
       structure and data from @a P. */
@@ -578,7 +578,7 @@ public:

       Presently, this method is used by SUNDIALS ARKStep integrator, for more
       details, see the ARKode User Guide. */
-   MFEM_DEPRECATED virtual int SUNMassSetup();
+   virtual int SUNMassSetup();

   /** @brief Solve the mass matrix linear system  M @a x = @a b, where M is
       defined by the method SUNMassSetup().
@@ -591,7 +591,7 @@ public:

       Presently, this method is used by SUNDIALS ARKStep integrator, for more
       details, see the ARKode User Guide. */
-   MFEM_DEPRECATED virtual int SUNMassSolve(const Vector &b, Vector &x, real_t tol);
+   virtual int SUNMassSolve(const Vector &b, Vector &x, real_t tol);

   /** @brief Compute the mass matrix-vector product @a v = M @a x, where M is
       defined by the method SUNMassSetup().
@@ -603,7 +603,7 @@ public:

       Presently, this method is used by SUNDIALS ARKStep integrator, for more
       details, see the ARKode User Guide. */
-   MFEM_DEPRECATED virtual int SUNMassMult(const Vector &x, Vector &v);
+   virtual int SUNMassMult(const Vector &x, Vector &v);

   virtual ~TimeDependentOperator() { }
 };
@@ -1367,84 +1367,6 @@ CVODESSolver::~CVODESSolver()
 // ARKStep interface
 // ---------------------------------------------------------------------------

-ARKStepSolver::TimeDependentOperatorWrapper::TimeDependentOperatorWrapper(
-   TimeDependentOperator *f)
-{
-   tdo = f;
-}
-
-int ARKStepSolver::TimeDependentOperatorWrapper::ARKSize() const
-{
-   return tdo->Height();
-}
-
-bool ARKStepSolver::TimeDependentOperatorWrapper::ARKInMassForm() const
-{
-   return (tdo->isExplicit() == false);
-}
-
-void ARKStepSolver::TimeDependentOperatorWrapper::ARKSetEvalMode(
-   const ARKEvalMode new_eval_mode)
-{
-   if (new_eval_mode == NORMAL)
-      tdo->SetEvalMode(tdo->NORMAL);
-   else if (new_eval_mode == ADDITIVE_TERM_1)
-      tdo->SetEvalMode(tdo->ADDITIVE_TERM_1);
-   else if (new_eval_mode == ADDITIVE_TERM_2)
-      tdo->SetEvalMode(tdo->ADDITIVE_TERM_2);
-   else
-      mfem_error("Unrecognized evaluation mode.");
-}
-
-void ARKStepSolver::TimeDependentOperatorWrapper::ARKEvaluateRHS(
-   const Vector &u, const real_t t, Vector &result) const
-{
-   tdo->SetTime(t);
-   if (ARKInMassForm())
-      tdo->Mult(u, result);
-   else
-      tdo->ExplicitMult(u, result);
-}
-
-int ARKStepSolver::TimeDependentOperatorWrapper::ARKImplicitSetup(
-   const Vector &u, const real_t t, const Vector &v, int jok, int *jcur,
-   real_t gamma)
-{
-   tdo->SetTime(t);
-   return tdo->SUNImplicitSetup(u, v, jok, jcur, gamma);
-}
-
-int ARKStepSolver::TimeDependentOperatorWrapper::ARKImplicitSolve(
-   const Vector &r, Vector &dk, real_t tol)
-{
-   return tdo->SUNImplicitSolve(r, dk, tol);
-}
-
-int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassSetup(const real_t t)
-{
-   tdo->SetTime(t);
-   return tdo->SUNMassSetup();
-}
-
-int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassSolve(const Vector &b,
-   Vector &x, real_t tol)
-{
-   return tdo->SUNMassSolve(b, x, tol);
-}
-
-int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassMult(const Vector &x,
-   const real_t t, Vector &v)
-{
-   tdo->SetTime(t);
-   return tdo->SUNMassMult(x, v);
-}
-
-int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassMult(const Vector &x,
-   Vector &v)
-{
-   return tdo->SUNMassMult(x, v);
-}
-
 int ARKStepSolver::RHS1(sunrealtype t, const N_Vector y, N_Vector result,
                        void *user_data)
 {
@@ -1459,11 +1381,19 @@ int ARKStepSolver::RHS1(sunrealtype t, const N_Vector y, N_Vector result,
   // or fe(t, y) in one of
   //   1. y' = fe(t, y) + fi(t, y)
   //   2. M y' = fe(t, y) + fi(t, y)
+   self->f->SetTime(t);
   if (self->rk_type == IMEX)
   {
-      self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_1);
+      self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_1);
+   }
+   if (self->f->isExplicit()) // ODE is in form 1
+   {
+      self->f->Mult(mfem_y, mfem_result);
+   }
+   else // ODE is in form 2
+   {
+      self->f->ExplicitMult(mfem_y, mfem_result);
   }
-   self->f_arkstep->ARKEvaluateRHS(mfem_y, t, mfem_result);

   // Return success
   return (0);
@@ -1480,8 +1410,16 @@ int ARKStepSolver::RHS2(sunrealtype t, const N_Vector y, N_Vector result,
   // Compute fi(t, y) in one of
   //   1. y' = fe(t, y) + fi(t, y)       (ODE is expressed in EXPLICIT form)
   //   2. M y' = fe(t, y) + fi(y, t)     (ODE is expressed in IMPLICIT form)
-   self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_2);
-   self->f_arkstep->ARKEvaluateRHS(mfem_y, t, mfem_result);
+   self->f->SetTime(t);
+   self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
+   if (self->f->isExplicit())
+   {
+      self->f->Mult(mfem_y, mfem_result);
+   }
+   else
+   {
+      self->f->ExplicitMult(mfem_y, mfem_result);
+   }

   // Return success
   return (0);
@@ -1498,11 +1436,12 @@ int ARKStepSolver::LinSysSetup(sunrealtype t, N_Vector y, N_Vector fy,
   ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(A));

   // Compute the linear system
+   self->f->SetTime(t);
   if (self->rk_type == IMEX)
   {
-      self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_2);
+      self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
   }
-   return (self->f_arkstep->ARKImplicitSetup(mfem_y, t, mfem_fy, jok, jcur, gamma));
+   return (self->f->SUNImplicitSetup(mfem_y, mfem_fy, jok, jcur, gamma));
 }

 int ARKStepSolver::LinSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
@@ -1515,9 +1454,9 @@ int ARKStepSolver::LinSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
   // Solve the linear system
   if (self->rk_type == IMEX)
   {
-      self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_2);
+      self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
   }
-   return (self->f_arkstep->ARKImplicitSolve(mfem_b, mfem_x, tol));
+   return (self->f->SUNImplicitSolve(mfem_b, mfem_x, tol));
 }

 int ARKStepSolver::MassSysSetup(sunrealtype t, SUNMatrix M,
@@ -1526,7 +1465,8 @@ int ARKStepSolver::MassSysSetup(sunrealtype t, SUNMatrix M,
   ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(M));

   // Compute the mass matrix system
-   return (self->f_arkstep->ARKMassSetup(t));
+   self->f->SetTime(t);
+   return (self->f->SUNMassSetup());
 }

 int ARKStepSolver::MassSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
@@ -1537,7 +1477,7 @@ int ARKStepSolver::MassSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
   ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(LS));

   // Solve the mass matrix system
-   return (self->f_arkstep->ARKMassSolve(mfem_b, mfem_x, tol));
+   return (self->f->SUNMassSolve(mfem_b, mfem_x, tol));
 }

 int ARKStepSolver::MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
@@ -1547,7 +1487,7 @@ int ARKStepSolver::MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
   ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(M));

   // Compute the mass matrix-vector product
-   return (self->f_arkstep->ARKMassMult(mfem_x, mfem_v));
+   return (self->f->SUNMassMult(mfem_x, mfem_v));
 }

 int ARKStepSolver::MassMult2(N_Vector x, N_Vector v, sunrealtype t,
@@ -1558,7 +1498,8 @@ int ARKStepSolver::MassMult2(N_Vector x, N_Vector v, sunrealtype t,
   ARKStepSolver *self = static_cast<ARKStepSolver*>(mtimes_data);

   // Compute the mass matrix-vector product
-   return (self->f_arkstep->ARKMassMult(mfem_x, t, mfem_v));
+   self->f->SetTime(t);
+   return (self->f->SUNMassMult(mfem_x, mfem_v));
 }

 ARKStepSolver::ARKStepSolver(Type type)
@@ -1577,12 +1518,13 @@ ARKStepSolver::ARKStepSolver(MPI_Comm comm, Type type)
 }
 #endif

-void ARKStepSolver::Init(ARKStepODE *f_ark_)
+void ARKStepSolver::Init(TimeDependentOperator &f_)
 {
-   f_arkstep = f_ark_;
+   // Initialize the base class
+   ODESolver::Init(f_);

   // Get the vector length
-   long local_size = f_arkstep->ARKSize();
+   long local_size = f_.Height();
 #ifdef MFEM_USE_MPI
   long global_size;
 #endif
@@ -1596,7 +1538,7 @@ void ARKStepSolver::Init(ARKStepODE *f_ark_)
   }

   // Get current time
-   double t = f ? f->GetTime() : 0.0;
+   double t = f_.GetTime();

   if (sundials_mem)
   {
@@ -1675,12 +1617,6 @@ void ARKStepSolver::Init(ARKStepODE *f_ark_)
   reinit = true;
 }

-void ARKStepSolver::Init(TimeDependentOperator &f_)
-{
-   f_tdo = std::make_unique<TimeDependentOperatorWrapper>(&f_);
-   Init(f_tdo.get());
-}
-
 void ARKStepSolver::Step(Vector &x, real_t &t, real_t &dt)
 {
   Y->MakeRef(x, 0, x.Size());
@@ -1773,9 +1709,6 @@ void ARKStepSolver::UseSundialsLinearSolver()

 void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
 {
-   // Check that the ODE is expressed in mass form
-   MFEM_VERIFY(f_arkstep->ARKInMassForm(), "ODE operator is not in mass form.")
-
   // Free any existing matrix and linear solver
   if (M != NULL)   { SUNMatDestroy(M); M = NULL; }
   if (LSM != NULL) { SUNLinSolFree(LSM); LSM = NULL; }
@@ -1806,13 +1739,13 @@ void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
   flag = MFEM_ARKode(SetMassFn)(sundials_mem, ARKStepSolver::MassSysSetup);
   MFEM_VERIFY(flag == ARK_SUCCESS,
               "error in " STR(MFEM_ARKode(SetMassFn)) "()");
+
+   // Check that the ODE is not expressed in EXPLICIT form
+   MFEM_VERIFY(!f->isExplicit(), "ODE operator is expressed in EXPLICIT form")
 }

 void ARKStepSolver::UseSundialsMassLinearSolver(int tdep)
 {
-   // Check that the ODE is expressed in mass form
-   MFEM_VERIFY(f_arkstep->ARKInMassForm(), "ODE operator is not in mass form.")
-
   // Free any existing matrix and linear solver
   if (M != NULL)   { SUNMatDestroy(A); M = NULL; }
   if (LSM != NULL) { SUNLinSolFree(LSM); LSM = NULL; }
@@ -1831,6 +1764,9 @@ void ARKStepSolver::UseSundialsMassLinearSolver(int tdep)
                                    ARKStepSolver::MassMult2, this);
   MFEM_VERIFY(flag == ARK_SUCCESS,
               "error in " STR(MFEM_ARKode(SetMassTimes)) "()");
+
+   // Check that the ODE is not expressed in EXPLICIT form
+   MFEM_VERIFY(!f->isExplicit(), "ODE operator is expressed in EXPLICIT form")
 }

 void ARKStepSolver::SetStepMode(int itask)
@@ -706,130 +706,9 @@ public:
 // Interface to ARKode's ARKStep module -- Additive Runge-Kutta methods
 // ---------------------------------------------------------------------------

-// Interface for defining ODE systems to be evolved using ARKStepSolver:
-//
-//  1)   du/dt = inv(M) f(u,t) ("MFEM" form)
-//  2) M dy/dt = f(u,t)        ("mass" form)
-//
-// where f(u,t) might be additively split, i.e., f(u,t) = f1(u,t) + f2(u,t)
-class ARKStepODE
-{
-
-public:
-
-   // the size of the ODE system
-   virtual int ARKSize() const = 0;
-
-   // return if the ODE system is of the form M du/dt = f(u,t), note the MFEM
-   // default is to use the form du/dt = int(M) f(u,t)
-   virtual bool ARKInMassForm() const { return false; };
-
-   // these flags are used by ARKStepSolver for switching between RK and ARK methods
-   enum ARKEvalMode
-   {  NORMAL,          // evaluate f(u,t)
-      ADDITIVE_TERM_1, // evaluate f1(u,t)
-      ADDITIVE_TERM_2  // evaluate f2(u,t)
-   };
-   virtual void ARKSetEvalMode(const ARKEvalMode new_eval_mode) {}
-
-   // evaluate either f(u,t) (mass form) or inv(M(t)) f(u,t) (MFEM form),
-   // which is necessary for solving ODEs with ERK or IMEX
-   virtual void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const
-   {
-      mfem_error("This function must be specified for ERK or IMEX methods.");
-   }
-
-  /** setup linear system for solving [M(t) - gamma Jf(u)] dk = f(u) - M(t) k,
-      which is necessary for solving ODEs with DIRK or IMEX methods
-      @param[in]  u     The state at which A(@a u,t) should be evaluated.
-      @param[in]  t     The time at which A(u,@a t) should be evaluated.
-      @param[in]  v     The value of inv(M) f(u,t) or f(u,t) for depending on form.
-      @param[in]  jok   Flag indicating if the Jacobian should be updated.
-      @param[out] jcur  Flag to signal if the Jacobian was updated.
-      @param[in]  gamma The scaled time step value.    */
-   virtual int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &v,
-                                int jok, int *jcur, real_t gamma)
-   {
-      mfem_error("This function must be specified for DIRK or IMEX methods.");
-   }
-
-  /** solve for dk in [M - gamma Jf(u)] dk = r, where r is either
-        inv(M) f(u,t) - k       (MFEM form)
-        f(u,t) - M k f(u) - M k (mass form)
-      when using DIRK or IMEX methods
-      @param[in]      r   inv(M) f(u,t) - k or f(u,t) - M k, depending on form.
-      @param[in,out]  dk  On input, the initial guess. On output, the solution.
-      @param[in]      tol Linear solve tolerance. */
-   virtual int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol)
-   {
-      mfem_error("This function must be specified for DIRK or IMEX methods.");
-   }
-
-  /** for mass form ODEs using an MFEM mass solver, setup the mass linear
-      system M(t) x = b
-      @param[in]  t  The time at which M(@a t) should be evaluated. */
-   virtual int ARKMassSetup(const real_t t)
-   {
-      mfem_error("This function must be specified to use MFEM mass solvers for mass form ODEs.");
-   }
-
-  /** for mass form ODEs using an MFEM mass solver, solve for x in M(t) x = b
-      @param[in]      b   The linear system right-hand side.
-      @param[in,out]  x   On input, the initial guess. On output, the solution.
-      @param[in]      tol Linear solve tolerance. */
-   virtual int ARKMassSolve(const Vector &b, Vector &x, real_t tol)
-   {
-      mfem_error("This function must be specified to use MFEM mass solver for mass form ODEs.");
-   }
-
-  /**  for mass form ODEs using an MFEM mass solver, evaluate M(t) x
-       @param[in]   x The vector to multiply.
-       @param[out]  v The result of the matrix-vector product. */
-   virtual int ARKMassMult(const Vector &x, Vector &v)
-   {
-      mfem_error("This function must be specified to use MFEM mass solver for mass form ODEs.");
-   }
-
-  /**  for mass form ODEs using a SUNDIALS mass solver, evaluate M(t) x
-       @param[in]   x The vector to multiply.
-       @param[in]   t The time at which M(@a t) should be evaluated.
-       @param[out]  v The result of the matrix-vector product. */
-   virtual int ARKMassMult(const Vector &x, const real_t t, Vector &v)
-   {
-      mfem_error("This function must be specified to use SUNDIALS mass solver for mass form ODEs.");
-   }
-
-};
-
 /// Interface to ARKode's ARKStep module -- additive Runge-Kutta methods.
 class ARKStepSolver : public ODESolver, public SundialsSolver
 {
-
-   // Wrapper class to provide backwards compatability with user code that
-   // derives from TimeDependentOperator instead of ARKStepODE
-   class TimeDependentOperatorWrapper : public ARKStepODE
-   {
-      TimeDependentOperator *tdo;
-
-   public:
-
-      TimeDependentOperatorWrapper(TimeDependentOperator *f);
-
-      int ARKSize() const override;
-      bool ARKInMassForm() const override;
-      void ARKSetEvalMode(const ARKEvalMode new_eval_mode) override;
-      void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const override;
-      int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &v,
-         int jok, int *jcur, real_t gamma) override;
-      int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
-      int ARKMassSetup(const real_t t) override;
-      int ARKMassSolve(const Vector &b, Vector &x, real_t tol) override;
-      int ARKMassMult(const Vector &x, Vector &v) override;
-      int ARKMassMult(const Vector &x, const real_t t, Vector &v) override;
-
-   };
-
-
 public:
   /// Types of ARKODE solvers.
   enum Type
@@ -843,8 +722,6 @@ protected:
   Type rk_type;      ///< Runge-Kutta type.
   int step_mode;     ///< ARKStep step mode (ARK_NORMAL or ARK_ONE_STEP).
   bool use_implicit; ///< True for implicit or imex integration.
-   ARKStepODE* f_arkstep;
-   std::unique_ptr<TimeDependentOperatorWrapper> f_tdo; // for backwards compatibility

   /** @name Wrappers to compute the ODE RHS functions.
       RHS1 is explicit RHS and RHS2 the implicit RHS for IMEX integration. When
@@ -907,19 +784,14 @@ public:
       then ARKStepReInit() will be called in the next call to Step(). If the
       problem size has changed, the ARKStep memory is freed and realloced
       for the new problem size. */
-   /** @param[in] f_ The ARKStepODE that defines the ODE system
+   /** @param[in] f_ The TimeDependentOperator that defines the ODE system

       @note All other methods must be called after Init().

       @note If this method is called a second time with a different problem
       size, then any non-default user-set options will be lost and will need
       to be set again. */
-   void Init(ARKStepODE *f_ark_);
-
-   // This method is provided for backwards compatibility with classes that
-   // derive TimeDependentOperator instead of ARKStepODE; however, those classes
-   // should be migrated.
-   MFEM_DEPRECATED void Init(TimeDependentOperator &f_) override;
+   void Init(TimeDependentOperator &f_) override;

   /// Integrate the ODE with ARKode using the specified step mode.
   /**
@@ -461,7 +461,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const

   Kmat.Mult(u, z);
   z.Neg(); // z = -z
-   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
+   K->EliminateVDofsInRHS(ess_tdof_list, u, z);

   M_solver.Mult(z, du_dt);
   du_dt.Print();
@@ -483,7 +483,7 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
   Kmat.Mult(u, z);
   z.Neg();
-   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
+   K->EliminateVDofsInRHS(ess_tdof_list, u, z);

   T_solver.Mult(z, du_dt);
   du_dt.SetSubVector(ess_tdof_list, 0.0);
@@ -89,7 +89,7 @@ The currently available test scripts are:
 - gitignore: Checks for compliance with MFEM's .gitignore rules

 - branch-history: Checks if the current branch history (the commits that will be
-  merged in master) contains unusually large files, unusually large number of
+  merged in main) contains unusually large files, unusually large number of
  changes in a commit, unusually large number of commits, etc.

 Adding new test scripts is easy and we welcome contributions from users and
@@ -27,7 +27,7 @@ sub usage {
          -h|--help   prints this usage information and exits

   This script checks if the current branch history, defined as the commits that
-   will be merged in master (those shown in a GitHub PR), contains unusually
+   will be merged in main (those shown in a GitHub PR), contains unusually
   large files, unusually large number of changes in a commit, unusually large
   number of commits, etc.

@@ -72,7 +72,7 @@ my $max_branch_kb = 1000;
 my $status = 0; # Return code

 # Get SHA hash of all commits in this branch
-my @commits = split /\n/, `git log --pretty=format:%H master..$branch`;
+my @commits = split /\n/, `git log --pretty=format:%H main..$branch`;

 # Check if total number of commits in this branch exceeds the maximum allowable
 my $ncommits = scalar @commits;
Author	SHA1	Message	Date
Tzanio Kolev	ae695de95a	Merge branch 'master' into main-dev	2025-10-16 07:06:45 -07:00
Veselin Dobrev	5436debf31	In config/githooks/pre-push, fetch 'main' instead of 'master'	2022-12-29 17:14:55 -08:00
Tzanio Kolev	2c875db251	Renamed the default branch from master to main	2022-12-29 12:03:23 -08:00