small updates prior to passing the code along for comparison

Eikonal solver for MFEM

.gitignore

@@ -57,6 +57,8 @@ examples/ex2[0-9]
 examples/ex2[0-9]p
 examples/ex3[0-9]
 examples/ex3[0-9]p
+examples/ex4[0-9]
+examples/ex4[0-9]p
 
 examples/refined.mesh
 examples/displaced.mesh

.gitlab-ci.yml

@@ -13,6 +13,9 @@
 # at Lawrence Livermore National Laboratory (LLNL). This entire pipeline is
 # LLNL-specific!
 
+include:
+  - project: 'lc-templates/id_tokens'
+    file: 'id_tokens.yml'
 
 # The pipeline is divided into stages. Usually, jobs in a given stage wait for
 # the preceding stages to complete before to start. However, we sometimes use

.gitlab/configs/common.yml

@@ -9,6 +9,10 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
+include:
+  - project: 'lc-templates/id_tokens'
+    file: 'id_tokens.yml'
+
 # We define the following GitLab pipeline variables:
 variables:
 

.gitlab/configs/lassen-config.yml

@@ -35,9 +35,8 @@ variables:
     - when: on_success
 
 # Lassen uses a different job scheduler (spectrum lsf) that does not allow
-# pre-allocation the same way slurm does.  We use pdebug queue on lassen
-# to speed-up the allocation.  However this would not be scalable to
-# multiple builds.
+# pre-allocation the same way slurm does.  We use the pci queue on lassen
+# to speed-up the allocation.
 .build_and_test_on_lassen:
   extends: [.on_lassen]
   stage: build_and_test
@@ -45,5 +44,5 @@ variables:
     - echo ${MFEM_DATA_DIR}
     - echo ${SPEC}
     # Next script uses 'THREADS': leaving it empty --> it uses 'make all -j'
-    - lalloc 1 -W 45 -q pdebug --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - lalloc 1 -W 45 -q pci --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
   needs: [setup]

.gitlab/configs/quartz-config.yml

@@ -52,4 +52,4 @@ variables:
     - echo ${JOBID}
     - echo ${MFEM_DATA_DIR}
     - echo ${SPEC}
-    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) --reservation=ci -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data

.gitlab/scripts/baseline

@@ -32,11 +32,11 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
 
 # run
 if [[ "${MACHINE_NAME}" == "quartz" || "${MACHINE_NAME}" == "ruby" ]]; then
-  salloc --nodes=1 -p pdebug ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  salloc --nodes=1 --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
   salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "lassen" ]]; then
-  lalloc 1 -q pdebug ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  lalloc 1 -q pci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
   echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
   exit 1

CHANGELOG

@@ -11,6 +11,9 @@
 Version 4.7.1 (development)
 ===========================
 
+- Added an MFEM example for the eikonal equation. This new solver is based on
+  the proximal Galerkin method introduced by Keith and Surowiec.
+
 
 Version 4.7, released on May 7, 2024
 ====================================

doc/CodeDocumentation.dox

@@ -114,7 +114,9 @@ namespace mfem {
  * - <a class="el" href="ex37p_8cpp_source.html">Example 37p</a>: parallel topology optimization
  * - <a class="el" href="ex38_8cpp_source.html">Example 38</a>: cut-surface and cut-volume integration
  * - <a class="el" href="ex39_8cpp_source.html">Example 39</a>: named mesh attributes
- * - <a class="el" href="ex39p_8cpp_source.html">Example 39</a>: parallel named mesh attributes
+ * - <a class="el" href="ex39p_8cpp_source.html">Example 39p</a>: parallel named mesh attributes
+ * - <a class="el" href="ex40_8cpp_source.html">Example 40</a>: eikonal equation
+ * - <a class="el" href="ex40p_8cpp_source.html">Example 40p</a>: parallel eikonal equation
  *
  * <H4>AmgX Examples</H4>
  * - Variants of Examples

examples/CMakeLists.txt

@@ -45,6 +45,7 @@ list(APPEND ALL_EXE_SRCS
   ex37.cpp
   ex38.cpp
   ex39.cpp
+  ex40.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -87,6 +88,7 @@ if (MFEM_USE_MPI)
     ex36p.cpp
     ex37p.cpp
     ex39p.cpp
+    ex40p.cpp
   )
 endif()
 
@@ -146,10 +148,10 @@ if (MFEM_ENABLE_TESTING)
   # Add CUDA/HIP tests.
   set(DEVICE_EXAMPLES
       # serial examples with device support:
-      ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26 ex34
+      ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
       # parallel examples with device support:
-      ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p ex24p ex25p ex26p
-      ex34p ex35p)
+      ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p ex22p ex24p ex25p
+      ex26p ex34p ex35p)
   set(MFEM_TEST_DEVICE)
   if (MFEM_USE_CUDA)
     set(MFEM_TEST_DEVICE "cuda")
@@ -159,6 +161,11 @@ if (MFEM_ENABLE_TESTING)
   if (MFEM_TEST_DEVICE)
     foreach(TEST_NAME ${DEVICE_EXAMPLES})
       set(THIS_TEST_OPTIONS "-no-vis" "-d" "${MFEM_TEST_DEVICE}")
+      if (${TEST_NAME} MATCHES "ex14p")
+        list(APPEND THIS_TEST_OPTIONS "-rs" "2" "-rp" "0" "-pa")
+      elseif (${TEST_NAME} MATCHES "ex14")
+        list(APPEND THIS_TEST_OPTIONS "-r" "2" "-pa")
+      endif()
       if (NOT (${TEST_NAME} MATCHES ".*p$"))
         add_test(NAME ${TEST_NAME}_${MFEM_TEST_DEVICE}_ser
           COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})

examples/ex40.cpp

@@ -0,0 +1,374 @@
+//                                MFEM Example 40
+//
+// Compile with: make ex40
+//
+// Sample runs: ex40 -step 10 -gr 2.0
+//              ex40 -step 10 -gr 2.0 -o 3 -r 1
+//              ex40 -step 10 -gr 2.0 -r 4 -m ../data/l-shape.mesh
+//              ex40 -step 10 -gr 2.0 -r 2 -m ../data/fichera.mesh
+//
+// Description: This example code demonstrates how to use MFEM to solve the
+//              eikonal equation,
+//
+//                      |∇𝑢| = 1 in Ω,  𝑢 = g on ∂Ω.
+//
+//              The solution of this problem coincides with the unique optimum of
+//              the nonlinear program
+//
+//                   maximize ∫_Ω 𝑢 d𝑥 subject to |∇𝑢| ≤ 1,   𝑢 = g on Ω,      (⋆)
+//
+//              which is the foundation for method implemented below.
+//
+//              Following the proximal Galerkin methodology [1] (see also Example
+//              36), we construct a Legendre function for the unit ball
+//              𝐵₁ := {𝑥 ∈ Rⁿ | |𝑥| < 1}. Our choice is the Hellinger entropy,
+//
+//                    h(𝑥) = −( 1 − |𝑥|² )^{1/2},
+//
+//              although other choices are possible, each leading to a slightly
+//              different algorithm. We then adaptively regularize the optimization
+//              problem (⋆) with the Bregman divergence of the Hellinger entropy,
+//
+//                 maximize  ∫_Ω 𝑢 d𝑥 - αₖ⁻¹ Dₕ(∇𝑢,∇𝑢ₖ₋₁)  subject to  𝑢 = g on Ω.
+//
+//              This results in a sequence of functions ( 𝜓ₖ , 𝑢ₖ ),
+//
+//                      𝑢ₖ → 𝑢,    𝜓ₖ/|𝜓ₖ| → ∇𝑢    as k → \infty,
+//
+//              defined by the nonlinear saddle-point problems
+//
+//               Find 𝜓ₖ ∈ H(div,Ω) and 𝑢ₖ ∈ L²(Ω) such that
+//               ( Zₖ(𝜓ₖ) , τ ) + ( 𝑢ₖ , ∇⋅τ ) = ⟨ g , τ⋅n ⟩         ∀ τ ∈ H(div,Ω)
+//               ( ∇⋅𝜓ₖ , v )                 = ( ∇⋅𝜓ₖ₋₁ - 1 , v )   ∀ v ∈ L²(Ω)
+//
+//              where Zₖ(𝜓) := ∇h⁻¹(αₖ 𝜓) = 𝜓 / ( αₖ⁻² + |𝜓|² )^{1/2} and step size
+//              αₖ > 0. These saddle-point problems are solved using a damped Newton's
+//              method. This example assumes that g = 0 and allows the step size to
+//              grow geometrically, αₖ = α₀rᵏ, where r ≥ 1 is the growth rate.
+//
+//              [1] Keith, B. and Surowiec, T. (2023) Proximal Galerkin: A structure-
+//                  preserving finite element method for pointwise bound constraints.
+//                  arXiv:2307.12444 [math.NA]
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   ZCoefficient(int vdim, GridFunction &psi_, real_t alpha_ = 1.0)
+      : VectorCoefficient(vdim), psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   DZCoefficient(int height, GridFunction &psi_, real_t alpha_ = 1.0)
+      : MatrixCoefficient(height),  psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_it = 5;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t growth_rate = 1.0;
+   real_t newton_scaling = 0.9;
+   real_t tichonov = 1e-1;
+   real_t tol = 1e-4;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Stopping criteria based on the difference between"
+                  "successive solution updates");
+   args.AddOption(&alpha, "-step", "--step",
+                  "Initial size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+                  "Growth rate of the step size alpha");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   MFEM_ASSERT(mesh.bdr_attributes.Size(),
+               "This example does not currently support meshes"
+               " without boundary attributes."
+              )
+
+   // 3. Postprocess the mesh.
+   // 3A. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 3B. Interpolate the geometry after refinement to control geometry error.
+   // NOTE: Minimum second-order interpolation is used to improve the accuracy.
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   RT_FECollection RTfec(order, dim);
+   FiniteElementSpace RTfes(&mesh, &RTfec);
+
+   L2_FECollection L2fec(order, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec);
+
+   cout << "Number of H(div) dofs: "
+        << RTfes.GetTrueVSize() << endl;
+   cout << "Number of L² dofs: "
+        << L2fes.GetTrueVSize() << endl;
+
+   // 5. Define the offsets for the block matrices
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = RTfes.GetVSize();
+   offsets[2] = L2fes.GetVSize();
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   // 6. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   GridFunction u_gf, delta_psi_gf;
+   delta_psi_gf.MakeRef(&RTfes,x,offsets[0]);
+   u_gf.MakeRef(&L2fes,x,offsets[1]);
+
+   GridFunction psi_old_gf(&RTfes);
+   GridFunction psi_gf(&RTfes);
+   GridFunction u_old_gf(&L2fes);
+
+   // 7. Define initial guesses for the solution variables.
+   delta_psi_gf = 0.0;
+   psi_gf = 0.0;
+   u_gf = 0.0;
+   psi_old_gf = psi_gf;
+   u_old_gf = u_gf;
+
+   // 8. Prepare for glvis output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   // 9. Coefficients to be used later.
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+   ConstantCoefficient tichonov_cf(tichonov);
+   ConstantCoefficient neg_tichonov_cf(-1.0*tichonov);
+   ZCoefficient Z(sdim, psi_gf, alpha);
+   DZCoefficient DZ(sdim, psi_gf, alpha);
+   ScalarVectorProductCoefficient neg_Z(-1.0, Z);
+   DivergenceGridFunctionCoefficient div_psi_cf(&psi_gf);
+   DivergenceGridFunctionCoefficient div_psi_old_cf(&psi_old_gf);
+   SumCoefficient psi_old_minus_psi(div_psi_old_cf, div_psi_cf, 1.0, -1.0);
+
+   // 10. Assemble constant matrices/vectors to avoid reassembly in the loop.
+   LinearForm b0, b1;
+   b0.MakeRef(&RTfes,rhs.GetBlock(0),0);
+   b1.MakeRef(&L2fes,rhs.GetBlock(1),0);
+
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(neg_Z));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(neg_one));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
+
+   BilinearForm a00(&RTfes);
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(DZ));
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(tichonov_cf));
+
+   MixedBilinearForm a10(&RTfes,&L2fes);
+   a10.AddDomainIntegrator(new VectorFEDivergenceIntegrator());
+   a10.Assemble();
+   a10.Finalize();
+   SparseMatrix &A10 = a10.SpMat();
+   SparseMatrix *A01 = Transpose(A10);
+
+   BilinearForm a11(&L2fes);
+   a11.AddDomainIntegrator(new MassIntegrator(neg_tichonov_cf));
+   a11.Assemble();
+   a11.Finalize();
+   SparseMatrix &A11 = a11.SpMat();
+
+   // 11. Iterate.
+   int k;
+   int total_iterations = 0;
+   real_t increment_u = 0.1;
+   GridFunction u_tmp(&L2fes);
+   for (k = 0; k < max_it; k++)
+   {
+      u_tmp = u_old_gf;
+      Z.SetAlpha(alpha);
+      DZ.SetAlpha(alpha);
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         b0.Assemble();
+         b1.Assemble();
+
+         a00.Assemble(false);
+         a00.Finalize(false);
+         SparseMatrix &A00 = a00.SpMat();
+
+         // Construct Schur-complement preconditioner
+         Vector A00_diag(a00.Height());
+         A00.GetDiag(A00_diag);
+         A00_diag.Reciprocal();
+         SparseMatrix *S = Mult_AtDA(*A01, A00_diag);
+
+         BlockDiagonalPreconditioner prec(offsets);
+         prec.SetDiagonalBlock(0,new DSmoother(A00));
+#ifndef MFEM_USE_SUITESPARSE
+         prec.SetDiagonalBlock(1,new GSSmoother(*S));
+#else
+         prec.SetDiagonalBlock(1,new UMFPackSolver(*S));
+#endif
+         prec.owns_blocks = 1;
+
+         BlockOperator A(offsets);
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,&A10);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,1,&A11);
+
+         GMRES(A,prec,rhs,x,0,2000,500,1e-12,0.0);
+         delete S;
+
+         u_tmp -= u_gf;
+         real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         // Damped Newton update
+         psi_gf.Add(newton_scaling, delta_psi_gf);
+         a00.Update();
+
+         if (visualization)
+         {
+            sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+         }
+
+         mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      alpha *= max(growth_rate, 1_r);
+
+   }
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+             << endl;
+
+   delete A01;
+   return 0;
+}
+
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(vdim);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   V = psi_vals;
+   V *= phi;
+}
+
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(height);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   K = 0.0;
+   for (int i = 0; i < height; i++)
+   {
+      K(i,i) = phi;
+      for (int j = 0; j < height; j++)
+      {
+         K(i,j) -= psi_vals(i) * psi_vals(j) * pow(phi, 3);
+      }
+   }
+}

examples/ex40p.cpp

@@ -0,0 +1,436 @@
+//                       MFEM Example 40 - Parallel Version
+//
+// Compile with: make ex40p
+//
+// Sample runs: mpirun -np 4 ex40p -step 10 -gr 2.0
+//              mpirun -np 4 ex40p -step 10 -gr 2.0 -o 3 -r 1
+//              mpirun -np 4 ex40p -step 10 -gr 2.0 -r 4 -m ../data/l-shape.mesh
+//              mpirun -np 4 ex40p -step 10 -gr 2.0 -r 2 -m ../data/fichera.mesh
+//
+// Description: This example code demonstrates how to use MFEM to solve the
+//              eikonal equation,
+//
+//                      |∇𝑢| = 1 in Ω,  𝑢 = g on ∂Ω.
+//
+//              The solution of this problem coincides with the unique optimum of
+//              the nonlinear program
+//
+//                   maximize ∫_Ω 𝑢 d𝑥 subject to |∇𝑢| ≤ 1,   𝑢 = g on Ω,      (⋆)
+//
+//              which is the foundation for method implemented below.
+//
+//              Following the proximal Galerkin methodology [1] (see also Example
+//              36), we construct a Legendre function for the unit ball
+//              𝐵₁ := {𝑥 ∈ Rⁿ | |𝑥| < 1}. Our choice is the Hellinger entropy,
+//
+//                    h(𝑥) = −( 1 − |𝑥|² )^{1/2},
+//
+//              although other choices are possible, each leading to a slightly
+//              different algorithm. We then adaptively regularize the optimization
+//              problem (⋆) with the Bregman divergence of the Hellinger entropy,
+//
+//                 maximize  ∫_Ω 𝑢 d𝑥 - αₖ⁻¹ Dₕ(∇𝑢,∇𝑢ₖ₋₁)  subject to  𝑢 = g on Ω.
+//
+//              This results in a sequence of functions ( 𝜓ₖ , 𝑢ₖ ),
+//
+//                      𝑢ₖ → 𝑢,    𝜓ₖ/|𝜓ₖ| → ∇𝑢    as k → \infty,
+//
+//              defined by the nonlinear saddle-point problems
+//
+//               Find 𝜓ₖ ∈ H(div,Ω) and 𝑢ₖ ∈ L²(Ω) such that
+//               ( Zₖ(𝜓ₖ) , τ ) + ( 𝑢ₖ , ∇⋅τ ) = ⟨ g , τ⋅n ⟩         ∀ τ ∈ H(div,Ω)
+//               ( ∇⋅𝜓ₖ , v )                 = ( ∇⋅𝜓ₖ₋₁ - 1 , v )   ∀ v ∈ L²(Ω)
+//
+//              where Zₖ(𝜓) := ∇h⁻¹(αₖ 𝜓) = 𝜓 / ( αₖ⁻² + |𝜓|² )^{1/2} and step size
+//              αₖ > 0. These saddle-point problems are solved using a damped Newton's
+//              method. This example assumes that g = 0 and allows the step size to
+//              grow geometrically, αₖ = α₀rᵏ, where r ≥ 1 is the growth rate.
+//
+//              [1] Keith, B. and Surowiec, T. (2023) Proximal Galerkin: A structure-
+//                  preserving finite element method for pointwise bound constraints.
+//                  arXiv:2307.12444 [math.NA]
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   ParGridFunction *psi;
+   real_t alpha;
+
+public:
+   ZCoefficient(int vdim, ParGridFunction &psi_, real_t alpha_ = 1.0)
+      : VectorCoefficient(vdim), psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   ParGridFunction *psi;
+   real_t alpha;
+
+public:
+   DZCoefficient(int height, ParGridFunction &psi_, real_t alpha_ = 1.0)
+      : MatrixCoefficient(height),  psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_it = 5;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t growth_rate = 1.0;
+   real_t newton_scaling = 0.9;
+   real_t tichonov = 1e-1;
+   real_t tol = 1e-4;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Stopping criteria based on the difference between"
+                  "successive solution updates");
+   args.AddOption(&alpha, "-step", "--step",
+                  "Initial size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+                  "Growth rate of the step size alpha");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. Read the mesh from the mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   MFEM_ASSERT(mesh.bdr_attributes.Size(),
+               "This example does not currently support meshes"
+               " without boundary attributes."
+              )
+
+   // 3. Postprocess the mesh.
+   // 3A. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 3B. Interpolate the geometry after refinement to control geometry error.
+   // NOTE: Minimum second-order interpolation is used to improve the accuracy.
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   RT_FECollection RTfec(order, dim);
+   ParFiniteElementSpace RTfes(&pmesh, &RTfec);
+
+   L2_FECollection L2fec(order, dim);
+   ParFiniteElementSpace L2fes(&pmesh, &L2fec);
+
+   int num_dofs_RT = RTfes.GlobalTrueVSize();
+   int num_dofs_L2 = L2fes.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of H(div) dofs: "
+           << num_dofs_RT << endl;
+      cout << "Number of L² dofs: "
+           << num_dofs_L2 << endl;
+   }
+
+   // 5. Define the offsets for the block matrices
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = RTfes.GetVSize();
+   offsets[2] = L2fes.GetVSize();
+   offsets.PartialSum();
+
+   Array<int> toffsets(3);
+   toffsets[0] = 0;
+   toffsets[1] = RTfes.GetTrueVSize();
+   toffsets[2] = L2fes.GetTrueVSize();
+   toffsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   BlockVector tx(toffsets), trhs(toffsets);
+   tx = 0.0; trhs = 0.0;
+
+   // 6. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   ParGridFunction u_gf, delta_psi_gf;
+   delta_psi_gf.MakeRef(&RTfes,x,offsets[0]);
+   u_gf.MakeRef(&L2fes,x,offsets[1]);
+
+   ParGridFunction psi_old_gf(&RTfes);
+   ParGridFunction psi_gf(&RTfes);
+   ParGridFunction u_old_gf(&L2fes);
+
+   // 7. Define initial guesses for the solution variables.
+   delta_psi_gf = 0.0;
+   psi_gf = 0.0;
+   u_gf = 0.0;
+   psi_old_gf = psi_gf;
+   u_old_gf = u_gf;
+
+   // 8. Prepare for glvis output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   // 9. Coefficients to be used later.
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+   ConstantCoefficient tichonov_cf(tichonov);
+   ConstantCoefficient neg_tichonov_cf(-1.0*tichonov);
+   ZCoefficient Z(sdim, psi_gf, alpha);
+   DZCoefficient DZ(sdim, psi_gf, alpha);
+   ScalarVectorProductCoefficient neg_Z(-1.0, Z);
+   DivergenceGridFunctionCoefficient div_psi_cf(&psi_gf);
+   DivergenceGridFunctionCoefficient div_psi_old_cf(&psi_old_gf);
+   SumCoefficient psi_old_minus_psi(div_psi_old_cf, div_psi_cf, 1.0, -1.0);
+
+   // 10. Assemble constant matrices/vectors to avoid reassembly in the loop.
+   ParLinearForm b0, b1;
+   b0.MakeRef(&RTfes,rhs.GetBlock(0),0);
+   b1.MakeRef(&L2fes,rhs.GetBlock(1),0);
+
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(neg_Z));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(neg_one));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
+
+   ParBilinearForm a00(&RTfes);
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(DZ));
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(tichonov_cf));
+
+   ParMixedBilinearForm a10(&RTfes,&L2fes);
+   a10.AddDomainIntegrator(new VectorFEDivergenceIntegrator());
+   a10.Assemble();
+   a10.Finalize();
+   HypreParMatrix *A10 = a10.ParallelAssemble();
+
+   HypreParMatrix *A01 = A10->Transpose();
+
+   ParBilinearForm a11(&L2fes);
+   a11.AddDomainIntegrator(new MassIntegrator(neg_tichonov_cf));
+   a11.Assemble();
+   a11.Finalize();
+   HypreParMatrix *A11 = a11.ParallelAssemble();
+
+   // 11. Iterate.
+   int k;
+   int total_iterations = 0;
+   real_t increment_u = 0.1;
+   ParGridFunction u_tmp(&L2fes);
+   for (k = 0; k < max_it; k++)
+   {
+      u_tmp = u_old_gf;
+      Z.SetAlpha(alpha);
+      DZ.SetAlpha(alpha);
+
+      if (myid == 0)
+      {
+         mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+      }
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         b0.Assemble();
+         b0.ParallelAssemble(trhs.GetBlock(0));
+
+         b1.Assemble();
+         b1.ParallelAssemble(trhs.GetBlock(1));
+
+         a00.Assemble(false);
+         a00.Finalize(false);
+         HypreParMatrix *A00 = a00.ParallelAssemble();
+
+         // Construct Schur-complement preconditioner
+         HypreParVector A00_diag(MPI_COMM_WORLD, A00->GetGlobalNumRows(),
+                                 A00->GetRowStarts());
+         A00->GetDiag(A00_diag);
+         HypreParMatrix S_tmp(*A01);
+         S_tmp.InvScaleRows(A00_diag);
+         HypreParMatrix *S = ParMult(A10, &S_tmp, true);
+
+         BlockDiagonalPreconditioner prec(toffsets);
+         HypreBoomerAMG P00(*A00);
+         P00.SetPrintLevel(0);
+         HypreBoomerAMG P11(*S);
+         P11.SetPrintLevel(0);
+         prec.SetDiagonalBlock(0,&P00);
+         prec.SetDiagonalBlock(1,&P11);
+
+         BlockOperator A(toffsets);
+         A.SetBlock(0,0,A00);
+         A.SetBlock(1,0,A10);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,1,A11);
+
+         GMRESSolver gmres(MPI_COMM_WORLD);
+         gmres.SetPrintLevel(-1);
+         gmres.SetRelTol(1e-8);
+         gmres.SetMaxIter(2000);
+         gmres.SetKDim(500);
+         gmres.SetOperator(A);
+         gmres.SetPreconditioner(prec);
+         gmres.Mult(trhs,tx);
+         delete S;
+         delete A00;
+
+         delta_psi_gf.SetFromTrueDofs(tx.GetBlock(0));
+         u_gf.SetFromTrueDofs(tx.GetBlock(1));
+
+         u_tmp -= u_gf;
+         real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         // Damped Newton update
+         psi_gf.Add(newton_scaling, delta_psi_gf);
+         a00.Update();
+
+         if (visualization)
+         {
+            sol_sock << "parallel " << num_procs << " " << myid << "\n";
+            sol_sock << "solution\n" << pmesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+         }
+
+         if (myid == 0)
+         {
+            mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+         }
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      if (myid == 0)
+      {
+         mfem::out << "Number of Newton iterations = " << j+1 << endl;
+         mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+      }
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      alpha *= max(growth_rate, 1_r);
+
+   }
+
+   // 12. Print stats.
+   if (myid == 0)
+   {
+      mfem::out << "\n Outer iterations: " << k+1
+                << "\n Total iterations: " << total_iterations
+                << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+                << endl;
+   }
+
+   // 13. Free the used memory.
+   delete A01;
+   delete A10;
+   delete A11;
+   return 0;
+}
+
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(vdim);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   V = psi_vals;
+   V *= phi;
+}
+
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(height);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   K = 0.0;
+   for (int i = 0; i < height; i++)
+   {
+      K(i,i) = phi;
+      for (int j = 0; j < height; j++)
+      {
+         K(i,j) -= psi_vals(i) * psi_vals(j) * pow(phi, 3);
+      }
+   }
+}

examples/makefile

@@ -23,14 +23,14 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
    ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-   ex31 ex33 ex34 ex36 ex37 ex38 ex39
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
    ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
    ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \
-   ex37p ex39p
-SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26 ex34
-PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p \
-   ex24p ex25p ex26p ex34p ex35p
+   ex37p ex39p ex40p
+SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
+PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
+   ex22p ex24p ex25p ex26p ex34p ex35p
 
 ifeq ($(MFEM_USE_LAPACK),YES)
    SEQ_EXAMPLES += ex38
@@ -138,6 +138,14 @@ ex10-test-seq: ex10
 	@$(call mfem-test,$<,, Serial example,-tf 5)
 ex10p-test-par: ex10p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-tf 5)
+ex14-test-seq-cuda: ex14
+	@$(call mfem-test,$<,, Serial CUDA example,-r 2 -pa -d cuda)
+ex14p-test-par-cuda: ex14p
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel CUDA example,-rs 2 -rp 0 -pa -d cuda)
+ex14-test-seq-hip: ex14
+	@$(call mfem-test,$<,, Serial HIP example,-r 2 -pa -d hip)
+ex14p-test-par-hip: ex14p
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel HIP example,-rs 2 -rp 0 -pa -d hip)
 ex15-test-seq: ex15
 	@$(call mfem-test,$<,, Serial example,-e 1)
 ex15p-test-par: ex15p

fem/bilininteg.hpp

@@ -1741,7 +1741,7 @@ public:
    { vector_fe.CalcPhysDShape(Trans, shape); }
 };
 
-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot  \nabla \cdot u, \nabla \cdot v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot \nabla u, \nabla \cdot v)$ in 2D
     or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ and $v$ is in $H(div)$. */
 class MixedGradDivIntegrator : public MixedScalarVectorIntegrator
 {
@@ -1780,7 +1780,7 @@ public:
    { scalar_fe.CalcPhysDivShape(Trans, shape); }
 };
 
-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla \cdot v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla v)$ in 2D
     or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H(div)$ and $v$ is in $H^1$. */
 class MixedDivGradIntegrator : public MixedScalarVectorIntegrator
 {
@@ -1820,7 +1820,7 @@ public:
    { scalar_fe.CalcPhysDivShape(Trans, shape); }
 };
 
-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla \cdot v)$ in 2D or 3D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla v)$ in 2D or 3D
     and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ or $L_2$ and $v$ is in $H^1$. */
 class MixedScalarWeakDivergenceIntegrator : public MixedScalarVectorIntegrator
 {

fem/coefficient.cpp

@@ -924,6 +924,75 @@ void MatrixArrayCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
    }
 }
 
+MatrixArrayVectorCoefficient::MatrixArrayVectorCoefficient (int dim)
+   : MatrixCoefficient (dim)
+{
+   Coeff.SetSize(height);
+   ownCoeff.SetSize(height);
+   for (int i = 0; i < height; i++)
+   {
+      Coeff[i] = NULL;
+      ownCoeff[i] = true;
+   }
+}
+
+void MatrixArrayVectorCoefficient::SetTime(real_t t)
+{
+   for (int i=0; i < height; i++)
+   {
+      if (Coeff[i]) { Coeff[i]->SetTime(t); }
+   }
+   this->MatrixCoefficient::SetTime(t);
+}
+
+void MatrixArrayVectorCoefficient::Set(int i, VectorCoefficient * c, bool own)
+{
+   MFEM_ASSERT(i < height && i >= 0, "Row "
+               << i << " does not exist. " <<
+               "Matrix height = " << height << ".");
+   if (ownCoeff[i]) { delete Coeff[i]; }
+   Coeff[i] = c;
+   ownCoeff[i] = own;
+}
+
+MatrixArrayVectorCoefficient::~MatrixArrayVectorCoefficient ()
+{
+   for (int i=0; i < height; i++)
+   {
+      if (ownCoeff[i]) { delete Coeff[i]; }
+   }
+}
+
+void MatrixArrayVectorCoefficient::Eval(int i, Vector &V,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(i < height && i >= 0, "Row "
+               << i << " does not exist. " <<
+               "Matrix height = " << height << ".");
+   if (Coeff[i])
+   {
+      Coeff[i] -> Eval(V, T, ip);
+   }
+   else
+   {
+      V = 0.0;
+   }
+}
+
+void MatrixArrayVectorCoefficient::Eval(DenseMatrix &K,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   K.SetSize(height, width);
+   Vector V(width);
+   for (int i = 0; i < height; i++)
+   {
+      this->Eval(i, V, T, ip);
+      K.SetRow(i, V);
+   }
+}
+
 void MatrixRestrictedCoefficient::SetTime(real_t t)
 {
    if (c) { c->SetTime(t); }
@@ -1041,6 +1110,27 @@ real_t DeterminantCoefficient::Eval(ElementTransformation &T,
    return ma.Det();
 }
 
+TraceCoefficient::TraceCoefficient(MatrixCoefficient &A)
+   : a(&A), ma(A.GetHeight(), A.GetWidth())
+{
+   MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
+               "TraceCoefficient:  "
+               "Argument must be a square matrix.");
+}
+
+void TraceCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+real_t TraceCoefficient::Eval(ElementTransformation &T,
+                              const IntegrationPoint &ip)
+{
+   a->Eval(ma, T, ip);
+   return ma.Trace();
+}
+
 VectorSumCoefficient::VectorSumCoefficient(int dim)
    : VectorCoefficient(dim),
      ACoef(NULL), BCoef(NULL),
@@ -1326,6 +1416,30 @@ void InverseMatrixCoefficient::Eval(DenseMatrix &M,
    M.Invert();
 }
 
+ExponentialMatrixCoefficient::ExponentialMatrixCoefficient(MatrixCoefficient &A)
+   : MatrixCoefficient(A.GetHeight(), A.GetWidth()), a(&A)
+{
+   MFEM_ASSERT(A.GetHeight() == A.GetWidth() && A.GetHeight() == 2,
+               "ExponentialMatrixCoefficient:  "
+               << "Argument must be a square 2x2 matrix."
+               << "  Height = " << A.GetHeight()
+               << ", Width = " << A.GetWidth());
+}
+
+void ExponentialMatrixCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
+void ExponentialMatrixCoefficient::Eval(DenseMatrix &M,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   a->Eval(M, T, ip);
+   M.Exponential();
+}
+
 OuterProductCoefficient::OuterProductCoefficient(VectorCoefficient &A,
                                                  VectorCoefficient &B)
    : MatrixCoefficient(A.GetVDim(), B.GetVDim()), a(&A), b(&B),

fem/coefficient.hpp

@@ -1334,6 +1334,46 @@ public:
    virtual ~MatrixArrayCoefficient();
 };
 
+/** @brief Matrix coefficient defined row-wise by an array of vector
+    coefficients. Rows that are not set will evaluate to zero. The
+    matrix coefficient is stored as an array indexing the rows of
+    the matrix. */
+class MatrixArrayVectorCoefficient : public MatrixCoefficient
+{
+private:
+   Array<VectorCoefficient *> Coeff;
+   Array<bool> ownCoeff;
+
+public:
+   /** @brief Construct a coefficient matrix of dimensions @a dim * @a dim. The
+       actual coefficients still need to be added with Set(). */
+   explicit MatrixArrayVectorCoefficient (int dim);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t) override;
+
+   /// Get the vector coefficient located at the i-th row of the matrix
+   VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }
+
+   /** @brief Set the coefficient located at the i-th row of the matrix.
+       By this will take ownership of the Coefficient passed in, but this
+       can be overridden with the @a own parameter. */
+   void Set(int i, VectorCoefficient * c, bool own=true);
+
+   using MatrixCoefficient::Eval;
+
+   /// Evaluate coefficient located at the i-th row of the matrix using integration
+   /// point @a ip.
+   void Eval(int i, Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+
+   /// Evaluate the matrix coefficient @a ip.
+   void Eval(DenseMatrix &K, ElementTransformation &T,
+             const IntegrationPoint &ip) override;
+
+   virtual ~MatrixArrayVectorCoefficient();
+};
+
 
 /** @brief Derived matrix coefficient that has the value of the parent matrix
     coefficient where it is active and is zero otherwise. */
@@ -1767,6 +1807,31 @@ public:
                        const IntegrationPoint &ip);
 };
 
+/// Scalar coefficient defined as the trace of a matrix coefficient
+class TraceCoefficient : public Coefficient
+{
+private:
+   MatrixCoefficient * a;
+
+   mutable DenseMatrix ma;
+
+public:
+   /// Construct with the matrix.
+   TraceCoefficient(MatrixCoefficient &A);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Evaluate the trace coefficient at @a ip.
+   virtual real_t Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+};
+
 /// Vector coefficient defined as the linear combination of two vectors
 class VectorSumCoefficient : public VectorCoefficient
 {
@@ -2118,7 +2183,7 @@ public:
                      const IntegrationPoint &ip);
 };
 
-/// Matrix coefficient defined as the transpose a matrix coefficient
+/// Matrix coefficient defined as the transpose of a matrix coefficient
 class TransposeMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2141,7 +2206,7 @@ public:
                      const IntegrationPoint &ip);
 };
 
-/// Matrix coefficient defined as the inverse a matrix coefficient.
+/// Matrix coefficient defined as the inverse of a matrix coefficient.
 class InverseMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2164,6 +2229,29 @@ public:
                      const IntegrationPoint &ip);
 };
 
+/// Matrix coefficient defined as the exponential of a matrix coefficient.
+class ExponentialMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   MatrixCoefficient * a;
+
+public:
+   /// Construct the matrix coefficient.  Result is $ \exp(A) $.
+   ExponentialMatrixCoefficient(MatrixCoefficient &A);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Evaluate the matrix coefficient at @a ip.
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
 /// Matrix coefficient defined as the outer product of two vector coefficients.
 class OuterProductCoefficient : public MatrixCoefficient
 {

fem/gridfunc.cpp

@@ -1321,9 +1321,9 @@ void GridFunction::ProjectVectorFieldOn(GridFunction &vec_field, int comp)
    }
 }
 
-void GridFunction::AccumulateAndCountDerivativeValues(int comp, int der_comp,
-                                                      GridFunction &der,
-                                                      Array<int> &zones_per_dof)
+void GridFunction::AccumulateAndCountDerivativeValues(
+   int comp, int der_comp, GridFunction &der,
+   Array<int> &zones_per_dof) const
 {
    FiniteElementSpace * der_fes = der.FESpace();
    ElementTransformation * transf;
@@ -1374,7 +1374,8 @@ void GridFunction::AccumulateAndCountDerivativeValues(int comp, int der_comp,
    }
 }
 
-void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der)
+void GridFunction::GetDerivative(int comp, int der_comp,
+                                 GridFunction &der) const
 {
    Array<int> overlap;
    AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);
@@ -2061,41 +2062,37 @@ void GridFunction::AccumulateAndCountBdrValues(
    Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
    Array<int> &values_counter)
 {
-   int i, j, fdof, d, ind, vdim;
-   real_t val;
-   const FiniteElement *fe;
-   ElementTransformation *transf;
    Array<int> vdofs;
    Vector vc;
 
    values_counter.SetSize(Size());
    values_counter = 0;
 
-   vdim = fes->GetVDim();
-
+   const int vdim = fes->GetVDim();
    HostReadWrite();
 
-   for (i = 0; i < fes->GetNBE(); i++)
+   for (int i = 0; i < fes->GetNBE(); i++)
    {
       if (attr[fes->GetBdrAttribute(i) - 1] == 0) { continue; }
 
-      fe = fes->GetBE(i);
-      fdof = fe->GetDof();
-      transf = fes->GetBdrElementTransformation(i);
+      const FiniteElement *fe = fes->GetBE(i);
+      const int fdof = fe->GetDof();
+      ElementTransformation *transf = fes->GetBdrElementTransformation(i);
       const IntegrationRule &ir = fe->GetNodes();
       fes->GetBdrElementVDofs(i, vdofs);
 
-      for (j = 0; j < fdof; j++)
+      for (int j = 0; j < fdof; j++)
       {
          const IntegrationPoint &ip = ir.IntPoint(j);
          transf->SetIntPoint(&ip);
          if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
-         for (d = 0; d < vdim; d++)
+         for (int d = 0; d < vdim; d++)
          {
             if (!vcoeff && !coeff[d]) { continue; }
 
-            val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
-            if ( (ind = vdofs[fdof*d+j]) < 0 )
+            real_t val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
+            int ind = vdofs[fdof*d+j];
+            if ( ind < 0 )
             {
                val = -val, ind = -1-ind;
             }
@@ -2117,10 +2114,11 @@ void GridFunction::AccumulateAndCountBdrValues(
    // iff A_ij != 0. It is sufficient to resolve just the first level of
    // dependency, since A is a projection matrix: A^n = A due to cR.cP = I.
    // Cases like these arise in 3D when boundary edges are constrained by
-   // (depend on) internal faces/elements. We use the virtual method
-   // GetBoundaryClosure from NCMesh to resolve the dependencies.
-
-   if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
+   // (depend on) internal faces/elements, or for internal boundaries in 2 or
+   // 3D. We use the virtual method GetBoundaryClosure from NCMesh to resolve
+   // the dependencies.
+   if (fes->Nonconforming() && (fes->GetMesh()->Dimension() == 2 ||
+                                fes->GetMesh()->Dimension() == 3))
    {
       Vector vals;
       Mesh *mesh = fes->GetMesh();
@@ -2128,26 +2126,19 @@ void GridFunction::AccumulateAndCountBdrValues(
       Array<int> bdr_edges, bdr_vertices, bdr_faces;
       ncmesh->GetBoundaryClosure(attr, bdr_vertices, bdr_edges, bdr_faces);
 
-      for (i = 0; i < bdr_edges.Size(); i++)
+      auto mark_dofs = [&](ElementTransformation &transf, const FiniteElement &fe)
       {
-         int edge = bdr_edges[i];
-         fes->GetEdgeVDofs(edge, vdofs);
-         if (vdofs.Size() == 0) { continue; }
-
-         transf = mesh->GetEdgeTransformation(edge);
-         transf->Attribute = -1; // TODO: set the boundary attribute
-         fe = fes->GetEdgeElement(edge);
          if (!vcoeff)
          {
-            vals.SetSize(fe->GetDof());
-            for (d = 0; d < vdim; d++)
+            vals.SetSize(fe.GetDof());
+            for (int d = 0; d < vdim; d++)
             {
                if (!coeff[d]) { continue; }
 
-               fe->Project(*coeff[d], *transf, vals);
+               fe.Project(*coeff[d], transf, vals);
                for (int k = 0; k < vals.Size(); k++)
                {
-                  ind = vdofs[d*vals.Size()+k];
+                  const int ind = vdofs[d*vals.Size()+k];
                   if (++values_counter[ind] == 1)
                   {
                      (*this)(ind) = vals(k);
@@ -2161,11 +2152,11 @@ void GridFunction::AccumulateAndCountBdrValues(
          }
          else // vcoeff != NULL
          {
-            vals.SetSize(vdim*fe->GetDof());
-            fe->Project(*vcoeff, *transf, vals);
+            vals.SetSize(vdim*fe.GetDof());
+            fe.Project(*vcoeff, transf, vals);
             for (int k = 0; k < vals.Size(); k++)
             {
-               ind = vdofs[k];
+               const int ind = vdofs[k];
                if (++values_counter[ind] == 1)
                {
                   (*this)(ind) = vals(k);
@@ -2176,6 +2167,26 @@ void GridFunction::AccumulateAndCountBdrValues(
                }
             }
          }
+      };
+
+      for (auto edge : bdr_edges)
+      {
+         fes->GetEdgeVDofs(edge, vdofs);
+         if (vdofs.Size() == 0) { continue; }
+
+         ElementTransformation *transf = mesh->GetEdgeTransformation(edge);
+         const FiniteElement *fe = fes->GetEdgeElement(edge);
+         mark_dofs(*transf, *fe);
+      }
+
+      for (auto face : bdr_faces)
+      {
+         fes->GetFaceVDofs(face, vdofs);
+         if (vdofs.Size() == 0) { continue; }
+
+         ElementTransformation *transf = mesh->GetFaceTransformation(face);
+         const FiniteElement *fe = fes->GetFaceElement(face);
+         mark_dofs(*transf, *fe);
       }
    }
 }
@@ -2228,26 +2239,37 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
       accumulate_dofs(dofs, lvec, *this, values_counter);
    }
 
-   if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
+   if (fes->Nonconforming() && (fes->GetMesh()->Dimension() == 2 ||
+                                fes->GetMesh()->Dimension() == 3))
    {
       Mesh *mesh = fes->GetMesh();
       NCMesh *ncmesh = mesh->ncmesh;
       Array<int> bdr_edges, bdr_vertices, bdr_faces;
       ncmesh->GetBoundaryClosure(bdr_attr, bdr_vertices, bdr_edges, bdr_faces);
 
-      for (int i = 0; i < bdr_edges.Size(); i++)
+      for (auto edge : bdr_edges)
       {
-         int edge = bdr_edges[i];
          fes->GetEdgeDofs(edge, dofs);
          if (dofs.Size() == 0) { continue; }
 
          T = mesh->GetEdgeTransformation(edge);
-         T->Attribute = -1; // TODO: set the boundary attribute
          fe = fes->GetEdgeElement(edge);
          lvec.SetSize(fe->GetDof());
          fe->Project(vcoeff, *T, lvec);
          accumulate_dofs(dofs, lvec, *this, values_counter);
       }
+
+      for (auto face : bdr_faces)
+      {
+         fes->GetFaceDofs(face, dofs);
+         if (dofs.Size() == 0) { continue; }
+
+         T = mesh->GetFaceTransformation(face);
+         fe = fes->GetFaceElement(face);
+         lvec.SetSize(fe->GetDof());
+         fe->Project(vcoeff, *T, lvec);
+         accumulate_dofs(dofs, lvec, *this, values_counter);
+      }
    }
 }
 

fem/gridfunc.hpp

@@ -321,7 +321,7 @@ public:
        @param[out] der       The resulting derivative (scalar function). The
                              FiniteElementSpace of this function must be set
                              before the call. */
-   void GetDerivative(int comp, int der_comp, GridFunction &der);
+   void GetDerivative(int comp, int der_comp, GridFunction &der) const;
 
    real_t GetDivergence(ElementTransformation &tr) const;
 
@@ -443,7 +443,7 @@ protected:
        GetDerivative() method; see its documentation. */
    void AccumulateAndCountDerivativeValues(int comp, int der_comp,
                                            GridFunction &der,
-                                           Array<int> &zones_per_dof);
+                                           Array<int> &zones_per_dof) const;
 
    void AccumulateAndCountBdrValues(Coefficient *coeff[],
                                     VectorCoefficient *vcoeff,

fem/gslib.cpp

@@ -1393,8 +1393,10 @@ GSOPGSLIB::~GSOPGSLIB()
 void GSOPGSLIB::UpdateIdentifiers(const Array<long long> &ids)
 {
    long long minval = ids.Min();
+#ifdef MFEM_USE_MPI
    MPI_Allreduce(MPI_IN_PLACE, &minval, 1, MPI_LONG_LONG_INT,
                  MPI_MIN, gsl_comm->c);
+#endif
    MFEM_VERIFY(minval >= 0, "Unique identifier cannot be negative.");
    if (gsl_data != NULL) { gslib_gs_free(gsl_data); }
    num_ids = ids.Size();

fem/pgridfunc.cpp

@@ -520,7 +520,7 @@ void ParGridFunction::CountElementsPerVDof(Array<int> &elem_per_vdof) const
 }
 
 void ParGridFunction::GetDerivative(int comp, int der_comp,
-                                    ParGridFunction &der)
+                                    ParGridFunction &der) const
 {
    Array<int> overlap;
    AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);

fem/pgridfunc.hpp

@@ -231,7 +231,7 @@ public:
    void CountElementsPerVDof(Array<int> &elem_per_vdof) const override;
 
    /// Parallel version of GridFunction::GetDerivative(); see its documentation.
-   void GetDerivative(int comp, int der_comp, ParGridFunction &der);
+   void GetDerivative(int comp, int der_comp, ParGridFunction &der) const;
 
    /** Sets the output vector @a dof_vals to the values of the degrees of
        freedom of element @a el. If @a el is greater than or equal to the number

fem/tmop.cpp

@@ -2949,6 +2949,15 @@ void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &s0,
    MFEM_VERIFY(surf_fit_pos == NULL,
                "Using both fitting approaches is not supported.");
 
+   const int dim = s0.FESpace()->GetMesh()->Dimension();
+   Mesh *mesh = s0.FESpace()->GetMesh();
+   MFEM_VERIFY(mesh->GetNodes()->Size() == dim*s0.Size(),
+               "Mesh and level-set polynomial order must be the same.");
+   const H1_FECollection *fec = dynamic_cast<const H1_FECollection *>
+                                (s0.FESpace()->FEColl());
+   MFEM_VERIFY(fec, "Only H1_FECollection is supported for the surface fitting "
+               "grid function.");
+
    delete surf_fit_gf;
    surf_fit_gf = new GridFunction(s0);
    surf_fit_gf->CountElementsPerVDof(surf_fit_dof_count);
@@ -2987,12 +2996,24 @@ void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &pos,
 void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
                                            const Array<bool> &smarker,
                                            Coefficient &coeff,
-                                           AdaptivityEvaluator &ae)
+                                           AdaptivityEvaluator &ae,
+                                           AdaptivityEvaluator *aegrad,
+                                           AdaptivityEvaluator *aehess)
 {
    // To have both we must duplicate the markers.
    MFEM_VERIFY(surf_fit_pos == NULL,
                "Using both fitting approaches is not supported.");
 
+   const int dim = s0.FESpace()->GetMesh()->Dimension();
+   ParMesh *pmesh = s0.ParFESpace()->GetParMesh();
+   MFEM_VERIFY(pmesh->GetNodes()->Size() == dim*s0.Size(),
+               "Mesh and level-set polynomial order must be the same.");
+   const H1_FECollection *fec = dynamic_cast<const H1_FECollection *>
+                                (s0.FESpace()->FEColl());
+   MFEM_VERIFY(fec, "Only H1_FECollection is supported for the surface fitting "
+               "grid function.");
+
+
    delete surf_fit_gf;
    surf_fit_gf = new GridFunction(s0);
    s0.CountElementsPerVDof(surf_fit_dof_count);
@@ -3000,11 +3021,80 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
    surf_fit_coeff = &coeff;
    surf_fit_eval = &ae;
 
-   surf_fit_eval->SetParMetaInfo(*s0.ParFESpace()->GetParMesh(),
-                                 *s0.ParFESpace());
+   surf_fit_eval->SetParMetaInfo(*pmesh, *s0.ParFESpace());
    surf_fit_eval->SetInitialField
    (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
-   surf_fit_gf_bg = false;
+
+   if (!aegrad) { return; }
+
+   MFEM_VERIFY(aehess, "AdaptivityEvaluator for Hessians must be provided too.");
+
+   ParFiniteElementSpace *fes = s0.ParFESpace();
+
+   // FE space for gradients.
+   delete surf_fit_grad;
+   H1_FECollection *fec_grad = new H1_FECollection(fec->GetOrder(), dim,
+                                                   fec->GetBasisType());
+   ParFiniteElementSpace *fes_grad = new ParFiniteElementSpace(pmesh, fec_grad,
+                                                               dim);
+   // Initial gradients.
+   surf_fit_grad = new GridFunction(fes_grad);
+   surf_fit_grad->MakeOwner(fec_grad);
+   for (int d = 0; d < dim; d++)
+   {
+      ParGridFunction surf_fit_grad_comp(fes, surf_fit_grad->GetData()+d*s0.Size());
+      s0.GetDerivative(1, d, surf_fit_grad_comp);
+   }
+   surf_fit_eval_grad = aegrad;
+   surf_fit_eval_grad->SetParMetaInfo(*pmesh, *fes_grad);
+   surf_fit_eval_grad->SetInitialField(*pmesh->GetNodes(), *surf_fit_grad);
+
+   // FE space for Hessians.
+   delete surf_fit_hess;
+   H1_FECollection *fec_hess = new H1_FECollection(fec->GetOrder(), dim,
+                                                   fec->GetBasisType());
+   ParFiniteElementSpace *fes_hess = new ParFiniteElementSpace(pmesh, fec_hess,
+                                                               dim*dim);
+   // Initial Hessians.
+   surf_fit_hess = new GridFunction(fes_hess);
+   surf_fit_hess->MakeOwner(fec_hess);
+   int id = 0;
+   for (int d = 0; d < dim; d++)
+   {
+      for (int idir = 0; idir < dim; idir++)
+      {
+         ParGridFunction surf_fit_grad_comp(fes,
+                                            surf_fit_grad->GetData()+d*s0.Size());
+         ParGridFunction surf_fit_hess_comp(fes,
+                                            surf_fit_hess->GetData()+id*s0.Size());
+         surf_fit_grad_comp.GetDerivative(1, idir, surf_fit_hess_comp);
+         id++;
+      }
+   }
+   surf_fit_eval_hess = aehess;
+   surf_fit_eval_hess->SetParMetaInfo(*pmesh, *fes_hess);
+   surf_fit_eval_hess->SetInitialField(*pmesh->GetNodes(), *surf_fit_hess);
+
+   // Store DOF indices that are marked for fitting. Used to reduce work for
+   // transferring information between source/background and current mesh.
+   surf_fit_marker_dof_index.SetSize(0);
+#ifdef MFEM_USE_GSLIB
+   if (dynamic_cast<InterpolatorFP *>(surf_fit_eval) &&
+       dynamic_cast<InterpolatorFP *>(surf_fit_eval_grad) &&
+       dynamic_cast<InterpolatorFP *>(surf_fit_eval_hess))
+   {
+      for (int i = 0; i < surf_fit_marker->Size(); i++)
+      {
+         if ((*surf_fit_marker)[i] == true)
+         {
+            surf_fit_marker_dof_index.Append(i);
+         }
+      }
+   }
+#endif
+
+   *surf_fit_grad = 0.0;
+   *surf_fit_hess = 0.0;
 }
 
 void TMOP_Integrator::EnableSurfaceFittingFromSource(
@@ -3022,12 +3112,10 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
    // Setup for level set function
    delete surf_fit_gf;
    surf_fit_gf = new GridFunction(s0);
-   *surf_fit_gf = 0.0;
    surf_fit_marker = &smarker;
    surf_fit_coeff = &coeff;
    surf_fit_eval = &ae;
 
-   surf_fit_gf_bg = true;
    surf_fit_eval->SetParMetaInfo(*s_bg.ParFESpace()->GetParMesh(),
                                  *s_bg.ParFESpace());
    surf_fit_eval->SetInitialField
@@ -3044,11 +3132,11 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
    delete surf_fit_grad;
    surf_fit_grad = new GridFunction(s0_grad);
    *surf_fit_grad = 0.0;
-   surf_fit_eval_bg_grad = &age;
-   surf_fit_eval_bg_hess = &ahe;
-   surf_fit_eval_bg_grad->SetParMetaInfo(*s_bg_grad.ParFESpace()->GetParMesh(),
-                                         *s_bg_grad.ParFESpace());
-   surf_fit_eval_bg_grad->SetInitialField
+   surf_fit_eval_grad = &age;
+   surf_fit_eval_hess = &ahe;
+   surf_fit_eval_grad->SetParMetaInfo(*s_bg_grad.ParFESpace()->GetParMesh(),
+                                      *s_bg_grad.ParFESpace());
+   surf_fit_eval_grad->SetInitialField
    (*s_bg_grad.FESpace()->GetMesh()->GetNodes(), s_bg_grad);
 
    // Setup for Hessian on background mesh
@@ -3059,9 +3147,9 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
    delete surf_fit_hess;
    surf_fit_hess = new GridFunction(s0_hess);
    *surf_fit_hess = 0.0;
-   surf_fit_eval_bg_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(),
-                                         *s_bg_hess.ParFESpace());
-   surf_fit_eval_bg_hess->SetInitialField
+   surf_fit_eval_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(),
+                                      *s_bg_hess.ParFESpace());
+   surf_fit_eval_hess->SetInitialField
    (*s_bg_hess.FESpace()->GetMesh()->GetNodes(), s_bg_hess);
 
    // Count number of zones that share each of the DOFs
@@ -3866,7 +3954,7 @@ void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x,
 
    Vector sigma_e(dof_s);
    DenseMatrix surf_fit_grad_e(dof_s, dim);
-   if (surf_fit_gf || surf_fit_gf_bg)
+   if (surf_fit_gf)
    {
       surf_fit_gf->GetSubVector(vdofs, sigma_e);
 
@@ -3874,7 +3962,7 @@ void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x,
       // The FE coefficients of the gradient go in surf_fit_grad_e.
       Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
       DenseMatrix grad_phys; // This will be (dof x dim, dof).
-      if (surf_fit_gf_bg)
+      if (surf_fit_grad)
       {
          surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
          surf_fit_grad->GetSubVector(dofs, grad_ptr);
@@ -3948,7 +4036,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
    Vector sigma_e(dof_s);
    DenseMatrix surf_fit_grad_e(dof_s, dim);
    DenseMatrix surf_fit_hess_e(dof_s, dim*dim);
-   if (surf_fit_gf || surf_fit_gf_bg)
+   if (surf_fit_gf)
    {
       surf_fit_gf->GetSubVector(vdofs, sigma_e);
 
@@ -3956,7 +4044,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
       // The FE coefficients of the gradient go in surf_fit_grad_e.
       Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
       DenseMatrix grad_phys; // This will be (dof x dim, dof).
-      if (surf_fit_gf_bg)
+      if (surf_fit_grad)
       {
          surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
          surf_fit_grad->GetSubVector(dofs, grad_ptr);
@@ -3970,7 +4058,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
       // Project the Hessian of sigma in the same space.
       // The FE coefficients of the Hessian go in surf_fit_hess_e.
       Vector hess_ptr(surf_fit_hess_e.GetData(), dof_s*dim*dim);
-      if (surf_fit_gf_bg)
+      if (surf_fit_hess)
       {
          surf_fit_hess->FESpace()->GetElementVDofs(el_id, dofs);
          surf_fit_hess->GetSubVector(dofs, hess_ptr);
@@ -3997,7 +4085,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
       Tpr.SetIntPoint(&ip);
       real_t w = surf_fit_normal * surf_fit_coeff->Eval(Tpr, ip);
 
-      if (surf_fit_gf || surf_fit_gf_bg)
+      if (surf_fit_gf)
       {
          Vector gg_ptr(surf_fit_hess_s.GetData(), dim * dim);
          surf_fit_hess_e.GetRow(s, gg_ptr);
@@ -4379,6 +4467,130 @@ void TMOP_Integrator::ComputeMinJac(const Vector &x,
    dx = detv_avg_min / dxscale;
 }
 
+void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x,
+                                                         int new_x_ordering)
+{
+   if (!surf_fit_gf) { return; }
+
+   if (surf_fit_marker_dof_index.Size())
+   {
+      // Interpolate information only at DOFs marked for fitting.
+      const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension();
+      const int cnt = surf_fit_marker_dof_index.Size();
+      const int total_cnt = new_x.Size()/dim;
+      Vector new_x_sorted(cnt*dim);
+      if (new_x_ordering == 0)
+      {
+         for (int d = 0; d < dim; d++)
+         {
+            for (int i = 0; i < cnt; i++)
+            {
+               int dof_index = surf_fit_marker_dof_index[i];
+               new_x_sorted(i + d*cnt) = new_x(dof_index + d*total_cnt);
+            }
+         }
+      }
+      else
+      {
+         for (int i = 0; i < cnt; i++)
+         {
+            int dof_index = surf_fit_marker_dof_index[i];
+            for (int d = 0; d < dim; d++)
+            {
+               new_x_sorted(d + i*dim) = new_x(d + dof_index*dim);
+            }
+         }
+      }
+
+      // Interpolate values of the LS.
+      Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int;
+      surf_fit_eval->ComputeAtNewPosition(new_x_sorted, surf_fit_gf_int,
+                                          new_x_ordering);
+      for (int i = 0; i < cnt; i++)
+      {
+         int dof_index = surf_fit_marker_dof_index[i];
+         (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i);
+      }
+
+      // Interpolate gradients of the LS.
+      surf_fit_eval_grad->ComputeAtNewPosition(new_x_sorted, surf_fit_grad_int,
+                                               new_x_ordering);
+      // Assumes surf_fit_grad and surf_fit_gf share the same space
+      const int grad_dim = surf_fit_grad->VectorDim();
+      const int grad_cnt = surf_fit_grad->Size()/grad_dim;
+      if (surf_fit_grad->FESpace()->GetOrdering() == Ordering::byNODES)
+      {
+         for (int d = 0; d < grad_dim; d++)
+         {
+            for (int i = 0; i < cnt; i++)
+            {
+               int dof_index = surf_fit_marker_dof_index[i];
+               (*surf_fit_grad)[dof_index + d*grad_cnt] =
+                  surf_fit_grad_int(i + d*cnt);
+            }
+         }
+      }
+      else
+      {
+         for (int i = 0; i < cnt; i++)
+         {
+            int dof_index = surf_fit_marker_dof_index[i];
+            for (int d = 0; d < grad_dim; d++)
+            {
+               (*surf_fit_grad)[dof_index*grad_dim + d] =
+                  surf_fit_grad_int(i*grad_dim + d);
+            }
+         }
+      }
+
+      // Interpolate Hessians of the LS.
+      surf_fit_eval_hess->ComputeAtNewPosition(new_x_sorted, surf_fit_hess_int,
+                                               new_x_ordering);
+      // Assumes surf_fit_hess and surf_fit_gf share the same space
+      const int hess_dim = surf_fit_hess->VectorDim();
+      const int hess_cnt = surf_fit_hess->Size()/hess_dim;
+      if (surf_fit_hess->FESpace()->GetOrdering() == Ordering::byNODES)
+      {
+         for (int d = 0; d < hess_dim; d++)
+         {
+            for (int i = 0; i < cnt; i++)
+            {
+               int dof_index = surf_fit_marker_dof_index[i];
+               (*surf_fit_hess)[dof_index + d*hess_cnt] =
+                  surf_fit_hess_int(i + d*cnt);
+            }
+         }
+      }
+      else
+      {
+         for (int i = 0; i < cnt; i++)
+         {
+            int dof_index = surf_fit_marker_dof_index[i];
+            for (int d = 0; d < hess_dim; d++)
+            {
+               (*surf_fit_hess)[dof_index*hess_dim + d] =
+                  surf_fit_hess_int(i*hess_dim + d);
+            }
+         }
+      }
+
+   }
+   else
+   {
+      surf_fit_eval->ComputeAtNewPosition(new_x, *surf_fit_gf, new_x_ordering);
+      if (surf_fit_eval_grad)
+      {
+         surf_fit_eval_grad->ComputeAtNewPosition(new_x, *surf_fit_grad,
+                                                  new_x_ordering);
+      }
+      if (surf_fit_eval_hess)
+      {
+         surf_fit_eval_hess->ComputeAtNewPosition(new_x, *surf_fit_hess,
+                                                  new_x_ordering);
+      }
+   }
+}
+
 void TMOP_Integrator::
 UpdateAfterMeshPositionChange(const Vector &x_new,
                               const FiniteElementSpace &x_fes)
@@ -4409,112 +4621,11 @@ UpdateAfterMeshPositionChange(const Vector &x_new,
       adapt_lim_eval->ComputeAtNewPosition(x_new, *adapt_lim_gf, ordering);
    }
 
-   // Update surf_fit_gf if surface fitting is enabled.
+   // Update surf_fit_gf (and optionally its gradients) if surface
+   // fitting is enabled.
    if (surf_fit_gf)
    {
-      if (surf_fit_gf_bg)
-      {
-         // Interpolate information for only DOFs marked for fitting.
-         const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension();
-         const int cnt = surf_fit_marker_dof_index.Size();
-         const int total_cnt = x_new.Size()/dim;
-         Vector new_x_sorted(cnt*dim);
-         if (ordering == 0)
-         {
-            for (int d = 0; d < dim; d++)
-            {
-               for (int i = 0; i < cnt; i++)
-               {
-                  int dof_index = surf_fit_marker_dof_index[i];
-                  new_x_sorted(i + d*cnt) = x_new(dof_index + d*total_cnt);
-               }
-            }
-         }
-         else
-         {
-            for (int i = 0; i < cnt; i++)
-            {
-               int dof_index = surf_fit_marker_dof_index[i];
-               for (int d = 0; d < dim; d++)
-               {
-                  new_x_sorted(d + i*dim) = x_new(d + dof_index*dim);
-               }
-            }
-         }
-
-         Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int;
-         surf_fit_eval->ComputeAtNewPosition(
-            new_x_sorted, surf_fit_gf_int, ordering);
-         for (int i = 0; i < cnt; i++)
-         {
-            int dof_index = surf_fit_marker_dof_index[i];
-            (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i);
-         }
-
-         surf_fit_eval_bg_grad->ComputeAtNewPosition(
-            new_x_sorted, surf_fit_grad_int, ordering);
-         // Assumes surf_fit_grad and surf_fit_gf share the same space
-         const int grad_dim = surf_fit_grad->VectorDim();
-         const int grad_cnt = surf_fit_grad->Size()/grad_dim;
-         if (surf_fit_grad->FESpace()->GetOrdering() == Ordering::byNODES)
-         {
-            for (int d = 0; d < grad_dim; d++)
-            {
-               for (int i = 0; i < cnt; i++)
-               {
-                  int dof_index = surf_fit_marker_dof_index[i];
-                  (*surf_fit_grad)[dof_index + d*grad_cnt] =
-                     surf_fit_grad_int(i + d*cnt);
-               }
-            }
-         }
-         else
-         {
-            for (int i = 0; i < cnt; i++)
-            {
-               int dof_index = surf_fit_marker_dof_index[i];
-               for (int d = 0; d < grad_dim; d++)
-               {
-                  (*surf_fit_grad)[dof_index*dim + d] =
-                     surf_fit_grad_int(i*dim + d);
-               }
-            }
-         }
-
-         surf_fit_eval_bg_hess->ComputeAtNewPosition(
-            new_x_sorted, surf_fit_hess_int, ordering);
-         // Assumes surf_fit_hess and surf_fit_gf share the same space
-         const int hess_dim = surf_fit_hess->VectorDim();
-         const int hess_cnt = surf_fit_hess->Size()/hess_dim;
-         if (surf_fit_hess->FESpace()->GetOrdering() == Ordering::byNODES)
-         {
-            for (int d = 0; d < hess_dim; d++)
-            {
-               for (int i = 0; i < cnt; i++)
-               {
-                  int dof_index = surf_fit_marker_dof_index[i];
-                  (*surf_fit_hess)[dof_index + d*hess_cnt] =
-                     surf_fit_hess_int(i + d*cnt);
-               }
-            }
-         }
-         else
-         {
-            for (int i = 0; i < cnt; i++)
-            {
-               int dof_index = surf_fit_marker_dof_index[i];
-               for (int d = 0; d < hess_dim; d++)
-               {
-                  (*surf_fit_hess)[dof_index*dim + d] =
-                     surf_fit_hess_int(i*dim + d);
-               }
-            }
-         }
-      }
-      else
-      {
-         surf_fit_eval->ComputeAtNewPosition(x_new, *surf_fit_gf, ordering);
-      }
+      RemapSurfaceFittingLevelSetAtNodes(x_new, ordering);
    }
 }
 

fem/tmop.hpp

@@ -1784,12 +1784,11 @@ protected:
    // Fitting to given physical positions.
    TMOP_QuadraticLimiter *surf_fit_limiter; // Owned. Created internally.
    const GridFunction *surf_fit_pos;        // Not owned. Positions to fit.
-   real_t surf_fit_normal;
-   bool surf_fit_gf_bg;
-   GridFunction *surf_fit_grad, *surf_fit_hess;
-   AdaptivityEvaluator *surf_fit_eval_bg_grad, *surf_fit_eval_bg_hess;
-   Array<int> surf_fit_dof_count;
-   Array<int> surf_fit_marker_dof_index;
+   real_t surf_fit_normal;                  // Normalization factor.
+   GridFunction *surf_fit_grad, *surf_fit_hess; // Owned. Created internally.
+   AdaptivityEvaluator *surf_fit_eval_grad, *surf_fit_eval_hess; // Not owned.
+   Array<int> surf_fit_dof_count;            // Number of dofs per node.
+   Array<int> surf_fit_marker_dof_index;     // Indices of nodes to fit.
 
    DiscreteAdaptTC *discr_tc;
 
@@ -1985,6 +1984,10 @@ protected:
    real_t ComputeUntanglerMaxMuBarrier(const Vector &x,
                                        const FiniteElementSpace &fes);
 
+   // Remaps the internal surface fitting gridfunction object at provided
+   // locations.
+   void RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x,
+                                           int new_x_ordering);
 public:
    /** @param[in] m    TMOP_QualityMetric for r-adaptivity (not owned).
        @param[in] tc   Target-matrix construction algorithm to use (not owned).
@@ -2000,9 +2003,8 @@ public:
         surf_fit_marker(NULL), surf_fit_coeff(NULL),
         surf_fit_gf(NULL), surf_fit_eval(NULL),
         surf_fit_limiter(NULL), surf_fit_pos(NULL),
-        surf_fit_normal(1.0),
-        surf_fit_gf_bg(false), surf_fit_grad(NULL), surf_fit_hess(NULL),
-        surf_fit_eval_bg_grad(NULL), surf_fit_eval_bg_hess(NULL),
+        surf_fit_normal(1.0), surf_fit_grad(NULL), surf_fit_hess(NULL),
+        surf_fit_eval_grad(NULL), surf_fit_eval_hess(NULL),
         discr_tc(dynamic_cast<DiscreteAdaptTC *>(tc)),
         fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false)
    { PA.enabled = false; }
@@ -2103,9 +2105,15 @@ public:
 
 #ifdef MFEM_USE_MPI
    /// Parallel support for surface fitting to the zero level set of a function.
+   /// Here, we add two optional inputs: @a aegrad and @a aehess. When provided,
+   /// the first and second derivative of the input level set are computed on
+   /// the initial mesh, and @a aegrad and @a aehess are used to remap grad_s(x)
+   /// from grad_s0(x0) and hess_s(x) from hess_s0(x0), respectively.
    void EnableSurfaceFitting(const ParGridFunction &s0,
                              const Array<bool> &smarker, Coefficient &coeff,
-                             AdaptivityEvaluator &ae);
+                             AdaptivityEvaluator &ae,
+                             AdaptivityEvaluator *aegrad = NULL,
+                             AdaptivityEvaluator *aehess = NULL);
 
    /** @brief Fitting of certain DOFs in the current mesh to the zero level set
        of a function defined on another (finer) source mesh.

general/socketstream.cpp

@@ -134,7 +134,7 @@ int socketbuf::open(const char hostname[], int port)
       {
          closesocket(socket_descriptor);
          socket_descriptor = -2;
-         return -1;
+         continue;
       }
 #endif
 
@@ -148,7 +148,7 @@ int socketbuf::open(const char hostname[], int port)
    }
 
    freeaddrinfo(res);
-   return 0;
+   return (socket_descriptor < 0) ? -1 : 0;
 }
 
 int socketbuf::close()

linalg/densemat.cpp

@@ -532,6 +532,69 @@ MatrixInverse *DenseMatrix::Inverse() const
    return new DenseMatrixInverse(*this);
 }
 
+void DenseMatrix::Exponential()
+{
+   MFEM_ASSERT(Height() == Width() && Height() <= 2,
+               "The matrix must be square and "
+               << "of size less than or equal to 2."
+               << "  Height() = " << Height()
+               << ", Width() = " << Width());
+
+   switch (Height())
+   {
+      case 1:
+      {
+         data[0] = std::exp(data[0]);
+         break;
+      }
+      case 2:
+      {
+         /// Formulas from Corollary 2.4 of doi:10.1109/9.233156
+         /// Note typo in the paper, in the prefactor in the equation under (i).
+         const real_t a = data[0];
+         const real_t b = data[1];
+         const real_t c = data[2];
+         const real_t d = data[3];
+         const real_t e = (a - d)*(a - d) + 4*b*c;
+         const real_t f = std::exp((a + d)/2.0);
+         const real_t g = std::sqrt(std::abs(e)) / 2.0;
+
+         if (e == 0)
+         {
+            data[0] = 1.0 + (a - d)/2.0;
+            data[3] = 1.0 - (a - d)/2.0;
+         }
+         else if (e > 0)
+         {
+            data[0] = std::cosh(g) + (a - d)/2 * std::sinh(g) / g;
+            data[1] = b * std::sinh(g) / g;
+            data[2] = c * std::sinh(g) / g;
+            data[3] = std::cosh(g) - (a - d)/2 * std::sinh(g) / g;
+         }
+         else
+         {
+            data[0] = std::cos(g) + (a - d)/2 * std::sin(g) / g;
+            data[1] = b * std::sin(g) / g;
+            data[2] = c * std::sin(g) / g;
+            data[3] = std::cos(g) - (a - d)/2 * std::sin(g) / g;
+         }
+         for (int i = 0; i < 4; i++)
+         {
+            data[i] *= f;
+         }
+         break;
+      }
+      case 3:
+      {
+         MFEM_ABORT("3x3 matrices are not currently supported");
+      }
+      default:
+      {
+         MFEM_ABORT("Only 1x1 and 2x2 matrices are currently supported");
+      }
+   }
+}
+
 real_t DenseMatrix::Det() const
 {
    MFEM_ASSERT(Height() == Width() && Height() > 0,

linalg/densemat.hpp

@@ -207,6 +207,10 @@ public:
    /// Replaces the current matrix with its square root inverse
    void SquareRootInverse();
 
+   /// Replaces the current matrix with its exponential
+   /// (currently only supports 2x2 matrices)
+   void Exponential();
+
    /// Calculates the determinant of the matrix
    /// (optimized for 2x2, 3x3, and 4x4 matrices)
    real_t Det() const;

mesh/mesh_readers.cpp

@@ -1151,15 +1151,24 @@ void Mesh::ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf,
    }
    if (cells_xml == NULL) { MFEM_ABORT(erstr); }
 
-   // Read the element attributes, which are stored as CellData named "material"
+   // Read the element attributes, which are stored as CellData named either
+   // "material" or "attribute". We prioritize "material" over "attribute" for
+   // backwards compatibility.
    Array<int> cell_attributes;
+   bool found_attributes = false;
    for (const XMLElement *cell_data_xml = piece->FirstChildElement();
         cell_data_xml != NULL;
         cell_data_xml = cell_data_xml->NextSiblingElement())
    {
-      if (StringCompare(cell_data_xml->Name(), "CellData")
-          && StringCompare(cell_data_xml->Attribute("Scalars"), "material"))
+      const bool is_cell_data =
+         StringCompare(cell_data_xml->Name(), "CellData");
+      const bool is_material =
+         StringCompare(cell_data_xml->Attribute("Scalars"), "material");
+      const bool is_attribute =
+         StringCompare(cell_data_xml->Attribute("Scalars"), "attribute");
+      if (is_cell_data && (is_material || (is_attribute && !found_attributes)))
       {
+         found_attributes = true;
          const XMLElement *data_xml = cell_data_xml->FirstChildElement();
          if (data_xml != NULL && StringCompare(data_xml->Name(), "DataArray"))
          {
@@ -1274,6 +1283,7 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
    // Read the cell materials
    // bool found_material = false;
    Array<int> cell_attributes;
+   bool found_attributes = false;
    while ((input.good()))
    {
       getline(input, buff);
@@ -1281,8 +1291,10 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
       {
          break; // We have entered the POINT_DATA block. Quit.
       }
-      else if (buff.rfind("SCALARS material") == 0)
+      else if (buff.rfind("SCALARS material") == 0 ||
+               (buff.rfind("SCALARS attribute") == 0 && !found_attributes))
       {
+         found_attributes = true;
          getline(input, buff); // LOOKUP_TABLE default
          if (buff.rfind("LOOKUP_TABLE default") != 0)
          {

miniapps/contact/contact.cpp

@@ -0,0 +1,637 @@
+//                               Parallel contact example
+// mpirun -np 4 ./contact -ls 2 -sr 1 -testno 4
+// CG iteration numbers            = 105 114 116 115 113 109 113 108 107 114 206 236 268 435 987
+
+// mpirun -np 4 ./contact -ls 2 -sr 0 -testno 5
+// CG iteration numbers            = 106 116 116 116 115 113 107 107 128 131 531 1437 1318
+
+// mpirun -np 4 ./contact -ls 2 -sr 0 -testno 6
+// CG iteration numbers            = 18 18 18 18 18 17 17 21 22 46 52 53
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "ipsolver/ParIPsolver.hpp"
+
+using namespace std;
+using namespace mfem;
+
+double GetBdrElementVolume(int i, Mesh & mesh)
+{
+   ElementTransformation *et = mesh.GetBdrElementTransformation(i);
+   const IntegrationRule &ir = IntRules.Get(mesh.GetBdrElementGeometry(i),
+                                            et->OrderJ());
+   double volume = 0.0;
+   for (int j = 0; j < ir.GetNPoints(); j++)
+   {
+      const IntegrationPoint &ip = ir.IntPoint(j);
+      et->SetIntPoint(&ip);
+      volume += ip.weight * et->Weight();
+   }
+
+   return volume;
+}
+
+
+double GetBdrArea(int bdrattr, Mesh&mesh)
+{
+   double area = 0.0;
+   for (int i = 0; i<mesh.GetNBE(); i++)
+   {
+      if (mesh.GetBdrAttribute(i) == bdrattr)
+      {
+         area += GetBdrElementVolume(i,mesh);
+      }
+   }
+
+   MPI_Allreduce(MPI_IN_PLACE,&area,1, MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+   return area;
+}
+
+void OutputData(ostringstream & file_name, double E0, double Ef, int dofs, int constr, int optit, const Array<int> & iters)
+{
+   file_name << ".csv";
+   std::ofstream outputfile(file_name.str().c_str());
+   //if (!outputfile.is_open()) 
+   //{
+   //   MFEM_ABORT("Failed to open file for writing.\n");
+   //}
+   outputfile << "Initial Energy objective        = " << E0 << endl;
+   outputfile << "Final Energy objective          = " << Ef << endl;
+   outputfile << "Global number of dofs           = " << dofs << endl;
+   outputfile << "Global number of constraints    = " << constr << endl;
+   outputfile << "Optimizer number of iterations  = " << optit << endl;
+   outputfile << "CG iteration numbers            = "; iters.Print(outputfile, iters.Size());
+   outputfile << "OptimizerIteration,CGIterations" << endl;
+   for (int i = 0; i< iters.Size(); i++)
+   {
+      outputfile << i+1 <<","<< iters[i] << endl;
+   }
+   outputfile.close();   
+   std::cout << " Data has been written to " << file_name.str().c_str() << endl;
+}
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
+   Hypre::Init();
+
+   int order = 1;
+   int sref = 1;
+   int pref = 0;
+   Array<int> attr;
+   Array<int> m_attr;
+   bool visualization = true;
+   bool paraview = false;
+   int paraview_plot_every = 1;
+   int SQPrepeat = 1;
+   double linsolverrtol = 1e-10;
+   double linsolveratol = 1e-12;
+   int relax_type = 8;
+   double optimizer_tol = 1e-6;
+   int optimizer_maxit = 20;
+   int linsolver = 2; // PCG  - AMG
+   bool elast = false;
+   bool nocontact = false;
+   int testNo = -1; // 0-6
+   int nsteps = 1;
+   bool outputfiles = false;
+   bool doublepass = false;
+   // 1. Parse command-line options.
+   OptionsParser args(argc, argv);
+   args.AddOption(&testNo, "-testno", "--test-number",
+                  "Choice of test problem:"
+                  "-1: default (original 2 block problem)"
+                  "0: not implemented yet"
+                  "1: not implemented yet"
+                  "2: not implemented yet"
+                  "3: not implemented yet"
+                  "4: two block problem - diablo"
+                  "41: two block problem - twisted"
+                  "5: ironing problem"
+                  "51: ironing problem extended"
+                  "6: nested spheres problem");
+   args.AddOption(&attr, "-at", "--attributes-surf",
+                  "Attributes of boundary faces on contact surface for mesh 2.");
+   args.AddOption(&sref, "-sr", "--serial-refinements",
+                  "Number of uniform refinements.");                  
+   args.AddOption(&nsteps, "-nsteps", "--nsteps",
+                  "Number of steps.");
+   args.AddOption(&pref, "-pr", "--parallel-refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&linsolverrtol, "-srtol", "--solver-rel-tol",
+                  "Linear Solver Relative Tolerance.");
+   args.AddOption(&linsolveratol, "-satol", "--solver-abs-tol",
+                  "Linear Solver Abs Tolerance.");
+   args.AddOption(&elast, "-elast", "--elast", "-no-elast",
+                  "--no-elast",
+                  "Enable or disable AMG Elasticity options.");
+   args.AddOption(&nocontact, "-nocontact", "--nocontact", "-no-nocontact",
+                  "--no-nocontact",
+                  "Enable or disable AMG solve with no contact for testing.");
+   args.AddOption(&doublepass, "-doublepass", "--double-pass", "-singlepass",
+                  "--single-pass",
+                  "Enable or disable double pass for contact constraints.");                  
+   args.AddOption(&optimizer_tol, "-otol", "--optimizer-tol",
+                  "Interior Point Solver Tolerance.");
+   args.AddOption(&optimizer_maxit, "-omaxit", "--optimizer-maxit",
+                  "Interior Point Solver maximum number of iterations.");
+   args.AddOption(&relax_type, "-rt", "--relax-type",
+                  "Selection of Smoother for AMG");
+   args.AddOption(&linsolver, "-ls", "--linear-solver",
+                  "Selection of inner linear solver:" 
+                  "0: mumps," 
+                  "1: mumps-reduced,"
+                  "2: PCG-AMG-reduced,"
+                  "3: PCG- with block-diag(AMG,direct solver)"
+                  "4: with static cond of contact dofs");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
+                  "--no-paraview",
+                  "Enable or disable ParaView visualization.");
+   args.AddOption(&paraview_plot_every, "-plot_every", "--plot-every",
+		   "Output every plot_every pseudotimesteps as a paraview file");
+   args.AddOption(&SQPrepeat, "-nSQPrepeat", "--nSQP-repeats", "Number of times to relinearize and resolve the SQP before incremenetally updating forcing and boundary terms");   
+   args.AddOption(&outputfiles, "-out", "--output", "-no-out",
+                  "--no-ouput",
+                  "Enable or disable ouput to files.");                  
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   if (Mpi::Root())
+   {
+      mfem::out << "Solving test problem number: " << testNo << endl;
+   }
+
+   const char *mesh_file = nullptr;
+
+   switch (testNo)
+   {
+      case -1:
+         mesh_file = "meshes/two-block.mesh";
+         break;
+      case 0:
+      case 1:
+      case 2:
+      case 3:
+      {
+         MFEM_ABORT("Problem not implemented yet");
+         break;
+      }
+      case 4:
+         mesh_file = "meshes/Test4.mesh";
+         break;
+      case 40:
+         mesh_file = "meshes/Test40.mesh";
+         break;   
+      case 41:
+         mesh_file = "meshes/Test41.mesh";
+         break;
+      case 42:
+         mesh_file = "meshes/Test42.mesh";
+         break;         
+      case 5:
+         mesh_file = "meshes/Test5.mesh";
+         break;
+      case 51:
+         mesh_file = "meshes/Test51.mesh";
+         break;
+      case 6:
+         mesh_file = "meshes/Test6.mesh";
+         break;
+      case 61:
+         // Something wrong with this mesh
+         mesh_file = "meshes/Test61.mesh";
+         break;
+      case 62:
+         mesh_file = "meshes/Test62.mesh";
+         break;         
+      default:
+         MFEM_ABORT("Should be unreachable");
+         break;
+   }
+
+   Mesh * mesh = new Mesh(mesh_file,1);
+   for (int i = 0; i<sref; i++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   ParMesh * pmesh = new ParMesh(MPI_COMM_WORLD,*mesh);
+
+   for (int i = 0; i<pref; i++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   Array<int> ess_bdr_attr;
+   Array<int> ess_bdr_attr_comp;
+   if (testNo == 6 || testNo == 61)
+   {
+      ess_bdr_attr.Append(1); ess_bdr_attr_comp.Append(1);
+      ess_bdr_attr.Append(2); ess_bdr_attr_comp.Append(2);
+      ess_bdr_attr.Append(4); ess_bdr_attr_comp.Append(0);
+      ess_bdr_attr.Append(5); ess_bdr_attr_comp.Append(-1);
+   }
+   else if (testNo == 62)
+   {
+      ess_bdr_attr.Append(4); ess_bdr_attr_comp.Append(0);
+      ess_bdr_attr.Append(5); ess_bdr_attr_comp.Append(-1);
+   }
+   else if (testNo == 40)
+   {
+      ess_bdr_attr.Append(1); ess_bdr_attr_comp.Append(-1);
+      ess_bdr_attr.Append(10); ess_bdr_attr_comp.Append(-1);
+   }
+   else
+   {
+      ess_bdr_attr.Append(2); ess_bdr_attr_comp.Append(-1);
+      ess_bdr_attr.Append(6); ess_bdr_attr_comp.Append(-1);
+   }
+   ParElasticityProblem * prob = new ParElasticityProblem(pmesh,
+                                                          ess_bdr_attr,ess_bdr_attr_comp,
+                                                          order);
+   Vector lambda(prob->GetMesh()->attributes.Max());
+   Vector mu(prob->GetMesh()->attributes.Max());
+
+   if (testNo == -1 )
+   {
+      lambda = 57.6923076923;
+      mu = 38.4615384615;
+   }
+   else if (testNo == 6 || testNo == 61 || testNo == 62)
+   {
+      lambda = (1000*0.3)/(1.3*0.4);
+      mu = 500/(1.3);
+   }
+   else
+   {
+      //lambda = 57.6923076923;
+      //mu = 38.4615384615;
+      //lambda = 0.499 / (1.499 * 0.002);
+      //mu = 1. / (2. * 1.499);
+      lambda[0] = 0.499/(1.499*0.002);
+      lambda[1] = 0.0;
+      mu[0]     = 1. / (2. * 1.499);
+      mu[1]     = 500.;
+   }
+
+   prob->SetLambda(lambda); prob->SetMu(mu);
+
+   int dim = pmesh->Dimension();
+   Vector ess_values(dim);
+   int essbdr_attr;
+   Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+
+   ess_values = 0.0;
+
+
+   double area = GetBdrArea(3,*mesh);
+
+   // ConstantCoefficient one(-area);
+   ConstantCoefficient one(-1.0);
+
+   std::set<int> mortar_attr;
+   std::set<int> nonmortar_attr;
+
+   if (testNo == 6 || testNo == 61)
+   {
+      ess_values = 0.0;
+      ess_bdr = 0;
+      ess_bdr[0] = 1;
+      ess_bdr[1] = 1;
+      ess_bdr[3] = 1;
+      ess_bdr[4] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      ess_bdr = 0;
+      ess_bdr[2] = 1;
+      // prob->SetNeumanPressureData(one,ess_bdr);
+      mortar_attr.insert(6);
+      mortar_attr.insert(9);
+      nonmortar_attr.insert(7);
+      nonmortar_attr.insert(8);
+   }
+   else if(testNo == 62)
+   {
+      ess_values = 0.0;
+      ess_bdr = 0;
+      ess_bdr[3] = 1;
+      ess_bdr[4] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      ess_bdr = 0;
+      ess_bdr[2] = 1;
+      // prob->SetNeumanPressureData(one,ess_bdr);
+      prob->SetNeumanData(0,3,-2.0);
+      mortar_attr.insert(6);
+      mortar_attr.insert(9);
+      nonmortar_attr.insert(7);
+      nonmortar_attr.insert(8);
+   }
+   else
+   {
+      if (testNo == -1 || testNo == 41)
+      {
+         ess_values[0] = 0.1/nsteps;
+      }
+      else
+      {
+         ess_values[2] = 1.0 / 1.4  / nsteps;
+	      //ess_values[2] = 0.25 / nsteps;//1.0/1.4/nsteps;
+         // ess_values[0] = -2.0/nsteps;
+      }
+      essbdr_attr = (testNo == 40) ? 1 : 2;
+      ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+      // prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      essbdr_attr = (testNo == 40) ? 10 : 6;
+      ess_values = 0.0; ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+      // prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      if (testNo == 40)
+      {
+         mortar_attr.insert(4);
+         nonmortar_attr.insert(7);
+      }
+      else
+      {
+         mortar_attr.insert(3);
+         nonmortar_attr.insert(4);
+      }
+   }
+
+   ParFiniteElementSpace * fes = prob->GetFESpace();
+   ParGridFunction x_gf(fes); x_gf = 0.0;
+   ParGridFunction xnew(fes); xnew = 0.0;
+   ParaViewDataCollection * paraview_dc = nullptr;
+   ParMesh pmesh_copy(*pmesh);
+   ParFiniteElementSpace fes_copy(*fes,pmesh_copy);
+   ParGridFunction xcopy_gf(&fes_copy); xcopy_gf = 0.0;
+   
+   if (paraview)
+   {
+      std::ostringstream paraview_file_name;
+      paraview_file_name << "QPContact-Test_" << testNo
+                         << "_par_ref_" << pref
+                         << "_ser_ref_" << sref;
+      paraview_dc = new ParaViewDataCollection(paraview_file_name.str(), &pmesh_copy);
+      paraview_dc->SetPrefixPath("ParaView");
+      paraview_dc->SetLevelsOfDetail(1);
+      paraview_dc->SetDataFormat(VTKFormat::BINARY);
+      paraview_dc->SetHighOrderOutput(true);
+      // paraview_dc->RegisterField("u", &x_gf);
+      paraview_dc->RegisterField("u", &xcopy_gf);
+      paraview_dc->SetCycle(0);
+      paraview_dc->SetTime(double(0));
+      paraview_dc->Save();
+   }
+   socketstream sol_sock;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      sol_sock.open(vishost, visport);
+      sol_sock.precision(8);
+   }
+   // ParGridFunction coords(prob->GetFESpace()); 
+   ParGridFunction ref_coords(prob->GetFESpace()); 
+   ParGridFunction new_coords(prob->GetFESpace()); 
+   pmesh->GetNodes(new_coords);
+   pmesh->GetNodes(ref_coords);
+   
+   Vector xref(x_gf.GetTrueVector().Size());
+
+   HypreParMatrix *dgdu;
+   double p = 1;
+   ConstantCoefficient f(p);
+   
+   // SQPrepeat solves on same problem (forcing/boundary conditions)
+   int Nsteps = nsteps * SQPrepeat;
+   
+   double pseudotime = 0.0;
+   double pseudotimestep = 1.0 / ((double) nsteps);
+   double paraview_time = 0.0;
+   double paraview_subtimestep = pseudotimestep / ((double) SQPrepeat);
+   int paraview_cycle = 1;
+
+   bool QPConverged;
+
+   std::ofstream numConstraintsStream;
+   std::ostringstream numConstraints_file_name;
+   numConstraints_file_name << "data/numConstraints_ref" << sref << ".dat";
+   if (Mpi::Root)
+   {
+      numConstraintsStream.open(numConstraints_file_name.str(), ios::out | ios::trunc);
+   }
+
+   for (int i = 0; i < nsteps; i++)
+   {
+      pseudotime = ((double) (i + 1)) / ((double) nsteps);
+      for (int j = 0; j < SQPrepeat; j++)
+      {
+         paraview_time = pseudotime + j * paraview_subtimestep;
+	 if (testNo == 6)
+         {
+            ess_bdr = 0;
+            ess_bdr[2] = 1;
+            f.constant = -p * pseudotime;
+            prob->SetNeumanPressureData(f,ess_bdr);
+            // prob->SetNeumanData(0,3,-p*(i+1)/nsteps);
+         }
+         else if (testNo == 4 || testNo == 40 || testNo == 5 || testNo == 51)
+         {
+            ess_bdr = 0;
+            essbdr_attr = (testNo == 40) ? 1 : 2;
+            ess_bdr[essbdr_attr-1] = 1;
+            ess_values = 0.0;
+            //ess_values[2] = 4.0 / 7.0 * pseudotime;
+	    //ess_values[2] = 0.25 * pseudotime; //1.0/1.4 * pseudotime;
+            ess_values[2] = 1.0 / 1.4 * pseudotime;
+	    prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+         }
+         else if (testNo == 41)
+         {
+            ess_values = 0.0;
+            ess_values[0] = 0.5 * pseudotime; //0.5/nsteps*(i+1);
+            // ess_values[0] = 0.0;
+            essbdr_attr =  2;
+            ess_bdr[essbdr_attr-1] = 1;
+            prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+            essbdr_attr = 6;
+            ess_values = 0.0; 
+            // ess_values[0] = -0.5/nsteps*(i+1);
+            if (myid == 0)
+            {
+               mfem::out << "ess_values[0] = " << ess_values[0] << endl;
+            }
+            ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+            prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+         }
+
+	 //xref.Set(1.0, x_gf.GetTrueVector());
+         xref = 0.0;
+	 ParContactProblem contact(prob, mortar_attr, nonmortar_attr, &new_coords, doublepass);
+         QPOptParContactProblem qpopt(&contact, xref);
+         int numconstr = contact.GetGlobalNumConstraints();
+         ParInteriorPointSolver optimizer(&qpopt);
+         optimizer.SetTol(optimizer_tol);
+         optimizer.SetMaxIter(optimizer_maxit);
+         optimizer.SetLinearSolver(linsolver);
+         optimizer.SetLinearSolveRelTol(linsolverrtol);
+         optimizer.SetLinearSolveAbsTol(linsolveratol);
+         optimizer.SetLinearSolveRelaxType(relax_type);
+         if (nocontact)
+         {
+            optimizer.EnableNoContactSolve();
+         }
+         if (elast)
+         {
+            optimizer.SetElasticityOptions(prob->GetFESpace());
+         }
+         // ParGridFunction x = prob->GetDisplacementGridFunction();
+         // x.SetTrueVector();
+         // Vector x0 = x.GetTrueVector();
+
+         x_gf.SetTrueVector();
+	 
+
+         Vector x0 = x_gf.GetTrueVector();
+         int ndofs = x0.Size();
+         Vector xf(ndofs); xf = 0.0;
+         optimizer.Mult(x0, xf);
+         QPConverged = optimizer.GetConverged();
+         
+	 /* exit if not converged */
+	 MFEM_VERIFY(QPConverged, "IPM not converged on QP contact problem");
+         
+	 
+         double Einitial = contact.E(x0);
+         double Efinal   = contact.E(xf);
+         Array<int> & CGiterations = optimizer.GetCGIterNumbers();
+         int gndofs = prob->GetGlobalNumDofs();
+         int gnconstraints = contact.GetGlobalNumConstraints();
+   
+	 //std::ofstream xfStream;
+	 //std::ostringstream xf_file_name;
+	 //xf_file_name << "data/xf_" << i << ".dat";
+	 //if (Mpi::Root())
+	 //{
+	 //   xfStream.open(xf_file_name.str(), ios::out | ios::trunc);
+	 //   for (int ii = 0; ii < xf.Size(); ii++)
+	 //   {
+	 //      xfStream << xf(ii) << "\n";
+	 //   }
+	 //   xfStream.close();
+	 //}
+   //if (Mpi::Root)
+   //{
+   //   numConstraintsStream.open(numConstraints_file_name.str(), ios::out | ios::trunc);
+   //}
+
+
+         if (Mpi::Root())
+         {
+            
+            mfem::out << endl;
+            mfem::out << " Initial Energy objective        = " << Einitial << endl;
+            mfem::out << " Final Energy objective          = " << Efinal << endl;
+            mfem::out << " Global number of dofs           = " << gndofs << endl;
+            mfem::out << " Global number of constraints    = " << numconstr << endl;
+            mfem::out << " Optimizer number of iterations  = " <<
+            optimizer.GetNumIterations() << endl;
+            if (linsolver == 2 || linsolver == 3 || linsolver == 4)
+            {
+               mfem::out << " CG iteration numbers            = " ;
+               CGiterations.Print(mfem::out, CGiterations.Size());
+            }
+            if (nocontact)
+            {
+               Array<int> & CGNoContactIterations = optimizer.GetCGNoContactIterNumbers();
+               mfem::out << " CG no Contact iteration numbers = " ;
+               CGNoContactIterations.Print(mfem::out, CGNoContactIterations.Size());
+            }
+            if (outputfiles)
+            {
+               ostringstream file_name;
+               file_name << "output/Testno-"<<testNo<<"-ref-"<<sref+pref << "-step-" << i; 
+               OutputData(file_name, Einitial, Efinal, gndofs,numconstr, optimizer.GetNumIterations(), CGiterations);
+            }
+	    numConstraintsStream << gnconstraints << endl;
+         }
+
+         // Vector X_new(xf.GetData(),fes->GetTrueVSize());
+         // xnew.SetFromTrueDofs(X_new);
+         // x_gf = xnew;
+         x_gf.SetFromTrueDofs(xf);
+         // mfem::out << "x_gf norm = " << x_gf.Norml2() << endl;
+         // cin.get();
+         // pmesh->MoveNodes(xnew);
+         // pmesh_copy.MoveNodes(xnew);
+         // pmesh_copy.MoveNodes(xnew);
+         add(ref_coords,x_gf,new_coords);
+         // mfem::out << " ref_coords norm " << ref_coords.Norml2() << endl;
+         // mfem::out << " x_gf norm " << x_gf.Norml2() << endl;
+         // mfem::out << " new_coords norm " << new_coords.Norml2() << endl;
+         // pmesh_copy.SetNodes(new_coords);
+         pmesh_copy.SetNodes(new_coords);
+         xcopy_gf = x_gf;
+         // pmesh_copy.MoveNodes(x_gf);
+         // pmesh_copy.SetNodes(x_gf);
+         if (paraview && ((i+1) % paraview_plot_every == 0 ))
+         {
+            paraview_cycle += 1;
+            paraview_dc->SetCycle(paraview_cycle) ;
+            paraview_dc->SetTime(paraview_time);
+            paraview_dc->Save();
+         }
+
+         if (visualization)
+         {
+            sol_sock << "parallel " << num_procs << " " << myid << "\n"
+                     << "solution\n" << pmesh_copy << x_gf << flush;
+         
+            if (i == nsteps - 1 && j == SQPrepeat - 1)
+            {
+               pmesh->MoveNodes(x_gf);
+               char vishost[] = "localhost";
+               int  visport   = 19916;
+               socketstream sol_sock1(vishost, visport);
+               sol_sock1 << "parallel " << num_procs << " " << myid << "\n";
+               sol_sock1.precision(8);
+               sol_sock1 << "solution\n" << *pmesh << x_gf << flush;
+            }
+         }
+         if (i == nsteps - 1 && j == SQPrepeat) break;
+
+         prob->UpdateStep();
+         if (testNo == 6 )
+         {
+            double area_new = GetBdrArea(3,*pmesh);
+            if (myid == 0)
+            {
+               mfem::out << "New area = " << area_new << endl;
+            }
+         }
+      }
+   }
+ 
+   if (Mpi::Root)
+   {
+      numConstraintsStream.close();
+   }
+   delete prob;
+   delete pmesh;
+   delete mesh;
+   return 0;
+}

miniapps/contact/ipsolver/ParIPsolver.hpp

@@ -0,0 +1,114 @@
+#include "mfem.hpp"
+#include "../problems/parproblems.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+#ifndef PARIPSOLVER 
+#define PARIPSOLVER
+
+class ParInteriorPointSolver
+{
+protected:
+    QPOptParContactProblem* problem = nullptr;
+    double OptTol;
+    int  max_iter;
+    int  iter=0;
+    double mu_k; // \mu_k
+    Vector lk, zlk;
+
+    double sMax, kSig, tauMin, eta, thetaMin, delta, sTheta, sPhi, kMu, thetaMu;
+    double thetaMax, kSoc, gTheta, gPhi, kEps;
+	
+    // filter
+    Array<double> F1, F2;
+	
+    // quantities computed in lineSearch
+    double alpha, alphaz;
+    double thx0, thxtrial;
+    double phx0, phxtrial;
+    bool descentDirection, switchCondition, sufficientDecrease, lineSearchSuccess, inFilterRegion;
+    double Dxphi0_xhat;
+
+    int dimU, dimM, dimC;
+    int gdimU, gdimM, gdimC;
+    Array<int> block_offsetsumlz, block_offsetsuml, block_offsetsx;
+    Vector ml;
+
+    Vector ckSoc;
+    HypreParMatrix * Huu = nullptr;
+    HypreParMatrix * Hum = nullptr;
+    HypreParMatrix * Hmu = nullptr;
+    HypreParMatrix * Hmm = nullptr;
+    HypreParMatrix * Wmm = nullptr;
+    HypreParMatrix * Ju = nullptr;
+    HypreParMatrix * Jm = nullptr;
+    HypreParMatrix * JuT = nullptr;
+    HypreParMatrix * JmT = nullptr;
+
+    Array<int> cgnum_iterations;
+    Array<int> cgnum_iterations_nocontact;
+    ParFiniteElementSpace *pfes = nullptr;
+    
+    int jOpt;
+    bool converged;
+    
+    int MyRank;
+    bool iAmRoot;
+
+    bool saveLogBarrierIterates = false;
+
+    int linSolver=0;
+    double linSolveAbsTol = 1e-12;
+    double linSolveRelTol = 1e-6;
+    int relax_type = 8;
+    bool nocontact = false;
+public:
+    ParInteriorPointSolver(QPOptParContactProblem*);
+    double MaxStepSize(Vector& , Vector& , Vector& , double);
+    double MaxStepSize(Vector& , Vector& , double);
+    void Mult(const BlockVector& , BlockVector&);
+    void Mult(const Vector&, Vector &); 
+    void FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
+    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , double, bool);
+    void lineSearch(BlockVector& , BlockVector& , double);
+    void projectZ(const Vector & , Vector &, double);
+    void filterCheck(double, double);
+    double E(const BlockVector &, const Vector &, const Vector &, double, bool);
+    double E(const BlockVector &, const Vector &, const Vector &, bool);
+    bool GetConverged() const;
+    Array<int> & GetCGIterNumbers() {return cgnum_iterations;}
+    Array<int> & GetCGNoContactIterNumbers() {return cgnum_iterations_nocontact;}
+    int GetNumIterations() {return iter;}
+    // TO DO: include Hessian of Lagrangian
+    double theta(const BlockVector &);
+    double phi(const BlockVector &, double);
+    void Dxphi(const BlockVector &, double, BlockVector &);
+    double L(const BlockVector &, const Vector &, const Vector &);
+    void DxL(const BlockVector &, const Vector &, const Vector &, BlockVector &);
+    void SetTol(double);
+    void SetMaxIter(int);
+    void SetBarrierParameter(double);    
+    void SaveLogBarrierHessianIterates(bool);
+    void SaveLambda(int);
+    void SaveZl(int);
+    void SetLinearSolver(int);
+    void SetLinearSolveAbsTol(double);
+    void SetLinearSolveRelTol(double);
+    void SetLinearSolveRelaxType(int);
+    
+    void SetElasticityOptions(ParFiniteElementSpace * pfes_)
+    {
+        pfes = pfes_;
+    }
+    void EnableNoContactSolve()
+    {
+        nocontact = true;
+    }
+    virtual ~ParInteriorPointSolver();
+};
+
+#endif

miniapps/contact/makefile

@@ -0,0 +1,109 @@
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/contact/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+# Include defaults.mk to get XLINKER
+#DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+#include $(DEFAULTS_MK)
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+FRAMEWORK_SRC = ipsolver/ParIPsolver.cpp problems/parproblems.cpp problems/parproblems_util.cpp
+CONTACT_SRC = contact.cpp $(FRAMEWORK_SRC)
+CONTACT_OBJ = $(CONTACT_SRC:.cpp=.o)
+
+CONTACT_FDCHECK_SRC = contactFDcheck.cpp $(FRAMEWORK_SRC)
+CONTACT_FDCHECK_OBJ = $(CONTACT_FDCHECK_SRC:.cpp=.o)
+
+SCRATCH_SRC = scratch.cpp $(FRAMEWORK_SRC)
+SCRATCH_OBJ = $(SCRATCH_SRC:.cpp=.o)
+
+
+SEQ_MINIAPPS = 
+PAR_MINIAPPS = scratch contact contactFDcheck
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
+
+# If MFEM_SHARED is set, add the ../common rpath
+COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all lib-common clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+%.o: %.cpp
+
+%.o: $(SRC)%.cpp $(wildcard $(SRC)%.hpp) $(MFEM_LIB_FILE)\
+ $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+problems/%.o: $(SRC)problems/%.cpp $(wildcard $(SRC)problems/%.hpp) $(MFEM_LIB_FILE)\
+ $(CONFIG_MK) | lib-common
+	mkdir -p $(@D)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+
+all: $(MINIAPPS)
+
+contact: $(CONTACT_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(CONTACT_OBJ) $(COMMON_LIB) $(MFEM_LIBS) \
+		-l$(patsubst lib%,%,$(basename $(notdir $(MFEM_LIB_FILE))))
+
+contactFDcheck: $(CONTACT_FDCHECK_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(CONTACT_FDCHECK_OBJ) $(COMMON_LIB) $(MFEM_LIBS) \
+		-l$(patsubst lib%,%,$(basename $(notdir $(MFEM_LIB_FILE))))
+
+scratch: $(SCRATCH_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(SCRATCH_OBJ) $(COMMON_LIB) $(MFEM_LIBS) \
+		-l$(patsubst lib%,%,$(basename $(notdir $(MFEM_LIB_FILE))))
+
+
+# Rule for building lib-common
+lib-common:
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Specific execution options
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+contact-test-par: contact
+	@$(call mfem-test,$<, $(RUN_MPI), pcontact miniapp,)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+	rm -f $(CONTACT_OBJ)  
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf ParaView 

miniapps/contact/meshes/Test41.mesh

@@ -0,0 +1,453 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+89
+1 5 0 1 5 4 40 41 45 44
+1 5 40 41 45 44 80 81 85 84
+1 5 44 45 49 48 84 85 89 88
+1 5 4 5 9 8 44 45 49 48
+1 5 5 6 10 9 45 46 50 49
+1 5 45 46 50 49 85 86 90 89
+1 5 41 42 46 45 81 82 86 85
+1 5 1 2 6 5 41 42 46 45
+1 5 2 3 7 6 42 43 47 46
+1 5 42 43 47 46 82 83 87 86
+1 5 6 7 11 10 46 47 51 50
+1 5 46 47 51 50 86 87 91 90
+1 5 86 87 91 90 126 127 131 130
+1 5 82 83 87 86 122 123 127 126
+1 5 81 82 86 85 121 122 126 125
+1 5 80 81 85 84 120 121 125 124
+1 5 84 85 89 88 124 125 129 128
+1 5 85 86 90 89 125 126 130 129
+1 5 89 90 94 93 129 130 134 133
+1 5 88 89 93 92 128 129 133 132
+1 5 92 93 97 96 132 133 137 136
+1 5 93 94 98 97 133 134 138 137
+1 5 94 95 99 98 134 135 139 138
+1 5 54 55 59 58 94 95 99 98
+1 5 90 91 95 94 130 131 135 134
+1 5 50 51 55 54 90 91 95 94
+1 5 10 11 15 14 50 51 55 54
+1 5 14 15 19 18 54 55 59 58
+1 5 13 14 18 17 53 54 58 57
+1 5 53 54 58 57 93 94 98 97
+1 5 49 50 54 53 89 90 94 93
+1 5 9 10 14 13 49 50 54 53
+1 5 8 9 13 12 48 49 53 52
+1 5 48 49 53 52 88 89 93 92
+1 5 52 53 57 56 92 93 97 96
+1 5 12 13 17 16 52 53 57 56
+1 5 16 17 21 20 56 57 61 60
+1 5 56 57 61 60 96 97 101 100
+1 5 57 58 62 61 97 98 102 101
+1 5 17 18 22 21 57 58 62 61
+1 5 18 19 23 22 58 59 63 62
+1 5 58 59 63 62 98 99 103 102
+1 5 98 99 103 102 138 139 143 142
+1 5 97 98 102 101 137 138 142 141
+1 5 96 97 101 100 136 137 141 140
+1 5 100 101 105 104 140 141 145 144
+1 5 101 102 106 105 141 142 146 145
+1 5 102 103 107 106 142 143 147 146
+1 5 62 63 67 66 102 103 107 106
+1 5 22 23 27 26 62 63 67 66
+1 5 21 22 26 25 61 62 66 65
+1 5 61 62 66 65 101 102 106 105
+1 5 60 61 65 64 100 101 105 104
+1 5 20 21 25 24 60 61 65 64
+1 5 24 25 29 28 64 65 69 68
+1 5 64 65 69 68 104 105 109 108
+1 5 68 69 73 72 108 109 113 112
+1 5 28 29 33 32 68 69 73 72
+1 5 29 30 34 33 69 70 74 73
+1 5 69 70 74 73 109 110 114 113
+1 5 65 66 70 69 105 106 110 109
+1 5 25 26 30 29 65 66 70 69
+1 5 26 27 31 30 66 67 71 70
+1 5 66 67 71 70 106 107 111 110
+1 5 30 31 35 34 70 71 75 74
+1 5 70 71 75 74 110 111 115 114
+1 5 110 111 115 114 150 151 155 154
+1 5 106 107 111 110 146 147 151 150
+1 5 105 106 110 109 145 146 150 149
+1 5 109 110 114 113 149 150 154 153
+1 5 104 105 109 108 144 145 149 148
+1 5 108 109 113 112 148 149 153 152
+1 5 112 113 117 116 152 153 157 156
+1 5 113 114 118 117 153 154 158 157
+1 5 114 115 119 118 154 155 159 158
+1 5 74 75 79 78 114 115 119 118
+1 5 34 35 39 38 74 75 79 78
+1 5 33 34 38 37 73 74 78 77
+1 5 73 74 78 77 113 114 118 117
+1 5 72 73 77 76 112 113 117 116
+1 5 32 33 37 36 72 73 77 76
+2 5 160 161 164 163 169 170 173 172
+2 5 163 164 167 166 172 173 176 175
+2 5 172 173 176 175 181 182 185 184
+2 5 169 170 173 172 178 179 182 181
+2 5 170 171 174 173 179 180 183 182
+2 5 173 174 177 176 182 183 186 185
+2 5 164 165 168 167 173 174 177 176
+2 5 161 162 165 164 170 171 174 173
+
+boundary
+150
+1 3 0 4 5 1
+1 3 1 5 6 2
+1 3 2 6 7 3
+1 3 4 8 9 5
+1 3 5 9 10 6
+1 3 6 10 11 7
+1 3 8 12 13 9
+1 3 9 13 14 10
+1 3 10 14 15 11
+1 3 12 16 17 13
+1 3 13 17 18 14
+1 3 14 18 19 15
+1 3 16 20 21 17
+1 3 17 21 22 18
+1 3 18 22 23 19
+1 3 20 24 25 21
+1 3 21 25 26 22
+1 3 22 26 27 23
+1 3 24 28 29 25
+1 3 25 29 30 26
+1 3 26 30 31 27
+1 3 28 32 33 29
+1 3 29 33 34 30
+1 3 30 34 35 31
+1 3 32 36 37 33
+1 3 33 37 38 34
+1 3 34 38 39 35
+1 3 120 121 125 124
+1 3 121 122 126 125
+1 3 122 123 127 126
+1 3 124 125 129 128
+1 3 125 126 130 129
+1 3 126 127 131 130
+1 3 128 129 133 132
+1 3 129 130 134 133
+1 3 130 131 135 134
+1 3 132 133 137 136
+1 3 133 134 138 137
+1 3 134 135 139 138
+1 3 136 137 141 140
+1 3 137 138 142 141
+1 3 138 139 143 142
+1 3 140 141 145 144
+1 3 141 142 146 145
+1 3 142 143 147 146
+1 3 144 145 149 148
+1 3 145 146 150 149
+1 3 146 147 151 150
+1 3 148 149 153 152
+1 3 149 150 154 153
+1 3 150 151 155 154
+1 3 152 153 157 156
+1 3 153 154 158 157
+1 3 154 155 159 158
+2 3 0 40 44 4
+2 3 4 44 48 8
+2 3 8 48 52 12
+2 3 12 52 56 16
+2 3 16 56 60 20
+2 3 20 60 64 24
+2 3 24 64 68 28
+2 3 28 68 72 32
+2 3 32 72 76 36
+2 3 40 80 84 44
+2 3 44 84 88 48
+2 3 48 88 92 52
+2 3 52 92 96 56
+2 3 56 96 100 60
+2 3 60 100 104 64
+2 3 64 104 108 68
+2 3 68 108 112 72
+2 3 72 112 116 76
+2 3 80 120 124 84
+2 3 84 124 128 88
+2 3 88 128 132 92
+2 3 92 132 136 96
+2 3 96 136 140 100
+2 3 100 140 144 104
+2 3 104 144 148 108
+2 3 108 148 152 112
+2 3 112 152 156 116
+3 3 3 7 47 43
+3 3 7 11 51 47
+3 3 11 15 55 51
+3 3 15 19 59 55
+3 3 19 23 63 59
+3 3 23 27 67 63
+3 3 27 31 71 67
+3 3 31 35 75 71
+3 3 35 39 79 75
+3 3 43 47 87 83
+3 3 47 51 91 87
+3 3 51 55 95 91
+3 3 55 59 99 95
+3 3 59 63 103 99
+3 3 63 67 107 103
+3 3 67 71 111 107
+3 3 71 75 115 111
+3 3 75 79 119 115
+3 3 83 87 127 123
+3 3 87 91 131 127
+3 3 91 95 135 131
+3 3 95 99 139 135
+3 3 99 103 143 139
+3 3 103 107 147 143
+3 3 107 111 151 147
+3 3 111 115 155 151
+3 3 115 119 159 155
+1 3 0 1 41 40
+1 3 40 41 81 80
+1 3 80 81 121 120
+1 3 1 2 42 41
+1 3 41 42 82 81
+1 3 81 82 122 121
+1 3 2 3 43 42
+1 3 42 43 83 82
+1 3 82 83 123 122
+1 3 36 76 77 37
+1 3 76 116 117 77
+1 3 116 156 157 117
+1 3 37 77 78 38
+1 3 77 117 118 78
+1 3 117 157 158 118
+1 3 38 78 79 39
+1 3 78 118 119 79
+1 3 118 158 159 119
+5 3 160 163 164 161
+5 3 161 164 165 162
+5 3 163 166 167 164
+5 3 164 167 168 165
+5 3 178 179 182 181
+5 3 179 180 183 182
+5 3 181 182 185 184
+5 3 182 183 186 185
+4 3 160 169 172 163
+4 3 163 172 175 166
+4 3 169 178 181 172
+4 3 172 181 184 175
+6 3 162 165 174 171
+6 3 165 168 177 174
+6 3 171 174 183 180
+6 3 174 177 186 183
+5 3 160 161 170 169
+5 3 169 170 179 178
+5 3 161 162 171 170
+5 3 170 171 180 179
+5 3 166 175 176 167
+5 3 175 184 185 176
+5 3 167 176 177 168
+5 3 176 185 186 177
+
+vertices
+187
+3
+-1 0 0
+-0.66666667 0 0
+-0.33333333 0 0
+0 0 0
+-1 0.33333333 0
+-0.66666667 0.33333333 0
+-0.33333333 0.33333333 0
+0 0.33333333 0
+-1 0.66666667 0
+-0.66666667 0.66666667 0
+-0.33333333 0.66666667 0
+0 0.66666667 0
+-1 1 0
+-0.66666667 1 0
+-0.33333333 1 0
+0 1 0
+-1 1.3333333 0
+-0.66666667 1.3333333 0
+-0.33333333 1.3333333 0
+0 1.3333333 0
+-1 1.6666667 0
+-0.66666667 1.6666667 0
+-0.33333333 1.6666667 0
+0 1.6666667 0
+-1 2 0
+-0.66666667 2 0
+-0.33333333 2 0
+0 2 0
+-1 2.3333333 0
+-0.66666667 2.3333333 0
+-0.33333333 2.3333333 0
+0 2.3333333 0
+-1 2.6666667 0
+-0.66666667 2.6666667 0
+-0.33333333 2.6666667 0
+0 2.6666667 0
+-1 3 0
+-0.66666667 3 0
+-0.33333333 3 0
+0 3 0
+-1 0 0.33333333
+-0.66666667 0 0.33333333
+-0.33333333 0 0.33333333
+0 0 0.33333333
+-1 0.33333333 0.33333333
+-0.66666667 0.33333333 0.33333333
+-0.33333333 0.33333333 0.33333333
+0 0.33333333 0.33333333
+-1 0.66666667 0.33333333
+-0.66666667 0.66666667 0.33333333
+-0.33333333 0.66666667 0.33333333
+0 0.66666667 0.33333333
+-1 1 0.33333333
+-0.66666667 1 0.33333333
+-0.33333333 1 0.33333333
+0 1 0.33333333
+-1 1.3333333 0.33333333
+-0.66666667 1.3333333 0.33333333
+-0.33333333 1.3333333 0.33333333
+0 1.3333333 0.33333333
+-1 1.6666667 0.33333333
+-0.66666667 1.6666667 0.33333333
+-0.33333333 1.6666667 0.33333333
+0 1.6666667 0.33333333
+-1 2 0.33333333
+-0.66666667 2 0.33333333
+-0.33333333 2 0.33333333
+0 2 0.33333333
+-1 2.3333333 0.33333333
+-0.66666667 2.3333333 0.33333333
+-0.33333333 2.3333333 0.33333333
+0 2.3333333 0.33333333
+-1 2.6666667 0.33333333
+-0.66666667 2.6666667 0.33333333
+-0.33333333 2.6666667 0.33333333
+0 2.6666667 0.33333333
+-1 3 0.33333333
+-0.66666667 3 0.33333333
+-0.33333333 3 0.33333333
+0 3 0.33333333
+-1 0 0.66666667
+-0.66666667 0 0.66666667
+-0.33333333 0 0.66666667
+0 0 0.66666667
+-1 0.33333333 0.66666667
+-0.66666667 0.33333333 0.66666667
+-0.33333333 0.33333333 0.66666667
+0 0.33333333 0.66666667
+-1 0.66666667 0.66666667
+-0.66666667 0.66666667 0.66666667
+-0.33333333 0.66666667 0.66666667
+0 0.66666667 0.66666667
+-1 1 0.66666667
+-0.66666667 1 0.66666667
+-0.33333333 1 0.66666667
+0 1 0.66666667
+-1 1.3333333 0.66666667
+-0.66666667 1.3333333 0.66666667
+-0.33333333 1.3333333 0.66666667
+0 1.3333333 0.66666667
+-1 1.6666667 0.66666667
+-0.66666667 1.6666667 0.66666667
+-0.33333333 1.6666667 0.66666667
+0 1.6666667 0.66666667
+-1 2 0.66666667
+-0.66666667 2 0.66666667
+-0.33333333 2 0.66666667
+0 2 0.66666667
+-1 2.3333333 0.66666667
+-0.66666667 2.3333333 0.66666667
+-0.33333333 2.3333333 0.66666667
+0 2.3333333 0.66666667
+-1 2.6666667 0.66666667
+-0.66666667 2.6666667 0.66666667
+-0.33333333 2.6666667 0.66666667
+0 2.6666667 0.66666667
+-1 3 0.66666667
+-0.66666667 3 0.66666667
+-0.33333333 3 0.66666667
+0 3 0.66666667
+-1 0 1
+-0.66666667 0 1
+-0.33333333 0 1
+0 0 1
+-1 0.33333333 1
+-0.66666667 0.33333333 1
+-0.33333333 0.33333333 1
+0 0.33333333 1
+-1 0.66666667 1
+-0.66666667 0.66666667 1
+-0.33333333 0.66666667 1
+0 0.66666667 1
+-1 1 1
+-0.66666667 1 1
+-0.33333333 1 1
+0 1 1
+-1 1.3333333 1
+-0.66666667 1.3333333 1
+-0.33333333 1.3333333 1
+0 1.3333333 1
+-1 1.6666667 1
+-0.66666667 1.6666667 1
+-0.33333333 1.6666667 1
+0 1.6666667 1
+-1 2 1
+-0.66666667 2 1
+-0.33333333 2 1
+0 2 1
+-1 2.3333333 1
+-0.66666667 2.3333333 1
+-0.33333333 2.3333333 1
+0 2.3333333 1
+-1 2.6666667 1
+-0.66666667 2.6666667 1
+-0.33333333 2.6666667 1
+0 2.6666667 1
+-1 3 1
+-0.66666667 3 1
+-0.33333333 3 1
+0 3 1
+0 1.5 0.25251263
+0.175 1.5 0.25251263
+0.35 1.5 0.25251263
+0 1.6237437 0.37625631
+0.175 1.6237437 0.37625631
+0.35 1.6237437 0.37625631
+0 1.7474874 0.5
+0.175 1.7474874 0.5
+0.35 1.7474874 0.5
+0 1.3762563 0.37625631
+0.175 1.3762563 0.37625631
+0.35 1.3762563 0.37625631
+0 1.5 0.5
+0.175 1.5 0.5
+0.35 1.5 0.5
+0 1.6237437 0.62374369
+0.175 1.6237437 0.62374369
+0.35 1.6237437 0.62374369
+0 1.2525126 0.5
+0.175 1.2525126 0.5
+0.35 1.2525126 0.5
+0 1.3762563 0.62374369
+0.175 1.3762563 0.62374369
+0.35 1.3762563 0.62374369
+0 1.5 0.74748737
+0.175 1.5 0.74748737
+0.35 1.5 0.74748737

miniapps/contact/meshes/Test42.mesh

@@ -0,0 +1,453 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+89
+1 5 0 1 5 4 40 41 45 44
+1 5 40 41 45 44 80 81 85 84
+1 5 44 45 49 48 84 85 89 88
+1 5 4 5 9 8 44 45 49 48
+1 5 5 6 10 9 45 46 50 49
+1 5 45 46 50 49 85 86 90 89
+1 5 41 42 46 45 81 82 86 85
+1 5 1 2 6 5 41 42 46 45
+1 5 2 3 7 6 42 43 47 46
+1 5 42 43 47 46 82 83 87 86
+1 5 6 7 11 10 46 47 51 50
+1 5 46 47 51 50 86 87 91 90
+1 5 86 87 91 90 126 127 131 130
+1 5 82 83 87 86 122 123 127 126
+1 5 81 82 86 85 121 122 126 125
+1 5 80 81 85 84 120 121 125 124
+1 5 84 85 89 88 124 125 129 128
+1 5 85 86 90 89 125 126 130 129
+1 5 89 90 94 93 129 130 134 133
+1 5 88 89 93 92 128 129 133 132
+1 5 92 93 97 96 132 133 137 136
+1 5 93 94 98 97 133 134 138 137
+1 5 94 95 99 98 134 135 139 138
+1 5 54 55 59 58 94 95 99 98
+1 5 90 91 95 94 130 131 135 134
+1 5 50 51 55 54 90 91 95 94
+1 5 10 11 15 14 50 51 55 54
+1 5 14 15 19 18 54 55 59 58
+1 5 13 14 18 17 53 54 58 57
+1 5 53 54 58 57 93 94 98 97
+1 5 49 50 54 53 89 90 94 93
+1 5 9 10 14 13 49 50 54 53
+1 5 8 9 13 12 48 49 53 52
+1 5 48 49 53 52 88 89 93 92
+1 5 52 53 57 56 92 93 97 96
+1 5 12 13 17 16 52 53 57 56
+1 5 16 17 21 20 56 57 61 60
+1 5 56 57 61 60 96 97 101 100
+1 5 57 58 62 61 97 98 102 101
+1 5 17 18 22 21 57 58 62 61
+1 5 18 19 23 22 58 59 63 62
+1 5 58 59 63 62 98 99 103 102
+1 5 98 99 103 102 138 139 143 142
+1 5 97 98 102 101 137 138 142 141
+1 5 96 97 101 100 136 137 141 140
+1 5 100 101 105 104 140 141 145 144
+1 5 101 102 106 105 141 142 146 145
+1 5 102 103 107 106 142 143 147 146
+1 5 62 63 67 66 102 103 107 106
+1 5 22 23 27 26 62 63 67 66
+1 5 21 22 26 25 61 62 66 65
+1 5 61 62 66 65 101 102 106 105
+1 5 60 61 65 64 100 101 105 104
+1 5 20 21 25 24 60 61 65 64
+1 5 24 25 29 28 64 65 69 68
+1 5 64 65 69 68 104 105 109 108
+1 5 68 69 73 72 108 109 113 112
+1 5 28 29 33 32 68 69 73 72
+1 5 29 30 34 33 69 70 74 73
+1 5 69 70 74 73 109 110 114 113
+1 5 65 66 70 69 105 106 110 109
+1 5 25 26 30 29 65 66 70 69
+1 5 26 27 31 30 66 67 71 70
+1 5 66 67 71 70 106 107 111 110
+1 5 30 31 35 34 70 71 75 74
+1 5 70 71 75 74 110 111 115 114
+1 5 110 111 115 114 150 151 155 154
+1 5 106 107 111 110 146 147 151 150
+1 5 105 106 110 109 145 146 150 149
+1 5 109 110 114 113 149 150 154 153
+1 5 104 105 109 108 144 145 149 148
+1 5 108 109 113 112 148 149 153 152
+1 5 112 113 117 116 152 153 157 156
+1 5 113 114 118 117 153 154 158 157
+1 5 114 115 119 118 154 155 159 158
+1 5 74 75 79 78 114 115 119 118
+1 5 34 35 39 38 74 75 79 78
+1 5 33 34 38 37 73 74 78 77
+1 5 73 74 78 77 113 114 118 117
+1 5 72 73 77 76 112 113 117 116
+1 5 32 33 37 36 72 73 77 76
+2 5 160 161 164 163 169 170 173 172
+2 5 163 164 167 166 172 173 176 175
+2 5 172 173 176 175 181 182 185 184
+2 5 169 170 173 172 178 179 182 181
+2 5 170 171 174 173 179 180 183 182
+2 5 173 174 177 176 182 183 186 185
+2 5 164 165 168 167 173 174 177 176
+2 5 161 162 165 164 170 171 174 173
+
+boundary
+150
+1 3 0 4 5 1
+1 3 1 5 6 2
+1 3 2 6 7 3
+1 3 4 8 9 5
+1 3 5 9 10 6
+1 3 6 10 11 7
+1 3 8 12 13 9
+1 3 9 13 14 10
+1 3 10 14 15 11
+1 3 12 16 17 13
+1 3 13 17 18 14
+1 3 14 18 19 15
+1 3 16 20 21 17
+1 3 17 21 22 18
+1 3 18 22 23 19
+1 3 20 24 25 21
+1 3 21 25 26 22
+1 3 22 26 27 23
+1 3 24 28 29 25
+1 3 25 29 30 26
+1 3 26 30 31 27
+1 3 28 32 33 29
+1 3 29 33 34 30
+1 3 30 34 35 31
+1 3 32 36 37 33
+1 3 33 37 38 34
+1 3 34 38 39 35
+1 3 120 121 125 124
+1 3 121 122 126 125
+1 3 122 123 127 126
+1 3 124 125 129 128
+1 3 125 126 130 129
+1 3 126 127 131 130
+1 3 128 129 133 132
+1 3 129 130 134 133
+1 3 130 131 135 134
+1 3 132 133 137 136
+1 3 133 134 138 137
+1 3 134 135 139 138
+1 3 136 137 141 140
+1 3 137 138 142 141
+1 3 138 139 143 142
+1 3 140 141 145 144
+1 3 141 142 146 145
+1 3 142 143 147 146
+1 3 144 145 149 148
+1 3 145 146 150 149
+1 3 146 147 151 150
+1 3 148 149 153 152
+1 3 149 150 154 153
+1 3 150 151 155 154
+1 3 152 153 157 156
+1 3 153 154 158 157
+1 3 154 155 159 158
+2 3 0 40 44 4
+2 3 4 44 48 8
+2 3 8 48 52 12
+2 3 12 52 56 16
+2 3 16 56 60 20
+2 3 20 60 64 24
+2 3 24 64 68 28
+2 3 28 68 72 32
+2 3 32 72 76 36
+2 3 40 80 84 44
+2 3 44 84 88 48
+2 3 48 88 92 52
+2 3 52 92 96 56
+2 3 56 96 100 60
+2 3 60 100 104 64
+2 3 64 104 108 68
+2 3 68 108 112 72
+2 3 72 112 116 76
+2 3 80 120 124 84
+2 3 84 124 128 88
+2 3 88 128 132 92
+2 3 92 132 136 96
+2 3 96 136 140 100
+2 3 100 140 144 104
+2 3 104 144 148 108
+2 3 108 148 152 112
+2 3 112 152 156 116
+3 3 3 7 47 43
+3 3 7 11 51 47
+3 3 11 15 55 51
+3 3 15 19 59 55
+3 3 19 23 63 59
+3 3 23 27 67 63
+3 3 27 31 71 67
+3 3 31 35 75 71
+3 3 35 39 79 75
+3 3 43 47 87 83
+3 3 47 51 91 87
+3 3 51 55 95 91
+3 3 55 59 99 95
+3 3 59 63 103 99
+3 3 63 67 107 103
+3 3 67 71 111 107
+3 3 71 75 115 111
+3 3 75 79 119 115
+3 3 83 87 127 123
+3 3 87 91 131 127
+3 3 91 95 135 131
+3 3 95 99 139 135
+3 3 99 103 143 139
+3 3 103 107 147 143
+3 3 107 111 151 147
+3 3 111 115 155 151
+3 3 115 119 159 155
+1 3 0 1 41 40
+1 3 40 41 81 80
+1 3 80 81 121 120
+1 3 1 2 42 41
+1 3 41 42 82 81
+1 3 81 82 122 121
+1 3 2 3 43 42
+1 3 42 43 83 82
+1 3 82 83 123 122
+1 3 36 76 77 37
+1 3 76 116 117 77
+1 3 116 156 157 117
+1 3 37 77 78 38
+1 3 77 117 118 78
+1 3 117 157 158 118
+1 3 38 78 79 39
+1 3 78 118 119 79
+1 3 118 158 159 119
+5 3 160 163 164 161
+5 3 161 164 165 162
+5 3 163 166 167 164
+5 3 164 167 168 165
+5 3 178 179 182 181
+5 3 179 180 183 182
+5 3 181 182 185 184
+5 3 182 183 186 185
+4 3 160 169 172 163
+4 3 163 172 175 166
+4 3 169 178 181 172
+4 3 172 181 184 175
+6 3 162 165 174 171
+6 3 165 168 177 174
+6 3 171 174 183 180
+6 3 174 177 186 183
+5 3 160 161 170 169
+5 3 169 170 179 178
+5 3 161 162 171 170
+5 3 170 171 180 179
+5 3 166 175 176 167
+5 3 175 184 185 176
+5 3 167 176 177 168
+5 3 176 185 186 177
+
+vertices
+187
+3
+-1 0 0
+-0.66666667 0 0
+-0.33333333 0 0
+0 0 0
+-1 0.33333333 0
+-0.66666667 0.33333333 0
+-0.33333333 0.33333333 0
+0 0.33333333 0
+-1 0.66666667 0
+-0.66666667 0.66666667 0
+-0.33333333 0.66666667 0
+0 0.66666667 0
+-1 1 0
+-0.66666667 1 0
+-0.33333333 1 0
+0 1 0
+-1 1.3333333 0
+-0.66666667 1.3333333 0
+-0.33333333 1.3333333 0
+0 1.3333333 0
+-1 1.6666667 0
+-0.66666667 1.6666667 0
+-0.33333333 1.6666667 0
+0 1.6666667 0
+-1 2 0
+-0.66666667 2 0
+-0.33333333 2 0
+0 2 0
+-1 2.3333333 0
+-0.66666667 2.3333333 0
+-0.33333333 2.3333333 0
+0 2.3333333 0
+-1 2.6666667 0
+-0.66666667 2.6666667 0
+-0.33333333 2.6666667 0
+0 2.6666667 0
+-1 3 0
+-0.66666667 3 0
+-0.33333333 3 0
+0 3 0
+-1 0 0.33333333
+-0.66666667 0 0.33333333
+-0.33333333 0 0.33333333
+0 0 0.33333333
+-1 0.33333333 0.33333333
+-0.66666667 0.33333333 0.33333333
+-0.33333333 0.33333333 0.33333333
+0 0.33333333 0.33333333
+-1 0.66666667 0.33333333
+-0.66666667 0.66666667 0.33333333
+-0.33333333 0.66666667 0.33333333
+0 0.66666667 0.33333333
+-1 1 0.33333333
+-0.66666667 1 0.33333333
+-0.33333333 1 0.33333333
+0 1 0.33333333
+-1 1.3333333 0.33333333
+-0.66666667 1.3333333 0.33333333
+-0.33333333 1.3333333 0.33333333
+0 1.3333333 0.33333333
+-1 1.6666667 0.33333333
+-0.66666667 1.6666667 0.33333333
+-0.33333333 1.6666667 0.33333333
+0 1.6666667 0.33333333
+-1 2 0.33333333
+-0.66666667 2 0.33333333
+-0.33333333 2 0.33333333
+0 2 0.33333333
+-1 2.3333333 0.33333333
+-0.66666667 2.3333333 0.33333333
+-0.33333333 2.3333333 0.33333333
+0 2.3333333 0.33333333
+-1 2.6666667 0.33333333
+-0.66666667 2.6666667 0.33333333
+-0.33333333 2.6666667 0.33333333
+0 2.6666667 0.33333333
+-1 3 0.33333333
+-0.66666667 3 0.33333333
+-0.33333333 3 0.33333333
+0 3 0.33333333
+-1 0 0.66666667
+-0.66666667 0 0.66666667
+-0.33333333 0 0.66666667
+0 0 0.66666667
+-1 0.33333333 0.66666667
+-0.66666667 0.33333333 0.66666667
+-0.33333333 0.33333333 0.66666667
+0 0.33333333 0.66666667
+-1 0.66666667 0.66666667
+-0.66666667 0.66666667 0.66666667
+-0.33333333 0.66666667 0.66666667
+0 0.66666667 0.66666667
+-1 1 0.66666667
+-0.66666667 1 0.66666667
+-0.33333333 1 0.66666667
+0 1 0.66666667
+-1 1.3333333 0.66666667
+-0.66666667 1.3333333 0.66666667
+-0.33333333 1.3333333 0.66666667
+0 1.3333333 0.66666667
+-1 1.6666667 0.66666667
+-0.66666667 1.6666667 0.66666667
+-0.33333333 1.6666667 0.66666667
+0 1.6666667 0.66666667
+-1 2 0.66666667
+-0.66666667 2 0.66666667
+-0.33333333 2 0.66666667
+0 2 0.66666667
+-1 2.3333333 0.66666667
+-0.66666667 2.3333333 0.66666667
+-0.33333333 2.3333333 0.66666667
+0 2.3333333 0.66666667
+-1 2.6666667 0.66666667
+-0.66666667 2.6666667 0.66666667
+-0.33333333 2.6666667 0.66666667
+0 2.6666667 0.66666667
+-1 3 0.66666667
+-0.66666667 3 0.66666667
+-0.33333333 3 0.66666667
+0 3 0.66666667
+-1 0 1
+-0.66666667 0 1
+-0.33333333 0 1
+0 0 1
+-1 0.33333333 1
+-0.66666667 0.33333333 1
+-0.33333333 0.33333333 1
+0 0.33333333 1
+-1 0.66666667 1
+-0.66666667 0.66666667 1
+-0.33333333 0.66666667 1
+0 0.66666667 1
+-1 1 1
+-0.66666667 1 1
+-0.33333333 1 1
+0 1 1
+-1 1.3333333 1
+-0.66666667 1.3333333 1
+-0.33333333 1.3333333 1
+0 1.3333333 1
+-1 1.6666667 1
+-0.66666667 1.6666667 1
+-0.33333333 1.6666667 1
+0 1.6666667 1
+-1 2 1
+-0.66666667 2 1
+-0.33333333 2 1
+0 2 1
+-1 2.3333333 1
+-0.66666667 2.3333333 1
+-0.33333333 2.3333333 1
+0 2.3333333 1
+-1 2.6666667 1
+-0.66666667 2.6666667 1
+-0.33333333 2.6666667 1
+0 2.6666667 1
+-1 3 1
+-0.66666667 3 1
+-0.33333333 3 1
+0 3 1
+0 0.83333333 0.25251263
+0.175 0.83333333 0.25251263
+0.35 0.83333333 0.25251263
+0 0.95707702 0.37625631
+0.175 0.95707702 0.37625631
+0.35 0.95707702 0.37625631
+0 1.0808207 0.5
+0.175 1.0808207 0.5
+0.35 1.0808207 0.5
+0 0.70958965 0.37625631
+0.175 0.70958965 0.37625631
+0.35 0.70958965 0.37625631
+0 0.83333333 0.5
+0.175 0.83333333 0.5
+0.35 0.83333333 0.5
+0 0.95707702 0.62374369
+0.175 0.95707702 0.62374369
+0.35 0.95707702 0.62374369
+0 0.58584596 0.5
+0.175 0.58584596 0.5
+0.35 0.58584596 0.5
+0 0.70958965 0.62374369
+0.175 0.70958965 0.62374369
+0.35 0.70958965 0.62374369
+0 0.83333333 0.74748737
+0.175 0.83333333 0.74748737
+0.35 0.83333333 0.74748737

miniapps/contact/meshes/two-block.mesh

@@ -0,0 +1,231 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+35
+1 5 0 1 5 4 16 17 21 20
+1 5 16 17 21 20 32 33 37 36
+1 5 17 18 22 21 33 34 38 37
+1 5 1 2 6 5 17 18 22 21
+1 5 5 6 10 9 21 22 26 25
+1 5 21 22 26 25 37 38 42 41
+1 5 20 21 25 24 36 37 41 40
+1 5 4 5 9 8 20 21 25 24
+1 5 8 9 13 12 24 25 29 28
+1 5 24 25 29 28 40 41 45 44
+1 5 9 10 14 13 25 26 30 29
+1 5 25 26 30 29 41 42 46 45
+1 5 41 42 46 45 57 58 62 61
+1 5 40 41 45 44 56 57 61 60
+1 5 36 37 41 40 52 53 57 56
+1 5 37 38 42 41 53 54 58 57
+1 5 32 33 37 36 48 49 53 52
+1 5 33 34 38 37 49 50 54 53
+1 5 34 35 39 38 50 51 55 54
+1 5 38 39 43 42 54 55 59 58
+1 5 42 43 47 46 58 59 63 62
+1 5 26 27 31 30 42 43 47 46
+1 5 10 11 15 14 26 27 31 30
+1 5 6 7 11 10 22 23 27 26
+1 5 22 23 27 26 38 39 43 42
+1 5 18 19 23 22 34 35 39 38
+1 5 2 3 7 6 18 19 23 22
+1 5 64 65 68 67 73 74 77 76
+1 5 67 68 71 70 76 77 80 79
+1 5 76 77 80 79 85 86 89 88
+1 5 73 74 77 76 82 83 86 85
+1 5 74 75 78 77 83 84 87 86
+1 5 77 78 81 80 86 87 90 89
+1 5 68 69 72 71 77 78 81 80
+1 5 65 66 69 68 74 75 78 77
+
+boundary
+78
+1 3 0 4 5 1
+1 3 1 5 6 2
+1 3 2 6 7 3
+1 3 4 8 9 5
+1 3 5 9 10 6
+1 3 6 10 11 7
+1 3 8 12 13 9
+1 3 9 13 14 10
+1 3 10 14 15 11
+1 3 48 49 53 52
+1 3 49 50 54 53
+1 3 50 51 55 54
+1 3 52 53 57 56
+1 3 53 54 58 57
+1 3 54 55 59 58
+1 3 56 57 61 60
+1 3 57 58 62 61
+1 3 58 59 63 62
+2 3 0 16 20 4
+2 3 4 20 24 8
+2 3 8 24 28 12
+2 3 16 32 36 20
+2 3 20 36 40 24
+2 3 24 40 44 28
+2 3 32 48 52 36
+2 3 36 52 56 40
+2 3 40 56 60 44
+3 3 3 7 23 19
+3 3 7 11 27 23
+3 3 11 15 31 27
+3 3 19 23 39 35
+3 3 23 27 43 39
+3 3 27 31 47 43
+3 3 35 39 55 51
+3 3 39 43 59 55
+3 3 43 47 63 59
+1 3 0 1 17 16
+1 3 16 17 33 32
+1 3 32 33 49 48
+1 3 1 2 18 17
+1 3 17 18 34 33
+1 3 33 34 50 49
+1 3 2 3 19 18
+1 3 18 19 35 34
+1 3 34 35 51 50
+1 3 12 28 29 13
+1 3 28 44 45 29
+1 3 44 60 61 45
+1 3 13 29 30 14
+1 3 29 45 46 30
+1 3 45 61 62 46
+1 3 14 30 31 15
+1 3 30 46 47 31
+1 3 46 62 63 47
+5 3 64 67 68 65
+5 3 65 68 69 66
+5 3 67 70 71 68
+5 3 68 71 72 69
+5 3 82 83 86 85
+5 3 83 84 87 86
+5 3 85 86 89 88
+5 3 86 87 90 89
+4 3 64 73 76 67
+4 3 67 76 79 70
+4 3 73 82 85 76
+4 3 76 85 88 79
+6 3 66 69 78 75
+6 3 69 72 81 78
+6 3 75 78 87 84
+6 3 78 81 90 87
+5 3 64 65 74 73
+5 3 73 74 83 82
+5 3 65 66 75 74
+5 3 74 75 84 83
+5 3 70 79 80 71
+5 3 79 88 89 80
+5 3 71 80 81 72
+5 3 80 89 90 81
+
+vertices
+91
+3
+-1 0 0
+-0.66666667 0 0
+-0.33333333 0 0
+0 0 0
+-1 0.33333333 0
+-0.66666667 0.33333333 0
+-0.33333333 0.33333333 0
+0 0.33333333 0
+-1 0.66666667 0
+-0.66666667 0.66666667 0
+-0.33333333 0.66666667 0
+0 0.66666667 0
+-1 1 0
+-0.66666667 1 0
+-0.33333333 1 0
+0 1 0
+-1 0 0.33333333
+-0.66666667 0 0.33333333
+-0.33333333 0 0.33333333
+0 0 0.33333333
+-1 0.33333333 0.33333333
+-0.66666667 0.33333333 0.33333333
+-0.33333333 0.33333333 0.33333333
+0 0.33333333 0.33333333
+-1 0.66666667 0.33333333
+-0.66666667 0.66666667 0.33333333
+-0.33333333 0.66666667 0.33333333
+0 0.66666667 0.33333333
+-1 1 0.33333333
+-0.66666667 1 0.33333333
+-0.33333333 1 0.33333333
+0 1 0.33333333
+-1 0 0.66666667
+-0.66666667 0 0.66666667
+-0.33333333 0 0.66666667
+0 0 0.66666667
+-1 0.33333333 0.66666667
+-0.66666667 0.33333333 0.66666667
+-0.33333333 0.33333333 0.66666667
+0 0.33333333 0.66666667
+-1 0.66666667 0.66666667
+-0.66666667 0.66666667 0.66666667
+-0.33333333 0.66666667 0.66666667
+0 0.66666667 0.66666667
+-1 1 0.66666667
+-0.66666667 1 0.66666667
+-0.33333333 1 0.66666667
+0 1 0.66666667
+-1 0 1
+-0.66666667 0 1
+-0.33333333 0 1
+0 0 1
+-1 0.33333333 1
+-0.66666667 0.33333333 1
+-0.33333333 0.33333333 1
+0 0.33333333 1
+-1 0.66666667 1
+-0.66666667 0.66666667 1
+-0.33333333 0.66666667 1
+0 0.66666667 1
+-1 1 1
+-0.66666667 1 1
+-0.33333333 1 1
+0 1 1
+0 0.5 0.14644661
+0.25 0.5 0.14644661
+0.5 0.5 0.14644661
+0 0.6767767 0.3232233
+0.25 0.6767767 0.3232233
+0.5 0.6767767 0.3232233
+0 0.85355339 0.5
+0.25 0.85355339 0.5
+0.5 0.85355339 0.5
+0 0.3232233 0.3232233
+0.25 0.3232233 0.3232233
+0.5 0.3232233 0.3232233
+0 0.5 0.5
+0.25 0.5 0.5
+0.5 0.5 0.5
+0 0.6767767 0.6767767
+0.25 0.6767767 0.6767767
+0.5 0.6767767 0.6767767
+0 0.14644661 0.5
+0.25 0.14644661 0.5
+0.5 0.14644661 0.5
+0 0.3232233 0.6767767
+0.25 0.3232233 0.6767767
+0.5 0.3232233 0.6767767
+0 0.5 0.85355339
+0.25 0.5 0.85355339
+0.5 0.5 0.85355339

miniapps/contact/problems/parproblems.cpp

@@ -0,0 +1,663 @@
+#include "parproblems.hpp"
+
+void ParElasticityProblem::Init()
+{
+   int dim = pmesh->Dimension();
+   fec = new H1_FECollection(order,dim);
+   fes = new ParFiniteElementSpace(pmesh,fec,dim,Ordering::byVDIM);
+   ndofs = fes->GetVSize();
+   ntdofs = fes->GetTrueVSize();
+   gndofs = fes->GlobalTrueVSize();
+   pmesh->SetNodalFESpace(fes);
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+   }
+   ess_bdr = 0; 
+   Array<int> ess_tdof_list_temp;
+   for (int i = 0; i < ess_bdr_attr.Size(); i++ )
+   {
+      ess_bdr[ess_bdr_attr[i]-1] = 1;
+      fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list_temp,ess_bdr_attr_comp[i]);
+      ess_tdof_list.Append(ess_tdof_list_temp);
+      ess_bdr[ess_bdr_attr[i]-1] = 0;
+   }
+   // Solution GridFunction
+   x.SetSpace(fes);  x = 0.0;
+   // RHS
+   b = new ParLinearForm(fes);
+
+   // Elasticity operator
+   lambda.SetSize(pmesh->attributes.Max()); lambda = 57.6923076923;
+   mu.SetSize(pmesh->attributes.Max()); mu = 38.4615384615;
+
+   lambda_cf.UpdateConstants(lambda);
+   mu_cf.UpdateConstants(mu);
+
+   a = new ParBilinearForm(fes);
+   a->AddDomainIntegrator(new ElasticityIntegrator(lambda_cf,mu_cf));
+}
+
+void ParElasticityProblem::FormLinearSystem()
+{
+   if (!formsystem) 
+   {
+      formsystem = true;
+      b->Assemble();
+      a->Assemble();
+      a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+   }
+}
+
+void ParElasticityProblem::UpdateLinearSystem()
+{
+   UpdateStep();
+   FormLinearSystem();
+}
+
+// #ifdef MFEM_USE_TRIBOL
+
+
+
+ParContactProblem::ParContactProblem(ParElasticityProblem * prob_, 
+                                                         const std::set<int> & mortar_attrs_, 
+                                                         const std::set<int> & nonmortar_attrs_,
+                                                         ParGridFunction * coords_,
+                                                         bool doublepass_)
+: prob(prob_), mortar_attrs(mortar_attrs_), nonmortar_attrs(nonmortar_attrs_), doublepass(doublepass_), coords(coords_)
+{
+   ParMesh* pmesh = prob->GetMesh();
+   comm = pmesh->GetComm();
+   MPI_Comm_rank(comm, &myid);
+   MPI_Comm_size(comm, &numprocs);
+ 
+   dim = pmesh->Dimension();
+   nodes0.SetSpace(pmesh->GetNodes()->FESpace());
+   nodes0 = *pmesh->GetNodes();
+   nodes1 = pmesh->GetNodes();
+
+   prob->FormLinearSystem();
+   K = new HypreParMatrix(prob->GetOperator());
+   B = new Vector(prob->GetRHS());
+   if (doublepass)
+   {
+      SetupTribolDoublePass();
+   }
+   else
+   {
+      SetupTribol();
+   }
+}
+
+void ParContactProblem::SetupTribol()
+{
+   axom::slic::SimpleLogger logger;
+   axom::slic::setIsRoot(mfem::Mpi::Root());
+
+   // Initialize Tribol contact library
+   tribol::initialize(3, MPI_COMM_WORLD);
+
+   int coupling_scheme_id = 0;
+   int mesh1_id = 0;
+   int mesh2_id = 1;
+   vfes = prob->GetFESpace();
+   ParMesh * pmesh = prob->GetMesh();
+   tribol::registerMfemCouplingScheme(
+      coupling_scheme_id, mesh1_id, mesh2_id,
+      *pmesh, *coords, mortar_attrs, nonmortar_attrs,
+      tribol::SURFACE_TO_SURFACE,
+      tribol::NO_SLIDING,
+      tribol::SINGLE_MORTAR,
+      tribol::FRICTIONLESS,
+      tribol::LAGRANGE_MULTIPLIER,
+      tribol::BINNING_GRID
+   );
+
+   // Access Tribol's pressure grid function (on the contact surface)
+   auto& pressure = tribol::getMfemPressure(coupling_scheme_id);
+   if (mfem::Mpi::Root())
+   {
+      std::cout << "Number of pressure unknowns: " <<
+                pressure.ParFESpace()->GlobalTrueVSize() << std::endl;
+   }
+
+   // Set Tribol options for Lagrange multiplier enforcement
+   tribol::setLagrangeMultiplierOptions(
+      coupling_scheme_id,
+      tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
+   );
+
+   // Update contact mesh decomposition
+   tribol::updateMfemParallelDecomposition();
+
+   // Update contact gaps, forces, and tangent stiffness
+   int cycle = 1;   // pseudo cycle
+   double t = 1.0;  // pseudo time
+   double dt = 1.0; // pseudo dt
+   tribol::update(cycle, t, dt);
+
+   // Return contact contribution to the tangent stiffness matrix
+   auto A_blk = tribol::getMfemBlockJacobian(coupling_scheme_id);
+   
+   HypreParMatrix * Mfull = (HypreParMatrix *)(&A_blk->GetBlock(1,0));
+   Mfull->EliminateCols(prob->GetEssentialDofs());
+   int h = Mfull->Height();
+   SparseMatrix merged;
+   Mfull->MergeDiagAndOffd(merged);
+   Array<int> nonzero_rows;
+   for (int i = 0; i<h; i++)
+   {
+      if (!merged.RowIsEmpty(i))
+      {
+         nonzero_rows.Append(i);
+      }
+   }
+
+   int hnew = nonzero_rows.Size();
+   SparseMatrix P(hnew,h);
+
+   for (int i = 0; i<hnew; i++)
+   {
+      int col = nonzero_rows[i];
+      P.Set(i,col,1.0);
+   }
+   P.Finalize();
+
+   SparseMatrix * reduced_merged = Mult(P,merged);
+
+   int rows[2];
+   int cols[2];
+   cols[0] = Mfull->ColPart()[0];
+   cols[1] = Mfull->ColPart()[1];
+   int nrows = reduced_merged->Height();
+
+   int row_offset;
+   MPI_Scan(&nrows,&row_offset,1,MPI_INT,MPI_SUM,Mfull->GetComm());
+
+   row_offset-=nrows;
+   rows[0] = row_offset;
+   rows[1] = row_offset+nrows;
+   int glob_nrows;
+   MPI_Allreduce(&nrows, &glob_nrows,1,MPI_INT,MPI_SUM,Mfull->GetComm());
+
+
+   int glob_ncols = reduced_merged->Width();
+   M = new HypreParMatrix(Mfull->GetComm(), nrows, glob_nrows,
+                          glob_ncols, reduced_merged->GetI(), reduced_merged->GetJ(),
+                          reduced_merged->GetData(), rows,cols); 
+
+   Vector gap;
+   tribol::getMfemGap(coupling_scheme_id, gap);
+   auto& P_submesh = *pressure.ParFESpace()->GetProlongationMatrix();
+   Vector gap_true;
+   gap_true.SetSize(P_submesh.Width());
+   P_submesh.MultTranspose(gap,gap_true);
+
+   gapv.SetSize(nrows);
+   for (int i = 0; i<nrows; i++)
+   {
+      gapv[i] = gap_true[nonzero_rows[i]];
+   }
+
+   constraints_starts.SetSize(2);
+   constraints_starts[0] = M->RowPart()[0];
+   constraints_starts[1] = M->RowPart()[1];
+
+   // find elast dofs in contact;
+   HypreParMatrix * Jt = (HypreParMatrix *)(&A_blk->GetBlock(0,1));
+   Jt->EliminateRows(prob->GetEssentialDofs());
+
+   int hJt = Jt->Height();
+   SparseMatrix mergedJt;
+   Jt->MergeDiagAndOffd(mergedJt);
+
+   Array<int> nonzerorows;
+   Array<int> zerorows;
+   for (int i = 0; i<hJt; i++)
+   {
+      if (!mergedJt.RowIsEmpty(i))
+      {
+         nonzerorows.Append(i);
+      }
+      else
+      {
+         zerorows.Append(i);
+      }
+   }
+
+   int hb = nonzerorows.Size();
+   SparseMatrix Pbt(hb,K->GetGlobalNumCols());
+
+   for (int i = 0; i<hb; i++)
+   {
+      int col = nonzerorows[i]+prob->GetFESpace()->GetMyTDofOffset();
+      Pbt.Set(i,col,1.0);
+   }
+   Pbt.Finalize();
+
+   int rows_b[2];
+   int cols_b[2];
+   int nrows_b = Pbt.Height();
+
+   int row_offset_b;
+   MPI_Scan(&nrows_b,&row_offset_b,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+
+   row_offset_b-=nrows_b;
+   rows_b[0] = row_offset_b;
+   rows_b[1] = row_offset_b+nrows_b;
+   cols_b[0] = K->ColPart()[0];
+   cols_b[1] = K->ColPart()[1];
+   int glob_nrows_b;
+   int glob_ncols_b = K->GetGlobalNumCols();
+   MPI_Allreduce(&nrows_b, &glob_nrows_b,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+
+   HypreParMatrix * P_bt = new HypreParMatrix(MPI_COMM_WORLD, nrows_b, glob_nrows_b,
+                            glob_ncols_b, Pbt.GetI(), Pbt.GetJ(),
+                            Pbt.GetData(), rows_b,cols_b); 
+
+   Pb = P_bt->Transpose();
+   delete P_bt;                         
+
+   int hi = zerorows.Size();
+   SparseMatrix Pit(hi,K->GetGlobalNumCols());
+
+   for (int i = 0; i<hi; i++)
+   {
+      int col = zerorows[i]+prob->GetFESpace()->GetMyTDofOffset();
+      Pit.Set(i,col,1.0);
+   }
+   Pit.Finalize();
+
+   int rows_i[2];
+   int cols_i[2];
+   int nrows_i = Pit.Height();
+
+   int row_offset_i;
+   MPI_Scan(&nrows_i,&row_offset_i,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+
+   row_offset_i-=nrows_i;
+   rows_i[0] = row_offset_i;
+   rows_i[1] = row_offset_i+nrows_i;
+   cols_i[0] = K->ColPart()[0];
+   cols_i[1] = K->ColPart()[1];
+   int glob_nrows_i;
+   int glob_ncols_i = K->GetGlobalNumCols();
+   MPI_Allreduce(&nrows_i, &glob_nrows_i,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+
+   HypreParMatrix * P_it = new HypreParMatrix(MPI_COMM_WORLD, nrows_i, glob_nrows_i,
+                            glob_ncols_i, Pit.GetI(), Pit.GetJ(),
+                            Pit.GetData(), rows_i,cols_i); 
+   
+   Pi = P_it->Transpose();
+   delete P_it;
+}
+
+void ParContactProblem::SetupTribolDoublePass()
+{
+   axom::slic::SimpleLogger logger1;
+   axom::slic::setIsRoot(mfem::Mpi::Root());
+
+   // Initialize Tribol contact library
+   tribol::initialize(3, MPI_COMM_WORLD);
+
+   int coupling_scheme_id1 = 0;
+   int mesh1_id1 = 0;
+   int mesh2_id1 = 1;
+   vfes = prob->GetFESpace();
+   ParGridFunction * coords1 = new ParGridFunction(vfes);
+   ParMesh * pmesh1 = prob->GetMesh();
+   pmesh1->SetNodalGridFunction(coords1);
+   tribol::registerMfemCouplingScheme(
+      coupling_scheme_id1, mesh1_id1, mesh2_id1,
+      *pmesh1, *coords1, mortar_attrs, nonmortar_attrs,
+      tribol::SURFACE_TO_SURFACE,
+      tribol::NO_SLIDING,
+      tribol::SINGLE_MORTAR,
+      tribol::FRICTIONLESS,
+      tribol::LAGRANGE_MULTIPLIER,
+      tribol::BINNING_GRID
+   );
+
+   // Access Tribol's pressure grid function (on the contact surface)
+   auto& pressure1 = tribol::getMfemPressure(coupling_scheme_id1);
+   if (mfem::Mpi::Root())
+   {
+      std::cout << "Number of pressure unknowns: " <<
+                pressure1.ParFESpace()->GlobalTrueVSize() << std::endl;
+   }
+
+   // Set Tribol options for Lagrange multiplier enforcement
+   tribol::setLagrangeMultiplierOptions(
+      coupling_scheme_id1,
+      tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
+   );
+
+   // Update contact mesh decomposition
+   tribol::updateMfemParallelDecomposition();
+
+   // Update contact gaps, forces, and tangent stiffness
+   int cycle1 = 1;   // pseudo cycle
+   double t1 = 1.0;  // pseudo time
+   double dt1 = 1.0; // pseudo dt
+   tribol::update(cycle1, t1, dt1);
+
+   // Return contact contribution to the tangent stiffness matrix
+   auto A_blk1 = tribol::getMfemBlockJacobian(coupling_scheme_id1);
+   
+   HypreParMatrix * Mfull1 = (HypreParMatrix *)(&A_blk1->GetBlock(1,0));
+   Mfull1->EliminateCols(prob->GetEssentialDofs());
+   int h1 = Mfull1->Height();
+   SparseMatrix merged1;
+   Mfull1->MergeDiagAndOffd(merged1);
+   Array<int> nonzero_rows1;
+   for (int i = 0; i<h1; i++)
+   {
+      if (!merged1.RowIsEmpty(i))
+      {
+         nonzero_rows1.Append(i);
+      }
+   }
+
+   int hnew1 = nonzero_rows1.Size();
+   SparseMatrix P1(hnew1,h1);
+
+   for (int i = 0; i<hnew1; i++)
+   {
+      int col = nonzero_rows1[i];
+      P1.Set(i,col,1.0);
+   }
+   P1.Finalize();
+
+   SparseMatrix * reduced_merged1 = Mult(P1,merged1);
+
+   int rows1[2];
+   int cols1[2];
+   cols1[0] = Mfull1->ColPart()[0];
+   cols1[1] = Mfull1->ColPart()[1];
+   int nrows1 = reduced_merged1->Height();
+
+   int row_offset1;
+   MPI_Scan(&nrows1,&row_offset1,1,MPI_INT,MPI_SUM,Mfull1->GetComm());
+
+   row_offset1-=nrows1;
+   rows1[0] = row_offset1;
+   rows1[1] = row_offset1+nrows1;
+   int glob_nrows1;
+   MPI_Allreduce(&nrows1, &glob_nrows1,1,MPI_INT,MPI_SUM,Mfull1->GetComm());
+
+
+   int glob_ncols1 = reduced_merged1->Width();
+   HypreParMatrix * M1 = new HypreParMatrix(Mfull1->GetComm(), nrows1, glob_nrows1,
+                          glob_ncols1, reduced_merged1->GetI(), reduced_merged1->GetJ(),
+                          reduced_merged1->GetData(), rows1,cols1); 
+
+   Vector gap1;
+   tribol::getMfemGap(coupling_scheme_id1, gap1);
+   auto& P_submesh1 = *pressure1.ParFESpace()->GetProlongationMatrix();
+   Vector gap_true1;
+   gap_true1.SetSize(P_submesh1.Width());
+   P_submesh1.MultTranspose(gap1,gap_true1);
+
+   tribol::finalize();
+
+   // ------------------------------ 
+   // second pass
+   // ------------------------------ 
+   // Initialize Tribol contact library
+   tribol::initialize(3, MPI_COMM_WORLD);
+
+   int coupling_scheme_id2 = 0;
+   int mesh1_id2 = 0;
+   int mesh2_id2 = 1;
+   ParGridFunction * coords2 = new ParGridFunction(vfes);
+   ParMesh * pmesh2 = prob->GetMesh();
+   pmesh2->SetNodalGridFunction(coords2);
+   tribol::registerMfemCouplingScheme(
+      coupling_scheme_id2, mesh1_id2, mesh2_id2,
+      *pmesh2, *coords2, nonmortar_attrs, mortar_attrs,
+      tribol::SURFACE_TO_SURFACE,
+      tribol::NO_SLIDING,
+      tribol::SINGLE_MORTAR,
+      tribol::FRICTIONLESS,
+      tribol::LAGRANGE_MULTIPLIER,
+      tribol::BINNING_GRID
+   );
+
+   // Access Tribol's pressure grid function (on the contact surface)
+   auto& pressure2 = tribol::getMfemPressure(coupling_scheme_id2);
+   if (mfem::Mpi::Root())
+   {
+      std::cout << "Number of pressure unknowns: " <<
+                pressure2.ParFESpace()->GlobalTrueVSize() << std::endl;
+   }
+
+   // Set Tribol options for Lagrange multiplier enforcement
+   tribol::setLagrangeMultiplierOptions(
+      coupling_scheme_id2,
+      tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
+   );
+
+   // Update contact mesh decomposition
+   tribol::updateMfemParallelDecomposition();
+
+   // Update contact gaps, forces, and tangent stiffness
+   int cycle2 = 1;   // pseudo cycle
+   double t2 = 1.0;  // pseudo time
+   double dt2 = 1.0; // pseudo dt
+   tribol::update(cycle2, t2, dt2);
+
+   // Return contact contribution to the tangent stiffness matrix
+   auto A_blk2 = tribol::getMfemBlockJacobian(coupling_scheme_id2);
+   
+   HypreParMatrix * Mfull2 = (HypreParMatrix *)(&A_blk2->GetBlock(1,0));
+   Mfull2->EliminateCols(prob->GetEssentialDofs());
+   int h2 = Mfull2->Height();
+   SparseMatrix merged2;
+   Mfull2->MergeDiagAndOffd(merged2);
+   Array<int> nonzero_rows2;
+   for (int i = 0; i<h2; i++)
+   {
+      if (!merged2.RowIsEmpty(i))
+      {
+         nonzero_rows2.Append(i);
+      }
+   }
+
+   int hnew2 = nonzero_rows2.Size();
+   SparseMatrix P2(hnew2,h2);
+
+   for (int i = 0; i<hnew2; i++)
+   {
+      int col = nonzero_rows2[i];
+      P2.Set(i,col,1.0);
+   }
+   P2.Finalize();
+
+   SparseMatrix * reduced_merged2 = Mult(P2,merged2);
+
+   int rows2[2];
+   int cols2[2];
+   cols2[0] = Mfull2->ColPart()[0];
+   cols2[1] = Mfull2->ColPart()[1];
+   int nrows2 = reduced_merged2->Height();
+
+   int row_offset2;
+   MPI_Scan(&nrows2,&row_offset2,1,MPI_INT,MPI_SUM,Mfull2->GetComm());
+
+   row_offset2-=nrows2;
+   rows2[0] = row_offset2;
+   rows2[1] = row_offset2+nrows2;
+   int glob_nrows2;
+   MPI_Allreduce(&nrows2, &glob_nrows2,1,MPI_INT,MPI_SUM,Mfull2->GetComm());
+
+
+   int glob_ncols2 = reduced_merged2->Width();
+   HypreParMatrix * M2 = new HypreParMatrix(Mfull2->GetComm(), nrows2, glob_nrows2,
+                          glob_ncols2, reduced_merged2->GetI(), reduced_merged2->GetJ(),
+                          reduced_merged2->GetData(), rows2,cols2); 
+
+   Vector gap2;
+   tribol::getMfemGap(coupling_scheme_id2, gap2);
+   auto& P_submesh2 = *pressure2.ParFESpace()->GetProlongationMatrix();
+   Vector gap_true2;
+   gap_true2.SetSize(P_submesh2.Width());
+   P_submesh2.MultTranspose(gap2,gap_true2);
+
+   tribol::finalize();
+
+
+   gapv.SetSize(nrows1+nrows2);
+   for (int i = 0; i<nrows1; i++)
+   {
+      gapv[i] = gap_true1[nonzero_rows1[i]];
+   }
+   for (int i = 0; i<nrows2; i++)
+   {
+      gapv[nrows1+i] = gap_true2[nonzero_rows2[i]];
+   }
+
+   Array2D<HypreParMatrix *> A_array(2,1);
+   A_array(0,0) = M1;
+   A_array(1,0) = M2;
+
+   M = HypreParMatrixFromBlocks(A_array);
+
+   constraints_starts.SetSize(2);
+   constraints_starts[0] = M->RowPart()[0];
+   constraints_starts[1] = M->RowPart()[1];
+
+}
+
+
+
+double ParContactProblem::E(const Vector & d)
+{
+   Vector kd(K->Height());
+   K->Mult(d,kd);
+   return 0.5 * InnerProduct(comm,d, kd) - InnerProduct(comm,d, *B);
+}
+
+void ParContactProblem::DdE(const Vector &d, Vector &gradE)
+{
+   gradE.SetSize(K->Height());
+   K->Mult(d, gradE);
+   gradE.Add(-1.0, *B); 
+}
+
+HypreParMatrix* ParContactProblem::DddE(const Vector &d)
+{
+   return K; 
+}
+
+void ParContactProblem::g(const Vector &d, Vector &gd)
+{
+   gd = GetGapFunction();
+}
+
+HypreParMatrix* ParContactProblem::Ddg(const Vector &d)
+{
+  return GetJacobian();
+}
+
+HypreParMatrix* ParContactProblem::lDddg(const Vector &d, const Vector &l)
+{
+   return nullptr; // for now
+}
+
+
+QPOptParContactProblem::QPOptParContactProblem(ParContactProblem * problem_, Vector &xref_)
+: problem(problem_)
+{
+   dimU = problem->GetNumDofs();
+   dimM = problem->GetNumConstraints();
+   dimC = problem->GetNumConstraints();
+   ml.SetSize(dimM); ml = 0.0;
+   Vector negone(dimM); negone = -1.0;
+   SparseMatrix diag(negone);
+
+   xref.SetSize(xref_.Size());
+   xref.Set(1.0, xref_);
+
+   int gsize = problem->GetGlobalNumConstraints();
+   int * rows = problem->GetConstraintsStarts().GetData();
+
+   NegId = new HypreParMatrix(problem->GetComm(),gsize, rows,&diag);
+   HypreStealOwnership(*NegId, diag);
+}
+
+int QPOptParContactProblem::GetDimU() { return dimU; }
+
+int QPOptParContactProblem::GetDimM() { return dimM; }
+
+int QPOptParContactProblem::GetDimC() { return dimC; }
+
+Vector & QPOptParContactProblem::Getml() { return ml; }
+
+HypreParMatrix * QPOptParContactProblem::Duuf(const BlockVector & x)
+{
+   return problem->DddE(x.GetBlock(0));
+}
+
+HypreParMatrix * QPOptParContactProblem::Dumf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+HypreParMatrix * QPOptParContactProblem::Dmuf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+HypreParMatrix * QPOptParContactProblem::Dmmf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+HypreParMatrix * QPOptParContactProblem::Duc(const BlockVector & x)
+{
+   return problem->Ddg(x.GetBlock(0));
+}
+
+HypreParMatrix * QPOptParContactProblem::Dmc(const BlockVector & x)
+{
+   return NegId;
+}
+
+HypreParMatrix * QPOptParContactProblem::lDuuc(const BlockVector & x, const Vector & l)
+{
+   return nullptr;
+}
+
+
+void QPOptParContactProblem::c(const BlockVector &x, Vector & y)
+{
+   Vector g0; // g(dref) 
+   problem->g(x.GetBlock(0), g0); // gap function
+
+   // temp = d - xref (expansion)
+   Vector temp(x.GetBlock(0).Size()); temp = 0.0;
+   temp.Set(1.0, x.GetBlock(0)); 
+   temp.Add(-1.0, xref); // displacement at previous time step  
+
+   problem->GetJacobian()->Mult(temp, y); // J * (d - xref)
+   y.Add(1.0, g0); // J * (d - xref) + g0 
+   y.Add(-1.0, x.GetBlock(1)); // J * (d - xref) + g0 - s
+}
+
+double QPOptParContactProblem::CalcObjective(const BlockVector & x)
+{
+   return problem->E(x.GetBlock(0));
+}
+
+void QPOptParContactProblem::CalcObjectiveGrad(const BlockVector & x, BlockVector & y)
+{
+   problem->DdE(x.GetBlock(0), y.GetBlock(0));
+   y.GetBlock(1) = 0.0;
+}
+
+QPOptParContactProblem::~QPOptParContactProblem()
+{
+   delete NegId;
+}
+
+// #endif

miniapps/contact/problems/parproblems.hpp

@@ -0,0 +1,335 @@
+
+#include "parproblems_util.hpp"
+
+class ParElasticityProblem
+{
+private:
+   MPI_Comm comm;
+   bool formsystem = false;
+   ParMesh * pmesh = nullptr;
+   Array<int> ess_bdr_attr, ess_bdr_attr_comp;
+   int order;
+   int ndofs;
+   int ntdofs;
+   int gndofs;
+   FiniteElementCollection * fec = nullptr;
+   ParFiniteElementSpace * fes = nullptr;
+   Vector lambda, mu;
+   PWConstCoefficient lambda_cf, mu_cf;
+   Array<int> ess_bdr, ess_tdof_list;
+   ParBilinearForm * a = nullptr;
+   ParLinearForm * b = nullptr;
+   ParGridFunction x;
+   HypreParMatrix A;
+   Vector B,X;
+   ConstantCoefficient pressure_cf;
+   VectorArrayCoefficient * bf = nullptr;
+   void Init();
+   bool own_mesh;
+public:
+   ParElasticityProblem(MPI_Comm comm_, const char *mesh_file , int sref, int pref, 
+                        Array<int> & ess_bdr_attr_, Array<int> & ess_bdr_attr_comp_, 
+                        int order_ = 1 ) 
+   : comm(comm_), ess_bdr_attr(ess_bdr_attr_),ess_bdr_attr_comp(ess_bdr_attr_comp_), order(order_)
+   {
+      own_mesh = true;
+      Mesh * mesh = new Mesh(mesh_file,1,1);
+      for (int i = 0; i<sref; i++)
+      {
+         mesh->UniformRefinement();
+      }
+      pmesh = new ParMesh(comm,*mesh);
+      MFEM_VERIFY(pmesh->GetNE(), "ParElasticityProblem::Empty partition");
+      delete mesh;
+      for (int i = 0; i<pref; i++)
+      {
+         pmesh->UniformRefinement();
+      }
+      Init();
+   }
+
+   ParElasticityProblem(ParMesh * pmesh_, Array<int> & ess_bdr_attr_, Array<int> & ess_bdr_attr_comp_,  int order_ = 1) 
+   :  pmesh(pmesh_), ess_bdr_attr(ess_bdr_attr_), ess_bdr_attr_comp(ess_bdr_attr_comp_), order(order_)
+   {
+      own_mesh = false;
+      comm = pmesh->GetComm();
+      Init();
+   }
+
+   ParMesh * GetMesh() { return pmesh; }
+   ParFiniteElementSpace * GetFESpace() { return fes; }
+   FiniteElementCollection * GetFECol() { return fec; }
+   int GetNumDofs() { return ndofs; }
+   int GetNumTDofs() { return ntdofs; }
+   int GetGlobalNumDofs() { return gndofs; }
+   HypreParMatrix & GetOperator() 
+   { 
+      MFEM_VERIFY(formsystem, "System not formed yet. Call FormLinearSystem()"); 
+      return A; 
+   }
+   Vector & GetRHS() 
+   { 
+      MFEM_VERIFY(formsystem, "System not formed yet. Call FormLinearSystem()"); 
+      return B; 
+   }
+
+   void SetLambda(const Vector & lambda_) 
+   { 
+      lambda = lambda_; 
+      lambda_cf.UpdateConstants(lambda);
+   }
+   void SetMu(const Vector & mu_) 
+   { 
+      mu = mu_; 
+      mu_cf.UpdateConstants(mu);
+   }
+
+   void SetNeumanPressureData(ConstantCoefficient &f, Array<int> & bdr_marker)
+   { 
+      pressure_cf.constant = f.constant;
+      b->AddBoundaryIntegrator(new VectorBoundaryFluxLFIntegrator(pressure_cf),bdr_marker);
+   }
+
+   void SetNeumanData(int comp, int bdrattr, double value)
+   { 
+      int dim = pmesh->Dimension();
+      bf = new VectorArrayCoefficient(dim);
+      for (int i = 0; i < dim; i++)
+      {
+         if (i == comp)
+         {
+            Vector pull_force(pmesh->bdr_attributes.Max());
+            pull_force = 0.0;
+            pull_force(bdrattr-1) = value;
+            bf->Set(i, new PWConstCoefficient(pull_force));
+         }
+         else
+         {
+            bf->Set(i, new ConstantCoefficient(0.0));
+         }
+      }
+      b->AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(*bf));
+   }
+
+   void UpdateEssentialBC(Array<int> & ess_bdr_attr_, Array<int> & ess_bdr_attr_comp_)
+   {
+      ess_bdr_attr = ess_bdr_attr_;
+      ess_bdr_attr_comp = ess_bdr_attr_comp_;
+      ess_tdof_list.SetSize(0);
+      if (pmesh->bdr_attributes.Size())
+      {
+         ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      }
+      ess_bdr = 0; 
+      Array<int> ess_tdof_list_temp;
+      for (int i = 0; i < ess_bdr_attr.Size(); i++ )
+      {
+         ess_bdr[ess_bdr_attr[i]-1] = 1;
+         fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list_temp,ess_bdr_attr_comp[i]);
+         ess_tdof_list.Append(ess_tdof_list_temp);
+         ess_bdr[ess_bdr_attr[i]-1] = 0;
+      }
+   }
+
+   void UpdateStep()
+   {
+      if (formsystem)
+      {
+         delete b;
+         b = new ParLinearForm(fes);
+         delete a;
+         a = new ParBilinearForm(fes);
+         a->AddDomainIntegrator(new ElasticityIntegrator(lambda_cf,mu_cf));
+
+
+         // a->Update();
+         formsystem = false;
+      }
+   }
+
+   void FormLinearSystem();
+   void UpdateLinearSystem();
+
+   void SetDisplacementDirichletData(const Vector & delta) 
+   {
+      VectorConstantCoefficient delta_cf(delta);
+      x.ProjectBdrCoefficient(delta_cf,ess_bdr);
+      bool vis = false;
+      if (vis)
+      {
+         int myid, num_procs;
+         MPI_Comm_rank(comm, &myid);
+         MPI_Comm_size(comm, &num_procs);
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock(vishost, visport);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n";
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << *pmesh << x << std::flush;
+         MFEM_ABORT("");
+      }
+   };
+
+   void ResetDisplacementDirichletData()
+   {
+      x = 0.0;
+   }
+
+   void SetDisplacementDirichletData(const Vector & delta, Array<int> essbdr) 
+   {
+      VectorConstantCoefficient delta_cf(delta);
+      x.ProjectBdrCoefficient(delta_cf,essbdr);
+      bool vis = false;
+      if (vis)
+      {
+         int myid, num_procs;
+         MPI_Comm_rank(comm, &myid);
+         MPI_Comm_size(comm, &num_procs);
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock(vishost, visport);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n";
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << *pmesh << x << std::flush;
+         MFEM_ABORT("");
+      }
+   };
+
+   ParGridFunction & GetDisplacementGridFunction() {return x;};
+   Array<int> & GetEssentialDofs() {return ess_tdof_list;};
+
+   ~ParElasticityProblem()
+   {
+      delete a;
+      delete b;
+      delete fes;
+      delete fec;
+      if (own_mesh)
+      {
+         delete pmesh;
+      }
+      delete bf;
+   }
+};
+
+
+// #ifdef MFEM_USE_TRIBOL
+
+class ParContactProblem
+{
+private:
+   MPI_Comm comm;
+   int numprocs;
+   int myid;
+   ParElasticityProblem * prob = nullptr;
+   ParFiniteElementSpace * vfes = nullptr;
+   int dim;
+   GridFunction nodes0;
+   GridFunction *nodes1 = nullptr;
+   std::set<int> contact_vertices;
+   std::vector<int> dof_offsets;
+   std::vector<int> vertex_offsets; 
+   std::vector<int> constraints_offsets; 
+   Array<int> tdof_offsets;
+   Array<int> constraints_starts;
+   Array<int> globalvertices;
+   Array<int> vertices;
+   ParGridFunction * coords = nullptr;
+   //ParGridFunction * xref = nullptr;
+
+protected:
+   int npoints=0;
+   int gnpoints=0;
+   int nv, gnv;
+   HypreParMatrix * K = nullptr;
+   HypreParMatrix * Pi = nullptr;
+   HypreParMatrix * Pb = nullptr;
+   Vector *B = nullptr;
+   Vector gapv;
+   HypreParMatrix * M=nullptr;
+   void SetupTribol();
+   void SetupTribolDoublePass();
+   std::set<int> mortar_attrs;
+   // plane of top block
+   std::set<int> nonmortar_attrs;
+   bool doublepass = false;
+
+public:
+   ParContactProblem(ParElasticityProblem * prob_, 
+                     const std::set<int> & mortar_attrs_, const std::set<int> & nonmortar_attrs_,
+                     ParGridFunction * coords_,
+                     bool doublepass = false);
+
+   ParElasticityProblem * GetElasticityProblem() {return prob;}
+   MPI_Comm GetComm() {return comm;}
+   int GetNumDofs() {return K->Height();}
+   int GetGlobalNumDofs() {return K->GetGlobalNumRows();}
+   int GetNumConstraints() {return M->Height();}
+   int GetGlobalNumConstraints() { return M->GetGlobalNumRows(); }
+   
+   std::vector<int> & GetDofOffets() { return dof_offsets; }
+   std::vector<int> & GetVertexOffsets() { return vertex_offsets; }
+   std::vector<int> & GetConstraintsOffsets() { return constraints_offsets; }
+   Array<int> & GetConstraintsStarts() { return constraints_starts; }
+   
+   Vector & GetGapFunction() {return gapv;}
+
+   HypreParMatrix * GetJacobian() {return M;}
+
+   double E(const Vector & d);
+   void DdE(const Vector &d, Vector &gradE);
+   HypreParMatrix* DddE(const Vector &d);
+   void g(const Vector &d, Vector &gd);
+   HypreParMatrix* Ddg(const Vector &d);
+   HypreParMatrix* lDddg(const Vector &d, const Vector &l);
+
+   HypreParMatrix * GetRestrictionToInteriorDofs() {return Pi;}
+   HypreParMatrix * GetRestrictionToContactDofs() {return Pb;}
+
+   ~ParContactProblem()
+   {
+      delete B;
+      delete K;
+      delete M;
+   }
+};
+
+
+class QPOptParContactProblem
+{
+private:
+   ParContactProblem * problem = nullptr;
+   int dimU, dimM, dimC;
+   Vector ml;
+   HypreParMatrix * NegId = nullptr;
+   Vector xref;
+public:
+   QPOptParContactProblem(ParContactProblem * problem_, Vector & xref_);
+   int GetDimU();
+   int GetDimM();
+   int GetDimC();
+   Vector & Getml();
+   MPI_Comm GetComm() {return problem->GetComm();}
+   int * GetConstraintsStarts() {return problem->GetConstraintsStarts().GetData();}
+   int GetGlobalNumConstraints() {return problem->GetGlobalNumConstraints();}
+
+   ParElasticityProblem * GetElasticityProblem() {return problem->GetElasticityProblem();}
+
+   HypreParMatrix * Duuf(const BlockVector &);
+   HypreParMatrix * Dumf(const BlockVector &);
+   HypreParMatrix * Dmuf(const BlockVector &);
+   HypreParMatrix * Dmmf(const BlockVector &);
+   HypreParMatrix * Duc(const BlockVector &);
+   HypreParMatrix * Dmc(const BlockVector &);
+   HypreParMatrix * lDuuc(const BlockVector &, const Vector &);
+
+   HypreParMatrix * GetRestrictionToInteriorDofs() {return problem->GetRestrictionToInteriorDofs();}
+   HypreParMatrix * GetRestrictionToContactDofs() {return problem->GetRestrictionToContactDofs();}
+
+   void c(const BlockVector &, Vector &);
+   double CalcObjective(const BlockVector &);
+   void CalcObjectiveGrad(const BlockVector &, BlockVector &);
+   ~QPOptParContactProblem();
+};
+
+// #endif

miniapps/contact/problems/parproblems_util.cpp

@@ -0,0 +1,115 @@
+#include "parproblems_util.hpp"
+
+int get_rank(int tdof, std::vector<int> & tdof_offsets)
+{
+   int size = tdof_offsets.size();
+   if (size == 1) { return 0; }
+   std::vector<int>::iterator up;
+   up=std::upper_bound(tdof_offsets.begin(), tdof_offsets.end(),tdof); //
+   return std::distance(tdof_offsets.begin(),up)-1;
+}
+
+void ComputeTdofOffsets(const ParFiniteElementSpace * pfes,
+                        std::vector<int> & tdof_offsets)
+{
+   MPI_Comm comm = pfes->GetComm();
+   int num_procs;
+   MPI_Comm_size(comm, &num_procs);
+   tdof_offsets.resize(num_procs);
+   int mytoffset = pfes->GetMyTDofOffset();
+   MPI_Allgather(&mytoffset,1,MPI_INT,&tdof_offsets[0],1,MPI_INT,comm);
+}
+
+void ComputeTdofOffsets(MPI_Comm comm, int mytoffset, std::vector<int> & tdof_offsets)
+{
+   int num_procs;
+   MPI_Comm_size(comm,&num_procs);
+   tdof_offsets.resize(num_procs);
+   MPI_Allgather(&mytoffset,1,MPI_INT,&tdof_offsets[0],1,MPI_INT,comm);
+}
+
+void ComputeTdofs(MPI_Comm comm, int mytoffs, std::vector<int> & tdofs)
+{
+   int num_procs;
+   MPI_Comm_size(comm,&num_procs);
+   tdofs.resize(num_procs);
+   MPI_Allgather(&mytoffs,1,MPI_INT,&tdofs,1,MPI_INT,comm);
+}         
+
+
+// Performs Pᵀ * A * P for BlockOperator  P (with blocks as HypreParMatrices)
+// and A a HypreParMatrix, i.e., this handles the special case 
+// where P = [P₁ P₂ ⋅⋅⋅ Pₙ] 
+// C = Pᵀ * A * P 
+void RAP(const HypreParMatrix & A, const BlockOperator & P, 
+         BlockOperator & C)
+{
+   int nblocks = P.NumColBlocks();
+
+   const HypreParMatrix * Pi = nullptr;
+   const HypreParMatrix * Pj = nullptr;
+   HypreParMatrix * PitAPj = nullptr; 
+
+   for (int i = 0; i< nblocks; i++)
+   {
+      if (P.IsZeroBlock(0,i)) continue;
+      Pi = dynamic_cast<const HypreParMatrix*>(&P.GetBlock(0,i));
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (P.IsZeroBlock(0,j)) continue;
+         Pj = dynamic_cast<const HypreParMatrix*>(&P.GetBlock(0,j));
+         if (i == j) 
+         {
+            PitAPj = RAP(&A, Pj); 
+         }
+         else
+         {
+            PitAPj = RAP(Pi, &A, Pj); 
+         }
+         C.SetBlock(i,j,PitAPj);
+      }
+   }
+}
+
+void ParAdd(const BlockOperator & A, const BlockOperator & B, BlockOperator & C)
+{
+   int n = A.NumRowBlocks();
+   int m = A.NumColBlocks();
+   MFEM_VERIFY(B.NumRowBlocks() == n, "Inconsistent number of row blocks");
+   MFEM_VERIFY(B.NumColBlocks() == m, "Inconsistent number of column blocks");
+
+   const HypreParMatrix * a;
+   const HypreParMatrix * b;
+   for (int i = 0; i<n; i++)
+   {
+      for (int j = 0; j<m; j++)
+      {
+         a = nullptr;
+         b = nullptr;
+         if (!A.IsZeroBlock(i,j))
+         {
+            a = dynamic_cast<const HypreParMatrix*>(&A.GetBlock(i,j));
+         }
+         if (!B.IsZeroBlock(i,j))
+         {
+            b = dynamic_cast<const HypreParMatrix*>(&B.GetBlock(i,j));
+         }
+         if (a && b)
+         {
+            C.SetBlock(i,j,ParAdd(a,b));
+         }
+         else if (a)
+         {
+            C.SetBlock(i,j,new HypreParMatrix(*a));
+         }
+         else if (b)
+         {
+            C.SetBlock(i,j,new HypreParMatrix(*b));
+         }
+         else
+         {
+            // do nothing
+         }
+      }
+   }
+}

miniapps/contact/problems/parproblems_util.hpp

@@ -0,0 +1,23 @@
+
+#include "mfem.hpp"
+
+
+using namespace std;
+using namespace mfem;
+
+#include "axom/slic.hpp"
+
+#include "tribol/interface/tribol.hpp"
+#include "tribol/interface/mfem_tribol.hpp"
+int get_rank(int tdof, std::vector<int> & tdof_offsets);
+void ComputeTdofOffsets(const ParFiniteElementSpace * pfes,
+                        std::vector<int> & tdof_offsets);
+void ComputeTdofOffsets(MPI_Comm comm, int mytoffset, std::vector<int> & tdof_offsets);
+void ComputeTdofs(MPI_Comm comm, int mytoffs, std::vector<int> & tdofs);
+
+
+// Performs Pᵀ * A * P for BlockOperator  P (with blocks as HypreParMatrices)
+// and A a HypreParMatrix, i.e., this handles the special case 
+// where P = [P₁ P₂ ⋅⋅⋅ Pₙ] 
+void RAP(const HypreParMatrix & A, const BlockOperator & P, BlockOperator & C);
+void ParAdd(const BlockOperator & A, const BlockOperator & B, BlockOperator & C);

miniapps/contact/scratch.cpp

@@ -0,0 +1,564 @@
+//                               Parallel contact example
+// mpirun -np 4 ./contact -ls 2 -sr 1 -testno 4
+// CG iteration numbers            = 105 114 116 115 113 109 113 108 107 114 206 236 268 435 987
+
+// mpirun -np 4 ./contact -ls 2 -sr 0 -testno 5
+// CG iteration numbers            = 106 116 116 116 115 113 107 107 128 131 531 1437 1318
+
+// mpirun -np 4 ./contact -ls 2 -sr 0 -testno 6
+// CG iteration numbers            = 18 18 18 18 18 17 17 21 22 46 52 53
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "ipsolver/ParIPsolver.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
+   Hypre::Init();
+
+   int order = 1;
+   int sref = 1;
+   int pref = 0;
+   Array<int> attr;
+   Array<int> m_attr;
+   bool visualization = true;
+   bool paraview = false;
+   bool elast = false;
+   bool nocontact = false;
+   int testNo = -1; // 0-6
+   // 1. Parse command-line options.
+   OptionsParser args(argc, argv);
+   args.AddOption(&testNo, "-testno", "--test-number",
+                  "Choice of test problem:"
+                  "-1: default (original 2 block problem)"
+                  "0: not implemented yet"
+                  "1: not implemented yet"
+                  "2: not implemented yet"
+                  "3: not implemented yet"
+                  "4: two block problem - diablo"
+                  "41: two block problem - twisted"
+                  "5: ironing problem"
+                  "51: ironing problem extended"
+                  "6: nested spheres problem");
+   args.AddOption(&attr, "-at", "--attributes-surf",
+                  "Attributes of boundary faces on contact surface for mesh 2.");
+   args.AddOption(&sref, "-sr", "--serial-refinements",
+                  "Number of uniform refinements.");                  
+   args.AddOption(&pref, "-pr", "--parallel-refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&elast, "-elast", "--elast", "-no-elast",
+                  "--no-elast",
+                  "Enable or disable AMG Elasticity options.");
+   args.AddOption(&nocontact, "-nocontact", "--nocontact", "-no-nocontact",
+                  "--no-nocontact",
+                  "Enable or disable AMG solve with no contact for testing.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
+                  "--no-paraview",
+                  "Enable or disable ParaView visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   if (Mpi::Root())
+   {
+      mfem::out << "Solving test problem number: " << testNo << endl;
+   }
+
+   const char *mesh_file = nullptr;
+
+   switch (testNo)
+   {
+      case -1:
+         mesh_file = "meshes/two-block.mesh";
+         break;
+      case 0:
+      case 1:
+      case 2:
+      case 3:
+      {
+         MFEM_ABORT("Problem not implemented yet");
+         break;
+      }
+      case 4:
+         mesh_file = "meshes/Test4.mesh";
+         break;
+      case 40:
+         mesh_file = "meshes/Test40.mesh";
+         break;   
+      case 41:
+         mesh_file = "meshes/Test41.mesh";
+         break;
+      case 42:
+         mesh_file = "meshes/Test42.mesh";
+         break;         
+      case 5:
+         mesh_file = "meshes/Test5.mesh";
+         break;
+      case 51:
+         mesh_file = "meshes/Test51.mesh";
+         break;
+      case 6:
+         mesh_file = "meshes/Test6.mesh";
+         break;
+      case 61:
+         // Something wrong with this mesh
+         mesh_file = "meshes/Test61.mesh";
+         break;
+      case 62:
+         mesh_file = "meshes/Test62.mesh";
+         break;         
+      default:
+         MFEM_ABORT("Should be unreachable");
+         break;
+   }
+
+   Mesh * mesh = new Mesh(mesh_file,1);
+   for (int i = 0; i<sref; i++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   ParMesh * pmesh = new ParMesh(MPI_COMM_WORLD,*mesh);
+
+   for (int i = 0; i<pref; i++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   Array<int> ess_bdr_attr;
+   Array<int> ess_bdr_attr_comp;
+   if (testNo == 6 || testNo == 61)
+   {
+      ess_bdr_attr.Append(1); ess_bdr_attr_comp.Append(1);
+      ess_bdr_attr.Append(2); ess_bdr_attr_comp.Append(2);
+      ess_bdr_attr.Append(4); ess_bdr_attr_comp.Append(0);
+      ess_bdr_attr.Append(5); ess_bdr_attr_comp.Append(-1);
+   }
+   else if (testNo == 62)
+   {
+      ess_bdr_attr.Append(4); ess_bdr_attr_comp.Append(0);
+      ess_bdr_attr.Append(5); ess_bdr_attr_comp.Append(-1);
+   }
+   else if (testNo == 40)
+   {
+      ess_bdr_attr.Append(1); ess_bdr_attr_comp.Append(-1);
+      ess_bdr_attr.Append(10); ess_bdr_attr_comp.Append(-1);
+   }
+   else
+   {
+      ess_bdr_attr.Append(2); ess_bdr_attr_comp.Append(-1);
+      ess_bdr_attr.Append(6); ess_bdr_attr_comp.Append(-1);
+   }
+   ParElasticityProblem * prob = new ParElasticityProblem(pmesh,
+                                                          ess_bdr_attr,ess_bdr_attr_comp,
+                                                          order);
+   Vector lambda(prob->GetMesh()->attributes.Max());
+   Vector mu(prob->GetMesh()->attributes.Max());
+
+   if (testNo == -1 )
+   {
+      lambda = 57.6923076923;
+      mu = 38.4615384615;
+   }
+   else if (testNo == 6 || testNo == 61 || testNo == 62)
+   {
+      lambda = (1000*0.3)/(1.3*0.4);
+      mu = 500/(1.3);
+   }
+   else
+   {
+      lambda[0] = 0.499/(1.499*0.002);
+      lambda[1] = 0.0;
+      mu[0] = 1./(2*1.499);
+      mu[1] = 500;
+   }
+
+   prob->SetLambda(lambda); prob->SetMu(mu);
+
+   int dim = pmesh->Dimension();
+   Vector ess_values(dim);
+   int essbdr_attr;
+   Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+
+   ess_values = 0.0;
+
+
+   // ConstantCoefficient one(-area);
+   ConstantCoefficient one(-1.0);
+
+   std::set<int> mortar_attr;
+   std::set<int> nonmortar_attr;
+
+   int nsteps = 100;
+   if (testNo == 6 || testNo == 61)
+   {
+      ess_values = 0.0;
+      ess_bdr = 0;
+      ess_bdr[0] = 1;
+      ess_bdr[1] = 1;
+      ess_bdr[3] = 1;
+      ess_bdr[4] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      ess_bdr = 0;
+      ess_bdr[2] = 1;
+      // prob->SetNeumanPressureData(one,ess_bdr);
+      mortar_attr.insert(6);
+      mortar_attr.insert(9);
+      nonmortar_attr.insert(7);
+      nonmortar_attr.insert(8);
+   }
+   else if(testNo == 62)
+   {
+      ess_values = 0.0;
+      ess_bdr = 0;
+      ess_bdr[3] = 1;
+      ess_bdr[4] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      ess_bdr = 0;
+      ess_bdr[2] = 1;
+      // prob->SetNeumanPressureData(one,ess_bdr);
+      prob->SetNeumanData(0,3,-2.0);
+      mortar_attr.insert(6);
+      mortar_attr.insert(9);
+      nonmortar_attr.insert(7);
+      nonmortar_attr.insert(8);
+   }
+   else
+   {
+      if (testNo == -1 || testNo == 41)
+      {
+         ess_values[0] = 0.1/nsteps;
+      }
+      else
+      {
+         ess_values[2] = 1.0/1.4/nsteps;
+         // ess_values[0] = -2.0/nsteps;
+      }
+      essbdr_attr = (testNo == 40) ? 1 : 2;
+      ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+      // prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      essbdr_attr = (testNo == 40) ? 10 : 6;
+      ess_values = 0.0; ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+      // prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      if (testNo == 40)
+      {
+         mortar_attr.insert(4);
+         nonmortar_attr.insert(7);
+      }
+      else
+      {
+         mortar_attr.insert(3);
+         nonmortar_attr.insert(4);
+      }
+   }
+
+   ParFiniteElementSpace * fes = prob->GetFESpace();
+   ParGridFunction x_gf(fes); x_gf = 0.0;
+   ParGridFunction xnew(fes); xnew = 0.0;
+   ParaViewDataCollection * paraview_dc = nullptr;
+   ParMesh pmesh_copy(*pmesh);
+   ParFiniteElementSpace fes_copy(*fes,pmesh_copy);
+   ParGridFunction xcopy_gf(&fes_copy); xcopy_gf = 0.0;
+   
+   if (paraview)
+   {
+      std::ostringstream paraview_file_name;
+      paraview_file_name << "QPContact-Test_" << testNo
+                         << "_par_ref_" << pref
+                         << "_ser_ref_" << sref;
+      paraview_dc = new ParaViewDataCollection(paraview_file_name.str(), &pmesh_copy);
+      paraview_dc->SetPrefixPath("ParaView");
+      paraview_dc->SetLevelsOfDetail(1);
+      paraview_dc->SetDataFormat(VTKFormat::BINARY);
+      paraview_dc->SetHighOrderOutput(true);
+      // paraview_dc->RegisterField("u", &x_gf);
+      paraview_dc->RegisterField("u", &xcopy_gf);
+      paraview_dc->SetCycle(0);
+      paraview_dc->SetTime(double(0));
+      paraview_dc->Save();
+   }
+   socketstream sol_sock;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      sol_sock.open(vishost, visport);
+      sol_sock.precision(8);
+   }
+   
+   
+   ParGridFunction ref_coords(prob->GetFESpace()); 
+   ParGridFunction new_coords(prob->GetFESpace()); 
+   pmesh->GetNodes(new_coords);
+   pmesh->GetNodes(ref_coords);
+   
+   Vector xref(x_gf.GetTrueVector().Size());
+
+   double p = 1;
+   ConstantCoefficient f(p);
+   
+   double pseudotime = 1.0 / ((double) nsteps);
+   if (testNo == 6)
+   {
+      ess_bdr = 0;
+      ess_bdr[2] = 1;
+      f.constant = -p * pseudotime;
+      prob->SetNeumanPressureData(f,ess_bdr);
+      // prob->SetNeumanData(0,3,-p*(i+1)/nsteps);
+   }
+   else if (testNo == 4 || testNo == 40 || testNo == 5 || testNo == 51)
+   {
+      ess_bdr = 0;
+      essbdr_attr = (testNo == 40) ? 1 : 2;
+      ess_bdr[essbdr_attr-1] = 1;
+      ess_values = 0.0;
+      //ess_values[2] = 4.0 / 7.0 * pseudotime;
+      ess_values[2] = 1.0/1.4 * pseudotime;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+   }
+   else if (testNo == 41)
+   {
+      ess_values = 0.0;
+      ess_values[0] = 0.5 * pseudotime; //0.5/nsteps*(i+1);
+      // ess_values[0] = 0.0;
+      essbdr_attr =  2;
+      ess_bdr[essbdr_attr-1] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      essbdr_attr = 6;
+      ess_values = 0.0; 
+      // ess_values[0] = -0.5/nsteps*(i+1);
+      if (myid == 0)
+      {
+         mfem::out << "ess_values[0] = " << ess_values[0] << endl;
+      }
+      ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+   }
+
+
+   /* ------- finite difference check -------- */
+   Vector x0(fes->GetTrueVSize()); x0 = 0.0;
+   //x0 = 2.0;
+   //x0.Randomize(); x0 *= 1.e-2;
+   Array<int> vdofs;
+   for (int i = 0; i < pmesh->GetNE(); i++)
+   {
+      cout << "attribute = " << pmesh->GetAttribute(i) << endl;
+      if (pmesh->GetAttribute(i) == 1)
+      {
+         continue;
+      }
+      fes->GetElementVDofs(i, vdofs);
+      for (int j = 0; j < vdofs.Size(); j++)
+      {
+         x0(vdofs[j]) = 0.01;
+      }
+   }
+   
+   
+   
+   
+   x_gf.SetFromTrueDofs(x0); 
+   add(ref_coords, x_gf, new_coords);
+
+   Vector x1(fes->GetTrueVSize()); x1 = 0.0;
+   Vector xdir(fes->GetTrueVSize()); xdir.Randomize();
+   Vector temp(fes->GetTrueVSize()); temp = 0.0;
+   xdir *= 1.e-2; // scale so as to avoid mesh tangling
+   double eps = 1.0;
+   
+   ParContactProblem ref_contact(prob, mortar_attr, nonmortar_attr, &new_coords);
+   int ndofs = ref_contact.GetNumDofs();
+   int nconstraints = ref_contact.GetNumConstraints();
+   Vector g0 = ref_contact.GetGapFunction();
+   g0.Print();
+   HypreParMatrix * J0 = ref_contact.GetJacobian();
+  
+   //for (int i = 0; i < 30; i++)
+   //{
+   //   x1.Set(1.0, x0);     // x1 = x0 + eps * xdir
+   //   x1.Add(eps, xdir);
+   //   x_gf.SetFromTrueDofs(x1);
+   //   add(ref_coords, x_gf, new_coords);
+   //   ParContactProblem new_contact(prob, mortar_attr, nonmortar_attr, &new_coords);
+   //   Vector g1 = new_contact.GetGapFunction(); // g1 = g(x0 + eps * xdir)
+   //   Vector fd_err(g1.Size());
+
+   //   // ||J0 * xdir - (g1 - g0) / eps||
+   //   J0->Mult(xdir, fd_err);
+   //   fd_err.Add(-1.0 / eps, g1);
+   //   fd_err.Add(1.0 / eps, g0);
+   //   cout << "fd err = " << fd_err.Norml2() << ", eps = " << eps << endl;
+   //   eps /= 2.0;
+   //}
+   
+   
+   
+   //for (int i = 0; i < 30; i++)
+   //{
+   ////      add(ref_coords,x_gf,new_coords);
+   //   
+   //}
+
+
+   //for (int i = 0; i < nsteps; i++)
+   //{
+   //   //pseudotime = ((double) (i) / ((double) SQPrepeat) + 1.) / ((double) nsteps);
+   //   pseudotime = ((double) (i)) / ((double) nsteps);
+   //   for (int j = 0; j < SQPrepeat; j++)
+   //   {
+   //      paraview_time = pseudotime + j * paraview_subtimestep;
+
+   //      //xref.Set(1.0, new_coords.GetTrueVector());
+   //      //xref.Add(-1.0, ref_coords.GetTrueVector());
+   //      xref.Set(1.0, x_gf.GetTrueVector());
+   //      ParContactProblem contact(prob, mortar_attr, nonmortar_attr, &new_coords, doublepass);
+   //      QPOptParContactProblem qpopt(&contact, xref);
+   //      int numconstr = contact.GetGlobalNumConstraints();
+   //      ParInteriorPointSolver optimizer(&qpopt);
+   //      optimizer.SetTol(optimizer_tol);
+   //      optimizer.SetMaxIter(optimizer_maxit);
+   //      optimizer.SetLinearSolver(linsolver);
+   //      optimizer.SetLinearSolveRelTol(linsolverrtol);
+   //      optimizer.SetLinearSolveAbsTol(linsolveratol);
+   //      optimizer.SetLinearSolveRelaxType(relax_type);
+   //      if (nocontact)
+   //      {
+   //         optimizer.EnableNoContactSolve();
+   //      }
+   //      if (elast)
+   //      {
+   //         optimizer.SetElasticityOptions(prob->GetFESpace());
+   //      }
+   //      // ParGridFunction x = prob->GetDisplacementGridFunction();
+   //      // x.SetTrueVector();
+   //      // Vector x0 = x.GetTrueVector();
+
+   //      x_gf.SetTrueVector();
+   //      Vector x0 = x_gf.GetTrueVector();
+   //      int ndofs = x0.Size();
+   //      Vector xf(ndofs); xf = 0.0;
+   //      optimizer.Mult(x0, xf);
+   //      QPConverged = optimizer.GetConverged();
+   //      MFEM_VERIFY(QPConverged, "IPM not converged on QP contact problem");
+   //      //optimizer.SaveLambda(i);
+   //      //optimizer.SaveZl(i);
+   //      Vector xf_copy(xf);
+   //      xf_copy+=x0;
+   //      double Einitial = contact.E(x0);
+   //      // double Efinal = contact.E(xf);
+   //      double Efinal = contact.E(xf_copy);
+   //      Array<int> & CGiterations = optimizer.GetCGIterNumbers();
+   //      int gndofs = prob->GetGlobalNumDofs();
+   //      //dgdu = contact.Ddg(xf_copy);
+   //      //std::ostringstream dgdu_file_name;
+   //      //dgdu_file_name << "Jacobians/J" << i;
+   //      //dgdu->Print(dgdu_file_name.str().c_str());
+   //      if (Mpi::Root())
+   //      {
+   //         mfem::out << endl;
+   //         mfem::out << " Initial Energy objective        = " << Einitial << endl;
+   //         mfem::out << " Final Energy objective          = " << Efinal << endl;
+   //         mfem::out << " Global number of dofs           = " << gndofs << endl;
+   //         mfem::out << " Global number of constraints    = " << numconstr << endl;
+   //         mfem::out << " Optimizer number of iterations  = " <<
+   //         optimizer.GetNumIterations() << endl;
+   //         if (linsolver == 2 || linsolver == 3 || linsolver == 4)
+   //         {
+   //            mfem::out << " CG iteration numbers            = " ;
+   //            CGiterations.Print(mfem::out, CGiterations.Size());
+   //         }
+   //         if (nocontact)
+   //         {
+   //            Array<int> & CGNoContactIterations = optimizer.GetCGNoContactIterNumbers();
+   //            mfem::out << " CG no Contact iteration numbers = " ;
+   //            CGNoContactIterations.Print(mfem::out, CGNoContactIterations.Size());
+   //         }
+   //         if (outputfiles)
+   //         {
+   //            ostringstream file_name;
+   //            file_name << "output/Testno-"<<testNo<<"-ref-"<<sref+pref << "-step-" << i; 
+   //            OutputData(file_name, Einitial, Efinal, gndofs,numconstr, optimizer.GetNumIterations(), CGiterations);
+   //         }
+   //      }
+
+   //      // Vector X_new(xf.GetData(),fes->GetTrueVSize());
+   //      // xnew.SetFromTrueDofs(X_new);
+   //      // x_gf = xnew;
+   //      x_gf.SetFromTrueDofs(xf);
+   //      // mfem::out << "x_gf norm = " << x_gf.Norml2() << endl;
+   //      // cin.get();
+   //      // pmesh->MoveNodes(xnew);
+   //      // pmesh_copy.MoveNodes(xnew);
+   //      // pmesh_copy.MoveNodes(xnew);
+   //      add(ref_coords,x_gf,new_coords);
+   //      // mfem::out << " ref_coords norm " << ref_coords.Norml2() << endl;
+   //      // mfem::out << " x_gf norm " << x_gf.Norml2() << endl;
+   //      // mfem::out << " new_coords norm " << new_coords.Norml2() << endl;
+   //      // pmesh_copy.SetNodes(new_coords);
+   //      pmesh_copy.SetNodes(new_coords);
+   //      xcopy_gf = x_gf;
+   //      // pmesh_copy.MoveNodes(x_gf);
+   //      // pmesh_copy.SetNodes(x_gf);
+   //      if (paraview && ((i+1) % paraview_plot_every == 0 ))
+   //      {
+   //         paraview_cycle += 1;
+   //         paraview_dc->SetCycle(paraview_cycle) ;
+   //         paraview_dc->SetTime(paraview_time);
+   //         paraview_dc->Save();
+   //      }
+
+   //      if (visualization)
+   //      {
+   //         sol_sock << "parallel " << num_procs << " " << myid << "\n"
+   //                  << "solution\n" << pmesh_copy << x_gf << flush;
+   //      
+   //         if (i == nsteps - 1 && j == SQPrepeat - 1)
+   //         {
+   //            pmesh->MoveNodes(x_gf);
+   //            char vishost[] = "localhost";
+   //            int  visport   = 19916;
+   //            socketstream sol_sock1(vishost, visport);
+   //            sol_sock1 << "parallel " << num_procs << " " << myid << "\n";
+   //            sol_sock1.precision(8);
+   //            sol_sock1 << "solution\n" << *pmesh << x_gf << flush;
+   //         }
+   //      }
+   //      if (i == nsteps - 1 && j == SQPrepeat) break;
+
+   //      prob->UpdateStep();
+   //      if (testNo == 6 )
+   //      {
+   //         double area_new = GetBdrArea(3,*pmesh);
+   //         if (myid == 0)
+   //         {
+   //            mfem::out << "New area = " << area_new << endl;
+   //         }
+   //      }
+   //   }
+   //}
+
+   delete prob;
+   delete pmesh;
+   delete mesh;
+   return 0;
+}

miniapps/meshing/pmesh-fitting.cpp

@@ -580,6 +580,22 @@ int main (int argc, char *argv[])
          }
       }
 
+      // Unify marker across processor boundary
+      surf_fit_mat_gf.ExchangeFaceNbrData();
+      {
+         GroupCommunicator &gcomm = surf_fit_mat_gf.ParFESpace()->GroupComm();
+         Array<real_t> gf_array(surf_fit_mat_gf.GetData(),
+                                surf_fit_mat_gf.Size());
+         gcomm.Reduce<real_t>(gf_array, GroupCommunicator::Max);
+         gcomm.Bcast(gf_array);
+      }
+      surf_fit_mat_gf.ExchangeFaceNbrData();
+
+      for (int i = 0; i < surf_fit_mat_gf.Size(); i++)
+      {
+         surf_fit_marker[i] = surf_fit_mat_gf(i) == 1.0;
+      }
+
       // Set AdaptivityEvaluators for transferring information from initial
       // mesh to current mesh as it moves during adaptivity.
       if (adapt_eval == 0)
@@ -591,11 +607,8 @@ int main (int argc, char *argv[])
       {
 #ifdef MFEM_USE_GSLIB
          adapt_surface = new InterpolatorFP;
-         if (surf_bg_mesh)
-         {
-            adapt_grad_surface = new InterpolatorFP;
-            adapt_hess_surface = new InterpolatorFP;
-         }
+         adapt_grad_surface = new InterpolatorFP;
+         adapt_hess_surface = new InterpolatorFP;
 #else
          MFEM_ABORT("MFEM is not built with GSLIB support!");
 #endif
@@ -605,7 +618,9 @@ int main (int argc, char *argv[])
       if (!surf_bg_mesh)
       {
          tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker,
-                                          surf_fit_coeff, *adapt_surface);
+                                          surf_fit_coeff, *adapt_surface,
+                                          adapt_grad_surface,
+                                          adapt_hess_surface);
       }
       else
       {
@@ -837,9 +852,13 @@ int main (int argc, char *argv[])
 
    if (surface_fit_const > 0.0)
    {
+      adapt_surface->ComputeAtNewPosition(x, surf_fit_gf0,
+                                          x.FESpace()->GetOrdering());
       if (visualization)
       {
-         socketstream vis2, vis3;
+         socketstream vis1, vis2, vis3;
+         common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0,
+                                "Level Set", 000, 400, 300, 300);
          common::VisualizeField(vis2, "localhost", 19916, mat,
                                 "Materials", 300, 400, 300, 300);
          common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,

miniapps/tribol/contact-patch-finite-difference-test.cpp

@@ -0,0 +1,453 @@
+// Copyright (c) 2010-2024, 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.
+//
+//                -----------------------------------------
+//                Tribol Miniapp: Mortar contact patch test
+//                -----------------------------------------
+//
+//
+// Command line options:
+//  - -r, --refine: number of uniform refinements of the mesh (default: 2)
+//
+
+#include "mfem.hpp"
+
+#include "axom/slic.hpp"
+
+#include "tribol/interface/tribol.hpp"
+#include "tribol/interface/mfem_tribol.hpp"
+
+// Define MPI_REAL_T
+#if defined(MFEM_USE_DOUBLE)
+#define MPI_REAL_T MPI_DOUBLE
+#else
+#error "Tribol requires MFEM built with double precision!"
+#endif
+
+
+using namespace mfem;
+
+class ContactObj
+{
+protected:
+   HypreParMatrix * Jacobian = nullptr;
+   mfem::Vector gap;
+   std::unique_ptr<mfem::BlockOperator> A_blk;
+   ParMesh * mesh = nullptr;
+   ParGridFunction * coords = nullptr;
+   std::set<int> mortar_attrs;
+   std::set<int> nonmortar_attrs;
+public:
+   ContactObj(ParMesh * mesh_, 
+	      const std::set<int> & mortar_attrs_,
+	      const std::set<int> & nonmortar_attrs_,
+	      ParGridFunction * coords_);
+   void GetGap(mfem::Vector & g) const;
+   mfem::HypreParMatrix * GetJacobian() const;
+   virtual ~ContactObj();
+};
+
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI
+   mfem::Mpi::Init();
+
+   // Initialize logging with axom::slic
+   axom::slic::SimpleLogger logger;
+   axom::slic::setIsRoot(mfem::Mpi::Root());
+
+   // Define command line options
+   int ref_levels = 2;   // number of times to uniformly refine the serial mesh
+   double u0shift = 0.0;
+   bool outputfiles = false;
+   // Parse command line options
+   mfem::OptionsParser args(argc, argv);
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&u0shift, "-u0shift", "--u0shift", "magnitude (inf norm) of random displacement where finite difference test is evaluated"); 
+   args.AddOption(&outputfiles, "-out", "--output", "-no-out",
+                  "--no-ouput",
+                  "Enable or disable ouput to files.");                  
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mfem::Mpi::Root())
+      {
+         args.PrintUsage(std::cout);
+      }
+      return EXIT_FAILURE;
+   }
+   if (mfem::Mpi::Root())
+   {
+      args.PrintOptions(std::cout);
+   }
+
+   // Fixed options
+   // two block mesh; bottom block = [0,1]^3 and top block = [0,1]x[0,1]x[0.99,1.99]
+   std::string mesh_file = "modified-two-hex.mesh";
+   // Problem dimension (NOTE: Tribol's mortar only works in 3D)
+   constexpr int dim = 3;
+   // FE polynomial degree (NOTE: only 1 works for now)
+   constexpr int order = 1;
+   // z=1 plane of bottom block (contact plane)
+   std::set<int> mortar_attrs({4});
+   // z=0.99 plane of top block (contact plane)
+   std::set<int> nonmortar_attrs({5});
+   // per-dimension sets of boundary attributes with homogeneous Dirichlet BCs.
+   // allows transverse deformation of the blocks while precluding rigid body
+   // rotations/translations.
+   std::vector<std::set<int>> fixed_attrs(dim);
+   fixed_attrs[0] = {1}; // x=0 plane of both blocks
+   fixed_attrs[1] = {2}; // y=0 plane of both blocks
+   fixed_attrs[2] = {3, 6}; // 3: z=0 plane of bottom block; 6: z=1.99 plane of top block
+
+   // Read the mesh, refine, and create a mfem::ParMesh
+   mfem::Mesh serial_mesh(mesh_file);
+   for (int i = 0; i < ref_levels; ++i)
+   {
+      serial_mesh.UniformRefinement();
+   }
+   mfem::ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   mfem::ParMesh mesh_copy(mesh);
+   
+   serial_mesh.Clear();
+
+   MFEM_ASSERT(dim == mesh.Dimension(),
+               "This miniapp must be run with the supplied two-hex.mesh file.");
+
+   // Create an H1 finite element space on the mesh for displacements/forces
+   mfem::H1_FECollection fec(order, dim);
+   mfem::ParFiniteElementSpace fespace(&mesh, &fec, dim);
+   auto n_displacement_dofs = fespace.GlobalTrueVSize();
+   if (mfem::Mpi::Root())
+   {
+      std::cout << "Number of displacement unknowns: " << n_displacement_dofs <<
+                std::endl;
+   }
+
+   // Create coordinate and displacement grid functions
+   mfem::ParGridFunction coords(&fespace);
+   mesh.SetNodalGridFunction(&coords);
+   mfem::ParGridFunction displacement(&fespace);
+   displacement = 0.0;
+
+   // Find true dofs with homogeneous Dirichlet BCs
+   mfem::Array<int> ess_tdof_list;
+   {
+      mfem::Array<int> ess_vdof_marker(fespace.GetVSize());
+      ess_vdof_marker = 0;
+      for (int i = 0; i < dim; ++i)
+      {
+         mfem::Array<int> ess_bdr(mesh.bdr_attributes.Max());
+         ess_bdr = 0;
+         for (auto xfixed_attr : fixed_attrs[i])
+         {
+            ess_bdr[xfixed_attr-1] = 1;
+         }
+         mfem::Array<int> new_ess_vdof_marker;
+         fespace.GetEssentialVDofs(ess_bdr, new_ess_vdof_marker, i);
+         for (int j = 0; j < new_ess_vdof_marker.Size(); ++j)
+         {
+            ess_vdof_marker[j] = ess_vdof_marker[j] || new_ess_vdof_marker[j];
+         }
+      }
+      mfem::Array<int> ess_tdof_marker;
+      fespace.GetRestrictionMatrix()->BooleanMult(ess_vdof_marker, ess_tdof_marker);
+      mfem::FiniteElementSpace::MarkerToList(ess_tdof_marker, ess_tdof_list);
+   }
+
+   // #1: Initialize Tribol contact library
+   tribol::initialize(dim, MPI_COMM_WORLD);
+
+   
+   /* Begin Tucker addition
+    * finite difference check of the gap function Jacobian at u = u0
+    * we evaluate the norm of the finite difference residual
+    * err(eps) = || (g(u0 + eps * udir) - g(u0)) / eps - J(u0) * udir ||_2
+    * which in the absence of finite-precision
+    * err(eps) = O(eps) when the gap is not linear
+    * err(eps) = 0, when the gap is linear
+    */
+   int dimU = fespace.GetTrueVSize(); 
+   Vector u0(dimU); u0 = 0.0; 
+   Vector u1(dimU); u1 = 0.0;
+   Vector udir(dimU); udir = 0.0; udir.Randomize(); udir *= 1.e-2;
+
+   Array<int> vdofs;
+   for (int i = 0; i < mesh.GetNBE(); i++)
+   {
+      const int attr = (mesh.GetBdrElement(i))->GetAttribute();
+      if (attr == 4)
+      {
+         fespace.GetBdrElementVDofs(i, vdofs);
+         for (int j = 0; j < vdofs.Size(); j++)
+         {
+             if (j / 4 == 2)
+             {      
+                u0(vdofs[j]) = -1.0 * u0shift;
+             }
+         }
+      }
+   }
+   
+   
+   ParGridFunction new_coords(&fespace);
+   mesh.GetNodes(new_coords);
+   
+   // evaluate the gap and gap Jacobian at u = u0   
+   u1.Set(1.0, u0);
+   displacement.SetFromTrueDofs(u1);
+   add(coords, displacement, new_coords);
+
+   ContactObj contact0(&mesh, mortar_attrs, nonmortar_attrs, &new_coords);
+   HypreParMatrix * J0 = contact0.GetJacobian();
+   int dimG = J0->Height();
+   
+   Vector g0(dimG); g0 = 0.0; contact0.GetGap(g0);
+   Vector g1(dimG); g1 = 0.0;
+
+   // finite difference residual
+   Vector fdres(dimG); fdres = 0.0;
+
+   // J0udir = J(u0) * udir
+   Vector J0udir(dimG); J0->Mult(udir, J0udir);
+
+   // output various configurations
+   // to visualize u = u0, u = u0 + eps * udir
+   // use linear adjustment for eps here  
+   std::ostringstream paraview_file_name;
+   paraview_file_name << "BlockConfigurations_ref_" << ref_levels << "shift" << u0shift;
+   ParaViewDataCollection * paraview_dc = new ParaViewDataCollection(paraview_file_name.str(), &mesh_copy);
+   paraview_dc->SetPrefixPath("ParaView");
+   paraview_dc->SetLevelsOfDetail(1);
+   paraview_dc->SetDataFormat(VTKFormat::BINARY);
+   paraview_dc->SetHighOrderOutput(true);
+   paraview_dc->SetCycle(0);
+   paraview_dc->SetTime(double(0));
+   paraview_dc->Save();
+   
+   std::ofstream fdepsStream;
+   std::ostringstream fdeps_file_name;
+   fdeps_file_name << "data/fdeps.dat";
+
+   std::ofstream fderrStream;
+   std::ostringstream fderr_file_name;
+   fderr_file_name << "data/fderr.dat";
+   
+   // write new configuration (reference coordinates + displacement u0) to file 
+   u1.Set(1.0, u0);
+   displacement.SetFromTrueDofs(u1);
+   add(coords, displacement, new_coords);
+   
+   Vector config(u0.Size()); config = 0.0;
+   new_coords.GetTrueDofs(config);
+   
+
+   if (mfem::Mpi::Root() && outputfiles)
+   {
+      fdepsStream.open(fdeps_file_name.str(), std::ios::out | std::ios::trunc);
+      fderrStream.open(fderr_file_name.str(), std::ios::out | std::ios::trunc);
+   }
+
+   double eps = 1.0;
+   int   neps = 40;
+   for (int i = 0; i < neps; i++) // eps_min = 0.5^(39) \approx 10^(-12)
+   {
+      // compute g1 = g(u1), u1 = u0 + eps * udir
+      u1.Set(1.0, u0);
+      u1.Add(eps, udir);
+      displacement.SetFromTrueDofs(u1);
+      add(coords, displacement, new_coords);
+      ContactObj contact1(&mesh, mortar_attrs, nonmortar_attrs, &new_coords);
+      contact1.GetGap(g1);
+
+      // determine finite difference residual: fdres = (g1 - g0) / eps - J0 * udir
+      fdres.Set(1. / eps, g1);
+      fdres.Add(-1. / eps, g0);
+      fdres.Add(-1, J0udir);
+      double fderr_l2norm = GlobalLpNorm(2, fdres.Norml2(), MPI_COMM_WORLD);
+      double udir_l2norm = GlobalLpNorm(2, udir.Norml2(), MPI_COMM_WORLD);
+
+      if (mfem::Mpi::Root())
+      {
+         std::cout << "--------------------------------------------\n\n";
+         std::cout << "||(g(u0 + eps * udir) - g(u0)) / eps - J(u0) * udir|| = " << fderr_l2norm << ", eps = " << eps << "\n\n";
+	 std::cout << "||(g(u0 + eps * udir) - g(u0)) / eps - J(u0) * udir||_2 / ||udir||_2 = " << fderr_l2norm / udir_l2norm << std::endl;
+      }
+      if (mfem::Mpi::Root() && outputfiles)
+      {
+	 fdepsStream << eps << std::endl;
+	 fderrStream << fderr_l2norm << std::endl;
+      }
+      eps /= 2.0;
+   }
+   if (mfem::Mpi::Root() && outputfiles)
+   {
+      fdepsStream.close();
+      fderrStream.close();
+   }
+
+
+   /* What follows we linearly modify epsilon
+    * output the gap, in order to check for discontinuities
+    * and also output the various states u0 + eps * udir to file
+    * in order to visualize the mesh configurations    *
+    * */
+   eps  = 1.0;
+   neps = 100;
+   double deps =  eps / ((double) neps);
+   std::ofstream epsStream;
+   std::ostringstream eps_file_name;
+   eps_file_name << "data/eps_ref_" << ref_levels << ".dat";
+   
+   std::ofstream gapStream;
+   std::ostringstream gap_file_name;
+   gap_file_name << "data/gap_ref_" << ref_levels << ".dat";
+   
+   if (mfem::Mpi::Root() && outputfiles)
+   {
+      epsStream.open(eps_file_name.str(), std::ios::out | std::ios::trunc);
+      gapStream.open(gap_file_name.str(), std::ios::out | std::ios::trunc);
+   }
+
+
+
+
+   for (int i = 0; i < neps; i++)
+   {
+      // compute g1 = g(u1), u1 = u0 + eps * udir
+      u1.Set(1.0, u0);
+      u1.Add(eps, udir);
+      displacement.SetFromTrueDofs(u1);
+      add(coords, displacement, new_coords);
+      ContactObj contact1(&mesh, mortar_attrs, nonmortar_attrs, &new_coords);
+      contact1.GetGap(g1);
+      double gap_l2norm = GlobalLpNorm(2, g1.Norml2(), MPI_COMM_WORLD);
+      if (mfem::Mpi::Root() && outputfiles)
+      {
+         epsStream << eps << std::endl;
+         gapStream << g1.Norml2() << std::endl;
+      } 
+      
+      // update mesh according to u1 and write to Paraview for visualization
+      mesh_copy.SetNodes(new_coords);
+      paraview_dc->SetCycle(i+1) ;
+      paraview_dc->SetTime((double) (i+1));
+      paraview_dc->Save();
+      
+      // linear update to eps: eps = eps - deps
+      eps -= deps;
+   }
+
+   if (mfem::Mpi::Root() && outputfiles)
+   {
+      epsStream.close();
+      gapStream.close();
+   }
+
+   // #7: Tribol cleanup: deletes coupling schemes and clears associated memory
+   tribol::finalize();
+
+   return 0;
+}
+
+ContactObj::ContactObj(ParMesh * mesh_, const std::set<int> & mortar_attrs_,
+	      const std::set<int> & nonmortar_attrs_,
+	      ParGridFunction * coords_) : 
+	mesh(mesh_), mortar_attrs(mortar_attrs_), 
+	nonmortar_attrs(nonmortar_attrs_),
+	coords(coords_)
+{
+   // #2: Create a Tribol coupling scheme: defines contact surfaces and enforcement
+   int coupling_scheme_id = 0;
+   // NOTE: While there is a single mfem ParMesh for this problem, Tribol
+   // defines a mortar and a nonmortar contact mesh, each with a unique mesh ID.
+   // The Tribol mesh IDs for each contact surface are defined here.
+   int mesh1_id = 0;
+   int mesh2_id = 1;
+   tribol::registerMfemCouplingScheme(
+      coupling_scheme_id, mesh1_id, mesh2_id,
+      *mesh, *coords, mortar_attrs, nonmortar_attrs,
+      tribol::SURFACE_TO_SURFACE,
+      tribol::NO_CASE,
+      tribol::SINGLE_MORTAR,
+      tribol::FRICTIONLESS,
+      tribol::LAGRANGE_MULTIPLIER,
+      tribol::BINNING_GRID
+   );
+   
+   // #3: Set additional options/access pressure grid function on contact surfaces
+   // Access Tribol's pressure grid function (on the contact surface). The
+   // pressure ParGridFunction is created upon calling
+   // registerMfemCouplingScheme(). It's lifetime coincides with the lifetime of
+   // the coupling scheme, so the host code can reference and update it as
+   // needed.
+   auto& pressure = tribol::getMfemPressure(coupling_scheme_id);
+
+   // Set Tribol options for Lagrange multiplier enforcement
+   tribol::setLagrangeMultiplierOptions(
+      coupling_scheme_id,
+      tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
+   );
+
+   // #4: Update contact mesh decomposition so the on-rank Tribol meshes
+   // coincide with the current configuration of the mesh. This must be called
+   // before tribol::update().
+   tribol::updateMfemParallelDecomposition();
+
+   // #5: Update contact gaps, forces, and tangent stiffness contributions
+   int cycle = 1;   // pseudo cycle
+   mfem::real_t t = 1.0;  // pseudo time
+   mfem::real_t dt = 1.0; // pseudo dt
+   tribol::update(cycle, t, dt);
+
+   
+   // #6a: Return contact contribution to the tangent stiffness matrix as a
+   // block operator. See documentation for getMfemBlockJacobian() for block
+   // definitions.
+   //auto A_blk = tribol::getMfemBlockJacobian(coupling_scheme_id);
+   A_blk = tribol::getMfemBlockJacobian(coupling_scheme_id);
+   Jacobian = (HypreParMatrix *)(& A_blk->GetBlock(1, 0));
+
+   
+    mfem::BlockVector B_blk(A_blk->RowOffsets());
+    B_blk = 0.0;
+
+   // Fill with initial nodal gaps.
+   // Note forces from contact are currently zero since pressure is zero prior
+   // to first solve.
+   mfem::Vector gap_temp;
+   // #6b: Return computed gap constraints on the contact surfaces
+   tribol::getMfemGap(coupling_scheme_id, gap_temp); // gap on ldofs
+   auto& P_submesh = *pressure.ParFESpace()->GetProlongationMatrix();
+   //auto& gap_true = B_blk.GetBlock(1); // gap tdof vectorParFESpace()
+   // gap is a dual vector, so (gap tdof vector) = P^T * (gap ldof vector)
+   gap.SetSize(P_submesh.Width()); gap = 0.0;
+   
+   P_submesh.MultTranspose(gap_temp, gap);
+}
+
+void ContactObj::GetGap(mfem::Vector & g) const
+{
+   g.SetSize(gap.Size());
+   g.Set(1.0, gap);
+}
+
+
+mfem::HypreParMatrix * ContactObj::GetJacobian() const
+{
+   return Jacobian;
+}
+
+ContactObj::~ContactObj()
+{
+}

miniapps/tribol/makefile

@@ -22,7 +22,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 ifeq ($(MFEM_USE_TRIBOL)$(MFEM_USE_MPI),YESYES)
-   MINIAPPS = contact-patch-test
+   MINIAPPS = contact-patch-test contact-patch-finite-difference-test
 else
    MINIAPPS =
 endif

miniapps/tribol/modified-two-hex.mesh

@@ -0,0 +1,57 @@
+MFEM mesh v1.0
+
+# two unit cubes occupying [0,1]^3 and [0,1]x[0,1]x[0.99,1.99]
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+2
+1 5 0 1 3 2 4 5 7 6
+1 5 8 9 11 10 12 13 15 14
+
+boundary
+12
+3 3 2 3 1 0
+2 3 0 1 5 4
+7 3 3 2 6 7
+1 3 2 0 4 6
+4 3 4 5 7 6
+7 3 1 3 7 5
+5 3 10 11 9 8
+2 3 8 9 13 12
+7 3 11 10 14 15
+1 3 10 8 12 14
+6 3 12 13 15 14
+7 3 9 11 15 13
+
+vertices
+16
+3
+0.25 0.25 0
+0.75 0.25 0
+0.25 0.75 0
+0.75 0.75 0
+0.25 0.25 1
+0.75 0.25 1
+0.25 0.75 1
+0.75 0.75 1
+0 0 1.00
+1 0 1.00
+0 1 1.00
+1 1 1.00
+0 0 2.00
+1 0 2.00
+0 1 2.00
+1 1 2.00

miniapps/tribol/two-hex.mesh

@@ -47,11 +47,11 @@ vertices
 1 0 1
 0 1 1
 1 1 1
-0 0 0.99
-1 0 0.99
-0 1 0.99
-1 1 0.99
+0 0 1.01
+1 0 1.01
+0 1 1.01 
+1 1 1.01
 0 0 1.99
 1 0 1.99
 0 1 1.99
-1 1 1.99
+1 1 1.99

tests/unit/data/quad_attribute.vtu

@@ -0,0 +1,58 @@
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="16" NumberOfCells="9">
+      <Points>
+        <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="0" RangeMax="4.242640687119285">
+          0 0 0 1 0 0
+          2 0 0 3 0 0
+          0 1 0 1 1 0
+          2 1 0 3 1 0
+          0 2 0 1 2 0
+          2 2 0 3 2 0
+          0 3 0 1 3 0
+          2 3 0 3 3 0
+          <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              0
+            </Value>
+            <Value index="1">
+              4.2426406871
+            </Value>
+          </InformationKey>
+          <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              0
+            </Value>
+            <Value index="1">
+              4.2426406871
+            </Value>
+          </InformationKey>
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="15">
+          0 1 5 4 1 2
+          6 5 2 3 7 6
+          4 5 9 8 5 6
+          10 9 6 7 11 10
+          8 9 13 12 9 10
+          14 13 10 11 15 14
+        </DataArray>
+        <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="4" RangeMax="36">
+          4 8 12 16 20 24
+          28 32 36
+        </DataArray>
+        <DataArray type="UInt8" Name="types" format="ascii" RangeMin="9" RangeMax="9">
+          9 9 9 9 9 9
+          9 9 9
+        </DataArray>
+      </Cells>
+      <CellData Scalars="attribute">
+        <DataArray type="UInt8" Name="attribute" format="ascii">
+          1 2 3 4 5 6
+          7 8 9
+        </DataArray>
+      </CellData>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>

tests/unit/data/quad_material_attribute.vtu

@@ -0,0 +1,64 @@
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="16" NumberOfCells="9">
+      <Points>
+        <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="0" RangeMax="4.242640687119285">
+          0 0 0 1 0 0
+          2 0 0 3 0 0
+          0 1 0 1 1 0
+          2 1 0 3 1 0
+          0 2 0 1 2 0
+          2 2 0 3 2 0
+          0 3 0 1 3 0
+          2 3 0 3 3 0
+          <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              0
+            </Value>
+            <Value index="1">
+              4.2426406871
+            </Value>
+          </InformationKey>
+          <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              0
+            </Value>
+            <Value index="1">
+              4.2426406871
+            </Value>
+          </InformationKey>
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="15">
+          0 1 5 4 1 2
+          6 5 2 3 7 6
+          4 5 9 8 5 6
+          10 9 6 7 11 10
+          8 9 13 12 9 10
+          14 13 10 11 15 14
+        </DataArray>
+        <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="4" RangeMax="36">
+          4 8 12 16 20 24
+          28 32 36
+        </DataArray>
+        <DataArray type="UInt8" Name="types" format="ascii" RangeMin="9" RangeMax="9">
+          9 9 9 9 9 9
+          9 9 9
+        </DataArray>
+      </Cells>
+      <CellData Scalars="attribute">
+        <DataArray type="UInt8" Name="attribute" format="ascii">
+          10 11 12 13 14 15
+          16 17 18
+        </DataArray>
+      </CellData>
+      <CellData Scalars="material">
+        <DataArray type="UInt8" Name="material" format="ascii">
+          1 2 3 4 5 6
+          7 8 9
+        </DataArray>
+      </CellData>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>

tests/unit/fem/test_coefficient.cpp

@@ -304,6 +304,66 @@ TEST_CASE("Piecewise Matrix Coefficient", "[Coefficient]")
    }
 }
 
+TEST_CASE("MatrixArrayVectorCoefficient", "[Coefficient]")
+{
+   Vector V1(2), V2(2);
+   V1(0) = 0.0; V1(1) = 1.0;
+   V2(0) = 2.0; V2(1) = 3.0;
+   VectorConstantCoefficient Coef1(V1), Coef2(V2);
+
+   IsoparametricTransformation T;
+   IntegrationPoint ip;
+
+   MatrixArrayVectorCoefficient mavc(2);
+   Vector V(2);
+
+   // Verify zeros for unset rows
+   int row = 0;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(0.0));
+
+   row = 1;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(0.0));
+
+   DenseMatrix K(2);
+   mavc.Eval(K, T, ip);
+   REQUIRE(K(0,0) == MFEM_Approx(0.0));
+   REQUIRE(K(0,1) == MFEM_Approx(0.0));
+   REQUIRE(K(1,0) == MFEM_Approx(0.0));
+   REQUIRE(K(1,1) == MFEM_Approx(0.0));
+
+   // Test setting individual rows
+   row = 0;
+   mavc.Set(row, &Coef1, false);
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(1.0));
+   row = 1;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(0.0));
+
+   mavc.Set(row, &Coef2, false);
+   row = 0;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(0.0));
+   REQUIRE(V(1) == MFEM_Approx(1.0));
+   row = 1;
+   mavc.Eval(row, V, T, ip);
+   REQUIRE(V(0) == MFEM_Approx(2.0));
+   REQUIRE(V(1) == MFEM_Approx(3.0));
+
+   mavc.Eval(K, T, ip);
+   REQUIRE(K(0,0) == MFEM_Approx(0.0));
+   REQUIRE(K(0,1) == MFEM_Approx(1.0));
+   REQUIRE(K(1,0) == MFEM_Approx(2.0));
+   REQUIRE(K(1,1) == MFEM_Approx(3.0));
+
+}
+
 TEST_CASE("Symmetric Matrix Coefficient", "[Coefficient]")
 {
    int d = 3;

tests/unit/linalg/test_matrix_dense.cpp

@@ -623,6 +623,67 @@ TEST_CASE("MatrixInverse", "[DenseMatrix]")
    }
 
 }
+
+TEST_CASE("Exponential", "[DenseMatrix]")
+{
+   // case 1
+   DenseMatrix A(2,2);
+   A(0,0) = 5.0;
+   A(0,1) = 3.0;
+   A(1,0) = 0.0;
+   A(1,1) = 5.0;
+   A.Exponential();
+
+   DenseMatrix expA(2,2);
+   expA(0,0) = std::exp(5.0);
+   expA(0,1) = 3.0 * std::exp(5.0);
+   expA(1,0) = 0.0;
+   expA(1,1) = std::exp(5.0);
+
+   A.Print();
+   expA.Print();
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+
+   // case 2
+   A(0,0) = 3.0;
+   A(0,1) = 5.0;
+   A(1,0) = 4.0;
+   A(1,1) = 2.0;
+   A.Exponential();
+
+   expA(0,0) = 4.0 / (9.0 * std::exp(2.0)) + (5.0 * std::exp(7.0)) / 9.0;
+   expA(0,1) = (5.0 * std::exp(7.0)) / 9.0 - 5.0 / (9.0 * std::exp(2.0));
+   expA(1,0) = (4.0 * std::exp(7.0)) / 9.0 - 4.0 / (9.0 * std::exp(2.0));
+   expA(1,1) = 5.0 / (9.0 * std::exp(2.0)) + (4.0 * std::exp(7.0)) / 9.0;
+
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+
+   // case 3
+   A(0,0) = 10.0;
+   A(0,1) = 2.0;
+   A(1,0) = -2.0;
+   A(1,1) = 8.0;
+   A.Exponential();
+
+   expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(3.0)
+                                + std::cos(std::sqrt(3.0)));
+   expA(0,1) = 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
+   expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
+   expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0))
+                                - std::sin(std::sqrt(3.0)) / std::sqrt(3.0));
+
+   REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
+   REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
+   REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
+   REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
+}
+
 #ifdef MFEM_USE_LAPACK
 
 enum class TestCase { GenEigSPD, GenEigGE, SVD};

tests/unit/mesh/test_ncmesh.cpp

@@ -2811,4 +2811,82 @@ TEST_CASE("RP=I", "[NCMesh]")
    }
 }
 
+
+TEST_CASE("InternalBoundaryProjectBdrCoefficient", "[NCMesh]")
+{
+   auto test_project_H1 = [](Mesh &mesh, int order, double coef)
+   {
+      MFEM_ASSERT(std::abs(coef) > 0,
+                  "Non zero coef value required for meaningful test.");
+      H1_FECollection fe_collection(order, mesh.SpaceDimension());
+      FiniteElementSpace fe_space(&mesh, &fe_collection);
+      GridFunction x(&fe_space);
+      x = -coef;
+      ConstantCoefficient c(coef);
+
+      // Check projecting on the internal face sets essential dof.
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 0;
+      ess_bdr.Last() = 1; // internal boundary
+      x.ProjectBdrCoefficient(c, ess_bdr);
+
+      Array<int> ess_vdofs_list, ess_vdofs_marker;
+      fe_space.GetEssentialVDofs(ess_bdr, ess_vdofs_marker);
+      fe_space.MarkerToList(ess_vdofs_marker, ess_vdofs_list);
+      for (auto ess_dof : ess_vdofs_list)
+      {
+         CHECK(x[ess_dof] == Approx(coef).epsilon(1e-8));
+      }
+
+      int iess = 0;
+      for (int i = 0; i < x.Size(); i++)
+      {
+         if (iess < ess_vdofs_list.Size() && i == ess_vdofs_list[iess])
+         {
+            iess++;
+            continue;
+         }
+         CHECK(x[i] == Approx(-coef).epsilon(1e-8));
+      }
+
+   };
+
+   auto OneSidedNCRefine = [](Mesh &mesh)
+   {
+      // Pick one element attached to the new boundary attribute and refine.
+      const auto interface_attr = mesh.bdr_attributes.Max();
+      Array<int> el_to_ref;
+      for (int nbe = 0; nbe < mesh.GetNBE(); nbe++)
+      {
+         if (mesh.GetBdrAttribute(nbe) == interface_attr)
+         {
+            int f, o, e1, e2;
+            mesh.GetBdrElementFace(nbe, &f, &o);
+            mesh.GetFaceElements(f, &e1, &e2);
+            el_to_ref.Append(e1);
+         }
+      }
+      mesh.GeneralRefinement(el_to_ref);
+      return;
+   };
+
+   SECTION("Hex")
+   {
+      auto smesh = DividingPlaneMesh(false, true);
+      smesh.EnsureNCMesh(true);
+      OneSidedNCRefine(smesh);
+      test_project_H1(smesh, 2, 0.25);
+   }
+
+   SECTION("Tet")
+   {
+      auto smesh = DividingPlaneMesh(true, true);
+      smesh.EnsureNCMesh(true);
+      OneSidedNCRefine(smesh);
+      test_project_H1(smesh, 3, 0.25);
+   }
+}
+
+
+
 } // namespace mfem

tests/unit/mesh/test_vtu.cpp

@@ -29,16 +29,35 @@ TEST_CASE("VTU XML Reader", "[Mesh][VTU][XML]")
       "quad_binary_compress.vtu"
    });
 #endif
-   for (const std::string &fname : mesh_filenames)
+
+   const auto fname = GENERATE_COPY(from_range(mesh_filenames));
+
+   Mesh mesh = Mesh::LoadFromFile("data/" + fname);
+   REQUIRE(mesh.Dimension() == 2);
+   REQUIRE(mesh.GetNE() == 9);
+   REQUIRE(mesh.GetNV() == 16);
+   REQUIRE(mesh.HasGeometry(Geometry::POINT));
+   REQUIRE(mesh.HasGeometry(Geometry::SEGMENT));
+   REQUIRE(mesh.HasGeometry(Geometry::SQUARE));
+   REQUIRE(mesh.GetNumGeometries(2) == 1);
+}
+
+TEST_CASE("VTU Attributes", "[VTU][XML]")
+{
+   // quad_attribute.vtu contains the attributes in a cell data array named
+   // "attribute"
+   Mesh mesh_1 = Mesh::LoadFromFile("data/quad_attribute.vtu");
+   // quad_material_attribute.vtu has cell data arrays named "material" and
+   // "attribute". The one named "material" should take precedence.
+   Mesh mesh_2 = Mesh::LoadFromFile("data/quad_material_attribute.vtu");
+
+   REQUIRE(mesh_1.GetNE() == 9);
+   REQUIRE(mesh_2.GetNE() == 9);
+
+   for (int i = 0; i < mesh_1.GetNE(); ++i)
    {
-      Mesh mesh = Mesh::LoadFromFile(("data/" + fname).c_str());
-      REQUIRE(mesh.Dimension() == 2);
-      REQUIRE(mesh.GetNE() == 9);
-      REQUIRE(mesh.GetNV() == 16);
-      REQUIRE(mesh.HasGeometry(Geometry::POINT));
-      REQUIRE(mesh.HasGeometry(Geometry::SEGMENT));
-      REQUIRE(mesh.HasGeometry(Geometry::SQUARE));
-      REQUIRE(mesh.GetNumGeometries(2) == 1);
+      REQUIRE(mesh_1.GetAttribute(i) == i+1);
+      REQUIRE(mesh_2.GetAttribute(i) == i+1);
    }
 }
 
@@ -52,8 +71,7 @@ TEST_CASE("VTU XML Compressed Blocks", "[VTU][XML][MFEMData]")
                       "bracket_inline_compressed.vtu"
                    );
 
-   std::string mesh_path = mfem_data_dir + "/vtk/" + filename;
-   Mesh mesh = Mesh::LoadFromFile(mesh_path.c_str());
+   Mesh mesh = Mesh::LoadFromFile(mfem_data_dir + "/vtk/" + filename);
 
    REQUIRE(mesh.Dimension() == 3);
    REQUIRE(mesh.GetNE() == 206208);