Merge branch 'master' into main-dev

.codecov.yml

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

CHANGELOG

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

CMakeLists.txt

@@ -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}

CONTRIBUTING.md

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

config/docker/Dockerfile.base

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

config/docker/spack.yaml

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

config/githooks/pre-push

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

examples/ex_sconv_p.cpp

@@ -1,503 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-static double freq = 0.5, kappa;
-static int dim;
-
-double u_func(const Vector &);
-
-enum SCA_TYPE {INVALID_SCA_TYPE = -1,
-               H1_TYPE = 0,
-               L2_TYPE,
-               L2I_TYPE,
-               NUM_SCA_TYPES
-              };
-enum CONV_TYPE {INVALID_CONV_TYPE = -1,
-                PROJECTION = 0,
-                INTERPOLATION_OP,
-                SOLVE,
-                SOLVE_W_DBC,
-                NUM_CONV_TYPES
-               };
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
-                                   FiniteElementCollection &fec);
-string GetTypeName(SCA_TYPE type);
-string GetConvTypeName(CONV_TYPE type);
-string GetConvTypeShortName(CONV_TYPE type);
-
-void Projection(const ParGridFunction &v0, ParGridFunction &v1);
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1);
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c);
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   MPI_Session mpi(argc, argv);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ser_ref_levels = 0;
-   int par_ref_levels = 0;
-   int order0 = 1;
-   int order1 = 1;
-   int type0 = 0;
-   int type1 = 1;
-   int conv_type = -1;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order0, "-o0", "--initial-order",
-                  "Finite element order (polynomial degree) "
-                  "for initial field.");
-   args.AddOption(&order1, "-o1", "--final-order",
-                  "Finite element order (polynomial degree) "
-                  "for final field.");
-   args.AddOption(&type0, "-t0", "--initial-type",
-                  "Set the basis type for the initial field: "
-                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
-   args.AddOption(&type1, "-t1", "--final-type",
-                  "Set the basis type for the final field: "
-                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
-   args.AddOption(&conv_type, "-c", "--conversion-type",
-                  "Set the conversion scheme: "
-                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
-                  "3-Least Squares with BC, -1 loop over all.");
-   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
-                  " solution.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-
-   args.Parse();
-   if (!args.Good())
-   {
-      if (mpi.Root()) { args.PrintUsage(cout); }
-      return 1;
-   }
-   if (mpi.Root()) { args.PrintOptions(cout); }
-   kappa = freq * M_PI;
-
-   // 3. Enable hardware devices such as GPUs, and programming models such as
-   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
-   Device device(device_config);
-   if (mpi.Root()) { device.Print(); }
-
-   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   dim = mesh->Dimension();
-
-   // 5. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
-   //    specified on the command line with -rs).
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution (1 time by
-   //    default, or specified on the command line with -rp). Once the parallel
-   //    mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh.UniformRefinement();
-   }
-
-   FunctionCoefficient uCoef(u_func);
-
-   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
-
-   if (mpi.Root())
-   {
-      cout << "L2 Errors:" << endl;
-   }
-   int t0a = (type0 == -1) ? 0 : type0;
-   int t0b = (type0 == -1) ? NUM_SCA_TYPES : (type0+1);
-   for (int t0 = t0a; t0 < t0b; t0++)
-   {
-      FiniteElementCollection *fec0 = GetFECollection((SCA_TYPE)t0, order0);
-      ParFiniteElementSpace   *fes0 = GetFESpace((SCA_TYPE)t0, pmesh, *fec0);
-
-      ParGridFunction x0(fes0);
-      x0.ProjectCoefficient(uCoef);
-
-      double err0 = x0.ComputeL2Error(uCoef);
-      if (mpi.Root())
-      {
-         cout << "Initial " << GetTypeName((SCA_TYPE)t0)
-              << ": \t\t" << err0 << endl;
-      }
-
-      // nn. Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")";
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock0(vishost, visport);
-         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
-                   << pmesh.GetMyRank() << '\n';
-         sol_sock0.precision(8);
-         sol_sock0 << "solution\n" << pmesh << x0
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
-                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
-                   << flush;
-      }
-
-      int t1a = (type1 == -1) ? 0 : type1;
-      int t1b = (type1 == -1) ? NUM_SCA_TYPES : (type1+1);
-      for (int t1 = t1a; t1 < t1b; t1++)
-      {
-         FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, order1);
-         ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, pmesh, *fec1);
-
-         ParGridFunction y1(fes1);
-
-         if (mpi.Root())
-         {
-            cout << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")"
-                 << " -> "
-                 << GetTypeName((SCA_TYPE)t1) << "(" << order1 << ")"
-                 << ":" << endl;
-         }
-
-         int c01a = (conv_type == -1) ? 0 : conv_type;
-         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
-         for (int c01 = c01a; c01 < c01b; c01++)
-         {
-            string cmnt = "";
-
-            switch ((CONV_TYPE)c01)
-            {
-               case PROJECTION:
-                  Projection(x0, y1);
-                  break;
-               case INTERPOLATION_OP:
-                  cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
-                         "(should match projection)" : "(not expected to succeed)";
-                  InterpolationOp(x0, y1);
-                  break;
-               case SOLVE:
-                  LeastSquares((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1);
-                  break;
-               case SOLVE_W_DBC:
-                  LeastSquaresBC((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1, uCoef);
-                  break;
-               default:
-                  y1 = 0.0;
-            }
-
-            double err1 = y1.ComputeL2Error(uCoef);
-            cout << GetConvTypeName((CONV_TYPE)c01)
-                 << "\t\t" << err1 << "\t" << cmnt << endl;
-
-            if (visualization)
-            {
-               ostringstream oss;
-               oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")" << " --"
-                   << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
-                   << GetTypeName((SCA_TYPE)t1)<< "(" << order1 << ")";
-               char vishost[] = "localhost";
-               int  visport   = 19916;
-               socketstream sol_sock1(vishost, visport);
-               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
-                         << pmesh.GetMyRank() << '\n';
-               sol_sock1.precision(8);
-               sol_sock1 << "solution\n" << pmesh << y1
-                         << "window_title '" << oss.str() << "'"
-                         << "window_geometry "
-                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
-                                  Ws * (t0 - t0a))
-                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
-                         << " " << Ww << " " << Wh
-                         << flush;
-            }
-         }
-         if (mpi.Root())
-         {
-            cout << endl;
-         }
-
-         delete fes1;
-         delete fec1;
-      }
-
-      delete fes0;
-      delete fec0;
-
-      if (t0 < t0b - 1)
-      {
-         char c;
-         if (mpi.Root())
-         {
-            cout << "press (q)uit or (c)ontinue --> " << flush;
-            cin >> c;
-         }
-         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
-
-         if (c != 'c')
-         {
-            break;
-         }
-      }
-      if (mpi.Root())
-      {
-         cout << endl;
-      }
-   }
-
-   return 0;
-}
-
-double u_func(const Vector &x)
-{
-   double kx = kappa * x[0];
-   double ky = kappa * x[1];
-   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
-
-   // Add the gradient of a scalar function
-   return cos(kx) * cos(ky) * cos(kz);
-}
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return new H1_FECollection(p, dim);
-      case L2_TYPE:
-         return new L2_FECollection(p-1, dim);
-      case L2I_TYPE:
-         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
-                                    FiniteElement::INTEGRAL);
-      default:
-         return NULL;
-   }
-}
-
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
-                                   ParMesh &pmesh,
-                                   FiniteElementCollection &fec)
-{
-   return new ParFiniteElementSpace(&pmesh, &fec);
-}
-
-string GetTypeName(SCA_TYPE type)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return "    H1";
-      case L2_TYPE:
-         return "    L2";
-      case L2I_TYPE:
-         return "    L2I";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Projection            ";
-      case INTERPOLATION_OP:
-         return "Interpolation Operator";
-      case SOLVE:
-         return "Least Squares         ";
-      case SOLVE_W_DBC:
-         return "Least Squares with BC ";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeShortName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Proj";
-      case INTERPOLATION_OP:
-         return "Interp";
-      case SOLVE:
-         return "LS";
-      case SOLVE_W_DBC:
-         return "LSwBC";
-      default:
-         return "--";
-   }
-}
-
-/** Perform a naive projection from one scalar field to another.
-
-   This scheme simply evaluates v0 at the interpolation points of v1.
-
-   If v0 has reduced continuity compared to v1 this can produce
-   results that depend on the order in which the elements are
-   traversed.
-
-   Suitable conversions:
-   H1 -> L2
-   H1 -> DG (same as L2)
- */
-void Projection(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   GridFunctionCoefficient v0Coef(&v0);
-   v1.ProjectCoefficient(v0Coef);
-}
-
-/** In theory this interpolation scheme should be equivalent to projection.
-
-    Building an interpolastion matrix could lead to computational
-    efficiency compared to simple projection if the operator will be
-    used several times.
-
-    Unfortunately this is broken for several combinations of source
-    and target fields.
-*/
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
-   op.AddDomainInterpolator(new IdentityInterpolator);
-   op.Assemble();
-   op.Finalize();
-
-   op.Mult(v0, v1);
-}
-
-/** Compute a least-squares best fit using the target basis functions.
-
-    This scheme is more difficult to setup and more computationally
-    expensive but the results can be significantly better than simple
-    projections.
-*/
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   HypreParMatrix * M = m.ParallelAssemble();
-
-   HypreDiagScale diag(*M);
-   HyprePCG pcg(*M);
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   Vector B, X;
-   b.ParallelAssemble(B);
-
-   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-
-   delete M;
-}
-
-/** Compute a least-squares best fit with boundary conditions.
-
-    This scheme is virtually identical to the previous one but it
-    makes use of boundary values, when available, to improve the
-    accuracy.  This scheme can produce significantly better results
-    when the normal derivative of the field is large near the
-    boundary.  This is particularly true when the field is
-    under-resolved near the boundary.
-*/
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
-   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
-      ess_bdr = 1;
-      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   if (t1 == H1_TYPE)
-   {
-      v1.ProjectBdrCoefficient(c, ess_bdr);
-   }
-
-   OperatorPtr M;
-   Vector B, X;
-   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
-
-   HypreDiagScale diag(*M.As<HypreParMatrix>());
-   HyprePCG pcg(*M.As<HypreParMatrix>());
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-}

examples/ex_vconv_p.cpp

@@ -1,586 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-static double freq = 0.5, kappa;
-static int dim;
-
-void u_func(const Vector &, Vector &);
-
-enum VEC_TYPE {INVALID_VEC_TYPE = -1,
-               H1V_TYPE = 0,
-               ND_TYPE,
-               RT_TYPE,
-               L2V_TYPE,
-               NUM_VEC_TYPES
-              };
-enum CONV_TYPE {INVALID_CONV_TYPE = -1,
-                PROJECTION = 0,
-                INTERPOLATION_OP,
-                SOLVE,
-                SOLVE_W_DBC,
-                NUM_CONV_TYPES
-               };
-
-FiniteElementCollection * GetFECollection(VEC_TYPE type, int p);
-ParFiniteElementSpace * GetFESpace(VEC_TYPE type, ParMesh &pmesh,
-                                   FiniteElementCollection &fec);
-string GetTypeName(VEC_TYPE type);
-string GetConvTypeName(CONV_TYPE type);
-
-void Projection(const ParGridFunction &v0, ParGridFunction &v1);
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
-void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
-                  VEC_TYPE t1, ParGridFunction &v1);
-void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
-                    VEC_TYPE t1, ParGridFunction &v1,
-                    VectorCoefficient &vc);
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   MPI_Session mpi(argc, argv);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ser_ref_levels = 0;
-   int par_ref_levels = 0;
-   int order0 = 1;
-   int order1 = 1;
-   int type0 = 0;
-   int type1 = 1;
-   int conv_type = -1;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order0, "-o0", "--initial-order",
-                  "Finite element order (polynomial degree) "
-                  "for initial field.");
-   args.AddOption(&order1, "-o1", "--final-order",
-                  "Finite element order (polynomial degree) "
-                  "for final field.");
-   args.AddOption(&type0, "-t0", "--initial-type",
-                  "Set the basis type for the initial field: "
-                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
-   args.AddOption(&type1, "-t1", "--final-type",
-                  "Set the basis type for the final field: "
-                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
-   args.AddOption(&conv_type, "-c", "--conversion-type",
-                  "Set the conversion scheme: "
-                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
-                  "3-Least Squares with BC, -1 loop over all.");
-   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
-                  " solution.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-
-   args.Parse();
-   if (!args.Good())
-   {
-      if (mpi.Root()) { args.PrintUsage(cout); }
-      return 1;
-   }
-   if (mpi.Root()) { args.PrintOptions(cout); }
-   kappa = freq * M_PI;
-
-   // 3. Enable hardware devices such as GPUs, and programming models such as
-   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
-   Device device(device_config);
-   if (mpi.Root()) { device.Print(); }
-
-   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   dim = mesh->Dimension();
-
-   // 5. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
-   //    specified on the command line with -rs).
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution (1 time by
-   //    default, or specified on the command line with -rp). Once the parallel
-   //    mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh.UniformRefinement();
-   }
-
-   VectorFunctionCoefficient uCoef(dim, u_func);
-
-   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
-
-   if (mpi.Root())
-   {
-      cout << "L2 Errors:" << endl;
-   }
-   int t0a = (type0 == -1) ? 0 : type0;
-   int t0b = (type0 == -1) ? NUM_VEC_TYPES : (type0+1);
-   for (int t0 = t0a; t0 < t0b; t0++)
-   {
-      FiniteElementCollection *fec0 = GetFECollection((VEC_TYPE)t0, order0);
-      ParFiniteElementSpace   *fes0 = GetFESpace((VEC_TYPE)t0, pmesh, *fec0);
-
-      ParGridFunction x0(fes0);
-      x0.ProjectCoefficient(uCoef);
-
-      double err0 = x0.ComputeL2Error(uCoef);
-      if (mpi.Root())
-      {
-         cout << "Initial " << GetTypeName((VEC_TYPE)t0)
-              << ": \t\t" << err0 << endl;
-      }
-
-      // nn. Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << GetTypeName((VEC_TYPE)t0);
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock0(vishost, visport);
-         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
-                   << pmesh.GetMyRank() << '\n';
-         sol_sock0.precision(8);
-         sol_sock0 << "solution\n" << pmesh << x0
-                   << "keys vvv "
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
-                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
-                   << flush;
-      }
-
-      int t1a = (type1 == -1) ? 0 : type1;
-      int t1b = (type1 == -1) ? NUM_VEC_TYPES : (type1+1);
-      for (int t1 = t1a; t1 < t1b; t1++)
-      {
-         FiniteElementCollection *fec1 = GetFECollection((VEC_TYPE)t1, order1);
-         ParFiniteElementSpace   *fes1 = GetFESpace((VEC_TYPE)t1, pmesh, *fec1);
-
-         ParGridFunction x1(fes1);
-
-         if (mpi.Root())
-         {
-            cout << GetTypeName((VEC_TYPE)t0) << " -> "
-                 << GetTypeName((VEC_TYPE)t1) << ":" << endl;
-         }
-
-         int c01a = (conv_type == -1) ? 0 : conv_type;
-         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
-         for (int c01 = c01a; c01 < c01b; c01++)
-         {
-            switch ((CONV_TYPE)c01)
-            {
-               case PROJECTION:
-                  Projection(x0, x1);
-                  break;
-               case INTERPOLATION_OP:
-                  // InterpolationOp(x0, x1);
-                  x1 = 0.0;
-                  break;
-               case SOLVE:
-                  LeastSquares((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1);
-                  break;
-               case SOLVE_W_DBC:
-                  LeastSquaresBC((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1, uCoef);
-                  break;
-               default:
-                  x1 = 0.0;
-            }
-
-            double err1 = x1.ComputeL2Error(uCoef);
-            cout << GetConvTypeName((CONV_TYPE)c01)
-                 << "\t\t" << err1 << endl;
-
-            if (visualization)
-            {
-               ostringstream oss;
-               oss << GetTypeName((VEC_TYPE)t0) << " --" << c01 << "--> "
-                   << GetTypeName((VEC_TYPE)t1);
-               char vishost[] = "localhost";
-               int  visport   = 19916;
-               socketstream sol_sock1(vishost, visport);
-               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
-                         << pmesh.GetMyRank() << '\n';
-               sol_sock1.precision(8);
-               sol_sock1 << "solution\n" << pmesh << x1
-                         << "keys vvv "
-                         << "window_title '" << oss.str() << "'"
-                         << "window_geometry "
-                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
-                                  Ws * (t0 - t0a))
-                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
-                         << " " << Ww << " " << Wh
-                         << flush;
-            }
-         }
-         if (mpi.Root())
-         {
-            cout << endl;
-         }
-
-         delete fes1;
-         delete fec1;
-      }
-
-      delete fes0;
-      delete fec0;
-
-      if (t0 < t0b - 1)
-      {
-         char c;
-         if (mpi.Root())
-         {
-            cout << "press (q)uit or (c)ontinue --> " << flush;
-            cin >> c;
-         }
-         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
-
-         if (c != 'c')
-         {
-            break;
-         }
-      }
-      if (mpi.Root())
-      {
-         cout << endl;
-      }
-   }
-
-   return 0;
-}
-
-void u_func(const Vector &x, Vector &u)
-{
-   u.SetSize(dim);
-
-   double kx = kappa * x[0];
-   double ky = kappa * x[1];
-   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
-
-   // Add the gradient of a scalar function
-   u(0) = sin(kx) * cos(ky);
-   u(1) = cos(kx) * sin(ky);
-   if (dim == 3)
-   {
-      u(0) *= cos(kz);
-      u(1) *= cos(kz);
-      u(2) = cos(kx) * cos(ky) * sin(kz);
-   }
-
-   // Add the curl of a vector function
-   u(0) -= cos(kx) * sin(ky);
-   u(1) += sin(kx) * cos(ky);
-   if (dim == 3)
-   {
-      u(0) += cos(kx) * sin(kz);
-      u(1) -= cos(ky) * sin(kz);
-      u(2) += (sin(ky) - sin(kx)) * cos(kz);
-   }
-}
-
-FiniteElementCollection * GetFECollection(VEC_TYPE type, int p)
-{
-   switch (type)
-   {
-      case H1V_TYPE:
-         return new H1_FECollection(p, dim);
-      case ND_TYPE:
-         return new ND_FECollection(p, dim);
-      case RT_TYPE:
-         return new RT_FECollection(p-1, dim);
-      case L2V_TYPE:
-         return new L2_FECollection(p-1, dim);
-      default:
-         return NULL;
-   }
-}
-
-ParFiniteElementSpace * GetFESpace(VEC_TYPE type,
-                                   ParMesh &pmesh,
-                                   FiniteElementCollection &fec)
-{
-   switch (type)
-   {
-      case H1V_TYPE:
-      case L2V_TYPE:
-         return new ParFiniteElementSpace(&pmesh, &fec, dim);
-      case ND_TYPE:
-      case RT_TYPE:
-         return new ParFiniteElementSpace(&pmesh, &fec);
-      default:
-         return NULL;
-   }
-}
-
-string GetTypeName(VEC_TYPE type)
-{
-   switch (type)
-   {
-      case H1V_TYPE:
-         return "    H1V";
-      case ND_TYPE:
-         return "H(Curl)";
-      case RT_TYPE:
-         return " H(Div)";
-      case L2V_TYPE:
-         return "    L2V";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Projection            ";
-      case INTERPOLATION_OP:
-         return "Interpolation Operator";
-      case SOLVE:
-         return "Least Squares         ";
-      case SOLVE_W_DBC:
-         return "Least Squares with BC ";
-      default:
-         return "--";
-   }
-}
-
-/** Perform a naive projection from one vector field to another.
-
-   This scheme simply evaluates v0 at the interpolation points of v1.
-
-   If v0 has reduced continuity compared to v1 this can produce
-   results that depend on the order in which the elements are
-   traversed.
-
-   Suitable conversions:
-   H1V     -> H(Curl), H(Div), or L2V
-   H(Curl) -> L2V
-   H(Div)  -> L2V
-
- */
-void Projection(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   VectorGridFunctionCoefficient v0Coef(&v0);
-   v1.ProjectCoefficient(v0Coef);
-}
-
-/** In theory this interpolation scheme should be equivalent to projection.
-
-    Building an interpolastion matrix could lead to computational
-    efficiency compared to simple projection if the operator will be
-    used several times.
-
-    Unfortunately this is broken for several combinations of source
-    and target fields.
-*/
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
-   op.AddDomainInterpolator(new IdentityInterpolator);
-   op.Assemble();
-   op.Finalize();
-
-   op.Mult(v0, v1);
-}
-
-/** Compute a least-squares best fit using the target basis functions.
-
-    This scheme is more difficult to setup and more computationally
-    expensive but the results can be significantly better than simple
-    projections.
-*/
-void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
-                  VEC_TYPE t1, ParGridFunction &v1)
-{
-   bool trans = false;
-   ParFiniteElementSpace *fes0, *fes1;
-   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
-       (t1 == ND_TYPE || t1 == RT_TYPE))
-   {
-      fes0 = v1.ParFESpace();
-      fes1 = v0.ParFESpace();
-      trans = true;
-   }
-   else
-   {
-      fes0 = v0.ParFESpace();
-      fes1 = v1.ParFESpace();
-   }
-   ParMixedBilinearForm op(fes0, fes1);
-   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      op.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      op.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   if (trans)
-   {
-      op.MultTranspose(v0, b);
-   }
-   else
-   {
-      op.Mult(v0, b);
-   }
-
-   ParBilinearForm m(v1.ParFESpace());
-   if (t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      m.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      m.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   m.Assemble();
-   m.Finalize();
-
-   HypreParMatrix * M = m.ParallelAssemble();
-
-   HypreDiagScale diag(*M);
-   HyprePCG pcg(*M);
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   Vector B, X;
-   b.ParallelAssemble(B);
-
-   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-
-   delete M;
-}
-
-/** Compute a least-squares best fit with boundary conditions.
-
-    This scheme is virtually identical to the previous one but it
-    makes use of boundary values, when available, to improve the
-    accuracy.  This scheme can produce significantly better results
-    when the normal derivative of the field is large near the
-    boundary.  This is particularly true when the field is
-    under-resolved near the boundary.
-*/
-void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
-                    VEC_TYPE t1, ParGridFunction &v1,
-                    VectorCoefficient &vc)
-{
-   bool trans = false;
-   ParFiniteElementSpace *fes0, *fes1;
-   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
-       (t1 == ND_TYPE || t1 == RT_TYPE))
-   {
-      fes0 = v1.ParFESpace();
-      fes1 = v0.ParFESpace();
-      trans = true;
-   }
-   else
-   {
-      fes0 = v0.ParFESpace();
-      fes1 = v1.ParFESpace();
-   }
-   ParMixedBilinearForm op(fes0, fes1);
-   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      op.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      op.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   if (trans)
-   {
-      op.MultTranspose(v0, b);
-   }
-   else
-   {
-      op.Mult(v0, b);
-   }
-
-   ParBilinearForm m(v1.ParFESpace());
-   if (t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      m.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      m.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   m.Assemble();
-   m.Finalize();
-
-   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
-   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
-      ess_bdr = 1;
-      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   if (t1 == H1V_TYPE)
-   {
-      v1.ProjectBdrCoefficient(vc, ess_bdr);
-   }
-   if (t1 == ND_TYPE)
-   {
-      v1.ProjectBdrCoefficientTangent(vc, ess_bdr);
-   }
-   else if (t1 == RT_TYPE)
-   {
-      v1.ProjectBdrCoefficientNormal(vc, ess_bdr);
-   }
-
-   OperatorPtr M;
-   Vector B, X;
-   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
-
-   HypreDiagScale diag(*M.As<HypreParMatrix>());
-   HyprePCG pcg(*M.As<HypreParMatrix>());
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-}

examples/jupyter/README.md

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

fem/nonlinearform.cpp

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

fem/nonlinearform.hpp

@@ -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().

fem/pbilinearform.cpp

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

fem/pbilinearform.hpp

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

fem/pnonlinearform.cpp

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

fem/pnonlinearform.hpp

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

linalg/hypre.cpp

@@ -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();
 

linalg/hypre.hpp

@@ -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. */

miniapps/navier/navier_cht.cpp

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

miniapps/tools/gf-conv.cpp

@@ -1,577 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-static int dim;
-
-enum SCA_TYPE {INVALID_SCA_TYPE = -1,
-               H1_TYPE = 0,
-               L2_TYPE,
-               L2I_TYPE,
-               NUM_SCA_TYPES
-              };
-enum CONV_TYPE {INVALID_CONV_TYPE = -1,
-                PROJECTION = 0,
-                INTERPOLATION_OP,
-                SOLVE,
-                SOLVE_W_DBC,
-                NUM_CONV_TYPES
-               };
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
-                                   FiniteElementCollection &fec);
-void parseFieldNames(const char * field_name_c_str,
-                     vector<string> &field_names);
-
-string GetTypeName(SCA_TYPE type);
-string GetTypeShortName(SCA_TYPE type);
-string GetConvTypeName(CONV_TYPE type);
-string GetConvTypeShortName(CONV_TYPE type);
-
-void Projection(const ParGridFunction &v0, ParGridFunction &v1);
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1);
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c);
-
-int main(int argc, char *argv[])
-{
-#ifdef MFEM_USE_MPI
-   Mpi::Init();
-   if (!Mpi::Root()) { mfem::out.Disable(); mfem::err.Disable(); }
-   Hypre::Init();
-#endif
-
-   // Parse command-line options.
-   const char *coll_name = NULL;
-   int cycle = 0;
-
-   const char *field_name_c_str = "ALL";
-
-   Array<int> orders;
-   Array<int> types;
-   Array<int> conv_types;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&coll_name, "-r", "--root-file",
-                  "Set the VisIt data collection root file prefix.", true);
-   args.AddOption(&cycle, "-c", "--cycle", "Set the cycle index to read.");
-   args.AddOption(&field_name_c_str, "-fn", "--field-names",
-                  "List of field names to get values from.");
-   args.AddOption(&orders, "-o", "--final-order",
-                  "Finite element orders for each final field "
-                  "(an array of integers for multiple fields).");
-   args.AddOption(&types, "-t", "--final-type",
-                  "Set the basis type for the final fields: "
-                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
-   args.AddOption(&conv_types, "-ct", "--conversion-type",
-                  "Set the conversion schemes: "
-                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
-                  "3-Least Squares with BC, -1 loop over all.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(mfem::out);
-      return 1;
-   }
-   args.PrintOptions(mfem::out);
-
-#ifdef MFEM_USE_MPI
-   VisItDataCollection dc(MPI_COMM_WORLD, coll_name);
-#else
-   VisItDataCollection dc(coll_name);
-#endif
-
-   dc.Load(cycle);
-
-   if (dc.Error() != DataCollection::NO_ERROR)
-   {
-      mfem::out << "Error loading VisIt data collection: " << coll_name << endl;
-      return 1;
-   }
-
-   dim = dc.GetMesh()->Dimension();
-   int spaceDim = dc.GetMesh()->SpaceDimension();
-
-   mfem::out << endl;
-   mfem::out << "Collection Name:    " << dc.GetCollectionName() << endl;
-   mfem::out << "Manifold Dimension: " << dim << endl;
-   mfem::out << "Space Dimension:    " << spaceDim << endl;
-   mfem::out << "Cycle:              " << dc.GetCycle() << endl;
-   mfem::out << "Time:               " << dc.GetTime() << endl;
-   mfem::out << "Time Step:          " << dc.GetTimeStep() << endl;
-   mfem::out << endl;
-
-   typedef DataCollection::FieldMapType fields_t;
-   const fields_t &fields = dc.GetFieldMap();
-   // Print the names of all fields.
-   mfem::out << "fields: [ ";
-   for (fields_t::const_iterator it = fields.begin(); it != fields.end(); ++it)
-   {
-      if (it != fields.begin()) { mfem::out << ", "; }
-      mfem::out << it->first;
-   }
-   mfem::out << " ]" << endl;
-
-   // Parsing desired field names
-   vector<string> field_names;
-   parseFieldNames(field_name_c_str, field_names);
-
-   if (field_names.size() == 1)
-   {
-      if (field_names[0] == "ALL")
-      {
-         fields_t::const_iterator it = fields.begin();
-         field_names[0] = it->first; it++;
-         for ( ; it != fields.end(); ++it)
-         {
-            field_names.push_back(it->first);
-         }
-      }
-   }
-
-   if (orders.Size() < field_names.size())
-   {
-      int size = orders.Size();
-      int order = (size > 0) ? orders[0] : 1;
-
-      orders.SetSize(field_names.size());
-      for (int i=size; i < field_names.size(); i++)
-      {
-         orders[i] = order;
-      }
-   }
-
-   if (types.Size() < field_names.size())
-   {
-      int size = types.Size();
-      int type = (size > 0) ? types[0] : 0;
-
-      types.SetSize(field_names.size());
-      for (int i=size; i < field_names.size(); i++)
-      {
-         types[i] = type;
-      }
-   }
-
-   if (conv_types.Size() < field_names.size())
-   {
-      int size = conv_types.Size();
-      int type = (size > 0) ? conv_types[0] : 0;
-
-      conv_types.SetSize(field_names.size());
-      for (int i=size; i < field_names.size(); i++)
-      {
-         conv_types[i] = type;
-      }
-   }
-
-   // Print field names to be extracted
-   mfem::out << "Extracting fields: ";
-   for (int i=0; i < field_names.size(); i++)
-   {
-      mfem::out << " \"" << field_names[i] << "\"";
-   }
-   mfem::out << endl;
-
-#ifdef MFEM_USE_MPI
-   ParMesh *mesh = dynamic_cast<ParMesh*>(dc.GetMesh());
-#else
-   Mesh *mesh = dc.GetMesh();
-#endif
-   if (mesh == NULL)
-   {
-      mfem::out << "Problem with mesh\n";
-      return 1;
-   }
-
-   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
-
-   // Loop over all requested fields.
-   for (int i=0; i < field_names.size(); i++)
-   {
-#ifdef MFEM_USE_MPI
-      ParGridFunction *x0 = dc.GetParField(field_names[i]);
-#else
-      GridFunction *x0 = dc.GetField(field_names[i]);
-#endif
-      if (x0 == NULL)
-      {
-         mfem::out << "Problem with x0 for field \"" << field_names[i] << "\"\n";
-         continue;
-      }
-
-      int t0 = 0;
-
-      // nn. Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << field_names[i];
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock0(vishost, visport);
-#ifdef MFEM_USE_MPI
-         sol_sock0 << "parallel " << mesh->GetNRanks() << ' '
-                   << mesh->GetMyRank() << '\n';
-#endif
-         sol_sock0.precision(8);
-         sol_sock0 << "solution\n" << *mesh << *x0
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << Ws * (t0) << " " << Ws * (t0) << " "
-                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
-                   << flush;
-      }
-
-      int t1 = types[i];
-      FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, orders[i]);
-      ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, *mesh, *fec1);
-
-      ParGridFunction *y1 = new ParGridFunction(fes1);
-
-      mfem::out << GetTypeName((SCA_TYPE)t1) << "(" << orders[i] << ")"
-                << ":" << endl;
-
-      int c01 = conv_types[i];
-      string cmnt = "";
-
-      switch ((CONV_TYPE)c01)
-      {
-         case PROJECTION:
-            Projection(*x0, *y1);
-            break;
-         case INTERPOLATION_OP:
-            cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
-                   "(should match projection)" : "(not expected to succeed)";
-            InterpolationOp(*x0, *y1);
-            break;
-         case SOLVE:
-            LeastSquares((SCA_TYPE)t0, *x0, (SCA_TYPE)t1, *y1);
-            break;
-         default:
-            *y1 = 0.0;
-      }
-
-      {
-         ostringstream oss;
-         oss << field_names[i] << "_" << GetConvTypeShortName((CONV_TYPE)c01)
-             << "_" << GetTypeShortName((SCA_TYPE)t1) << "_o" << orders[i];
-
-         dc.RegisterField(oss.str(), y1);
-      }
-
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
-             << GetTypeName((SCA_TYPE)t1)<< "(" << orders[i] << ")";
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock1(vishost, visport);
-#ifdef MFEM_USE_MPI
-         sol_sock1 << "parallel " << mesh->GetNRanks() << ' '
-                   << mesh->GetMyRank() << '\n';
-#endif
-         sol_sock1.precision(8);
-         sol_sock1 << "solution\n" << *mesh << y1
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << (int)((Ww + Fw) * (1.5 + c01) +
-                            Ws * (t0))
-                   << " " << (Wh + Fh) * (t1) + Ws * (t0)
-                   << " " << Ww << " " << Wh
-                   << flush;
-      }
-      mfem::out << endl;
-
-      // delete fes1;
-      // delete fec1;
-   }
-   dc.Save();
-
-   return 0;
-}
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return new H1_FECollection(p, dim);
-      case L2_TYPE:
-         return new L2_FECollection(p-1, dim);
-      case L2I_TYPE:
-         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
-                                    FiniteElement::INTEGRAL);
-      default:
-         return NULL;
-   }
-}
-
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
-                                   ParMesh &pmesh,
-                                   FiniteElementCollection &fec)
-{
-   return new ParFiniteElementSpace(&pmesh, &fec);
-}
-
-string GetTypeName(SCA_TYPE type)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return "    H1";
-      case L2_TYPE:
-         return "    L2";
-      case L2I_TYPE:
-         return "    L2I";
-      default:
-         return "--";
-   }
-}
-
-string GetTypeShortName(SCA_TYPE type)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return "H1";
-      case L2_TYPE:
-         return "L2";
-      case L2I_TYPE:
-         return "L2I";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Projection            ";
-      case INTERPOLATION_OP:
-         return "Interpolation Operator";
-      case SOLVE:
-         return "Least Squares         ";
-      case SOLVE_W_DBC:
-         return "Least Squares with BC ";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeShortName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Proj";
-      case INTERPOLATION_OP:
-         return "Interp";
-      case SOLVE:
-         return "LS";
-      case SOLVE_W_DBC:
-         return "LSwBC";
-      default:
-         return "--";
-   }
-}
-
-void parseFieldNames(const char * field_name_c_str, vector<string> &field_names)
-{
-   string field_name_str(field_name_c_str);
-   string field_name;
-
-   for (string::iterator it=field_name_str.begin();
-        it!=field_name_str.end(); it++)
-   {
-      if (*it == '\\')
-      {
-         it++;
-         field_name.push_back(*it);
-      }
-      else if (*it == ' ')
-      {
-         if (!field_name.empty())
-         {
-            field_names.push_back(field_name);
-         }
-         field_name.clear();
-      }
-      else if (it == field_name_str.end() - 1)
-      {
-         field_name.push_back(*it);
-         field_names.push_back(field_name);
-      }
-      else
-      {
-         field_name.push_back(*it);
-      }
-   }
-   if (field_names.size() == 0)
-   {
-      field_names.push_back("ALL");
-   }
-}
-
-/** Perform a naive projection from one scalar field to another.
-
-   This scheme simply evaluates v0 at the interpolation points of v1.
-
-   If v0 has reduced continuity compared to v1 this can produce
-   results that depend on the order in which the elements are
-   traversed.
-
-   Suitable conversions:
-   H1 -> L2
-   H1 -> DG (same as L2)
- */
-void Projection(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   GridFunctionCoefficient v0Coef(&v0);
-   v1.ProjectCoefficient(v0Coef);
-}
-
-/** In theory this interpolation scheme should be equivalent to projection.
-
-    Building an interpolastion matrix could lead to computational
-    efficiency compared to simple projection if the operator will be
-    used several times.
-
-    Unfortunately this is broken for several combinations of source
-    and target fields.
-*/
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
-   op.AddDomainInterpolator(new IdentityInterpolator);
-   op.Assemble();
-   op.Finalize();
-
-   op.Mult(v0, v1);
-}
-
-/** Compute a least-squares best fit using the target basis functions.
-
-    This scheme is more difficult to setup and more computationally
-    expensive but the results can be significantly better than simple
-    projections.
-*/
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   HypreParMatrix * M = m.ParallelAssemble();
-
-   HypreDiagScale diag(*M);
-   HyprePCG pcg(*M);
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   Vector B, X;
-   b.ParallelAssemble(B);
-
-   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-
-   delete M;
-}
-
-/** Compute a least-squares best fit with boundary conditions.
-
-    This scheme is virtually identical to the previous one but it
-    makes use of boundary values, when available, to improve the
-    accuracy.  This scheme can produce significantly better results
-    when the normal derivative of the field is large near the
-    boundary.  This is particularly true when the field is
-    under-resolved near the boundary.
-*/
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
-   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
-      ess_bdr = 1;
-      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   if (t1 == H1_TYPE)
-   {
-      v1.ProjectBdrCoefficient(c, ess_bdr);
-   }
-
-   OperatorPtr M;
-   Vector B, X;
-   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
-
-   HypreDiagScale diag(*M.As<HypreParMatrix>());
-   HyprePCG pcg(*M.As<HypreParMatrix>());
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-}

tests/scripts/README

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

tests/scripts/branch-history

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