Extending a unit test to include pyramids

vsx128 GNU fix [vsx]

.gitlab/configs/lassen-config.yml

@@ -45,5 +45,5 @@ variables:
     - echo ${MFEM_DATA_DIR}
     - echo ${SPEC}
     # Next script uses 'THREADS': leaving it empty --> it uses 'make all -j'
-    - lalloc 1 -W 30 -q pdebug tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - lalloc 1 -W 30 -q pdebug --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
   needs: [setup]

CHANGELOG

@@ -10,10 +10,16 @@
 
 Version 4.3.1 (development)
 ===========================
-- Added support for automatic differentiation. Users can select between 
-  native implementation and external library implementation at the 
-  configuration phase. A parallel and two serial examples are implemented
-  in the autodiff miniapp directory.
+- More explicit and consistent formating of the output of iterative solvers
+  with the new IterativeSolver::PrintLevel options. See linalg/solvers.hpp.
+
+- Added support for automatic differentiation. Users can select between native
+  implementation and external library implementation during configuration. One
+  parallel and two serial examples are implemented in the miniapps/autodiff/
+  directory.
+
+- GridFunctionCoefficient (and the related vector, gradient, divergence, and
+  curl classes) now work properly with LORDiscretization and LORSolver.
 
 - Added support for mesh preprocessing to resolve fine scale problem data
   before simulation. This feature uses adaptive mesh refinement to control the
@@ -71,13 +77,15 @@ Version 4.3.1 (development)
 
 - Added initial TMOP-based capabilities for surface fitting and tangential
   relaxation in the mesh-optimizer and pmesh-optimizer miniapps.
- 
+
 - Added ParMesh Adjaceny Set (adjset) creation support to the Conduit Mesh
   Blueprint MFEM wrapper functions in ConduitDataCollection.
 
 - `HypreParVector` and `Vector` now support move semantics, and the copy
   constructor for `HypreParVector` now copies the local vector data.
 
+- The HPC versions of ex1 and ex1p (in miniapps/performance) now support
+  runtime selection of either 2D or 3D meshes.
 
 Version 4.3, released on July 29, 2021
 ======================================

CODE_OF_CONDUCT.md

@@ -0,0 +1,133 @@
+# MFEM Code of Conduct
+
+## Our Pledge
+
+We as members, contributors, and leaders pledge to make participation in our
+community a harassment-free experience for everyone, regardless of age, body
+size, visible or invisible disability, ethnicity, sex characteristics, gender
+identity and expression, level of experience, education, socio-economic status,
+nationality, personal appearance, race, caste, color, religion, or sexual
+identity and orientation.
+
+We pledge to act and interact in ways that contribute to an open, welcoming,
+diverse, inclusive, and healthy community.
+
+## Our Standards
+
+Examples of behavior that contributes to a positive environment for our
+community include:
+
+* Demonstrating empathy and kindness toward other people
+* Being respectful of differing opinions, viewpoints, and experiences
+* Giving and gracefully accepting constructive feedback
+* Accepting responsibility and apologizing to those affected by our mistakes,
+  and learning from the experience
+* Focusing on what is best not just for us as individuals, but for the overall
+  community
+
+Examples of unacceptable behavior include:
+
+* The use of sexualized language or imagery, and sexual attention or advances of
+  any kind
+* Trolling, insulting or derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or email address,
+  without their explicit permission
+* Other conduct which could reasonably be considered inappropriate in a
+  professional setting
+
+## Enforcement Responsibilities
+
+Community leaders are responsible for clarifying and enforcing our standards of
+acceptable behavior and will take appropriate and fair corrective action in
+response to any behavior that they deem inappropriate, threatening, offensive,
+or harmful.
+
+Community leaders have the right and responsibility to remove, edit, or reject
+comments, commits, code, wiki edits, issues, and other contributions that are
+not aligned to this Code of Conduct, and will communicate reasons for moderation
+decisions when appropriate.
+
+## Scope
+
+This Code of Conduct applies within all community spaces, and also applies when
+an individual is officially representing the community in public spaces.
+Examples of representing our community include using an official e-mail address,
+posting via an official social media account, or acting as an appointed
+representative at an online or offline event.
+
+## Enforcement
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be
+reported to the community leaders responsible for enforcement at mfem@llnl.gov.
+All complaints will be reviewed and investigated promptly and fairly.
+
+All community leaders are obligated to respect the privacy and security of the
+reporter of any incident. Anyone involved in the reported behavior will recuse
+themselves from the investigation and decision making about the resolution of
+the complaint.
+
+## Enforcement Guidelines
+
+Community leaders will follow these Community Impact Guidelines in determining
+the consequences for any action they deem in violation of this Code of Conduct:
+
+### 1. Correction
+
+**Community Impact**: Use of inappropriate language or other behavior deemed
+unprofessional or unwelcome in the community.
+
+**Consequence**: A private, written warning from community leaders, providing
+clarity around the nature of the violation and an explanation of why the
+behavior was inappropriate. A public apology may be requested.
+
+### 2. Warning
+
+**Community Impact**: A violation through a single incident or series of
+actions.
+
+**Consequence**: A warning with consequences for continued behavior. No
+interaction with the people involved, including unsolicited interaction with
+those enforcing the Code of Conduct, for a specified period of time. This
+includes avoiding interactions in community spaces as well as external channels
+like social media. Violating these terms may lead to a temporary or permanent
+ban.
+
+### 3. Temporary Ban
+
+**Community Impact**: A serious violation of community standards, including
+sustained inappropriate behavior.
+
+**Consequence**: A temporary ban from any sort of interaction or public
+communication with the community for a specified period of time. No public or
+private interaction with the people involved, including unsolicited interaction
+with those enforcing the Code of Conduct, is allowed during this period.
+Violating these terms may lead to a permanent ban.
+
+### 4. Permanent Ban
+
+**Community Impact**: Demonstrating a pattern of violation of community
+standards, including sustained inappropriate behavior, harassment of an
+individual, or aggression toward or disparagement of classes of individuals.
+
+**Consequence**: A permanent ban from any sort of public interaction within the
+community.
+
+## Attribution
+
+This Code of Conduct is adapted from the [Contributor Covenant][homepage],
+version 2.1, available at
+[https://www.contributor-covenant.org/version/2/1/code_of_conduct.html][v2.1].
+
+Community Impact Guidelines were inspired by
+[Mozilla's code of conduct enforcement ladder][Mozilla CoC].
+
+For answers to common questions about this code of conduct, see the FAQ at
+[https://www.contributor-covenant.org/faq][FAQ]. Translations are available at
+[https://www.contributor-covenant.org/translations][translations].
+
+[homepage]: https://www.contributor-covenant.org
+[v2.1]: https://www.contributor-covenant.org/version/2/1/code_of_conduct.html
+[Mozilla CoC]: https://github.com/mozilla/diversity
+[FAQ]: https://www.contributor-covenant.org/faq
+[translations]: https://www.contributor-covenant.org/translations

CONTRIBUTING.md

@@ -21,6 +21,9 @@ documentation; new examples and miniapps; HPC performance improvements; etc.
 MFEM is distributed under the terms of the BSD-3 license. All new contributions
 must be made under this license.
 
+Note also that MFEM has a [Code of Conduct](CODE_OF_CONDUCT.md). By participating
+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

README

@@ -12,6 +12,9 @@ to enable high-performance scalable finite element discretization research and
 application development on a wide variety of platforms, ranging from laptops to
 supercomputers.
 
+We welcome contributions and feedback from the community. Please see the file
+CONTRIBUTING.md for additional details about our development process.
+
 * For building instructions, see the file INSTALL, or type "make help".
 
 * Copyright and licensing information can be found in files LICENSE and NOTICE.
@@ -19,9 +22,6 @@ supercomputers.
 * The best starting point for new users interested in MFEM's features is to
   review the examples and miniapps at https://mfem.org/examples.
 
-* Developers interested in contributing to the library, should read the
-  instructions and documentation in the CONTRIBUTING.md file.
-
 Conceptually, MFEM can be viewed as a finite element toolbox that provides the
 building blocks for developing finite element algorithms in a manner similar to
 that of MATLAB for linear algebra methods. In particular, MFEM provides support

examples/ex10.cpp

@@ -149,7 +149,7 @@ void InitialDeformation(const Vector &x, Vector &y);
 
 void InitialVelocity(const Vector &x, Vector &v);
 
-void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
+void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
                GridFunction *field, const char *field_name = NULL,
                bool init_vis = false);
 
@@ -376,10 +376,10 @@ int main(int argc, char *argv[])
 }
 
 
-void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
+void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
                GridFunction *field, const char *field_name, bool init_vis)
 {
-   if (!out)
+   if (!os)
    {
       return;
    }
@@ -389,24 +389,25 @@ void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
 
    mesh->SwapNodes(nodes, owns_nodes);
 
-   out << "solution\n" << *mesh << *field;
+   os << "solution\n" << *mesh << *field;
 
    mesh->SwapNodes(nodes, owns_nodes);
 
    if (init_vis)
    {
-      out << "window_size 800 800\n";
-      out << "window_title '" << field_name << "'\n";
+      os << "window_size 800 800\n";
+      os << "window_title '" << field_name << "'\n";
       if (mesh->SpaceDimension() == 2)
       {
-         out << "view 0 0\n"; // view from top
-         out << "keys jl\n";  // turn off perspective and light
+         os << "view 0 0\n"; // view from top
+         os << "keys jl\n";  // turn off perspective and light
       }
-      out << "keys cm\n";         // show colorbar and mesh
-      out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
-      out << "pause\n";
+      os << "keys cm\n";         // show colorbar and mesh
+      // update value-range; keep mesh-extents fixed
+      os << "autoscale value\n";
+      os << "pause\n";
    }
-   out << flush;
+   os << flush;
 }
 
 

examples/ex10p.cpp

@@ -154,7 +154,8 @@ void InitialDeformation(const Vector &x, Vector &y);
 
 void InitialVelocity(const Vector &x, Vector &v);
 
-void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
+void visualize(ostream &os, ParMesh *mesh,
+               ParGridFunction *deformed_nodes,
                ParGridFunction *field, const char *field_name = NULL,
                bool init_vis = false);
 
@@ -438,10 +439,11 @@ int main(int argc, char *argv[])
    return 0;
 }
 
-void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
+void visualize(ostream &os, ParMesh *mesh,
+               ParGridFunction *deformed_nodes,
                ParGridFunction *field, const char *field_name, bool init_vis)
 {
-   if (!out)
+   if (!os)
    {
       return;
    }
@@ -451,25 +453,27 @@ void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
 
    mesh->SwapNodes(nodes, owns_nodes);
 
-   out << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
-   out << "solution\n" << *mesh << *field;
+   os << "parallel " << mesh->GetNRanks()
+      << " " << mesh->GetMyRank() << "\n";
+   os << "solution\n" << *mesh << *field;
 
    mesh->SwapNodes(nodes, owns_nodes);
 
    if (init_vis)
    {
-      out << "window_size 800 800\n";
-      out << "window_title '" << field_name << "'\n";
+      os << "window_size 800 800\n";
+      os << "window_title '" << field_name << "'\n";
       if (mesh->SpaceDimension() == 2)
       {
-         out << "view 0 0\n"; // view from top
-         out << "keys jl\n";  // turn off perspective and light
+         os << "view 0 0\n"; // view from top
+         os << "keys jl\n";  // turn off perspective and light
       }
-      out << "keys cm\n";         // show colorbar and mesh
-      out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
-      out << "pause\n";
+      os << "keys cm\n";         // show colorbar and mesh
+      // update value-range; keep mesh-extents fixed
+      os << "autoscale value\n";
+      os << "pause\n";
    }
-   out << flush;
+   os << flush;
 }
 
 

examples/ex18.hpp

@@ -32,7 +32,7 @@ private:
    mutable DenseTensor flux;
    mutable Vector z;
 
-   void GetFlux(const DenseMatrix &state, DenseTensor &flux) const;
+   void GetFlux(const DenseMatrix &state_, DenseTensor &flux_) const;
 
 public:
    FE_Evolution(FiniteElementSpace &vfes_,
@@ -256,26 +256,26 @@ inline double ComputeMaxCharSpeed(const Vector &state, const int dim)
 }
 
 // Compute the flux at solution nodes.
-void FE_Evolution::GetFlux(const DenseMatrix &x, DenseTensor &flux) const
+void FE_Evolution::GetFlux(const DenseMatrix &x_, DenseTensor &flux_) const
 {
-   const int dof = flux.SizeI();
-   const int dim = flux.SizeJ();
+   const int flux_dof = flux_.SizeI();
+   const int flux_dim = flux_.SizeJ();
 
-   for (int i = 0; i < dof; i++)
+   for (int i = 0; i < flux_dof; i++)
    {
-      for (int k = 0; k < num_equation; k++) { state(k) = x(i, k); }
-      ComputeFlux(state, dim, f);
+      for (int k = 0; k < num_equation; k++) { state(k) = x_(i, k); }
+      ComputeFlux(state, flux_dim, f);
 
-      for (int d = 0; d < dim; d++)
+      for (int d = 0; d < flux_dim; d++)
       {
          for (int k = 0; k < num_equation; k++)
          {
-            flux(i, d, k) = f(k, d);
+            flux_(i, d, k) = f(k, d);
          }
       }
 
       // Update max char speed
-      const double mcs = ComputeMaxCharSpeed(state, dim);
+      const double mcs = ComputeMaxCharSpeed(state, flux_dim);
       if (mcs > max_char_speed) { max_char_speed = mcs; }
    }
 }

examples/ex19.cpp

@@ -171,7 +171,7 @@ public:
 };
 
 // Visualization driver
-void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
+void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
                GridFunction *field, const char *field_name = NULL,
                bool init_vis = false);
 
@@ -542,10 +542,10 @@ RubberOperator::~RubberOperator()
 
 
 // Inline visualization
-void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
+void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
                GridFunction *field, const char *field_name, bool init_vis)
 {
-   if (!out)
+   if (!os)
    {
       return;
    }
@@ -555,23 +555,25 @@ void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
 
    mesh->SwapNodes(nodes, owns_nodes);
 
-   out << "solution\n" << *mesh << *field;
+   os << "solution\n" << *mesh << *field;
 
    mesh->SwapNodes(nodes, owns_nodes);
 
    if (init_vis)
    {
-      out << "window_size 800 800\n";
-      out << "window_title '" << field_name << "'\n";
+      os << "window_size 800 800\n";
+      os << "window_title '" << field_name << "'\n";
       if (mesh->SpaceDimension() == 2)
       {
-         out << "view 0 0\n"; // view from top
-         out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
+         os << "view 0 0\n"; // view from top
+         // turn off perspective and light, +anti-aliasing
+         os << "keys jlA\n";
       }
-      out << "keys cmA\n";        // show colorbar and mesh, +anti-aliasing
-      out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
+      os << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
+      // update value-range; keep mesh-extents fixed
+      os << "autoscale value\n";
    }
-   out << flush;
+   os << flush;
 }
 
 void ReferenceConfiguration(const Vector &x, Vector &y)

examples/ex19p.cpp

@@ -185,7 +185,8 @@ public:
 };
 
 // Visualization driver
-void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
+void visualize(ostream &os, ParMesh *mesh,
+               ParGridFunction *deformed_nodes,
                ParGridFunction *field, const char *field_name = NULL,
                bool init_vis = false);
 
@@ -617,10 +618,11 @@ RubberOperator::~RubberOperator()
 
 
 // Inline visualization
-void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
+void visualize(ostream &os, ParMesh *mesh,
+               ParGridFunction *deformed_nodes,
                ParGridFunction *field, const char *field_name, bool init_vis)
 {
-   if (!out)
+   if (!os)
    {
       return;
    }
@@ -630,24 +632,27 @@ void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
 
    mesh->SwapNodes(nodes, owns_nodes);
 
-   out << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
-   out << "solution\n" << *mesh << *field;
+   os << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() <<
+      "\n";
+   os << "solution\n" << *mesh << *field;
 
    mesh->SwapNodes(nodes, owns_nodes);
 
    if (init_vis)
    {
-      out << "window_size 800 800\n";
-      out << "window_title '" << field_name << "'\n";
+      os << "window_size 800 800\n";
+      os << "window_title '" << field_name << "'\n";
       if (mesh->SpaceDimension() == 2)
       {
-         out << "view 0 0\n"; // view from top
-         out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
+         os << "view 0 0\n"; // view from top
+         // turn off perspective and light, +anti-aliasing
+         os << "keys jlA\n";
       }
-      out << "keys cmA\n";        // show colorbar and mesh, +anti-aliasing
-      out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
+      os << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
+      // update value-range; keep mesh-extents fixed
+      os << "autoscale value\n";
    }
-   out << flush;
+   os << flush;
 }
 
 void ReferenceConfiguration(const Vector &x, Vector &y)

examples/ex27.cpp

@@ -75,7 +75,7 @@ Mesh * GenerateSerialMesh(int ref);
 // alpha*n.Grad(sol) + beta*sol - gamma over the same boundary.
 double IntegrateBC(const GridFunction &sol, const Array<int> &bdr_marker,
                    double alpha, double beta, double gamma,
-                   double &err);
+                   double &error);
 
 int main(int argc, char *argv[])
 {
@@ -302,26 +302,26 @@ int main(int argc, char *argv[])
    {
       // Integrate the solution on the Dirichlet boundary and compare to the
       // expected value.
-      double err, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, err);
+      double error, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, error);
 
       bool hom_dbc = (dbc_val == 0.0);
-      err /=  hom_dbc ? 1.0 : fabs(dbc_val);
+      error /=  hom_dbc ? 1.0 : fabs(dbc_val);
       mfem::out << "Average of solution on Gamma_dbc:\t"
                 << avg << ", \t"
                 << (hom_dbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
    {
       // Integrate n.Grad(u) on the inhomogeneous Neumann boundary and compare
       // to the expected value.
-      double err, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, err);
+      double error, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, error);
 
       bool hom_nbc = (nbc_val == 0.0);
-      err /=  hom_nbc ? 1.0 : fabs(nbc_val);
+      error /=  hom_nbc ? 1.0 : fabs(nbc_val);
       mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
                 << avg << ", \t"
                 << (hom_nbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
    {
       // Integrate n.Grad(u) on the homogeneous Neumann boundary and compare to
@@ -330,25 +330,26 @@ int main(int argc, char *argv[])
       nbc0_bdr = 0;
       nbc0_bdr[3] = 1;
 
-      double err, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, err);
+      double error, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, error);
 
       bool hom_nbc = true;
       mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
                 << avg << ", \t"
                 << (hom_nbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
    {
       // Integrate n.Grad(u) + a * u on the Robin boundary and compare to the
       // expected value.
-      double err, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, err);
+      double error;
+      double avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, error);
 
       bool hom_rbc = (rbc_b_val == 0.0);
-      err /=  hom_rbc ? 1.0 : fabs(rbc_b_val);
+      error /=  hom_rbc ? 1.0 : fabs(rbc_b_val);
       mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
                 << avg << ", \t"
                 << (hom_rbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
 
    // 14. Save the refined mesh and the solution. This output can be viewed
@@ -637,11 +638,11 @@ Mesh * GenerateSerialMesh(int ref)
 
 double IntegrateBC(const GridFunction &x, const Array<int> &bdr,
                    double alpha, double beta, double gamma,
-                   double &err)
+                   double &error)
 {
    double nrm = 0.0;
    double avg = 0.0;
-   err = 0.0;
+   error = 0.0;
 
    const bool a_is_zero = alpha == 0.0;
    const bool b_is_zero = beta == 0.0;
@@ -705,20 +706,20 @@ double IntegrateBC(const GridFunction &x, const Array<int> &bdr,
 
          // Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
          val -= gamma;
-         err += (val*val) * ip.weight * face_weight;
+         error += (val*val) * ip.weight * face_weight;
       }
    }
 
    // Normalize by the length of the boundary
    if (std::abs(nrm) > 0.0)
    {
-      err /= nrm;
+      error /= nrm;
       avg /= nrm;
    }
 
    // Compute l2 norm of the error in the boundary condition (negative
-   // quadrature weights may produce negative 'err')
-   err = (err >= 0.0) ? sqrt(err) : -sqrt(-err);
+   // quadrature weights may produce negative 'error')
+   error = (error >= 0.0) ? sqrt(error) : -sqrt(-error);
 
    // Return the average value of alpha * n.Grad(x) + beta * x
    return avg;

examples/ex27p.cpp

@@ -75,7 +75,7 @@ Mesh * GenerateSerialMesh(int ref);
 // alpha*n.Grad(sol) + beta*sol - gamma over the same boundary.
 double IntegrateBC(const ParGridFunction &sol, const Array<int> &bdr_marker,
                    double alpha, double beta, double gamma,
-                   double &err);
+                   double &error);
 
 int main(int argc, char *argv[])
 {
@@ -321,26 +321,26 @@ int main(int argc, char *argv[])
    {
       // Integrate the solution on the Dirichlet boundary and compare to the
       // expected value.
-      double err, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, err);
+      double error, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, error);
 
       bool hom_dbc = (dbc_val == 0.0);
-      err /=  hom_dbc ? 1.0 : fabs(dbc_val);
+      error /=  hom_dbc ? 1.0 : fabs(dbc_val);
       mfem::out << "Average of solution on Gamma_dbc:\t"
                 << avg << ", \t"
                 << (hom_dbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
    {
       // Integrate n.Grad(u) on the inhomogeneous Neumann boundary and compare
       // to the expected value.
-      double err, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, err);
+      double error, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, error);
 
       bool hom_nbc = (nbc_val == 0.0);
-      err /=  hom_nbc ? 1.0 : fabs(nbc_val);
+      error /=  hom_nbc ? 1.0 : fabs(nbc_val);
       mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
                 << avg << ", \t"
                 << (hom_nbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
    {
       // Integrate n.Grad(u) on the homogeneous Neumann boundary and compare to
@@ -349,25 +349,26 @@ int main(int argc, char *argv[])
       nbc0_bdr = 0;
       nbc0_bdr[3] = 1;
 
-      double err, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, err);
+      double error, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, error);
 
       bool hom_nbc = true;
       mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
                 << avg << ", \t"
                 << (hom_nbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
    {
       // Integrate n.Grad(u) + a * u on the Robin boundary and compare to the
       // expected value.
-      double err, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, err);
+      double error, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val,
+                                      error);
 
       bool hom_rbc = (rbc_b_val == 0.0);
-      err /=  hom_rbc ? 1.0 : fabs(rbc_b_val);
+      error /=  hom_rbc ? 1.0 : fabs(rbc_b_val);
       mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
                 << avg << ", \t"
                 << (hom_rbc ? "absolute" : "relative")
-                << " error " << err << endl;
+                << " error " << error << endl;
    }
 
    // 15. Save the refined mesh and the solution in parallel. This output can be
@@ -667,11 +668,11 @@ double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
    double loc_vals[3];
    double &nrm = loc_vals[0];
    double &avg = loc_vals[1];
-   double &err = loc_vals[2];
+   double &error = loc_vals[2];
 
    nrm = 0.0;
    avg = 0.0;
-   err = 0.0;
+   error = 0.0;
 
    const bool a_is_zero = alpha == 0.0;
    const bool b_is_zero = beta == 0.0;
@@ -735,7 +736,7 @@ double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
 
          // Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
          val -= gamma;
-         err += (val*val) * ip.weight * face_weight;
+         error += (val*val) * ip.weight * face_weight;
       }
    }
 
@@ -754,7 +755,7 @@ double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
    }
 
    // Compute l2 norm of the error in the boundary condition (negative
-   // quadrature weights may produce negative 'err')
+   // quadrature weights may produce negative 'error')
    glb_err = (glb_err >= 0.0) ? sqrt(glb_err) : -sqrt(-glb_err);
 
    // Return the average value of alpha * n.Grad(x) + beta * x

examples/ex29.cpp

@@ -167,9 +167,9 @@ int main(int argc, char *argv[])
 
    // 13. Compute error in the solution and its flux
    FunctionCoefficient uCoef(uExact);
-   double err = x.ComputeL2Error(uCoef);
+   double error = x.ComputeL2Error(uCoef);
 
-   cout << "|u - u_h|_2 = " << err << endl;
+   cout << "|u - u_h|_2 = " << error << endl;
 
    FiniteElementSpace flux_fespace(mesh, &fec, 3);
    GridFunction flux(&flux_fespace);

examples/ex29p.cpp

@@ -197,9 +197,9 @@ int main(int argc, char *argv[])
 
    // 15. Compute error in the solution and its flux
    FunctionCoefficient uCoef(uExact);
-   double err = x.ComputeL2Error(uCoef);
+   double error = x.ComputeL2Error(uCoef);
 
-   if (myid == 0) { cout << "|u - u_h|_2 = " << err << endl; }
+   if (myid == 0) { cout << "|u - u_h|_2 = " << error << endl; }
 
    ParFiniteElementSpace flux_fespace(&pmesh, &fec, 3);
    ParGridFunction flux(&flux_fespace);

examples/ex3p.cpp

@@ -255,10 +255,10 @@ int main(int argc, char *argv[])
 
    // 15. Compute and print the L^2 norm of the error.
    {
-      double err = x.ComputeL2Error(E);
+      double error = x.ComputeL2Error(E);
       if (myid == 0)
       {
-         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+         cout << "\n|| E_h - E ||_{L^2} = " << error << '\n' << endl;
       }
    }
 

examples/ex4p.cpp

@@ -256,10 +256,10 @@ int main(int argc, char *argv[])
 
    // 15. Compute and print the L^2 norm of the error.
    {
-      double err = x.ComputeL2Error(F);
+      double error = x.ComputeL2Error(F);
       if (myid == 0)
       {
-         cout << "\n|| F_h - F ||_{L^2} = " << err << '\n' << endl;
+         cout << "\n|| F_h - F ||_{L^2} = " << error << '\n' << endl;
       }
    }
 

examples/ex7p.cpp

@@ -282,10 +282,10 @@ int main(int argc, char *argv[])
    delete b;
 
    // 12. Compute and print the L^2 norm of the error.
-   double err = x.ComputeL2Error(sol_coef);
+   double error = x.ComputeL2Error(sol_coef);
    if (myid == 0)
    {
-      cout << "\nL2 norm of error: " << err << endl;
+      cout << "\nL2 norm of error: " << error << endl;
    }
 
    // 13. Save the refined mesh and the solution. This output can be viewed

fem/coefficient.cpp

@@ -21,6 +21,33 @@ namespace mfem
 
 using namespace std;
 
+// Given an ElementTransformation and IntegrationPoint in a refined mesh,
+// return the ElementTransformation of the parent coarse element, and set
+// coarse_ip to the location of the original ip within the coarse element.
+ElementTransformation *RefinedToCoarse(
+   Mesh &coarse_mesh, const ElementTransformation &T,
+   const IntegrationPoint &ip, IntegrationPoint &coarse_ip)
+{
+   Mesh &fine_mesh = *T.mesh;
+   // Get the element transformation of the coarse element containing the
+   // fine element.
+   int fine_element = T.ElementNo;
+   const CoarseFineTransformations &cf = fine_mesh.GetRefinementTransforms();
+   int coarse_element = cf.embeddings[fine_element].parent;
+   ElementTransformation *coarse_T = coarse_mesh.GetElementTransformation(
+                                        coarse_element);
+   // Transform the integration point from fine element coordinates to coarse
+   // element coordinates.
+   Geometry::Type geom = T.GetGeometryType();
+   IntegrationPointTransformation fine_to_coarse;
+   IsoparametricTransformation &emb_tr = fine_to_coarse.Transf;
+   emb_tr.SetIdentityTransformation(geom);
+   emb_tr.SetPointMat(cf.point_matrices[geom](cf.embeddings[fine_element].matrix));
+   fine_to_coarse.Transform(ip, coarse_ip);
+   coarse_T->SetIntPoint(&coarse_ip);
+   return coarse_T;
+}
+
 double PWConstCoefficient::Eval(ElementTransformation & T,
                                 const IntegrationPoint & ip)
 {
@@ -95,7 +122,17 @@ double FunctionCoefficient::Eval(ElementTransformation & T,
 double GridFunctionCoefficient::Eval (ElementTransformation &T,
                                       const IntegrationPoint &ip)
 {
-   return GridF -> GetValue (T, ip, Component);
+   Mesh *gf_mesh = GridF->FESpace()->GetMesh();
+   if (T.mesh == gf_mesh)
+   {
+      return GridF->GetValue(T, ip, Component);
+   }
+   else
+   {
+      IntegrationPoint coarse_ip;
+      ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
+      return GridF->GetValue(*coarse_T, coarse_ip, Component);
+   }
 }
 
 void TransformedCoefficient::SetTime(double t)
@@ -305,13 +342,30 @@ void VectorGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
 void VectorGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
                                          const IntegrationPoint &ip)
 {
-   GridFunc->GetVectorValue(T, ip, V);
+   Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
+   if (T.mesh == gf_mesh)
+   {
+      GridFunc->GetVectorValue(T, ip, V);
+   }
+   else
+   {
+      IntegrationPoint coarse_ip;
+      ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
+      GridFunc->GetVectorValue(*coarse_T, coarse_ip, V);
+   }
 }
 
 void VectorGridFunctionCoefficient::Eval(
    DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
 {
-   GridFunc->GetVectorValues(T, ir, M);
+   if (T.mesh == GridFunc->FESpace()->GetMesh())
+   {
+      GridFunc->GetVectorValues(T, ir, M);
+   }
+   else
+   {
+      VectorCoefficient::Eval(M, T, ir);
+   }
 }
 
 GradientGridFunctionCoefficient::GradientGridFunctionCoefficient (
@@ -331,13 +385,30 @@ void GradientGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
 void GradientGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
                                            const IntegrationPoint &ip)
 {
-   GridFunc->GetGradient(T, V);
+   Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
+   if (T.mesh == gf_mesh)
+   {
+      GridFunc->GetGradient(T, V);
+   }
+   else
+   {
+      IntegrationPoint coarse_ip;
+      ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
+      GridFunc->GetGradient(*coarse_T, V);
+   }
 }
 
 void GradientGridFunctionCoefficient::Eval(
    DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
 {
-   GridFunc->GetGradients(T, ir, M);
+   if (T.mesh == GridFunc->FESpace()->GetMesh())
+   {
+      GridFunc->GetGradients(T, ir, M);
+   }
+   else
+   {
+      VectorCoefficient::Eval(M, T, ir);
+   }
 }
 
 CurlGridFunctionCoefficient::CurlGridFunctionCoefficient(
@@ -363,7 +434,17 @@ void CurlGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
 void CurlGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
                                        const IntegrationPoint &ip)
 {
-   GridFunc->GetCurl(T, V);
+   Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
+   if (T.mesh == gf_mesh)
+   {
+      GridFunc->GetCurl(T, V);
+   }
+   else
+   {
+      IntegrationPoint coarse_ip;
+      ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
+      GridFunc->GetCurl(*coarse_T, V);
+   }
 }
 
 DivergenceGridFunctionCoefficient::DivergenceGridFunctionCoefficient (
@@ -375,7 +456,17 @@ DivergenceGridFunctionCoefficient::DivergenceGridFunctionCoefficient (
 double DivergenceGridFunctionCoefficient::Eval(ElementTransformation &T,
                                                const IntegrationPoint &ip)
 {
-   return GridFunc->GetDivergence(T);
+   Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
+   if (T.mesh == gf_mesh)
+   {
+      return GridFunc->GetDivergence(T);
+   }
+   else
+   {
+      IntegrationPoint coarse_ip;
+      ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
+      return GridFunc->GetDivergence(*coarse_T);
+   }
 }
 
 void VectorDeltaCoefficient::SetTime(double t)

fem/coefficient.hpp

@@ -574,6 +574,7 @@ public:
    /// Evaluate the coefficient.
    virtual void Eval(Vector &V, ElementTransformation &T,
                      const IntegrationPoint &ip);
+   using VectorCoefficient::Eval;
 };
 
 /// A general vector function coefficient

fem/conduitdatacollection.cpp

@@ -330,7 +330,7 @@ ConduitDataCollection::BlueprintMeshToMesh(const Node &n_mesh,
          }
          else
          {
-            Node &(n_bndry_conn_conv) =
+            Node &n_bndry_conn_conv =
                n_conv["topologies"][bndry_topo_name]["elements/connectivity"];
             n_bndry_conn.to_int_array(n_bndry_conn_conv);
             bndry_indices = (n_bndry_conn_conv).value();

fem/doftrans.hpp

@@ -243,6 +243,7 @@ public:
    void TransformDual(double *v) const;
 
    void InvTransformDual(double *v) const;
+   using DofTransformation::InvTransformDual;
 };
 
 /// DoF transformation implementation for the Nedelec basis on tetrahedra

fem/eltrans.cpp

@@ -21,7 +21,8 @@ ElementTransformation::ElementTransformation()
      EvalState(0),
      geom(Geometry::INVALID),
      Attribute(-1),
-     ElementNo(-1)
+     ElementNo(-1),
+     mesh(nullptr)
 { }
 
 double ElementTransformation::EvalWeight()

fem/eltrans.hpp

@@ -75,6 +75,12 @@ public:
 
    int Attribute, ElementNo, ElementType;
 
+   /// The Mesh object containing the element.
+   /** If the element transformation belongs to a mesh, this will point to the
+       containing Mesh object. ElementNo will be the number of the element in
+       this Mesh. This will be NULL if the element does not belong to a mesh. */
+   class Mesh *mesh;
+
    ElementTransformation();
 
    /** @brief Force the reevaluation of the Jacobian in the next call. */

fem/fe/fe_base.cpp

@@ -1993,6 +1993,57 @@ void Poly_1D::CalcDBinomTerms(const int p, const double x, const double y,
    }
 }
 
+void Poly_1D::CalcJacobi(const int p, const double alpha, const double beta,
+                         const double x, double *u)
+{
+   // use the recursive definition for [-1,1]:
+   // 2(n+1)(n+a+b+1)(2n+a+b)P^{(a,b)}_{n+1}(z) =
+   //    (2n+a+b+1)((2n+a+b+2)(2n+a+b)z+a^2-b^2)P^{(a,b)}_n(z)
+   //    -2(n+a)(n+b)(2n+a+b+2)P^{(a,b)}_{n-1}(z)
+   u[0] = 1.;
+   if (p == 0) { return; }
+   double z = 2.*x - 1.;
+   u[1] = 0.5 * ((alpha + beta + 2.) * z + alpha - beta);
+   for (int n = 1; n < p; n++)
+   {
+      double c2 = 2.*(1. + n)*(1. + n + alpha + beta)*(2.*n + alpha + beta);
+      double c1 = (1. + 2.*n + alpha + beta) *
+                  ((2. + 2.*n + alpha + beta)*
+                   (2.*n + alpha + beta)*z + pow(alpha,2) - pow(beta,2));
+      double c0 = 2.*(alpha + n)*(beta + n)*(2. + 2.*n + alpha + beta);
+
+      u[n+1] = (c1 * u[n] - c0 * u[n-1]) / c2;
+   }
+}
+
+void Poly_1D::CalcJacobi(const int p, const double alpha, const double beta,
+                         const double x, double *u, double *d)
+{
+   // use the recursive definition for [-1,1]:
+   // 2(n+1)(n+a+b+1)(2n+a+b)P^{(a,b)}_{n+1}(z) =
+   //    (2n+a+b+1)((2n+a+b+2)(2n+a+b)z+a^2-b^2)P^{(a,b)}_n(z)
+   //    -2(n+a)(n+b)(2n+a+b+2)P^{(a,b)}_{n-1}(z)
+   u[0] = 1.;
+   d[0] = 0.;
+   if (p == 0) { return; }
+   double z = 2.*x - 1.;
+   u[1] = 0.5 * ((alpha + beta + 2.) * z + alpha - beta);
+   d[1] = alpha + beta + 2.;
+   for (int n = 1; n < p; n++)
+   {
+      double c2 = 2.*(1. + n)*(1. + n + alpha + beta)*(2.*n + alpha + beta);
+      double c1 = (1. + 2.*n + alpha + beta) *
+                  ((2. + 2.*n + alpha + beta)*
+                   (2.*n + alpha + beta)*z + pow(alpha,2) - pow(beta,2));
+      double c0 = 2.*(alpha + n)*(beta + n)*(2. + 2.*n + alpha + beta);
+      double dc1 = 2.*(1. + 2.*n + alpha + beta) *
+                   (2. + 2.*n + alpha + beta)*(2.*n + alpha + beta);
+
+      u[n+1] = (c1 * u[n] - c0 * u[n-1]) / c2;
+      d[n+1] = (c1 * d[n] + dc1 * u[n] - c0 * d[n-1]) / c2;
+   }
+}
+
 void Poly_1D::CalcLegendre(const int p, const double x, double *u)
 {
    // use the recursive definition for [-1,1]:

fem/fe/fe_base.hpp

@@ -1112,7 +1112,28 @@ public:
    static void CalcBernstein(const int p, const double x, double *u, double *d)
    { CalcBinomTerms(p, x, 1. - x, u, d); }
 
+   // Evaluate the values of the Jacobi polynomials, P^{(\alpha,\beta)}_n,
+   // on [0,1] with degrees ranging from 0 to p and parameters alpha and beta
+   // at a point x.  Internally the point is transformed to the traditional
+   // domain of [-1,1].
+   static void CalcJacobi(const int p, const double alpha, const double beta,
+                          const double x, double *u);
+
+   // Evaluate the values and derivatives of the Jacobi polynomials,
+   // P^{(\alpha,\beta)}_n, on [0,1] with degrees ranging from 0 to p and
+   // parameters alpha and beta at a point x.  Internally the point is
+   // transformed to the traditional domain of [-1,1].
+   static void CalcJacobi(const int p, const double alpha, const double beta,
+                          const double x, double *u, double *d);
+
+   // Evaluate the values of the Legendre polynomials, P_n, on [0,1] with
+   // degrees ranging from 0 to p at a point x.  Internally the point is
+   // transformed to the traditional domain of [-1,1].
    static void CalcLegendre(const int p, const double x, double *u);
+
+   // Evaluate the values and derivatives of the Legendre polynomials, P_n,
+   // on [0,1] with degrees ranging from 0 to p at a point x.  Internally the
+   // point is transformed to the traditional domain of [-1,1].
    static void CalcLegendre(const int p, const double x, double *u, double *d);
 
    ~Poly_1D();

fem/fe/fe_h1.cpp

@@ -1040,4 +1040,315 @@ void H1_WedgeElement::CalcDShape(const IntegrationPoint &ip,
    }
 }
 
+H1_PyramidElement::H1_PyramidElement(const int p, const int btype)
+   : NodalFiniteElement(3, Geometry::PYRAMID,
+                        (p + 1) * (p + 2) * (2 * p + 3) / 6,
+                        p, FunctionSpace::Qk)
+{
+   const double *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+
+   const double **mcp = new const double*[p];
+   mcp[0] = NULL;
+   for (int k=1; k<p; k++)
+   {
+      mcp[k] = poly1d.ClosedPoints(k, VerifyNodal(VerifyClosed(btype)));
+   }
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(p + 1);
+   shape_y.SetSize(p + 1);
+   shape_z.SetSize(p + 1);
+   dshape_x.SetSize(p + 1);
+   dshape_y.SetSize(p + 1);
+   dshape_z.SetSize(p + 1);
+   ddshape_x.SetSize(p + 1);
+   ddshape_y.SetSize(p + 1);
+   ddshape_z.SetSize(p + 1);
+   u.SetSize(dof);
+   du.SetSize(dof, dim);
+   ddu.SetSize(dof, (dim * (dim + 1)) / 2);
+#else
+   Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
+   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
+   Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
+   Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
+   Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
+
+   // edges
+   int o = 5;
+   for (int i = 1; i < p; i++)  // (0,1)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (1,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (3,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (1,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (2,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[p-i], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (3,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
+   }
+
+   // quadrilateral face
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[0]);
+      }
+   }
+
+   // triangular faces
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (0,1,4)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
+      }
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (1,2,4)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(1.0 - cp[j]/w, cp[i]/w, cp[j]/w);
+      }
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (3,4,2)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[j]/w, 1.0 - cp[i]/w, cp[i]/w);
+      }
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (0,4,3)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
+      }
+   // interior
+   for (int k = 1; k < p - 1; k++)
+   {
+      double wk = 1.0 - cp[k];
+      for (int j = 1; j < p - k; j++)
+      {
+         for (int i = 1; i < p - k; i++)
+         {
+            Nodes.IntPoint(o++).Set3(mcp[p-k][i] * wk,
+                                     mcp[p-k][j] * wk,
+                                     cp[k]);
+         }
+      }
+   }
+
+   MFEM_ASSERT(o == dof,
+               "Number of nodes does not match the "
+               "number of degrees of freedom");
+   DenseMatrix T(dof);
+
+   for (int m = 0; m < dof; m++)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(m);
+
+      double oz = 1.0 - ip.z;
+
+      double tol = 1e-6;
+
+      double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+      double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+      double z = ip.z;
+
+      poly1d.CalcLegendre(p, x, shape_x);
+      poly1d.CalcLegendre(p, y, shape_y);
+
+      o = 0;
+      for (int i = 0; i <= p; i++)
+      {
+         for (int j = 0; j <= p; j++)
+         {
+            int maxij = std::max(i, j);
+            poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z, shape_z);
+
+            for (int k = 0; k <= p - maxij; k++)
+            {
+               if (oz <= tol)
+               {
+                  if (maxij == 0)
+                  {
+                     T(o++, m) = shape_z(k);
+                  }
+                  else
+                  {
+                     T(o++, m) = 0.0;
+                  }
+               }
+               else
+               {
+                  T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
+                              pow(1.0 - ip.z, maxij);
+               }
+            }
+         }
+      }
+   }
+   Ti.Factor(T);
+
+   delete [] mcp;
+}
+
+void H1_PyramidElement::CalcShape(const IntegrationPoint &ip,
+                                  Vector &shape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(order+1);
+   Vector shape_y(order+1);
+   Vector shape_z(order+1);
+   Vector u(dof);
+#endif
+
+   double oz = 1.0 - ip.z;
+
+   double tol = 1e-6;
+
+   double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+   double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+   double z = ip.z;
+
+   poly1d.CalcLegendre(p, x, shape_x);
+   poly1d.CalcLegendre(p, y, shape_y);
+
+   for (int o = 0, i = 0; i <= p; i++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z, shape_z);
+         for (int k = 0; k <= p - maxij; k++)
+         {
+            if (oz <= tol)
+            {
+               if (maxij == 0)
+               {
+                  u(o++) = shape_z(k);
+               }
+               else
+               {
+                  u(o++) = 0.0;
+               }
+            }
+            else
+            {
+               u(o++) =
+                  shape_x(i) * shape_y(j) * shape_z(k) * pow(1.0 - ip.z, maxij);
+            }
+         }
+      }
+   }
+
+   Ti.Mult(u, shape);
+}
+
+void H1_PyramidElement::CalcDShape(const IntegrationPoint &ip,
+                                   DenseMatrix &dshape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector  shape_x(p + 1),  shape_y(p + 1),  shape_z(p + 1);
+   Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1);
+   DenseMatrix du(dof, dim);
+#endif
+
+   double oz = 1.0 - ip.z;
+
+   double tol = 1e-6;
+
+   if (oz <= tol)
+   {
+      mfem::out << "H1_PyramidElement::CalcDShape near apex" << std::endl;
+   }
+
+   double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+   double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+   double z = ip.z;
+
+   poly1d.CalcLegendre(p, x, shape_x, dshape_x);
+   poly1d.CalcLegendre(p, y, shape_y, dshape_y);
+
+   for (int o = 0, i = 0; i <= p; i++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z,
+                           shape_z, dshape_z);
+         for (int k = 0; k <= p - maxij; k++)
+         {
+            if (oz > tol)
+            {
+               du(o, 0) = dshape_x(i) * shape_y(j) * shape_z(k) *
+                          pow(1.0 - ip.z, maxij - 1);
+               du(o, 1) = shape_x(i) * dshape_y(j) * shape_z(k) *
+                          pow(1.0 - ip.z, maxij - 1);
+               du(o, 2) = (ip.x * dshape_x(i) * shape_y(j) * shape_z(k) +
+                           ip.y * shape_x(i) * dshape_y(j) * shape_z(k) -
+                           shape_x(i) * shape_y(j) * (1.0 - ip.z) *
+                           (shape_z(k) * maxij -
+                            (1.0 - ip.z) * dshape_z(k))
+                          ) * pow(1.0 - ip.z, maxij - 2);
+            }
+            else
+            {
+               // The gradient is multivalued at the apex. The following
+               // values were produced by setting x and y equal to zero
+               // and computing the limit as z->1.
+
+               du(o, 0) = 0.0;
+               du(o, 1) = 0.0;
+               du(o, 2) = 0.0;
+
+               if (i == 0 && j == 0)
+               {
+                  du(o, 2) = k * (k + 1) * (k + 2) * (k + 3) / 6;
+               }
+               else if (i == 0 && j == 1)
+               {
+                  du(o, 1) = -(k + 1) * (k + 2) * (k + 3) * (k + 4) / 24;
+               }
+               else if (i == 1 && j == 0)
+               {
+                  du(o, 0) = -(k + 1) * (k + 2) * (k + 3) * (k + 4) / 24;
+               }
+            }
+            o++;
+         }
+      }
+   }
+
+   Ti.Mult(du, dshape);
+}
+
 }

fem/fe/fe_h1.hpp

@@ -148,6 +148,31 @@ public:
                            DenseMatrix &dshape) const;
 };
 
+class H1_PyramidElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   // mutable Vector shape_0, shape_1, shape_2;
+   // mutable Vector dshape_0_0, dshape_1_0, dshape_2_0;
+   // mutable Vector dshape_0_1, dshape_1_1, dshape_2_1;
+   // mutable Vector u;
+   // mutable DenseMatrix du;
+
+   mutable Vector shape_x, shape_y, shape_z;
+   mutable Vector dshape_x, dshape_y, dshape_z, u;
+   mutable Vector ddshape_x, ddshape_y, ddshape_z;
+   mutable DenseMatrix du, ddu;
+#endif
+   DenseMatrixInverse Ti;
+
+public:
+   H1_PyramidElement(const int p,
+                     const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+};
+
 } // namespace mfem
 
 #endif

fem/fe_coll.cpp

@@ -1749,7 +1749,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
          H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
          H1_dof[Geometry::CUBE] = QuadDof*pm1;
          H1_dof[Geometry::PRISM] = TriDof*pm1;
-         H1_dof[Geometry::PYRAMID] = 0;
+         H1_dof[Geometry::PYRAMID] = pm1*pm2*(2*p-3)/6;
          if (b_type == BasisType::Positive)
          {
             H1_Elements[Geometry::TETRAHEDRON] = new H1Pos_TetrahedronElement(p);
@@ -1763,7 +1763,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
             H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p, btype);
             H1_Elements[Geometry::PRISM] = new H1_WedgeElement(p, btype);
          }
-         H1_Elements[Geometry::PYRAMID] = new LinearPyramidFiniteElement;
+         H1_Elements[Geometry::PYRAMID] = new H1_PyramidElement(p, btype);
 
          const int &TetDof = H1_dof[Geometry::TETRAHEDRON];
          TetDofOrd[0] = new int[24*TetDof];
@@ -1861,16 +1861,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
 const FiniteElement *
 H1_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
 {
-   if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
-   {
-      return H1_Elements[GeomType];
-   }
-   else
-   {
-      MFEM_ABORT("H1 Pyramid basis functions are not yet supported "
-                 "for order > 1.");
-      return NULL;
-   }
+   return H1_Elements[GeomType];
 }
 
 const int *H1_FECollection::DofOrderForOrientation(Geometry::Type GeomType,

fem/fmsconvert.cpp

@@ -112,7 +112,7 @@ FmsFieldToGridFunction(FmsMesh fms_mesh, FmsField f, Mesh *mesh,
    // NOTE: transplanted from the FmsMeshToMesh function
    //       We should do this work once and save it.
    //--------------------------------------------------
-   FmsInt dim, n_vert, n_elem, space_dim;
+   FmsInt dim, n_elem, space_dim;
 
    // Find the first component that has coordinates - that will be the new mfem
    // mesh.
@@ -144,7 +144,6 @@ FmsFieldToGridFunction(FmsMesh fms_mesh, FmsField f, Mesh *mesh,
          n_ents[et] += num_ents;
       }
    }
-   n_vert = n_ents[FMS_VERTEX];
    //--------------------------------------------------
 
    // Interrogate the field.

fem/intrules.cpp

@@ -1675,7 +1675,7 @@ IntegrationRule *IntegrationRules::PyramidIntegrationRule(int Order)
       ipp.x = ipc.x * (1.0 - ipc.z);
       ipp.y = ipc.y * (1.0 - ipc.z);
       ipp.z = ipc.z;
-      ipp.weight = ipc.weight / 3.0;
+      ipp.weight = 4.0 * ipc.weight * pow(1.0 - ipc.z, 2);
    }
    return PyramidIntRules[Order];
 }

fem/pfespace.cpp

@@ -1010,12 +1010,9 @@ void ParFiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
 {
    FiniteElementSpace::GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component);
 
-   if (Conforming())
-   {
-      // Make sure that processors without boundary elements mark
-      // their boundary dofs (if they have any).
-      Synchronize(ess_dofs);
-   }
+   // Make sure that processors without boundary elements mark
+   // their boundary dofs (if they have any).
+   Synchronize(ess_dofs);
 }
 
 void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
@@ -1041,7 +1038,8 @@ void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
    {
       if (bool(ted[i]) != bool(true_ess_dofs2[i])) { counter++; }
    }
-   MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter);
+   MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter
+               << ", rank = " << MyRank);
 #endif
 
    MarkerToList(true_ess_dofs, ess_tdof_list);

fem/qinterp/eval.hpp

@@ -65,13 +65,13 @@ static void Values2D(const int NE,
       MFEM_SHARED double sm0[NBZ][MDQ*MDQ];
       MFEM_SHARED double sm1[NBZ][MDQ*MDQ];
 
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
+
       ConstDeviceMatrix B(sB, D1D,Q1D);
       DeviceMatrix DD(sm0[tidz], MD1, MD1);
       DeviceMatrix DQ(sm1[tidz], MD1, MQ1);
       DeviceMatrix QQ(sm0[tidz], MQ1, MQ1);
 
-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
-
       for (int c = 0; c < VDIM; c++)
       {
          kernels::internal::LoadX(e,D1D,c,x,DD);
@@ -126,14 +126,14 @@ static void Values3D(const int NE,
       MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
       MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
 
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
+
       ConstDeviceMatrix B(sB, D1D,Q1D);
       DeviceCube DDD(sm0, MD1,MD1,MD1);
       DeviceCube DDQ(sm1, MD1,MD1,MQ1);
       DeviceCube DQQ(sm0, MD1,MQ1,MQ1);
       DeviceCube QQQ(sm1, MQ1,MQ1,MQ1);
 
-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
-
       for (int c = 0; c < VDIM; c++)
       {
          kernels::internal::LoadX(e,D1D,c,x,DDD);

fem/tmop.hpp

@@ -1133,7 +1133,7 @@ protected:
 
    // Evaluation of the discrete target specification on different meshes.
    // Owned.
-   AdaptivityEvaluator *adapt_eval;
+   AdaptivityEvaluator *adapt_lim_eval;
 
    void SetDiscreteTargetBase(const GridFunction &tspec_);
    void SetTspecAtIndex(int idx, const GridFunction &tspec_);
@@ -1156,7 +1156,7 @@ public:
 #endif
         amr_el(-1), lim_min_size(-0.1),
         good_tspec(false), good_tspec_grad(false), good_tspec_hess(false),
-        adapt_eval(NULL) { }
+        adapt_lim_eval(NULL) { }
 
    virtual ~DiscreteAdaptTC();
 
@@ -1232,8 +1232,8 @@ public:
 
    void SetAdaptivityEvaluator(AdaptivityEvaluator *ae)
    {
-      if (adapt_eval) { delete adapt_eval; }
-      adapt_eval = ae;
+      if (adapt_lim_eval) { delete adapt_lim_eval; }
+      adapt_lim_eval = ae;
    }
 
    const Vector &GetTspecPert1H()   { return tspec_pert1h; }
@@ -1315,15 +1315,15 @@ protected:
    int integ_order;
 
    // Weight Coefficient multiplying the quality metric term.
-   Coefficient *coeff1; // not owned, if NULL -> coeff1 is 1.
+   Coefficient *metric_coeff; // not owned, if NULL -> metric_coeff is 1.
    // Normalization factor for the metric term.
    double metric_normal;
 
    // Nodes and weight Coefficient used for "limiting" the TMOP_Integrator.
    // These are both NULL when there is no limiting.
-   // The class doesn't own nodes0 and coeff0.
-   const GridFunction *nodes0;
-   Coefficient *coeff0;
+   // The class doesn't own lim_nodes0 and lim_coeff.
+   const GridFunction *lim_nodes0;
+   Coefficient *lim_coeff;
    // Limiting reference distance. Not owned.
    const GridFunction *lim_dist;
    // Limiting function. Owned.
@@ -1332,20 +1332,21 @@ protected:
    double lim_normal;
 
    // Adaptive limiting.
-   const GridFunction *zeta_0;       // Not owned.
+   const GridFunction *adapt_lim_gf0;    // Not owned.
 #ifdef MFEM_USE_MPI
-   const ParGridFunction *pzeta_0;
+   const ParGridFunction *adapt_lim_pgf0;
 #endif
-   GridFunction *zeta;               // Owned. Updated by adapt_eval.
-   Coefficient *coeff_zeta;          // Not owned.
-   AdaptivityEvaluator *adapt_eval;  // Not owned.
+   GridFunction *adapt_lim_gf;           // Owned. Updated by adapt_lim_eval.
+   Coefficient *adapt_lim_coeff;         // Not owned.
+   AdaptivityEvaluator *adapt_lim_eval;  // Not owned.
 
    // Surface fitting.
-   GridFunction *sigma, *sigma_bar;  // Owned. Updated by sigma_eval.
-   const Array<bool> *sigma_marker;  // Not owned.
-   Coefficient *coeff_sigma;         // Not owned.
-   AdaptivityEvaluator *sigma_eval;  // Not owned.
-   double sigma_normal;
+   GridFunction *surf_fit_gf,
+                *surf_fit_gf_bar;       // Owned, Updated by surf_fit_eval.
+   const Array<bool> *surf_fit_marker;  // Not owned.
+   Coefficient *surf_fit_coeff;         // Not owned.
+   AdaptivityEvaluator *surf_fit_eval;  // Not owned.
+   double surf_fit_normal;
 
    DiscreteAdaptTC *discr_tc;
 
@@ -1416,7 +1417,7 @@ protected:
 
    void ComputeNormalizationEnergies(const GridFunction &x,
                                      double &metric_energy, double &lim_energy,
-                                     double &sigma_energy);
+                                     double &surf_fit_gf_energy);
 
    void AssembleElementVectorExact(const FiniteElement &el,
                                    ElementTransformation &T,
@@ -1471,7 +1472,7 @@ protected:
 
    void DisableLimiting()
    {
-      nodes0 = NULL; coeff0 = NULL; lim_dist = NULL;
+      lim_nodes0 = NULL; lim_coeff = NULL; lim_dist = NULL;
       delete lim_func; lim_func = NULL;
    }
 
@@ -1531,12 +1532,14 @@ public:
    TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc,
                    TMOP_QualityMetric *hm)
       : h_metric(hm), metric(m), targetC(tc), IntegRules(NULL),
-        integ_order(-1), coeff1(NULL), metric_normal(1.0),
-        nodes0(NULL), coeff0(NULL),
+        integ_order(-1), metric_coeff(NULL), metric_normal(1.0),
+        lim_nodes0(NULL), lim_coeff(NULL),
         lim_dist(NULL), lim_func(NULL), lim_normal(1.0),
-        zeta_0(NULL), zeta(NULL), coeff_zeta(NULL), adapt_eval(NULL),
-        sigma(NULL), sigma_bar(NULL), sigma_marker(NULL), coeff_sigma(NULL),
-        sigma_eval(NULL), sigma_normal(1.0),
+        adapt_lim_gf0(NULL), adapt_lim_gf(NULL), adapt_lim_coeff(NULL),
+        adapt_lim_eval(NULL),
+        surf_fit_gf(NULL), surf_fit_gf_bar(NULL), surf_fit_marker(NULL),
+        surf_fit_coeff(NULL),
+        surf_fit_eval(NULL), surf_fit_normal(1.0),
         discr_tc(dynamic_cast<DiscreteAdaptTC *>(tc)),
         fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false)
    { PA.enabled = false; }
@@ -1564,7 +1567,7 @@ public:
 
        Note that the Coefficient is evaluated in the physical configuration and
        not in the target configuration which may be undefined. */
-   void SetCoefficient(Coefficient &w1) { coeff1 = &w1; }
+   void SetCoefficient(Coefficient &w1) { metric_coeff = &w1; }
 
    /** @brief Limiting of the mesh displacements (general version).
 
@@ -1631,7 +1634,7 @@ public:
    void GetSurfaceFittingErrors(double &err_avg, double &err_max);
 
    /// Update the original/reference nodes used for limiting.
-   void SetLimitingNodes(const GridFunction &n0) { nodes0 = &n0; }
+   void SetLimitingNodes(const GridFunction &n0) { lim_nodes0 = &n0; }
 
    /** @brief Computes the integral of W(Jacobian(Trt)) over a target zone.
        @param[in] el     Type of FiniteElement.

fem/tmop/tmop_pa.cpp

@@ -38,13 +38,13 @@ void TMOP_Integrator::AssembleGradPA(const Vector &xe,
    if (PA.dim == 2)
    {
       AssembleGradPA_2D(xe);
-      if (coeff0) { AssembleGradPA_C0_2D(xe); }
+      if (lim_coeff) { AssembleGradPA_C0_2D(xe); }
    }
 
    if (PA.dim == 3)
    {
       AssembleGradPA_3D(xe);
-      if (coeff0) { AssembleGradPA_C0_3D(xe); }
+      if (lim_coeff) { AssembleGradPA_C0_3D(xe); }
    }
 }
 
@@ -53,8 +53,8 @@ void TMOP_Integrator::AssemblePA_Limiting()
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                          Device::GetDeviceMemoryType() : pa_mt;
    // Return immediately if limiting is not enabled
-   if (coeff0 == nullptr) { return; }
-   MFEM_VERIFY(nodes0, "internal error");
+   if (lim_coeff == nullptr) { return; }
+   MFEM_VERIFY(lim_nodes0, "internal error");
 
    MFEM_VERIFY(PA.enabled, "AssemblePA_Limiting but PA is not enabled!");
    MFEM_VERIFY(lim_func, "No TMOP_LimiterFunction specification!")
@@ -68,14 +68,14 @@ void TMOP_Integrator::AssemblePA_Limiting()
 
    const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
 
-   // H0 for coeff0, (dim x dim) Q-vector
+   // H0 for lim_coeff, (dim x dim) Q-vector
    PA.H0.UseDevice(true);
    PA.H0.SetSize(PA.dim * PA.dim * PA.nq * NE, mt);
 
-   // coeff0 -> PA.C0 (Q-vector)
+   // lim_coeff -> PA.C0 (Q-vector)
    PA.C0.UseDevice(true);
    if (ConstantCoefficient* cQ =
-          dynamic_cast<ConstantCoefficient*>(coeff0))
+          dynamic_cast<ConstantCoefficient*>(lim_coeff))
    {
       PA.C0.SetSize(1, Device::GetMemoryType());
       PA.C0.HostWrite();
@@ -90,17 +90,17 @@ void TMOP_Integrator::AssemblePA_Limiting()
          ElementTransformation& T = *fes->GetElementTransformation(e);
          for (int q = 0; q < ir.GetNPoints(); ++q)
          {
-            C0(q,e) = coeff0->Eval(T, ir.IntPoint(q));
+            C0(q,e) = lim_coeff->Eval(T, ir.IntPoint(q));
          }
       }
    }
 
-   // nodes0 -> PA.X0 (E-vector)
-   MFEM_VERIFY(nodes0->FESpace() == fes, "");
+   // lim_nodes0 -> PA.X0 (E-vector)
+   MFEM_VERIFY(lim_nodes0->FESpace() == fes, "");
    const Operator *n0_R = fes->GetElementRestriction(ordering);
    PA.X0.SetSize(n0_R->Height(), Device::GetMemoryType());
    PA.X0.UseDevice(true);
-   n0_R->Mult(*nodes0, PA.X0);
+   n0_R->Mult(*lim_nodes0, PA.X0);
 
    // Limiting distances: lim_dist -> PA.LD (E-vector)
    // TODO: remove the hack for the case lim_dist == NULL.
@@ -217,8 +217,8 @@ void TMOP_Integrator::AssemblePA(const FiniteElementSpace &fes)
    PA.Jtr_needs_update = true;
    PA.Jtr_debug_grad = false;
 
-   // Limiting: coeff0 -> PA.C0, nodes0 -> PA.X0, lim_dist -> PA.LD, PA.H0
-   if (coeff0) { AssemblePA_Limiting(); }
+   // Limiting: lim_coeff -> PA.C0, lim_nodes0 -> PA.X0, lim_dist -> PA.LD, PA.H0
+   if (lim_coeff) { AssemblePA_Limiting(); }
 }
 
 void TMOP_Integrator::AssembleGradDiagonalPA(Vector &de) const
@@ -236,13 +236,13 @@ void TMOP_Integrator::AssembleGradDiagonalPA(Vector &de) const
    if (PA.dim == 2)
    {
       AssembleDiagonalPA_2D(de);
-      if (coeff0) { AssembleDiagonalPA_C0_2D(de); }
+      if (lim_coeff) { AssembleDiagonalPA_C0_2D(de); }
    }
 
    if (PA.dim == 3)
    {
       AssembleDiagonalPA_3D(de);
-      if (coeff0) { AssembleDiagonalPA_C0_3D(de); }
+      if (lim_coeff) { AssembleDiagonalPA_C0_3D(de); }
    }
 }
 
@@ -258,13 +258,13 @@ void TMOP_Integrator::AddMultPA(const Vector &xe, Vector &ye) const
    if (PA.dim == 2)
    {
       AddMultPA_2D(xe,ye);
-      if (coeff0) { AddMultPA_C0_2D(xe,ye); }
+      if (lim_coeff) { AddMultPA_C0_2D(xe,ye); }
    }
 
    if (PA.dim == 3)
    {
       AddMultPA_3D(xe,ye);
-      if (coeff0) { AddMultPA_C0_3D(xe,ye); }
+      if (lim_coeff) { AddMultPA_C0_3D(xe,ye); }
    }
 }
 
@@ -283,13 +283,13 @@ void TMOP_Integrator::AddMultGradPA(const Vector &re, Vector &ce) const
    if (PA.dim == 2)
    {
       AddMultGradPA_2D(re,ce);
-      if (coeff0) { AddMultGradPA_C0_2D(re,ce); }
+      if (lim_coeff) { AddMultGradPA_C0_2D(re,ce); }
    }
 
    if (PA.dim == 3)
    {
       AddMultGradPA_3D(re,ce);
-      if (coeff0) { AddMultGradPA_C0_3D(re,ce); }
+      if (lim_coeff) { AddMultGradPA_C0_3D(re,ce); }
    }
 }
 
@@ -307,13 +307,13 @@ double TMOP_Integrator::GetLocalStateEnergyPA(const Vector &xe) const
    if (PA.dim == 2)
    {
       energy = GetLocalStateEnergyPA_2D(xe);
-      if (coeff0) { energy += GetLocalStateEnergyPA_C0_2D(xe); }
+      if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_2D(xe); }
    }
 
    if (PA.dim == 3)
    {
       energy = GetLocalStateEnergyPA_3D(xe);
-      if (coeff0) { energy += GetLocalStateEnergyPA_C0_3D(xe); }
+      if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_3D(xe); }
    }
 
    return energy;

fem/tmop/tmop_pa.hpp

@@ -181,9 +181,9 @@ template<int T_D1D, int T_Q1D, int T_MAX> return_t kernel(__VA_ARGS__)
 if (K##kernel.Find(id)) { return K##kernel.At(id)(__VA_ARGS__,0,0); }\
 else {\
    constexpr int T_MAX = 4;\
-   const int D1D = (id>>4)&0xF, Q1D = id&0xF;\
-   MFEM_VERIFY(D1D <= MAX_D1D && Q1D <= MAX_Q1D, "Max size error!");\
-   return kernel<0,0,T_MAX>(__VA_ARGS__,D1D,Q1D); }
+   const int d1d = (id>>4)&0xF, q1d = id&0xF;\
+   MFEM_VERIFY(d1d <= MAX_D1D && q1d <= MAX_Q1D, "Max size error!");\
+   return kernel<0,0,T_MAX>(__VA_ARGS__,d1d,q1d); }
 
 } // namespace kernels
 

fem/tmop/tmop_pa_da3.cpp

@@ -55,6 +55,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, DatcSize,
       MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
       MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
 
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
+
       ConstDeviceMatrix B(sB, D1D, Q1D);
       DeviceCube DDD(sm0, MD1,MD1,MD1);
       DeviceCube DDQ(sm1, MD1,MD1,MQ1);
@@ -88,7 +90,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, DatcSize,
       }
       min = min_size[0];
 
-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
       kernels::internal::EvalX(D1D,Q1D,B,DDD,DDQ);
       kernels::internal::EvalY(D1D,Q1D,B,DDQ,DQQ);
       kernels::internal::EvalZ(D1D,Q1D,B,DQQ,QQQ);

fem/tmop/tmop_pa_h2d.cpp

@@ -97,9 +97,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, AssembleDiagonalPA_Kernel_2D,
                   const double *Jtr = &J(0,0,qx,qy,e);
 
                   // Jrt = Jtr^{-1}
-                  double j[4];
-                  ConstDeviceMatrix Jrt(j,2,2);
-                  kernels::CalcInverse<2>(Jtr, j);
+                  double jrt_data[4];
+                  ConstDeviceMatrix Jrt(jrt_data,2,2);
+                  kernels::CalcInverse<2>(Jtr, jrt_data);
 
                   const double gg = G(qy,dy) * G(qy,dy);
                   const double gb = G(qy,dy) * B(qy,dy);

fem/tmop/tmop_pa_h2m_c0.cpp

@@ -68,8 +68,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_2D,
             double Xh[2];
             kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ,Xh);
 
-            double B[4];
-            DeviceMatrix H(B,2,2);
+            double H_data[4];
+            DeviceMatrix H(H_data,2,2);
             for (int i = 0; i < DIM; i++)
             {
                for (int j = 0; j < DIM; j++)
@@ -78,9 +78,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_2D,
                }
             }
 
-            // p2 = B . Xh
+            // p2 = H . Xh
             double p2[2];
-            kernels::Mult(2,2,B,Xh,p2);
+            kernels::Mult(2,2,H_data,Xh,p2);
             kernels::internal::PushEval<MQ1,NBZ>(Q1D,qx,qy,p2,QQ);
          }
       }

fem/tmop/tmop_pa_h3m_c0.cpp

@@ -70,8 +70,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_3D,
                double Xh[3];
                kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ,Xh);
 
-               double B[9];
-               DeviceMatrix H(B,3,3);
+               double H_data[9];
+               DeviceMatrix H(H_data,3,3);
                for (int i = 0; i < DIM; i++)
                {
                   for (int j = 0; j < DIM; j++)
@@ -80,9 +80,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_3D,
                   }
                }
 
-               // p2 = B . Xh
+               // p2 = H . Xh
                double p2[3];
-               kernels::Mult(3,3,B,Xh,p2);
+               kernels::Mult(3,3,H_data,Xh,p2);
                kernels::internal::PushEval<MQ1>(Q1D,qx,qy,qz,p2,QQQ);
             }
          }

fem/tmop/tmop_pa_h3s_c0.cpp

@@ -55,6 +55,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
       constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
 
       MFEM_SHARED double sBLD[MQ1*MD1];
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
       ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
 
       MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
@@ -66,8 +67,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
 
       kernels::internal::LoadX(e,D1D,LD,DDD);
 
-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
-
       kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
       kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);
       kernels::internal::EvalZ(D1D,Q1D,BLD,DQQ,QQQ);

fem/tmop/tmop_pa_p3_c0.cpp

@@ -62,6 +62,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
 
       MFEM_SHARED double B[MQ1*MD1];
       MFEM_SHARED double sBLD[MQ1*MD1];
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
       ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
 
       MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
@@ -86,7 +87,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
       kernels::internal::LoadX<MD1>(e,D1D,X1,DDD1);
 
       kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
 
       kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
       kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);

fem/tmop/tmop_pa_w3_c0.cpp

@@ -63,6 +63,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
 
       MFEM_SHARED double B[MQ1*MD1];
       MFEM_SHARED double sBLD[MQ1*MD1];
+      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
       ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
 
       MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
@@ -87,7 +88,6 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
       kernels::internal::LoadX<MD1>(e,D1D,X1,DDD1);
 
       kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
 
       kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
       kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);

fem/tmop_amr.cpp

@@ -837,26 +837,26 @@ void TMOPHRSolver::ParUpdate()
 }
 #endif
 
-void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
+void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh_, NonlinearForm *nlf_)
 {
-   const FiniteElementSpace &fes = *mesh->GetNodalFESpace();
+   const FiniteElementSpace &fes = *mesh_->GetNodalFESpace();
 
    // Update Nonlinear form and Set Essential BC
-   nlf->Update();
+   nlf_->Update();
    const int dim = fes.GetFE(0)->GetDim();
    if (move_bnd == false)
    {
-      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      Array<int> ess_bdr(mesh_->bdr_attributes.Max());
       ess_bdr = 1;
-      nlf->SetEssentialBC(ess_bdr);
+      nlf_->SetEssentialBC(ess_bdr);
    }
    else
    {
       const int nd  = fes.GetBE(0)->GetDof();
       int n = 0;
-      for (int i = 0; i < mesh->GetNBE(); i++)
+      for (int i = 0; i < mesh_->GetNBE(); i++)
       {
-         const int attr = mesh->GetBdrElement(i)->GetAttribute();
+         const int attr = mesh_->GetBdrElement(i)->GetAttribute();
          MFEM_VERIFY(!(dim == 2 && attr == 3),
                      "Boundary attribute 3 must be used only for 3D meshes. "
                      "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
@@ -866,9 +866,9 @@ void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
       }
       Array<int> ess_vdofs(n), vdofs;
       n = 0;
-      for (int i = 0; i < mesh->GetNBE(); i++)
+      for (int i = 0; i < mesh_->GetNBE(); i++)
       {
-         const int attr = mesh->GetBdrElement(i)->GetAttribute();
+         const int attr = mesh_->GetBdrElement(i)->GetAttribute();
          fes.GetBdrElementVDofs(i, vdofs);
          if (attr == 1) // Fix x components.
          {
@@ -891,7 +891,7 @@ void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
             { ess_vdofs[n++] = vdofs[j]; }
          }
       }
-      nlf->SetEssentialVDofs(ess_vdofs);
+      nlf_->SetEssentialVDofs(ess_vdofs);
    }
 }
 

fem/tmop_tools.cpp

@@ -444,15 +444,19 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       if (untangling == false && min_detT_out < 0.0)
       {
          // No untangling, and detJ got negative -- no good.
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Neg det(J) found.\n"; }
+         if (print_options.iterations)
+         {
+            mfem::out << "Scale = " << scale << " Neg det(J) found.\n";
+         }
          scale *= detJ_factor; continue;
       }
       if (untangling == true && min_detT_out < *min_det_ptr)
       {
          // Untangling, and detJ got even more negative -- no good.
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Neg det(J) decreased.\n"; }
+         if (print_options.iterations)
+         {
+            mfem::out << "Scale = " << scale << " Neg det(J) decreased.\n";
+         }
          scale *= detJ_factor; continue;
       }
 
@@ -477,7 +481,7 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       if (energy_out > energy_in + 0.2*fabs(energy_in) ||
           std::isnan(energy_out) != 0)
       {
-         if (print_level >= 0)
+         if (print_options.iterations)
          {
             mfem::out << "Scale = " << scale << " Increasing energy: "
                       << energy_in << " --> " << energy_out << '\n';
@@ -492,7 +496,7 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
 
       if (norm_out > 1.2*norm_in)
       {
-         if (print_level >= 0)
+         if (print_options.iterations)
          {
             mfem::out << "Scale = " << scale << " Norm increased: "
                       << norm_in << " --> " << norm_out << '\n';
@@ -508,13 +512,15 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       if (min_detT_out > 0.0)
       {
          *min_det_ptr = 0.0;
-         if (print_level >= 0)
+         if (print_options.summary || print_options.iterations ||
+             print_options.first_and_last)
          { mfem::out << "The mesh has been untangled at the used points!\n"; }
       }
       else { *min_det_ptr = untangle_factor * min_detT_out; }
    }
 
-   if (print_level >= 0)
+   if (print_options.summary || print_options.iterations ||
+       print_options.first_and_last)
    {
       if (untangling)
       {
@@ -568,9 +574,9 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
    if (parallel)
    {
 #ifdef MFEM_USE_MPI
-      const ParNonlinearForm *nlf =
+      const ParNonlinearForm *pnlf =
          dynamic_cast<const ParNonlinearForm *>(oper);
-      const ParFiniteElementSpace *pfesc = nlf->ParFESpace();
+      const ParFiniteElementSpace *pfesc = pnlf->ParFESpace();
       Vector x_loc(pfesc->GetVSize());
       pfesc->GetProlongationMatrix()->Mult(x, x_loc);
       for (int i = 0; i < integs.Size(); i++)

general/device.cpp

@@ -177,7 +177,7 @@ Device::~Device()
    Get().device_mem_class = MemoryClass::HOST;
 }
 
-void Device::Configure(const std::string &device, const int dev)
+void Device::Configure(const std::string &device, const int device_id)
 {
    // If a device was configured via the environment, skip the configuration,
    // and avoid the 'singleton_device' to destroy the mm.
@@ -240,7 +240,7 @@ void Device::Configure(const std::string &device, const int dev)
 #endif
 
    // Perform setup.
-   Get().Setup(dev);
+   Get().Setup(device_id);
 
    // Enable the device
    Enable();
@@ -276,35 +276,35 @@ void Device::SetMemoryTypes(MemoryType h_mt, MemoryType d_mt)
    // the call mm.Configure(...) in UpdateMemoryTypeAndClass()
 }
 
-void Device::Print(std::ostream &out)
+void Device::Print(std::ostream &os)
 {
-   out << "Device configuration: ";
+   os << "Device configuration: ";
    bool add_comma = false;
    for (int i = 0; i < Backend::NUM_BACKENDS; i++)
    {
       if (backends & internal::backend_list[i])
       {
-         if (add_comma) { out << ','; }
+         if (add_comma) { os << ','; }
          add_comma = true;
-         out << internal::backend_name[i];
+         os << internal::backend_name[i];
       }
    }
-   out << '\n';
+   os << '\n';
 #ifdef MFEM_USE_CEED
    if (Allows(Backend::CEED_MASK))
    {
       const char *ceed_backend;
       CeedGetResource(internal::ceed, &ceed_backend);
-      out << "libCEED backend: " << ceed_backend << '\n';
+      os << "libCEED backend: " << ceed_backend << '\n';
    }
 #endif
-   out << "Memory configuration: "
-       << MemoryTypeName[static_cast<int>(host_mem_type)];
+   os << "Memory configuration: "
+      << MemoryTypeName[static_cast<int>(host_mem_type)];
    if (Device::Allows(Backend::DEVICE_MASK))
    {
-      out << ',' << MemoryTypeName[static_cast<int>(device_mem_type)];
+      os << ',' << MemoryTypeName[static_cast<int>(device_mem_type)];
    }
-   out << std::endl;
+   os << std::endl;
 }
 
 void Device::UpdateMemoryTypeAndClass()
@@ -502,12 +502,12 @@ static void CeedDeviceSetup(const char* ceed_spec)
 #endif
 }
 
-void Device::Setup(const int device)
+void Device::Setup(const int device_id)
 {
    MFEM_VERIFY(ngpu == -1, "the mfem::Device is already configured!");
 
    ngpu = 0;
-   dev = device;
+   dev = device_id;
 #ifndef MFEM_USE_CUDA
    MFEM_VERIFY(!Allows(Backend::CUDA_MASK),
                "the CUDA backends require MFEM built with MFEM_USE_CUDA=YES");

general/device.hpp

@@ -150,7 +150,7 @@ private:
    static Device& Get() { return device_singleton; }
 
    /// Setup switcher based on configuration settings
-   void Setup(const int dev = 0);
+   void Setup(const int device_id = 0);
 
    void MarkBackend(Backend::Id b) { backends |= b; }
 

general/gecko.cpp

@@ -1228,12 +1228,12 @@ Graph::reweight(uint k)
 
 // Linearly order graph.
 void
-Graph::order(Functional* functional, uint iterations, uint window, uint period,
-             uint seed, Progress* progress)
+Graph::order(Functional* functional_, uint iterations, uint window, uint period,
+             uint seed, Progress* progress_)
 {
    // Initialize graph.
-   this->functional = functional;
-   progress = this->progress = progress ? progress : new Progress;
+   this->functional = functional_;
+   progress_ = this->progress = progress_ ? progress_ : new Progress;
    for (level = 0; (1u << level) < nodes(); level++);
    place();
    Float mincost = cost();

general/mem_manager.hpp

@@ -849,13 +849,14 @@ inline void Memory<T>::New(int size, MemoryType mt)
 }
 
 template <typename T>
-inline void Memory<T>::New(int size, MemoryType h_mt, MemoryType d_mt)
+inline void Memory<T>::New(int size, MemoryType host_mt, MemoryType device_mt)
 {
    capacity = size;
    const size_t bytes = size*sizeof(T);
-   this->h_mt = h_mt;
-   T *h_tmp = (h_mt == MemoryType::HOST) ? NewHOST(size) : nullptr;
-   h_ptr = (T*)MemoryManager::New_(h_tmp, bytes, h_mt, d_mt, VALID_HOST, flags);
+   this->h_mt = host_mt;
+   T *h_tmp = (host_mt == MemoryType::HOST) ? NewHOST(size) : nullptr;
+   h_ptr = (T*)MemoryManager::New_(h_tmp, bytes, host_mt, device_mt,
+                                   VALID_HOST, flags);
 }
 
 template <typename T>

general/tic_toc.hpp

@@ -42,7 +42,7 @@ public:
    /// Clear the elapsed time on the stopwatch and restart it if it's running.
    void Clear();
 
-   /// Clear the elapsed time and start the stopwatch.
+   /// Start the stopwatch. The elapsed time is @b not cleared.
    void Start();
 
    /// Stop the stopwatch.

linalg/constraints.cpp

@@ -327,7 +327,7 @@ void EliminationSolver::Mult(const Vector& rhs, Vector& sol) const
    krylov->SetMaxIter(max_iter);
    krylov->SetRelTol(rel_tol);
    krylov->SetAbsTol(abs_tol);
-   krylov->SetPrintLevel(print_level);
+   krylov->SetPrintLevel(print_options);
 
    Vector rtilde(rhs.Size());
    if (constraint_rhs.Size() > 0)
@@ -447,7 +447,7 @@ void PenaltyConstrainedSolver::Mult(const Vector& b, Vector& x) const
    krylov->SetRelTol(rel_tol);
    krylov->SetAbsTol(abs_tol);
    krylov->SetMaxIter(max_iter);
-   krylov->SetPrintLevel(print_level);
+   krylov->SetPrintLevel(print_options);
    krylov->Mult(penalized_rhs, x);
    final_iter = krylov->GetNumIterations();
    final_norm = krylov->GetFinalNorm();
@@ -574,7 +574,7 @@ void SchurConstrainedSolver::LagrangeSystemMult(const Vector& x,
    gmres->SetRelTol(rel_tol);
    gmres->SetAbsTol(abs_tol);
    gmres->SetMaxIter(max_iter);
-   gmres->SetPrintLevel(print_level);
+   gmres->SetPrintLevel(print_options);
    gmres->SetPreconditioner(
       const_cast<BlockDiagonalPreconditioner&>(*block_pc));
 

linalg/operator.hpp

@@ -701,11 +701,15 @@ private:
 public:
    /// Create an operator which is a scalar multiple of A.
    explicit ScaledOperator(const Operator *A, double a)
-      : Operator(A->Width(), A->Height()), A_(*A), a_(a) { }
+      : Operator(A->Height(), A->Width()), A_(*A), a_(a) { }
 
    /// Operator application
    virtual void Mult(const Vector &x, Vector &y) const
    { A_.Mult(x, y); y *= a_; }
+
+   /// Application of the transpose.
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   { A_.MultTranspose(x, y); y *= a_; }
 };
 
 

linalg/simd/auto.hpp

@@ -29,6 +29,10 @@ struct MFEM_AUTOSIMD_ALIGN__ AutoSIMD
 
    scalar_t vec[size];
 
+   AutoSIMD() = default;
+
+   AutoSIMD(const AutoSIMD &) = default;
+
    inline MFEM_ALWAYS_INLINE scalar_t &operator[](int i)
    {
       return vec[i];
@@ -117,6 +121,11 @@ struct MFEM_AUTOSIMD_ALIGN__ AutoSIMD
       return r;
    }
 
+   inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
+   {
+      return *this;
+   }
+
    inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
    {
       AutoSIMD r;

linalg/simd/m128.hpp

@@ -38,6 +38,10 @@ template <> struct AutoSIMD<double,2,16>
       double vec[size];
    };
 
+   AutoSIMD() = default;
+
+   AutoSIMD(const AutoSIMD &) = default;
+
    inline MFEM_ALWAYS_INLINE double &operator[](int i)
    {
       return vec[i];
@@ -115,6 +119,11 @@ template <> struct AutoSIMD<double,2,16>
       return r;
    }
 
+   inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
+   {
+      return *this;
+   }
+
    inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
    {
       AutoSIMD r;
@@ -251,4 +260,3 @@ AutoSIMD<double,2,16> operator/(const double &e,
 #endif // __SSE2__
 
 #endif // MFEM_SIMD_M128_HPP
-

linalg/simd/m256.hpp

@@ -38,6 +38,10 @@ template <> struct AutoSIMD<double,4,32>
       double vec[size];
    };
 
+   AutoSIMD() = default;
+
+   AutoSIMD(const AutoSIMD &) = default;
+
    inline MFEM_ALWAYS_INLINE double &operator[](int i)
    {
       return vec[i];
@@ -115,6 +119,11 @@ template <> struct AutoSIMD<double,4,32>
       return r;
    }
 
+   inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
+   {
+      return *this;
+   }
+
    inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
    {
       AutoSIMD r;

linalg/simd/m512.hpp

@@ -39,6 +39,10 @@ template <> struct AutoSIMD<double,8,64>
       double vec[size];
    };
 
+   AutoSIMD() = default;
+
+   AutoSIMD(const AutoSIMD &) = default;
+
    inline MFEM_ALWAYS_INLINE double &operator[](int i)
    {
       return vec[i];
@@ -120,6 +124,11 @@ template <> struct AutoSIMD<double,8,64>
       return r;
    }
 
+   inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
+   {
+      return *this;
+   }
+
    inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
    {
       AutoSIMD r;

linalg/simd/qpx256.hpp

@@ -34,6 +34,10 @@ template <> struct AutoSIMD<double,4,32>
       double vec[size];
    };
 
+   AutoSIMD() = default;
+
+   AutoSIMD(const AutoSIMD &) = default;
+
    inline __ATTRS_ai double &operator[](int i) { return vec[i]; }
 
    inline __ATTRS_ai const double &operator[](int i) const { return vec[i]; }
@@ -105,6 +109,11 @@ template <> struct AutoSIMD<double,4,32>
       return r;
    }
 
+   inline __ATTRS_ai AutoSIMD operator+() const
+   {
+      return *this;
+   }
+
    inline __ATTRS_ai AutoSIMD operator+(const AutoSIMD &v) const
    {
       AutoSIMD r;

linalg/simd/sve.hpp

@@ -33,6 +33,10 @@ template <> struct MFEM_AUTOSIMD_ALIGN_SVE AutoSIMD<double,8,64>
 
    double vec[size];
 
+   AutoSIMD() = default;
+
+   AutoSIMD(const AutoSIMD &) = default;
+
    inline MFEM_ALWAYS_INLINE double &operator[](int i)
    {
       return vec[i];
@@ -123,6 +127,11 @@ template <> struct MFEM_AUTOSIMD_ALIGN_SVE AutoSIMD<double,8,64>
       return r;
    }
 
+   inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
+   {
+      return *this;
+   }
+
    inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
    {
       AutoSIMD r;

linalg/simd/vsx128.hpp

@@ -17,6 +17,10 @@
 #include "../../config/tconfig.hpp"
 #include <altivec.h>
 
+#ifdef __GNUC__
+#undef bool
+#endif
+
 namespace mfem
 {
 
@@ -34,6 +38,10 @@ template <> struct AutoSIMD<double,2,16>
       double vec[size];
    };
 
+   AutoSIMD() = default;
+
+   AutoSIMD(const AutoSIMD &) = default;
+
    inline MFEM_ALWAYS_INLINE double &operator[](int i)
    {
       return vec[i];
@@ -107,10 +115,19 @@ template <> struct AutoSIMD<double,2,16>
    inline MFEM_ALWAYS_INLINE AutoSIMD operator-() const
    {
       AutoSIMD r;
+#ifndef __GNUC__
       r.vd = vec_neg(vd);
+#else
+      r.vd = vec_splats(0.0) - vd;
+#endif
       return r;
    }
 
+   inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
+   {
+      return *this;
+   }
+
    inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
    {
       AutoSIMD r;

linalg/solvers.hpp

@@ -56,7 +56,7 @@ public:
    {
    }
 
-   /** @brief This method is invoked by ItertiveSolver::SetMonitor, informing
+   /** @brief This method is invoked by IterativeSolver::SetMonitor, informing
        the monitor which IterativeSolver is using it. */
    void SetIterativeSolver(const IterativeSolver &solver)
    { iter_solver = &solver; }
@@ -65,10 +65,55 @@ public:
 /// Abstract base class for iterative solver
 class IterativeSolver : public Solver
 {
+public:
+   /** @brief Settings for the output behavior of the IterativeSolver.
+
+       By default, all output is suppressed. The construction of the desired
+       print level can be achieved through a builder pattern, for example
+
+           PrintLevel().Errors().Warnings()
+
+       constructs the print level with only errors and warnings enabled.
+     */
+   struct PrintLevel
+   {
+      /** @brief If a fatal problem has been detected the failure will be
+          reported to @ref mfem::err. */
+      bool errors = false;
+      /** @brief If a non-fatal problem has been detected some context-specific
+          information will be reported to @ref mfem::out */
+      bool warnings = false;
+      /** @brief Detailed information about each iteration will be reported to
+          @ref mfem::out */
+      bool iterations = false;
+      /** @brief A summary of the solver process will be reported after the last
+          iteration to @ref mfem::out */
+      bool summary = false;
+      /** @brief Information about the first and last iteration will be printed
+          to @ref mfem::out */
+      bool first_and_last = false;
+
+      /// Initializes the print level to suppress
+      PrintLevel() = default;
+
+      /** @name Builder
+         These methods are utilized to construct PrintLevel objects through a
+         builder approach by chaining the function calls in this group. */
+      ///@{
+      PrintLevel &None() { *this = PrintLevel(); return *this; }
+      PrintLevel &Warnings() { warnings=true; return *this; }
+      PrintLevel &Errors() { errors=true; return *this; }
+      PrintLevel &Iterations() { iterations=true; return *this; }
+      PrintLevel &FirstAndLast() { first_and_last=true; return *this; }
+      PrintLevel &Summary() { summary=true; return *this; }
+      PrintLevel &All() { return Warnings().Errors().Iterations().FirstAndLast().Summary(); }
+      ///@}
+   };
+
 #ifdef MFEM_USE_MPI
 private:
    int dot_prod_type; // 0 - local, 1 - global over 'comm'
-   MPI_Comm comm;
+   MPI_Comm comm = MPI_COMM_NULL;
 #endif
 
 protected:
@@ -76,13 +121,52 @@ protected:
    Solver *prec;
    IterativeSolverMonitor *monitor = nullptr;
 
-   int max_iter, print_level;
-   double rel_tol, abs_tol;
+   /// @name Reporting (protected attributes and member functions)
+   ///@{
 
-   // stats
-   mutable int final_iter, converged;
+   /** @brief (DEPRECATED) Legacy print level definition, which is left for
+       compatibility with custom iterative solvers.
+       @deprecated #print_options should be used instead. */
+   int print_level = -1;
+
+   /** @brief Output behavior for the iterative solver.
+
+       This primarily controls the output behavior of the iterative solvers
+       provided by this library. This member must be synchronized with
+       #print_level to ensure compatibility with custom iterative solvers. */
+   PrintLevel print_options;
+
+   /// Convert a legacy print level integer to a PrintLevel object
+   PrintLevel FromLegacyPrintLevel(int);
+
+   /// @brief Use some heuristics to guess a legacy print level corresponding to
+   /// the given PrintLevel.
+   static int GuessLegacyPrintLevel(PrintLevel);
+   ///@}
+
+   /// @name Convergence (protected attributes)
+   ///@{
+
+   /// Limit for the number of iterations the solver is allowed to do
+   int max_iter;
+
+   /// Relative tolerance.
+   double rel_tol;
+
+   /// Absolute tolerance.
+   double abs_tol;
+
+   ///@}
+
+   /// @name Solver statistics (protected attributes)
+   ///@{
+
+   mutable int final_iter;
+   mutable bool converged;
    mutable double final_norm;
 
+   ///@}
+
    double Dot(const Vector &x, const Vector &y) const;
    double Norm(const Vector &x) const { return sqrt(Dot(x, x)); }
    void Monitor(int it, double norm, const Vector& r, const Vector& x,
@@ -95,20 +179,80 @@ public:
    IterativeSolver(MPI_Comm comm_);
 #endif
 
+   /** @name Convergence
+       @brief Termination criteria for the iterative solvers.
+
+       @details While the convergence criterion is solver specific, most of the
+       provided iterative solvers use one of the following criteria
+
+       \f$ ||r||_X \leq tol_{rel}||r_0||_X \f$,
+
+       \f$ ||r||_X \leq tol_{abs} \f$,
+
+       \f$ ||r||_X \leq \max\{ tol_{abs}, tol_{rel} ||r_0||_X \} \f$,
+
+       where X denotes the space in which the norm is measured. The choice of
+       X depends on the specific iterative solver.
+      */
+   ///@{
    void SetRelTol(double rtol) { rel_tol = rtol; }
    void SetAbsTol(double atol) { abs_tol = atol; }
    void SetMaxIter(int max_it) { max_iter = max_it; }
-   void SetPrintLevel(int print_lvl);
+   ///@}
 
+   /** @name Reporting
+       These options control the internal reporting behavior into ::mfem::out
+       and ::mfem::err of the iterative solvers.
+    */
+   ///@{
+
+   /// @brief Legacy method to set the level of verbosity of the solver output.
+   /** This is the old way to control what information will be printed to
+       ::mfem::out and ::mfem::err. The behavior for the print level for all
+       iterative solvers is:
+
+       - -1: Suppress all outputs.
+       -  0: Print information about all detected issues (e.g. no convergence).
+       -  1: Same as level 0, but with detailed information about each
+             iteration.
+       -  2: Print detected issues and a summary when the solver terminates.
+       -  3: Same as 2, but print also the first and last iterations.
+       - >3: Custom print options which are dependent on the specific solver.
+
+       In parallel, only rank 0 produces output.
+
+       @note It is recommended to use @ref SetPrintLevel(PrintLevel) instead.
+
+       @note Some derived classes, like KINSolver, redefine this method and use
+       their own set of print level constants. */
+   virtual void SetPrintLevel(int print_lvl);
+
+   /// @brief Set the level of verbosity of the solver output.
+   /** In parallel, only rank 0 produces outputs. Errors are output to
+       ::mfem::err and all other information to ::mfem::out.
+
+       @note Not all subclasses of IterativeSolver support all possible options.
+
+       @note Some derived classes, like KINSolver, disable this method in favor
+       of SetPrintLevel(int).
+
+       @sa PrintLevel for possible options.
+   */
+   virtual void SetPrintLevel(PrintLevel);
+   ///@}
+
+   /// @name Solver statistics
+   ///@{
    int GetNumIterations() const { return final_iter; }
-   int GetConverged() const { return converged; }
+   bool GetConverged() const { return converged; }
    double GetFinalNorm() const { return final_norm; }
+   ///@}
 
    /// This should be called before SetOperator
    virtual void SetPreconditioner(Solver &pr);
 
    /// Also calls SetOperator for the preconditioner if there is one
-   virtual void SetOperator(const Operator &op);
+   virtual void SetOperator(const Operator &op) override;
 
    /// Set the iterative solver monitor
    void SetMonitor(IterativeSolverMonitor &m)
@@ -163,6 +307,9 @@ public:
 
    ~OperatorJacobiSmoother() {}
 
+   /// Replace diagonal entries with their absolute values.
+   void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
+
    void Mult(const Vector &x, Vector &y) const;
    void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
 
@@ -184,6 +331,8 @@ private:
    const double damping;
    const Array<int> *ess_tdof_list; // not owned; may be NULL
    mutable Vector residual;
+   /// Uses absolute values of the diagonal entries.
+   bool use_abs_diag = false;
 
    const Operator *oper; // not owned
 

linalg/sparsemat.cpp

@@ -2382,7 +2382,7 @@ double SparseMatrix::GetJacobiScaling() const
 }
 
 void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
-                          double sc) const
+                          double sc, bool use_abs_diag) const
 {
    MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
 
@@ -2403,7 +2403,8 @@ void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
       }
       if (d >= 0 && A[d] != 0.0)
       {
-         x1(i) = sc * (sum / A[d]) + (1.0 - sc) * x0(i);
+         const double diag = (use_abs_diag) ? fabs(A[d]) : A[d];
+         x1(i) = sc * (sum / diag) + (1.0 - sc) * x0(i);
       }
       else
       {
@@ -2412,7 +2413,8 @@ void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
    }
 }
 
-void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
+void SparseMatrix::DiagScale(const Vector &b, Vector &x,
+                             double sc, bool use_abs_diag) const
 {
    MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
 
@@ -2438,11 +2440,12 @@ void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
          }
          if (Jp[j] == i)
          {
-            if (!(std::abs(Ap[j]) > 0.0))
+            const double diag = (use_abs_diag) ? fabs(Ap[j]) : Ap[j];
+            if (diag == 0.0)
             {
                MFEM_ABORT_KERNEL("Zero diagonal in SparseMatrix::DiagScale");
             }
-            xp[i] = sc * bp[i] / Ap[j];
+            xp[i] = sc * bp[i] / diag;
             break;
          }
       }

linalg/sparsemat.hpp

@@ -446,10 +446,14 @@ public:
    /// Determine appropriate scaling for Jacobi iteration
    double GetJacobiScaling() const;
    /** One scaled Jacobi iteration for the system A x = b.
-       x1 = x0 + sc D^{-1} (b - A x0)  where D is the diag of A. */
-   void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const;
+       x1 = x0 + sc D^{-1} (b - A x0)  where D is the diag of A.
+       Absolute values of D are used when use_abs_diag = true. */
+   void Jacobi(const Vector &b, const Vector &x0, Vector &x1,
+               double sc, bool use_abs_diag = false) const;
 
-   void DiagScale(const Vector &b, Vector &x, double sc = 1.0) const;
+   /// x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
+   void DiagScale(const Vector &b, Vector &x,
+                  double sc = 1.0, bool use_abs_diag = false) const;
 
    /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j |A_{ij}| \f$. */
    void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,

linalg/sparsesmoothers.cpp

@@ -65,7 +65,7 @@ void DSmoother::Mult(const Vector &x, Vector &y) const
 {
    if (!iterative_mode && type == 0 && iterations == 1)
    {
-      oper->DiagScale(x, y, scale);
+      oper->DiagScale(x, y, scale, use_abs_diag);
       return;
    }
 
@@ -90,7 +90,7 @@ void DSmoother::Mult(const Vector &x, Vector &y) const
    {
       if (type == 0)
       {
-         oper->Jacobi(x, *p, *r, scale);
+         oper->Jacobi(x, *p, *r, scale, use_abs_diag);
       }
       else if (type == 1)
       {

linalg/sparsesmoothers.hpp

@@ -58,6 +58,8 @@ protected:
    int type; // 0, 1, 2 - scaled Jacobi, scaled l1-Jacobi, scaled lumped-Jacobi
    double scale;
    int iterations;
+   /// Uses abs values of the diagonal entries. Relevant only when type = 0.
+   bool use_abs_diag = false;
 
    mutable Vector z;
 
@@ -69,6 +71,9 @@ public:
    /// Create Jacobi smoother.
    DSmoother(const SparseMatrix &a, int t = 0, double s = 1., int it = 1);
 
+   /// Replace diag entries with their abs values. Relevant only when type = 0.
+   void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
+
    /// Matrix vector multiplication with Jacobi smoother.
    virtual void Mult(const Vector &x, Vector &y) const;
 };

linalg/sundials.cpp

@@ -1976,6 +1976,11 @@ void KINSolver::SetMAA(int m_aa)
    }
 }
 
+void KINSolver::SetPrintLevel(PrintLevel)
+{
+   MFEM_ABORT("this method is not supported! Use SetPrintLevel(int) instead.");
+}
+
 // Compute the scaling vectors and solve nonlinear system
 void KINSolver::Mult(const Vector&, Vector &x) const
 {

linalg/sundials.hpp

@@ -833,6 +833,12 @@ public:
    /** @note Only valid in combination with JFNK */
    void SetLSMaxRestarts(int m) { maxlrs = m; }
 
+   /// Set the print level for the KINSetPrintLevel function.
+   virtual void SetPrintLevel(int print_lvl) { print_level = print_lvl; }
+
+   /// This method is not supported and will throw an error.
+   virtual void SetPrintLevel(PrintLevel);
+
    /// Solve the nonlinear system \f$ F(x) = 0 \f$.
    /** This method computes the x_scale and fx_scale vectors and calls the
        other Mult(Vector&, Vector&, Vector&) const method. The x_scale vector

mesh/mesh.cpp

@@ -350,6 +350,7 @@ void Mesh::GetElementTransformation(int i, IsoparametricTransformation *ElTr)
    ElTr->Attribute = GetAttribute(i);
    ElTr->ElementNo = i;
    ElTr->ElementType = ElementTransformation::ELEMENT;
+   ElTr->mesh = this;
    ElTr->Reset();
    if (Nodes == NULL)
    {
@@ -382,6 +383,7 @@ void Mesh::GetElementTransformation(int i, const Vector &nodes,
    ElTr->Attribute = GetAttribute(i);
    ElTr->ElementNo = i;
    ElTr->ElementType = ElementTransformation::ELEMENT;
+   ElTr->mesh = this;
    DenseMatrix &pm = ElTr->GetPointMat();
    ElTr->Reset();
    nodes.HostRead();
@@ -437,6 +439,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
    ElTr->Attribute = GetBdrAttribute(i);
    ElTr->ElementNo = i; // boundary element number
    ElTr->ElementType = ElementTransformation::BDR_ELEMENT;
+   ElTr->mesh = this;
    DenseMatrix &pm = ElTr->GetPointMat();
    ElTr->Reset();
    if (Nodes == NULL)
@@ -481,6 +484,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
                      "Mesh requires nodal Finite Element.");
          IntegrationRule eir(face_el->GetDof());
          FaceElemTr.Loc1.Transf.ElementNo = elem_id;
+         FaceElemTr.Loc1.Transf.mesh = this;
          FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
          FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
          Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
@@ -495,6 +499,7 @@ void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
    FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
    FTr->ElementNo = FaceNo;
    FTr->ElementType = ElementTransformation::FACE;
+   FTr->mesh = this;
    DenseMatrix &pm = FTr->GetPointMat();
    FTr->Reset();
    if (Nodes == NULL)
@@ -551,6 +556,7 @@ void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
          IntegrationRule eir(face_el->GetDof());
          FaceElemTr.Loc1.Transf.ElementNo = face_info.Elem1No;
          FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
+         FaceElemTr.Loc1.Transf.mesh = this;
          FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
          Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
 
@@ -580,6 +586,7 @@ void Mesh::GetEdgeTransformation(int EdgeNo, IsoparametricTransformation *EdTr)
    EdTr->Attribute = 1;
    EdTr->ElementNo = EdgeNo;
    EdTr->ElementType = ElementTransformation::EDGE;
+   EdTr->mesh = this;
    DenseMatrix &pm = EdTr->GetPointMat();
    EdTr->Reset();
    if (Nodes == NULL)
@@ -1098,6 +1105,7 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
    tr->Attribute = boundary[BdrElemNo]->GetAttribute();
    tr->ElementNo = BdrElemNo;
    tr->ElementType = ElementTransformation::BDR_FACE;
+   tr->mesh = this;
    return tr;
 }
 

mesh/ncmesh.hpp

@@ -737,6 +737,8 @@ protected: // implementation
 
       Point() { dim = 0; }
 
+      Point(const Point &) = default;
+
       Point(double x)
       { dim = 1; coord[0] = x; }
 

mesh/pmesh.cpp

@@ -1890,6 +1890,7 @@ void ParMesh::GetFaceNbrElementTransformation(
    ElTr->Attribute = elem->GetAttribute();
    ElTr->ElementNo = NumOfElements + i;
    ElTr->ElementType = ElementTransformation::ELEMENT;
+   ElTr->mesh = this;
    ElTr->Reset();
 
    if (Nodes == NULL)
@@ -2711,6 +2712,7 @@ STable3D *ParMesh::GetFaceNbrElementToFaceTable(int ret_ftbl)
    }
    face_nbr_el_to_face->Finalize();
 
+   delete sfaces_tbl;
    if (ret_ftbl)
    {
       return faces_tbl;
@@ -6178,6 +6180,9 @@ void ParMesh::Destroy()
       FreeElement(shared_edges[i]);
    }
    shared_edges.DeleteAll();
+
+   delete face_nbr_el_to_face;
+   face_nbr_el_to_face = NULL;
 }
 
 ParMesh::~ParMesh()

miniapps/autodiff/example.hpp

@@ -804,11 +804,11 @@ public:
          mfem::ad::ADFloatType aa;
          mfem::ad::ADVectorType lvec(ndof);
 
-         for (int i = 0; i < ir.GetNPoints(); i++)
+         for (int q = 0; q < ir.GetNPoints(); q++)
          {
             lvec=0.0;
 
-            const IntegrationPoint &ip = ir.IntPoint(i);
+            const IntegrationPoint &ip = ir.IntPoint(q);
             trans.SetIntPoint(&ip);
             w = trans.Weight();
             detJ = (square ? w : w * w);

miniapps/autodiff/par_example.cpp

@@ -485,7 +485,7 @@ int main(int argc, char *argv[])
                    << std::endl;
       }
       // Compute the energy
-      double energy = nr->GetEnergy(*sv);
+      energy = nr->GetEnergy(*sv);
       if (myrank==0)
       {
          std::cout << "[pp="<<i<<"] The total energy of the system is E=" << energy
@@ -517,7 +517,7 @@ int main(int argc, char *argv[])
                    << std::endl;
       }
       // Compute the energy
-      double energy = nr->GetEnergy(*sv);
+      energy = nr->GetEnergy(*sv);
       if (myrank==0)
       {
          std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy

miniapps/electromagnetics/joule_solver.cpp

@@ -678,22 +678,22 @@ void MagneticDiffusionEOperator::buildA1(double muInv,
    dt_A1 = dt;
 }
 
-void MagneticDiffusionEOperator::buildA2(MeshDependentCoefficient &InvTcond,
-                                         MeshDependentCoefficient &InvTcap,
-                                         double dt)
+void MagneticDiffusionEOperator::buildA2(MeshDependentCoefficient &InvTcond_,
+                                         MeshDependentCoefficient &InvTcap_,
+                                         double dt_)
 {
    if ( a2 != NULL ) { delete a2; }
 
-   InvTcap.SetScaleFactor(dt);
+   InvTcap_.SetScaleFactor(dt_);
    a2 = new ParBilinearForm(&HDivFESpace);
-   a2->AddDomainIntegrator(new VectorFEMassIntegrator(InvTcond));
-   a2->AddDomainIntegrator(new DivDivIntegrator(InvTcap));
+   a2->AddDomainIntegrator(new VectorFEMassIntegrator(InvTcond_));
+   a2->AddDomainIntegrator(new DivDivIntegrator(InvTcap_));
    if (STATIC_COND == 1) { a2->EnableStaticCondensation(); }
    a2->Assemble();
 
    // Don't finalize or parallel assemble this is done in FormLinearSystem.
 
-   dt_A2 = dt;
+   dt_A2 = dt_;
 }
 
 void MagneticDiffusionEOperator::buildM1(MeshDependentCoefficient &Sigma)
@@ -719,13 +719,13 @@ void MagneticDiffusionEOperator::buildM2(MeshDependentCoefficient &Alpha)
    // Don't finalize or parallel assemble this is done in FormLinearSystem.
 }
 
-void MagneticDiffusionEOperator::buildM3(MeshDependentCoefficient &Tcapacity)
+void MagneticDiffusionEOperator::buildM3(MeshDependentCoefficient &Tcapacity_)
 {
    if ( m3 != NULL ) { delete m3; }
 
    // ConstantCoefficient Sigma(sigma);
    m3 = new ParBilinearForm(&L2FESpace);
-   m3->AddDomainIntegrator(new MassIntegrator(Tcapacity));
+   m3->AddDomainIntegrator(new MassIntegrator(Tcapacity_));
    m3->Assemble();
    m3->Finalize();
    M3 = m3->ParallelAssemble();
@@ -741,13 +741,13 @@ void MagneticDiffusionEOperator::buildS1(double muInv)
    s1->Assemble();
 }
 
-void MagneticDiffusionEOperator::buildS2(MeshDependentCoefficient &InvTcap)
+void MagneticDiffusionEOperator::buildS2(MeshDependentCoefficient &InvTcap_)
 {
    if ( s2 != NULL ) { delete s2; }
 
    // ConstantCoefficient param(a);
    s2 = new ParBilinearForm(&HDivFESpace);
-   s2->AddDomainIntegrator(new DivDivIntegrator(InvTcap));
+   s2->AddDomainIntegrator(new DivDivIntegrator(InvTcap_));
    s2->Assemble();
 }
 
@@ -768,13 +768,13 @@ void MagneticDiffusionEOperator::buildCurl(double muInv)
    // no ParallelAssemble since this will be applied to GridFunctions
 }
 
-void MagneticDiffusionEOperator::buildDiv(MeshDependentCoefficient &InvTcap)
+void MagneticDiffusionEOperator::buildDiv(MeshDependentCoefficient &InvTcap_)
 {
    if ( weakDiv != NULL ) { delete weakDiv; }
    if ( weakDivC != NULL ) { delete weakDivC; }
 
    weakDivC = new ParMixedBilinearForm(&HDivFESpace, &L2FESpace);
-   weakDivC->AddDomainIntegrator(new VectorFEDivergenceIntegrator(InvTcap));
+   weakDivC->AddDomainIntegrator(new VectorFEDivergenceIntegrator(InvTcap_));
    weakDivC->Assemble();
 
    weakDiv = new ParMixedBilinearForm(&HDivFESpace, &L2FESpace);

miniapps/meshing/mesh-explorer.cpp

@@ -1143,8 +1143,8 @@ int main (int argc, char *argv[])
 
       if (mk == 'S')
       {
-         const char mesh_file[] = "mesh-explorer.mesh";
-         ofstream omesh(mesh_file);
+         const char omesh_file[] = "mesh-explorer.mesh";
+         ofstream omesh(omesh_file);
          omesh.precision(14);
          mesh->Print(omesh);
          cout << "New mesh file: " << mesh_file << endl;
@@ -1152,8 +1152,8 @@ int main (int argc, char *argv[])
 
       if (mk == 'V')
       {
-         const char mesh_file[] = "mesh-explorer.vtk";
-         ofstream omesh(mesh_file);
+         const char omesh_file[] = "mesh-explorer.vtk";
+         ofstream omesh(omesh_file);
          omesh.precision(14);
          mesh->PrintVTK(omesh);
          cout << "New VTK mesh file: " << mesh_file << endl;

miniapps/meshing/mesh-optimizer.cpp

@@ -326,7 +326,7 @@ int main(int argc, char *argv[])
    //    In addition, compute average mesh size and total volume.
    Vector h0(fespace->GetNDofs());
    h0 = infinity();
-   double volume = 0.0;
+   double mesh_volume = 0.0;
    Array<int> dofs;
    for (int i = 0; i < mesh->GetNE(); i++)
    {
@@ -338,9 +338,9 @@ int main(int argc, char *argv[])
       {
          h0(dofs[j]) = min(h0(dofs[j]), hi);
       }
-      volume += mesh->GetElementVolume(i);
+      mesh_volume += mesh->GetElementVolume(i);
    }
-   const double small_phys_size = pow(volume, 1.0 / dim) / 100.0;
+   const double small_phys_size = pow(mesh_volume, 1.0 / dim) / 100.0;
 
    // 8. Add a random perturbation to the nodes in the interior of the domain.
    //    We define a random grid function of fespace and make sure that it is
@@ -413,6 +413,10 @@ int main(int argc, char *argv[])
       case 315: metric = new TMOP_Metric_315; break;
       case 316: metric = new TMOP_Metric_316; break;
       case 321: metric = new TMOP_Metric_321; break;
+      case 328: metric = new TMOP_Metric_328(0.5); break;
+      case 332: metric = new TMOP_Metric_332(0.5); break;
+      case 333: metric = new TMOP_Metric_333(0.5); break;
+      case 334: metric = new TMOP_Metric_334(0.5); break;
       // case 352: metric = new TMOP_Metric_352(tauval); break;
       // A-metrics
       case 11: metric = new TMOP_AMetric_011; break;
@@ -461,7 +465,7 @@ int main(int argc, char *argv[])
    H1_FECollection ind_fec(mesh_poly_deg, dim);
    FiniteElementSpace ind_fes(mesh, &ind_fec);
    FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
-   GridFunction size(&ind_fes), aspr(&ind_fes), disc(&ind_fes), ori(&ind_fes);
+   GridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
    GridFunction aspr3d(&ind_fesv);
 
    const AssemblyLevel al =
@@ -499,13 +503,13 @@ int main(int argc, char *argv[])
          }
          if (dim == 2)
          {
-            FunctionCoefficient ind_coeff(discrete_size_2d);
-            size.ProjectCoefficient(ind_coeff);
+            FunctionCoefficient size_coeff(discrete_size_2d);
+            size.ProjectCoefficient(size_coeff);
          }
          else if (dim == 3)
          {
-            FunctionCoefficient ind_coeff(discrete_size_3d);
-            size.ProjectCoefficient(ind_coeff);
+            FunctionCoefficient size_coeff(discrete_size_3d);
+            size.ProjectCoefficient(size_coeff);
          }
          tc->SetSerialDiscreteTargetSize(size);
          target_c = tc;
@@ -513,12 +517,12 @@ int main(int argc, char *argv[])
       }
       case 6: // Discrete size + aspect ratio - 2D
       {
-         GridFunction d_x(&ind_fes), d_y(&ind_fes);
+         GridFunction d_x(&ind_fes), d_y(&ind_fes), disc(&ind_fes);
 
          target_t = TargetConstructor::GIVEN_SHAPE_AND_SIZE;
          DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
-         FunctionCoefficient ind_coeff(material_indicator_2d);
-         disc.ProjectCoefficient(ind_coeff);
+         FunctionCoefficient mat_coeff(material_indicator_2d);
+         disc.ProjectCoefficient(mat_coeff);
          if (adapt_eval == 0)
          {
             tc->SetAdaptivityEvaluator(new AdvectorCG(al));
@@ -648,8 +652,8 @@ int main(int argc, char *argv[])
 
          if (metric_id == 14 || metric_id == 36)
          {
-            ConstantCoefficient ind_coeff(0.1*0.1);
-            size.ProjectCoefficient(ind_coeff);
+            ConstantCoefficient size_coeff(0.1*0.1);
+            size.ProjectCoefficient(size_coeff);
             tc->SetSerialDiscreteTargetSize(size);
          }
 
@@ -686,16 +690,16 @@ int main(int argc, char *argv[])
       target_c = new TargetConstructor(target_t);
    }
    target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c,
-                                                       h_metric);
+   TMOP_Integrator *tmop_integ = new TMOP_Integrator(metric, target_c,
+                                                     h_metric);
 
    // Finite differences for computations of derivatives.
    if (fdscheme)
    {
       MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
-      he_nlf_integ->EnableFiniteDifferences(x);
+      tmop_integ->EnableFiniteDifferences(x);
    }
-   he_nlf_integ->SetExactActionFlag(exactaction);
+   tmop_integ->SetExactActionFlag(exactaction);
 
    // Setup the quadrature rules for the TMOP integrator.
    IntegrationRules *irules = NULL;
@@ -706,7 +710,7 @@ int main(int argc, char *argv[])
       case 3: irules = &IntRulesCU; break;
       default: cout << "Unknown quad_type: " << quad_type << endl; return 3;
    }
-   he_nlf_integ->SetIntegrationRules(*irules, quad_order);
+   tmop_integ->SetIntegrationRules(*irules, quad_order);
    if (dim == 2)
    {
       cout << "Triangle quadrature points: "
@@ -732,49 +736,50 @@ int main(int argc, char *argv[])
    // The small_phys_size is relevant only with proper normalization.
    if (normalization) { dist = small_phys_size; }
    ConstantCoefficient lim_coeff(lim_const);
-   if (lim_const != 0.0) { he_nlf_integ->EnableLimiting(x0, dist, lim_coeff); }
+   if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
 
    // Adaptive limiting.
-   GridFunction zeta_0(&ind_fes);
-   ConstantCoefficient coef_zeta(adapt_lim_const);
-   AdaptivityEvaluator *adapt_evaluator = NULL;
+   GridFunction adapt_lim_gf0(&ind_fes);
+   ConstantCoefficient adapt_lim_coeff(adapt_lim_const);
+   AdaptivityEvaluator *adapt_lim_eval = NULL;
    if (adapt_lim_const > 0.0)
    {
       MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
 
-      FunctionCoefficient alim_coeff(adapt_lim_fun);
-      zeta_0.ProjectCoefficient(alim_coeff);
+      FunctionCoefficient adapt_lim_gf0_coeff(adapt_lim_fun);
+      adapt_lim_gf0.ProjectCoefficient(adapt_lim_gf0_coeff);
 
-      if (adapt_eval == 0) { adapt_evaluator = new AdvectorCG(al); }
+      if (adapt_eval == 0) { adapt_lim_eval = new AdvectorCG(al); }
       else if (adapt_eval == 1)
       {
 #ifdef MFEM_USE_GSLIB
-         adapt_evaluator = new InterpolatorFP;
+         adapt_lim_eval = new InterpolatorFP;
 #else
          MFEM_ABORT("MFEM is not built with GSLIB support!");
 #endif
       }
       else { MFEM_ABORT("Bad interpolation option."); }
 
-      he_nlf_integ->EnableAdaptiveLimiting(zeta_0, coef_zeta, *adapt_evaluator);
+      tmop_integ->EnableAdaptiveLimiting(adapt_lim_gf0, adapt_lim_coeff,
+                                         *adapt_lim_eval);
       if (visualization)
       {
          socketstream vis1;
-         common::VisualizeField(vis1, "localhost", 19916, zeta_0, "Zeta 0",
+         common::VisualizeField(vis1, "localhost", 19916, adapt_lim_gf0, "Zeta 0",
                                 300, 600, 300, 300);
       }
    }
 
    // Surface fitting.
    L2_FECollection mat_coll(0, dim);
-   H1_FECollection sigma_fec(mesh_poly_deg, dim);
-   FiniteElementSpace sigma_fes(mesh, &sigma_fec);
+   H1_FECollection surf_fit_fec(mesh_poly_deg, dim);
+   FiniteElementSpace surf_fit_fes(mesh, &surf_fit_fec);
    FiniteElementSpace mat_fes(mesh, &mat_coll);
    GridFunction mat(&mat_fes);
-   GridFunction marker_gf(&sigma_fes);
-   GridFunction ls_0(&sigma_fes);
-   Array<bool> marker(ls_0.Size());
-   ConstantCoefficient coef_ls(surface_fit_const);
+   GridFunction surf_fit_mat_gf(&surf_fit_fes);
+   GridFunction surf_fit_gf0(&surf_fit_fes);
+   Array<bool> surf_fit_marker(surf_fit_gf0.Size());
+   ConstantCoefficient surf_fit_coeff(surface_fit_const);
    AdaptivityEvaluator *adapt_surface = NULL;
    if (surface_fit_const > 0.0)
    {
@@ -784,27 +789,27 @@ int main(int argc, char *argv[])
                   "Surface fitting with PA is not implemented yet.");
 
       FunctionCoefficient ls_coeff(surface_level_set);
-      ls_0.ProjectCoefficient(ls_coeff);
+      surf_fit_gf0.ProjectCoefficient(ls_coeff);
 
       for (int i = 0; i < mesh->GetNE(); i++)
       {
-         mat(i) = material_id(i, ls_0);
+         mat(i) = material_id(i, surf_fit_gf0);
          mesh->SetAttribute(i, mat(i) + 1);
       }
 
-      GridFunctionCoefficient coeff_mat(&mat);
-      marker_gf.ProjectDiscCoefficient(coeff_mat, GridFunction::ARITHMETIC);
-      for (int j = 0; j < marker.Size(); j++)
+      GridFunctionCoefficient mat_coeff(&mat);
+      surf_fit_mat_gf.ProjectDiscCoefficient(mat_coeff, GridFunction::ARITHMETIC);
+      for (int j = 0; j < surf_fit_marker.Size(); j++)
       {
-         if (marker_gf(j) > 0.1 && marker_gf(j) < 0.9)
+         if (surf_fit_mat_gf(j) > 0.1 && surf_fit_mat_gf(j) < 0.9)
          {
-            marker[j] = true;
-            marker_gf(j) = 1.0;
+            surf_fit_marker[j] = true;
+            surf_fit_mat_gf(j) = 1.0;
          }
          else
          {
-            marker[j] = false;
-            marker_gf(j) = 0.0;
+            surf_fit_marker[j] = false;
+            surf_fit_mat_gf(j) = 0.0;
          }
       }
 
@@ -819,22 +824,24 @@ int main(int argc, char *argv[])
       }
       else { MFEM_ABORT("Bad interpolation option."); }
 
-      he_nlf_integ->EnableSurfaceFitting(ls_0, marker, coef_ls, *adapt_surface);
+      tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker,
+                                       surf_fit_coeff, *adapt_surface);
       if (visualization)
       {
          socketstream vis1, vis2, vis3;
-         common::VisualizeField(vis1, "localhost", 19916, ls_0, "Level Set 0",
+         common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0, "Level Set 0",
                                 300, 600, 300, 300);
          common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
                                 600, 600, 300, 300);
-         common::VisualizeField(vis3, "localhost", 19916, marker_gf, "Dofs to Move",
+         common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
+                                "Dofs to Move",
                                 900, 600, 300, 300);
       }
    }
 
    // Has to be after the enabling of the limiting / alignment, as it computes
    // normalization factors for these terms as well.
-   if (normalization) { he_nlf_integ->EnableNormalization(x0); }
+   if (normalization) { tmop_integ->EnableNormalization(x0); }
 
    // 12. Setup the final NonlinearForm (which defines the integral of interest,
    //     its first and second derivatives). Here we can use a combination of
@@ -844,39 +851,39 @@ int main(int argc, char *argv[])
    //     metric; one should update those in the code.
    NonlinearForm a(fespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   ConstantCoefficient *coeff1 = NULL;
+   ConstantCoefficient *metric_coeff1 = NULL;
    TMOP_QualityMetric *metric2 = NULL;
    TargetConstructor *target_c2 = NULL;
-   FunctionCoefficient coeff2(weight_fun);
+   FunctionCoefficient metric_coeff2(weight_fun);
 
    // Explicit combination of metrics.
    if (combomet > 0)
    {
       // First metric.
-      coeff1 = new ConstantCoefficient(1.0);
-      he_nlf_integ->SetCoefficient(*coeff1);
+      metric_coeff1 = new ConstantCoefficient(1.0);
+      tmop_integ->SetCoefficient(*metric_coeff1);
 
       // Second metric.
       if (dim == 2) { metric2 = new TMOP_Metric_077; }
       else          { metric2 = new TMOP_Metric_315; }
-      TMOP_Integrator *he_nlf_integ2 = NULL;
+      TMOP_Integrator *tmop_integ2 = NULL;
       if (combomet == 1)
       {
          target_c2 = new TargetConstructor(
             TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE);
          target_c2->SetVolumeScale(0.01);
          target_c2->SetNodes(x0);
-         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
-         he_nlf_integ2->SetCoefficient(coeff2);
+         tmop_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
+         tmop_integ2->SetCoefficient(metric_coeff2);
       }
-      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
-      he_nlf_integ2->SetIntegrationRules(*irules, quad_order);
-      if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
-      he_nlf_integ2->SetExactActionFlag(exactaction);
+      else { tmop_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
+      tmop_integ2->SetIntegrationRules(*irules, quad_order);
+      if (fdscheme) { tmop_integ2->EnableFiniteDifferences(x); }
+      tmop_integ2->SetExactActionFlag(exactaction);
 
       TMOPComboIntegrator *combo = new TMOPComboIntegrator;
-      combo->AddTMOPIntegrator(he_nlf_integ);
-      combo->AddTMOPIntegrator(he_nlf_integ2);
+      combo->AddTMOPIntegrator(tmop_integ);
+      combo->AddTMOPIntegrator(tmop_integ2);
       if (normalization) { combo->EnableNormalization(x0); }
       if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
 
@@ -884,7 +891,7 @@ int main(int argc, char *argv[])
    }
    else
    {
-      a.AddDomainIntegrator(he_nlf_integ);
+      a.AddDomainIntegrator(tmop_integ);
    }
 
    if (pa) { a.Setup(); }
@@ -930,13 +937,13 @@ int main(int argc, char *argv[])
    if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
    {
       lim_coeff.constant = 0.0;
-      coef_zeta.constant = 0.0;
-      coef_ls.constant   = 0.0;
+      adapt_lim_coeff.constant = 0.0;
+      surf_fit_coeff.constant   = 0.0;
       init_metric_energy = a.GetGridFunctionEnergy(x) /
                            (hradaptivity ? mesh->GetNE() : 1);
       lim_coeff.constant = lim_const;
-      coef_zeta.constant = adapt_lim_const;
-      coef_ls.constant   = surface_fit_const;
+      adapt_lim_coeff.constant = adapt_lim_const;
+      surf_fit_coeff.constant   = surface_fit_const;
    }
 
    // Visualize the starting mesh and metric values.
@@ -972,7 +979,7 @@ int main(int argc, char *argv[])
          if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
          if (attr == 4) { n += nd * dim; }
       }
-      Array<int> ess_vdofs(n), vdofs;
+      Array<int> ess_vdofs(n);
       n = 0;
       for (int i = 0; i < mesh->GetNBE(); i++)
       {
@@ -1033,11 +1040,15 @@ int main(int argc, char *argv[])
          if (pa)
          {
             MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
-            S_prec = new OperatorJacobiSmoother;
+            auto js = new OperatorJacobiSmoother;
+            js->SetPositiveDiagonal(true);
+            S_prec = js;
          }
          else
          {
-            S_prec = new DSmoother((lin_solver == 3) ? 0 : 1, 1.0, 1);
+            auto ds = new DSmoother((lin_solver == 3) ? 0 : 1, 1.0, 1);
+            ds->SetPositiveDiagonal(true);
+            S_prec = ds;
          }
          minres->SetPreconditioner(*S_prec);
       }
@@ -1080,7 +1091,7 @@ int main(int argc, char *argv[])
    hr_solver.AddGridFunctionForUpdate(&x0);
    if (adapt_lim_const > 0.)
    {
-      hr_solver.AddGridFunctionForUpdate(&zeta_0);
+      hr_solver.AddGridFunctionForUpdate(&adapt_lim_gf0);
       hr_solver.AddFESpaceForUpdate(&ind_fes);
    }
    hr_solver.Mult();
@@ -1099,13 +1110,13 @@ int main(int argc, char *argv[])
    if (lim_const > 0.0 || adapt_lim_const > 0.0)
    {
       lim_coeff.constant = 0.0;
-      coef_zeta.constant = 0.0;
-      coef_ls.constant   = 0.0;
+      adapt_lim_coeff.constant = 0.0;
+      surf_fit_coeff.constant  = 0.0;
       fin_metric_energy = a.GetGridFunctionEnergy(x) /
                           (hradaptivity ? mesh->GetNE() : 1);
       lim_coeff.constant = lim_const;
-      coef_zeta.constant = adapt_lim_const;
-      coef_ls.constant   = surface_fit_const;
+      adapt_lim_coeff.constant = adapt_lim_const;
+      surf_fit_coeff.constant  = surface_fit_const;
    }
    std::cout << std::scientific << std::setprecision(4);
    cout << "Initial strain energy: " << init_energy
@@ -1127,7 +1138,7 @@ int main(int argc, char *argv[])
    if (adapt_lim_const > 0.0 && visualization)
    {
       socketstream vis0;
-      common::VisualizeField(vis0, "localhost", 19916, zeta_0, "Xi 0",
+      common::VisualizeField(vis0, "localhost", 19916, adapt_lim_gf0, "Xi 0",
                              600, 600, 300, 300);
    }
 
@@ -1138,11 +1149,11 @@ int main(int argc, char *argv[])
          socketstream vis2, vis3;
          common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
                                 600, 900, 300, 300);
-         common::VisualizeField(vis3, "localhost", 19916, marker_gf, "Surface dof",
+         common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf, "Surface dof",
                                 900, 900, 300, 300);
       }
       double err_avg, err_max;
-      he_nlf_integ->GetSurfaceFittingErrors(err_avg, err_max);
+      tmop_integ->GetSurfaceFittingErrors(err_avg, err_max);
       std::cout << "Avg fitting error: " << err_avg << std::endl
                 << "Max fitting error: " << err_max << std::endl;
    }
@@ -1166,8 +1177,8 @@ int main(int argc, char *argv[])
    delete S_prec;
    delete target_c2;
    delete metric2;
-   delete coeff1;
-   delete adapt_evaluator;
+   delete metric_coeff1;
+   delete adapt_lim_eval;
    delete adapt_surface;
    delete target_c;
    delete hr_adapt_coeff;

miniapps/meshing/minimal-surface.cpp

@@ -392,16 +392,17 @@ public:
       void Amr()
       {
          MFEM_VERIFY(opt.amr_threshold >= 0.0 && opt.amr_threshold <= 1.0, "");
-         Mesh *mesh = S.mesh;
+         Mesh *smesh = S.mesh;
          Array<Refinement> amr;
-         const int NE = mesh->GetNE();
+         const int NE = smesh->GetNE();
          DenseMatrix Jadjt, Jadj(DIM, SDIM);
          for (int e = 0; e < NE; e++)
          {
             double minW = +NL_DMAX;
             double maxW = -NL_DMAX;
-            ElementTransformation *eTr = mesh->GetElementTransformation(e);
-            const Geometry::Type &type = mesh->GetElement(e)->GetGeometryType();
+            ElementTransformation *eTr = smesh->GetElementTransformation(e);
+            const Geometry::Type &type =
+               smesh->GetElement(e)->GetGeometryType();
             const IntegrationRule *ir = &IntRules.Get(type, opt.order);
             const int NQ = ir->GetNPoints();
             for (int q = 0; q < NQ; q++)
@@ -424,8 +425,8 @@ public:
          }
          if (amr.Size()>0)
          {
-            mesh->GetNodes()->HostReadWrite();
-            mesh->GeneralRefinement(amr);
+            smesh->GetNodes()->HostReadWrite();
+            smesh->GeneralRefinement(amr);
             S.fes->Update();
             x.HostReadWrite();
             x.Update();

miniapps/meshing/pmesh-optimizer.cpp

@@ -358,9 +358,9 @@ int main (int argc, char *argv[])
       }
       vol_loc += pmesh->GetElementVolume(i);
    }
-   double volume;
-   MPI_Allreduce(&vol_loc, &volume, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-   const double small_phys_size = pow(volume, 1.0 / dim) / 100.0;
+   double vol_glb;
+   MPI_Allreduce(&vol_loc, &vol_glb, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   const double small_phys_size = pow(vol_glb, 1.0 / dim) / 100.0;
 
    // 9. Add a random perturbation to the nodes in the interior of the domain.
    //    We define a random grid function of fespace and make sure that it is
@@ -438,6 +438,10 @@ int main (int argc, char *argv[])
       case 315: metric = new TMOP_Metric_315; break;
       case 316: metric = new TMOP_Metric_316; break;
       case 321: metric = new TMOP_Metric_321; break;
+      case 328: metric = new TMOP_Metric_328(0.5); break;
+      case 332: metric = new TMOP_Metric_332(0.5); break;
+      case 333: metric = new TMOP_Metric_333(0.5); break;
+      case 334: metric = new TMOP_Metric_334(0.5); break;
       // case 352: metric = new TMOP_Metric_352(tauval); break;
       // A-metrics
       case 11: metric = new TMOP_AMetric_011; break;
@@ -486,7 +490,7 @@ int main (int argc, char *argv[])
    H1_FECollection ind_fec(mesh_poly_deg, dim);
    ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
    ParFiniteElementSpace ind_fesv(pmesh, &ind_fec, dim);
-   ParGridFunction size(&ind_fes), aspr(&ind_fes), disc(&ind_fes), ori(&ind_fes);
+   ParGridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
    ParGridFunction aspr3d(&ind_fesv);
 
    const AssemblyLevel al =
@@ -524,13 +528,13 @@ int main (int argc, char *argv[])
          }
          if (dim == 2)
          {
-            FunctionCoefficient ind_coeff(discrete_size_2d);
-            size.ProjectCoefficient(ind_coeff);
+            FunctionCoefficient size_coeff(discrete_size_2d);
+            size.ProjectCoefficient(size_coeff);
          }
          else if (dim == 3)
          {
-            FunctionCoefficient ind_coeff(discrete_size_3d);
-            size.ProjectCoefficient(ind_coeff);
+            FunctionCoefficient size_coeff(discrete_size_3d);
+            size.ProjectCoefficient(size_coeff);
          }
          tc->SetParDiscreteTargetSize(size);
          target_c = tc;
@@ -538,12 +542,12 @@ int main (int argc, char *argv[])
       }
       case 6: // material indicator 2D
       {
-         ParGridFunction d_x(&ind_fes), d_y(&ind_fes);
+         ParGridFunction d_x(&ind_fes), d_y(&ind_fes), disc(&ind_fes);
 
          target_t = TargetConstructor::GIVEN_SHAPE_AND_SIZE;
          DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
-         FunctionCoefficient ind_coeff(material_indicator_2d);
-         disc.ProjectCoefficient(ind_coeff);
+         FunctionCoefficient mat_coeff(material_indicator_2d);
+         disc.ProjectCoefficient(mat_coeff);
          if (adapt_eval == 0)
          {
             tc->SetAdaptivityEvaluator(new AdvectorCG(al));
@@ -678,8 +682,8 @@ int main (int argc, char *argv[])
 
          if (metric_id == 14 || metric_id == 36)
          {
-            ConstantCoefficient ind_coeff(0.1*0.1);
-            size.ProjectCoefficient(ind_coeff);
+            ConstantCoefficient size_coeff(0.1*0.1);
+            size.ProjectCoefficient(size_coeff);
             tc->SetParDiscreteTargetSize(size);
          }
 
@@ -719,16 +723,16 @@ int main (int argc, char *argv[])
       target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
    }
    target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c,
-                                                       h_metric);
+   TMOP_Integrator *tmop_integ = new TMOP_Integrator(metric, target_c,
+                                                     h_metric);
 
    // Finite differences for computations of derivatives.
    if (fdscheme)
    {
       MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
-      he_nlf_integ->EnableFiniteDifferences(x);
+      tmop_integ->EnableFiniteDifferences(x);
    }
-   he_nlf_integ->SetExactActionFlag(exactaction);
+   tmop_integ->SetExactActionFlag(exactaction);
 
    // Setup the quadrature rules for the TMOP integrator.
    IntegrationRules *irules = NULL;
@@ -741,7 +745,7 @@ int main (int argc, char *argv[])
          if (myid == 0) { cout << "Unknown quad_type: " << quad_type << endl; }
          return 3;
    }
-   he_nlf_integ->SetIntegrationRules(*irules, quad_order);
+   tmop_integ->SetIntegrationRules(*irules, quad_order);
    if (myid == 0 && dim == 2)
    {
       cout << "Triangle quadrature points: "
@@ -767,49 +771,50 @@ int main (int argc, char *argv[])
    // The small_phys_size is relevant only with proper normalization.
    if (normalization) { dist = small_phys_size; }
    ConstantCoefficient lim_coeff(lim_const);
-   if (lim_const != 0.0) { he_nlf_integ->EnableLimiting(x0, dist, lim_coeff); }
+   if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
 
    // Adaptive limiting.
-   ParGridFunction zeta_0(&ind_fes);
-   ConstantCoefficient coef_zeta(adapt_lim_const);
-   AdaptivityEvaluator *adapt_evaluator = NULL;
+   ParGridFunction adapt_lim_gf0(&ind_fes);
+   ConstantCoefficient adapt_lim_coeff(adapt_lim_const);
+   AdaptivityEvaluator *adapt_lim_eval = NULL;
    if (adapt_lim_const > 0.0)
    {
       MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
 
-      FunctionCoefficient alim_coeff(adapt_lim_fun);
-      zeta_0.ProjectCoefficient(alim_coeff);
+      FunctionCoefficient adapt_lim_gf0_coeff(adapt_lim_fun);
+      adapt_lim_gf0.ProjectCoefficient(adapt_lim_gf0_coeff);
 
-      if (adapt_eval == 0) { adapt_evaluator = new AdvectorCG(al); }
+      if (adapt_eval == 0) { adapt_lim_eval = new AdvectorCG(al); }
       else if (adapt_eval == 1)
       {
 #ifdef MFEM_USE_GSLIB
-         adapt_evaluator = new InterpolatorFP;
+         adapt_lim_eval = new InterpolatorFP;
 #else
          MFEM_ABORT("MFEM is not built with GSLIB support!");
 #endif
       }
       else { MFEM_ABORT("Bad interpolation option."); }
 
-      he_nlf_integ->EnableAdaptiveLimiting(zeta_0, coef_zeta, *adapt_evaluator);
+      tmop_integ->EnableAdaptiveLimiting(adapt_lim_gf0, adapt_lim_coeff,
+                                         *adapt_lim_eval);
       if (visualization)
       {
          socketstream vis1;
-         common::VisualizeField(vis1, "localhost", 19916, zeta_0, "Zeta 0",
+         common::VisualizeField(vis1, "localhost", 19916, adapt_lim_gf0, "Zeta 0",
                                 300, 600, 300, 300);
       }
    }
 
    // Surface fitting.
    L2_FECollection mat_coll(0, dim);
-   H1_FECollection sigma_fec(mesh_poly_deg, dim);
-   ParFiniteElementSpace sigma_fes(pmesh, &sigma_fec);
+   H1_FECollection surf_fit_fec(mesh_poly_deg, dim);
+   ParFiniteElementSpace surf_fit_fes(pmesh, &surf_fit_fec);
    ParFiniteElementSpace mat_fes(pmesh, &mat_coll);
    ParGridFunction mat(&mat_fes);
-   ParGridFunction marker_gf(&sigma_fes);
-   ParGridFunction ls_0(&sigma_fes);
-   Array<bool> marker(ls_0.Size());
-   ConstantCoefficient coef_ls(surface_fit_const);
+   ParGridFunction surf_fit_mat_gf(&surf_fit_fes);
+   ParGridFunction surf_fit_gf0(&surf_fit_fes);
+   Array<bool> surf_fit_marker(surf_fit_gf0.Size());
+   ConstantCoefficient surf_fit_coeff(surface_fit_const);
    AdaptivityEvaluator *adapt_surface = NULL;
    if (surface_fit_const > 0.0)
    {
@@ -819,27 +824,27 @@ int main (int argc, char *argv[])
                   "Surface fitting with PA is not implemented yet.");
 
       FunctionCoefficient ls_coeff(surface_level_set);
-      ls_0.ProjectCoefficient(ls_coeff);
+      surf_fit_gf0.ProjectCoefficient(ls_coeff);
 
       for (int i = 0; i < pmesh->GetNE(); i++)
       {
-         mat(i) = material_id(i, ls_0);
+         mat(i) = material_id(i, surf_fit_gf0);
          pmesh->SetAttribute(i, mat(i) + 1);
       }
 
       GridFunctionCoefficient coeff_mat(&mat);
-      marker_gf.ProjectDiscCoefficient(coeff_mat, GridFunction::ARITHMETIC);
-      for (int j = 0; j < marker.Size(); j++)
+      surf_fit_mat_gf.ProjectDiscCoefficient(coeff_mat, GridFunction::ARITHMETIC);
+      for (int j = 0; j < surf_fit_marker.Size(); j++)
       {
-         if (marker_gf(j) > 0.1 && marker_gf(j) < 0.9)
+         if (surf_fit_mat_gf(j) > 0.1 && surf_fit_mat_gf(j) < 0.9)
          {
-            marker[j] = true;
-            marker_gf(j) = 1.0;
+            surf_fit_marker[j] = true;
+            surf_fit_mat_gf(j) = 1.0;
          }
          else
          {
-            marker[j] = false;
-            marker_gf(j) = 0.0;
+            surf_fit_marker[j] = false;
+            surf_fit_mat_gf(j) = 0.0;
          }
       }
 
@@ -854,22 +859,24 @@ int main (int argc, char *argv[])
       }
       else { MFEM_ABORT("Bad interpolation option."); }
 
-      he_nlf_integ->EnableSurfaceFitting(ls_0, marker, coef_ls, *adapt_surface);
+      tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker, surf_fit_coeff,
+                                       *adapt_surface);
       if (visualization)
       {
          socketstream vis1, vis2, vis3;
-         common::VisualizeField(vis1, "localhost", 19916, ls_0, "Level Set 0",
+         common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0, "Level Set 0",
                                 300, 600, 300, 300);
          common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
                                 600, 600, 300, 300);
-         common::VisualizeField(vis3, "localhost", 19916, marker_gf, "Dofs to Move",
+         common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
+                                "Dofs to Move",
                                 900, 600, 300, 300);
       }
    }
 
    // Has to be after the enabling of the limiting / alignment, as it computes
    // normalization factors for these terms as well.
-   if (normalization) { he_nlf_integ->ParEnableNormalization(x0); }
+   if (normalization) { tmop_integ->ParEnableNormalization(x0); }
 
    // 13. Setup the final NonlinearForm (which defines the integral of interest,
    //     its first and second derivatives). Here we can use a combination of
@@ -879,39 +886,39 @@ int main (int argc, char *argv[])
    //     metric; one should update those in the code.
    ParNonlinearForm a(pfespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   ConstantCoefficient *coeff1 = NULL;
+   ConstantCoefficient *metric_coeff1 = NULL;
    TMOP_QualityMetric *metric2 = NULL;
    TargetConstructor *target_c2 = NULL;
-   FunctionCoefficient coeff2(weight_fun);
+   FunctionCoefficient metric_coeff2(weight_fun);
 
    // Explicit combination of metrics.
    if (combomet > 0)
    {
       // First metric.
-      coeff1 = new ConstantCoefficient(1.0);
-      he_nlf_integ->SetCoefficient(*coeff1);
+      metric_coeff1 = new ConstantCoefficient(1.0);
+      tmop_integ->SetCoefficient(*metric_coeff1);
 
       // Second metric.
       if (dim == 2) { metric2 = new TMOP_Metric_077; }
       else          { metric2 = new TMOP_Metric_315; }
-      TMOP_Integrator *he_nlf_integ2 = NULL;
+      TMOP_Integrator *tmop_integ2 = NULL;
       if (combomet == 1)
       {
          target_c2 = new TargetConstructor(
             TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE, MPI_COMM_WORLD);
          target_c2->SetVolumeScale(0.01);
          target_c2->SetNodes(x0);
-         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
-         he_nlf_integ2->SetCoefficient(coeff2);
+         tmop_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
+         tmop_integ2->SetCoefficient(metric_coeff2);
       }
-      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
-      he_nlf_integ2->SetIntegrationRules(*irules, quad_order);
-      if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
-      he_nlf_integ2->SetExactActionFlag(exactaction);
+      else { tmop_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
+      tmop_integ2->SetIntegrationRules(*irules, quad_order);
+      if (fdscheme) { tmop_integ2->EnableFiniteDifferences(x); }
+      tmop_integ2->SetExactActionFlag(exactaction);
 
       TMOPComboIntegrator *combo = new TMOPComboIntegrator;
-      combo->AddTMOPIntegrator(he_nlf_integ);
-      combo->AddTMOPIntegrator(he_nlf_integ2);
+      combo->AddTMOPIntegrator(tmop_integ);
+      combo->AddTMOPIntegrator(tmop_integ2);
       if (normalization) { combo->ParEnableNormalization(x0); }
       if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
 
@@ -919,7 +926,7 @@ int main (int argc, char *argv[])
    }
    else
    {
-      a.AddDomainIntegrator(he_nlf_integ);
+      a.AddDomainIntegrator(tmop_integ);
    }
 
    if (pa) { a.Setup(); }
@@ -971,13 +978,13 @@ int main (int argc, char *argv[])
    if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
    {
       lim_coeff.constant = 0.0;
-      coef_zeta.constant = 0.0;
-      coef_ls.constant   = 0.0;
+      adapt_lim_coeff.constant = 0.0;
+      surf_fit_coeff.constant   = 0.0;
       init_metric_energy = a.GetParGridFunctionEnergy(x) /
                            (hradaptivity ? pmesh->GetGlobalNE() : 1);
       lim_coeff.constant = lim_const;
-      coef_zeta.constant = adapt_lim_const;
-      coef_ls.constant   = surface_fit_const;
+      adapt_lim_coeff.constant = adapt_lim_const;
+      surf_fit_coeff.constant  = surface_fit_const;
    }
 
    // Visualize the starting mesh and metric values.
@@ -990,9 +997,8 @@ int main (int argc, char *argv[])
 
    // 14. Fix all boundary nodes, or fix only a given component depending on the
    //     boundary attributes of the given mesh.  Attributes 1/2/3 correspond to
-   //     fixed x/y/z components of the node.  Attribute 4 corresponds to an
-   //     entirely fixed node.  Other boundary attributes do not affect the node
-   //     movement boundary conditions.
+   //     fixed x/y/z components of the node.  Attribute dim+1 corresponds to
+   //     an entirely fixed node.
    if (move_bnd == false)
    {
       Array<int> ess_bdr(pmesh->bdr_attributes.Max());
@@ -1013,7 +1019,7 @@ int main (int argc, char *argv[])
          if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
          if (attr == 4) { n += nd * dim; }
       }
-      Array<int> ess_vdofs(n), vdofs;
+      Array<int> ess_vdofs(n);
       n = 0;
       for (int i = 0; i < pmesh->GetNBE(); i++)
       {
@@ -1074,13 +1080,16 @@ int main (int argc, char *argv[])
          if (pa)
          {
             MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
-            S_prec = new OperatorJacobiSmoother;
+            auto js = new OperatorJacobiSmoother;
+            js->SetPositiveDiagonal(true);
+            S_prec = js;
          }
          else
          {
-            HypreSmoother *hs = new HypreSmoother;
+            auto hs = new HypreSmoother;
             hs->SetType((lin_solver == 3) ? HypreSmoother::Jacobi
-                        : HypreSmoother::l1Jacobi, 1);
+                        /* */             : HypreSmoother::l1Jacobi, 1);
+            hs->SetPositiveDiagonal(true);
             S_prec = hs;
          }
          minres->SetPreconditioner(*S_prec);
@@ -1124,7 +1133,7 @@ int main (int argc, char *argv[])
    hr_solver.AddGridFunctionForUpdate(&x0);
    if (adapt_lim_const > 0.)
    {
-      hr_solver.AddGridFunctionForUpdate(&zeta_0);
+      hr_solver.AddGridFunctionForUpdate(&adapt_lim_gf0);
       hr_solver.AddFESpaceForUpdate(&ind_fes);
    }
    hr_solver.Mult();
@@ -1146,13 +1155,13 @@ int main (int argc, char *argv[])
    if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
    {
       lim_coeff.constant = 0.0;
-      coef_zeta.constant = 0.0;
-      coef_ls.constant   = 0.0;
+      adapt_lim_coeff.constant = 0.0;
+      surf_fit_coeff.constant  = 0.0;
       fin_metric_energy  = a.GetParGridFunctionEnergy(x) /
                            (hradaptivity ? pmesh->GetGlobalNE() : 1);
       lim_coeff.constant = lim_const;
-      coef_zeta.constant = adapt_lim_const;
-      coef_ls.constant   = surface_fit_const;
+      adapt_lim_coeff.constant = adapt_lim_const;
+      surf_fit_coeff.constant  = surface_fit_const;
    }
    if (myid == 0)
    {
@@ -1177,7 +1186,7 @@ int main (int argc, char *argv[])
    if (adapt_lim_const > 0.0 && visualization)
    {
       socketstream vis0;
-      common::VisualizeField(vis0, "localhost", 19916, zeta_0, "Xi 0",
+      common::VisualizeField(vis0, "localhost", 19916, adapt_lim_gf0, "Xi 0",
                              600, 600, 300, 300);
    }
 
@@ -1188,11 +1197,11 @@ int main (int argc, char *argv[])
          socketstream vis2, vis3;
          common::VisualizeField(vis2, "localhost", 19916, mat,
                                 "Materials", 600, 900, 300, 300);
-         common::VisualizeField(vis3, "localhost", 19916, marker_gf,
+         common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
                                 "Surface dof", 900, 900, 300, 300);
       }
       double err_avg, err_max;
-      he_nlf_integ->GetSurfaceFittingErrors(err_avg, err_max);
+      tmop_integ->GetSurfaceFittingErrors(err_avg, err_max);
       if (myid == 0)
       {
          std::cout << "Avg fitting error: " << err_avg << std::endl
@@ -1226,8 +1235,8 @@ int main (int argc, char *argv[])
    delete S_prec;
    delete target_c2;
    delete metric2;
-   delete coeff1;
-   delete adapt_evaluator;
+   delete metric_coeff1;
+   delete adapt_lim_eval;
    delete adapt_surface;
    delete target_c;
    delete hr_adapt_coeff;

miniapps/meshing/pminimal-surface.cpp

@@ -399,16 +399,18 @@ public:
       void Amr()
       {
          MFEM_VERIFY(opt.amr_threshold >= 0.0 && opt.amr_threshold <= 1.0, "");
-         Mesh *mesh = S.mesh;
+         Mesh *smesh = S.mesh;
          Array<Refinement> amr;
-         const int NE = mesh->GetNE();
+         const int NE = smesh->GetNE();
          DenseMatrix Jadjt, Jadj(DIM, SDIM);
          for (int e = 0; e < NE; e++)
          {
             double minW = +NL_DMAX;
             double maxW = -NL_DMAX;
-            ElementTransformation *eTr = mesh->GetElementTransformation(e);
-            const Geometry::Type &type = mesh->GetElement(e)->GetGeometryType();
+            ElementTransformation *eTr = smesh->GetElementTransformation(e);
+            const Geometry::Type &type =
+               smesh->GetElement(e)->GetGeometryType();
+
             const IntegrationRule *ir = &IntRules.Get(type, opt.order);
             const int NQ = ir->GetNPoints();
             for (int q = 0; q < NQ; q++)
@@ -431,8 +433,8 @@ public:
          }
          if (amr.Size()>0)
          {
-            mesh->GetNodes()->HostReadWrite();
-            mesh->GeneralRefinement(amr);
+            smesh->GetNodes()->HostReadWrite();
+            smesh->GeneralRefinement(amr);
             S.fes->Update();
             x.HostReadWrite();
             x.Update();

miniapps/meshing/polar-nc.cpp

@@ -185,9 +185,9 @@ Mesh* Make2D(int nsteps, double rstep, double phi, double aspect, int order,
 
          for (int j = 0, k = 0; j < coords.Size(); k++, j += 2)
          {
-            int sfc = ((i & 1) ? coords[j] : (width-1 - coords[j]))
-                      + coords[j+1]*width;
-            int old_index = beg + sfc;
+            int sfc_index = ((i & 1) ? coords[j] : (width-1 - coords[j]))
+                            + coords[j+1]*width;
+            int old_index = beg + sfc_index;
 
             ordering[old_index] = beg + k;
             new_params[beg + k] = params[old_index];
@@ -221,7 +221,7 @@ Mesh* Make2D(int nsteps, double rstep, double phi, double aspect, int order,
 
          for (int j = 0; j < dofs.Size(); j++)
          {
-            double r, a;
+            double a;
             if (geom == Geometry::SQUARE)
             {
                r = par.r + ir[j].x * par.dr;
@@ -449,7 +449,7 @@ Mesh* Make3D(int nsteps, double rstep, double aspect, int order, bool sfc)
          {
             const IntegrationPoint &ip = ir[j];
 
-            double u, v, w, r;
+            double u, v, w;
             if (geom == Geometry::PRISM)
             {
                double l1 = 1.0 - ip.x - ip.y;

miniapps/meshing/shaper.cpp

@@ -158,22 +158,22 @@ int main(int argc, char *argv[])
                const int s = sd+1;
                if (dim == 2)
                {
-                  for (int j = 0; j <= sd; j++)
-                     for (int i = 0; i < sd; i++)
+                  for (int jj = 0; jj <= sd; jj++)
+                     for (int ii = 0; ii < sd; ii++)
                      {
-                        dx += abs(mat[j*s + i+1] - mat[j*s + i]);
-                        dy += abs(mat[(i+1)*s + j] - mat[i*s + j]);
+                        dx += abs(mat[jj*s + ii+1] - mat[jj*s + ii]);
+                        dy += abs(mat[(ii+1)*s + jj] - mat[ii*s + jj]);
                      }
                }
                else if (dim == 3)
                {
-                  for (int k = 0; k <= sd; k++)
-                     for (int j = 0; j <= sd; j++)
-                        for (int i = 0; i < sd; i++)
+                  for (int kk = 0; kk <= sd; kk++)
+                     for (int jj = 0; jj <= sd; jj++)
+                        for (int ii = 0; ii < sd; ii++)
                         {
-                           dx += abs(mat[(k*s + j)*s + i+1] - mat[(k*s + j)*s + i]);
-                           dy += abs(mat[(k*s + i+1)*s + j] - mat[(k*s + i)*s + j]);
-                           dz += abs(mat[((i+1)*s + j)*s + k] - mat[(i*s + j)*s + k]);
+                           dx += abs(mat[(kk*s + jj)*s + ii+1] - mat[(kk*s + jj)*s + ii]);
+                           dy += abs(mat[(kk*s + ii+1)*s + jj] - mat[(kk*s + ii)*s + jj]);
+                           dz += abs(mat[((ii+1)*s + jj)*s + kk] - mat[(ii*s + jj)*s + kk]);
                         }
                }
                type = 0;

miniapps/meshing/toroid.cpp

@@ -127,18 +127,20 @@ int main(int argc, char *argv[])
    }
 
    // Add Elements of the desired type
-   int v[8];
-   for (int i=0; i < nphi_; i++)
    {
-      if (el_type_ == Element::WEDGE)
+      int v[8];
+      for (int i=0; i < nphi_; i++)
       {
-         for (int j = 0; j < 6; j++) { v[j] = 3*i+j; }
-         mesh->AddWedge(v);
-      }
-      else
-      {
-         for (int j = 0; j < 8; j++) { v[j] = 4*i+j; }
-         mesh->AddHex(v);
+         if (el_type_ == Element::WEDGE)
+         {
+            for (int j = 0; j < 6; j++) { v[j] = 3*i+j; }
+            mesh->AddWedge(v);
+         }
+         else
+         {
+            for (int j = 0; j < 8; j++) { v[j] = 4*i+j; }
+            mesh->AddHex(v);
+         }
       }
    }
    mesh->FinalizeTopology();

miniapps/navier/navier_solver.cpp

@@ -341,11 +341,12 @@ void NavierSolver::UpdateTimestepHistory(double dt)
    un_gf.SetFromTrueDofs(un);
 }
 
-void NavierSolver::Step(double &time, double dt, int cur_step, bool provisional)
+void NavierSolver::Step(double &time, double dt, int current_step,
+                        bool provisional)
 {
    sw_step.Start();
 
-   SetTimeIntegrationCoefficients(cur_step);
+   SetTimeIntegrationCoefficients(current_step);
 
    // Set current time for velocity Dirichlet boundary conditions.
    for (auto &vel_dbc : vel_dbcs)
@@ -850,7 +851,7 @@ void NavierSolver::ComputeCurl2D(ParGridFunction &u,
 
 double NavierSolver::ComputeCFL(ParGridFunction &u, double dt)
 {
-   ParMesh *pmesh = u.ParFESpace()->GetParMesh();
+   ParMesh *pmesh_u = u.ParFESpace()->GetParMesh();
    FiniteElementSpace *fes = u.FESpace();
    int vdim = fes->GetVDim();
 
@@ -879,7 +880,7 @@ double NavierSolver::ComputeCFL(ParGridFunction &u, double dt)
          ut.SetSize(uz.Size());
       }
 
-      double hmin = pmesh->GetElementSize(e, 1) /
+      double hmin = pmesh_u->GetElementSize(e, 1) /
                     (double) fes->GetElementOrder(0);
 
       for (int i = 0; i < ir.GetNPoints(); ++i)
@@ -924,7 +925,7 @@ double NavierSolver::ComputeCFL(ParGridFunction &u, double dt)
                  1,
                  MPI_DOUBLE,
                  MPI_MAX,
-                 pmesh->GetComm());
+                 pmesh_u->GetComm());
 
    return cflmax_global;
 }

miniapps/navier/navier_tgv.cpp

@@ -366,11 +366,11 @@ int main(int argc, char *argv[])
          pvdc.Save();
       }
 
-      double u_inf_loc = u_gf->Normlinf();
-      double p_inf_loc = p_gf->Normlinf();
-      double u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
-      double p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
-      double ke = kin_energy.ComputeKineticEnergy(*u_gf);
+      u_inf_loc = u_gf->Normlinf();
+      p_inf_loc = p_gf->Normlinf();
+      u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
+      p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
+      ke = kin_energy.ComputeKineticEnergy(*u_gf);
       if (mpi.Root())
       {
          printf("%.5E %.5E %.5E %.5E %.5E\n", t, dt, u_inf, p_inf, ke);

miniapps/nurbs/nurbs_ex1.cpp

@@ -109,11 +109,11 @@ public:
             w *= Q->Eval(Trans, ip);
          }
 
-         for (int j = 0; j < nd; j++)
+         for (int jj = 0; jj < nd; jj++)
          {
-            for (int i = 0; i < nd; i++)
+            for (int ii = 0; ii < nd; ii++)
             {
-               elmat(i, j) += w*shape(i)*laplace(j);
+               elmat(ii, jj) += w*shape(ii)*laplace(jj);
             }
          }
       }

miniapps/nurbs/nurbs_ex1p.cpp

@@ -121,11 +121,11 @@ public:
             w *= Q->Eval(Trans, ip);
          }
 
-         for (int j = 0; j < nd; j++)
+         for (int jj = 0; jj < nd; jj++)
          {
-            for (int i = 0; i < nd; i++)
+            for (int ii = 0; ii < nd; ii++)
             {
-               elmat(i, j) += w*shape(i)*laplace(j);
+               elmat(ii, jj) += w*shape(ii)*laplace(jj);
             }
          }
       }

miniapps/performance/ex1.cpp

@@ -13,6 +13,15 @@
 //               ex1 -m ../../data/ball-nurbs.mesh -std  -asm -pc ho  -sc
 //               ex1 -m ../../data/pipe-nurbs.mesh -perf -mf  -pc lor
 //               ex1 -m ../../data/pipe-nurbs.mesh -std  -asm -pc ho  -sc
+//               ex1 -m ../../data/star.mesh -perf -mf  -pc lor
+//               ex1 -m ../../data/star.mesh -perf -asm -pc ho
+//               ex1 -m ../../data/star.mesh -perf -asm -pc ho -sc
+//               ex1 -m ../../data/star.mesh -std  -asm -pc ho
+//               ex1 -m ../../data/star.mesh -std  -asm -pc ho -sc
+//               ex1 -m ../../data/amr-quad.mesh -perf -asm -pc ho -sc
+//               ex1 -m ../../data/amr-quad.mesh -std  -asm -pc ho -sc
+//               ex1 -m ../../data/disc-nurbs.mesh -perf -asm -pc ho  -sc
+//               ex1 -m ../../data/disc-nurbs.mesh -std  -asm -pc ho  -sc
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Laplace problem
@@ -36,29 +45,45 @@
 using namespace std;
 using namespace mfem;
 
+enum class PCType { NONE, LOR, HO };
+
 // Define template parameters for optimized build.
-const Geometry::Type geom     = Geometry::CUBE; // mesh elements  (default: hex)
-const int            mesh_p   = 3;              // mesh curvature (default: 3)
-const int            sol_p    = 3;              // solution order (default: 3)
-const int            rdim     = Geometry::Constants<geom>::Dimension;
-const int            ir_order = 2*sol_p+rdim-1;
+template <int dim> struct geom_t { };
+template <>
+struct geom_t<2> { static const Geometry::Type value = Geometry::SQUARE; };
+template <>
+struct geom_t<3> { static const Geometry::Type value = Geometry::CUBE; };
 
-// Static mesh type
-typedef H1_FiniteElement<geom,mesh_p>         mesh_fe_t;
-typedef H1_FiniteElementSpace<mesh_fe_t>      mesh_fes_t;
-typedef TMesh<mesh_fes_t>                     mesh_t;
+const int mesh_p   = 3; // mesh curvature (default: 3)
+const int sol_p    = 3; // solution order (default: 3)
 
-// Static solution finite element space type
-typedef H1_FiniteElement<geom,sol_p>          sol_fe_t;
-typedef H1_FiniteElementSpace<sol_fe_t>       sol_fes_t;
+template <int dim>
+struct ex1_t
+{
+   static const Geometry::Type geom = geom_t<dim>::value;
+   static const int            rdim     = Geometry::Constants<geom>::Dimension;
+   static const int            ir_order = 2*sol_p+rdim-1;
 
-// Static quadrature, coefficient and integrator types
-typedef TIntegrationRule<geom,ir_order>       int_rule_t;
-typedef TConstantCoefficient<>                coeff_t;
-typedef TIntegrator<coeff_t,TDiffusionKernel> integ_t;
+   // Static mesh type
+   using mesh_fe_t = H1_FiniteElement<geom,mesh_p>;
+   using mesh_fes_t = H1_FiniteElementSpace<mesh_fe_t>;
+   using mesh_t = TMesh<mesh_fes_t>;
 
-// Static bilinear form type, combining the above types
-typedef TBilinearForm<mesh_t,sol_fes_t,int_rule_t,integ_t> HPCBilinearForm;
+   // Static solution finite element space type
+   using sol_fe_t = H1_FiniteElement<geom,sol_p>;
+   using sol_fes_t = H1_FiniteElementSpace<sol_fe_t>;
+
+   // Static quadrature, coefficient and integrator types
+   using int_rule_t = TIntegrationRule<geom,ir_order>;
+   using coeff_t = TConstantCoefficient<>;
+   using integ_t = TIntegrator<coeff_t,TDiffusionKernel>;
+
+   using HPCBilinearForm = TBilinearForm<mesh_t,sol_fes_t,int_rule_t,integ_t>;
+
+   static int run(Mesh *mesh, int ref_levels, int order, int basis,
+                  bool static_cond, PCType pc_choice, bool perf,
+                  bool matrix_free, bool visualization);
+};
 
 int main(int argc, char *argv[])
 {
@@ -113,11 +138,10 @@ int main(int argc, char *argv[])
                "--standard-version is not compatible with --matrix-free");
    args.PrintOptions(cout);
 
-   enum PCType { NONE, LOR, HO };
    PCType pc_choice;
-   if (!strcmp(pc, "ho")) { pc_choice = HO; }
-   else if (!strcmp(pc, "lor")) { pc_choice = LOR; }
-   else if (!strcmp(pc, "none")) { pc_choice = NONE; }
+   if (!strcmp(pc, "ho")) { pc_choice = PCType::HO; }
+   else if (!strcmp(pc, "lor")) { pc_choice = PCType::LOR; }
+   else if (!strcmp(pc, "none")) { pc_choice = PCType::NONE; }
    else
    {
       mfem_error("Invalid Preconditioner specified");
@@ -137,6 +161,29 @@ int main(int argc, char *argv[])
    Mesh *mesh = new Mesh(mesh_file, 1, 1);
    int dim = mesh->Dimension();
 
+   if (dim == 2)
+   {
+      return ex1_t<2>::run(mesh, ref_levels, order, basis, static_cond,
+                           pc_choice, perf, matrix_free, visualization);
+   }
+   else if (dim == 3)
+   {
+      return ex1_t<3>::run(mesh, ref_levels, order, basis, static_cond,
+                           pc_choice, perf, matrix_free, visualization);
+   }
+   else
+   {
+      MFEM_ABORT("Dimension must be 2 or 3.")
+   }
+
+   return 0;
+}
+
+template <int dim>
+int ex1_t<dim>::run(Mesh *mesh, int ref_levels, int order, int basis,
+                    bool static_cond, PCType pc_choice, bool perf,
+                    bool matrix_free, bool visualization)
+{
    // 3. Check if the optimized version matches the given mesh
    if (perf)
    {
@@ -172,8 +219,8 @@ int main(int argc, char *argv[])
    }
    if (mesh->MeshGenerator() & 1) // simplex mesh
    {
-      MFEM_VERIFY(pc_choice != LOR, "triangle and tet meshes do not support"
-                  " the LOR preconditioner yet");
+      MFEM_VERIFY(pc_choice != PCType::LOR, "triangle and tet meshes do not "
+                  " support the LOR preconditioner yet");
    }
 
    // 5. Define a finite element space on the mesh. Here we use continuous
@@ -202,7 +249,7 @@ int main(int argc, char *argv[])
    Mesh mesh_lor;
    FiniteElementCollection *fec_lor = NULL;
    FiniteElementSpace *fespace_lor = NULL;
-   if (pc_choice == LOR)
+   if (pc_choice == PCType::LOR)
    {
       int basis_lor = basis;
       if (basis == BasisType::Positive) { basis_lor=BasisType::ClosedUniform; }
@@ -253,8 +300,8 @@ int main(int argc, char *argv[])
    //     Optionally setup a form to be assembled for preconditioning (a_pc).
    BilinearForm *a = new BilinearForm(fespace);
    BilinearForm *a_pc = NULL;
-   if (pc_choice == LOR) { a_pc = new BilinearForm(fespace_lor); }
-   if (pc_choice == HO)  { a_pc = new BilinearForm(fespace); }
+   if (pc_choice == PCType::LOR) { a_pc = new BilinearForm(fespace_lor); }
+   if (pc_choice == PCType::HO)  { a_pc = new BilinearForm(fespace); }
 
    // 11. Assemble the bilinear form and the corresponding linear system,
    //     applying any necessary transformations such as: eliminating boundary
@@ -263,7 +310,7 @@ int main(int argc, char *argv[])
    if (static_cond)
    {
       a->EnableStaticCondensation();
-      MFEM_VERIFY(pc_choice != LOR,
+      MFEM_VERIFY(pc_choice != PCType::LOR,
                   "cannot use LOR preconditioner with static condensation");
    }
 
@@ -322,7 +369,7 @@ int main(int argc, char *argv[])
    tic_toc.Start();
 
    SparseMatrix A_pc;
-   if (pc_choice == LOR)
+   if (pc_choice == PCType::LOR)
    {
       // TODO: assemble the LOR matrix using the performance code
       a_pc->AddDomainIntegrator(new DiffusionIntegrator(one));
@@ -330,7 +377,7 @@ int main(int argc, char *argv[])
       a_pc->Assemble();
       a_pc->FormSystemMatrix(ess_tdof_list, A_pc);
    }
-   else if (pc_choice == HO)
+   else if (pc_choice == PCType::HO)
    {
       if (!matrix_free)
       {
@@ -348,7 +395,7 @@ int main(int argc, char *argv[])
    cout << " done, " << tic_toc.RealTime() << "s." << endl;
 
    // Solve with CG or PCG, depending if the matrix A_pc is available
-   if (pc_choice != NONE)
+   if (pc_choice != PCType::NONE)
    {
       GSSmoother M(A_pc);
       PCG(*a_oper, M, B, X, 1, 500, 1e-12, 0.0);

miniapps/performance/ex1p.cpp

@@ -13,6 +13,15 @@
 //               mpirun -np 4 ex1p -m ../../data/ball-nurbs.mesh -std  -asm -pc ho  -sc
 //               mpirun -np 4 ex1p -m ../../data/pipe-nurbs.mesh -perf -mf  -pc lor
 //               mpirun -np 4 ex1p -m ../../data/pipe-nurbs.mesh -std  -asm -pc ho  -sc
+//               mpirun -np 4 ex1p -m ../../data/star.mesh -perf -mf  -pc lor
+//               mpirun -np 4 ex1p -m ../../data/star.mesh -perf -asm -pc ho
+//               mpirun -np 4 ex1p -m ../../data/star.mesh -perf -asm -pc ho -sc
+//               mpirun -np 4 ex1p -m ../../data/star.mesh -std  -asm -pc ho
+//               mpirun -np 4 ex1p -m ../../data/star.mesh -std  -asm -pc ho -sc
+//               mpirun -np 4 ex1p -m ../../data/amr-quad.mesh -perf -asm -pc ho -sc
+//               mpirun -np 4 ex1p -m ../../data/amr-quad.mesh -std  -asm -pc ho -sc
+//               mpirun -np 4 ex1p -m ../../data/disc-nurbs.mesh -perf -asm -pc ho  -sc
+//               mpirun -np 4 ex1p -m ../../data/disc-nurbs.mesh -std  -asm -pc ho  -sc
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Laplace problem
@@ -36,29 +45,45 @@
 using namespace std;
 using namespace mfem;
 
+enum class PCType { NONE, LOR, HO };
+
 // Define template parameters for optimized build.
-const Geometry::Type geom     = Geometry::CUBE; // mesh elements  (default: hex)
-const int            mesh_p   = 3;              // mesh curvature (default: 3)
-const int            sol_p    = 3;              // solution order (default: 3)
-const int            rdim     = Geometry::Constants<geom>::Dimension;
-const int            ir_order = 2*sol_p+rdim-1;
+template <int dim> struct geom_t { };
+template <>
+struct geom_t<2> { static const Geometry::Type value = Geometry::SQUARE; };
+template <>
+struct geom_t<3> { static const Geometry::Type value = Geometry::CUBE; };
 
-// Static mesh type
-typedef H1_FiniteElement<geom,mesh_p>         mesh_fe_t;
-typedef H1_FiniteElementSpace<mesh_fe_t>      mesh_fes_t;
-typedef TMesh<mesh_fes_t>                     mesh_t;
+const int mesh_p   = 3; // mesh curvature (default: 3)
+const int sol_p    = 3; // solution order (default: 3)
 
-// Static solution finite element space type
-typedef H1_FiniteElement<geom,sol_p>          sol_fe_t;
-typedef H1_FiniteElementSpace<sol_fe_t>       sol_fes_t;
+template <int dim>
+struct ex1_t
+{
+   static const Geometry::Type geom = geom_t<dim>::value;
+   static const int            rdim     = Geometry::Constants<geom>::Dimension;
+   static const int            ir_order = 2*sol_p+rdim-1;
 
-// Static quadrature, coefficient and integrator types
-typedef TIntegrationRule<geom,ir_order>       int_rule_t;
-typedef TConstantCoefficient<>                coeff_t;
-typedef TIntegrator<coeff_t,TDiffusionKernel> integ_t;
+   // Static mesh type
+   using mesh_fe_t = H1_FiniteElement<geom,mesh_p>;
+   using mesh_fes_t = H1_FiniteElementSpace<mesh_fe_t>;
+   using mesh_t = TMesh<mesh_fes_t>;
 
-// Static bilinear form type, combining the above types
-typedef TBilinearForm<mesh_t,sol_fes_t,int_rule_t,integ_t> HPCBilinearForm;
+   // Static solution finite element space type
+   using sol_fe_t = H1_FiniteElement<geom,sol_p>;
+   using sol_fes_t = H1_FiniteElementSpace<sol_fe_t>;
+
+   // Static quadrature, coefficient and integrator types
+   using int_rule_t = TIntegrationRule<geom,ir_order>;
+   using coeff_t = TConstantCoefficient<>;
+   using integ_t = TIntegrator<coeff_t,TDiffusionKernel>;
+
+   using HPCBilinearForm = TBilinearForm<mesh_t,sol_fes_t,int_rule_t,integ_t>;
+
+   static int run(Mesh *mesh, int ser_ref_levels, int par_ref_levels, int order,
+                  int basis, bool static_cond, PCType pc_choice, bool perf,
+                  bool matrix_free, bool visualization);
+};
 
 int main(int argc, char *argv[])
 {
@@ -69,7 +94,11 @@ int main(int argc, char *argv[])
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
 
    // 2. Parse command-line options.
+#ifdef MFEM_HPC_EX1_2D
+   const char *mesh_file = "../../data/star.mesh";
+#else
    const char *mesh_file = "../../data/fichera.mesh";
+#endif
    int ser_ref_levels = -1;
    int par_ref_levels = 1;
    int order = sol_p;
@@ -133,11 +162,10 @@ int main(int argc, char *argv[])
       args.PrintOptions(cout);
    }
 
-   enum PCType { NONE, LOR, HO };
    PCType pc_choice;
-   if (!strcmp(pc, "ho")) { pc_choice = HO; }
-   else if (!strcmp(pc, "lor")) { pc_choice = LOR; }
-   else if (!strcmp(pc, "none")) { pc_choice = NONE; }
+   if (!strcmp(pc, "ho")) { pc_choice = PCType::HO; }
+   else if (!strcmp(pc, "lor")) { pc_choice = PCType::LOR; }
+   else if (!strcmp(pc, "none")) { pc_choice = PCType::NONE; }
    else
    {
       mfem_error("Invalid Preconditioner specified");
@@ -163,6 +191,35 @@ int main(int argc, char *argv[])
    Mesh *mesh = new Mesh(mesh_file, 1, 1);
    int dim = mesh->Dimension();
 
+   if (dim == 2)
+   {
+      return ex1_t<2>::run(mesh, ser_ref_levels, par_ref_levels, order, basis,
+                           static_cond, pc_choice, perf, matrix_free,
+                           visualization);
+   }
+   else if (dim == 3)
+   {
+      return ex1_t<3>::run(mesh, ser_ref_levels, par_ref_levels, order,
+                           basis, static_cond, pc_choice, perf, matrix_free,
+                           visualization);
+   }
+   else
+   {
+      MFEM_ABORT("Dimension must be 2 or 3.")
+   }
+
+   return 0;
+}
+
+template <int dim>
+int ex1_t<dim>::run(Mesh *mesh, int ser_ref_levels, int par_ref_levels,
+                    int order, int basis, bool static_cond, PCType pc_choice,
+                    bool perf, bool matrix_free, bool visualization)
+{
+   int num_procs, myid;
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
    // 4. Check if the optimized version matches the given mesh
    if (perf)
    {
@@ -230,8 +287,8 @@ int main(int argc, char *argv[])
    }
    if (pmesh->MeshGenerator() & 1) // simplex mesh
    {
-      MFEM_VERIFY(pc_choice != LOR, "triangle and tet meshes do not support"
-                  " the LOR preconditioner yet");
+      MFEM_VERIFY(pc_choice != PCType::LOR, "triangle and tet meshes do not "
+                  "support the LOR preconditioner yet");
    }
 
    // 7. Define a parallel finite element space on the parallel mesh. Here we
@@ -264,7 +321,7 @@ int main(int argc, char *argv[])
    ParMesh pmesh_lor;
    FiniteElementCollection *fec_lor = NULL;
    ParFiniteElementSpace *fespace_lor = NULL;
-   if (pc_choice == LOR)
+   if (pc_choice == PCType::LOR)
    {
       int basis_lor = basis;
       if (basis == BasisType::Positive) { basis_lor=BasisType::ClosedUniform; }
@@ -318,8 +375,8 @@ int main(int argc, char *argv[])
    //     that will hold the matrix corresponding to the Laplacian operator.
    ParBilinearForm *a = new ParBilinearForm(fespace);
    ParBilinearForm *a_pc = NULL;
-   if (pc_choice == LOR) { a_pc = new ParBilinearForm(fespace_lor); }
-   if (pc_choice == HO)  { a_pc = new ParBilinearForm(fespace); }
+   if (pc_choice == PCType::LOR) { a_pc = new ParBilinearForm(fespace_lor); }
+   if (pc_choice == PCType::HO)  { a_pc = new ParBilinearForm(fespace); }
 
    // 13. Assemble the parallel bilinear form and the corresponding linear
    //     system, applying any necessary transformations such as: parallel
@@ -328,7 +385,7 @@ int main(int argc, char *argv[])
    if (static_cond)
    {
       a->EnableStaticCondensation();
-      MFEM_VERIFY(pc_choice != LOR,
+      MFEM_VERIFY(pc_choice != PCType::LOR,
                   "cannot use LOR preconditioner with static condensation");
    }
 
@@ -404,7 +461,7 @@ int main(int argc, char *argv[])
    tic_toc.Start();
 
    HypreParMatrix A_pc;
-   if (pc_choice == LOR)
+   if (pc_choice == PCType::LOR)
    {
       // TODO: assemble the LOR matrix using the performance code
       a_pc->AddDomainIntegrator(new DiffusionIntegrator(one));
@@ -412,7 +469,7 @@ int main(int argc, char *argv[])
       a_pc->Assemble();
       a_pc->FormSystemMatrix(ess_tdof_list, A_pc);
    }
-   else if (pc_choice == HO)
+   else if (pc_choice == PCType::HO)
    {
       if (!matrix_free)
       {
@@ -441,7 +498,7 @@ int main(int argc, char *argv[])
    HypreSolver *amg = NULL;
 
    pcg->SetOperator(*a_oper);
-   if (pc_choice != NONE)
+   if (pc_choice != PCType::NONE)
    {
       amg = new HypreBoomerAMG(A_pc);
       pcg->SetPreconditioner(*amg);

miniapps/shifted/CMakeLists.txt

@@ -17,6 +17,7 @@ if (MFEM_USE_MPI)
   list(APPEND DIST_COMMON_HEADERS
         dist_solver.hpp
         sbm_solver.hpp
+	sbm_aux.hpp
         marking.hpp)
 
 convert_filenames_to_full_paths(DIST_COMMON_SOURCES)

miniapps/shifted/diffusion.cpp

@@ -604,7 +604,7 @@ int main(int argc, char *argv[])
    if (dirichlet_level_set_type == 2 || dirichlet_level_set_type == 3 ||
        (dirichlet_level_set_type == -1 && neumann_level_set_type == 2))
    {
-      ParGridFunction err(x);
+      ParGridFunction error(x);
       Vector pxyz(dim);
       pxyz(0) = 0.;
       for (int i = 0; i < nodes_cnt; i++)
@@ -620,7 +620,7 @@ int main(int argc, char *argv[])
          {
             exact_val = dirichlet_velocity_xy_sinusoidal(pxyz);
          }
-         err(i) = std::fabs(x(i) - exact_val);
+         error(i) = std::fabs(x(i) - exact_val);
       }
 
       if (visualization)
@@ -628,7 +628,7 @@ int main(int argc, char *argv[])
          char vishost[] = "localhost";
          int  visport   = 19916, s = 350;
          socketstream sol_sock;
-         common::VisualizeField(sol_sock, vishost, visport, err,
+         common::VisualizeField(sol_sock, vishost, visport, error,
                                 "Error", 2*s, 0, s, s, "Rj");
       }
 

miniapps/shifted/dist_solver.cpp

@@ -211,13 +211,13 @@ void HeatDistanceSolver::ComputeScalarDistance(Coefficient &zero_level_set,
       a_n.RecoverFEMSolution(X, b, u_neumann);
       delete prec2;
 
-      for (int i = 0; i < diffused_source.Size(); i++)
+      for (int ii = 0; ii < diffused_source.Size(); ii++)
       {
          // This assumes that the magnitudes of the two solutions are somewhat
          // similar; otherwise one of the solutions would dominate and the BC
          // won't look correct. To avoid this, it's good to have the source
          // away from the boundary (i.e. have more resolution).
-         diffused_source(i) = 0.5 * (u_neumann(i) + u_dirichlet(i));
+         diffused_source(ii) = 0.5 * (u_neumann(ii) + u_dirichlet(ii));
       }
       source = diffused_source;
    }

miniapps/shifted/makefile

@@ -53,15 +53,15 @@ COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
 %: %.cpp
 %.o: %.cpp
 
-%.o: $(SRC)%.cpp $(wildcard $(SRC)%.hpp) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+%.o: $(SRC)%.cpp sbm_aux.hpp $(wildcard $(SRC)%.hpp) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
 
 all: $(MINIAPPS)
 
-distance: sbm_aux.hpp $(DISTANCE_OBJ)
+distance: $(DISTANCE_OBJ)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(DISTANCE_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
 
-diffusion: sbm_aux.hpp $(DIFFUSION_OBJ)
+diffusion: $(DIFFUSION_OBJ)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(DIFFUSION_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
 
 # Rule for building lib-common

miniapps/shifted/sbm_solver.cpp

@@ -233,7 +233,7 @@ void SBM2DirichletIntegrator::AssembleFaceMatrix(
    Vector D(vD->GetVDim());
    // Assemble: -< \nabla u.n, w >
    //           -< u + \nabla u.d + h.o.t, \nabla w.n>
-   //           -<alpha h^{-1} (u + \nabla u.d + h.o.t), w + \nabla w.d + h.o.t>
+   //           +<alpha h^{-1} (u + \nabla u.d + h.o.t), w + \nabla w.d + h.o.t>
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);
@@ -324,10 +324,9 @@ void SBM2DirichletIntegrator::AssembleFaceMatrix(
       w = ip.weight*alpha*hinvdx;
       // + <alpha * hinv * u + grad u.d + h.o.t, w + grad w.d + h.o.t> - Term 4
       AddMult_a_VVt(w, wrk, temp_elmat);
-
-      int offset = elem1f ? 0 : ndof1;
-      elmat.CopyMN(temp_elmat, offset, offset);
    } // p < ir->GetNPoints()
+   int offset = elem1f ? 0 : ndof1;
+   elmat.CopyMN(temp_elmat, offset, offset);
 
    for (int i = 0; i < dkphi_dxk.Size(); i++)
    {
@@ -436,8 +435,8 @@ void SBM2DirichletLFIntegrator::AssembleRHSElementVect(
       }
    }
 
-   temp_elvect.SetSize(ndof);
-   temp_elvect = 0.0;
+   int offset = elem1f ? 0 : ndof1;
+   temp_elvect.SetDataAndSize(elvect.GetData()+offset, ndof);
 
    nor.SetSize(dim);
    nh.SetSize(dim);
@@ -541,10 +540,9 @@ void SBM2DirichletLFIntegrator::AssembleRHSElementVect(
    Vector D(vD->GetVDim());
    Vector wrk = shape;
    // Assemble: -< u_D, \nabla w.n >
-   //           -<alpha h^{-1} u_D, w + \nabla w.d + h.o.t>
+   //           +<alpha h^{-1} u_D, w + \nabla w.d + h.o.t>
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
-
       const IntegrationPoint &ip = ir->IntPoint(p);
 
       // Set the integration point in the face and the neighboring element
@@ -636,12 +634,6 @@ void SBM2DirichletLFIntegrator::AssembleRHSElementVect(
       wrk += dshape_dd; // \grad w .d
       wrk += q_hess_dot_d;
       temp_elvect.Add(w, wrk); // <u, gradw.d>
-
-      int offset = elem1f ? 0 : ndof1;
-      for (int i = 0; i < temp_elvect.Size(); i++)
-      {
-         elvect(i+offset) = temp_elvect(i);
-      }
    }
 
    for (int i = 0; i < dkphi_dxk.Size(); i++)
@@ -927,9 +919,9 @@ void SBM2NeumannIntegrator::AssembleFaceMatrix(
       wrk *= ip.weight * n_dot_ntilde;
 
       AddMult_a_VWt(1., shape, wrk, temp_elmat);
-      int offset = elem1f ? 0 : ndof1;
-      elmat.CopyMN(temp_elmat, offset, offset);
    } //p < ir->GetNPoints()
+   int offset = elem1f ? 0 : ndof1;
+   elmat.CopyMN(temp_elmat, offset, offset);
 
    for (int i = 0; i < dkphi_dxk.Size(); i++)
    {
@@ -1036,8 +1028,8 @@ void SBM2NeumannLFIntegrator::AssembleRHSElementVect(
       }
    }
 
-   temp_elvect.SetSize(ndof);
-   temp_elvect = 0.0;
+   int offset = elem1f ? 0 : ndof1;
+   temp_elvect.SetDataAndSize(elvect.GetData()+offset, ndof);
 
    nor.SetSize(dim);
    shape.SetSize(ndof);
@@ -1094,12 +1086,6 @@ void SBM2NeumannLFIntegrator::AssembleRHSElementVect(
       wrk.Set(n_dot_ntilde*w, shape);
       //<w, (nhat.n)t_n)
       temp_elvect.Add(1., wrk);
-
-      int offset = elem1f ? 0 : ndof1;
-      for (int i = 0; i < temp_elvect.Size(); i++)
-      {
-         elvect(i+offset) = temp_elvect(i);
-      }
    }
 }
 

miniapps/shifted/sbm_solver.hpp

@@ -78,7 +78,7 @@ public:
 /// method.
 /// A(u, w) = -<nabla u.n, w>
 ///           -<u + nabla u.d + h.o.t, nabla w.n>
-///           -<alpha h^{-1} (u + nabla u.d + h.o.t), w + nabla w.d + h.o.t>
+///           +<alpha h^{-1} (u + nabla u.d + h.o.t), w + nabla w.d + h.o.t>
 /// where h.o.t include higher-order derivatives (nabla^k u) due to Taylor
 /// expansion. Since this interior face integrator is applied to the surrogate
 /// boundary (see marking.hpp for notes on how the surrogate faces are
@@ -133,7 +133,7 @@ public:
 /// LinearFormIntegrator for the high-order extension of shifted boundary
 /// method.
 /// (u, w) = -<u_D, nabla w.n >
-///          -<alpha h^{-1} u_D, w + nabla w.d + h.o.t>
+///          +<alpha h^{-1} u_D, w + nabla w.d + h.o.t>
 /// where h.o.t include higher-order derivatives (nabla^k u) due to Taylor
 /// expansion. Since this interior face integrator is applied to the surrogate
 /// boundary (see marking.hpp for notes on how the surrogate faces are

miniapps/tools/convert-dc.cpp

@@ -157,15 +157,15 @@ int main(int argc, char *argv[])
    }
    args.PrintOptions(mfem::out);
 
-   DataCollection *src = create_data_collection(std::string(src_coll_name),
-                                                std::string(src_coll_type));
+   DataCollection *src_dc = create_data_collection(std::string(src_coll_name),
+                                                   std::string(src_coll_type));
 
-   DataCollection *out = create_data_collection(std::string(out_coll_name),
-                                                std::string(out_coll_type));
+   DataCollection *out_dc = create_data_collection(std::string(out_coll_name),
+                                                   std::string(out_coll_type));
 
-   src->Load(src_cycle);
+   src_dc->Load(src_cycle);
 
-   if (src->Error() != DataCollection::NO_ERROR)
+   if (src_dc->Error() != DataCollection::NO_ERROR)
    {
       mfem::out << "Error loading data collection: "
                 << src_coll_name
@@ -176,33 +176,33 @@ int main(int argc, char *argv[])
       return 1;
    }
 
-   out->SetOwnData(false);
+   out_dc->SetOwnData(false);
 
    // add mesh from source dc to output dc
 #ifdef MFEM_USE_MPI
-   out->SetMesh(MPI_COMM_WORLD,src->GetMesh());
+   out_dc->SetMesh(MPI_COMM_WORLD,src_dc->GetMesh());
 #else
-   out->SetMesh(src->GetMesh());
+   out_dc->SetMesh(src_dc->GetMesh());
 #endif
 
    // propagate the basics
-   out->SetCycle(src->GetCycle());
-   out->SetTime(src->GetTime());
-   out->SetTimeStep(src->GetTimeStep());
+   out_dc->SetCycle(src_dc->GetCycle());
+   out_dc->SetTime(src_dc->GetTime());
+   out_dc->SetTimeStep(src_dc->GetTimeStep());
 
    // loop over all fields in the source dc, and add them to the output dc
-   const DataCollection::FieldMapType &src_fields = src->GetFieldMap();
+   const DataCollection::FieldMapType &src_fields = src_dc->GetFieldMap();
 
    for (DataCollection::FieldMapType::const_iterator it = src_fields.begin();
         it != src_fields.end();
         ++it)
    {
-      out->RegisterField(it->first,it->second);
+      out_dc->RegisterField(it->first,it->second);
    }
 
-   out->Save();
+   out_dc->Save();
 
-   if (out->Error() != DataCollection::NO_ERROR)
+   if (out_dc->Error() != DataCollection::NO_ERROR)
    {
       mfem::out << "Error saving data collection: "
                 << out_coll_name
@@ -214,8 +214,8 @@ int main(int argc, char *argv[])
    }
 
    // cleanup
-   delete src;
-   delete out;
+   delete src_dc;
+   delete out_dc;
 
    return 0;
 }

miniapps/toys/lissajous.cpp

@@ -148,7 +148,7 @@ double u_function(const Vector &x)
 
 // Tubular Lissajous curve with the given parameters (a, b, theta)
 void lissajous_trans(const Vector &x, Vector &p,
-                     double a, double b, double delta)
+                     double a_, double b_, double delta_)
 {
    p.SetSize(3);
 
@@ -156,13 +156,13 @@ void lissajous_trans(const Vector &x, Vector &p,
    double theta = x[1];
    double t = phi;
 
-   double A = b; // Scaling of the curve along the x-axis
-   double B = a; // Scaling of the curve along the y-axis
+   double A = b_; // Scaling of the curve along the x-axis
+   double B = a_; // Scaling of the curve along the y-axis
 
    // Lissajous curve on a 3D cylinder
-   p[0] = B*cos(b*t);
-   p[1] = B*sin(b*t); // Y
-   p[2] = A*sin(a*t + delta); // X
+   p[0] = B*cos(b_*t);
+   p[1] = B*sin(b_*t); // Y
+   p[2] = A*sin(a_*t + delta_); // X
 
    // Turn the curve into a tubular surface
    {
@@ -170,14 +170,14 @@ void lissajous_trans(const Vector &x, Vector &p,
       double R = 0.02*(A+B);
 
       // normal to the cylinder at p(t)
-      double normal[3] = { cos(b*t), sin(b*t), 0 };
+      double normal[3] = { cos(b_*t), sin(b_*t), 0 };
 
       // tangent to the curve, dp/dt(t)
-      // double tangent[3] = { -b*B*sin(b*t), b*B*cos(b*t), A*a*cos(a*t+delta) };
+      // double tangent[3] = { -b_*B*sin(b_*t), b_*B*cos(b_*t), A*a_*cos(a_*t+delta_) };
 
       // normalized cross product of tangent and normal at p(t)
       double cn = 1e-128;
-      double cross[3] = { A*a*sin(b*t)*cos(a*t+delta), -A*a*cos(b*t)*cos(a*t+delta), b*B };
+      double cross[3] = { A*a_*sin(b_*t)*cos(a_*t+delta_), -A*a_*cos(b_*t)*cos(a_*t+delta_), b_*B };
       for (int i = 0; i < 3; i++) { cn += cross[i]*cross[i]; }
       for (int i = 0; i < 3; i++) { cross[i] /= sqrt(cn); }
 

miniapps/toys/mandel.cpp

@@ -36,15 +36,15 @@ using namespace std;
 //
 // In this particular miniapp, the material value is based on the number of
 // iterations for the point from the definition of the Mandelbrot set.
-int material(Vector &x, Vector &xmin, Vector &xmax)
+int material(Vector &p, Vector &pmin, Vector &pmax)
 {
    // Rescaling to [0,1]^sdim
-   for (int i = 0; i < x.Size(); i++)
+   for (int i = 0; i < p.Size(); i++)
    {
-      x(i) = (x(i)-xmin(i))/(xmax(i)-xmin(i));
+      p(i) = (p(i)-pmin(i))/(pmax(i)-pmin(i));
    }
-   x(0) -= 0.1;
-   double col = x(0), row = x(1);
+   p(0) -= 0.1;
+   double col = p(0), row = p(1);
    {
       int width = 1080, height = 1080;
       col *= width;
@@ -131,14 +131,14 @@ int main(int argc, char *argv[])
    for (int iter = 0; 1; iter++)
    {
       Array<Refinement> refs;
-      for (int i = 0; i < mesh.GetNE(); i++)
+      for (int e = 0; e < mesh.GetNE(); e++)
       {
          bool refine = false;
 
          // Sample materials in each element using "sd" sub-divisions
          Vector pt;
-         Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-         ElementTransformation *T = mesh.GetElementTransformation(i);
+         Geometry::Type geom = mesh.GetElementBaseGeometry(e);
+         ElementTransformation *T = mesh.GetElementTransformation(e);
          RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
          IntegrationRule &ir = RefG->RefPts;
 
@@ -160,8 +160,8 @@ int main(int argc, char *argv[])
          }
 
          // Set the element attribute as the "average". Other choices are
-         // possible here too, e.g. attr(i) = mat;
-         attr(i) = round(matsum/ir.GetNPoints());
+         // possible here too, e.g. attr(e) = mat;
+         attr(e) = round(matsum/ir.GetNPoints());
 
          // Mark the element for refinement
          if (refine)
@@ -200,7 +200,7 @@ int main(int argc, char *argv[])
                if (!type) { type = 7; } // because of tol
             }
 
-            refs.Append(Refinement(i, type));
+            refs.Append(Refinement(e, type));
          }
       }
 

miniapps/toys/mondrian.cpp

@@ -180,22 +180,22 @@ int main(int argc, char *argv[])
                const int s = sd+1;
                if (dim == 2)
                {
-                  for (int j = 0; j <= sd; j++)
-                     for (int i = 0; i < sd; i++)
+                  for (int jj = 0; jj <= sd; jj++)
+                     for (int ii = 0; ii < sd; ii++)
                      {
-                        dx += abs(mat[j*s + i+1] - mat[j*s + i]);
-                        dy += abs(mat[(i+1)*s + j] - mat[i*s + j]);
+                        dx += abs(mat[jj*s + ii+1] - mat[jj*s + ii]);
+                        dy += abs(mat[(ii+1)*s + jj] - mat[ii*s + jj]);
                      }
                }
                else if (dim == 3)
                {
-                  for (int k = 0; k <= sd; k++)
-                     for (int j = 0; j <= sd; j++)
-                        for (int i = 0; i < sd; i++)
+                  for (int kk = 0; kk <= sd; kk++)
+                     for (int jj = 0; jj <= sd; jj++)
+                        for (int ii = 0; ii < sd; ii++)
                         {
-                           dx += abs(mat[(k*s + j)*s + i+1] - mat[(k*s + j)*s + i]);
-                           dy += abs(mat[(k*s + i+1)*s + j] - mat[(k*s + i)*s + j]);
-                           dz += abs(mat[((i+1)*s + j)*s + k] - mat[(i*s + j)*s + k]);
+                           dx += abs(mat[(kk*s + jj)*s + ii+1] - mat[(kk*s + jj)*s + ii]);
+                           dy += abs(mat[(kk*s + ii+1)*s + jj] - mat[(kk*s + ii)*s + jj]);
+                           dz += abs(mat[((ii+1)*s + jj)*s + kk] - mat[(ii*s + jj)*s + kk]);
                         }
                }
                type = 0;

miniapps/toys/rubik.cpp

@@ -873,30 +873,30 @@ update_state(char axis, int tier, int incr)
 }
 
 void
-print_state(ostream & out)
+print_state(ostream & out_stream)
 {
-   out << "Rubik's Cube State:\n";
-   out << "  Centers: ";
+   out_stream << "Rubik's Cube State:\n";
+   out_stream << "  Centers: ";
    for (int i=0; i<6; i++)
    {
-      out << " " << rubik.cent_[i];
+      out_stream << " " << rubik.cent_[i];
    }
-   out << "\n";
-   out << "  Edges:   ";
+   out_stream << "\n";
+   out_stream << "  Edges:   ";
    for (int i=0; i<12; i++)
    {
-      out << " " << rubik.edge_[2 * i + 0]
-          << ":" << rubik.edge_[2 * i + 1];
+      out_stream << " " << rubik.edge_[2 * i + 0]
+                 << ":" << rubik.edge_[2 * i + 1];
    }
-   out << "\n";
-   out << "  Corners: ";
+   out_stream << "\n";
+   out_stream << "  Corners: ";
    for (int i=0; i<8; i++)
    {
-      out << " " << rubik.corn_[3 * i + 0]
-          << ":" << rubik.corn_[3 * i + 1]
-          << ":" << rubik.corn_[3 * i + 2];
+      out_stream << " " << rubik.corn_[3 * i + 0]
+                 << ":" << rubik.corn_[3 * i + 1]
+                 << ":" << rubik.corn_[3 * i + 2];
    }
-   out << "\n";
+   out_stream << "\n";
 }
 
 void repaint_cube(Mesh & mesh, GridFunction & color, socketstream & sock)

miniapps/toys/snake.cpp

@@ -366,12 +366,12 @@ rotate(double * x)
 }
 
 void
-trans(const int * conf, Mesh & mesh)
+trans(const int * new_conf, Mesh & mesh)
 {
    for (int i=0; i<23; i++)
    {
       joint_ = i;
-      notch_ = conf[i];
+      notch_ = new_conf[i];
 
       if (notch_ != 0)
       {
@@ -444,12 +444,12 @@ rotate_step(double * x)
 }
 
 bool
-anim_step(const int * conf, Mesh & mesh)
+anim_step(const int * new_conf, Mesh & mesh)
 {
    if (notch_ == 2 && step_ == 2 * nstep_) { joint_++; step_ = 0; }
    if (notch_ != 2 && step_ == nstep_) { joint_++; step_ = 0; }
    if (joint_ == 23) { return false; }
-   notch_ = conf[joint_];
+   notch_ = new_conf[joint_];
 
    if (notch_ == 0)
    {

tests/unit/CMakeLists.txt

@@ -70,6 +70,7 @@ set(UNIT_TESTS_SRCS
   fem/test_lexicographic_ordering.cpp
   fem/test_lin_interp.cpp
   fem/test_linear_fes.cpp
+  fem/test_lor.cpp
   fem/test_operatorjacobismoother.cpp
   fem/test_pa_coeff.cpp
   fem/test_pa_grad.cpp