small updates prior to passing the code along for comparison

miniapps/contact/contact.cpp

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

miniapps/contact/ipsolver/ParIPsolver.hpp

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

miniapps/contact/makefile

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

miniapps/contact/meshes/Test41.mesh

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

miniapps/contact/meshes/Test42.mesh

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

miniapps/contact/meshes/two-block.mesh

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

miniapps/contact/problems/parproblems.cpp

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

miniapps/contact/problems/parproblems.hpp

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

miniapps/contact/problems/parproblems_util.cpp

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

miniapps/contact/problems/parproblems_util.hpp

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

miniapps/contact/scratch.cpp

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

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

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

miniapps/tribol/makefile

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

miniapps/tribol/modified-two-hex.mesh

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

miniapps/tribol/two-hex.mesh

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