Few changes

.gitignore

@@ -57,8 +57,6 @@ 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,9 +13,6 @@
 # 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,10 +9,6 @@
 # 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,8 +35,9 @@ 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 the pci queue on lassen
-# to speed-up the allocation.
+# 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.
 .build_and_test_on_lassen:
   extends: [.on_lassen]
   stage: build_and_test
@@ -44,5 +45,5 @@ variables:
     - echo ${MFEM_DATA_DIR}
     - echo ${SPEC}
     # Next script uses 'THREADS': leaving it empty --> it uses 'make all -j'
-    - lalloc 1 -W 45 -q pci --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - lalloc 1 -W 45 -q pdebug --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}" ) --reservation=ci -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -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 --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  salloc --nodes=1 -p pdebug ../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 pci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  lalloc 1 -q pdebug ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
   echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
   exit 1

CHANGELOG

@@ -11,9 +11,6 @@
 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,9 +114,7 @@ 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 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
+ * - <a class="el" href="ex39p_8cpp_source.html">Example 39</a>: parallel named mesh attributes
  *
  * <H4>AmgX Examples</H4>
  * - Variants of Examples

examples/CMakeLists.txt

@@ -45,7 +45,6 @@ list(APPEND ALL_EXE_SRCS
   ex37.cpp
   ex38.cpp
   ex39.cpp
-  ex40.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -88,7 +87,6 @@ if (MFEM_USE_MPI)
     ex36p.cpp
     ex37p.cpp
     ex39p.cpp
-    ex40p.cpp
   )
 endif()
 
@@ -148,10 +146,10 @@ if (MFEM_ENABLE_TESTING)
   # Add CUDA/HIP tests.
   set(DEVICE_EXAMPLES
       # serial examples with device support:
-      ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
+      ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26 ex34
       # parallel examples with device support:
-      ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p ex22p ex24p ex25p
-      ex26p ex34p ex35p)
+      ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p ex24p ex25p ex26p
+      ex34p ex35p)
   set(MFEM_TEST_DEVICE)
   if (MFEM_USE_CUDA)
     set(MFEM_TEST_DEVICE "cuda")
@@ -161,11 +159,6 @@ 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

@@ -1,374 +0,0 @@
-//                                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

@@ -1,436 +0,0 @@
-//                       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 ex40
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39
 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 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
+   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
 
 ifeq ($(MFEM_USE_LAPACK),YES)
    SEQ_EXAMPLES += ex38
@@ -138,14 +138,6 @@ 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 u, \nabla \cdot v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot  \nabla \cdot 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 v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla \cdot 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 v)$ in 2D or 3D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla \cdot 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,75 +924,6 @@ 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); }
@@ -1110,27 +1041,6 @@ 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),
@@ -1416,30 +1326,6 @@ 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,46 +1334,6 @@ 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. */
@@ -1807,31 +1767,6 @@ 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
 {
@@ -2183,7 +2118,7 @@ public:
                      const IntegrationPoint &ip);
 };
 
-/// Matrix coefficient defined as the transpose of a matrix coefficient
+/// Matrix coefficient defined as the transpose a matrix coefficient
 class TransposeMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2206,7 +2141,7 @@ public:
                      const IntegrationPoint &ip);
 };
 
-/// Matrix coefficient defined as the inverse of a matrix coefficient.
+/// Matrix coefficient defined as the inverse a matrix coefficient.
 class InverseMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2229,29 +2164,6 @@ 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) const
+void GridFunction::AccumulateAndCountDerivativeValues(int comp, int der_comp,
+                                                      GridFunction &der,
+                                                      Array<int> &zones_per_dof)
 {
    FiniteElementSpace * der_fes = der.FESpace();
    ElementTransformation * transf;
@@ -1374,8 +1374,7 @@ void GridFunction::AccumulateAndCountDerivativeValues(
    }
 }
 
-void GridFunction::GetDerivative(int comp, int der_comp,
-                                 GridFunction &der) const
+void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der)
 {
    Array<int> overlap;
    AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);
@@ -2062,37 +2061,41 @@ 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;
 
-   const int vdim = fes->GetVDim();
+   vdim = fes->GetVDim();
+
    HostReadWrite();
 
-   for (int i = 0; i < fes->GetNBE(); i++)
+   for (i = 0; i < fes->GetNBE(); i++)
    {
       if (attr[fes->GetBdrAttribute(i) - 1] == 0) { continue; }
 
-      const FiniteElement *fe = fes->GetBE(i);
-      const int fdof = fe->GetDof();
-      ElementTransformation *transf = fes->GetBdrElementTransformation(i);
+      fe = fes->GetBE(i);
+      fdof = fe->GetDof();
+      transf = fes->GetBdrElementTransformation(i);
       const IntegrationRule &ir = fe->GetNodes();
       fes->GetBdrElementVDofs(i, vdofs);
 
-      for (int j = 0; j < fdof; j++)
+      for (j = 0; j < fdof; j++)
       {
          const IntegrationPoint &ip = ir.IntPoint(j);
          transf->SetIntPoint(&ip);
          if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
-         for (int d = 0; d < vdim; d++)
+         for (d = 0; d < vdim; d++)
          {
             if (!vcoeff && !coeff[d]) { continue; }
 
-            real_t val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
-            int ind = vdofs[fdof*d+j];
-            if ( ind < 0 )
+            val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
+            if ( (ind = vdofs[fdof*d+j]) < 0 )
             {
                val = -val, ind = -1-ind;
             }
@@ -2114,11 +2117,10 @@ 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, 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))
+   // (depend on) internal faces/elements. We use the virtual method
+   // GetBoundaryClosure from NCMesh to resolve the dependencies.
+
+   if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
    {
       Vector vals;
       Mesh *mesh = fes->GetMesh();
@@ -2126,19 +2128,26 @@ void GridFunction::AccumulateAndCountBdrValues(
       Array<int> bdr_edges, bdr_vertices, bdr_faces;
       ncmesh->GetBoundaryClosure(attr, bdr_vertices, bdr_edges, bdr_faces);
 
-      auto mark_dofs = [&](ElementTransformation &transf, const FiniteElement &fe)
+      for (i = 0; i < bdr_edges.Size(); i++)
       {
+         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 (int d = 0; d < vdim; d++)
+            vals.SetSize(fe->GetDof());
+            for (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++)
                {
-                  const int ind = vdofs[d*vals.Size()+k];
+                  ind = vdofs[d*vals.Size()+k];
                   if (++values_counter[ind] == 1)
                   {
                      (*this)(ind) = vals(k);
@@ -2152,11 +2161,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++)
             {
-               const int ind = vdofs[k];
+               ind = vdofs[k];
                if (++values_counter[ind] == 1)
                {
                   (*this)(ind) = vals(k);
@@ -2167,26 +2176,6 @@ 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);
       }
    }
 }
@@ -2239,37 +2228,26 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
       accumulate_dofs(dofs, lvec, *this, values_counter);
    }
 
-   if (fes->Nonconforming() && (fes->GetMesh()->Dimension() == 2 ||
-                                fes->GetMesh()->Dimension() == 3))
+   if (fes->Nonconforming() && 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 (auto edge : bdr_edges)
+      for (int i = 0; i < bdr_edges.Size(); i++)
       {
+         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) const;
+   void GetDerivative(int comp, int der_comp, GridFunction &der);
 
    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) const;
+                                           Array<int> &zones_per_dof);
 
    void AccumulateAndCountBdrValues(Coefficient *coeff[],
                                     VectorCoefficient *vcoeff,

fem/gslib.cpp

@@ -1393,10 +1393,8 @@ 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) const
+                                    ParGridFunction &der)
 {
    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) const;
+   void GetDerivative(int comp, int der_comp, ParGridFunction &der);
 
    /** 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,15 +2949,6 @@ 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);
@@ -2996,24 +2987,12 @@ void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &pos,
 void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
                                            const Array<bool> &smarker,
                                            Coefficient &coeff,
-                                           AdaptivityEvaluator &ae,
-                                           AdaptivityEvaluator *aegrad,
-                                           AdaptivityEvaluator *aehess)
+                                           AdaptivityEvaluator &ae)
 {
    // 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);
@@ -3021,80 +3000,11 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
    surf_fit_coeff = &coeff;
    surf_fit_eval = &ae;
 
-   surf_fit_eval->SetParMetaInfo(*pmesh, *s0.ParFESpace());
+   surf_fit_eval->SetParMetaInfo(*s0.ParFESpace()->GetParMesh(),
+                                 *s0.ParFESpace());
    surf_fit_eval->SetInitialField
    (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
-
-   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;
+   surf_fit_gf_bg = false;
 }
 
 void TMOP_Integrator::EnableSurfaceFittingFromSource(
@@ -3112,10 +3022,12 @@ 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
@@ -3132,11 +3044,11 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
    delete surf_fit_grad;
    surf_fit_grad = new GridFunction(s0_grad);
    *surf_fit_grad = 0.0;
-   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
+   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
    (*s_bg_grad.FESpace()->GetMesh()->GetNodes(), s_bg_grad);
 
    // Setup for Hessian on background mesh
@@ -3147,9 +3059,9 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
    delete surf_fit_hess;
    surf_fit_hess = new GridFunction(s0_hess);
    *surf_fit_hess = 0.0;
-   surf_fit_eval_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(),
-                                      *s_bg_hess.ParFESpace());
-   surf_fit_eval_hess->SetInitialField
+   surf_fit_eval_bg_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(),
+                                         *s_bg_hess.ParFESpace());
+   surf_fit_eval_bg_hess->SetInitialField
    (*s_bg_hess.FESpace()->GetMesh()->GetNodes(), s_bg_hess);
 
    // Count number of zones that share each of the DOFs
@@ -3954,7 +3866,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)
+   if (surf_fit_gf || surf_fit_gf_bg)
    {
       surf_fit_gf->GetSubVector(vdofs, sigma_e);
 
@@ -3962,7 +3874,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_grad)
+      if (surf_fit_gf_bg)
       {
          surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
          surf_fit_grad->GetSubVector(dofs, grad_ptr);
@@ -4036,7 +3948,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)
+   if (surf_fit_gf || surf_fit_gf_bg)
    {
       surf_fit_gf->GetSubVector(vdofs, sigma_e);
 
@@ -4044,7 +3956,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_grad)
+      if (surf_fit_gf_bg)
       {
          surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
          surf_fit_grad->GetSubVector(dofs, grad_ptr);
@@ -4058,7 +3970,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_hess)
+      if (surf_fit_gf_bg)
       {
          surf_fit_hess->FESpace()->GetElementVDofs(el_id, dofs);
          surf_fit_hess->GetSubVector(dofs, hess_ptr);
@@ -4085,7 +3997,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)
+      if (surf_fit_gf || surf_fit_gf_bg)
       {
          Vector gg_ptr(surf_fit_hess_s.GetData(), dim * dim);
          surf_fit_hess_e.GetRow(s, gg_ptr);
@@ -4467,130 +4379,6 @@ 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)
@@ -4621,11 +4409,112 @@ UpdateAfterMeshPositionChange(const Vector &x_new,
       adapt_lim_eval->ComputeAtNewPosition(x_new, *adapt_lim_gf, ordering);
    }
 
-   // Update surf_fit_gf (and optionally its gradients) if surface
-   // fitting is enabled.
+   // Update surf_fit_gf if surface fitting is enabled.
    if (surf_fit_gf)
    {
-      RemapSurfaceFittingLevelSetAtNodes(x_new, ordering);
+      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);
+      }
    }
 }
 

fem/tmop.hpp

@@ -1784,11 +1784,12 @@ 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;                  // 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.
+   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;
 
    DiscreteAdaptTC *discr_tc;
 
@@ -1984,10 +1985,6 @@ 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).
@@ -2003,8 +2000,9 @@ 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_grad(NULL), surf_fit_hess(NULL),
-        surf_fit_eval_grad(NULL), surf_fit_eval_hess(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),
         discr_tc(dynamic_cast<DiscreteAdaptTC *>(tc)),
         fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false)
    { PA.enabled = false; }
@@ -2105,15 +2103,9 @@ 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 *aegrad = NULL,
-                             AdaptivityEvaluator *aehess = NULL);
+                             AdaptivityEvaluator &ae);
 
    /** @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;
-         continue;
+         return -1;
       }
 #endif
 
@@ -148,7 +148,7 @@ int socketbuf::open(const char hostname[], int port)
    }
 
    freeaddrinfo(res);
-   return (socket_descriptor < 0) ? -1 : 0;
+   return 0;
 }
 
 int socketbuf::close()

linalg/densemat.cpp

@@ -532,69 +532,6 @@ 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,10 +207,6 @@ 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,24 +1151,15 @@ 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 either
-   // "material" or "attribute". We prioritize "material" over "attribute" for
-   // backwards compatibility.
+   // Read the element attributes, which are stored as CellData named "material"
    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())
    {
-      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)))
+      if (StringCompare(cell_data_xml->Name(), "CellData")
+          && StringCompare(cell_data_xml->Attribute("Scalars"), "material"))
       {
-         found_attributes = true;
          const XMLElement *data_xml = cell_data_xml->FirstChildElement();
          if (data_xml != NULL && StringCompare(data_xml->Name(), "DataArray"))
          {
@@ -1283,7 +1274,6 @@ 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);
@@ -1291,10 +1281,8 @@ 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 ||
-               (buff.rfind("SCALARS attribute") == 0 && !found_attributes))
+      else if (buff.rfind("SCALARS material") == 0)
       {
-         found_attributes = true;
          getline(input, buff); // LOOKUP_TABLE default
          if (buff.rfind("LOOKUP_TABLE default") != 0)
          {

miniapps/contact/contact.cpp

@@ -1,637 +0,0 @@
-//                               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

@@ -1,114 +0,0 @@
-#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

@@ -1,109 +0,0 @@
-# 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

@@ -1,453 +0,0 @@
-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

@@ -1,453 +0,0 @@
-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

@@ -1,231 +0,0 @@
-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

@@ -1,663 +0,0 @@
-#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

@@ -1,335 +0,0 @@
-
-#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

@@ -1,115 +0,0 @@
-#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

@@ -1,23 +0,0 @@
-
-#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

@@ -1,564 +0,0 @@
-//                               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,22 +580,6 @@ 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)
@@ -607,8 +591,11 @@ int main (int argc, char *argv[])
       {
 #ifdef MFEM_USE_GSLIB
          adapt_surface = new InterpolatorFP;
-         adapt_grad_surface = new InterpolatorFP;
-         adapt_hess_surface = new InterpolatorFP;
+         if (surf_bg_mesh)
+         {
+            adapt_grad_surface = new InterpolatorFP;
+            adapt_hess_surface = new InterpolatorFP;
+         }
 #else
          MFEM_ABORT("MFEM is not built with GSLIB support!");
 #endif
@@ -618,9 +605,7 @@ int main (int argc, char *argv[])
       if (!surf_bg_mesh)
       {
          tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker,
-                                          surf_fit_coeff, *adapt_surface,
-                                          adapt_grad_surface,
-                                          adapt_hess_surface);
+                                          surf_fit_coeff, *adapt_surface);
       }
       else
       {
@@ -852,13 +837,9 @@ 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 vis1, vis2, vis3;
-         common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0,
-                                "Level Set", 000, 400, 300, 300);
+         socketstream vis2, vis3;
          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

@@ -1,453 +0,0 @@
-// 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 contact-patch-finite-difference-test
+   MINIAPPS = contact-patch-test
 else
    MINIAPPS =
 endif

miniapps/tribol/modified-two-hex.mesh

@@ -1,57 +0,0 @@
-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 1.01
-1 0 1.01
-0 1 1.01 
-1 1 1.01
+0 0 0.99
+1 0 0.99
+0 1 0.99
+1 1 0.99
 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

@@ -1,58 +0,0 @@
-<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

@@ -1,64 +0,0 @@
-<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,66 +304,6 @@ 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,67 +623,6 @@ 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,82 +2811,4 @@ 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,35 +29,16 @@ TEST_CASE("VTU XML Reader", "[Mesh][VTU][XML]")
       "quad_binary_compress.vtu"
    });
 #endif
-
-   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)
+   for (const std::string &fname : mesh_filenames)
    {
-      REQUIRE(mesh_1.GetAttribute(i) == i+1);
-      REQUIRE(mesh_2.GetAttribute(i) == i+1);
+      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);
    }
 }
 
@@ -71,7 +52,8 @@ TEST_CASE("VTU XML Compressed Blocks", "[VTU][XML][MFEMData]")
                       "bracket_inline_compressed.vtu"
                    );
 
-   Mesh mesh = Mesh::LoadFromFile(mfem_data_dir + "/vtk/" + filename);
+   std::string mesh_path = mfem_data_dir + "/vtk/" + filename;
+   Mesh mesh = Mesh::LoadFromFile(mesh_path.c_str());
 
    REQUIRE(mesh.Dimension() == 3);
    REQUIRE(mesh.GetNE() == 206208);