Removing unit tests of Full operations in parallel since these are not supported (apparently...)

Positive diagonal in Jacobi smoothers

.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,14 @@
 
 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.
+
+- 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,7 +75,7 @@ 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.
 

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

fem/bilinearform.cpp

@@ -1766,8 +1766,13 @@ void MixedBilinearForm::FormRectangularSystemMatrix(
    mat_e = new SparseMatrix(mat->Height(), mat->Width());
    mat->EliminateCols(ess_trial_tdof_marker, *mat_e);
 
+   Array<int> cols;
+   Vector srow;
+
    for (int i=0; i<test_tdof_list.Size(); ++i)
    {
+      mat->GetRow(test_tdof_list[i], cols, srow);
+      mat_e->AddRow(test_tdof_list[i], cols, srow);
       mat->EliminateRow(test_tdof_list[i]);
    }
    mat_e->Finalize();

fem/bilinearform.hpp

@@ -276,9 +276,19 @@ public:
 
    /** @brief Add the original uneliminated matrix vector multiple to a vector.
        The original matrix is \f$ M + Me \f$ so we have:
-       \f$ y += M x + M_e x \f$ */
-   void FullAddMult(const Vector &x, Vector &y) const
-   { mat->AddMult(x, y); mat_e->AddMult(x, y); }
+       \f$ y += a * M x + a * M_e x \f$ */
+   void FullAddMult(const Vector &x, Vector &y, const double a = 1.0) const
+   { mat->AddMult(x, y, a); mat_e->AddMult(x, y, a); }
+
+   /// Matrix transpose vector multiplication:  \f$ y = M^T x \f$
+   virtual void MultTranspose(const Vector & x, Vector & y) const
+   { y = 0.0; AddMultTranspose (x, y); }
+
+   /** @brief Matrix transpose vector multiplication with the original
+       uneliminated matrix.  The original matrix is \f$ M + M_e \f$ so we have:
+       \f$ y = M^T x + {M_e}^T x \f$ */
+   void FullMultTranspose(const Vector &x, Vector &y) const
+   { mat->MultTranspose(x, y); mat_e->AddMultTranspose(x, y); }
 
    /// Add the matrix transpose vector multiplication:  \f$ y += a M^T x \f$
    virtual void AddMultTranspose(const Vector & x, Vector & y,
@@ -287,18 +297,19 @@ public:
 
    /** @brief Add the original uneliminated matrix transpose vector
        multiple to a vector. The original matrix is \f$ M + M_e \f$
-       so we have: \f$ y += M^T x + {M_e}^T x \f$ */
-   void FullAddMultTranspose(const Vector & x, Vector & y) const
-   { mat->AddMultTranspose(x, y); mat_e->AddMultTranspose(x, y); }
-
-   /// Matrix transpose vector multiplication:  \f$ y = M^T x \f$
-   virtual void MultTranspose(const Vector & x, Vector & y) const
-   { y = 0.0; AddMultTranspose (x, y); }
+       so we have: \f$ y += a * M^T x + a * {M_e}^T x \f$ */
+   void FullAddMultTranspose(const Vector & x, Vector & y,
+                             const double a = 1.0) const
+   { mat->AddMultTranspose(x, y, a); mat_e->AddMultTranspose(x, y, a); }
 
    /// Compute \f$ y^T M x \f$
    double InnerProduct(const Vector &x, const Vector &y) const
    { return mat->InnerProduct (x, y); }
 
+   /// Compute inner product for full uneliminated matrix \f$ y^T M x + y^T M_e x \f$
+   double FullInnerProduct(const Vector &x, const Vector &y) const
+   { return mat->InnerProduct(x, y) + mat_e->InnerProduct(x, y); }
+
    /// Returns a pointer to (approximation) of the matrix inverse:  \f$ M^{-1} \f$
    virtual MatrixInverse *Inverse() const;
 
@@ -434,8 +445,14 @@ public:
        recovered by calling RecoverFEMSolution() (with the same vectors @a X,
        @a b, and @a x).
 
-       NOTE: If there are no transformations, @a X simply reuses the data of
-             @a x. */
+       @note If there are no transformations, @a X simply reuses the data of
+             @a x.
+
+       @note This method does modify the bilinear form operator. For example,
+             calls to Mult() will produce different results before and after
+             use of this method. Use FullMult() to obtain the original behavior.
+             Similar methods exist for AddMult(), MultTranspose(), etc..
+   */
    virtual void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
                                  Vector &b, OperatorHandle &A, Vector &X,
                                  Vector &B, int copy_interior = 0);
@@ -590,10 +607,6 @@ public:
    void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x,
                             Vector &b);
 
-   /// Compute inner product for full uneliminated matrix \f$ y^T M x + y^T M_e x \f$
-   double FullInnerProduct(const Vector &x, const Vector &y) const
-   { return mat->InnerProduct(x, y) + mat_e->InnerProduct(x, y); }
-
    /// Update the @a FiniteElementSpace and delete all data associated with the old one.
    virtual void Update(FiniteElementSpace *nfes = NULL);
 
@@ -714,13 +727,42 @@ public:
    /// Matrix multiplication: \f$ y = M x \f$
    virtual void Mult(const Vector & x, Vector & y) const;
 
+   /** @brief Matrix vector multiplication with the original uneliminated
+       matrix.  The original matrix is \f$ M + M_e \f$ so we have:
+       \f$ y = M x + M_e x \f$ */
+   void FullMult(const Vector &x, Vector &y) const
+   { mat->Mult(x, y); mat_e->AddMult(x, y); }
+
+   /// Add the matrix vector multiple to a vector:  \f$ y += a M x \f$
    virtual void AddMult(const Vector & x, Vector & y,
                         const double a = 1.0) const;
 
+   /** @brief Add the original uneliminated matrix vector multiple to a vector.
+       The original matrix is \f$ M + Me \f$ so we have:
+       \f$ y += a * M x + a * M_e x \f$ */
+   void FullAddMult(const Vector &x, Vector &y, const double a = 1.0) const
+   { mat->AddMult(x, y, a); mat_e->AddMult(x, y, a); }
+
+   /// Matrix transpose vector multiplication:  \f$ y = M^T x \f$
    virtual void MultTranspose(const Vector & x, Vector & y) const;
+
+   /** @brief Matrix transpose vector multiplication with the original
+       uneliminated matrix.  The original matrix is \f$ M + M_e \f$ so we have:
+       \f$ y = M^T x + {M_e}^T x \f$ */
+   void FullMultTranspose(const Vector &x, Vector &y) const
+   { mat->MultTranspose(x, y); mat_e->AddMultTranspose(x, y); }
+
+   /// Add the matrix transpose vector multiplication:  \f$ y += a M^T x \f$
    virtual void AddMultTranspose(const Vector & x, Vector & y,
                                  const double a = 1.0) const;
 
+   /** @brief Add the original uneliminated matrix transpose vector
+       multiple to a vector. The original matrix is \f$ M + M_e \f$
+       so we have: \f$ y += a * M^T x + a * {M_e}^T x \f$ */
+   void FullAddMultTranspose(const Vector & x, Vector & y,
+                             const double a = 1.0) const
+   { mat->AddMultTranspose(x, y, a); mat_e->AddMultTranspose(x, y, a); }
+
    virtual MatrixInverse *Inverse() const;
 
    /// Finalizes the matrix initialization.
@@ -904,12 +946,17 @@ public:
       A.MakeRef(*A_ptr);
    }
 
-   /** @brief Form the linear system A X = B, corresponding to this mixed bilinear
-       form and the linear form @a b(.).
+   /** @brief Form the linear system A X = B, corresponding to this mixed
+       bilinear form and the linear form @a b(.). */
+   /** Return in @a A a *reference* to the system matrix that is
+       column-constrained. The reference will be invalidated when
+       SetOperatorType(), Update(), or the destructor is called.
 
-       Return in @a A a *reference* to the system matrix that is column-constrained.
-       The reference will be invalidated when SetOperatorType(), Update(), or the
-       destructor is called. */
+       @note This method does modify the bilinear form operator. For example,
+             calls to Mult() will produce different results before and after
+             use of this method. Use FullMult() to obtain the original behavior.
+             Similar methods exist for AddMult(), MultTranspose(), etc..
+   */
    virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
                                             const Array<int> &test_tdof_list,
                                             Vector &x, Vector &b,

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

linalg/solvers.cpp

@@ -199,7 +199,16 @@ void OperatorJacobiSmoother::Setup(const Vector &diag)
    const double delta = damping;
    auto D = diag.Read();
    auto DI = dinv.Write();
-   MFEM_FORALL(i, height, DI[i] = delta / D[i]; );
+   const bool use_abs_diag_ = use_abs_diag;
+   MFEM_FORALL(i, height,
+   {
+      if (D[i] == 0.0)
+      {
+         MFEM_ABORT_KERNEL("Zero diagonal entry in OperatorJacobiSmoother");
+      }
+      if (!use_abs_diag_) { DI[i] = delta / D[i]; }
+      else                { DI[i] = delta / std::abs(D[i]); }
+   });
    if (ess_tdof_list && ess_tdof_list->Size() > 0)
    {
       auto I = ess_tdof_list->Read();
@@ -229,7 +238,10 @@ void OperatorJacobiSmoother::Mult(const Vector &x, Vector &y) const
    auto DI = dinv.Read();
    auto R = residual.Read();
    auto Y = y.ReadWrite();
-   MFEM_FORALL(i, height, Y[i] += DI[i] * R[i]; );
+   MFEM_FORALL(i, height,
+   {
+      Y[i] += DI[i] * R[i];
+   });
 }
 
 OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,

linalg/solvers.hpp

@@ -163,6 +163,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 +187,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;
 };

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/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/meshing/mesh-optimizer.cpp

@@ -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.
@@ -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/pmesh-optimizer.cpp

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

tests/unit/CMakeLists.txt

@@ -47,6 +47,7 @@ set(UNIT_TESTS_SRCS
   mesh/test_pmesh.cpp
   mesh/test_periodic_mesh.cpp
   mesh/test_vtu.cpp
+  fem/common_get_mesh.cpp
   fem/test_1d_bilininteg.cpp
   fem/test_2d_bilininteg.cpp
   fem/test_3d_bilininteg.cpp
@@ -70,6 +71,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

tests/unit/fem/common_get_mesh.cpp

@@ -0,0 +1,416 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "common_get_mesh.hpp"
+
+using namespace mfem;
+
+namespace mfem_test_fem
+{
+
+Mesh * GetMesh(MeshType type, double lx, double ly, double lz)
+{
+   Mesh * mesh = NULL;
+   double c[3];
+   int    v[8];
+
+   switch (type)
+   {
+      case SEGMENT:
+         mesh = new Mesh(1, 2, 1);
+         c[0] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx;
+         mesh->AddVertex(c);
+         v[0] = 0; v[1] = 1;
+         mesh->AddSegment(v);
+         {
+            Element * el = mesh->NewElement(Geometry::POINT);
+            el->SetAttribute(1);
+            el->SetVertices(&v[0]);
+            mesh->AddBdrElement(el);
+         }
+         {
+            Element * el = mesh->NewElement(Geometry::POINT);
+            el->SetAttribute(2);
+            el->SetVertices(&v[1]);
+            mesh->AddBdrElement(el);
+         }
+         break;
+      case QUADRILATERAL:
+         mesh = new Mesh(2, 4, 1);
+         c[0] = 0.0; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
+         mesh->AddQuad(v);
+         break;
+      case TRIANGLE2A:
+         mesh = new Mesh(2, 4, 2);
+         c[0] = 0.0; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 2;
+         mesh->AddTri(v);
+         v[0] = 2; v[1] = 3; v[2] = 0;
+         mesh->AddTri(v);
+         break;
+      case TRIANGLE2B:
+         mesh = new Mesh(2, 4, 2);
+         c[0] = 0.0; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly;
+         mesh->AddVertex(c);
+
+         v[0] = 1; v[1] = 2; v[2] = 0;
+         mesh->AddTri(v);
+         v[0] = 3; v[1] = 0; v[2] = 2;
+         mesh->AddTri(v);
+         break;
+      case TRIANGLE2C:
+         mesh = new Mesh(2, 4, 2);
+         c[0] = 0.0; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly;
+         mesh->AddVertex(c);
+
+         v[0] = 2; v[1] = 0; v[2] = 1;
+         mesh->AddTri(v);
+         v[0] = 0; v[1] = 2; v[2] = 3;
+         mesh->AddTri(v);
+         break;
+      case TRIANGLE4:
+         mesh = new Mesh(2, 5, 4);
+         c[0] = 0.0; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lx; c[1] = 0.5 * ly;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 4;
+         mesh->AddTri(v);
+         v[0] = 1; v[1] = 2; v[2] = 4;
+         mesh->AddTri(v);
+         v[0] = 2; v[1] = 3; v[2] = 4;
+         mesh->AddTri(v);
+         v[0] = 3; v[1] = 0; v[2] = 4;
+         mesh->AddTri(v);
+         break;
+      case MIXED2D:
+         mesh = new Mesh(2, 6, 4);
+         c[0] = 0.0; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * ly; c[1] = 0.5 * ly;
+         mesh->AddVertex(c);
+         c[0] = lx - 0.5 * ly; c[1] = 0.5 * ly;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 5; v[3] = 4;
+         mesh->AddQuad(v);
+         v[0] = 1; v[1] = 2; v[2] = 5;
+         mesh->AddTri(v);
+         v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 5;
+         mesh->AddQuad(v);
+         v[0] = 3; v[1] = 0; v[2] = 4;
+         mesh->AddTri(v);
+         break;
+      case HEXAHEDRON:
+         mesh = new Mesh(3, 8, 1);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
+         v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
+         mesh->AddHex(v);
+         break;
+      case HEXAHEDRON2A:
+      case HEXAHEDRON2B:
+      case HEXAHEDRON2C:
+      case HEXAHEDRON2D:
+         mesh = new Mesh(3, 12, 2);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 5; v[2] = 11; v[3] = 6;
+         v[4] = 1; v[5] = 4; v[6] = 10; v[7] = 7;
+         mesh->AddHex(v);
+
+         switch (type)
+         {
+            case HEXAHEDRON2A: // Face Orientation 1
+               v[0] = 4; v[1] = 10; v[2] = 7; v[3] = 1;
+               v[4] = 3; v[5] = 9; v[6] = 8; v[7] = 2;
+               mesh->AddHex(v);
+               break;
+            case HEXAHEDRON2B: // Face Orientation 3
+               v[0] = 10; v[1] = 7; v[2] = 1; v[3] = 4;
+               v[4] = 9; v[5] = 8; v[6] = 2; v[7] = 3;
+               mesh->AddHex(v);
+               break;
+            case HEXAHEDRON2C: // Face Orientation 5
+               v[0] = 7; v[1] = 1; v[2] = 4; v[3] = 10;
+               v[4] = 8; v[5] = 2; v[6] = 3; v[7] = 9;
+               mesh->AddHex(v);
+               break;
+            case HEXAHEDRON2D: // Face Orientation 7
+               v[0] = 1; v[1] = 4; v[2] = 10; v[3] = 7;
+               v[4] = 2; v[5] = 3; v[6] = 9; v[7] = 8;
+               mesh->AddHex(v);
+               break;
+            default:
+               // Cannot happen
+               break;
+         }
+         break;
+      case WEDGE2:
+         mesh = new Mesh(3, 8, 2);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 4; v[4] = 5; v[5] = 6;
+         mesh->AddWedge(v);
+         v[0] = 0; v[1] = 2; v[2] = 3; v[3] = 4; v[4] = 6; v[5] = 7;
+         mesh->AddWedge(v);
+         break;
+      case TETRAHEDRA:
+         mesh = new Mesh(3, 8, 5);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 2; v[2] = 7; v[3] = 5;
+         mesh->AddTet(v);
+         v[0] = 6; v[1] = 7; v[2] = 2; v[3] = 5;
+         mesh->AddTet(v);
+         v[0] = 4; v[1] = 7; v[2] = 5; v[3] = 0;
+         mesh->AddTet(v);
+         v[0] = 1; v[1] = 0; v[2] = 5; v[3] = 2;
+         mesh->AddTet(v);
+         v[0] = 3; v[1] = 7; v[2] = 0; v[3] = 2;
+         mesh->AddTet(v);
+         break;
+      case WEDGE4:
+         mesh = new Mesh(3, 10, 4);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lx; c[1] = 0.5 * ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lx; c[1] = 0.5 * ly; c[2] = lz;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 4; v[3] = 5; v[4] = 6; v[5] = 9;
+         mesh->AddWedge(v);
+         v[0] = 1; v[1] = 2; v[2] = 4; v[3] = 6; v[4] = 7; v[5] = 9;
+         mesh->AddWedge(v);
+         v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 7; v[4] = 8; v[5] = 9;
+         mesh->AddWedge(v);
+         v[0] = 3; v[1] = 0; v[2] = 4; v[3] = 8; v[4] = 5; v[5] = 9;
+         mesh->AddWedge(v);
+         break;
+      case MIXED3D6:
+         mesh = new Mesh(3, 12, 6);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lz; c[1] = 0.5 * lz; c[2] = 0.5 * lz;
+         mesh->AddVertex(c);
+         c[0] = lx - 0.5 * lz; c[1] = 0.5 * lz; c[2] = 0.5 * lz;
+         mesh->AddVertex(c);
+         c[0] = lx - 0.5 * lz; c[1] = ly - 0.5 * lz; c[2] = 0.5 * lz;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * lz; c[1] = ly - 0.5 * lz; c[2] = 0.5 * lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
+         v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
+         mesh->AddHex(v);
+         v[0] = 0; v[1] = 4; v[2] = 8; v[3] = 1; v[4] = 5; v[5] = 9;
+         mesh->AddWedge(v);
+         v[0] = 1; v[1] = 5; v[2] = 9; v[3] = 2; v[4] = 6; v[5] = 10;
+         mesh->AddWedge(v);
+         v[0] = 2; v[1] = 6; v[2] = 10; v[3] = 3; v[4] = 7; v[5] = 11;
+         mesh->AddWedge(v);
+         v[0] = 3; v[1] = 7; v[2] = 11; v[3] = 0; v[4] = 4; v[5] = 8;
+         mesh->AddWedge(v);
+         v[0] = 4; v[1] = 5; v[2] = 6; v[3] = 7;
+         v[4] = 8; v[5] = 9; v[6] = 10; v[7] = 11;
+         mesh->AddHex(v);
+         break;
+      case MIXED3D8:
+         mesh = new Mesh(3, 10, 8);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = 0.0;
+         mesh->AddVertex(c);
+
+         c[0] = 0.25 * lx; c[1] = 0.5 * ly; c[2] = 0.5 * lz;
+         mesh->AddVertex(c);
+         c[0] = 0.75 * lx; c[1] = 0.5 * ly; c[2] = 0.5 * lz;
+         mesh->AddVertex(c);
+
+         c[0] = 0.0; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = 0.0; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = lx; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = ly; c[2] = lz;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 3; v[2] = 4; v[3] = 1; v[4] = 2; v[5] = 5;
+         mesh->AddWedge(v);
+         v[0] = 3; v[1] = 9; v[2] = 4; v[3] = 2; v[4] = 8; v[5] = 5;
+         mesh->AddWedge(v);
+         v[0] = 9; v[1] = 6; v[2] = 4; v[3] = 8; v[4] = 7; v[5] = 5;
+         mesh->AddWedge(v);
+         v[0] = 6; v[1] = 0; v[2] = 4; v[3] = 7; v[4] = 1; v[5] = 5;
+         mesh->AddWedge(v);
+         v[0] = 0; v[1] = 3; v[2] = 9; v[3] = 4;
+         mesh->AddTet(v);
+         v[0] = 0; v[1] = 9; v[2] = 6; v[3] = 4;
+         mesh->AddTet(v);
+         v[0] = 1; v[1] = 7; v[2] = 2; v[3] = 5;
+         mesh->AddTet(v);
+         v[0] = 8; v[1] = 2; v[2] = 7; v[3] = 5;
+         mesh->AddTet(v);
+         break;
+   }
+   mesh->FinalizeTopology();
+
+   return mesh;
+}
+
+} // namespace mfem_test_fem

tests/unit/fem/common_get_mesh.hpp

@@ -0,0 +1,41 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+
+namespace mfem_test_fem
+{
+
+enum MeshType
+{
+   SEGMENT = 0,
+   QUADRILATERAL = 1,
+   TRIANGLE2A = 2,
+   TRIANGLE2B = 3,
+   TRIANGLE2C = 4,
+   TRIANGLE4 = 5,
+   MIXED2D = 6,
+   HEXAHEDRON = 7,
+   HEXAHEDRON2A = 8,
+   HEXAHEDRON2B = 9,
+   HEXAHEDRON2C = 10,
+   HEXAHEDRON2D = 11,
+   WEDGE2 = 12,
+   TETRAHEDRA = 13,
+   WEDGE4 = 14,
+   MIXED3D6 = 15,
+   MIXED3D8 = 16
+};
+
+mfem::Mesh * GetMesh(MeshType type,
+                     double lx = 1.0, double ly = 1.0, double lz = 1.0);
+
+}

tests/unit/fem/test_bilinearform.cpp

@@ -11,10 +11,19 @@
 
 #include "mfem.hpp"
 #include "unit_tests.hpp"
+#include "common_get_mesh.hpp"
 
 #include <iostream>
 
 using namespace mfem;
+using namespace mfem_test_fem;
+
+namespace bilinearform
+{
+
+static double a_ = 5.0;
+static double b_ = 3.0;
+static double c_ = 2.0;
 
 TEST_CASE("Test order of boundary integrators",
           "[BilinearForm]")
@@ -142,3 +151,293 @@ TEST_CASE("FormLinearSystem/SolutionScope",
       REQUIRE(AsConst(sol)(bdr_dof) == 0.0);
    }
 }
+
+enum FEType
+{
+   H1_FEC = 0,
+   ND_FEC,
+   RT_FEC,
+   L2V_FEC,
+   L2I_FEC,
+};
+
+TEST_CASE("BilinearForm Full Ops",
+          "[BilinearForm]")
+{
+   int order = 2;
+   double alpha = M_E;
+
+   for (int mt = (int)MeshType::SEGMENT;
+        mt <= (int)MeshType::MIXED3D8; mt++)
+   {
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
+      int  dim = mesh->Dimension();
+      mesh->UniformRefinement();
+
+      Vector oneVec(dim); oneVec = 1.0;
+
+      ConstantCoefficient oneCoef(1.0);
+      VectorConstantCoefficient oneVecCoef(oneVec);
+
+      for (int ft = (int)FEType::H1_FEC; ft <= (int)FEType::RT_FEC; ft++)
+      {
+         // if (ft == (int)FEType::ND_FEC || ft == (int)FEType::RT_FEC)
+         // { continue; }
+         bool vec = (ft == (int)FEType::ND_FEC || ft == (int)FEType::RT_FEC);
+
+         if (dim == 1 && vec) { continue; }
+         if (vec && (mt == (int)MeshType::WEDGE2 ||
+                     mt == (int)MeshType::WEDGE4 ||
+                     mt == (int)MeshType::MIXED3D6 ||
+                     mt == (int)MeshType::MIXED3D8))
+         { continue; }
+
+         SECTION("Integral of field " + std::to_string(ft) +
+                 " on mesh type " + std::to_string(mt) )
+         {
+
+            FiniteElementCollection *fec = NULL;
+            switch ((FEType)ft)
+            {
+               case FEType::H1_FEC:
+                  fec = new H1_FECollection(order, dim);
+                  break;
+               case FEType::ND_FEC:
+                  fec = new ND_FECollection(order, dim);
+                  break;
+               case FEType::RT_FEC:
+                  fec = new RT_FECollection(order-1, dim);
+                  break;
+               case FEType::L2V_FEC:
+                  fec = new L2_FECollection(order-1, dim);
+                  break;
+               case FEType::L2I_FEC:
+                  fec = new L2_FECollection(order, dim,
+                                            BasisType::GaussLegendre,
+                                            FiniteElement::INTEGRAL);
+                  break;
+               default:
+                  MFEM_ABORT("Invalid vector FE type");
+            }
+            FiniteElementSpace fespace(mesh, fec);
+
+            Array<int> ess_tdof_list;
+            if (mesh->bdr_attributes.Size())
+            {
+               Array<int> ess_bdr(mesh->bdr_attributes.Max());
+               ess_bdr = 1;
+               fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+            }
+
+            GridFunction u(&fespace);
+            if (!vec)
+            {
+               u.ProjectCoefficient(oneCoef);
+            }
+            else
+            {
+               u.ProjectCoefficient(oneVecCoef);
+            }
+
+            BilinearForm a(&fespace);
+            if (!vec)
+            {
+               a.AddDomainIntegrator(new MassIntegrator(oneCoef));
+            }
+            else
+            {
+               a.AddDomainIntegrator(new VectorFEMassIntegrator(oneCoef));
+            }
+            a.Assemble();
+
+            LinearForm Au(&fespace);
+            LinearForm ATu(&fespace);
+            LinearForm aAu(&fespace);
+            LinearForm aATu(&fespace);
+            LinearForm b(&fespace);
+
+            a.Mult(u, Au);
+            a.MultTranspose(u, ATu);
+
+            aAu = Au;
+            aATu = ATu;
+
+            a.AddMult(u, aAu, alpha);
+            a.AddMultTranspose(u, aATu, alpha);
+
+            // Modify the Bilinear Form
+            OperatorPtr A;
+            a.FormSystemMatrix(ess_tdof_list, A);
+
+            a.FullMult(u, b);
+            b -= Au;
+
+            REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
+
+            a.FullMultTranspose(u, b);
+            b -= ATu;
+
+            REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
+
+            b = Au;
+            a.FullAddMult(u, b, alpha);
+            b -= aAu;
+
+            REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
+
+            b = ATu;
+            a.FullAddMultTranspose(u, b, alpha);
+            b -= aATu;
+
+            REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
+
+            delete fec;
+         }
+      }
+
+      delete mesh;
+   }
+}
+
+TEST_CASE("MixedBilinearform Full Ops",
+          "[MixedBilinearForm]")
+{
+   int order = 2;
+   double alpha = M_E;
+
+   for (int mt = (int)MeshType::SEGMENT;
+        mt <= (int)MeshType::MIXED3D8; mt++)
+   {
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
+      int  dim = mesh->Dimension();
+      mesh->UniformRefinement();
+
+      Vector oneVec(dim); oneVec = 1.0;
+
+      ConstantCoefficient oneCoef(1.0);
+      VectorConstantCoefficient oneVecCoef(oneVec);
+
+      for (int ft = (int)FEType::H1_FEC; ft <= (int)FEType::RT_FEC; ft++)
+      {
+         bool vec = (ft == (int)FEType::ND_FEC || ft == (int)FEType::RT_FEC);
+
+         if (dim == 1 && vec) { continue; }
+         if (vec && (mt == (int)MeshType::WEDGE2 ||
+                     mt == (int)MeshType::WEDGE4 ||
+                     mt == (int)MeshType::MIXED3D6 ||
+                     mt == (int)MeshType::MIXED3D8))
+         { continue; }
+
+         SECTION("Integral of field " + std::to_string(ft) +
+                 " on mesh type " + std::to_string(mt) )
+         {
+
+            FiniteElementCollection *fec_dom = NULL;
+            FiniteElementCollection *fec_ran = NULL;
+            switch ((FEType)ft)
+            {
+               case FEType::H1_FEC:
+                  fec_dom = new H1_FECollection(order, dim);
+                  fec_ran = new H1_FECollection(order-1, dim);
+                  break;
+               case FEType::ND_FEC:
+                  fec_dom = new ND_FECollection(order, dim);
+                  fec_ran = new RT_FECollection(order-1, dim);
+                  break;
+               case FEType::RT_FEC:
+                  fec_dom = new RT_FECollection(order-1, dim);
+                  fec_ran = new ND_FECollection(order, dim);
+                  break;
+               default:
+                  MFEM_ABORT("Invalid vector FE type");
+            }
+            FiniteElementSpace fespace_dom(mesh, fec_dom);
+            FiniteElementSpace fespace_ran(mesh, fec_ran);
+
+            Array<int> ess_tdof_list_dom;
+            Array<int> ess_tdof_list_ran;
+            if (mesh->bdr_attributes.Size())
+            {
+               Array<int> ess_bdr(mesh->bdr_attributes.Max());
+               ess_bdr = 1;
+               fespace_dom.GetEssentialTrueDofs(ess_bdr, ess_tdof_list_dom);
+               fespace_ran.GetEssentialTrueDofs(ess_bdr, ess_tdof_list_ran);
+            }
+
+            GridFunction u_dom(&fespace_dom);
+            GridFunction u_ran(&fespace_ran);
+            if (!vec)
+            {
+               u_dom.ProjectCoefficient(oneCoef);
+               u_ran.ProjectCoefficient(oneCoef);
+            }
+            else
+            {
+               u_dom.ProjectCoefficient(oneVecCoef);
+               u_ran.ProjectCoefficient(oneVecCoef);
+            }
+
+            MixedBilinearForm a(&fespace_dom, &fespace_ran);
+            if (!vec)
+            {
+               a.AddDomainIntegrator(new MassIntegrator(oneCoef));
+            }
+            else
+            {
+               a.AddDomainIntegrator(new VectorFEMassIntegrator(oneCoef));
+            }
+            a.Assemble();
+
+            LinearForm Au(&fespace_ran);
+            LinearForm ATu(&fespace_dom);
+            LinearForm aAu(&fespace_ran);
+            LinearForm aATu(&fespace_dom);
+            LinearForm b_ran(&fespace_ran);
+            LinearForm b_dom(&fespace_dom);
+
+            a.Mult(u_dom, Au);
+            a.MultTranspose(u_ran, ATu);
+
+            aAu = Au;
+            aATu = ATu;
+
+            a.AddMult(u_dom, aAu, alpha);
+            a.AddMultTranspose(u_ran, aATu, alpha);
+
+            // Modify the Bilinear Form
+            OperatorPtr A;
+            a.FormRectangularSystemMatrix(ess_tdof_list_dom,
+                                          ess_tdof_list_ran, A);
+
+            a.FullMult(u_dom, b_ran);
+            b_ran -= Au;
+
+            REQUIRE(b_ran.Norml2() == MFEM_Approx( 0.0));
+
+            a.FullMultTranspose(u_ran, b_dom);
+            b_dom -= ATu;
+
+            REQUIRE(b_dom.Norml2() == MFEM_Approx( 0.0));
+
+            b_ran = Au;
+            a.FullAddMult(u_dom, b_ran, alpha);
+            b_ran -= aAu;
+
+            REQUIRE(b_ran.Norml2() == MFEM_Approx( 0.0));
+
+            b_dom = ATu;
+            a.FullAddMultTranspose(u_ran, b_dom, alpha);
+            b_dom -= aATu;
+
+            REQUIRE(b_dom.Norml2() == MFEM_Approx( 0.0));
+
+            delete fec_dom;
+            delete fec_ran;
+         }
+      }
+
+      delete mesh;
+   }
+}
+
+} // namespace bilinearform

tests/unit/fem/test_build_dof_to_arrays.cpp

@@ -11,6 +11,9 @@
 
 #include "mfem.hpp"
 #include "unit_tests.hpp"
+#include "common_get_mesh.hpp"
+
+using namespace mfem_test_fem;
 
 namespace mfem
 {
@@ -22,29 +25,6 @@ static double a_ = M_PI;
 static double b_ = M_PI / sqrt(2.0);
 static double c_ = M_PI / 2.0;
 
-enum MeshType
-{
-   SEGMENT = 0,
-   QUADRILATERAL = 1,
-   TRIANGLE2A = 2,
-   TRIANGLE2B = 3,
-   TRIANGLE2C = 4,
-   TRIANGLE4 = 5,
-   MIXED2D = 6,
-   HEXAHEDRON = 7,
-   HEXAHEDRON2A = 8,
-   HEXAHEDRON2B = 9,
-   HEXAHEDRON2C = 10,
-   HEXAHEDRON2D = 11,
-   WEDGE2 = 12,
-   TETRAHEDRA = 13,
-   WEDGE4 = 14,
-   MIXED3D6 = 15,
-   MIXED3D8 = 16
-};
-
-Mesh * GetMesh(MeshType type);
-
 enum BasisType
 {
    H1 = 0, ND = 1, RT = 2, L2 = 3
@@ -61,7 +41,7 @@ TEST_CASE("Build Dof To Arrays",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int  dim = mesh->Dimension();
       if (dim < 3 ||
           mt == MeshType::HEXAHEDRON ||
@@ -156,7 +136,7 @@ TEST_CASE("Build Dof To Arrays (Parallel)",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int  dim = mesh->Dimension();
       if (dim < 3 ||
           mt == MeshType::HEXAHEDRON ||
@@ -238,239 +218,6 @@ TEST_CASE("Build Dof To Arrays (Parallel)",
 }
 #endif // MFEM_USE_MPI
 
-Mesh * GetMesh(MeshType type)
-{
-   Mesh * mesh = NULL;
-
-   switch (type)
-   {
-      case SEGMENT:
-         mesh = new Mesh(1, 2, 1);
-         mesh->AddVertex(0.0);
-         mesh->AddVertex(a_);
-
-         mesh->AddSegment(0, 1);
-
-         mesh->AddBdrPoint(0);
-         mesh->AddBdrPoint(1);
-         break;
-      case QUADRILATERAL:
-         mesh = new Mesh(2, 4, 1);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddQuad(0, 1, 2, 3);
-         break;
-      case TRIANGLE2A:
-         mesh = new Mesh(2, 4, 2);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddTriangle(0, 1, 2);
-         mesh->AddTriangle(2, 3, 0);
-         break;
-      case TRIANGLE2B:
-         mesh = new Mesh(2, 4, 2);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddTriangle(1, 2, 0);
-         mesh->AddTriangle(3, 0, 2);
-         break;
-      case TRIANGLE2C:
-         mesh = new Mesh(2, 4, 2);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddTriangle(2, 0, 1);
-         mesh->AddTriangle(0, 2, 3);
-         break;
-      case TRIANGLE4:
-         mesh = new Mesh(2, 5, 4);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-         mesh->AddVertex(0.5 * a_, 0.5 * b_);
-
-         mesh->AddTriangle(0, 1, 4);
-         mesh->AddTriangle(1, 2, 4);
-         mesh->AddTriangle(2, 3, 4);
-         mesh->AddTriangle(3, 0, 4);
-         break;
-      case MIXED2D:
-         mesh = new Mesh(2, 6, 4);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-         mesh->AddVertex(0.5 * b_, 0.5 * b_);
-         mesh->AddVertex(a_ - 0.5 * b_, 0.5 * b_);
-
-         mesh->AddQuad(0, 1, 5, 4);
-         mesh->AddTriangle(1, 2, 5);
-         mesh->AddQuad(2, 3, 4, 5);
-         mesh->AddTriangle(3, 0, 4);
-         break;
-      case HEXAHEDRON:
-         mesh = new Mesh(3, 8, 1);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
-         break;
-      case HEXAHEDRON2A:
-      case HEXAHEDRON2B:
-      case HEXAHEDRON2C:
-      case HEXAHEDRON2D:
-         mesh = new Mesh(3, 12, 2);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(0.5 * a_, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.5 * a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(0.5 * a_, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.5 * a_, b_, c_);
-         mesh->AddVertex(0.0,b_, c_);
-
-         mesh->AddHex(0, 5, 11, 6, 1, 4, 10, 7);
-
-         switch (type)
-         {
-            case HEXAHEDRON2A: // Face Orientation 1
-               mesh->AddHex(4, 10, 7, 1, 3, 9, 8, 2);
-               break;
-            case HEXAHEDRON2B: // Face Orientation 3
-               mesh->AddHex(10, 7, 1, 4, 9, 8, 2, 3);
-               break;
-            case HEXAHEDRON2C: // Face Orientation 5
-               mesh->AddHex(7, 1, 4, 10, 8, 2, 3, 9);
-               break;
-            case HEXAHEDRON2D: // Face Orientation 7
-               mesh->AddHex(1, 4, 10, 7, 2, 3, 9, 8);
-               break;
-            default:
-               // Cannot happen
-               break;
-         }
-         break;
-      case WEDGE2:
-         mesh = new Mesh(3, 8, 2);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddWedge(0, 1, 2, 4, 5, 6);
-         mesh->AddWedge(0, 2, 3, 4, 6, 7);
-         break;
-      case TETRAHEDRA:
-         mesh = new Mesh(3, 8, 5);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddTet(0, 2, 7, 5);
-         mesh->AddTet(6, 7, 2, 5);
-         mesh->AddTet(4, 7, 5, 0);
-         mesh->AddTet(1, 0, 5, 2);
-         mesh->AddTet(3, 7, 0, 2);
-         break;
-      case WEDGE4:
-         mesh = new Mesh(3, 10, 4);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.5 * a_, 0.5 * b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-         mesh->AddVertex(0.5 * a_, 0.5 * b_, c_);
-
-         mesh->AddWedge(0, 1, 4, 5, 6, 9);
-         mesh->AddWedge(1, 2, 4, 6, 7, 9);
-         mesh->AddWedge(2, 3, 4, 7, 8, 9);
-         mesh->AddWedge(3, 0, 4, 8, 5, 9);
-         break;
-      case MIXED3D6:
-         mesh = new Mesh(3, 12, 6);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.5 * c_, 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(a_ - 0.5 * c_, 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(a_ - 0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
-         mesh->AddWedge(0, 4, 8, 1, 5, 9);
-         mesh->AddWedge(1, 5, 9, 2, 6, 10);
-         mesh->AddWedge(2, 6, 10, 3, 7, 11);
-         mesh->AddWedge(3, 7, 11, 0, 4, 8);
-         mesh->AddHex(4, 5, 6, 7, 8, 9, 10, 11);
-         break;
-      case MIXED3D8:
-         mesh = new Mesh(3, 10, 8);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.25 * a_, 0.5 * b_, 0.5 * c_);
-         mesh->AddVertex(0.75 * a_, 0.5 * b_, 0.5 * c_);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddWedge(0, 3, 4, 1, 2, 5);
-         mesh->AddWedge(3, 9, 4, 2, 8, 5);
-         mesh->AddWedge(9, 6, 4, 8, 7, 5);
-         mesh->AddWedge(6, 0, 4, 7, 1, 5);
-         mesh->AddTet(0, 3, 9, 4);
-         mesh->AddTet(0, 9, 6, 4);
-         mesh->AddTet(1, 7, 2, 5);
-         mesh->AddTet(8, 2, 7, 5);
-         break;
-   }
-   mesh->FinalizeTopology();
-
-   return mesh;
-}
-
 } // namespace build_dof_to_arrays
 
 } // namespace mfem

tests/unit/fem/test_domain_int.cpp

@@ -11,8 +11,10 @@
 
 #include "mfem.hpp"
 #include "unit_tests.hpp"
+#include "common_get_mesh.hpp"
 
 using namespace mfem;
+using namespace mfem_test_fem;
 
 namespace domain_int
 {
@@ -46,29 +48,6 @@ enum FEType
    L2I_FEC,
 };
 
-enum MeshType
-{
-   SEGMENT = 0,
-   QUADRILATERAL = 1,
-   TRIANGLE2A = 2,
-   TRIANGLE2B = 3,
-   TRIANGLE2C = 4,
-   TRIANGLE4 = 5,
-   MIXED2D = 6,
-   HEXAHEDRON = 7,
-   HEXAHEDRON2A = 8,
-   HEXAHEDRON2B = 9,
-   HEXAHEDRON2C = 10,
-   HEXAHEDRON2D = 11,
-   WEDGE2 = 12,
-   TETRAHEDRA = 13,
-   WEDGE4 = 14,
-   MIXED3D6 = 15,
-   MIXED3D8 = 16
-};
-
-Mesh * GetMesh(MeshType type);
-
 TEST_CASE("Domain Integration (Scalar Field)",
           "[H1_FECollection]"
           "[L2_FECollection]"
@@ -80,7 +59,7 @@ TEST_CASE("Domain Integration (Scalar Field)",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int  dim = mesh->Dimension();
       mesh->UniformRefinement();
 
@@ -155,7 +134,7 @@ TEST_CASE("Domain Integration (Vector Field)",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int  dim = mesh->Dimension();
       int sdim = mesh->SpaceDimension();
       mesh->UniformRefinement();
@@ -258,7 +237,7 @@ TEST_CASE("Domain Integration in Parallel (Scalar Field)",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int dim = mesh->Dimension();
       while (mesh->GetNE() < num_procs)
       {
@@ -340,7 +319,7 @@ TEST_CASE("Domain Integration in Parallel (Vector Field)",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int  dim = mesh->Dimension();
       int sdim = mesh->SpaceDimension();
       while (mesh->GetNE() < num_procs)
@@ -427,237 +406,4 @@ TEST_CASE("Domain Integration in Parallel (Vector Field)",
 
 #endif // MFEM_USE_MPI
 
-Mesh * GetMesh(MeshType type)
-{
-   Mesh * mesh = NULL;
-
-   switch (type)
-   {
-      case SEGMENT:
-         mesh = new Mesh(1, 2, 1);
-         mesh->AddVertex(0.0);
-         mesh->AddVertex(a_);
-
-         mesh->AddSegment(0, 1);
-
-         mesh->AddBdrPoint(0);
-         mesh->AddBdrPoint(1);
-         break;
-      case QUADRILATERAL:
-         mesh = new Mesh(2, 4, 1);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddQuad(0, 1, 2, 3);
-         break;
-      case TRIANGLE2A:
-         mesh = new Mesh(2, 4, 2);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddTriangle(0, 1, 2);
-         mesh->AddTriangle(2, 3, 0);
-         break;
-      case TRIANGLE2B:
-         mesh = new Mesh(2, 4, 2);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddTriangle(1, 2, 0);
-         mesh->AddTriangle(3, 0, 2);
-         break;
-      case TRIANGLE2C:
-         mesh = new Mesh(2, 4, 2);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-
-         mesh->AddTriangle(2, 0, 1);
-         mesh->AddTriangle(0, 2, 3);
-         break;
-      case TRIANGLE4:
-         mesh = new Mesh(2, 5, 4);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-         mesh->AddVertex(0.5 * a_, 0.5 * b_);
-
-         mesh->AddTriangle(0, 1, 4);
-         mesh->AddTriangle(1, 2, 4);
-         mesh->AddTriangle(2, 3, 4);
-         mesh->AddTriangle(3, 0, 4);
-         break;
-      case MIXED2D:
-         mesh = new Mesh(2, 6, 4);
-         mesh->AddVertex(0.0, 0.0);
-         mesh->AddVertex(a_, 0.0);
-         mesh->AddVertex(a_, b_);
-         mesh->AddVertex(0.0, b_);
-         mesh->AddVertex(0.5 * b_, 0.5 * b_);
-         mesh->AddVertex(a_ - 0.5 * b_, 0.5 * b_);
-
-         mesh->AddQuad(0, 1, 5, 4);
-         mesh->AddTriangle(1, 2, 5);
-         mesh->AddQuad(2, 3, 4, 5);
-         mesh->AddTriangle(3, 0, 4);
-         break;
-      case HEXAHEDRON:
-         mesh = new Mesh(3, 8, 1);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
-         break;
-      case HEXAHEDRON2A:
-      case HEXAHEDRON2B:
-      case HEXAHEDRON2C:
-      case HEXAHEDRON2D:
-         mesh = new Mesh(3, 12, 2);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(0.5 * a_, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.5 * a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(0.5 * a_, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.5 * a_, b_, c_);
-         mesh->AddVertex(0.0,b_, c_);
-
-         mesh->AddHex(0, 5, 11, 6, 1, 4, 10, 7);
-
-         switch (type)
-         {
-            case HEXAHEDRON2A: // Face Orientation 1
-               mesh->AddHex(4, 10, 7, 1, 3, 9, 8, 2);
-               break;
-            case HEXAHEDRON2B: // Face Orientation 3
-               mesh->AddHex(10, 7, 1, 4, 9, 8, 2, 3);
-               break;
-            case HEXAHEDRON2C: // Face Orientation 5
-               mesh->AddHex(7, 1, 4, 10, 8, 2, 3, 9);
-               break;
-            case HEXAHEDRON2D: // Face Orientation 7
-               mesh->AddHex(1, 4, 10, 7, 2, 3, 9, 8);
-               break;
-            default:
-               // Cannot happen
-               break;
-         }
-         break;
-      case WEDGE2:
-         mesh = new Mesh(3, 8, 2);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddWedge(0, 1, 2, 4, 5, 6);
-         mesh->AddWedge(0, 2, 3, 4, 6, 7);
-         break;
-      case TETRAHEDRA:
-         mesh = new Mesh(3, 8, 5);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddTet(0, 2, 7, 5);
-         mesh->AddTet(6, 7, 2, 5);
-         mesh->AddTet(4, 7, 5, 0);
-         mesh->AddTet(1, 0, 5, 2);
-         mesh->AddTet(3, 7, 0, 2);
-         break;
-      case WEDGE4:
-         mesh = new Mesh(3, 10, 4);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.5 * a_, 0.5 * b_, 0.0);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-         mesh->AddVertex(0.5 * a_, 0.5 * b_, c_);
-
-         mesh->AddWedge(0, 1, 4, 5, 6, 9);
-         mesh->AddWedge(1, 2, 4, 6, 7, 9);
-         mesh->AddWedge(2, 3, 4, 7, 8, 9);
-         mesh->AddWedge(3, 0, 4, 8, 5, 9);
-         break;
-      case MIXED3D6:
-         mesh = new Mesh(3, 12, 6);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.5 * c_, 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(a_ - 0.5 * c_, 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(a_ - 0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
-         mesh->AddWedge(0, 4, 8, 1, 5, 9);
-         mesh->AddWedge(1, 5, 9, 2, 6, 10);
-         mesh->AddWedge(2, 6, 10, 3, 7, 11);
-         mesh->AddWedge(3, 7, 11, 0, 4, 8);
-         mesh->AddHex(4, 5, 6, 7, 8, 9, 10, 11);
-         break;
-      case MIXED3D8:
-         mesh = new Mesh(3, 10, 8);
-         mesh->AddVertex(0.0, 0.0, 0.0);
-         mesh->AddVertex(a_, 0.0, 0.0);
-         mesh->AddVertex(a_, b_, 0.0);
-         mesh->AddVertex(0.0, b_, 0.0);
-         mesh->AddVertex(0.25 * a_, 0.5 * b_, 0.5 * c_);
-         mesh->AddVertex(0.75 * a_, 0.5 * b_, 0.5 * c_);
-         mesh->AddVertex(0.0, 0.0, c_);
-         mesh->AddVertex(a_, 0.0, c_);
-         mesh->AddVertex(a_, b_, c_);
-         mesh->AddVertex(0.0, b_, c_);
-
-         mesh->AddWedge(0, 3, 4, 1, 2, 5);
-         mesh->AddWedge(3, 9, 4, 2, 8, 5);
-         mesh->AddWedge(9, 6, 4, 8, 7, 5);
-         mesh->AddWedge(6, 0, 4, 7, 1, 5);
-         mesh->AddTet(0, 3, 9, 4);
-         mesh->AddTet(0, 9, 6, 4);
-         mesh->AddTet(1, 7, 2, 5);
-         mesh->AddTet(8, 2, 7, 5);
-         break;
-   }
-   mesh->FinalizeTopology();
-
-   return mesh;
-}
-
 } // namespace domain_int

tests/unit/fem/test_eigs.cpp

@@ -11,38 +11,21 @@
 
 #include "mfem.hpp"
 #include "unit_tests.hpp"
+#include "common_get_mesh.hpp"
 
 using namespace mfem;
+using namespace mfem_test_fem;
 
 namespace eigs
 {
 
+#if defined MFEM_USE_LAPACK || defined MFEM_USE_MPI
+
 static double a_ = M_PI;
 static double b_ = M_PI / sqrt(2.0);
 static double c_ = M_PI / 2.0;
 
-enum MeshType
-{
-   SEGMENT = 0,
-   QUADRILATERAL = 1,
-   TRIANGLE2A = 2,
-   TRIANGLE2B = 3,
-   TRIANGLE2C = 4,
-   TRIANGLE4 = 5,
-   MIXED2D = 6,
-   HEXAHEDRON = 7,
-   HEXAHEDRON2A = 8,
-   HEXAHEDRON2B = 9,
-   HEXAHEDRON2C = 10,
-   HEXAHEDRON2D = 11,
-   WEDGE2 = 12,
-   TETRAHEDRA = 13,
-   WEDGE4 = 14,
-   MIXED3D6 = 15,
-   MIXED3D8 = 16
-};
-
-Mesh * GetMesh(MeshType type);
+#endif
 
 int eigs[21] =
 {
@@ -63,7 +46,7 @@ TEST_CASE("Laplacian Eigenvalues",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int  dim = mesh->Dimension();
       if (dim < 3 ||
           mt == MeshType::HEXAHEDRON ||
@@ -178,7 +161,7 @@ TEST_CASE("Laplacian Eigenvalues in Parallel",
    for (int mt = (int)MeshType::SEGMENT;
         mt <= (int)MeshType::MIXED3D8; mt++)
    {
-      Mesh *mesh = GetMesh((MeshType)mt);
+      Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
       int  dim = mesh->Dimension();
       if (dim < 3 ||
           mt == MeshType::HEXAHEDRON ||
@@ -273,401 +256,4 @@ TEST_CASE("Laplacian Eigenvalues in Parallel",
 
 #endif // MFEM_USE_MPI
 
-Mesh * GetMesh(MeshType type)
-{
-   Mesh * mesh = NULL;
-   double c[3];
-   int    v[8];
-
-   switch (type)
-   {
-      case SEGMENT:
-         mesh = new Mesh(1, 2, 1);
-         c[0] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_;
-         mesh->AddVertex(c);
-         v[0] = 0; v[1] = 1;
-         mesh->AddSegment(v);
-         {
-            Element * el = mesh->NewElement(Geometry::POINT);
-            el->SetAttribute(1);
-            el->SetVertices(&v[0]);
-            mesh->AddBdrElement(el);
-         }
-         {
-            Element * el = mesh->NewElement(Geometry::POINT);
-            el->SetAttribute(2);
-            el->SetVertices(&v[1]);
-            mesh->AddBdrElement(el);
-         }
-         break;
-      case QUADRILATERAL:
-         mesh = new Mesh(2, 4, 1);
-         c[0] = 0.0; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
-         mesh->AddQuad(v);
-         break;
-      case TRIANGLE2A:
-         mesh = new Mesh(2, 4, 2);
-         c[0] = 0.0; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 2;
-         mesh->AddTri(v);
-         v[0] = 2; v[1] = 3; v[2] = 0;
-         mesh->AddTri(v);
-         break;
-      case TRIANGLE2B:
-         mesh = new Mesh(2, 4, 2);
-         c[0] = 0.0; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_;
-         mesh->AddVertex(c);
-
-         v[0] = 1; v[1] = 2; v[2] = 0;
-         mesh->AddTri(v);
-         v[0] = 3; v[1] = 0; v[2] = 2;
-         mesh->AddTri(v);
-         break;
-      case TRIANGLE2C:
-         mesh = new Mesh(2, 4, 2);
-         c[0] = 0.0; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_;
-         mesh->AddVertex(c);
-
-         v[0] = 2; v[1] = 0; v[2] = 1;
-         mesh->AddTri(v);
-         v[0] = 0; v[1] = 2; v[2] = 3;
-         mesh->AddTri(v);
-         break;
-      case TRIANGLE4:
-         mesh = new Mesh(2, 5, 4);
-         c[0] = 0.0; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * a_; c[1] = 0.5 * b_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 4;
-         mesh->AddTri(v);
-         v[0] = 1; v[1] = 2; v[2] = 4;
-         mesh->AddTri(v);
-         v[0] = 2; v[1] = 3; v[2] = 4;
-         mesh->AddTri(v);
-         v[0] = 3; v[1] = 0; v[2] = 4;
-         mesh->AddTri(v);
-         break;
-      case MIXED2D:
-         mesh = new Mesh(2, 6, 4);
-         c[0] = 0.0; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * b_; c[1] = 0.5 * b_;
-         mesh->AddVertex(c);
-         c[0] = a_ - 0.5 * b_; c[1] = 0.5 * b_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 5; v[3] = 4;
-         mesh->AddQuad(v);
-         v[0] = 1; v[1] = 2; v[2] = 5;
-         mesh->AddTri(v);
-         v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 5;
-         mesh->AddQuad(v);
-         v[0] = 3; v[1] = 0; v[2] = 4;
-         mesh->AddTri(v);
-         break;
-      case HEXAHEDRON:
-         mesh = new Mesh(3, 8, 1);
-         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
-         v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
-         mesh->AddHex(v);
-         break;
-      case HEXAHEDRON2A:
-      case HEXAHEDRON2B:
-      case HEXAHEDRON2C:
-      case HEXAHEDRON2D:
-         mesh = new Mesh(3, 12, 2);
-         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 5; v[2] = 11; v[3] = 6;
-         v[4] = 1; v[5] = 4; v[6] = 10; v[7] = 7;
-         mesh->AddHex(v);
-
-         switch (type)
-         {
-            case HEXAHEDRON2A: // Face Orientation 1
-               v[0] = 4; v[1] = 10; v[2] = 7; v[3] = 1;
-               v[4] = 3; v[5] = 9; v[6] = 8; v[7] = 2;
-               mesh->AddHex(v);
-               break;
-            case HEXAHEDRON2B: // Face Orientation 3
-               v[0] = 10; v[1] = 7; v[2] = 1; v[3] = 4;
-               v[4] = 9; v[5] = 8; v[6] = 2; v[7] = 3;
-               mesh->AddHex(v);
-               break;
-            case HEXAHEDRON2C: // Face Orientation 5
-               v[0] = 7; v[1] = 1; v[2] = 4; v[3] = 10;
-               v[4] = 8; v[5] = 2; v[6] = 3; v[7] = 9;
-               mesh->AddHex(v);
-               break;
-            case HEXAHEDRON2D: // Face Orientation 7
-               v[0] = 1; v[1] = 4; v[2] = 10; v[3] = 7;
-               v[4] = 2; v[5] = 3; v[6] = 9; v[7] = 8;
-               mesh->AddHex(v);
-               break;
-            default:
-               // Cannot happen
-               break;
-         }
-         break;
-      case WEDGE2:
-         mesh = new Mesh(3, 8, 2);
-         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 4; v[4] = 5; v[5] = 6;
-         mesh->AddWedge(v);
-         v[0] = 0; v[1] = 2; v[2] = 3; v[3] = 4; v[4] = 6; v[5] = 7;
-         mesh->AddWedge(v);
-         break;
-      case TETRAHEDRA:
-         mesh = new Mesh(3, 8, 5);
-         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 2; v[2] = 7; v[3] = 5;
-         mesh->AddTet(v);
-         v[0] = 6; v[1] = 7; v[2] = 2; v[3] = 5;
-         mesh->AddTet(v);
-         v[0] = 4; v[1] = 7; v[2] = 5; v[3] = 0;
-         mesh->AddTet(v);
-         v[0] = 1; v[1] = 0; v[2] = 5; v[3] = 2;
-         mesh->AddTet(v);
-         v[0] = 3; v[1] = 7; v[2] = 0; v[3] = 2;
-         mesh->AddTet(v);
-         break;
-      case WEDGE4:
-         mesh = new Mesh(3, 10, 4);
-         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * a_; c[1] = 0.5 * b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * a_; c[1] = 0.5 * b_; c[2] = c_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 4; v[3] = 5; v[4] = 6; v[5] = 9;
-         mesh->AddWedge(v);
-         v[0] = 1; v[1] = 2; v[2] = 4; v[3] = 6; v[4] = 7; v[5] = 9;
-         mesh->AddWedge(v);
-         v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 7; v[4] = 8; v[5] = 9;
-         mesh->AddWedge(v);
-         v[0] = 3; v[1] = 0; v[2] = 4; v[3] = 8; v[4] = 5; v[5] = 9;
-         mesh->AddWedge(v);
-         break;
-      case MIXED3D6:
-         mesh = new Mesh(3, 12, 6);
-         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * c_; c[1] = 0.5 * c_; c[2] = 0.5 * c_;
-         mesh->AddVertex(c);
-         c[0] = a_ - 0.5 * c_; c[1] = 0.5 * c_; c[2] = 0.5 * c_;
-         mesh->AddVertex(c);
-         c[0] = a_ - 0.5 * c_; c[1] = b_ - 0.5 * c_; c[2] = 0.5 * c_;
-         mesh->AddVertex(c);
-         c[0] = 0.5 * c_; c[1] = b_ - 0.5 * c_; c[2] = 0.5 * c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
-         v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
-         mesh->AddHex(v);
-         v[0] = 0; v[1] = 4; v[2] = 8; v[3] = 1; v[4] = 5; v[5] = 9;
-         mesh->AddWedge(v);
-         v[0] = 1; v[1] = 5; v[2] = 9; v[3] = 2; v[4] = 6; v[5] = 10;
-         mesh->AddWedge(v);
-         v[0] = 2; v[1] = 6; v[2] = 10; v[3] = 3; v[4] = 7; v[5] = 11;
-         mesh->AddWedge(v);
-         v[0] = 3; v[1] = 7; v[2] = 11; v[3] = 0; v[4] = 4; v[5] = 8;
-         mesh->AddWedge(v);
-         v[0] = 4; v[1] = 5; v[2] = 6; v[3] = 7;
-         v[4] = 8; v[5] = 9; v[6] = 10; v[7] = 11;
-         mesh->AddHex(v);
-         break;
-      case MIXED3D8:
-         mesh = new Mesh(3, 10, 8);
-         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
-         mesh->AddVertex(c);
-
-         c[0] = 0.25 * a_; c[1] = 0.5 * b_; c[2] = 0.5 * c_;
-         mesh->AddVertex(c);
-         c[0] = 0.75 * a_; c[1] = 0.5 * b_; c[2] = 0.5 * c_;
-         mesh->AddVertex(c);
-
-         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = 0.0; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = a_; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-         c[0] = 0.0; c[1] = b_; c[2] = c_;
-         mesh->AddVertex(c);
-
-         v[0] = 0; v[1] = 3; v[2] = 4; v[3] = 1; v[4] = 2; v[5] = 5;
-         mesh->AddWedge(v);
-         v[0] = 3; v[1] = 9; v[2] = 4; v[3] = 2; v[4] = 8; v[5] = 5;
-         mesh->AddWedge(v);
-         v[0] = 9; v[1] = 6; v[2] = 4; v[3] = 8; v[4] = 7; v[5] = 5;
-         mesh->AddWedge(v);
-         v[0] = 6; v[1] = 0; v[2] = 4; v[3] = 7; v[4] = 1; v[5] = 5;
-         mesh->AddWedge(v);
-         v[0] = 0; v[1] = 3; v[2] = 9; v[3] = 4;
-         mesh->AddTet(v);
-         v[0] = 0; v[1] = 9; v[2] = 6; v[3] = 4;
-         mesh->AddTet(v);
-         v[0] = 1; v[1] = 7; v[2] = 2; v[3] = 5;
-         mesh->AddTet(v);
-         v[0] = 8; v[1] = 2; v[2] = 7; v[3] = 5;
-         mesh->AddTet(v);
-         break;
-   }
-   mesh->FinalizeTopology();
-
-   return mesh;
-}
-
 } // namespace eigs

tests/unit/fem/test_lor.cpp

@@ -0,0 +1,120 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+#include <fstream>
+#include <sstream>
+
+using namespace mfem;
+
+// Tests the use of refined/LOR grid function coefficients.
+//
+// Given a space fes, and a refined space fes_refined (either vector or scalar
+// spaces), projects coeff_1 onto a grid function in fes, and then creates the
+// corresponding grid function coefficient (could be scalar, vector, grad, div,
+// or curl grid function coefficients). Then, this grid function coefficient is
+// projected onto fes_refined, and compared with the result of projecting
+// coeff_2 onto fes_refined.
+//
+// If coeff_1 can be represented exactly in fes, then these two projections
+// should be identical.
+template <typename GridFunctionCoeffType=GridFunctionCoefficient,
+          typename CoeffType1, typename CoeffType2>
+void TestRefinedGridFunctionCoefficient(
+   FiniteElementSpace &fes, FiniteElementSpace &fes_refined,
+   CoeffType1 &coeff_1, CoeffType2 &coeff_2)
+{
+   GridFunction gf(&fes);
+   gf.ProjectCoefficient(coeff_1);
+   GridFunctionCoeffType gf_coeff(&gf);
+
+   GridFunction gf_refined_1(&fes_refined), gf_refined_2(&fes_refined);
+   gf_refined_1.ProjectCoefficient(coeff_2);
+   gf_refined_2.ProjectCoefficient(gf_coeff);
+
+   gf_refined_2 -= gf_refined_1;
+   REQUIRE(gf_refined_2.Normlinf() == MFEM_Approx(0.0));
+}
+
+// Forward declarations for functions defined in test_lin_interp.cpp
+namespace lin_interp
+{
+double f2(const Vector & x);
+void F2(const Vector & x, Vector & v);
+void Grad_f2(const Vector & x, Vector & df);
+double curlF2(const Vector & x);
+double DivF2(const Vector & x);
+double f3(const Vector & x);
+void F3(const Vector & x, Vector & v);
+void Grad_f3(const Vector & x, Vector & df);
+void CurlF3(const Vector & x, Vector & df);
+double DivF3(const Vector & x);
+}
+
+namespace detail
+{
+Mesh MakeCartesian(int dim, int nx)
+{
+   if (dim == 1) { return Mesh::MakeCartesian1D(nx); }
+   else if (dim == 2) { return Mesh::MakeCartesian2D(nx, nx, Element::QUADRILATERAL); }
+   else { return Mesh::MakeCartesian3D(nx, nx, nx, Element::HEXAHEDRON); }
+}
+}
+
+TEST_CASE("LOR GridFunction Coefficient", "[LOR][GridFunctionCoefficient]")
+{
+   auto dim = GENERATE(2, 3);
+   Mesh mesh = detail::MakeCartesian(dim, 2);
+   Mesh mesh_refined = Mesh::MakeRefined(mesh, 3, Quadrature1D::GaussLobatto);
+
+   int order = 1;
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+   FiniteElementSpace fes_refined(&mesh_refined, &fec);
+
+   FiniteElementSpace vec_fes(&mesh, &fec, dim);
+   FiniteElementSpace vec_fes_refined(&mesh_refined, &fec, dim);
+
+   auto f = (dim == 2)? lin_interp::f2 : lin_interp::f3;
+   auto F = (dim == 2)? lin_interp::F2 : lin_interp::F3;
+   auto grad = (dim == 2)? lin_interp::Grad_f2 : lin_interp::Grad_f3;
+   auto div = (dim == 2)? lin_interp::DivF2 : lin_interp::DivF3;
+
+   FunctionCoefficient f_coeff(f);
+   VectorFunctionCoefficient vec_coeff(dim, F);
+   VectorFunctionCoefficient grad_coeff(dim, grad);
+   FunctionCoefficient div_coeff(div);
+
+   TestRefinedGridFunctionCoefficient<GridFunctionCoefficient>(
+      fes, fes_refined, f_coeff, f_coeff);
+   TestRefinedGridFunctionCoefficient<VectorGridFunctionCoefficient>(
+      vec_fes, vec_fes_refined, vec_coeff, vec_coeff);
+   TestRefinedGridFunctionCoefficient<DivergenceGridFunctionCoefficient>(
+      vec_fes, fes_refined, vec_coeff, div_coeff);
+   TestRefinedGridFunctionCoefficient<GradientGridFunctionCoefficient>(
+      fes, vec_fes_refined, f_coeff, grad_coeff);
+
+   // Curl is treated differently for dim = 2 (where it is a scalar quantity)
+   // and dim = 3 (where it is a vector quantity)
+   if (dim == 2)
+   {
+      FunctionCoefficient curl_coeff(lin_interp::curlF2);
+      TestRefinedGridFunctionCoefficient<CurlGridFunctionCoefficient>(
+         vec_fes, fes_refined, vec_coeff, curl_coeff);
+   }
+   else if (dim == 3)
+   {
+      VectorFunctionCoefficient curl_coeff(dim, lin_interp::CurlF3);
+      TestRefinedGridFunctionCoefficient<CurlGridFunctionCoefficient>(
+         vec_fes, vec_fes_refined, vec_coeff, curl_coeff);
+   }
+}