make style

(cherry picked from commit 4c111d6f9198f418d32df98fd5e711455f73170a)

.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

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
 ====================================
@@ -38,6 +35,9 @@ Meshing improvements
 
 - Added support for internal boundary elements in nonconforming meshes.
 
+- Added ExodusII output capability. The writer can handle first-order (Pyramid5, 
+  Wedge6, Hex8, Tet4) and second-order FE types (Pyramid14, Wedge18, Hex27, Tet10).
+
 - The ReadCubit Genesis mesh importer has been rewritten to improve readability.
 
 Discretization improvements

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

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,11 +23,11 @@ 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
+   ex37p ex39p
 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

mesh/CMakeLists.txt

@@ -12,6 +12,7 @@
 set(SRCS
   attribute_sets.cpp
   element.cpp
+  exodus_writer.cpp
   face_nbr_geom.cpp
   gmsh.cpp
   hexahedron.cpp

mesh/mesh.hpp

@@ -2331,6 +2331,11 @@ public:
                     bool high_order_output=false,
                     int compression_level=0);
 
+#ifdef MFEM_USE_NETCDF
+   /// @brief Export a mesh to an Exodus II file.
+   void PrintExodusII(const std::string fpath);
+#endif
+
    /** @brief Prints the mesh with boundary elements given by the boundary of
        the subdomains, so that the boundary of subdomain i has boundary
        attribute i+1. */

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

@@ -386,6 +386,9 @@ int main (int argc, char *argv[])
            "S) Save in MFEM serial format\n"
            "T) Save in MFEM parallel format using the current partitioning\n"
            "V) Save in VTK format (only linear and quadratic meshes)\n"
+#ifdef MFEM_USE_NETCDF
+           "X) Save in Exodus II format (only linear and quadratic meshes)\n"
+#endif
            "D) Save as a DataCollection\n"
            "q) Quit\n"
 #ifdef MFEM_USE_ZLIB
@@ -1288,6 +1291,15 @@ int main (int argc, char *argv[])
          cout << "New VTK mesh file: " << omesh_file << endl;
       }
 
+#ifdef MFEM_USE_NETCDF
+      if (mk == 'X')
+      {
+         const char omesh_file[] = "mesh-explorer.e";
+         mesh->PrintExodusII(omesh_file);
+         cout << "New Exodus II mesh file: " << omesh_file << endl;
+      }
+#endif
+
       if (mk == 'D')
       {
          cout << "What type of DataCollection?\n"

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/CMakeLists.txt

@@ -56,6 +56,7 @@ set(UNIT_TESTS_SRCS
   mesh/test_submesh.cpp
   mesh/test_vtu.cpp
   mesh/test_nurbs.cpp
+  mesh/test_exodus_writer.cpp
   fem/test_1d_bilininteg.cpp
   fem/test_2d_bilininteg.cpp
   fem/test_3d_bilininteg.cpp

tests/unit/mesh/test_exodus_writer.cpp

@@ -0,0 +1,273 @@
+// 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.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+#ifdef MFEM_USE_NETCDF
+static void CompareMeshes(Mesh & mesh1, Mesh & mesh2)
+{
+   REQUIRE(mesh1.GetNE() == mesh2.GetNE());
+   REQUIRE(mesh1.GetNV() == mesh2.GetNV());
+   REQUIRE(mesh1.GetNBE() == mesh2.GetNBE());
+   REQUIRE(mesh1.GetNFaces() == mesh2.GetNFaces());
+
+   const FiniteElementSpace * fespace1 = mesh1.GetNodalFESpace();
+   const FiniteElementSpace * fespace2 = mesh2.GetNodalFESpace();
+
+   // Check elements.
+   Array<int> element_faces1, element_faces2;
+   Array<int> element_orient1, element_orient2;
+   Array<int> dofs1, dofs2;
+
+   for (int ielement = 0; ielement < mesh1.GetNE(); ielement++)
+   {
+      int attr1 = mesh1.GetAttribute(ielement);
+      int attr2 = mesh2.GetAttribute(ielement);
+
+      REQUIRE(attr1 == attr2);
+
+      Element::Type type1 = mesh1.GetElementType(ielement);
+      Element::Type type2 = mesh2.GetElementType(ielement);
+
+      REQUIRE(type1 == type2);
+
+      mesh1.GetElementFaces(ielement, element_faces1, element_orient1);
+      mesh2.GetElementFaces(ielement, element_faces2, element_orient2);
+
+      REQUIRE(element_faces1 == element_faces2);
+      REQUIRE(element_orient1 == element_orient2);
+
+      if (fespace1 && fespace2)
+      {
+         fespace1->GetElementDofs(ielement, dofs1);
+         fespace2->GetElementDofs(ielement, dofs2);
+      }
+      else
+      {
+         mesh1.GetElementVertices(ielement, dofs1);
+         mesh2.GetElementVertices(ielement, dofs2);
+      }
+
+      REQUIRE(dofs1 == dofs2);
+   }
+
+   // Check bdr elements.
+   for (int ibdr_element = 0; ibdr_element < mesh1.GetNBE(); ibdr_element++)
+   {
+      int attr1 = mesh1.GetBdrAttribute(ibdr_element);
+      int attr2 = mesh2.GetBdrAttribute(ibdr_element);
+
+      REQUIRE(attr1 == attr2);
+
+      Element::Type type1 = mesh1.GetBdrElementType(ibdr_element);
+      Element::Type type2 = mesh2.GetBdrElementType(ibdr_element);
+
+      REQUIRE(type1 == type2);
+
+      int face_index1 = mesh1.GetBdrElementFaceIndex(ibdr_element);
+      int face_index2 = mesh2.GetBdrElementFaceIndex(ibdr_element);
+
+      REQUIRE(face_index1 == face_index2);
+   }
+
+   // Check face vertices.
+   Array<int> face_vertices1, face_vertices2;
+   for (int iface_index = 0; iface_index < mesh1.GetNFaces(); iface_index++)
+   {
+      mesh1.GetFaceVertices(iface_index, face_vertices1);
+      mesh2.GetFaceVertices(iface_index, face_vertices2);
+
+      REQUIRE(face_vertices1 == face_vertices2);
+   }
+}
+#endif
+
+TEST_CASE("ExodusII Write Hex8", "[Mesh]")
+{
+#ifdef MFEM_USE_NETCDF
+   // Load Exodus II mesh from file.
+   std::string fpath_original = "data/simple-cube-hex8.e";
+   Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+   // Write generated Exodus II mesh to file.
+   std::string fpath_generated = "simple-cube-hex8-out.e";
+   original_mesh.PrintExodusII(fpath_generated);
+
+   // Load generated Exodus II mesh.
+   Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+
+   CompareMeshes(original_mesh, generated_mesh);
+
+   // Remove temporary file.
+   REQUIRE(remove(fpath_generated.c_str()) == 0);
+#endif
+}
+
+TEST_CASE("ExodusII Write Hex27", "[Mesh]")
+{
+#ifdef MFEM_USE_NETCDF
+   // Load Exodus II mesh from file.
+   std::string fpath_original = "data/simple-cube-hex27.e";
+   Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+   // Write generated Exodus II mesh to file.
+   std::string fpath_generated = "simple-cube-hex27-out.e";
+   original_mesh.PrintExodusII(fpath_generated);
+
+   // Load generated Exodus II mesh.
+   Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+
+   CompareMeshes(original_mesh, generated_mesh);
+
+   REQUIRE(remove(fpath_generated.c_str()) == 0);
+#endif
+}
+
+TEST_CASE("ExodusII Write Tet4", "[Mesh]")
+{
+#ifdef MFEM_USE_NETCDF
+   // Load Exodus II mesh from file. NB: - Do NOT refine as this changes vertex ordering!
+   std::string fpath_original = "data/simple-cube-tet4.e";
+   Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+   // Write generated Exodus II mesh to file.
+   std::string fpath_generated = "simple-cube-tet4-out.e";
+   original_mesh.PrintExodusII(fpath_generated);
+
+   // Load generated Exodus II mesh.
+   Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+
+   CompareMeshes(original_mesh, generated_mesh);
+
+   REQUIRE(remove(fpath_generated.c_str()) == 0);
+#endif
+}
+
+TEST_CASE("ExodusII Write Tet10", "[Mesh]")
+{
+#ifdef MFEM_USE_NETCDF
+   // Load Exodus II mesh from file.
+   std::string fpath_original = "data/simple-cube-tet10.e";
+   Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+   // Write generated Exodus II mesh to file.
+   std::string fpath_generated = "simple-cube-tet10-out.e";
+   original_mesh.PrintExodusII(fpath_generated);
+
+   // Load generated Exodus II mesh.
+   Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+
+   CompareMeshes(original_mesh, generated_mesh);
+
+   REQUIRE(remove(fpath_generated.c_str()) == 0);
+#endif
+}
+
+// TEST_CASE("ExodusII Write Wedge6", "[Mesh]")
+// {
+// #ifdef MFEM_USE_NETCDF
+//    std::string fpath_original = "data/simple-cube-wedge6.e";
+//    Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+//    std::string fpath_generated = "simple-cube-wedge6-out.e";
+//    original_mesh.PrintExodusII(fpath_generated);
+
+//    Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+//    CompareMeshes(original_mesh, generated_mesh);
+//
+//    REQUIRE(remove(fpath_generated.c_str()) == 0);
+// #endif
+// }
+
+// TEST_CASE("ExodusII Write Wedge18", "[Mesh]")
+// {
+// #ifdef MFEM_USE_NETCDF
+//    std::string fpath_original = "data/simple-cube-wedge18.e";
+//    Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+//    std::string fpath_generated = "simple-cube-wedge18-out.e";
+//    original_mesh.PrintExodusII(fpath_generated);
+
+//    Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+//    CompareMeshes(original_mesh, generated_mesh);
+
+//    REQUIRE(remove(fpath_generated.c_str()) == 0);
+// #endif
+// }
+
+// TEST_CASE("ExodusII Write Pyramid5", "[Mesh]")
+// {
+// #ifdef MFEM_USE_NETCDF
+//    std::string fpath_original = "data/simple-cube-pyramid5.e";
+//    Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+//    std::string fpath_generated = "simple-cube-pyramid5-out.e";
+//    original_mesh.PrintExodusII(fpath_generated);
+
+//    Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+//    CompareMeshes(original_mesh, generated_mesh);
+
+//    REQUIRE(remove(fpath_generated.c_str()) == 0);
+// #endif
+// }
+
+// TEST_CASE("ExodusII Write Pyramid14", "[Mesh]")
+// {
+// #ifdef MFEM_USE_NETCDF
+//    std::string fpath_original = "data/simple-cube-pyramid14.e";
+//    Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+//    std::string fpath_generated = "simple-cube-pyramid14-out.e";
+//    original_mesh.PrintExodusII(fpath_generated);
+
+//    Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+//    CompareMeshes(original_mesh, generated_mesh);
+
+//    REQUIRE(remove(fpath_generated.c_str()) == 0);
+// #endif
+// }
+
+// TEST_CASE("ExodusII Write Mixed First-Order", "[Mesh]")
+// {
+// #ifdef MFEM_USE_NETCDF
+//    // Contains Hex8, Tet4, Wedge6, Pyramid5 elements.
+//    std::string fpath_original = "data/simple-cube-multi-element-order1.e";
+//    Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+//    std::string fpath_generated = "simple-cube-multi-element-order1-out.e";
+//    original_mesh.PrintExodusII(fpath_generated);
+
+//    Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+//    CompareMeshes(original_mesh, generated_mesh);
+
+//    REQUIRE(remove(fpath_generated.c_str()) == 0);
+// #endif
+// }
+
+// TEST_CASE("ExodusII Write Mixed Second-Order", "[Mesh]")
+// {
+// #ifdef MFEM_USE_NETCDF
+//    // Contains Hex27 and Tet10 elements.
+//    std::string fpath_original = "data/simple-cube-multi-element-order2.e";
+//    Mesh original_mesh = Mesh::LoadFromFile(fpath_original, 0, 0, true);
+
+//    std::string fpath_generated = "simple-cube-multi-element-order2-out.e";
+//    original_mesh.PrintExodusII(fpath_generated);
+
+//    Mesh generated_mesh = Mesh::LoadFromFile(fpath_generated, 0, 0, true);
+//    CompareMeshes(original_mesh, generated_mesh);
+
+//    REQUIRE(remove(fpath_generated.c_str()) == 0);
+// #endif
+// }