setting up test problem 6

config/config.hpp.in

@@ -198,6 +198,9 @@
 // Enable functionality based on the Google Benchmark library.
 // #define MFEM_USE_BENCHMARK
 
+// Enable interface to Tribol
+// #define MFEM_USE_TRIBOL
+
 // Enable the Enzyme LLVM plugin
 // #define MFEM_USE_ENZYME
 

linalg/operator.cpp

@@ -365,6 +365,44 @@ void SecondOrderTimeDependentOperator::ImplicitSolve(const double dt0,
    mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
 }
 
+SumOperator::SumOperator(const Operator *A, const double alpha,
+                         const Operator *B, const double beta,
+                         bool ownA, bool ownB)
+   : Operator(A->Height(), A->Width()),
+     A(A), B(B), alpha(alpha), beta(beta), ownA(ownA), ownB(ownB),
+     z(A->Height())
+{
+   MFEM_VERIFY(A->Width() == B->Width(),
+               "incompatible Operators: different widths\n"
+               << "A->Width() = " << A->Width()
+               << ", B->Width() = " << B->Width() );
+   MFEM_VERIFY(A->Height() == B->Height(),
+               "incompatible Operators: different heights\n"
+               << "A->Height() = " << A->Height()
+               << ", B->Height() = " << B->Height() );
+
+   {
+      const Solver* SolverA = dynamic_cast<const Solver*>(A);
+      const Solver* SolverB = dynamic_cast<const Solver*>(B);
+      if (SolverA)
+      {
+         MFEM_VERIFY(!(SolverA->iterative_mode),
+                     "Operator A of a SumOperator should not be in iterative mode");
+      }
+      if (SolverB)
+      {
+         MFEM_VERIFY(!(SolverB->iterative_mode),
+                     "Operator B of a SumOperator should not be in iterative mode");
+      }
+   }
+
+}
+
+SumOperator::~SumOperator()
+{
+   if (ownA) { delete A; }
+   if (ownB) { delete B; }
+}
 
 ProductOperator::ProductOperator(const Operator *A, const Operator *B,
                                  bool ownA, bool ownB)

linalg/operator.hpp

@@ -769,6 +769,28 @@ public:
    { A.Mult(x, y); }
 };
 
+/// General linear combination operator: x -> a A(x) + b B(x).
+class SumOperator : public Operator
+{
+   const Operator *A, *B;
+   const double alpha, beta;
+   bool ownA, ownB;
+   mutable Vector z;
+
+public:
+   SumOperator(
+      const Operator *A, const double alpha,
+      const Operator *B, const double beta,
+      bool ownA, bool ownB);
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   { z.SetSize(A->Height()); A->Mult(x, z); B->Mult(x, y); add(alpha, z, beta, y, y); }
+
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   { z.SetSize(A->Width()); A->MultTranspose(x, z); B->MultTranspose(x, y); add(alpha, z, beta, y, y); }
+
+   virtual ~SumOperator();
+};
 
 /// General product operator: x -> (A*B)(x) = A(B(x)).
 class ProductOperator : public Operator

makefile

@@ -742,7 +742,7 @@ ASTYLE_VER = "Artistic Style Version 3.1"
 FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
 FORMAT_FILES += tests/unit/*.?pp
 UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed
-MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners hooke/materials hooke/kernels
+MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners hooke/materials hooke/kernels tribol/contact
 FORMAT_FILES += $(foreach dir,$(UNIT_TESTS_SUBDIRS),tests/unit/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(MINIAPPS_SUBDIRS),miniapps/$(dir)/*.?pp)
 FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp

miniapps/tribol/contact/contact_driver.cpp

@@ -0,0 +1,151 @@
+//                                Contact example
+//
+// Compile with: make contact
+//
+// Sample runs:  ./contact -m1 block1.mesh -m2 block2.mesh -at "5 6 7 8"
+// Sample runs:  ./contact -m1 block1_d.mesh -m2 block2_d.mesh -at "5 6 7 8"
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "ipsolver/IPsolver.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file1 = "meshes/block1.mesh";
+   const char *mesh_file2 = "meshes/rotatedblock2.mesh";
+   int order = 1;
+   int ref = 0;
+   Array<int> attr;
+   Array<int> m_attr;
+   int linSolver = 2;
+   bool paraview = false;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file1, "-m1", "--mesh1",
+                  "First mesh file to use.");
+   args.AddOption(&mesh_file2, "-m2", "--mesh2",
+                  "Second mesh file to use.");
+   args.AddOption(&attr, "-at", "--attributes-surf",
+                  "Attributes of boundary faces on contact surface for mesh 2.");
+   args.AddOption(&ref, "-r", "--refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
+                  "--no-paraview",
+                  "Enable or disable ParaView visualization.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   ElasticityProblem prob1(mesh_file1,ref,order);
+   ElasticityProblem prob2(mesh_file2,ref,order);
+
+   ContactProblem contact(&prob1, &prob2);
+   QPOptContactProblem qpopt(&contact);
+   int numconstr = contact.GetNumConstraints();
+
+   InteriorPointSolver optimizer(&qpopt);
+   optimizer.SetTol(1e-6);
+   optimizer.SetMaxIter(50);
+   optimizer.SetLinearSolver(linSolver);
+   optimizer.SetLinearSolveTol(1e-10);
+
+   GridFunction x1 = prob1.GetDisplacementGridFunction();
+   GridFunction x2 = prob2.GetDisplacementGridFunction();
+
+   int ndofs1 = prob1.GetNumDofs();
+   int ndofs2 = prob2.GetNumDofs();
+   int ndofs = ndofs1 + ndofs2;
+
+   Vector x0(ndofs); x0 = 0.0;
+   x0.SetVector(x1,0);
+   x0.SetVector(x2,x1.Size());
+
+   Vector xf(ndofs); xf = 0.0;
+   optimizer.Mult(x0, xf);
+   Array<int> & CGiterations = optimizer.GetCGIterNumbers();
+
+   double Einitial = contact.E(x0);
+   double Efinal = contact.E(xf);
+
+   mfem::out << endl;
+   mfem::out << " Initial Energy objective     = " << Einitial << endl;
+   mfem::out << " Final Energy objective       = " << Efinal << endl;
+   mfem::out << " Global number of dofs        = " << ndofs1 + ndofs2 << endl;
+   mfem::out << " Global number of constraints = " << numconstr << endl;
+   mfem::out << " CG iteration numbers         = " ;
+   CGiterations.Print(mfem::out, CGiterations.Size());
+
+   MFEM_VERIFY(optimizer.GetConverged(),
+               "Interior point solver did not converge.");
+
+   if (visualization || paraview)
+   {
+      FiniteElementSpace * fes1 = prob1.GetFESpace();
+      FiniteElementSpace * fes2 = prob2.GetFESpace();
+
+      Mesh * mesh1 = fes1->GetMesh();
+      Mesh * mesh2 = fes2->GetMesh();
+
+      GridFunction x1_gf(fes1,xf.GetData());
+      GridFunction x2_gf(fes2,&xf.GetData()[fes1->GetTrueVSize()]);
+
+      mesh1->MoveNodes(x1_gf);
+      mesh2->MoveNodes(x2_gf);
+
+      if (paraview)
+      {
+         ParaViewDataCollection paraview_dc1("QPContactBody1", mesh1);
+         paraview_dc1.SetPrefixPath("ParaView");
+         paraview_dc1.SetLevelsOfDetail(1);
+         paraview_dc1.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc1.SetHighOrderOutput(true);
+         paraview_dc1.SetCycle(0);
+         paraview_dc1.SetTime(0.0);
+         paraview_dc1.RegisterField("Body1", &x1_gf);
+         paraview_dc1.Save();
+
+         ParaViewDataCollection paraview_dc2("QPContactBody2", mesh2);
+         paraview_dc2.SetPrefixPath("ParaView");
+         paraview_dc2.SetLevelsOfDetail(1);
+         paraview_dc2.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc2.SetHighOrderOutput(true);
+         paraview_dc2.SetCycle(0);
+         paraview_dc2.SetTime(0.0);
+         paraview_dc2.RegisterField("Body2", &x2_gf);
+         paraview_dc2.Save();
+      }
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int visport = 19916;
+         {
+            socketstream sol_sock(vishost, visport);
+            sol_sock.precision(8);
+            sol_sock << "parallel " << 2 << " " << 0 << "\n"
+                     << "solution\n" << *mesh1 << x1_gf << flush;
+         }
+         {
+            socketstream sol_sock(vishost, visport);
+            sol_sock.precision(8);
+            sol_sock << "parallel " << 2 << " " << 1 << "\n"
+                     << "solution\n" << *mesh2 << x2_gf << flush;
+         }
+      }
+   }
+
+   return 0;
+}

miniapps/tribol/contact/ipsolver/IPsolver.cpp

@@ -0,0 +1,818 @@
+#include "mfem.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+#include <cstdlib>
+
+using namespace std;
+using namespace mfem;
+
+
+
+InteriorPointSolver::InteriorPointSolver(QPOptContactProblem * Problem) 
+: optProblem(Problem), block_offsetsumlz(5), block_offsetsuml(4), block_offsetsx(3),
+saveLogBarrierIterates(false)
+{
+   rel_tol  = 1.e-2;
+   max_iter = 20;
+   mu_k     = 1.0;
+
+   sMax     = 1.e2;
+   kSig     = 1.e10;   // control deviation from primal Hessian
+   tauMin   = 0.8;     // control rate at which iterates can approach the boundary
+   eta      = 1.e-4;   // backtracking constant
+   thetaMin = 1.e-4;   // allowed violation of the equality constraints
+
+   // constants in line-step A-5.4
+   delta    = 1.0;
+   sTheta   = 1.1;
+   sPhi     = 2.3;
+
+   // control the rate at which the penalty parameter is decreased
+   kMu     = 0.2;
+   thetaMu = 1.5;
+
+   // TO DO -- include the filter
+
+   thetaMax = 1.e6; // maximum constraint violation
+   // data for the second order correction
+   kSoc     = 0.99;
+
+   // equation (18)
+   gTheta = 1.e-5;
+   gPhi   = 1.e-5;
+
+   kEps   = 1.e1;
+
+   dimU = optProblem->GetDimU();
+   dimM = optProblem->GetDimM();
+   dimC = optProblem->GetDimC();
+   ckSoc.SetSize(dimC);
+   
+   block_offsetsumlz[0] = 0;
+   block_offsetsumlz[1] = dimU; // u
+   block_offsetsumlz[2] = dimM; // m
+   block_offsetsumlz[3] = dimC; // lambda
+   block_offsetsumlz[4] = dimM; // zl
+   block_offsetsumlz.PartialSum();
+
+   for(int i = 0; i < block_offsetsuml.Size(); i++)  { block_offsetsuml[i] = block_offsetsumlz[i]; }
+   for(int i = 0; i < block_offsetsx.Size(); i++)    { block_offsetsx[i]   = block_offsetsuml[i] ; }
+
+   // lower-bound for the inequality constraint m >= ml
+   ml = optProblem->Getml();
+   
+   lk.SetSize(dimC);  lk  = 0.0;
+   zlk.SetSize(dimM); zlk = 0.0;
+
+   linSolver = 0;
+   MyRank = 0;
+   iAmRoot = MyRank == 0 ? true : false;
+}
+
+double InteriorPointSolver::MaxStepSize(Vector &x, Vector &xl, Vector &xhat, double tau)
+{
+   double alphaMaxloc = 1.0;
+   double alphaTmp;
+   for(int i = 0; i < x.Size(); i++)
+   {
+      if( xhat(i) < 0. )
+      {
+         alphaTmp = -1. * tau * (x(i) - xl(i)) / xhat(i);
+         alphaMaxloc = min(alphaMaxloc, alphaTmp);
+      } 
+   }
+
+   // alphaMaxloc is the local maximum step size which is
+   // distinct on each MPI process. Need to compute
+   // the global maximum step size 
+   double alphaMaxglb;
+   alphaMaxglb = alphaMaxloc;
+   return alphaMaxglb;
+}
+
+double InteriorPointSolver::MaxStepSize(Vector &x, Vector &xhat, double tau)
+{
+   Vector zero(x.Size()); zero = 0.0;
+   return MaxStepSize(x, zero, xhat, tau);
+}
+
+
+void InteriorPointSolver::Mult(const Vector &x0, Vector &xf)
+{
+   BlockVector x0block(block_offsetsx); x0block = 0.0;
+   x0block.GetBlock(0).Set(1.0, x0);
+   // To do: give options for user specificiation of initialization m0
+   x0block.GetBlock(1) = 1.0;
+   x0block.GetBlock(1).Add(1.0, ml);
+   BlockVector xfblock(block_offsetsx); xfblock = 0.0;
+   Mult(x0block, xfblock);
+   xf.Set(1.0, xfblock.GetBlock(0));
+}
+
+void InteriorPointSolver::Mult(const BlockVector &x0, BlockVector &xf)
+{
+   converged = false;
+   
+   BlockVector xk(block_offsetsx), xhat(block_offsetsx); xk = 0; xhat = 0.0;
+   BlockVector Xk(block_offsetsumlz), Xhat(block_offsetsumlz); Xk = 0.0; Xhat = 0.0;
+   BlockVector Xhatuml(block_offsetsuml); Xhatuml = 0.0;
+   Vector zlhat(dimM); zlhat = 0.0;
+
+   xk.GetBlock(0).Set(1.0, x0.GetBlock(0));
+   xk.GetBlock(1).Set(1.0, x0.GetBlock(1));
+   // running estimate of the final values of the Lagrange multipliers
+   lk  = 0.0;
+   zlk = 0.0;
+
+   for(int i = 0; i < dimM; i++)
+   {
+      zlk(i) = 1.e1 * mu_k / (xk(i+dimU) - ml(i));
+   }
+
+   Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
+   Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
+   Xk.GetBlock(2).Set(1.0, lk);
+   Xk.GetBlock(3).Set(1.0, zlk);
+
+   /* set theta0 = theta(x0)
+      *     thetaMin
+      *     thetaMax
+      * when theta(xk) < thetaMin and the switching condition holds
+      * then we ask for the Armijo sufficient decrease of the barrier
+      * objective to be satisfied, in order to accept the trial step length alphakl
+      * 
+      * thetaMax controls how the filter is initialized for each log-barrier subproblem
+      * F0 = {(th, phi) s.t. th > thetaMax}
+      * that is the filter does not allow for iterates where the constraint violation
+      * is larger than that of thetaMax
+      */
+   double theta0 = theta(xk);
+   thetaMin = 1.e-4 * max(1.0, theta0);
+   thetaMax = 1.e8  * thetaMin;
+
+   double Eeval, maxBarrierSolves, Eevalmu0;
+   bool printOptimalityError; // control optimality error print to console for log-barrier subproblems
+   
+   maxBarrierSolves = 10;
+
+   for(jOpt = 0; jOpt < max_iter; jOpt++)
+   {
+      mfem::out << "interior-point solve step " << jOpt << endl;
+      // A-2. Check convergence of overall optimization problem
+      printOptimalityError = false;
+      Eevalmu0 = E(xk, lk, zlk, printOptimalityError);
+      if(Eevalmu0 < rel_tol)
+      {
+         converged = true;
+         mfem::out << "solved optimization problem :)\n";
+         break;
+      }
+      
+      if(jOpt > 0) { maxBarrierSolves = 1; }
+      
+      for(int i = 0; i < maxBarrierSolves; i++)
+      {
+         // A-3. Check convergence of the barrier subproblem
+         printOptimalityError = true;
+         Eeval = E(xk, lk, zlk, mu_k, printOptimalityError);
+         if(Eeval < kEps * mu_k)
+         {
+            mfem::out << "solved barrier subproblem, for mu = " << mu_k << endl;
+            // A-3.1. Recompute the barrier parameter
+            mu_k  = max(rel_tol / 10., min(kMu * mu_k, pow(mu_k, thetaMu)));
+            // A-3.2. Re-initialize the filter
+            F1.DeleteAll();
+            F2.DeleteAll();
+         }
+         else
+         {
+            break;
+         }
+      }
+      // A-4. Compute the search direction
+      // solve for (uhat, mhat, lhat)
+      mfem::out << "\n** A-4. IP-Newton solve **\n";
+      zlhat = 0.0; Xhatuml = 0.0;
+      // why do we have Xhatuml ....???
+      // TO DO: remove Xhatuml in favor of passing Xhat
+      IPNewtonSolve(xk, lk, zlk, zlhat, Xhatuml, mu_k, false); 
+
+      // assign data stack, X = (u, m, l, zl)
+      Xk = 0.0;
+      Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
+      Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
+      Xk.GetBlock(2).Set(1.0, lk);
+      Xk.GetBlock(3).Set(1.0, zlk);
+
+      // assign data stack, Xhat = (uhat, mhat, lhat, zlhat)
+      Xhat = 0.0;
+      for(int i = 0; i < 3; i++)
+      {
+         Xhat.GetBlock(i).Set(1.0, Xhatuml.GetBlock(i));
+      }
+      Xhat.GetBlock(3).Set(1.0, zlhat);
+
+      // A-5. Backtracking line search.
+      mfem::out << "\n** A-5. Linesearch **\n";
+      mfem::out << "mu = " << mu_k << endl;
+      
+      lineSearch(Xk, Xhat, mu_k);
+      if(lineSearchSuccess)
+      {
+         if(!switchCondition || !sufficientDecrease)
+         {
+            F1.Append( (1. - gTheta) * thx0);
+            F2.Append( phx0 - gPhi * thx0);
+         }
+         // ----- A-6: Accept the trial point
+         // print info regarding zl...
+         xk.GetBlock(0).Add(alpha, Xhat.GetBlock(0));
+         xk.GetBlock(1).Add(alpha, Xhat.GetBlock(1));
+         lk.Add(alpha,   Xhat.GetBlock(2));
+         zlk.Add(alphaz, Xhat.GetBlock(3));
+         projectZ(xk, zlk, mu_k);
+      }
+      else
+      {
+         mfem::out << "lineSearch not successful :(\n";
+         mfem::out << "attempting feasibility restoration with theta = " << thx0 << endl;
+         mfem::out << "no feasibility restoration implemented, exiting now \n";
+         break;
+      }
+      //
+      if(jOpt + 1 == max_iter) 
+      { 
+         mfem::out << "maximum optimization iterations :(\n";
+      }
+   }
+   // done with optimization routine, just reassign data to xf reference so
+   // that the application code has access to the optimal point
+   xf = 0.0;
+   xf.GetBlock(0).Set(1.0, xk.GetBlock(0));
+   xf.GetBlock(1).Set(1.0, xk.GetBlock(1));
+}
+
+void InteriorPointSolver::FormIPNewtonMat(BlockVector & x, Vector & l, Vector &zl, BlockOperator &Ak)
+{
+   // WARNING: Huu, Hum, Hmu, Hmm should all be Hessian terms of the Lagrangian, currently we 
+   //          them by Hessian terms of the objective function and neglect the Hessian of l^T c
+
+   Huu = optProblem->Duuf(x);
+   Hum = optProblem->Dumf(x);
+   Hmu = optProblem->Dmuf(x);
+   Hmm = optProblem->Dmmf(x);
+
+   Vector DiagLogBar(dimM); DiagLogBar = 0.0;
+   for(int ii = 0; ii < dimM; ii++)
+   {
+      DiagLogBar(ii) = zl(ii) / (x(ii+dimU) - ml(ii));
+   }
+
+   if(saveLogBarrierIterates)
+   {
+      std::ofstream diagStream;
+      char diagString[100];
+      snprintf(diagString, 100, "logBarrierHessiandata/D%d.dat", jOpt);
+      diagStream.open(diagString, ios::out | ios::trunc);
+      for(int ii = 0; ii < dimM; ii++)
+      {
+         diagStream << setprecision(30) << DiagLogBar(ii) << endl;
+      }
+      diagStream.close();
+   } 
+
+   delete Wmm;
+   if(Hmm != nullptr)
+   {
+      SparseMatrix * D = new SparseMatrix(DiagLogBar);
+      Wmm = Add(*Hmm, *D);
+      delete D;
+   }
+   else
+   {
+      Wmm = new SparseMatrix(DiagLogBar);
+   }
+
+   delete JuT;
+   delete JmT;
+   Ju = optProblem->Duc(x); JuT = Transpose(*Ju);
+   Jm = optProblem->Dmc(x); JmT = Transpose(*Jm);
+   
+   Huucl = optProblem->lDuuc(x, l);
+   if(Huucl != nullptr)
+   {
+      delete HLuucl;
+      HLuucl = Add(*Huucl, *Huu);
+      Ak.SetBlock(0, 0, HLuucl);
+   } 
+   else
+   {
+      Ak.SetBlock(0, 0, Huu);
+   }
+
+   //         IP-Newton system matrix
+   //    Ak = [[H_(u,u)  H_(u,m)   J_u^T]
+   //          [H_(m,u)  W_(m,m)   J_m^T]
+   //          [ J_u      J_m       0  ]]
+
+   Ak.SetBlock(0, 0, Huu);
+   Ak.SetBlock(0, 2, JuT);
+   Ak.SetBlock(1, 1, Wmm);   
+   Ak.SetBlock(1, 2, JmT);
+   Ak.SetBlock(2, 0, Ju); 
+   Ak.SetBlock(2, 1, Jm); 
+
+   if(Hum != nullptr) { Ak.SetBlock(0, 1, Hum); Ak.SetBlock(1, 0, Hmu); }
+}
+
+// perturbed KKT system solve
+// determine the search direction
+void InteriorPointSolver::IPNewtonSolve(BlockVector &x, Vector &l, Vector &zl, Vector &zlhat, BlockVector &Xhat, double mu, bool socSolve)
+{
+   // solve A x = b, where A is the IP-Newton matrix
+   BlockOperator A(block_offsetsuml, block_offsetsuml); BlockVector b(block_offsetsuml); b = 0.0;
+   FormIPNewtonMat(x, l, zl, A);
+
+   //       [grad_u phi + Ju^T l]
+   // b = - [grad_m phi + Jm^T l]
+   //       [          c        ]
+   BlockVector gradphi(block_offsetsx); gradphi = 0.0;
+   BlockVector JTl(block_offsetsx); JTl = 0.0;
+   Dxphi(x, mu, gradphi);
+  
+   (A.GetBlock(0,2)).Mult(l, JTl.GetBlock(0));
+   (A.GetBlock(1,2)).Mult(l, JTl.GetBlock(1));
+
+   
+   for(int ii = 0; ii < 2; ii++)
+   {
+      b.GetBlock(ii).Set(1.0, gradphi.GetBlock(ii));
+      b.GetBlock(ii).Add(1.0, JTl.GetBlock(ii));
+   }
+   if(!socSolve) 
+   {
+      optProblem->c(x, b.GetBlock(2));
+   }
+   else
+   {
+      b.GetBlock(2).Set(1.0, ckSoc);
+   }
+   b *= -1.0; 
+   Xhat = 0.0;
+
+
+#ifdef MFEM_USE_SUITESPARSE
+   // Direct solve for IP-Newton saddle-point system
+   // A = [ [ Huu   0    Ju^T]
+   //       [  0    D     -I ]
+   //       [ Ju   -I      0 ]]
+   // if(linSolver == 0)
+   // {
+   //    BlockMatrix ABlockMatrix(block_offsetsuml, block_offsetsuml);
+   //    for(int ii = 0; ii < 3; ii++)
+   //    {
+   //       for(int jj = 0; jj < 3; jj++)
+   //       {
+   //          if(!A.IsZeroBlock(ii, jj))
+   //          {
+   //             ABlockMatrix.SetBlock(ii, jj, dynamic_cast<SparseMatrix *>(&(A.GetBlock(ii, jj))));
+   //          }
+   //       }
+   //    }
+   //    /* direct solve of the 3x3 IP-Newton linear system */
+   //    UMFPackSolver ASolver;
+   //    SparseMatrix *ASparse = ABlockMatrix.CreateMonolithic();
+   //    ASolver.SetOperator(*ASparse);
+   //    ASolver.Mult(b, Xhat);
+
+   //    Vector residual(Xhat.Size());
+   //    ASparse->Mult(Xhat, residual);
+   //    residual.Add(-1.0, b);
+   //    delete ASparse;
+   // }
+   // else if(linSolver == 1)
+   // {
+   //    // Direct solve for 0,0 Schur complement of IP-Newton system, Huu + Ju^T Wmm Ju,
+   //    // where Wmm = D for contact problems
+   //    SparseMatrix * Huuloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 0)));
+   //    SparseMatrix * Wmmloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(1, 1)));
+   //    SparseMatrix * Juloc  = dynamic_cast<SparseMatrix *>(&(A.GetBlock(2, 0)));
+   //    SparseMatrix * JuTloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 2)));
+   //    Vector DVec(dimM); DVec = 0.0;
+   //    Vector one(dimM); one = 1.0;
+   //    D->Mult(one, DVec);
+   //    SparseMatrix *JuTDJu = Mult_AtDA(*Juloc, DVec);     // Ju^T D Ju
+   //    SparseMatrix *Areduced = Add(*Huuloc, *JuTDJu);  // Huu + Ju^T D Ju
+      
+   //    /* prepare the reduced rhs */
+   //    // breduced = bu + Ju^T (bm + Wmm bl)
+   //    Vector breduced(dimU); breduced = 0.0;
+   //    Vector tempVec(dimM); tempVec = 0.0;
+   //    Wmmloc->Mult(b.GetBlock(2), tempVec);
+   //    tempVec.Add(1.0, b.GetBlock(1));
+   //    JuTloc->Mult(tempVec, breduced);
+   //    breduced.Add(1.0, b.GetBlock(0));
+
+   //    // solve the reduced linear system
+   //    UMFPackSolver AreducedSolver;
+   //    AreducedSolver.SetOperator(*Areduced);
+   //    AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
+
+   //    // now propagate solved uhat to obtain mhat and lhat
+   //    // xm = Ju xu - bl
+   //    Juloc->Mult(Xhat.GetBlock(0), Xhat.GetBlock(1));
+   //    Xhat.GetBlock(1).Add(-1.0, b.GetBlock(2));
+
+   //    // xl = Wmm xm - bm
+   //    Wmmloc->Mult(Xhat.GetBlock(1), Xhat.GetBlock(2));
+   //    Xhat.GetBlock(2).Add(-1.0, b.GetBlock(1));
+
+   //    delete JuTDJu;
+   //    delete Areduced;
+   // }
+#else
+   MFEM_VERIFY(linSolver > 1, "linSolver = 0, 1 require MFEM_USE_SUITESPARSE=YES");
+#endif
+   // if(linSolver ==2)
+   { 
+      // Iterative solve for 0,0 Schur complement of IP-Newton system, Huu + Ju^T Wmm Ju,
+      // where Wmm = D for contact problems
+      // here the iterative solver is a Jacobi-preconditioned CG-solve
+      SparseMatrix * Huuloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 0)));
+      SparseMatrix * Wmmloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(1, 1)));
+      SparseMatrix * Juloc  = dynamic_cast<SparseMatrix *>(&(A.GetBlock(2, 0)));
+      SparseMatrix * JuTloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 2)));
+
+      SparseMatrix *JuTDJu   = RAP(*Juloc,*Wmmloc,*Juloc);     // Ju^T D Ju
+      SparseMatrix *Areduced = Add(*Huuloc, *JuTDJu);  // Huu + Ju^T D Ju
+
+      /* prepare the reduced rhs */
+      // breduced = bu + Ju^T (bm + Wmm bl)
+      Vector breduced(dimU); breduced = 0.0;
+      Vector tempVec(dimM); tempVec = 0.0;
+      Wmmloc->SortColumnIndices();
+
+      Wmmloc->Mult(b.GetBlock(2), tempVec);
+      tempVec.Add(1.0, b.GetBlock(1));
+      JuTloc->Mult(tempVec, breduced);
+
+
+      breduced.Add(1.0, b.GetBlock(0));
+      int globalNumRows = dimU;
+      HYPRE_BigInt rowStarts[2];
+      rowStarts[0] = 0;
+      rowStarts[1] = dimU;
+
+      HypreParMatrix Ahypre(MPI_COMM_WORLD, globalNumRows, rowStarts, Areduced);
+      HypreBoomerAMG Aprec(Ahypre);
+      Aprec.SetPrintLevel(0);
+      Aprec.SetSystemsOptions(3,false);
+      HyprePCG AreducedSolver(MPI_COMM_WORLD);
+      AreducedSolver.SetOperator(Ahypre);
+      // AreducedSolver.SetRelTol(linSolveTol);
+      // AreducedSolver.SetRelTol(1e-6);
+      AreducedSolver.SetTol(1e-6);
+      AreducedSolver.SetMaxIter(1000);
+      AreducedSolver.SetPreconditioner(Aprec);
+      // AreducedSolver.SetResidualConvergenceOptions();
+      AreducedSolver.SetPrintLevel(2);
+
+      AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
+      int num_iterations;
+      AreducedSolver.GetNumIterations(num_iterations);
+      cgnum_iterations.Append(num_iterations);
+    
+      // now propagate solved uhat to obtain mhat and lhat
+      // xm = Ju xu - bl
+      Juloc->Mult(Xhat.GetBlock(0), Xhat.GetBlock(1));
+      Xhat.GetBlock(1).Add(-1.0, b.GetBlock(2));
+      // xl = Wmm xm - bm
+      Wmmloc->Mult(Xhat.GetBlock(1), Xhat.GetBlock(2));
+      Xhat.GetBlock(2).Add(-1.0, b.GetBlock(1));
+
+      delete JuTDJu;
+      delete Areduced;
+   }
+    
+   /* backsolve to determine zlhat */
+   for(int ii = 0; ii < dimM; ii++)
+   {
+      zlhat(ii) = -1.*(zl(ii) + (zl(ii) * Xhat(ii + dimU) - mu) / (x(ii + dimU) - ml(ii)) );
+   }
+}
+
+// here Xhat, X will be BlockVectors w.r.t. the 4 partitioning X = (u, m, l, zl)
+
+void InteriorPointSolver::lineSearch(BlockVector& X0, BlockVector& Xhat, double mu)
+{
+  double tau  = max(tauMin, 1.0 - mu);
+  Vector u0   = X0.GetBlock(0);
+  Vector m0   = X0.GetBlock(1);
+  Vector l0   = X0.GetBlock(2);
+  Vector z0   = X0.GetBlock(3);
+  Vector uhat = Xhat.GetBlock(0);
+  Vector mhat = Xhat.GetBlock(1);
+  Vector lhat = Xhat.GetBlock(2);
+  Vector zhat = Xhat.GetBlock(3);
+  double alphaMax  = MaxStepSize(m0, ml, mhat, tau);
+  double alphaMaxz = MaxStepSize(z0, zhat, tau);
+  alphaz = alphaMaxz;
+
+
+  BlockVector x0(block_offsetsx); x0 = 0.0;
+  x0.GetBlock(0).Set(1.0, u0);
+  x0.GetBlock(1).Set(1.0, m0);
+  
+  BlockVector xhat(block_offsetsx); xhat = 0.0;
+  xhat.GetBlock(0).Set(1.0, uhat);
+  xhat.GetBlock(1).Set(1.0, mhat);
+  
+  BlockVector xtrial(block_offsetsx); xtrial = 0.0;
+  BlockVector Dxphi0(block_offsetsx); Dxphi0 = 0.0;
+  int maxBacktrack = 20;
+  alpha = alphaMax;
+
+
+  Vector ck0(dimC); ck0 = 0.0;
+  Vector zhatsoc(dimM); zhatsoc = 0.0;
+  BlockVector Xhatumlsoc(block_offsetsuml); Xhatumlsoc = 0.0;
+  BlockVector xhatsoc(block_offsetsx); xhatsoc = 0.0;
+  Vector uhatsoc(dimU); uhatsoc = 0.0;
+  Vector mhatsoc(dimM); mhatsoc = 0.0;
+
+  Dxphi(x0, mu, Dxphi0);
+  Dxphi0_xhat = InnerProduct(Dxphi0, xhat);
+  descentDirection = Dxphi0_xhat < 0. ? true : false;
+  if(descentDirection)
+  {
+    mfem::out << "is a descent direction for the log-barrier objective\n";
+  }
+  else
+  {
+    mfem::out << "is not a descent direction for the log-barrier objective\n";
+  }
+  mfem::out << "Dxphi^T xhat / (|| Dxphi||_2 * || xhat ||_2) = " << Dxphi0_xhat / (xhat.Norml2() * Dxphi0.Norml2()) << endl;
+  thx0 = theta(x0);
+  phx0 = phi(x0, mu);
+
+  lineSearchSuccess = false;
+  for(int i = 0; i < maxBacktrack; i++)
+  {
+    mfem::out << "\n--------- alpha = " << alpha << " ---------\n";
+
+    // ----- A-5.2. Compute trial point: xtrial = x0 + alpha_i xhat
+    xtrial.Set(1.0, x0);
+    xtrial.Add(alpha, xhat);
+
+    // ------ A-5.3. if not in filter region go to A.5.4 otherwise go to A-5.5.
+    thxtrial = theta(xtrial);
+    phxtrial = phi(xtrial, mu);
+
+    filterCheck(thxtrial, phxtrial);    
+    if(!inFilterRegion)
+    {
+      mfem::out << "not in filter region :)\n";
+      // ------ A.5.4: Check sufficient decrease
+      if(!descentDirection)
+      {
+        switchCondition = false;
+      }
+      else
+      {
+        switchCondition = (alpha * pow(abs(Dxphi0_xhat), sPhi) > delta * pow(thx0, sTheta)) ? true : false;
+      }
+      mfem::out << "theta(x0) = "     << thx0     << ", thetaMin = "                  << thetaMin             << endl;
+      mfem::out << "theta(xtrial) = " << thxtrial << ", (1-gTheta) *theta(x0) = "     << (1. - gTheta) * thx0 << endl;
+      mfem::out << "phi(xtrial) = "   << phxtrial << ", phi(x0) - gPhi *theta(x0) = " << phx0 - gPhi * thx0   << endl;
+      
+      // Case I      
+      if(thx0 <= thetaMin && switchCondition)
+      {
+        sufficientDecrease = phxtrial <= phx0 + eta * alpha * Dxphi0_xhat ? true : false;
+        if(sufficientDecrease)
+        {
+          mfem::out << "Accepted step length -- sufficient decrease in log-barrier objective.\n"; 
+          // accept the trial step
+          lineSearchSuccess = true;
+          break;
+        }
+      }
+      else
+      {
+        if(thxtrial <= (1. - gTheta) * thx0 || phxtrial <= phx0 - gPhi * thx0)
+        {
+          mfem::out << "Accepted step length -- decrease in either constraint violation or log-barrier objective.\n"; 
+          // accept the trial step
+          lineSearchSuccess = true;
+          break;
+        }
+      }
+      // A-5.5: Initialize the second-order correction
+      if((!(thx0 < thxtrial)) && i == 0)
+      {
+        mfem::out << "second order correction\n";
+        optProblem->c(xtrial, ckSoc);
+        optProblem->c(x0, ck0);
+        ckSoc.Add(alphaMax, ck0);
+        // A-5.6 Compute the second-order correction.
+        IPNewtonSolve(x0, l0, z0, zhatsoc, Xhatumlsoc, mu, true);
+        mhatsoc.Set(1.0, Xhatumlsoc.GetBlock(1));
+        // alphasoc = MaxStepSize(m0, ml, mhatsoc, tau);
+        //WARNING: not complete but currently solver isn't entering this region
+      }
+    }
+    else
+    {
+      mfem::out << "in filter region\n";
+    }
+
+    // include more if needed
+    alpha *= 0.5;
+
+  } 
+}
+
+
+void InteriorPointSolver::projectZ(const Vector &x, Vector &z, double mu)
+{
+  double zi;
+  double mudivmml;
+  for(int i = 0; i < dimM; i++)
+  {
+    zi = z(i);
+    mudivmml = mu / (x(i + dimU) - ml(i));
+    z(i) = max(min(zi, kSig * mudivmml), mudivmml / kSig);
+  }
+}
+
+void InteriorPointSolver::filterCheck(double th, double ph)
+{
+  inFilterRegion = false;
+  if(th > thetaMax)
+  {
+    inFilterRegion = true;
+  }
+  else
+  {
+    for(int i = 0; i < F1.Size(); i++)
+    {
+      if(th >= F1[i] && ph >= F2[i])
+      {
+        inFilterRegion = true;
+        break;
+      }
+    }
+  }
+}
+
+double InteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, double mu, bool print)
+{
+  double E1, E2, E3;
+  double sc, sd;
+  BlockVector gradL(block_offsetsx); gradL = 0.0; // stationarity grad L = grad f + J^T l - z
+  Vector cx(dimC); cx = 0.0;     // feasibility c = c(x)
+  Vector comp(dimM); comp = 0.0; // complementarity M Z - mu 1
+
+  DxL(x, l, zl, gradL);
+  E1 = gradL.Normlinf();
+
+  optProblem->c(x, cx);
+  E2 = cx.Normlinf();
+
+  for(int ii = 0; ii < dimM; ii++) 
+  { 
+    comp(ii) = x(dimU + ii) * zl(ii) - mu;
+  }
+  E3 = comp.Normlinf();
+
+  double ll1, zl1;
+  zl1 = zl.Norml1() / double(dimC + dimM);
+  ll1 =  l.Norml1();
+  sc = max(sMax, zl1 / (double(dimM)) ) / sMax;
+  sd = max(sMax, (ll1 + zl1) / (double(dimC + dimM))) / sMax;
+  if(print)
+  {
+    mfem::out << "evaluating optimality error for mu = " << mu << endl;
+    mfem::out << "stationarity measure = "    << E1 / sd << endl;
+    mfem::out << "feasibility measure  = "    << E2      << endl;
+    mfem::out << "complimentarity measure = " << E3 / sc << endl;
+  }
+  return max(max(E1 / sd, E2), E3 / sc);
+}
+
+double InteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, bool print)
+{
+  return E(x, l, zl, 0.0, print);
+}
+
+double InteriorPointSolver::theta(const BlockVector &x)
+{
+  Vector cx(dimC); cx = 0.0;
+  optProblem->c(x, cx);
+  return cx.Norml2();
+}
+
+// log-barrier objective
+double InteriorPointSolver::phi(const BlockVector &x, double mu)
+{
+  double fx = optProblem->CalcObjective(x); 
+  double logBarrierLoc = 0.0;
+  for(int i = 0; i < dimM; i++) 
+  { 
+    logBarrierLoc += log(x(dimU+i)-ml(i));
+  }
+  double logBarrierGlb = 0.0;
+  logBarrierGlb = logBarrierLoc;
+  return fx - mu * logBarrierGlb;
+}
+
+// gradient of log-barrier objective with respect to x = (u, m)
+void InteriorPointSolver::Dxphi(const BlockVector &x, double mu, BlockVector &y)
+{
+  optProblem->CalcObjectiveGrad(x, y);
+  for(int i = 0; i < dimM; i++) 
+  { 
+    y(dimU + i) -= mu / (x(dimU + i));
+  } 
+}
+
+// Lagrangian function evaluation
+// L(x, l, zl) = f(x) + l^T c(x) - zl^T m
+double InteriorPointSolver::L(const BlockVector &x, const Vector &l, const Vector &zl)
+{
+  double fx = optProblem->CalcObjective(x);
+  Vector cx(dimC); optProblem->c(x, cx);
+  return (fx + InnerProduct(cx, l) - InnerProduct(x.GetBlock(1), zl));
+}
+
+void InteriorPointSolver::DxL(const BlockVector &x, const Vector &l, const Vector &zl, BlockVector &y)
+{
+  // evaluate the gradient of the objective with respect to the primal variables x = (u, m)
+  BlockVector gradxf(block_offsetsx); gradxf = 0.0;
+  optProblem->CalcObjectiveGrad(x, gradxf);
+  
+  SparseMatrix *Jacu, *Jacm, *JacuT, *JacmT;
+  Jacu = optProblem->Duc(x); Jacm = optProblem->Dmc(x);
+  JacuT = Transpose(*Jacu);
+  JacmT = Transpose(*Jacm);
+  JacuT->Mult(l, y.GetBlock(0));
+  JacmT->Mult(l, y.GetBlock(1));
+  delete JacuT; 
+  delete JacmT;
+  y.Add(1.0, gradxf);
+  (y.GetBlock(1)).Add(-1.0, zl);
+}
+
+
+bool InteriorPointSolver::GetConverged() const
+{
+  return converged;
+}
+
+void InteriorPointSolver::SetTol(double Tol)
+{
+  rel_tol = Tol;
+}
+
+void InteriorPointSolver::SetMaxIter(int max_it)
+{
+  max_iter = max_it;
+}
+
+void InteriorPointSolver::SetBarrierParameter(double mu_0)
+{
+  mu_k = mu_0;
+}
+
+void InteriorPointSolver::SaveLogBarrierHessianIterates(bool save)
+{
+   MFEM_ASSERT(MyRank == 0 || save == false, "currently can only save logbarrier hessian in serial codes");
+   saveLogBarrierIterates = save;
+}
+
+void InteriorPointSolver::SetLinearSolver(int LinSolver)
+{
+   linSolver = LinSolver;
+}
+
+void InteriorPointSolver::SetLinearSolveTol(double Tol)
+{
+   linSolveTol = Tol;
+}
+
+
+InteriorPointSolver::~InteriorPointSolver() 
+{
+   delete HLuucl;
+   delete JuT;
+   delete JmT; 
+   delete Wmm;
+
+   F1.DeleteAll();
+   F2.DeleteAll();
+   block_offsetsx.DeleteAll();
+   block_offsetsumlz.DeleteAll();
+   block_offsetsuml.DeleteAll();
+   ml.SetSize(0);
+}

miniapps/tribol/contact/ipsolver/IPsolver.hpp

@@ -0,0 +1,94 @@
+#include "mfem.hpp"
+#include "../problems/problems.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+#ifndef IPSOLVER 
+#define IPSOLVER
+
+class InteriorPointSolver
+{
+protected:
+    QPOptContactProblem * optProblem;
+    double rel_tol;
+    int  max_iter;
+    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;
+    Array<int> block_offsetsumlz, block_offsetsuml, block_offsetsx;
+    Vector ml;
+
+    Vector ckSoc;
+    SparseMatrix * Huu = nullptr; 
+    SparseMatrix * Hum = nullptr; 
+    SparseMatrix * Hmu = nullptr; 
+    SparseMatrix * Hmm = nullptr; 
+    SparseMatrix * Wmm = nullptr; 
+    SparseMatrix * Ju  = nullptr; 
+    SparseMatrix * Jm  = nullptr; 
+    SparseMatrix * JmT  = nullptr; 
+    SparseMatrix * JuT = nullptr; 
+    SparseMatrix * Huucl = nullptr;
+    SparseMatrix * HLuucl = nullptr;
+
+    int jOpt;
+    bool converged;
+    
+    int MyRank;
+    bool iAmRoot;
+
+    bool saveLogBarrierIterates;
+
+    int linSolver;
+    double linSolveTol;
+    Array<int> cgnum_iterations;
+
+public:
+    InteriorPointSolver(QPOptContactProblem*);
+    void Mult(const BlockVector& , BlockVector&); // used when the user wants to be aware of bound-constrained variable m >= ml
+    void Mult(const Vector&, Vector &); // useful when the user doesn't need to know about bound-constrained variable   m >= ml, e.g., when m is a slack variable
+    double MaxStepSize(Vector& , Vector& , Vector& , double);
+    double MaxStepSize(Vector& , Vector& , double);
+    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);
+    Array<int> & GetCGIterNumbers() {return cgnum_iterations;}
+    bool GetConverged() const;
+    // 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 SetLinearSolver(int);
+    void SetLinearSolveTol(double);
+    virtual ~InteriorPointSolver();
+};
+
+#endif

miniapps/tribol/contact/ipsolver/ParIPsolver.hpp

@@ -0,0 +1,115 @@
+#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;
+    // QPOptParContactProblemTribol* problem = nullptr;
+    QPOptParContactProblemSingleMesh* 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*);
+    // ParInteriorPointSolver(QPOptParContactProblemTribol*);
+    ParInteriorPointSolver(QPOptParContactProblemSingleMesh*);
+    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 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/tribol/contact/makefile

@@ -0,0 +1,122 @@
+# 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/tribol/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)
+
+CONTACT_SEQ_SRC = problems/problems.cpp problems/problems_util.cpp util/util.cpp ipsolver/IPsolver.cpp
+CONTACT_SEC_OBJ = $(CONTACT_PAR_SRC:.cpp=.o)
+CONTACT_PAR_SRC = $(CONTACT_SEQ_SRC) ipsolver/ParIPsolver.cpp problems/parproblems.cpp problems/parproblems_util.cpp util/mpicomm.cpp 
+CONTACT_PAR_OBJ = $(CONTACT_PAR_SRC:.cpp=.o)
+
+CONTACT_SRC = contact_driver.cpp $(CONTACT_SEQ_SRC)
+CONTACT_OBJ = $(CONTACT_SRC:.cpp=.o)
+
+PCONTACT_SRC = pcontact_driver.cpp $(CONTACT_PAR_SRC)
+PCONTACT_OBJ = $(PCONTACT_SRC:.cpp=.o)
+
+PCONTACT_TRIBOL_SRC = pcontact-tribol.cpp $(CONTACT_PAR_SRC)
+PCONTACT_TRIBOL_OBJ = $(PCONTACT_TRIBOL_SRC:.cpp=.o)
+
+PCONTACT_SINGLEMESH_SRC = pcontact_single_mesh.cpp $(CONTACT_PAR_SRC)
+PCONTACT_SINGLEMESH_OBJ = $(PCONTACT_SINGLEMESH_SRC:.cpp=.o)
+
+SEQ_MINIAPPS = contact_driver
+PAR_MINIAPPS = pcontact_driver pcontact-tribol pcontact_single_mesh
+
+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 $@
+
+util/%.o: $(SRC)util/%.cpp $(wildcard $(SRC)util/%.hpp) $(MFEM_LIB_FILE)\
+ $(CONFIG_MK) | lib-common
+	mkdir -p $(@D)
+	$(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_driver: $(CONTACT_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(CONTACT_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+pcontact_driver: $(PCONTACT_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(PCONTACT_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+pcontact-tribol: $(PCONTACT_TRIBOL_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(PCONTACT_TRIBOL_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+pcontact_single_mesh: $(PCONTACT_SINGLEMESH_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(PCONTACT_SINGLEMESH_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+
+# 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-seq: diffusion
+	@$(call mfem-test,$<,, contact miniapp,)
+pcontact-test-par: pcontact
+	@$(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) $(PCONTACT_OBJ) $(PCONTACT_TRIBOL_OBJ) $(PCONTACT_SINGLEMESH_OBJ) 
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf ParaView 

miniapps/tribol/contact/meshes/block1.mesh

@@ -0,0 +1,103 @@
+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
+#
+
+dimension
+3
+
+elements
+9
+1   5     0      1     3     2     8    9     11    10
+1   5     2      3     5     4    10    11    13    12
+1   5     4      5     7     6    12    13    15    14
+1   5     8      9    11    10    16    17    19    18
+1   5     10    11    13    12    18    19    21    20
+1   5     12    13    15    14    20    21    23    22
+1   5     16    17    19    18    24    25    27    26
+1   5     18    19    21    20    26    27    29    28
+1   5     20    21    23    22    28    29    31    30
+
+
+
+# 0 nothing
+# 1 dirichlet bc
+# 2 contact
+boundary
+30
+1   3    1     0     2     3
+1   3    3     2     4     5
+1   3    5     4     6     7
+1   3    24    25    27    26
+1   3    26    27    29    28
+1   3    28    29    31    30
+2   3    2     0     8    10
+2   3    4     2    10    12
+2   3    6     4    12    14
+2  3    10    8    16    18
+2  3    12    10    18    20
+2  3    14    12    20    22
+2  3    18    16    24    26
+2  3    20    18    26    28
+2  3    22    20    28    30
+3  3    1     3    11    9
+3  3     3     5    13    11
+3  3     5     7    15    13
+3  3    9    11    19    17
+3  3    11    13    21    19
+3  3    13    15    23    21
+3  3    17    19    27    25
+3  3    19    21    29    27
+3  3    21    23    31    29
+1  3     8     0     1    9
+1  3    16     8    9    17
+1  3    24    16    17    25
+1  3     6    14    15     7
+1  3    14    22    23    15
+1  3    22    30    31    23
+
+
+vertices
+32
+3
+   -1.0000         0         0
+         0         0         0
+   -1.0000    0.3000         0
+         0    0.3000         0
+   -1.0000    0.6500         0
+         0    0.6500         0
+   -1.0000    1.0000         0
+         0    1.0000         0
+   -1.0000         0    0.3000
+         0         0    0.3000
+   -1.0000    0.3000    0.3500
+         0    0.3000    0.3500
+   -1.0000    0.6500    0.3000
+         0    0.6500    0.3000
+   -1.0000    1.0000    0.3000
+         0    1.0000    0.3000
+   -1.0000         0    0.6500
+         0         0    0.6500
+   -1.0000    0.3000    0.6500
+         0    0.3000    0.6500
+   -1.0000    0.6500    0.6500
+         0    0.6500    0.6500
+   -1.0000    1.0000    0.6500
+         0    1.0000    0.6500
+   -1.0000         0    1.0000
+         0         0    1.0000
+   -1.0000    0.3000    1.0000
+         0    0.3000    1.0000
+   -1.0000    0.6500    1.0000
+         0    0.6500    1.0000
+   -1.0000    1.0000    1.0000
+         0    1.0000    1.0000

miniapps/tribol/contact/meshes/block2.mesh

@@ -0,0 +1,68 @@
+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
+#
+
+dimension
+3
+
+# 1 nothing
+elements
+4
+1   5     0     1     3     2     6     7    9     8
+1   5     2     3     5     4     8     9    11    10
+1   5     6     7     9     8    12    13    15    14
+1   5     8     9    11    10    14    15    17    16
+
+# 0 nothing
+# 1 dirichlet bc
+# 2 contact
+boundary
+16
+1   3     1     0     2     3
+1   3     3     2     4     5
+1   3    12    13    15    14
+1   3    14    15    17    16
+3   3     2     0     6     8
+3   3     4     2     8    10
+3   3     8     6    12    14
+3   3    10     8    14    16
+2   3     1     3     9     7
+2   3     3     5    11    9
+2   3     7    9    15    13
+2   3     9    11    17    15
+1   3     6     0     1     7
+1   3    12     6     7    13
+1   3     4    10    11     5
+1   3    10    16    17    11
+
+vertices
+18
+3
+         0    0.2464    0.2464
+    0.5071    0.2464    0.2464
+         0    0.5000    0.2464
+    0.5071    0.5000    0.2464
+         0    0.7536    0.2464
+    0.5071    0.7536    0.2464
+         0    0.2464    0.5000
+    0.5071    0.2464    0.5000
+         0    0.5000    0.5000
+    0.5071    0.5000    0.5000
+         0    0.7536    0.5000
+    0.5071    0.7536    0.5000
+         0    0.2464    0.7536
+    0.5071    0.2464    0.7536
+         0    0.5000    0.7536
+    0.5071    0.5000    0.7536
+         0    0.7536    0.7536
+    0.5071    0.7536    0.7536

miniapps/tribol/contact/meshes/merged.mesh

@@ -0,0 +1,136 @@
+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
+13
+1 5 0 1 3 2 8 9 11 10
+1 5 2 3 5 4 10 11 13 12
+1 5 4 5 7 6 12 13 15 14
+1 5 8 9 11 10 16 17 19 18
+1 5 10 11 13 12 18 19 21 20
+1 5 12 13 15 14 20 21 23 22
+1 5 16 17 19 18 24 25 27 26
+1 5 18 19 21 20 26 27 29 28
+1 5 20 21 23 22 28 29 31 30
+2 5 32 33 35 34 38 39 41 40
+2 5 34 35 37 36 40 41 43 42
+2 5 38 39 41 40 44 45 47 46
+2 5 40 41 43 42 46 47 49 48
+
+boundary
+46
+1 3 1 0 2 3
+1 3 3 2 4 5
+1 3 5 4 6 7
+1 3 24 25 27 26
+1 3 26 27 29 28
+1 3 28 29 31 30
+2 3 2 0 8 10
+2 3 4 2 10 12
+2 3 6 4 12 14
+2 3 10 8 16 18
+2 3 12 10 18 20
+2 3 14 12 20 22
+2 3 18 16 24 26
+2 3 20 18 26 28
+2 3 22 20 28 30
+3 3 1 3 11 9
+3 3 3 5 13 11
+3 3 5 7 15 13
+3 3 9 11 19 17
+3 3 11 13 21 19
+3 3 13 15 23 21
+3 3 17 19 27 25
+3 3 19 21 29 27
+3 3 21 23 31 29
+1 3 8 0 1 9
+1 3 16 8 9 17
+1 3 24 16 17 25
+1 3 6 14 15 7
+1 3 14 22 23 15
+1 3 22 30 31 23
+5 3 33 32 34 35
+5 3 35 34 36 37
+5 3 44 45 47 46
+5 3 46 47 49 48
+4 3 34 32 38 40
+4 3 36 34 40 42
+4 3 40 38 44 46
+4 3 42 40 46 48
+6 3 33 35 41 39
+6 3 35 37 43 41
+6 3 39 41 47 45
+6 3 41 43 49 47
+5 3 38 32 33 39
+5 3 44 38 39 45
+5 3 36 42 43 37
+5 3 42 48 49 43
+
+vertices
+50
+3
+-1 0 0
+0 0 0
+-1 0.3 0
+0 0.3 0
+-1 0.65 0
+0 0.65 0
+-1 1 0
+0 1 0
+-1 0 0.3
+0 0 0.3
+-1 0.3 0.35
+0 0.3 0.35
+-1 0.65 0.3
+0 0.65 0.3
+-1 1 0.3
+0 1 0.3
+-1 0 0.65
+0 0 0.65
+-1 0.3 0.65
+0 0.3 0.65
+-1 0.65 0.65
+0 0.65 0.65
+-1 1 0.65
+0 1 0.65
+-1 0 1
+0 0 1
+-1 0.3 1
+0 0.3 1
+-1 0.65 1
+0 0.65 1
+-1 1 1
+0 1 1
+0 0.14577095 0.44389563
+0.5071 0.14577095 0.44389563
+0 0.35093766 0.29483329
+0.5071 0.35093766 0.29483329
+0 0.55610437 0.14577095
+0.5071 0.55610437 0.14577095
+0 0.29483329 0.64906234
+0.5071 0.29483329 0.64906234
+0 0.5 0.5
+0.5071 0.5 0.5
+0 0.70516671 0.35093766
+0.5071 0.70516671 0.35093766
+0 0.44389563 0.85422905
+0.5071 0.44389563 0.85422905
+0 0.64906234 0.70516671
+0.5071 0.64906234 0.70516671
+0 0.85422905 0.55610437
+0.5071 0.85422905 0.55610437

miniapps/tribol/contact/meshes/merged_new.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/tribol/contact/meshes/newmesh1.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/tribol/contact/meshes/newmesh2.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/tribol/contact/meshes/newmesh3.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.525 0.325
+0.175 0.525 0.325
+0.35 0.525 0.325
+0 0.7 0.325
+0.175 0.7 0.325
+0.35 0.7 0.325
+0 0.875 0.325
+0.175 0.875 0.325
+0.35 0.875 0.325
+0 0.525 0.5
+0.175 0.525 0.5
+0.35 0.525 0.5
+0 0.7 0.5
+0.175 0.7 0.5
+0.35 0.7 0.5
+0 0.875 0.5
+0.175 0.875 0.5
+0.35 0.875 0.5
+0 0.525 0.675
+0.175 0.525 0.675
+0.35 0.525 0.675
+0 0.7 0.675
+0.175 0.7 0.675
+0.35 0.7 0.675
+0 0.875 0.675
+0.175 0.875 0.675
+0.35 0.875 0.675

miniapps/tribol/contact/meshes/rotatedblock2.mesh

@@ -0,0 +1,70 @@
+
+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
+#
+
+dimension
+3
+
+# 1 nothing
+elements
+4
+1   5     0     1     3     2     6     7    9     8
+1   5     2     3     5     4     8     9    11    10
+1   5     6     7     9     8    12    13    15    14
+1   5     8     9    11    10    14    15    17    16
+
+# 0 nothing
+# 1 dirichlet bc
+# 2 contact
+boundary
+16
+1   3     1     0     2     3
+1   3     3     2     4     5
+1   3    12    13    15    14
+1   3    14    15    17    16
+3   3     2     0     6     8
+3   3     4     2     8    10
+3   3     8     6    12    14
+3   3    10     8    14    16
+2   3     1     3     9     7
+2   3     3     5    11    9
+2   3     7    9    15    13
+2   3     9    11    17    15
+1   3     6     0     1     7
+1   3    12     6     7    13
+1   3     4    10    11     5
+1   3    10    16    17    11
+
+vertices
+18
+3
+
+0.000000000000 0.145770950245 0.443895630208
+0.507100000000 0.145770950245 0.443895630208
+0.000000000000 0.350937660019 0.294833290227
+0.507100000000 0.350937660019 0.294833290227
+0.000000000000 0.556104369792 0.145770950245
+0.507100000000 0.556104369792 0.145770950245
+0.000000000000 0.294833290227 0.649062339981
+0.507100000000 0.294833290227 0.649062339981
+0.000000000000 0.500000000000 0.500000000000
+0.507100000000 0.500000000000 0.500000000000
+0.000000000000 0.705166709773 0.350937660019
+0.507100000000 0.705166709773 0.350937660019
+0.000000000000 0.443895630208 0.854229049755
+0.507100000000 0.443895630208 0.854229049755
+0.000000000000 0.649062339981 0.705166709773
+0.507100000000 0.649062339981 0.705166709773
+0.000000000000 0.854229049755 0.556104369792
+0.507100000000 0.854229049755 0.556104369792

miniapps/tribol/contact/pcontact-tribol.cpp

@@ -0,0 +1,180 @@
+//                               Parallel contact example
+//
+// Compile with: make pcontact_driver
+// sample run
+// mpirun -np 6 ./pcontact_driver -sr 2 -pr 2
+
+#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();
+   // 1. Parse command-line options.
+   const char *mesh_file = "meshes/merged.mesh";
+   int order = 1;
+   int sref = 0;
+   int pref = 0;
+   Array<int> attr;
+   Array<int> m_attr;
+   bool visualization = true;
+   bool paraview = false;
+   double linsolvertol = 1e-6;
+   int relax_type = 8;
+   double optimizer_tol = 1e-6;
+   int optimizer_maxit = 10;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   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(&linsolvertol, "-stol", "--solver-tol",
+                  "Linear Solver Tolerance.");
+   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(&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);
+   }
+
+   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();
+   }
+
+   MFEM_VERIFY(pmesh->GetNE(), "Empty partition pmesh");
+
+   Array<int> ess_bdr_attr;
+   ess_bdr_attr.Append(2);
+   ess_bdr_attr.Append(6);
+   ParElasticityProblem * prob = new ParElasticityProblem(pmesh,ess_bdr_attr,
+                                                          order);
+
+   Vector lambda(prob->GetMesh()->attributes.Max()); lambda = 57.6923076923;
+   Vector mu(prob->GetMesh()->attributes.Max()); mu = 38.4615384615;
+   prob->SetLambda(lambda); prob->SetMu(mu);
+
+#ifdef MFEM_USE_TRIBOL
+   ParContactProblemTribol contact_tribol(prob);
+   QPOptParContactProblemTribol qpopt(&contact_tribol);
+   int numconstr = contact_tribol.GetGlobalNumConstraints();
+   ParInteriorPointSolver optimizer(&qpopt);
+   optimizer.SetTol(optimizer_tol);
+   optimizer.SetMaxIter(optimizer_maxit);
+
+   int linsolver = 2;
+   optimizer.SetLinearSolver(linsolver);
+   optimizer.SetLinearSolveTol(linsolvertol);
+   optimizer.SetLinearSolveRelaxType(relax_type);
+   ParGridFunction x = prob->GetDisplacementGridFunction();
+   Vector x0 = x.GetTrueVector();
+   int ndofs = x0.Size();
+   Vector xf(ndofs); xf = 0.0;
+   optimizer.Mult(x0, xf);
+   double Einitial = contact_tribol.E(x0);
+   double Efinal = contact_tribol.E(xf);
+   Array<int> & CGiterations = optimizer.GetCGIterNumbers();
+   int gndofs = prob->GetGlobalNumDofs();
+   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 = " << CGiterations.Size() << endl
+                ;
+      mfem::out << " CG iteration numbers           = " ;
+      CGiterations.Print(mfem::out, CGiterations.Size());
+   }
+
+   // MFEM_VERIFY(optimizer.GetConverged(),
+   //             "Interior point solver did not converge.");
+
+
+   if (visualization || paraview)
+   {
+      ParFiniteElementSpace * fes = prob->GetFESpace();
+      ParMesh * pmesh = fes->GetParMesh();
+
+      Vector X_new(xf.GetData(),fes->GetTrueVSize());
+
+      ParGridFunction x_gf(fes);
+
+      x_gf.SetFromTrueDofs(X_new);
+
+      pmesh->MoveNodes(x_gf);
+
+      if (paraview)
+      {
+         ParaViewDataCollection paraview_dc("QPContactBodyTribol", pmesh);
+         paraview_dc.SetPrefixPath("ParaView");
+         paraview_dc.SetLevelsOfDetail(1);
+         paraview_dc.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc.SetHighOrderOutput(true);
+         paraview_dc.SetCycle(0);
+         paraview_dc.SetTime(0.0);
+         paraview_dc.RegisterField("Body", &x_gf);
+         paraview_dc.Save();
+      }
+
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int visport = 19916;
+
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n"
+                  << "solution\n" << *pmesh << x_gf << flush;
+      }
+   }
+#endif
+
+   delete prob;
+   delete pmesh;
+   delete mesh;
+   return 0;
+}

miniapps/tribol/contact/pcontact_driver.cpp

@@ -0,0 +1,260 @@
+//                               Parallel contact example
+//
+// Compile with: make pcontact_driver
+// sample run
+// mpirun -np 6 ./pcontact_driver -sr 2 -pr 2
+
+#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();
+   // 1. Parse command-line options.
+   const char *mesh_file = "meshes/merged.mesh";
+   int order = 1;
+   int sref = 0;
+   int pref = 0;
+   Array<int> attr;
+   Array<int> m_attr;
+   bool visualization = true;
+   bool paraview = false;
+   double linsolvertol = 1e-6;
+   int relax_type = 8;
+   double optimizer_tol = 1e-6;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   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(&linsolvertol, "-stol", "--solver-tol",
+                  "Linear Solver Tolerance.");
+   args.AddOption(&optimizer_tol, "-otol", "--optimizer-tol",
+                  "Interior Point Solver Tolerance.");
+   args.AddOption(&relax_type, "-rt", "--relax-type",
+                  "Selection of Smoother for AMG");
+   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);
+   }
+
+   Mesh * merged_mesh = new Mesh(mesh_file,1);
+
+
+   Array<int> attr1; attr1.Append(1);
+   Array<int> attr2; attr2.Append(2);
+   Mesh * mesh1 = new Mesh(SubMesh::CreateFromDomain(*merged_mesh,attr1));
+   Mesh * mesh2 = new Mesh(SubMesh::CreateFromDomain(*merged_mesh,attr2));
+
+   for (int i = 0; i<sref; i++)
+   {
+      mesh1->UniformRefinement();
+      mesh2->UniformRefinement();
+   }
+   for (int i = 0; i<mesh1->GetNE(); i++)
+   {
+      mesh1->SetAttribute(i,1);
+   }
+   mesh1->SetAttributes();
+   for (int i = 0; i<mesh2->GetNE(); i++)
+   {
+      mesh2->SetAttribute(i,2);
+   }
+   mesh2->SetAttributes();
+
+   ParMesh * pmesh1 = new ParMesh(MPI_COMM_WORLD,*mesh1);
+   ParMesh * pmesh2 = new ParMesh(MPI_COMM_WORLD,*mesh2);
+
+   for (int i = 0; i<pref; i++)
+   {
+      pmesh1->UniformRefinement();
+      pmesh2->UniformRefinement();
+   }
+
+   MFEM_VERIFY(pmesh1->GetNE(), "Empty partition mesh1");
+   MFEM_VERIFY(pmesh2->GetNE(), "Empty partition mesh2");
+
+
+   Array<int> ess_bdr_attr1; ess_bdr_attr1.Append(2);
+   Array<int> ess_bdr_attr2; ess_bdr_attr2.Append(6);
+   ParElasticityProblem * prob1 = new ParElasticityProblem(pmesh1,ess_bdr_attr1,
+                                                           order);
+   ParElasticityProblem * prob2 = new ParElasticityProblem(pmesh2,ess_bdr_attr2,
+                                                           order);
+
+
+   Vector lambda1(prob1->GetMesh()->attributes.Max()); lambda1 = 57.6923076923;
+   Vector mu1(prob1->GetMesh()->attributes.Max()); mu1 = 38.4615384615;
+   Vector lambda2(prob2->GetMesh()->attributes.Max()); lambda2 = 57.6923076923;
+   Vector mu2(prob2->GetMesh()->attributes.Max()); mu2 = 38.4615384615;
+
+   prob1->SetLambda(lambda1); prob1->SetMu(mu1);
+   prob2->SetLambda(lambda2); prob2->SetMu(mu2);
+
+   ParContactProblem contact(prob1,prob2);
+   QPOptParContactProblem qpopt(&contact);
+   int numconstr = contact.GetGlobalNumConstraints();
+
+   ParInteriorPointSolver optimizer(&qpopt);
+
+   optimizer.SetTol(optimizer_tol);
+   optimizer.SetMaxIter(50);
+
+   int linsolver = 2;
+   optimizer.SetLinearSolver(linsolver);
+   optimizer.SetLinearSolveTol(linsolvertol);
+   optimizer.SetLinearSolveRelaxType(relax_type);
+
+   ParGridFunction x1 = prob1->GetDisplacementGridFunction();
+   ParGridFunction x2 = prob2->GetDisplacementGridFunction();
+
+   int ndofs1 = prob1->GetNumTDofs();
+   int ndofs2 = prob2->GetNumTDofs();
+   int gndofs1 = prob1->GetGlobalNumDofs();
+   int gndofs2 = prob2->GetGlobalNumDofs();
+   int ndofs = ndofs1 + ndofs2;
+
+   Vector X1 = x1.GetTrueVector();
+   Vector X2 = x2.GetTrueVector();
+
+   Vector x0(ndofs); x0 = 0.0;
+   x0.SetVector(X1,0);
+   x0.SetVector(X2,X1.Size());
+
+   Vector xf(ndofs); xf = 0.0;
+   optimizer.Mult(x0, xf);
+   Vector & g = contact.GetGapFunction();
+
+   double Einitial = contact.E(x0);
+   double Efinal = contact.E(xf);
+   Array<int> & CGiterations = optimizer.GetCGIterNumbers();
+   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        = " << gndofs1 + gndofs2 << endl;
+      mfem::out << " Global number of constraints = " << numconstr << endl;
+      mfem::out << " CG iteration numbers         = " ;
+      CGiterations.Print(mfem::out, CGiterations.Size());
+   }
+
+   MFEM_VERIFY(optimizer.GetConverged(),
+               "Interior point solver did not converge.");
+
+
+   if (visualization || paraview)
+   {
+      ParFiniteElementSpace * fes1 = prob1->GetFESpace();
+      ParFiniteElementSpace * fes2 = prob2->GetFESpace();
+
+      ParMesh * pmesh_1 = fes1->GetParMesh();
+      ParMesh * pmesh_2 = fes2->GetParMesh();
+
+      Vector X1_new(xf.GetData(),fes1->GetTrueVSize());
+      Vector X2_new(&xf.GetData()[fes1->GetTrueVSize()],fes2->GetTrueVSize());
+
+      ParGridFunction x1_gf(fes1);
+      ParGridFunction x2_gf(fes2);
+
+      x1_gf.SetFromTrueDofs(X1_new);
+      x2_gf.SetFromTrueDofs(X2_new);
+
+      pmesh_1->MoveNodes(x1_gf);
+      pmesh_2->MoveNodes(x2_gf);
+
+      if (paraview)
+      {
+         ParaViewDataCollection paraview_dc1("QPContactBody1", pmesh_1);
+         paraview_dc1.SetPrefixPath("ParaView");
+         paraview_dc1.SetLevelsOfDetail(1);
+         paraview_dc1.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc1.SetHighOrderOutput(true);
+         paraview_dc1.SetCycle(0);
+         paraview_dc1.SetTime(0.0);
+         paraview_dc1.RegisterField("Body1", &x1_gf);
+         paraview_dc1.Save();
+
+         ParaViewDataCollection paraview_dc2("QPContactBody2", pmesh_2);
+         paraview_dc2.SetPrefixPath("ParaView");
+         paraview_dc2.SetLevelsOfDetail(1);
+         paraview_dc2.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc2.SetHighOrderOutput(true);
+         paraview_dc2.SetCycle(0);
+         paraview_dc2.SetTime(0.0);
+         paraview_dc2.RegisterField("Body2", &x2_gf);
+         paraview_dc2.Save();
+      }
+
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int visport = 19916;
+
+         // {
+         //    socketstream sol_sock1(vishost, visport);
+         //    sol_sock1.precision(8);
+         //    sol_sock1 << "parallel " << num_procs << " " << myid << "\n"
+         //              << "solution\n" << *pmesh_1 << x1_gf << flush;
+         // }
+         // {
+         //    socketstream sol_sock2(vishost, visport);
+         //    sol_sock2.precision(8);
+         //    sol_sock2 << "parallel " << num_procs << " " << myid << "\n"
+         //              << "solution\n" << *pmesh_2 << x2_gf << flush;
+         // }
+
+         {
+            socketstream sol_sock(vishost, visport);
+            sol_sock.precision(8);
+            sol_sock << "parallel " << 2*num_procs << " " << myid << "\n"
+                     << "solution\n" << *pmesh_1 << x1_gf << flush;
+         }
+         {
+            socketstream sol_sock(vishost, visport);
+            sol_sock.precision(8);
+            sol_sock << "parallel " << 2*num_procs << " " << myid+num_procs << "\n"
+                     << "solution\n" << *pmesh_2 << x2_gf << flush;
+         }
+      }
+   }
+
+   delete prob2;
+   delete prob1;
+   delete pmesh2;
+   delete pmesh1;
+   // delete mesh1;
+   // delete mesh2;
+
+   return 0;
+}

miniapps/tribol/contact/pcontact_single_mesh.cpp

@@ -0,0 +1,379 @@
+//                               Parallel contact example
+//
+// Compile with: make pcontact_driver
+// sample run
+// mpirun -np 6 ./pcontact_driver -sr 2 -pr 2
+// mpirun -np 8 ./pcontact_single_mesh -ls 2 -sr 5 -omaxit 20  -otol 1e-6 -tribol -rt 18
+#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 = 0;
+   int pref = 0;
+   Array<int> attr;
+   Array<int> m_attr;
+   bool visualization = true;
+   bool paraview = false;
+   double linsolverrtol = 1e-10;
+   double linsolveratol = 1e-12;
+   int relax_type = 8;
+   double optimizer_tol = 1e-6;
+   int optimizer_maxit = 20;
+   bool enable_tribol = true;
+   int linsolver = 2; // PCG  - AMG
+   bool elast = false;
+   bool nocontact = false;
+   int testNo = 4; // 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(&linsolverrtol, "-srtol", "--solver-rel-tol",
+                  "Linear Solver Relative Tolerance.");
+   args.AddOption(&linsolveratol, "-satol", "--solver-abs-tol",
+                  "Linear Solver Abs Tolerance.");
+   args.AddOption(&enable_tribol, "-tribol", "--tribol", "-no-tribol",
+                  "--no-tribol",
+                  "Enable or disable Tribol interface.");
+   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(&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 cont 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.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/merged_new.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 41:
+         mesh_file = "meshes/newmesh1.mesh";
+         break;
+      case 5:
+         mesh_file = "meshes/Test5.mesh";
+         break;
+      case 51:
+         mesh_file = "meshes/Test5mod.mesh";
+         break;
+      case 6:
+         mesh_file = "meshes/Test6.mesh";
+         break;
+      case 61:
+         // Something wrong with this mesh
+         mesh_file = "meshes/Test6mod.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();
+   }
+
+   Array<int> part;
+   Array<int> attr1, attr2;
+   attr1.Append(1);
+   attr2.Append(2);
+   // SubMesh mesh1 = SubMesh::CreateFromDomain(*mesh,attr1);
+   // SubMesh mesh2 = SubMesh::CreateFromDomain(*mesh,attr2);
+
+   // Array<int> part1(mesh1.GeneratePartitioning(num_procs),mesh1.GetNE());
+   // Array<int> part2(mesh2.GeneratePartitioning(num_procs),mesh2.GetNE());
+
+   // part.Append(part1);
+   // part.Append(part2);
+
+   // ParMesh * pmesh = new ParMesh(MPI_COMM_WORLD,*mesh,part.GetData());
+   ParMesh * pmesh = new ParMesh(MPI_COMM_WORLD,*mesh);
+
+   for (int i = 0; i<pref; i++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // MFEM_VERIFY(pmesh->GetNE(), "Empty partition pmesh");
+
+   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
+   {
+      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());
+   lambda[1] = 0.0;
+   lambda[0] = 0.499/(1.499*0.002);
+
+   Vector mu(prob->GetMesh()->attributes.Max());
+   mu[1] = 500;
+   mu[0] = 1./(2*1.499);
+
+   if (testNo == 6 || testNo == 61 )
+   {
+      lambda = (1000*0.3)/(1.3*0.4);
+      mu = 500/(1.3);
+   }
+
+   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(-10.0);
+
+   std::set<int> mortar_attr;
+   std::set<int> nonmortar_attr;
+
+   if (testNo == 6 || testNo == 61)
+   {
+      /* material 1: e = 1000
+                     nu = 0.3
+         material 2: e = 1000
+                     nu = 0.3
+         material 3: e = 1000
+                     nu = 0.3
+         Attr    value
+         1      Dirichlet y = 0
+         2      Dirichlet z = 0
+         3      Neuman = 1.0
+         4      Dirichlet x = 0
+         5      Dirichlet x=y=z=0
+         6      Mortar
+         7      NonMortar
+         8      NonMortar
+         9      Mortar
+      */
+      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 == -1 || testNo == 41)
+      {
+         ess_values[0] = 0.1;
+      }
+      else
+      {
+         ess_values[2] = 0.7;
+      }
+      essbdr_attr = 2;
+      ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      essbdr_attr = 6;
+      ess_values = 0.0; ess_bdr = 0; ess_bdr[essbdr_attr - 1] = 1;
+      prob->SetDisplacementDirichletData(ess_values, ess_bdr);
+      mortar_attr.insert(3);
+      nonmortar_attr.insert(4);
+   }
+
+
+
+   ParContactProblemSingleMesh contact(prob, mortar_attr, nonmortar_attr,
+                                       enable_tribol);
+   QPOptParContactProblemSingleMesh qpopt(&contact);
+   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();
+   Vector x0 = x.GetTrueVector();
+   int ndofs = x0.Size();
+   Vector xf(ndofs); xf = 0.0;
+   optimizer.Mult(x0, xf);
+   double Einitial = contact.E(x0);
+   double Efinal = contact.E(xf);
+   Array<int> & CGiterations = optimizer.GetCGIterNumbers();
+   int gndofs = prob->GetGlobalNumDofs();
+   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());
+      }
+
+   }
+
+   // MFEM_VERIFY(optimizer.GetConverged(),
+   //             "Interior point solver did not converge.");
+
+
+   if (visualization || paraview)
+   {
+      ParFiniteElementSpace * fes = prob->GetFESpace();
+      ParMesh * pmesh = fes->GetParMesh();
+
+      Vector X_new(xf.GetData(),fes->GetTrueVSize());
+
+      ParGridFunction x_gf(fes);
+
+      x_gf.SetFromTrueDofs(X_new);
+      // x_gf*=-1.0;
+
+      pmesh->MoveNodes(x_gf);
+
+      if (paraview)
+      {
+         std::ostringstream paraview_file_name;
+         paraview_file_name << "QPContactBody"
+                            <<"_Tribol_" << (int)enable_tribol
+                            << "_par_ref_" << pref
+                            << "_ser_ref_" << sref;
+         ParaViewDataCollection paraview_dc(paraview_file_name.str(), pmesh);
+         paraview_dc.SetPrefixPath("ParaView");
+         paraview_dc.SetLevelsOfDetail(1);
+         paraview_dc.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc.SetHighOrderOutput(true);
+         paraview_dc.SetCycle(0);
+         paraview_dc.SetTime(0.0);
+         paraview_dc.RegisterField("Body", &x_gf);
+         paraview_dc.Save();
+      }
+
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int visport = 19916;
+
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n"
+                  << "solution\n" << *pmesh << x_gf << flush;
+      }
+   }
+
+   delete prob;
+   delete pmesh;
+   delete mesh;
+   return 0;
+}

miniapps/tribol/contact/problems/parproblems.hpp

@@ -0,0 +1,492 @@
+
+#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;
+   ParGridFunction x;
+   HypreParMatrix A;
+   Vector B,X;
+   ConstantCoefficient pressure_cf;
+   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 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 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 fes;
+      delete fec;
+      if (own_mesh)
+      {
+         delete pmesh;
+      }
+   }
+};
+
+
+class ParContactProblem
+{
+private:
+   MPI_Comm comm;
+   int numprocs;
+   int myid;
+   ParElasticityProblem * prob1 = nullptr;
+   ParElasticityProblem * prob2 = nullptr;
+   ParFiniteElementSpace * vfes1 = nullptr;
+   ParFiniteElementSpace * vfes2 = nullptr;
+   int dim;
+   GridFunction nodes0;
+   GridFunction *nodes1 = nullptr;
+   std::set<int> contact_vertices;
+   bool recompute = true;
+   bool compute_hessians = true;
+   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> globalvertices1;
+   Array<int> globalvertices2;
+   Array<int> vertices2;
+   Array<int> vertices1;
+
+protected:
+   int npoints=0;
+   int gnpoints=0;
+   int nv, gnv;
+   HypreParMatrix * K = nullptr;
+   BlockVector *B = nullptr;
+   Vector gapv;
+   HypreParMatrix * M=nullptr;
+   Array<HypreParMatrix*> dM;
+   void ComputeContactVertices();
+
+public:
+   ParContactProblem(ParElasticityProblem * prob1_, ParElasticityProblem * prob2_);
+
+   ParElasticityProblem * GetElasticityProblem1() {return prob1;}
+   ParElasticityProblem * GetElasticityProblem2() {return prob2;}
+   MPI_Comm GetComm() {return comm;}
+   int GetNumDofs() {return K->Height();}
+   int GetGlobalNumDofs() {return K->GetGlobalNumRows();}
+   int GetNumContraints() {return npoints;}
+   int GetGlobalNumConstraints() {return gnpoints;}
+   
+   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;}
+   Array<HypreParMatrix*> & GetHessian() {return dM;}
+   void ComputeGapFunctionAndDerivatives(const Vector & displ1, const Vector &displ2);
+
+   double E(const Vector & d);
+   void DdE(const Vector &d, Vector &gradE);
+   HypreParMatrix* DddE(const Vector &d);
+   void g(const Vector &d, Vector &gd, bool compute_hessians_ = true);
+   HypreParMatrix* Ddg(const Vector &d);
+   HypreParMatrix* lDddg(const Vector &d, const Vector &l);
+
+   ~ParContactProblem()
+   {
+      delete B;
+      delete K;
+      delete M;
+      for (int i = 0; i<dM.Size(); i++)
+      {
+         delete dM[i];
+      }
+      delete vfes1;
+      delete vfes2;
+   }
+};
+
+
+class QPOptParContactProblem
+{
+private:
+   ParContactProblem * problem = nullptr;
+   int dimU, dimM, dimC;
+   // Array<int> block_offsets;
+   Vector ml;
+   HypreParMatrix * NegId = nullptr;
+public:
+   QPOptParContactProblem(ParContactProblem * problem_);
+   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 * GetElasticityProblem1() {return problem->GetElasticityProblem1();}
+   ParElasticityProblem * GetElasticityProblem2() {return problem->GetElasticityProblem2();}
+
+   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 &);
+   void c(const BlockVector &, Vector &);
+   double CalcObjective(const BlockVector &);
+   void CalcObjectiveGrad(const BlockVector &, BlockVector &);
+   ~QPOptParContactProblem();
+};
+
+
+
+#ifdef MFEM_USE_TRIBOL
+class ParContactProblemTribol
+{
+private:
+   MPI_Comm comm;
+   int numprocs;
+   int myid;
+   ParElasticityProblem * prob = nullptr;
+   ParMesh * pmesh = nullptr;
+   ParFiniteElementSpace * vfes = nullptr;
+   int dim;
+   ParGridFunction *coords = nullptr;
+   Array<int> constraints_starts;
+   // gap function
+   Vector gapv;
+
+   void SetupTribol();
+
+protected:
+   HypreParMatrix * K = nullptr;
+   Vector *B = nullptr;
+   // Gap Jacobian
+   HypreParMatrix * J=nullptr;
+
+public:
+// for now we work on 1 (merged mesh).
+   ParContactProblemTribol(ParElasticityProblem * prob_);
+   MPI_Comm GetComm() {return comm;}
+   int GetNumDofs() {return K->Height();}
+   int GetGlobalNumDofs() {return K->GetGlobalNumRows();}
+   int GetNumContraints() {return J->Height();}
+   int GetGlobalNumConstraints() {return J->GetGlobalNumRows();}
+   Array<int> & GetConstraintsStarts() { return constraints_starts; }
+   Vector & GetGapFunction() {return gapv;}
+   HypreParMatrix * GetJacobian() {return J;}
+   double E(const Vector & d);
+   void DdE(const Vector &d, Vector &gradE);
+   HypreParMatrix* DddE();
+   void g(Vector &gd);
+   HypreParMatrix* Ddg();
+
+   ~ParContactProblemTribol()
+   {
+      delete B;
+      delete K;
+   }
+};
+
+
+class QPOptParContactProblemTribol
+{
+private:
+   ParContactProblemTribol * problem = nullptr;
+   int dimU, dimM, dimC;
+   // Array<int> block_offsets;
+   Vector ml;
+   HypreParMatrix * NegId = nullptr;
+public:
+   QPOptParContactProblemTribol(ParContactProblemTribol * problem_);
+   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();}
+
+   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 &);
+   void c(const BlockVector &, Vector &);
+   double CalcObjective(const BlockVector &);
+   void CalcObjectiveGrad(const BlockVector &, BlockVector &);
+   ~QPOptParContactProblemTribol();
+};
+
+#endif
+
+
+class ParContactProblemSingleMesh
+{
+private:
+   MPI_Comm comm;
+   int numprocs;
+   int myid;
+   ParElasticityProblem * prob = nullptr;
+   ParFiniteElementSpace * vfes = nullptr;
+   int dim;
+   bool recompute = true;
+   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;
+
+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;
+   Array<HypreParMatrix*> dM;
+   void ComputeContactVertices();
+   void SetupTribol();
+   bool enable_tribol = false;
+   std::set<int> mortar_attrs;
+   // plane of top block
+   std::set<int> nonmortar_attrs;
+
+public:
+   ParContactProblemSingleMesh(ParElasticityProblem * prob_, 
+                               const std::set<int> & mortar_attrs_, const std::set<int> & nonmortar_attrs_,
+                               bool enable_tribol_ = false);
+
+   ParElasticityProblem * GetElasticityProblem() {return prob;}
+   MPI_Comm GetComm() {return comm;}
+   int GetNumDofs() {return K->Height();}
+   int GetGlobalNumDofs() {return K->GetGlobalNumRows();}
+   int GetNumContraints() 
+   {
+      if (enable_tribol)
+      {
+         return M->Height();
+      }
+      else
+      {
+         return npoints;
+      }
+   }
+   int GetGlobalNumConstraints() 
+   {
+      if (enable_tribol)
+      {
+         return M->GetGlobalNumRows();
+      }
+      else
+      {
+         return gnpoints;         
+      }
+   }
+   
+   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;}
+   Array<HypreParMatrix*> & GetHessian() {return dM;}
+   void ComputeGapFunctionAndDerivatives(const Vector & displ);
+
+   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;}
+
+   ~ParContactProblemSingleMesh()
+   {
+      delete B;
+      delete K;
+      delete M;
+      for (int i = 0; i<dM.Size(); i++)
+      {
+         delete dM[i];
+      }
+      if (!enable_tribol) delete vfes;
+   }
+};
+
+
+class QPOptParContactProblemSingleMesh
+{
+private:
+   ParContactProblemSingleMesh * problem = nullptr;
+   int dimU, dimM, dimC;
+   Vector ml;
+   HypreParMatrix * NegId = nullptr;
+public:
+   QPOptParContactProblemSingleMesh(ParContactProblemSingleMesh * problem_);
+   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 &);
+   ~QPOptParContactProblemSingleMesh();
+};

miniapps/tribol/contact/problems/parproblems_util.cpp

@@ -0,0 +1,579 @@
+#include "parproblems_util.hpp"
+
+void FindPointsInMesh(Mesh & mesh, const Array<int> & gvert, const Vector & xyz, const Array<int> & s_conn, Array<int>& conn,
+                      Vector & xyz2, Array<int> & s_conn2, Vector& xi, DenseMatrix & coords)
+{
+   const int dim = mesh.Dimension();
+   const int np = xyz.Size() / dim;
+
+   MFEM_VERIFY(np * dim == xyz.Size(), "");
+
+   mesh.EnsureNodes();
+
+   FindPointsGSLIB finder(MPI_COMM_WORLD);
+
+   finder.SetDistanceToleranceForPointsFoundOnBoundary(0.5);
+
+   const double bb_t = 0.5;
+   finder.Setup(mesh, bb_t);
+
+   finder.FindPoints(xyz,mfem::Ordering::byVDIM);
+
+   Array<unsigned int> procs = finder.GetProc();
+
+   /// Return code for each point searched by FindPoints: inside element (0), on
+   /// element boundary (1), or not found (2).
+   Array<unsigned int> codes = finder.GetCode();
+
+   /// Return element number for each point found by FindPoints.
+   Array<unsigned int> elems = finder.GetElem();
+
+   /// Return reference coordinates for each point found by FindPoints.
+   Vector refcrd = finder.GetReferencePosition();
+
+   /// Return distance between the sought and the found point in physical space,
+   /// for each point found by FindPoints.
+   Vector dist = finder.GetDist();
+
+   finder.FreeData();
+
+   MFEM_VERIFY(dist.Size() == np, "");
+   MFEM_VERIFY(refcrd.Size() == np * dim, "");
+   MFEM_VERIFY(elems.Size() == np, "");
+   MFEM_VERIFY(codes.Size() == np, "");
+
+   bool allfound = true;
+   for (auto code : codes)
+      if (code == 2) { allfound = false; }
+
+   MFEM_VERIFY(allfound, "A point was not found");
+
+   // cout << "Maximum distance of projected points: " << dist.Max() << endl;
+
+
+   Array<unsigned int> elems_recv, proc_recv;
+   Vector ref_recv;
+   Vector xyz_recv;
+   Array<int> s_conn_recv;
+
+   MPICommunicator mycomm(MPI_COMM_WORLD, procs);
+   mycomm.Communicate(xyz,xyz_recv,3,mfem::Ordering::byNODES);
+   mycomm.Communicate(elems,elems_recv,1,mfem::Ordering::byVDIM);
+   mycomm.Communicate(refcrd,ref_recv,3,mfem::Ordering::byVDIM);
+   mycomm.Communicate(s_conn,s_conn_recv,1,mfem::Ordering::byVDIM);
+
+   proc_recv = mycomm.GetOriginProcs();
+
+   int np_loc = elems_recv.Size();
+   Array<int> conn_loc(np_loc*4);
+   Vector xi_send(np_loc*(dim-1));
+   for (int i=0; i<np_loc; ++i)
+   {
+      int refFace, refNormal;
+      // int refNormalSide;
+      bool is_interior = -1;
+
+      Vector normal = GetNormalVector(mesh, elems_recv[i],
+                                      ref_recv.GetData() + (i*dim),
+                                      refFace, refNormal, is_interior);
+
+      // continue;
+      int phyFace;
+      if (is_interior)
+      {
+         phyFace = -1; // the id of the face that has the closest point
+         FindSurfaceToProject(mesh, elems_recv[i], phyFace); // seems that this works
+
+         Array<int> cbdrVert;
+         mesh.GetFaceVertices(phyFace, cbdrVert);
+         Vector xs(dim);
+         xs[0] = xyz_recv[i + 0*np_loc];
+         xs[1] = xyz_recv[i + 1*np_loc];
+         xs[2] = xyz_recv[i + 2*np_loc];
+
+         Vector xi_tmp(dim-1);
+         // get nodes!
+
+         GridFunction *nodes = mesh.GetNodes();
+         DenseMatrix coord(4,3);
+         for (int j=0; j<4; j++)
+         {
+            for (int k=0; k<3; k++)
+            {
+               coord(j,k) = (*nodes)[cbdrVert[j]*3+k];
+            }
+         }
+         SlaveToMaster(coord, xs, xi_tmp);
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi_send[i*(dim-1)+j] = xi_tmp[j];
+         }
+         // now get get the projection to the surface
+      }
+      else
+      {
+         Vector faceRefCrd(dim-1);
+         {
+            int fd = 0;
+            for (int j=0; j<dim; ++j)
+            {
+               if (j == refNormal)
+               {
+                  // refNormalSide = (ref_recv[(i*dim) + j] > 0.5); // not used
+               }
+               else
+               {
+                  faceRefCrd[fd] = ref_recv[(i*dim) + j];
+                  fd++;
+               }
+            }
+            MFEM_VERIFY(fd == dim-1, "");
+         }
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi_send[i*(dim-1)+j] = faceRefCrd[j]*2.0 - 1.0;
+         }
+      }
+      // Get the element face
+      Array<int> faces;
+      Array<int> ori;
+      int face;
+
+      if (is_interior)
+      {
+         face = phyFace;
+      }
+      else
+      {
+         mesh.GetElementFaces(elems_recv[i], faces, ori);
+         face = faces[refFace];
+      }
+
+      Array<int> faceVert;
+      mesh.GetFaceVertices(face, faceVert);
+
+      for (int p=0; p<4; p++)
+      {
+         conn_loc[4*i+p] = faceVert[p];
+      }
+   }
+
+   if (0) // for debugging
+   {
+      int sz = xi_send.Size()/2;
+
+      for (int i = 0; i<sz; i++)
+      {
+         mfem::out << "("<<xi_send[i*(dim-1)]<<","<<xi_send[i*(dim-1)+1]<<"): -> ";
+         for (int j = 0; j<4; j++)
+         {
+            double * vc = mesh.GetVertex(conn_loc[4*i+j]);
+            if (j<3)
+            {
+               mfem::out << "("<<vc[0]<<","<<vc[1]<<","<<vc[2]<<"), ";
+            }
+            else
+            {
+               mfem::out << "("<<vc[0]<<","<<vc[1]<<","<<vc[2]<<") \n " << endl;
+            }
+         }
+      }
+   }
+
+   int sz = xi_send.Size()/2;
+   DenseMatrix coordsm(sz*4, dim);
+   for (int i = 0; i<sz; i++)
+   {
+      for (int j = 0; j<4; j++)
+      {
+         for (int k=0; k<dim; k++)
+         {
+            coordsm(i*4+j,k) = mesh.GetVertex(conn_loc[i*4+j])[k];
+         }
+      }
+   }
+
+   // pass global indices for conn_loc
+   for (int i = 0; i<conn_loc.Size(); i++)
+   {
+      conn_loc[i] = gvert[conn_loc[i]];
+   }
+
+   mycomm.UpdateDestinationProcs();
+   mycomm.Communicate(xyz_recv,xyz2,3,mfem::Ordering::byNODES);
+   mycomm.Communicate(xi_send,xi,2,mfem::Ordering::byVDIM);
+   mycomm.Communicate(s_conn_recv,s_conn2,1,mfem::Ordering::byVDIM);
+   mycomm.Communicate(conn_loc,conn,4,mfem::Ordering::byVDIM);
+   mycomm.Communicate(coordsm,coords,4,mfem::Ordering::byVDIM);
+}
+
+
+void FindPointsInMesh(ParMesh & pmesh, const Array<int> & gvert, Array<int> & s_conn, const Vector &x1, Vector & xyz, Array<int>& conn,
+                      Vector& xi, DenseMatrix & coords, bool singlemesh)
+{
+   const int dim = pmesh.Dimension();
+   const int np = xyz.Size() / dim;
+   MFEM_VERIFY(np * dim == xyz.Size(), "");
+
+   pmesh.EnsureNodes();
+
+   ParSubMesh * psubmesh = nullptr;
+   ParMesh * mesh = nullptr;
+   Array<int> elem_map;
+   if (singlemesh)
+   {
+      Array<int> attr; attr.Append(1);
+      psubmesh = new ParSubMesh(ParSubMesh::CreateFromDomain(pmesh,attr));
+      mesh = (ParMesh *)psubmesh;
+   }
+   else
+   {
+      mesh = &pmesh;
+   }
+
+   FindPointsGSLIB finder(MPI_COMM_WORLD);
+
+   finder.SetDistanceToleranceForPointsFoundOnBoundary(0.5);
+
+   const double bb_t = 0.5;
+   MFEM_VERIFY(mesh->GetNE(), "FindPointsGSLIB does not support empty partition");
+   finder.Setup(*mesh, bb_t);
+
+   finder.FindPoints(xyz,mfem::Ordering::byVDIM);
+
+   Array<unsigned int> procs = finder.GetProc();
+
+   /// Return code for each point searched by FindPoints: inside element (0), on
+   /// element boundary (1), or not found (2).
+   Array<unsigned int> codes = finder.GetCode();
+
+   /// Return element number for each point found by FindPoints.
+   Array<unsigned int> elems = finder.GetElem();
+
+   /// Return reference coordinates for each point found by FindPoints.
+   Vector refcrd = finder.GetReferencePosition();
+
+   /// Return distance between the sought and the found point in physical space,
+   /// for each point found by FindPoints.
+   Vector dist = finder.GetDist();
+
+   finder.FreeData();
+
+   MFEM_VERIFY(dist.Size() == np, "");
+   MFEM_VERIFY(refcrd.Size() == np * dim, "");
+   MFEM_VERIFY(elems.Size() == np, "");
+   MFEM_VERIFY(codes.Size() == np, "");
+
+   bool allfound = true;
+   for (auto code : codes)
+      if (code == 2) { allfound = false; }
+
+   MFEM_VERIFY(allfound, "A point was not found");
+
+   // reorder data so that the procs are in ascending order
+   // sort procs and save the permutation
+   std::vector<unsigned int> procs_index(np);
+   std::iota(procs_index.begin(),procs_index.end(),0); //Initializing
+   sort( procs_index.begin(),procs_index.end(), [&](int i,int j){return procs[i]<procs[j];} );
+
+   // map to sorted
+   Array<unsigned int> procs_sorted(np);
+   Array<unsigned int> elems_sorted(np);
+   Vector xyz_sorted(np*dim);
+   Vector refcrd_sorted(np*dim);
+   Array<int> s_conn_sorted(np);
+   for (int i = 0; i<np; i++)
+   {
+      int j = procs_index[i];
+      procs_sorted[i] = procs[j];
+      elems_sorted[i] = elems[j];
+      s_conn_sorted[i] = s_conn[j];
+      for (int d = 0; d<dim; d++)
+      {
+         xyz_sorted(i*dim+d) = xyz(j*dim+d);
+         refcrd_sorted(i*dim+d) = refcrd(j*dim+d);
+      }
+   }
+
+   Array<unsigned int> elems_recv, proc_recv;
+   xyz = xyz_sorted;
+   s_conn = s_conn_sorted;
+   Vector ref_recv;
+   Vector xyz_recv;
+
+   MPICommunicator mycomm(MPI_COMM_WORLD, procs_sorted);
+   mycomm.Communicate(xyz_sorted,xyz_recv,3,mfem::Ordering::byVDIM);
+   mycomm.Communicate(elems_sorted,elems_recv,1,mfem::Ordering::byVDIM);
+   mycomm.Communicate(refcrd_sorted,ref_recv,3,mfem::Ordering::byVDIM);
+
+   proc_recv = mycomm.GetOriginProcs();
+   int np_loc = elems_recv.Size();
+   if (singlemesh)
+   {
+      elem_map = psubmesh->GetParentElementIDMap();
+      for (int i = 0; i< np_loc; i++)
+      {
+         elems_recv[i] = elem_map[elems_recv[i]];
+      }
+      delete mesh;
+   }
+
+
+   Array<int> conn_loc(np_loc*4);
+   Vector xi_send(np_loc*(dim-1));
+   for (int i=0; i<np_loc; ++i)
+   {
+      int refFace, refNormal; 
+      // int refNormalSide;
+      bool is_interior = -1;
+      Vector normal = GetNormalVector(pmesh, elems_recv[i],
+                                      ref_recv.GetData() + (i*dim),
+                                      refFace, refNormal, is_interior);
+
+      // continue;
+      int phyFace;
+      if (is_interior)
+      {
+         phyFace = -1; // the id of the face that has the closest point
+         FindSurfaceToProject(pmesh, elems_recv[i], phyFace); // seems that this works
+
+         Array<int> cbdrVert;
+         pmesh.GetFaceVertices(phyFace, cbdrVert);
+         Vector xs(dim);
+         xs[0] = xyz_recv[i*dim + 0];
+         xs[1] = xyz_recv[i*dim + 1];
+         xs[2] = xyz_recv[i*dim + 2];
+
+         Vector xi_tmp(dim-1);
+         // get nodes!
+
+         GridFunction *nodes = pmesh.GetNodes();
+         DenseMatrix coord(4,3);
+         for (int j=0; j<4; j++)
+         {
+            for (int k=0; k<3; k++)
+            {
+               coord(j,k) = (*nodes)[cbdrVert[j]*3+k];
+            }
+         }
+         SlaveToMaster(coord, xs, xi_tmp);
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi_send[i*(dim-1)+j] = xi_tmp[j];
+         }
+         // now get the projection to the surface
+      }
+      else
+      {
+         Vector faceRefCrd(dim-1);
+         {
+            int fd = 0;
+            for (int j=0; j<dim; ++j)
+            {
+               if (j == refNormal)
+               {
+                  // refNormalSide = (ref_recv[(i*dim) + j] > 0.5); // not used
+               }
+               else
+               {
+                  faceRefCrd[fd] = ref_recv[(i*dim) + j];
+                  fd++;
+               }
+            }
+            MFEM_VERIFY(fd == dim-1, "");
+         }
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi_send[i*(dim-1)+j] = faceRefCrd[j]*2.0 - 1.0;
+         }
+      }
+      // Get the element face
+      Array<int> faces;
+      Array<int> ori;
+      int face;
+
+      if (is_interior)
+      {
+         face = phyFace;
+      }
+      else
+      {
+         pmesh.GetElementFaces(elems_recv[i], faces, ori);
+         face = faces[refFace];
+      }
+
+      Array<int> faceVert;
+      pmesh.GetFaceVertices(face, faceVert);
+
+      for (int p=0; p<4; p++)
+      {
+         conn_loc[4*i+p] = faceVert[p];
+      }
+   }
+
+   if (0) // for debugging
+   {
+      int sz = xi_send.Size()/2;
+
+      for (int i = 0; i<sz; i++)
+      {
+         mfem::out << "("<<xi_send[i*(dim-1)]<<","<<xi_send[i*(dim-1)+1]<<"): -> ";
+         for (int j = 0; j<4; j++)
+         {
+            double * vc = pmesh.GetVertex(conn_loc[4*i+j]);
+            if (j<3)
+            {
+               mfem::out << "("<<vc[0]<<","<<vc[1]<<","<<vc[2]<<"), ";
+            }
+            else
+            {
+               mfem::out << "("<<vc[0]<<","<<vc[1]<<","<<vc[2]<<") \n " << endl;
+            }
+         }
+      }
+   }
+
+   int sz = xi_send.Size()/2;
+   DenseMatrix coordsm(sz*4, dim);
+   for (int i = 0; i<sz; i++)
+   {
+      for (int j = 0; j<4; j++)
+      {
+         for (int k=0; k<dim; k++)
+         {
+            coordsm(i*4+j,k) = pmesh.GetVertex(conn_loc[i*4+j])[k]+x1[dim*conn_loc[i*4+j]+k];
+         }
+      }
+   }
+
+   // pass global indices for conn_loc
+   for (int i = 0; i<conn_loc.Size(); i++)
+   {
+      conn_loc[i] = gvert[conn_loc[i]];
+   }
+
+   mycomm.UpdateDestinationProcs();
+   mycomm.Communicate(xi_send,xi,2,mfem::Ordering::byVDIM);
+   mycomm.Communicate(conn_loc,conn,4,mfem::Ordering::byVDIM);
+   mycomm.Communicate(coordsm,coords,4,mfem::Ordering::byVDIM);
+
+
+}
+
+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/tribol/contact/problems/parproblems_util.hpp

@@ -0,0 +1,31 @@
+
+#include "mfem.hpp"
+#include "problems_util.hpp"
+#include "../util/mpicomm.hpp"
+#include "axom/slic.hpp"
+
+#include "tribol/interface/tribol.hpp"
+#include "tribol/interface/mfem_tribol.hpp"
+// Coordinates in xyz are assumed to be ordered as [X, Y, Z]
+// where X is the list of x-coordinates for all points and so on.
+// conn: connectivity of the target surface elements
+// xi: surface reference cooridnates for the cloest point, involves a linear transformation from [0,1] to [-1,1]
+void FindPointsInMesh(Mesh & mesh, const Array<int> & gvert, const Vector & xyz, const Array<int> & s_conn, Array<int>& conn,
+                      Vector & xyz2, Array<int> & s_conn2, Vector& xi, DenseMatrix & coords);
+
+// somewhat simplified version of the above
+void FindPointsInMesh(ParMesh & mesh, const Array<int> & gvert, Array<int> & s_conn, const Vector &x1, Vector & xyz, Array<int>& conn,
+                      Vector& xi, DenseMatrix & coords, bool singlemesh = false);                   
+
+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/tribol/contact/problems/problems.cpp

@@ -0,0 +1,367 @@
+#include "problems.hpp"
+
+
+void ElasticityProblem::Init()
+{
+   int dim = mesh->Dimension();
+   fec = new H1_FECollection(order,dim);
+   fes = new FiniteElementSpace(mesh,fec,dim,Ordering::byVDIM);
+   ndofs = fes->GetTrueVSize();
+   mesh->SetNodalFESpace(fes);
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh->bdr_attributes.Max());
+   }
+   ess_bdr = 0; ess_bdr[1] = 1;
+   fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+   // Solution GridFunction
+   x.SetSpace(fes);  x = 0.0;
+   // RHS
+   b.Update(fes); 
+   // Elasticity operator
+   lambda.SetSize(mesh->attributes.Max()); lambda = 57.6923076923;
+   mu.SetSize(mesh->attributes.Max()); mu = 38.4615384615;
+
+   lambda_cf.UpdateConstants(lambda);
+   mu_cf.UpdateConstants(mu);
+   a = new BilinearForm(fes);
+   a->SetDiagonalPolicy(mfem::Operator::DIAG_ONE);
+   a->AddDomainIntegrator(new ElasticityIntegrator(lambda_cf,mu_cf));
+}
+
+void ElasticityProblem::FormLinearSystem()
+{
+   if (!formsystem) 
+   {
+      formsystem = true;
+      b.Assemble();
+      a->Assemble();
+      a->FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+   }
+}
+void ElasticityProblem::UpdateLinearSystem()
+{
+   if (formsystem)
+   {
+      b.Update();
+      a->Update();
+      formsystem = false;
+   }
+   FormLinearSystem();
+}
+
+ContactProblem::ContactProblem(ElasticityProblem * prob1_, ElasticityProblem * prob2_)
+: prob1(prob1_), prob2(prob2_)
+{
+   // 1. Set up block system
+   Mesh* mesh1 = prob1->GetMesh();
+   int dim = mesh1->Dimension();
+
+   nodes0.SetSpace(mesh1->GetNodes()->FESpace());
+   nodes0 = *mesh1->GetNodes();
+   nodes1 = mesh1->GetNodes();
+
+   Vector delta1(dim);
+   delta1 = 0.0; delta1[0] = 0.1;
+   prob1->SetDisplacementDirichletData(delta1);
+   prob1->FormLinearSystem();
+   
+   Vector delta2(dim);
+   delta2 = 0.0; 
+   prob2->SetDisplacementDirichletData(delta2);
+   prob2->FormLinearSystem();
+
+   int ndof1 = prob1->GetNumDofs();
+   int ndof2 = prob2->GetNumDofs();
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = ndof1;
+   offsets[2] = ndof2;
+   offsets.PartialSum();
+
+   BlockMatrix Kb(offsets);
+   SparseMatrix A1 = prob1->GetOperator();
+   SparseMatrix A2 = prob2->GetOperator();
+
+   Kb.SetBlock(0,0,&A1);
+   Kb.SetBlock(1,1,&A2);
+
+   K = Kb.CreateMonolithic();
+   K->Threshold(0.0);
+   K->SortColumnIndices();
+
+   B = new BlockVector(offsets);
+   B->GetBlock(0).Set(1.0, prob1->GetRHS());
+   B->GetBlock(1).Set(1.0, prob2->GetRHS());
+
+   ComputeContactVertrices();
+}
+
+void ContactProblem::ComputeContactVertrices()
+{
+   if (npoints>0) return;
+   Mesh * mesh2 = prob2->GetMesh();
+   Array<int> vert;
+   for (int b=0; b<mesh2->GetNBE(); b++)
+   {
+      if (mesh2->GetBdrAttribute(b) == 3)
+      {
+         mesh2->GetBdrElementVertices(b, vert);
+         for (auto v : vert)
+         {
+            contact_vertices.insert(v);
+         }
+      }
+   }
+   npoints = contact_vertices.size();
+}
+
+void ContactProblem::ComputeGapFunctionAndDerivatives(const Vector &displ1, 
+                                                      const Vector & displ2)
+{
+   ComputeContactVertrices();
+
+   Mesh * mesh1 = prob1->GetMesh();
+   int dim = mesh1->Dimension();
+   Mesh * mesh2 = prob2->GetMesh();
+
+   int ndof1 = prob1->GetNumDofs();
+   int ndof2 = prob2->GetNumDofs();
+   int ndofs = ndof1 + ndof2;
+
+   int nv1 = mesh1->GetNV();
+   // connectivity of the second mesh
+
+   Array<int> conn2(npoints); 
+   // mesh2->MoveNodes(displ2);
+   Vector xyz(dim * npoints);
+
+   int cnt = 0;
+   for (auto v : contact_vertices)
+   {
+      for (int d = 0; d<dim; d++)
+      {
+         xyz(cnt*dim + d) = mesh2->GetVertex(v)[d]+displ2[v*dim+d];
+      }
+      conn2[cnt] = v + nv1; 
+      cnt++;
+   }
+   
+   MFEM_VERIFY(cnt == npoints, "");
+   gapv.SetSize(npoints*dim);
+
+   // segment reference coordinates of the closest point
+   Vector xi1(npoints*(dim-1));
+   Array<int> conn1(npoints*4);
+   
+   // add(nodes0, displ1, *nodes1);
+   FindPointsInMesh(*mesh1, xyz, conn1, xi1);
+
+   DenseMatrix coordsm(npoints*4, dim);
+   for (int i=0; i<npoints; i++)
+   {
+      for (int j=0; j<4; j++)
+      {
+         for (int k=0; k<dim; k++)
+         {
+            coordsm(i*4+j,k) = mesh1->GetVertex(conn1[i*4+j])[k]+displ1[dim*conn1[i*4+j]+k];
+         }
+      }
+   }
+
+   if (M)
+   {
+      delete M;
+      for (int i = 0; i<dM.Size(); i++)
+      {
+         delete dM[i];
+      }
+      dM.SetSize(0);
+   }
+
+   int h = npoints;
+   M = new SparseMatrix(h,ndofs);
+   dM.SetSize(npoints);
+   for (int i = 0; i<npoints; i++)
+   {
+      dM[i] = new SparseMatrix(ndofs,ndofs);
+   }
+   Assemble_Contact(xyz, xi1, coordsm, conn2, conn1, gapv, *M, dM);
+}
+
+
+double ContactProblem::E(const Vector & d)
+{
+   return 0.5 * K->InnerProduct(d, d) - InnerProduct(d, *B);
+}
+
+void ContactProblem::DdE(const Vector &d, Vector &gradE)
+{
+   gradE.SetSize(K->Height());
+   K->Mult(d, gradE);
+   gradE.Add(-1.0, *B); 
+}
+
+SparseMatrix* ContactProblem::DddE(const Vector &d)
+{
+  return K; 
+}
+
+void ContactProblem::g(const Vector &d, Vector &gd)
+{
+   int ndof1 = prob1->GetNumDofs();
+   int ndof2 = prob2->GetNumDofs();
+   double * data = d.GetData();
+   Vector displ1(data,ndof1);
+   Vector displ2(&data[ndof1],ndof2);
+   if (recompute)
+   {
+      ComputeGapFunctionAndDerivatives(displ1, displ2);
+      recompute = false;
+   }
+
+   gd = GetGapFunction();
+}
+
+SparseMatrix* ContactProblem::Ddg(const Vector &d)
+{
+  return GetJacobian();
+}
+
+SparseMatrix* ContactProblem::lDddg(const Vector &d, const Vector &l)
+{
+   return nullptr; // for now
+}
+
+QPContactProblem::QPContactProblem(ElasticityProblem * prob1_, ElasticityProblem * prob2_)
+: ContactProblem(prob1_,prob2_) 
+{ 
+   ContactProblem::ComputeContactVertrices();
+   dimS = npoints;
+   dimD = K->Height();
+}
+
+// E(d) = 1 / 2 d^T K d + f^T d
+double QPContactProblem::E(const Vector &d)
+{
+   return ContactProblem::E(d);
+}
+
+// gradient(E) = K d + f
+void QPContactProblem::DdE(const Vector &d, Vector &gradE)
+{
+   ContactProblem::DdE(d,gradE);
+}
+
+// Hessian(E) = K
+SparseMatrix* QPContactProblem::DddE(const Vector &d)
+{
+  return ContactProblem::DddE(d);
+}
+
+// g(d) = J * d + g0 >= 0
+void QPContactProblem::g(const Vector &d, Vector &gd)
+{
+  Vector g0;
+  ContactProblem::g(d,g0);
+  M->Mult(d, gd);
+  gd.Add(1.0, g0);
+}
+
+// Jacobian(g) = J
+SparseMatrix* QPContactProblem::Ddg(const Vector &d)
+{
+  return M;
+}
+
+SparseMatrix* QPContactProblem::lDddg(const Vector &d, const Vector &l)
+{
+   return ContactProblem::lDddg(d,l);
+}
+
+
+QPOptContactProblem::QPOptContactProblem(ContactProblem * problem_)
+: problem(problem_)
+{
+   dimU = problem->GetNumDofs();
+   dimM = problem->GetNumConstraints();
+   dimC = problem->GetNumConstraints();
+   block_offsets.SetSize(3);
+   block_offsets[0] = 0;
+   block_offsets[1] = dimU;
+   block_offsets[2] = dimM;
+   block_offsets.PartialSum();
+   ml.SetSize(dimM); ml = 0.0;
+   Vector negone(dimM); negone = -1.0;
+   NegId = new SparseMatrix(negone);
+}
+
+int QPOptContactProblem::GetDimU() { return dimU; }
+
+int QPOptContactProblem::GetDimM() { return dimM; }
+
+int QPOptContactProblem::GetDimC() { return dimC; }
+
+Vector & QPOptContactProblem::Getml() { return ml; }
+
+SparseMatrix * QPOptContactProblem::Duuf(const BlockVector & x)
+{
+   return problem->DddE(x.GetBlock(0));
+}
+
+SparseMatrix * QPOptContactProblem::Dumf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+SparseMatrix * QPOptContactProblem::Dmuf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+SparseMatrix * QPOptContactProblem::Dmmf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+SparseMatrix * QPOptContactProblem::Duc(const BlockVector & x)
+{
+   return problem->Ddg(x.GetBlock(0));
+}
+
+SparseMatrix * QPOptContactProblem::Dmc(const BlockVector & x)
+{
+   return NegId;
+}
+
+SparseMatrix * QPOptContactProblem::lDuuc(const BlockVector & x, const Vector & l)
+{
+   return nullptr;
+}
+
+void QPOptContactProblem::c(const BlockVector &x, Vector & y)
+{
+   Vector g0;
+   problem->g(x.GetBlock(0),g0); // gap function
+   g0.Add(-1.0, x.GetBlock(1));  
+
+   problem->GetJacobian()->Mult(x.GetBlock(0),y);
+   y.Add(1.0, g0);
+}
+
+double QPOptContactProblem::CalcObjective(const BlockVector & x)
+{
+   return problem->E(x.GetBlock(0));
+}
+
+void QPOptContactProblem::CalcObjectiveGrad(const BlockVector & x, BlockVector & y)
+{
+   problem->DdE(x.GetBlock(0), y.GetBlock(0));
+   y.GetBlock(1) = 0.0;
+}
+
+QPOptContactProblem::~QPOptContactProblem()
+{
+   delete NegId;
+}

miniapps/tribol/contact/problems/problems.hpp

@@ -0,0 +1,169 @@
+#include "problems_util.hpp"
+
+
+class ElasticityProblem
+{
+private:
+   bool formsystem = false;
+   Mesh * mesh = nullptr;
+   int order;
+   int ndofs;
+   FiniteElementCollection * fec = nullptr;
+   FiniteElementSpace * fes = nullptr;
+   Vector lambda, mu;
+   PWConstCoefficient lambda_cf, mu_cf;
+   Array<int> ess_bdr, ess_tdof_list;
+   BilinearForm *a=nullptr;
+   LinearForm b;
+   GridFunction x;
+   SparseMatrix A;
+   Vector B,X;
+   void Init();
+public:
+   ElasticityProblem(const char *mesh_file , int ref, int order_ = 1) : order(order_) 
+   {
+      mesh = new Mesh(mesh_file,1,1);
+      for (int i = 0; i<ref; i++)
+      {
+         mesh->UniformRefinement();
+      }
+      Init();
+   }
+
+   Mesh * GetMesh() { return mesh; }
+   FiniteElementSpace * GetFESpace() { return fes; }
+   int GetNumDofs() { return ndofs; }
+   SparseMatrix & 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 FormLinearSystem();
+   void UpdateLinearSystem();
+
+   void SetDisplacementDirichletData(const Vector & delta) 
+   {
+      VectorConstantCoefficient delta_cf(delta);
+      x.ProjectBdrCoefficient(delta_cf,ess_bdr);
+   };
+
+   void UpdateDisplacement(const Vector & x_) 
+   {
+      // x = x_;
+      // mesh->MoveVertices(x);
+      // mesh->NodesUpdated();
+   };
+
+   GridFunction & GetDisplacementGridFunction() {return x;};
+   Array<int> & GetEssentialDofs() {return ess_tdof_list;};
+
+   ~ElasticityProblem()
+   {
+      delete a;
+      delete fes;
+      delete fec;
+      delete mesh;
+   }
+};
+
+
+class ContactProblem
+{
+private:
+   ElasticityProblem * prob1 = nullptr;
+   ElasticityProblem * prob2 = nullptr;
+   GridFunction nodes0;
+   GridFunction *nodes1 = nullptr;
+   std::set<int> contact_vertices;
+   bool recompute = true;
+
+protected:
+   int npoints=0;
+   SparseMatrix *K =nullptr;
+   BlockVector *B = nullptr;
+   Vector gapv;
+   Array<SparseMatrix*> dM;
+   SparseMatrix * M=nullptr;
+   void ComputeContactVertrices();
+public:
+   ContactProblem(ElasticityProblem * prob1_, ElasticityProblem * prob2_);
+
+   ElasticityProblem * GetElasticityProblem1() {return prob1;}
+   ElasticityProblem * GetElasticityProblem2() {return prob2;}
+
+   int GetNumDofs() {return K->Height();}
+   int GetNumConstraints() {return npoints;}
+   Vector & GetGapFunction() {return gapv;}
+   SparseMatrix * GetJacobian() {return M;}
+   Array<SparseMatrix*> & GetHessian() {return dM;}
+   void ComputeGapFunctionAndDerivatives(const Vector & displ1, const Vector &displ2);
+
+   virtual double E(const Vector & d);
+   virtual void DdE(const Vector &d, Vector &gradE);
+   virtual SparseMatrix* DddE(const Vector &d);
+   void g(const Vector &d, Vector &gd);
+   virtual SparseMatrix* Ddg(const Vector &d);
+   virtual SparseMatrix* lDddg(const Vector &d, const Vector &l);
+
+   ~ContactProblem()
+   {
+      delete B;
+      delete K;
+      delete M;
+      for (int i = 0; i<dM.Size(); i++)
+      {
+         delete dM[i];
+      }
+   }
+};
+
+
+class QPContactProblem : public ContactProblem
+{
+private:
+   int dimD, dimS;
+public:
+   QPContactProblem(ElasticityProblem * prob1_, ElasticityProblem * prob2_);
+
+   double E(const Vector & d);
+   void DdE(const Vector &d, Vector &gradE);
+   SparseMatrix* DddE(const Vector &d);
+   void g(const Vector &d, Vector &gd);
+   SparseMatrix* Ddg(const Vector &d);
+   SparseMatrix* lDddg(const Vector &d, const Vector &l);
+};
+
+
+class QPOptContactProblem
+{
+private:
+   ContactProblem * problem = nullptr;
+   int dimU, dimM, dimC;
+   Array<int> block_offsets;
+   Vector ml;
+   SparseMatrix * NegId = nullptr;
+public:
+   QPOptContactProblem(ContactProblem * problem_);
+   int GetDimU();
+   int GetDimM();
+   int GetDimC();
+   Vector & Getml();
+   SparseMatrix * Duuf(const BlockVector &);
+   SparseMatrix * Dumf(const BlockVector &);
+   SparseMatrix * Dmuf(const BlockVector &);
+   SparseMatrix * Dmmf(const BlockVector &);
+   SparseMatrix * Duc(const BlockVector &);
+   SparseMatrix * Dmc(const BlockVector &);
+   SparseMatrix * lDuuc(const BlockVector &, const Vector &);
+   void c(const BlockVector &, Vector &);
+   double CalcObjective(const BlockVector &);
+   void CalcObjectiveGrad(const BlockVector &, BlockVector &);
+   ~QPOptContactProblem();
+};

miniapps/tribol/contact/problems/problems_util.hpp

@@ -0,0 +1,65 @@
+
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+void BasisEval(const Vector xi, Vector &N, DenseMatrix &dNdxi); // dNdxi is 2*4
+void BasisEvalDerivs(const Vector xi, Vector& N, DenseMatrix& dNdxi,
+                     DenseMatrix& dN2dxi);
+// returns the vector and matrix form of the shape functions and its derivative
+void BasisVectorDerivs(const Vector xi, DenseMatrix& N, DenseMatrix& dNdxi,
+                       DenseMatrix& ddNdxi);
+void cross(const Vector a, const Vector b, Vector& c);
+// a outer b
+void outer(const Vector a, const Vector b, DenseMatrix& c);
+// dphidxi 2*4
+// coords 4*3
+void ComputeNormal(const DenseMatrix& dphidxi, const DenseMatrix& coords,
+                   Vector& normal, double& nnorm);
+void SlaveToMaster(const DenseMatrix& m_coords, const Vector& s_x, Vector& xi);
+
+// m_coords is expected to be 4 * 3
+void  ComputeGapJacobian(const Vector x_s, const Vector xi,
+                         const DenseMatrix m_coords,
+                         double& gap, Vector& normal, Vector& dgdxm, Vector& dgdxs);
+
+void ComputeGapHessian(const Vector x_s, const Vector xi,
+                       const DenseMatrix m_coords,
+                       DenseMatrix& dg2dx);
+void NodeSegConPairs(const Vector x1, const Vector xi2,
+                     const DenseMatrix coords2,
+                     double& node_g, Vector& node_dg, DenseMatrix& node_dg2);
+// coordsm : (npoints*4, 3) use what class?
+// m_conn: (npoints*4)
+void Assemble_Contact(const Vector x_s,
+                      const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector& g, SparseMatrix& M,
+                      Array<SparseMatrix *> & dM);
+
+void Assemble_Contact(const Vector x_s, const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector& g, SparseMatrix& M,
+                      const Array<int> & points_map);
+
+void Assemble_Contact(const Vector x_s,
+                      const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector & g, SparseMatrix & M1, SparseMatrix & M2, 
+                      Array<SparseMatrix *> & dM11,
+                      Array<SparseMatrix *> & dM12,
+                      Array<SparseMatrix *> & dM21,
+                      Array<SparseMatrix *> & dM22);
+void Assemble_Contact(const Vector x_s,
+                      const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector & g, SparseMatrix & M1, SparseMatrix & M2,const Array<int> & points_map); 
+
+void FindSurfaceToProject(Mesh& mesh, const int elem, int& cbdrface);
+
+Vector GetNormalVector(Mesh & mesh, const int elem, const double *ref,
+                       int & refFace, int & refNormal, bool & interior);
+int GetHexVertex(int cdim, int c, int fa, int fb, Vector & refCrd);
+
+// Coordinates in xyz are assumed to be ordered as [X, Y, Z]
+// where X is the list of x-coordinates for all points and so on.
+// conn: connectivity of the target surface elements
+// xi: surface reference cooridnates for the cloest point, involves a linear transformation from [0,1] to [-1,1]
+void FindPointsInMesh(Mesh & mesh, Vector const& xyz, Array<int>& conn, Vector& xi);

miniapps/tribol/contact/util/mpicomm.cpp

@@ -0,0 +1,530 @@
+#include "mpicomm.hpp"
+#include "util.hpp"
+
+
+MPICommunicator::MPICommunicator(MPI_Comm comm_, int offset_, int gsize)
+: comm(comm_), offset(offset_)
+{
+   MPI_Comm_size(comm,&num_procs);
+   MPI_Comm_rank(comm,&myid);
+   offsets.resize(num_procs);
+   MPI_Allgather(&offset,1,MPI_INT,&offsets[0],1,MPI_INT,comm);
+   lsize = (myid == num_procs-1) ? gsize - offsets[myid]
+                                : offsets[myid+1]-offsets[myid];
+
+   send_count.SetSize(num_procs); send_count = 0;
+   send_displ.SetSize(num_procs); send_displ = 0;
+   recv_count.SetSize(num_procs); recv_count = 0;
+   recv_displ.SetSize(num_procs); recv_displ = 0;
+} 
+
+MPICommunicator::MPICommunicator(MPI_Comm comm_, Array<unsigned int> & destination_procs_)
+: comm(comm_), destination_procs(destination_procs_)
+{
+   MPI_Comm_size(comm,&num_procs);
+   MPI_Comm_rank(comm,&myid);   
+   send_count.SetSize(num_procs); 
+   send_displ.SetSize(num_procs); 
+   recv_count.SetSize(num_procs); 
+   recv_displ.SetSize(num_procs); 
+   resetcounts();
+}
+
+
+int MPICommunicator::get_rank(int dof)
+{
+   if (num_procs == 1) { return 0; }
+   std::vector<int>::iterator up;
+   up=std::upper_bound(offsets.begin(), offsets.end(),dof);
+   return std::distance(offsets.begin(),up)-1;
+}
+
+
+void MPICommunicator::Communicate(const Vector & x_s, Vector & x_r, int vdim, int ordering)
+{
+   int npts = x_s.Size()/vdim;
+   MFEM_VERIFY(npts == destination_procs.Size(), "Inconsistent number of points to be send");
+
+   // construct send count
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      send_count[rank] += vdim + 1; // including the sending processor id
+   }
+
+   // 2. Compute recv_count
+   MPI_Alltoall(&send_count[0],1,MPI_INT,&recv_count[0],1,MPI_INT,comm);
+
+   // 3. Compute displacements
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   int sbuff_size = send_count.Sum();
+   int rbuff_size = recv_count.Sum();
+
+   // 4. Allocate memory and fill in send buffers
+   Array<double> sendvals(sbuff_size);  sendvals = 0.0;
+   Array<int> sendoffs(num_procs); sendoffs = 0;
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      int j = send_displ[rank] + sendoffs[rank];
+      sendoffs[rank] += vdim+1;
+      sendvals[j] = (double)myid;
+      for (int k = 0; k<vdim; k++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? k*npts+i : i*vdim + k;
+         sendvals[j+k+1] = x_s(kk);
+      }
+   }
+
+   // 5. Communication
+   Array<double> recvvals(rbuff_size);
+
+   double * sendvals_ptr = nullptr;
+   double * recvvals_ptr = nullptr;
+   if (sbuff_size !=0 ) { sendvals_ptr = &sendvals[0]; }   
+   if (rbuff_size !=0 ) { recvvals_ptr = &recvvals[0]; }
+
+   MPI_Alltoallv(sendvals_ptr, send_count, send_displ, MPI_DOUBLE, recvvals_ptr,
+                 recv_count, recv_displ, MPI_DOUBLE, comm);
+
+   // 6. Unpack
+   int n = rbuff_size/(vdim+1);
+   origin_procs.SetSize(n);
+   x_r.SetSize(vdim*n);
+   for (int i = 0; i<n; i++)
+   {
+      origin_procs[i] = (unsigned int)recvvals[(vdim+1)*i];
+      for (int j=0; j<vdim; j++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? j*n+i : i*vdim + j;
+         x_r(kk) = recvvals[(vdim+1)*i + j+1];
+      }
+   }
+   resetcounts();
+}
+
+void MPICommunicator::Communicate(const Array<unsigned int> & x_s, Array<unsigned int> & x_r, int vdim, int ordering)
+{
+   int npts = x_s.Size()/vdim;
+   MFEM_VERIFY(npts == destination_procs.Size(), "Inconsistent number of points to be send");
+
+   // construct send count
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      send_count[rank] += vdim + 1; // including the sending processor id
+   }
+
+   // 2. Compute recv_count
+   MPI_Alltoall(&send_count[0],1,MPI_INT,&recv_count[0],1,MPI_INT,comm);
+
+   // 3. Compute displacements
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   int sbuff_size = send_count.Sum();
+   int rbuff_size = recv_count.Sum();
+
+   // 4. Allocate memory and fill in send buffers
+   Array<unsigned int> sendvals(sbuff_size);  sendvals = 0.0;
+   Array<int> sendoffs(num_procs); sendoffs = 0;
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      int j = send_displ[rank] + sendoffs[rank];
+      sendoffs[rank] += vdim+1;
+      sendvals[j] = myid;
+      for (int k = 0; k<vdim; k++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? k*npts+i : i*vdim + k;
+         sendvals[j+k+1] = x_s[kk];
+      }
+   }
+
+   // 5. Communication
+   Array<unsigned int> recvvals(rbuff_size);
+
+   unsigned int * sendvals_ptr = nullptr;
+   unsigned int * recvvals_ptr = nullptr;
+   if (sbuff_size !=0 ) { sendvals_ptr = &sendvals[0]; }   
+   if (rbuff_size !=0 ) { recvvals_ptr = &recvvals[0]; }
+
+   MPI_Alltoallv(sendvals_ptr, send_count, send_displ, MPI_UNSIGNED, recvvals_ptr,
+                 recv_count, recv_displ, MPI_UNSIGNED, comm);
+
+   // 6. Unpack
+   int n = rbuff_size/(vdim+1);
+   origin_procs.SetSize(n);
+   x_r.SetSize(vdim*n);
+   for (int i = 0; i<n; i++)
+   {
+      origin_procs[i] = recvvals[(vdim+1)*i];
+      for (int j=0; j<vdim; j++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? j*n+i : i*vdim + j;
+         x_r[kk] = recvvals[(vdim+1)*i + j+1];
+      }
+   }
+   resetcounts();
+}
+
+void MPICommunicator::Communicate(const Array<int> & x_s, Array<int> & x_r, int vdim, int ordering)
+{
+   int npts = x_s.Size()/vdim;
+   MFEM_VERIFY(npts == destination_procs.Size(), "Inconsistent number of points to be send");
+
+   // construct send count
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      send_count[rank] += vdim + 1; // including the sending processor id
+   }
+
+   // 2. Compute recv_count
+   MPI_Alltoall(&send_count[0],1,MPI_INT,&recv_count[0],1,MPI_INT,comm);
+
+   // 3. Compute displacements
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   int sbuff_size = send_count.Sum();
+   int rbuff_size = recv_count.Sum();
+
+   // 4. Allocate memory and fill in send buffers
+   Array<int> sendvals(sbuff_size);  sendvals = 0.0;
+   Array<int> sendoffs(num_procs); sendoffs = 0;
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      int j = send_displ[rank] + sendoffs[rank];
+      sendoffs[rank] += vdim+1;
+      sendvals[j] = myid;
+      for (int k = 0; k<vdim; k++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? k*npts+i : i*vdim + k;
+         sendvals[j+k+1] = x_s[kk];
+      }
+   }
+
+   // 5. Communication
+   Array<int> recvvals(rbuff_size);
+
+   int * sendvals_ptr = nullptr;
+   int * recvvals_ptr = nullptr;
+   if (sbuff_size !=0 ) { sendvals_ptr = &sendvals[0]; }   
+   if (rbuff_size !=0 ) { recvvals_ptr = &recvvals[0]; }
+
+   MPI_Alltoallv(sendvals_ptr, send_count, send_displ, MPI_INT, recvvals_ptr,
+                 recv_count, recv_displ, MPI_INT, comm);
+
+   // 6. Unpack
+   int n = rbuff_size/(vdim+1);
+   origin_procs.SetSize(n);
+   x_r.SetSize(vdim*n);
+   for (int i = 0; i<n; i++)
+   {
+      origin_procs[i] = (unsigned int)recvvals[(vdim+1)*i];
+      for (int j=0; j<vdim; j++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? j*n+i : i*vdim + j;
+         x_r[kk] = recvvals[(vdim+1)*i + j+1];
+      }
+   }
+   resetcounts();
+}
+
+void MPICommunicator::Communicate(const DenseMatrix & A_s, DenseMatrix & A_r, int vdim, int ordering)
+{
+   // matrix width corresponds to dim coordinates
+   // matrix rows might include vdim copies
+   int npts = A_s.Height()/vdim;
+   int dim = A_s.Width();
+   MFEM_VERIFY(npts == destination_procs.Size(), "Inconsistent number of points to be send");
+
+   // construct send count
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      send_count[rank] += dim*vdim + 1; // including the sending processor id
+   }
+
+   // 2. Compute recv_count
+   MPI_Alltoall(&send_count[0],1,MPI_INT,&recv_count[0],1,MPI_INT,comm);
+
+   // 3. Compute displacements
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   int sbuff_size = send_count.Sum();
+   int rbuff_size = recv_count.Sum();
+
+   // 4. Allocate memory and fill in send buffers
+   Array<double> sendvals(sbuff_size);  sendvals = 0.0;
+   Array<int> sendoffs(num_procs); sendoffs = 0;
+   for (int i = 0; i<npts; i++)
+   {
+      int rank = destination_procs[i];
+      int j = send_displ[rank] + sendoffs[rank];
+      sendoffs[rank] += dim*vdim+1;
+      sendvals[j] = myid;
+      for (int k = 0; k<vdim; k++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? k*npts+i : i*vdim + k;
+         for (int d=0; d<dim; d++)
+         {
+            sendvals[j+k*dim+d+1] = A_s(kk,d);
+         }
+      }
+   }
+   // 5. Communication
+   Array<double> recvvals(rbuff_size);
+
+   double * sendvals_ptr = nullptr;
+   double * recvvals_ptr = nullptr;
+   if (sbuff_size !=0 ) { sendvals_ptr = &sendvals[0]; }   
+   if (rbuff_size !=0 ) { recvvals_ptr = &recvvals[0]; }
+
+   MPI_Alltoallv(sendvals_ptr, send_count, send_displ, MPI_DOUBLE, recvvals_ptr,
+                 recv_count, recv_displ, MPI_DOUBLE, comm);
+
+   // 6. Unpack
+   int n = rbuff_size/(dim*vdim+1);
+   origin_procs.SetSize(n);
+   A_r.SetSize(vdim*n,dim);
+
+   for (int i = 0; i<n; i++)
+   {
+      origin_procs[i] = (unsigned int)recvvals[(dim*vdim+1)*i];
+      for (int j=0; j<vdim; j++)
+      {
+         int kk = (ordering == mfem::Ordering::byNODES) ? j*n+i : i*vdim + j;
+         for (int d=0; d<dim; d++)
+         {
+            A_r(kk,d) = recvvals[(dim*vdim+1)*i + j*dim + d+1];
+         }
+      }
+   }
+   resetcounts();
+
+}
+
+
+void MPICommunicator::Communicate(const SparseMatrix & mat_s , SparseMatrix & mat_r)
+{
+   // 1. Compute send_count 
+   int n = mat_s.NumRows();
+   for (int i = 0; i<n; i++)
+   {
+      int rsize = mat_s.RowSize(i);
+      if (rsize == 0) continue;
+      int rank = get_rank(i);
+      send_count[rank] += rsize+2;
+   }
+   // 2. Compute recv_count
+   MPI_Alltoall(&send_count[0],1,MPI_INT,&recv_count[0],1,MPI_INT,comm);
+
+   // 3. Compute displacements
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   int sbuff_size = send_count.Sum();
+   int rbuff_size = recv_count.Sum();
+
+   // 4. Allocate memory and fill in send buffers
+   Array<double> sendvals(sbuff_size);  sendvals = 0.0;
+   Array<int> sendcols(sbuff_size);  sendcols = 0;
+   Array<int> sendoffs(num_procs); sendoffs = 0;
+   Array<int> cols;
+   Vector vals;
+   for (int i = 0; i<n; i++)
+   {
+      int rsize = mat_s.RowSize(i);
+      if (rsize == 0) continue;
+      int rank = get_rank(i);
+      int j = send_displ[rank] + sendoffs[rank];
+      mat_s.GetRow(i,cols,vals);
+      sendoffs[rank] += rsize+2;
+      sendvals[j] = (double)i;
+      sendvals[j+1] = (double)rsize;
+      sendcols[j] = i;
+      sendcols[j+1] = rsize;
+      for (int l=0; l<rsize ; l++)
+      {
+         sendvals[j+l+2] = vals[l];
+         sendcols[j+l+2] = cols[l];
+      }
+   }
+
+   // 5. Communication
+   Array<double> recvvals(rbuff_size);
+   Array<int> recvcols(rbuff_size);
+
+   double * sendvals_ptr = nullptr;
+   double * recvvals_ptr = nullptr;
+   int * sendcols_ptr = nullptr;
+   int * recvcols_ptr = nullptr;
+   if (sbuff_size !=0 ) 
+   {
+      sendvals_ptr = &sendvals[0]; 
+      sendcols_ptr = &sendcols[0]; 
+   }   
+   if (rbuff_size !=0 ) 
+   {
+      recvvals_ptr = &recvvals[0]; 
+      recvcols_ptr = &recvcols[0]; 
+   }
+
+   MPI_Alltoallv(sendvals_ptr, send_count, send_displ, MPI_DOUBLE, recvvals_ptr,
+                 recv_count, recv_displ, MPI_DOUBLE, comm);
+
+   MPI_Alltoallv(sendcols_ptr, send_count, send_displ, MPI_INT, recvcols_ptr,
+                 recv_count, recv_displ, MPI_INT, comm);
+
+   // 6. Unpack and store to the output SparseMatrix
+   MFEM_VERIFY(mat_r.Height() == lsize, "Inconsistent row size of output SparseMatrix");
+   MFEM_VERIFY(mat_r.Width() == mat_s.Width(), "Inconsistent column size of output SparseMatrix");
+
+   int counter = 0;
+   while (counter < rbuff_size)
+   {
+      int row = recvcols[counter] - offset;
+      int size = recvcols[counter+1];
+      vals.SetSize(size);
+      cols.SetSize(size);
+      for (int i = 0; i<size; i++)
+      {
+         vals[i] = recvvals[counter+2 + i];
+         cols[i] = recvcols[counter+2 + i];
+      }
+      mat_r.AddRow(row,cols,vals);
+      counter += size+2; 
+   }
+   MFEM_VERIFY(counter == rbuff_size, "inconsistent rbuff size");
+   mat_r.Finalize();
+   mat_r.SortColumnIndices();
+   resetcounts();
+}
+
+void MPICommunicator::Communicate(const Array<SparseMatrix*> & vmat_s, Array<SparseMatrix*> & vmat_r)
+{
+   // 1. Compute send_count 
+   for (int k = 0; k<vmat_s.Size(); k++)
+   {
+      if (!vmat_s[k]) continue;
+      if (vmat_s[k]->NumNonZeroElems() == 0) continue;
+      int nrows = vmat_s[k]->NumRows();
+      for (int i = 0; i<nrows; i++)
+      {
+         int rsize = vmat_s[k]->RowSize(i);
+         if (rsize == 0) continue;
+         int rank = get_rank(i);
+         send_count[rank] += rsize+3;
+      }
+   } 
+
+   // 2. Compute recv_count
+   MPI_Alltoall(&send_count[0],1,MPI_INT,&recv_count[0],1,MPI_INT,comm);
+   
+   // 3. Compute displacements
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   int sbuff_size = send_count.Sum();
+   int rbuff_size = recv_count.Sum();
+
+   // 4. Allocate memory and fill in send buffers
+   Array<double> sendvals(sbuff_size);  sendvals = 0.0;
+   Array<int> sendcols(sbuff_size);  sendcols = 0;
+   Array<int> sendoffs(num_procs); sendoffs = 0;
+   for (int k = 0; k<vmat_s.Size(); k++)
+   {
+      if (!vmat_s[k]) continue;
+      if (vmat_s[k]->NumNonZeroElems() == 0) continue;
+      int nrows = vmat_s[k]->NumRows();
+      for (int i = 0; i<nrows; i++)
+      {
+         int rsize = vmat_s[k]->RowSize(i);
+         if (rsize == 0) continue;
+         int rank = get_rank(i);
+         int j = send_displ[rank] + sendoffs[rank];
+         Array<int> cols;
+         Vector vals;
+         vmat_s[k]->GetRow(i,cols,vals);
+         sendoffs[rank] += rsize+3;
+         sendvals[j] = (double)k;
+         sendvals[j+1] = (double)i;
+         sendvals[j+2] = (double)rsize;
+         sendcols[j] = k;
+         sendcols[j+1] = i;
+         sendcols[j+2] = rsize;
+         for (int l=0; l<rsize ; l++)
+         {
+            sendvals[j+l+3] = vals[l];
+            sendcols[j+l+3] = cols[l];
+         }
+      }
+   }
+   // 5. Communication
+   Array<double> recvvals(rbuff_size);
+   Array<int> recvcols(rbuff_size);
+   double * sendvals_ptr = nullptr;
+   double * recvvals_ptr = nullptr;
+   int * sendcols_ptr = nullptr;
+   int * recvcols_ptr = nullptr;
+   if (sbuff_size !=0 ) 
+   {
+      sendvals_ptr = &sendvals[0]; 
+      sendcols_ptr = &sendcols[0]; 
+   }   
+   if (rbuff_size !=0 ) 
+   {
+      recvvals_ptr = &recvvals[0]; 
+      recvcols_ptr = &recvcols[0]; 
+   }
+
+   MPI_Alltoallv(sendvals_ptr, send_count, send_displ, MPI_DOUBLE, recvvals_ptr,
+                 recv_count, recv_displ, MPI_DOUBLE,comm);
+
+   MPI_Alltoallv(sendcols_ptr, send_count, send_displ, MPI_INT, recvcols_ptr,
+                 recv_count, recv_displ, MPI_INT,comm);
+
+   // 6. Unpack and store to the output SparseMatrix
+   int counter = 0;
+   while (counter < rbuff_size)
+   {
+      int npt = recvcols[counter];
+      int row = recvcols[counter+1] - offset;
+      int size = recvcols[counter+2];
+      Vector vals(size);
+      Array<int> cols(size);
+      for (int i = 0; i<size; i++)
+      {
+         vals[i] = recvvals[counter+3 + i];
+         cols[i] = recvcols[counter+3 + i];
+      }
+      vmat_r[npt]->AddRow(row,cols,vals);
+      counter += size+3; 
+   }
+   MFEM_VERIFY(counter == rbuff_size, "inconsistent size");
+
+   for (int i = 0; i<vmat_r.Size(); i++)
+   {
+      vmat_r[i]->Finalize();
+      vmat_r[i]->SortColumnIndices();
+   }
+   resetcounts();
+}

miniapps/tribol/contact/util/mpicomm.hpp

@@ -0,0 +1,47 @@
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+class MPICommunicator
+{
+private:
+    MPI_Comm comm;
+    int myid, num_procs;
+    Array<unsigned int > origin_procs;
+    Array<unsigned int > destination_procs;
+    int offset, lsize;
+    std::vector<int> offsets; 
+    Array<int> send_count;
+    Array<int> send_displ;
+    Array<int> recv_count;
+    Array<int> recv_displ;
+    void resetcounts()
+    {
+        send_count = 0;
+        send_displ = 0;
+        recv_count = 0;
+        recv_displ = 0;
+    }
+
+public:
+    MPICommunicator(MPI_Comm comm_, int offset_, int gsize);
+    MPICommunicator(MPI_Comm comm_, Array<unsigned int> & destination_procs_);
+
+    int get_rank(int dof);
+
+    Array<unsigned int> & GetOriginProcs() {return origin_procs;}
+    void UpdateDestinationProcs() 
+    {   
+        destination_procs.SetSize(origin_procs.Size());
+        destination_procs = origin_procs;
+        resetcounts();
+    }
+    void Communicate(const Vector & x_s, Vector & x_r, int vdim, int ordering);
+    void Communicate(const Array<int> & x_s, Array<int> & x_r, int vdim, int ordering);
+    void Communicate(const DenseMatrix & A_s, DenseMatrix & A_r, int vdim, int ordering);
+    void Communicate(const Array<unsigned int> & x_s, Array<unsigned int> & x_r, int vdim, int ordering);
+    void Communicate(const SparseMatrix & mat_s , SparseMatrix & mat_r);
+    void Communicate(const Array<SparseMatrix*> & vmat_s, Array<SparseMatrix*> & vmat_r);
+};

miniapps/tribol/contact/util/util.cpp

@@ -0,0 +1,171 @@
+#include "util.hpp"
+
+
+void PrintVertex(Mesh * mesh, int vertex)
+{
+   Array<int> vertices;
+   mfem::out << "vertex: " << vertex << ":  ";
+   double *coords = mesh->GetVertex(vertex);
+   mfem::out << "(" << coords[0] << ", " << coords[1] << ", " << coords[2] << ") \n";
+}
+
+void PrintElementVertices(Mesh * mesh, int elem)
+{
+   Array<int> vertices;
+   mfem::out <<  "elem: " << elem << ". Vertices = \n" ;
+   mesh->GetElementVertices(elem,vertices);
+   for (int i = 0; i<vertices.Size(); i++)
+   {
+      PrintVertex(mesh,vertices[i]);
+   }
+   mfem::out << endl;
+}
+
+void PrintFaceVertices(Mesh * mesh, int face)
+{
+   Array<int> vertices;
+   mfem::out << "face: " << face << ". Vertices = \n" ;
+   mesh->GetFaceVertices(face,vertices);
+   for (int i = 0; i<vertices.Size(); i++)
+   {
+      PrintVertex(mesh,vertices[i]);
+   }
+   mfem::out << endl;
+}
+
+void PrintSet(const std::set<int> & a, const char *aname)
+{
+   mfem::out << aname << " = " ;
+   for (std::set<int>::iterator it = a.begin(); it!= a.end(); it++)
+   {
+      mfem::out << *it << "  ";
+   }
+   mfem::out << endl;
+}
+
+void PrintVector(const Vector & a, const char *aname)
+{
+   int sz = a.Size();
+   mfem::out << aname << " = " ;
+   for (int i = 0; i<sz; i++)
+   {
+      mfem::out << a[i] << "  ";
+   }
+   mfem::out << endl;
+}
+
+void PrintVertex(Mesh * mesh, int vertex,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   if (myid == printid)
+   {
+      mfem::out << "myid = " << myid <<":   " <<  "vertex: " << vertex << ":  ";
+      double *coords = mesh->GetVertex(vertex);
+      mfem::out << "(" << coords[0] << ", " << coords[1] << ", " << coords[2] << ")\n";
+   }
+}
+
+void PrintElementVertices(Mesh * mesh, int elem,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   Array<int> vertices;
+   if (myid == printid)
+   {
+      mfem::out << "myid = " << myid <<":   " <<  "elem: " << elem <<
+                ". Vertices = \n" ;
+      mesh->GetElementVertices(elem,vertices);
+      for (int i = 0; i<vertices.Size(); i++)
+      {
+         PrintVertex(mesh,vertices[i],printid);
+      }
+      mfem::out << endl;
+   }
+}
+
+void PrintFaceVertices(Mesh * mesh, int face,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   Array<int> vertices;
+   if (myid == printid)
+   {
+      mfem::out << "myid = " << myid <<":   " <<  "face: " << face <<
+                ". Vertices = \n" ;
+      mesh->GetFaceVertices(face,vertices);
+      for (int i = 0; i<vertices.Size(); i++)
+      {
+         PrintVertex(mesh,vertices[i],printid);
+      }
+      mfem::out << endl;
+   }
+}
+
+
+void PrintSet(const std::set<int> & a, const char *aname,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   if (myid == printid)
+   {
+      mfem::out << "myid = " << myid <<":   " << aname << " = " ;
+      for (std::set<int>::iterator it = a.begin(); it!= a.end(); it++)
+      {
+         mfem::out << *it << "  ";
+      }
+      mfem::out << endl;
+   }
+}
+
+void PrintVector(const Vector & a, const char *aname,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   if (myid == printid)
+   {
+      int sz = a.Size();
+      mfem::out << "myid = " << myid <<":   " << aname << " = " ;
+      for (int i = 0; i<sz; i++)
+      {
+         mfem::out << a[i] << "  ";
+      }
+      mfem::out << endl;
+   }
+}
+
+void PrintVector(const std::vector<int> & a, const char *aname,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   if (myid == printid)
+   {
+      int sz = a.size();
+      mfem::out << "myid = " << myid <<":   " << aname << " = " ;
+      for (int i = 0; i<sz; i++)
+      {
+         mfem::out << a[i] << "  ";
+      }
+      mfem::out << endl;
+   }
+}
+
+void PrintVector(const std::vector<unsigned int> & a, const char *aname,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   if (myid == printid)
+   {
+      int sz = a.size();
+      mfem::out << "myid = " << myid <<":   " << aname << " = " ;
+      for (int i = 0; i<sz; i++)
+      {
+         mfem::out << a[i] << "  ";
+      }
+      mfem::out << endl;
+   }
+}
+
+void PrintSparseMatrix(const SparseMatrix & a, const char *aname,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   if (myid == printid)
+   {
+      mfem::out << "myid = " << myid <<":   " << aname << " = " ;
+      a.PrintMatlab(mfem::out);
+   }
+   mfem::out << endl;
+}

miniapps/tribol/contact/util/util.hpp

@@ -0,0 +1,46 @@
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+void PrintVertex(Mesh * mesh, int vertex);
+void PrintElementVertices(Mesh * mesh, int elem);
+void PrintFaceVertices(Mesh * mesh, int face);
+template <class T>
+void PrintArray(const Array<T> & a, const char *aname)
+{
+   int sz = a.Size();
+   mfem::out << aname << " = " ;
+   for (int i = 0; i<sz; i++)
+   {
+      mfem::out << a[i] << "  ";
+   }
+   mfem::out << endl;
+}
+void PrintSet(const std::set<int> & a, const char *aname);
+void PrintVector(const Vector & a, const char *aname);
+
+// for parallel 
+void PrintVertex(Mesh * mesh, int vertex,  int printid);
+void PrintElementVertices(Mesh * mesh, int elem,  int printid);
+void PrintFaceVertices(Mesh * mesh, int face,  int printid);
+template <class T>
+void PrintArray(const Array<T> & a, const char *aname,  int printid)
+{
+   int myid = Mpi::WorldRank();
+   if (myid == printid)
+   {
+      int sz = a.Size();
+      mfem::out << "myid = " << myid <<":   " << aname << " = " ;
+      for (int i = 0; i<sz; i++)
+      {
+         mfem::out << a[i] << "  ";
+      }
+      mfem::out << endl;
+   }
+}
+void PrintSet(const std::set<int> & a, const char *aname,  int printid);
+void PrintVector(const Vector & a, const char *aname,  int printid);
+void PrintVector(const std::vector<int> & a, const char *aname,  int printid);
+void PrintVector(const std::vector<unsigned int> & a, const char *aname,  int printid);
+void PrintSparseMatrix(const SparseMatrix & a, const char *aname,  int printid);

miniapps/tribol/makefile

@@ -67,5 +67,7 @@ clean: clean-build clean-exec
 clean-build:
 	rm -f *.o *~ $(MINIAPPS)
 	rm -rf *.dSYM *.TVD.*breakpoints
+	rm -rf ContactPatchTestSurface_000000*
+	rm -rf ContactPatchTestVolume_000000*
 
 clean-exec: