adding no contact AMG solves for every step

linalg/blockoperator.cpp

@@ -357,4 +357,154 @@ BlockLowerTriangularPreconditioner::~BlockLowerTriangularPreconditioner()
    }
 }
 
+
+BlockTriangularSymmetricPreconditioner::BlockTriangularSymmetricPreconditioner(
+   const Array<int> & offsets_)
+   : Solver(offsets_.Last()),
+     owns_blocks(0),
+     nBlocks(offsets_.Size() - 1),
+     offsets(0),
+     ops(nBlocks, nBlocks)
+{
+   ops = static_cast<Operator *>(NULL);
+   offsets.MakeRef(offsets_);
+}
+
+void BlockTriangularSymmetricPreconditioner::SetDiagonalBlock(int iblock,
+                                                              Operator *op)
+{
+   MFEM_VERIFY(offsets[iblock+1] - offsets[iblock] == op->Height() &&
+               offsets[iblock+1] - offsets[iblock] == op->Width(),
+               "incompatible Operator dimensions");
+
+   SetBlock(iblock, iblock, op);
+}
+
+void BlockTriangularSymmetricPreconditioner::SetBlock(int iRow, int iCol,
+                                                      Operator *op)
+{
+   MFEM_VERIFY(offsets[iRow+1] - offsets[iRow] == op->NumRows() &&
+               offsets[iCol+1] - offsets[iCol] == op->NumCols(),
+               "incompatible Operator dimensions");
+
+   ops(iRow, iCol) = op;
+}
+
+// Operator application
+
+void BlockTriangularSymmetricPreconditioner::ForwardPass(const Vector & x,
+                                                         Vector & y) const
+{
+   // Forward sweep: Solve for y1, then y2
+   for (int iRow = 0; iRow < nBlocks; ++iRow)
+   {
+      tmp.SetSize(offsets[iRow + 1] - offsets[iRow]);
+      tmp2.SetSize(offsets[iRow + 1] - offsets[iRow]);
+      tmp2 = 0.0;
+      tmp2 += xblock.GetBlock(iRow);  // tmp2 = xblock(iRow)
+
+      // Process the lower triangular part (jCol < iRow)
+      for (int jCol = 0; jCol < iRow; ++jCol)
+      {
+         if (ops(iRow, jCol))
+         {
+            ops(iRow, jCol)->Mult(yblock.GetBlock(jCol),
+                                  tmp);  // tmp = A(iRow,jCol) * yblock(jCol)
+            tmp2 -= tmp;  // tmp2 -= A(iRow, jCol) * yblock(jCol)
+         }
+      }
+
+      // Apply the diagonal block
+      if (ops(iRow, iRow))
+      {
+         ops(iRow, iRow)->Mult(tmp2,
+                               yblock.GetBlock(iRow));  // yblock(iRow) = A(iRow,iRow)^-1 * tmp2
+      }
+      else
+      {
+         yblock.GetBlock(iRow) = tmp2;  // If no diagonal operator, set yblock directly
+      }
+   }
+}
+
+void BlockTriangularSymmetricPreconditioner::BackwardPass(const Vector & x,
+                                                          Vector & y) const
+{
+   // Backward sweep: Adjust y1 based on y2
+   for (int iRow = nBlocks - 1; iRow >= 0; --iRow)
+   {
+      tmp.SetSize(offsets[iRow + 1] - offsets[iRow]);
+      tmp2.SetSize(offsets[iRow + 1] - offsets[iRow]);
+      tmp2 = 0.0;
+      tmp2 += xblock.GetBlock(iRow);  // tmp2 = yblock(iRow) from forward sweep
+
+      // Process the upper triangular part (jCol > iRow)
+      for (int jCol = iRow + 1; jCol < nBlocks; ++jCol)
+      {
+         if (ops(iRow, jCol))
+         {
+            ops(iRow, jCol)->Mult(yblock.GetBlock(jCol),
+                                  tmp);  // tmp = A(iRow,jCol) * yblock(jCol)
+            tmp2 -= tmp;  // tmp2 -= A(iRow,jCol) * yblock(jCol)
+         }
+      }
+
+      // Reapply diagonal block to correct y1
+      if (ops(iRow, iRow))
+      {
+         ops(iRow, iRow)->Mult(tmp2,
+                               yblock.GetBlock(iRow));  // Final correction for yblock(iRow)
+      }
+      else
+      {
+         yblock.GetBlock(iRow) = tmp2;  // If no diagonal operator, set yblock directly
+      }
+   }
+}
+
+
+
+void BlockTriangularSymmetricPreconditioner::Mult (const Vector & x,
+                                                   Vector & y) const
+{
+   MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
+   MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
+
+   // Update block views of the vectors y and x using offsets
+   yblock.Update(y.GetData(), offsets);
+   xblock.Update(x.GetData(), offsets);
+
+   // Initialize y to zero
+   y = 0.0;
+   ForwardPass(x,y);
+
+   // Update Residual
+   r.SetSize(x.Size());
+   r = 0.0; r+=x;
+   Op->AddMult(y,r,-1.0);
+
+   Vector y1(y);
+   yblock.Update(y1.GetData(), offsets);
+   xblock.Update(r.GetData(), offsets);
+   BackwardPass(r,y1);
+   y+=y1;
+}
+
+BlockTriangularSymmetricPreconditioner::~BlockTriangularSymmetricPreconditioner()
+{
+   if (owns_blocks)
+   {
+      for (int iRow=0; iRow < nBlocks; ++iRow)
+      {
+         for (int jCol=0; jCol < nBlocks; ++jCol)
+         {
+            delete ops(jCol,iRow);
+         }
+      }
+   }
+}
+
+
+
+
 }

linalg/blockoperator.hpp

@@ -288,6 +288,78 @@ private:
    mutable Vector tmp2;
 };
 
+
+class BlockTriangularSymmetricPreconditioner : public Solver
+{
+private:
+   const Operator * Op;
+public:
+   //! Constructor for BlockTriangularSymmetricPreconditioners with the same
+   //! block-structure for rows and columns.
+   /**
+    *  @param offsets  Offsets that mark the start of each row/column block
+    *                  (size nBlocks+1).
+    *
+    *  @note BlockTriangularSymmetricPreconditioner will not own/copy the data
+    *  contained in @a offsets.
+    */
+   BlockTriangularSymmetricPreconditioner(const Array<int> & offsets);
+
+   //! Add block op in the block-entry (iblock, iblock).
+   /**
+    * @param iblock  The block will be inserted in location (iblock, iblock).
+    * @param op      The Operator to be inserted.
+    */
+   void SetDiagonalBlock(int iblock, Operator *op);
+   //! Add a block opt in the block-entry (iblock, jblock).
+   /**
+    * @param iRow, iCol  The block will be inserted in location (iRow, iCol).
+    * @param op          The Operator to be inserted.
+    */
+   void SetBlock(int iRow, int iCol, Operator *op);
+   //! This method is present since required by the abstract base class Solver
+   virtual void SetOperator(const Operator &op) {Op = &op;}
+
+   //! Return the number of blocks
+   int NumBlocks() const { return nBlocks; }
+
+   //! Return a reference to block i,j.
+   Operator & GetBlock(int iblock, int jblock)
+   { MFEM_VERIFY(ops(iblock,jblock), ""); return *ops(iblock,jblock); }
+
+   //! Return the offsets for block starts
+   Array<int> & Offsets() { return offsets; }
+
+   /// Operator application
+   virtual void Mult (const Vector & x, Vector & y) const;
+
+   ~BlockTriangularSymmetricPreconditioner();
+
+   //! Controls the ownership of the blocks: if nonzero,
+   //! BlockTriangularSymmetricPreconditioner will delete all blocks that are set
+   //! (non-NULL); the default value is zero.
+   int owns_blocks;
+
+private:
+   //! Number of block rows/columns
+   int nBlocks;
+   //! Offsets for the starting position of each block
+   Array<int> offsets;
+   //! 2D array that stores each block of the operator.
+   Array2D<Operator *> ops;
+
+   //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
+   //! methods.
+   mutable BlockVector xblock;
+   mutable BlockVector yblock;
+   mutable Vector tmp;
+   mutable Vector tmp2;
+   mutable Vector r;
+   void ForwardPass(const Vector & x, Vector & y) const;
+   void BackwardPass(const Vector & x, Vector & y) const;
+};
+
+
 }
 
 #endif /* MFEM_BLOCKOPERATOR */

linalg/solvers.cpp

@@ -574,6 +574,7 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
    {
       oper->Mult(x, r);
       subtract(b, r, r); // r = b - A x
+      Monitor(0,r.Norml2(),r,x,false);
    }
    else
    {
@@ -623,6 +624,7 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
 
       oper->Mult(x, r);
       subtract(b, r, r); // r = b - A x
+      Monitor(i,nom,r,x,false);
 
       if (prec)
       {
@@ -634,6 +636,8 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
          nom = sqrt(Dot(r, r));
       }
 
+
+
       cf = nom/nomold;
       nomold = nom;
 

mesh/pmesh.cpp

@@ -5602,6 +5602,7 @@ Mesh ParMesh::GetSerialMesh(int save_rank) const
    }
 
    MPI_Barrier(MyComm);
+   serialmesh.SetAttributes();
    return serialmesh;
 }
 

miniapps/contact/contact-vis.cpp

@@ -0,0 +1,562 @@
+//                               contact-visualization test
+
+// srun -n 8 ./contact-vis -testno 6 -paraview -sr 0 -pr 0 -tr 1.0 -sn 1.0
+// srun -n 8 ./contact-vis -testno 6 -paraview -sr 1 -pr 0 -tr 16.0 -sn 1.0
+// srun -n 8 ./contact-vis -testno 6 -paraview -sr 2 -pr 0 -tr 4.0 -sn 1.0
+// srun -n 8 ./contact-vis -testno 6 -paraview -sr 3 -pr 0 -tr 4.0 -sn 1.0
+
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "axom/slic.hpp"
+
+#include "tribol/interface/tribol.hpp"
+#include "tribol/interface/mfem_tribol.hpp"
+#include "tribol/mesh/CouplingScheme.hpp"
+
+using namespace std;
+using namespace mfem;
+
+void GetEssentialTdofs(const ParFiniteElementSpace & fes, const Array<int> & ess_bdr_attr, 
+                       const Array<int> & ess_bdr_attr_comp, Array<int> & ess_tdof_list);
+HypreParMatrix * SetupTribol(ParMesh * pmesh, ParGridFunction * coords, 
+                             const Array<int> & ess_tdofs, const std::set<int> & mortar_attrs, 
+                             const std::set<int> & non_mortar_attrs,
+                             Vector &gap,  double ratio, int tribol_nprocs);
+
+HypreParMatrix * GetContactProlongation(ParFiniteElementSpace * fes, 
+                                        HypreParMatrix *J, 
+                                        Array<int> ess_tdof_list);
+
+int main(int argc, char *argv[])
+{
+    Mpi::Init();
+    int myid = Mpi::WorldRank();
+    int num_procs = Mpi::WorldSize();
+    Hypre::Init();
+
+    int sref = 1;
+    int pref = 0;
+    bool visualization = true;
+    bool paraview = false;
+    bool visit = false;
+    int testNo = -1; // 0-6
+    double tribol_ratio = 8.0;
+    double separation = 0.1;
+    bool disable_essbdr = false;
+    double scale_nodes = 1.0;
+   int tribol_nprocs = num_procs;
+
+    // 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(&sref, "-sr", "--serial-refinements",
+                   "Number of uniform refinements.");                  
+    args.AddOption(&pref, "-pr", "--parallel-refinements",
+                   "Number of uniform refinements.");
+    args.AddOption(&separation, "-s", "--mesh-separation",
+                    "Mesh separation distance.");                   
+    args.AddOption(&tribol_ratio, "-tr", "--tribol-proximity-parameter",
+                   "Tribol-proximity-parameter.");
+    args.AddOption(&scale_nodes, "-sn", "--scale-nodes",
+                   "Scale the nodes of the mesh.");                     
+    args.AddOption(&disable_essbdr, "-noessbdr", "--no-ess-bdr", "-ess-bdr",
+                   "--ess-bdr",
+                   "Enable or disable essential boundary.");                  
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                   "--no-visualization",
+                   "Enable or disable GLVis visualization.");
+    args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
+                   "--no-paraview",
+                   "Enable or disable ParaView visualization.");
+    args.AddOption(&visit, "-visit", "--visit", "-no-visit",
+                    "--no-visit",
+                    "Enable or disable VisIT visualization.");     
+    args.AddOption(&tribol_nprocs, "-tn", "--tribol-nprocs",
+                        "Number of ranks used in tribol redecomposition" );                                  
+    args.Parse();
+    if (!args.Good())
+    {
+        if (myid == 0)
+        {
+            args.PrintUsage(cout);
+        }
+        return 1;
+    }
+    if (myid == 0)
+    {
+        args.PrintOptions(cout);
+    }
+
+    if (Mpi::Root())
+    {
+        mfem::out << "Visualizing test problem number: " << testNo << endl;
+    }
+
+    const char *mesh_file = nullptr;
+
+    switch (testNo)
+    {
+        case 4:
+            mesh_file = "meshes/Test4.mesh";
+            break;
+        case 40:
+            mesh_file = "meshes/Test40.mesh";
+            break;   
+        case 5:
+            mesh_file = "meshes/Test5.mesh";
+            break;
+        case 51:
+            mesh_file = "meshes/Test51.mesh";
+            break;
+        case 6:
+            mesh_file = "meshes/Test6mod2.mesh";
+            break;
+        case -1:
+            mesh_file = "meshes/two-block2.mesh";
+            break;    
+        default:
+            MFEM_ABORT("Should be unreachable");
+            break;
+    }
+
+    Mesh mesh(mesh_file,1);
+
+    for (int i = 0; i<sref; i++)
+    {
+        mesh.UniformRefinement();
+    }
+
+    // move nodes of mesh elements with attr = 2
+    // mesh.EnsureNodes();
+    // GridFunction * nodes = mesh.GetNodes();
+    // Array<int> nodes_marker(nodes->Size());
+    // nodes_marker = 0;
+    // const FiniteElementSpace * meshfes = mesh.GetNodalFESpace();
+    // Array<int> vdofs;
+    // for (int i=0; i<mesh.GetNE(); i++)
+    // {
+    //     int attr = mesh.GetAttribute(i);
+    //     if (attr == 2)
+    //     {
+    //         // get element size
+    //         double h = mesh.GetElementSize(i,0);
+    //         mfem::out << "attr 2 h = " << h << endl;
+    //         meshfes->GetElementVDofs(i,vdofs);
+    //         for (int j = 0; j<vdofs.Size()/3; j++)
+    //         {
+
+    //             int xdof = vdofs[j];
+    //             if (!nodes_marker[xdof])
+    //             {
+    //                 (*nodes)[xdof] += separation;
+    //                 nodes_marker[xdof] = 1;    
+    //             }
+    //         }
+    //     }
+    //     else
+    //     {
+    //         double h = mesh.GetElementSize(i,0);
+    //         mfem::out << "attr 1 h = " << h << endl;
+    //     }
+    // }
+
+
+    ParMesh pmesh(MPI_COMM_WORLD,mesh);
+    mesh.Clear();
+    for (int i = 0; i<pref; i++)
+    {
+        pmesh.UniformRefinement();
+    }
+
+    int dim = pmesh.Dimension();
+
+    std::set<int> mortar_attr;
+    std::set<int> nonmortar_attr;
+    Array<int> ess_bdr_attr;
+    Array<int> ess_bdr_attr_comp;
+
+    // count faces of 9 and 8
+
+    if (testNo == 6)
+    {
+        mortar_attr.insert(6);
+        mortar_attr.insert(9);
+        nonmortar_attr.insert(7);
+        nonmortar_attr.insert(8);
+        
+        // mortar_attr.insert(8);
+        // nonmortar_attr.insert(9);
+
+
+        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
+    {
+        mortar_attr.insert(3);
+        nonmortar_attr.insert(4);
+        ess_bdr_attr.Append(2); ess_bdr_attr_comp.Append(-1);
+        ess_bdr_attr.Append(6); ess_bdr_attr_comp.Append(-1);
+    }
+
+    H1_FECollection fec(1,dim);
+    ParFiniteElementSpace fes(&pmesh,&fec,dim,Ordering::byVDIM);   
+    pmesh.SetNodalFESpace(&fes);
+
+    GridFunction * pnodes = pmesh.GetNodes();
+    for (int i = 0; i < pnodes->Size(); i++)
+    {
+        (*pnodes)(i) *= scale_nodes;
+    }
+
+    int gndofs = fes.GlobalTrueVSize();
+    if (myid == 0)
+    {
+        mfem::out << "--------------------------------------" << endl;
+        mfem::out << "Global number of dofs = " << gndofs << endl;
+        mfem::out << "--------------------------------------" << endl;
+    }
+    ParGridFunction contact_gf(&fes); contact_gf = 0.0;
+    Vector contact_tdofs(fes.GetTrueVSize()); contact_tdofs = 0.0;
+    ParaViewDataCollection * paraview_dc = nullptr;
+    VisItDataCollection * visit_dc = nullptr;
+
+    // Get the essential true dofs
+    Array<int> ess_tdof_list;
+    if (!disable_essbdr)
+    {
+        GetEssentialTdofs(fes, ess_bdr_attr, ess_bdr_attr_comp, ess_tdof_list);
+    }
+    // Set up the Tribol contact problem
+    ParGridFunction ref_coords(&fes);
+    pmesh.GetNodes(ref_coords);
+    Vector gap;
+    HypreParMatrix * J = SetupTribol(&pmesh, &ref_coords, ess_tdof_list, mortar_attr, nonmortar_attr, gap, tribol_ratio, tribol_nprocs);
+
+    // Get the contact prolongation operator
+    HypreParMatrix * Pc = GetContactProlongation(&fes, J, ess_tdof_list);
+
+    if (paraview)
+    {
+        std::ostringstream paraview_file_name;
+        paraview_file_name << "ContactTribolVis_TestNo_" << testNo
+                           << "_par_ref_" << pref
+                           << "_ser_ref_" << sref
+                           << "_tribol-scale_"  << tribol_ratio
+                           << "_mesh_scale_" << scale_nodes;
+        paraview_dc = new ParaViewDataCollection(paraview_file_name.str(), &pmesh);
+        paraview_dc->SetPrefixPath("ParaView");
+        paraview_dc->SetLevelsOfDetail(1);
+        paraview_dc->SetDataFormat(VTKFormat::BINARY32);
+        paraview_dc->SetHighOrderOutput(true);
+        paraview_dc->RegisterField("u_c", &contact_gf);
+        paraview_dc->SetCycle(0);
+        paraview_dc->SetTime(0.0);
+        paraview_dc->Save();
+    }
+    if (visit)
+    {
+        std::ostringstream visit_file_name;
+        visit_file_name << "ContactTribolVis_TestNo_" << testNo
+                           << "_par_ref_" << pref
+                           << "_ser_ref_" << sref
+                           << "_tribol-scale_"  << tribol_ratio
+                           << "_mesh_scale_" << scale_nodes;
+        visit_dc = new VisItDataCollection(visit_file_name.str(), &pmesh);
+        visit_dc->SetPrefixPath("VisIT");
+        visit_dc->RegisterField("u_c", &contact_gf);
+        visit_dc->SetCycle(0);
+        visit_dc->SetTime(0.0);
+        visit_dc->Save();
+    }
+
+
+
+    socketstream sol_sock;
+    if (visualization)
+    {
+        char vishost[] = "localhost";
+        int visport = 19916;
+        sol_sock.open(vishost, visport);
+        sol_sock.precision(8);
+    }
+   
+    int gncols = Pc->GetGlobalNumCols();
+    int gnrows = Pc->GetGlobalNumRows();
+
+    if (Mpi::Root())
+    {
+        mfem::out << "--------------------------------------" << endl;
+        mfem::out << "Global number of contact dofs   = " << gncols << endl;
+        mfem::out << "Global number of contact rows   = " << gnrows << endl;
+        mfem::out << "--------------------------------------" << endl;
+    }
+
+    contact_tdofs = 0.0;
+    Vector Ptc(Pc->Width()); Ptc = 1.0;
+    Pc->Mult(Ptc, contact_tdofs);
+    contact_gf.SetFromTrueDofs(contact_tdofs);
+
+    if (visualization)
+    {
+        sol_sock << "parallel " << num_procs << " " << myid << "\n"
+                 << "solution\n" << pmesh << contact_gf << flush;
+    }
+
+    if (paraview)
+    {
+        paraview_dc->SetCycle(0);
+        paraview_dc->SetTime(0.0);
+        paraview_dc->Save();
+    }
+    if (visit)
+    {
+        visit_dc->SetCycle(0);
+        visit_dc->SetTime(0.0);
+        visit_dc->Save();
+    }
+
+    if (paraview_dc) 
+    {
+        delete paraview_dc;
+    }
+    if (visit_dc) 
+    {
+        delete visit_dc;
+    }
+    return 0;
+}
+
+void GetEssentialTdofs(const ParFiniteElementSpace & fes, const Array<int> & ess_bdr_attr, 
+                       const Array<int> & ess_bdr_attr_comp, Array<int> & ess_tdof_list)
+{
+    ess_tdof_list.SetSize(0);
+    ParMesh * pmesh = fes.GetParMesh();
+    Array<int> ess_bdr;
+    if (pmesh->bdr_attributes.Size())
+    {
+        ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+    }
+    ess_bdr = 0; 
+    Array<int> ess_tdof_list_temp;
+    for (int i = 0; i < ess_bdr_attr.Size(); i++ )
+    {
+        ess_bdr[ess_bdr_attr[i]-1] = 1;
+        fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list_temp, ess_bdr_attr_comp[i]);
+        ess_tdof_list.Append(ess_tdof_list_temp);
+        ess_bdr[ess_bdr_attr[i]-1] = 0;
+    }
+}
+
+HypreParMatrix * SetupTribol(ParMesh * pmesh, ParGridFunction * coords, 
+                             const Array<int> & ess_tdofs, const std::set<int> & mortar_attrs, 
+                             const std::set<int> & non_mortar_attrs,
+                             Vector &gap,  double ratio, int tribol_nprocs)
+{
+    axom::slic::SimpleLogger logger;
+    axom::slic::setIsRoot(mfem::Mpi::Root());
+
+    int coupling_scheme_id = 0;
+    int mesh1_id = 0; int mesh2_id = 1;
+
+    tribol::registerMfemCouplingScheme(
+        coupling_scheme_id, mesh1_id, mesh2_id,
+        *pmesh, *coords, mortar_attrs, non_mortar_attrs,
+        tribol::SURFACE_TO_SURFACE,
+        tribol::NO_SLIDING,
+        tribol::SINGLE_MORTAR,
+        tribol::FRICTIONLESS,
+        tribol::LAGRANGE_MULTIPLIER,
+        tribol::BINNING_GRID
+    );
+
+    tribol::setBinningProximityScale(coupling_scheme_id, ratio);
+    tribol::CouplingSchemeManager::getInstance().findData( coupling_scheme_id )->getParameters().gap_separation_ratio = ratio;
+    // tribol::setMfemLORFactor(coupling_scheme_id, 8);
+    // tribol::setContactAreaFrac( coupling_scheme_id, 1e-8);
+
+    // Access Tribol's pressure grid function (on the contact surface)
+    auto& pressure = tribol::getMfemPressure(coupling_scheme_id);
+    int vsize = pressure.ParFESpace()->GlobalTrueVSize();
+    if (mfem::Mpi::Root())
+    {
+        std::cout << "Number of pressure unknowns: " <<
+        vsize << std::endl;
+    }
+
+    // Set Tribol options for Lagrange multiplier enforcement
+    tribol::setLagrangeMultiplierOptions(
+        coupling_scheme_id,
+        tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
+    );
+
+    // Update contact mesh decomposition
+    tribol::updateMfemParallelDecomposition(tribol_nprocs);
+
+    // Update contact gaps, forces, and tangent stiffness
+    int cycle = 1;   // pseudo cycle
+    double t = 1.0;  // pseudo time
+    double dt = 1.0; // pseudo dt
+    tribol::update(cycle, t, dt);
+
+    // tribol::saveRedecompMesh(0);
+
+    // Return contact contribution to the tangent stiffness matrix
+    auto A_blk = tribol::getMfemBlockJacobian(coupling_scheme_id);
+
+
+
+    HypreParMatrix * Mfull = (HypreParMatrix *)(&A_blk->GetBlock(1,0));
+    
+    // mfem::out << "Mfull size  = " << Mfull->GetGlobalNumRows() << " x " << Mfull->GetGlobalNumCols() << std::endl;
+    
+    HypreParMatrix * Me = Mfull->EliminateCols(ess_tdofs);
+    delete Me;
+
+    int h = Mfull->Height();
+    SparseMatrix merged;
+    Mfull->MergeDiagAndOffd(merged);
+    Array<int> nonzero_rows;
+    for (int i = 0; i<h; i++)
+    {
+        if (!merged.RowIsEmpty(i))
+        {
+            nonzero_rows.Append(i);
+        }
+    }
+
+    int hnew = nonzero_rows.Size();
+
+
+    SparseMatrix P(hnew,h);
+
+    for (int i = 0; i<hnew; i++)
+    {
+        int col = nonzero_rows[i];
+        P.Set(i,col,1.0);
+    }
+    P.Finalize();
+
+    SparseMatrix * reduced_merged = Mult(P,merged);
+
+    int rows[2];
+    int cols[2];
+    cols[0] = Mfull->ColPart()[0];
+    cols[1] = Mfull->ColPart()[1];
+    int nrows = reduced_merged->Height();
+
+    int row_offset;
+    MPI_Scan(&nrows,&row_offset,1,MPI_INT,MPI_SUM,Mfull->GetComm());
+
+    row_offset-=nrows;
+    rows[0] = row_offset;
+    rows[1] = row_offset+nrows;
+    int glob_nrows;
+    MPI_Allreduce(&nrows, &glob_nrows,1,MPI_INT,MPI_SUM,Mfull->GetComm());
+
+    int glob_ncols = reduced_merged->Width();
+    HypreParMatrix * M = new HypreParMatrix(Mfull->GetComm(), nrows, glob_nrows,
+                                            glob_ncols, reduced_merged->GetI(), reduced_merged->GetJ(),
+                                            reduced_merged->GetData(), rows,cols); 
+    delete reduced_merged;                          
+
+    Vector gap_full;
+    tribol::getMfemGap(coupling_scheme_id, gap_full);
+
+    // mfem::out << "gapsize = " << gap_full.Size() << endl;
+    // mfem::out << "gap norm  = " << gap_full.Norml1() << endl;
+
+    //  count zeros in gap
+    int gap_nonzeros=0;
+    for (int i = 0; i<gap_full.Size(); i++)
+    {
+        if (gap_full[i] > 1e-15)
+        {
+            gap_nonzeros++;
+        }
+    }
+    // mfem::out << "gap_nonzeros = " << gap_nonzeros << endl;
+
+    auto& P_submesh = *pressure.ParFESpace()->GetProlongationMatrix();
+    Vector gap_true(P_submesh.Width());
+    P_submesh.MultTranspose(gap_full,gap_true);
+    gap.SetSize(nrows);
+
+    for (int i = 0; i<nrows; i++) 
+    {
+        gap[i] = gap_true[nonzero_rows[i]];
+    }
+    tribol::finalize();
+    return M;
+}
+
+HypreParMatrix * GetContactProlongation(ParFiniteElementSpace * fes, 
+                                        HypreParMatrix *J, 
+                                        Array<int> ess_tdof_list)
+{
+    HypreParMatrix * Jt = J->Transpose();
+    Jt->EliminateRows(ess_tdof_list);
+    int hJt = Jt->Height();
+    SparseMatrix mergedJt;
+    Jt->MergeDiagAndOffd(mergedJt);
+    Array<int> nonzerorows;
+    for (int i = 0; i<hJt; i++)
+    {
+        if (!mergedJt.RowIsEmpty(i))
+        {
+            nonzerorows.Append(i);
+        }
+    }
+    int hc = nonzerorows.Size();
+    SparseMatrix Pct(hc,fes->GlobalTrueVSize());
+
+    // mfem::out << "number of nonzero cols (contact dofs) = " << hc << std::endl;
+
+    for (int i = 0; i<hc; i++)
+    {
+        int col = nonzerorows[i]+fes->GetMyTDofOffset();
+        Pct.Set(i,col,1.0);
+    }
+    Pct.Finalize();
+
+    int rows_c[2];
+    int cols_c[2];
+    int nrows_c = Pct.Height();
+
+    int row_offset_c;
+    MPI_Scan(&nrows_c,&row_offset_c,1,MPI_INT,MPI_SUM,J->GetComm());
+
+    row_offset_c-=nrows_c;
+    rows_c[0] = row_offset_c;
+    rows_c[1] = row_offset_c+nrows_c;
+    for (int i = 0; i < 2; i++)
+    {
+        cols_c[i] = fes->GetTrueDofOffsets()[i];
+    }
+    int glob_nrows_c;
+    int glob_ncols_c = fes->GlobalTrueVSize();
+    MPI_Allreduce(&nrows_c, &glob_nrows_c,1,MPI_INT,MPI_SUM,J->GetComm());
+
+    HypreParMatrix * P_ct = new HypreParMatrix(J->GetComm(), nrows_c, glob_nrows_c,
+                                               glob_ncols_c, Pct.GetI(), Pct.GetJ(),
+                                               Pct.GetData(), rows_c,cols_c); 
+    HypreParMatrix * Pc = P_ct->Transpose();
+    delete P_ct;                         
+    return Pc;
+}

miniapps/contact/ipsolver/ParIPsolver.hpp

@@ -0,0 +1,160 @@
+#include "mfem.hpp"
+#include "../problems/parproblems.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+#ifndef PARIPSOLVER 
+#define PARIPSOLVER
+
+class ParInteriorPointSolver
+{
+protected:
+    OptContactProblem* problem = nullptr;
+    int numActiveConstraints = -1;
+    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 dimG; // num of gap constraints
+    Array<int> constraint_offsets;
+    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 * Wuu = nullptr;
+    HypreParMatrix * Wmm = nullptr;
+    HypreParMatrix * Ju = nullptr;
+    HypreParMatrix * Jm = nullptr;
+    HypreParMatrix * JuT = nullptr;
+    HypreParMatrix * JmT = nullptr;
+
+    Vector Mcslump;
+    Vector Mvlump;
+
+
+    double alphaCurvatureTest;
+    double deltaRegLast;
+    double deltaRegMin;
+    double deltaRegMax;
+    double deltaReg0;
+
+    double kRegMinus;
+    double kRegBarPlus;
+    double kRegPlus;
+
+
+    Array<int> cgnum_iterations;
+    Array<int> amg_num_iterations;
+    Array<double> dmaxmin_ratio;
+    Array<double> jtdj_ratio;
+    Array<double> Adiag_ratio;
+    bool no_contact_solve = false;
+    bool amg_contact_solve = false;
+    Array<int> cgnum_iterations_nocontact;
+    ParFiniteElementSpace *pfes = nullptr;
+    
+    int jOpt;
+    bool converged;
+    
+    int MyRank;
+    bool iAmRoot;
+
+    bool saveLogBarrierIterates = false;
+
+    int linSolver=0;
+    bool dynamicsolver=false;
+    int dynamiclinSolver=0;
+    double linSolveAbsTol = 1e-12;
+    double linSolveRelTol = 1e-6;
+    int relax_type = 88;
+    bool monitor = false;
+    bool save_matrix_data = false;
+    bool useMassWeights = false;
+    int label = -1;
+    MPI_Comm comm;
+public:
+    ParInteriorPointSolver(OptContactProblem*);
+    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 &, double delta = 0.0);
+    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , bool &, double, bool, double delta = 0.0);
+    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> & GetAMGIterNumbers() {return amg_num_iterations;};
+    Array<double> & GetDMaxMinRatios() {return dmaxmin_ratio;};
+    Array<double> & GetJtDJMaxMinRatios() {return jtdj_ratio;};
+    Array<double> & GetAdiagMaxMinRatios() {return Adiag_ratio;};
+    Array<int> & GetCGNoContactIterNumbers() {return cgnum_iterations_nocontact;};
+    int GetNumIterations() {return iter;};
+    double theta(const BlockVector &);
+    double phi(const BlockVector &, double);
+    double phi(const BlockVector &, double, int &);
+    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 SetUsingMassWeights(bool); 
+    void SaveLogBarrierHessianIterates(bool);
+    void SetLinearSolver(int);
+    void SetLinearSolveAbsTol(double);
+    void SetLinearSolveRelTol(double);
+    void SetLinearSolveRelaxType(int);
+    void SetElasticityOptions(ParFiniteElementSpace * pfes_)
+    {
+        pfes = pfes_;
+    };
+    void EnableDynamicSolverChoice() { dynamicsolver = true;};
+    void DisableDynamicSolverChoice() { dynamicsolver = false;};
+    bool CurvatureTest(const BlockOperator & A, const BlockVector & Xhat, const Vector &l, const BlockVector & b, const double & delta);
+    void EnableMonitor() { monitor = true;};
+    void DisableMonitor() { monitor = false;};
+    void EnableSaveMatrix() { save_matrix_data = true;};
+    void DisableSaveMatrix() { save_matrix_data = false;};
+    void EnableNoContactSolve() {no_contact_solve = true;};
+    void EnableAMGContactSolve() {amg_contact_solve = true;};
+    void SetProblemLabel(int label_) { label = label_;};
+    double GetNumActiveConstraints() { return numActiveConstraints;};
+    void Clear() 
+    {
+       F1.DeleteAll();
+       F2.DeleteAll();
+       mu_k = 1.0;
+    };
+    virtual ~ParInteriorPointSolver();
+};
+
+#endif

miniapps/contact/ipsolver/two-level-solver.cpp

@@ -0,0 +1,130 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "two-level-solver.hpp"
+
+namespace mfem
+{
+
+TwoLevelAMGSolver::TwoLevelAMGSolver(MPI_Comm comm_): Solver()
+{
+    Init(comm_);
+}
+
+TwoLevelAMGSolver::TwoLevelAMGSolver(const Operator & Op, const Operator & P_)
+: Solver()
+{
+   auto APtr = dynamic_cast<const HypreParMatrix *>(&Op);
+   MFEM_VERIFY(APtr, "Operator: Not a compatible matrix type");
+   Init(APtr->GetComm());
+   auto PPtr = dynamic_cast<const HypreParMatrix *>(&P_);
+   MFEM_VERIFY(PPtr, "Transfer Map: not a compatible matrix type");
+   
+   SetOperator(Op);
+   SetContactTransferMap(P_);
+}
+
+void TwoLevelAMGSolver::Init(MPI_Comm comm_)
+{
+    comm=comm_;
+    MPI_Comm_size(comm, &numProcs);
+    MPI_Comm_rank(comm, &myid);
+}
+
+void TwoLevelAMGSolver::SetOperator(const Operator & Op)
+{
+    A = dynamic_cast<const HypreParMatrix *>(&Op);
+    height = A->Height();
+    width = A->Width();
+    InitAMG();
+}
+
+void TwoLevelAMGSolver::SetContactTransferMap(const Operator & P)
+{
+    Pc = dynamic_cast<const HypreParMatrix *>(&P);
+    InitCoarseSolver();
+}
+
+void TwoLevelAMGSolver::SetNonContactTransferMap(const Operator & P)
+{
+    Pnc = dynamic_cast<const HypreParMatrix *>(&P);
+}
+
+void TwoLevelAMGSolver::InitAMG()
+{
+    amg = new HypreBoomerAMG(*A);
+    amg->SetPrintLevel(0);
+    amg->SetSystemsOptions(3);
+    amg->SetRelaxType(relax_type);
+}
+
+void TwoLevelAMGSolver::InitCoarseSolver()
+{
+    Ac = RAP(A, Pc);
+#ifdef MFEM_USE_MUMPS
+    Mcoarse = new MUMPSSolver(comm);
+    auto M = dynamic_cast<MUMPSSolver *>(Mcoarse);
+    M->SetPrintLevel(0);
+    M->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE);
+    M->SetOperator(*Ac);
+#else
+#ifdef MFEM_USE_MKL_CPARDISO
+    Mcoarse = new CPardisoSolver(comm);
+    auto M = dynamic_cast<CPardisoSolver *>(Mcoarse);
+    M->SetMatrixType(CPardisoSolver::MatType::REAL_NONSYMMETRIC);
+    M->SetOperator(*Ac);
+#else
+    MFEM_VERIFY(false, "TwoLevelSolver will only work for an mfem build that uses mumps or mkl_cpardiso");
+#endif
+#endif
+}
+
+
+void TwoLevelAMGSolver::Mult(const Vector & b, Vector & x) const
+{
+    MFEM_VERIFY(b.Size() == x.Size(), "Inconsistent x and y size");
+
+    x = 0.0;
+    Vector z(x);
+    amg->Mult(b, z);
+    x+=z;
+    Vector rc(Pc->Width());
+    Vector xc(Pc->Width());
+    if (additive)
+    {
+        Pc->MultTranspose(b,rc);
+        Mcoarse->Mult(rc,xc);
+        Pc->Mult(xc,z);
+    }
+    else
+    {
+        Vector r(b.Size());
+        // 2. Compute Residual r = b - A x
+        A->Mult(x,r);
+        r.Neg(); r+=b;
+        // 3. Restrict to subspace
+        Pc->MultTranspose(r,rc);
+        // 4. Solve on the subspace
+        Mcoarse->Mult(rc,xc);
+        // 5. Transfer to fine space
+        Pc->Mult(xc,z);
+        // 6. Update Correction
+        x+=z;
+        // 7. Compute Residual r = b - A x
+        A->Mult(x,r);
+        r.Neg(); r+=b;
+        // 8. Post V-Cycle 
+        amg->Mult(r, z);
+    }
+    x+= z;
+}
+
+} // namespace

miniapps/contact/ipsolver/two-level-solver.hpp

@@ -0,0 +1,53 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+#ifndef TWOLEVELSOLVER 
+#define TWOLEVELSOLVER
+
+namespace mfem
+{
+
+class TwoLevelAMGSolver : public Solver
+{
+private:
+    MPI_Comm comm;
+    int numProcs, myid;
+    const HypreParMatrix * A = nullptr;
+    const HypreParMatrix * Pc = nullptr;
+    const HypreParMatrix * Pnc = nullptr;
+    HypreBoomerAMG * amg = nullptr;
+    HypreParMatrix * Ac = nullptr;
+    Solver * Mcoarse = nullptr; // previously a mumps solver
+    bool additive = false;
+    int relax_type = 88;
+    void Init(MPI_Comm comm_);
+    void InitAMG();
+    void InitCoarseSolver();
+public:
+    TwoLevelAMGSolver(MPI_Comm comm_);
+    TwoLevelAMGSolver(const Operator & Op, const Operator & P_);
+    void SetOperator(const Operator &op);
+    void SetContactTransferMap(const Operator & P);
+    void SetNonContactTransferMap(const Operator & P);
+    void EnableAdditiveCoupling() { additive = true; }
+    void EnableMultiplicativeCoupling() { additive = false; }
+    void SetAMGRelaxType(int relax_type_) { relax_type = relax_type_;  }
+
+    virtual void Mult(const Vector & y, Vector & x) const; 
+
+    ~TwoLevelAMGSolver()
+    {
+        delete amg;
+        delete Ac;
+        delete Mcoarse;
+    }
+};
+
+}
+
+#endif

miniapps/contact/makefile

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

miniapps/contact/meshes/Test41.mesh

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

miniapps/contact/meshes/Test42.mesh

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

miniapps/contact/meshes/two-block.mesh

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

miniapps/contact/problems/parproblems.cpp

@@ -0,0 +1,818 @@
+#include "parproblems.hpp"
+
+
+ElasticityOperator::ElasticityOperator(ParMesh * pmesh_, Array<int> & ess_bdr_attr_, Array<int> & ess_bdr_attr_comp_,
+                       const Vector & E, const Vector & nu, bool nonlinear_)
+: nonlinear(nonlinear_), pmesh(pmesh_), ess_bdr_attr(ess_bdr_attr_), ess_bdr_attr_comp(ess_bdr_attr_comp_) 
+{
+   comm = pmesh->GetComm();
+   SetParameters(E,nu);
+   Init();
+} 
+
+void ElasticityOperator::SetParameters(const Vector & E, const Vector & nu)
+{
+   int n = (pmesh->attributes.Size()) ?  pmesh->attributes.Max() : 0;
+   MFEM_VERIFY(E.Size() == n, "Incorrect parameter size E");
+   MFEM_VERIFY(nu.Size() == n, "Incorrect parameter size nu");
+   c1.SetSize(n);
+   c2.SetSize(n);
+   if (nonlinear)
+   {
+      for (int i = 0; i<n; i++)
+      {
+         c1(i) = 0.5*E(i) / (1+nu(i));
+         c2(i) = E(i)/(1-2*nu(i))/3;
+      }
+   }
+   else
+   {
+      for (int i = 0; i<n; i++)
+      {
+         c1(i) = E(i) * nu(i) / ( (1+nu(i)) * (1-2*nu(i)) );
+         c2(i) = 0.5 * E(i)/(1+nu(i));
+      }
+   }
+   c1_cf.UpdateConstants(c1);
+   c2_cf.UpdateConstants(c2);
+}
+  
+
+void ElasticityOperator::Init()
+{
+   int dim = pmesh->Dimension();
+   fec = new H1_FECollection(order,dim);
+   fes = new ParFiniteElementSpace(pmesh,fec,dim,Ordering::byVDIM);
+   ndofs = fes->GetVSize();
+   ntdofs = fes->GetTrueVSize();
+   gndofs = fes->GlobalTrueVSize();
+   pmesh->SetNodalFESpace(fes);
+
+   auto ref_func = [](const Vector & x, Vector & y) { y = x; };
+   VectorFunctionCoefficient ref_cf(dim,ref_func);
+   ParGridFunction xr(fes); xr.ProjectCoefficient(ref_cf);
+   xr.GetTrueDofs(xref);
+   SetEssentialBC();
+   SetUpOperator();
+}
+
+void ElasticityOperator::SetEssentialBC()
+{
+   ess_tdof_list.SetSize(0);
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+   }
+   ess_bdr = 0; 
+   Array<int> ess_tdof_list_temp;
+   for (int i = 0; i < ess_bdr_attr.Size(); i++ )
+   {
+      ess_bdr[ess_bdr_attr[i]-1] = 1;
+      fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list_temp,ess_bdr_attr_comp[i]);
+      ess_tdof_list.Append(ess_tdof_list_temp);
+      ess_bdr[ess_bdr_attr[i]-1] = 0;
+   }
+}
+
+void ElasticityOperator::SetUpOperator()
+{
+   x.SetSpace(fes);  x = 0.0;
+   b = new ParLinearForm(fes);
+   if (nonlinear)
+   {
+      material_model = new NeoHookeanModel(c1_cf, c2_cf);
+      op = new ParNonlinearForm(fes);
+      dynamic_cast<ParNonlinearForm*>(op)->AddDomainIntegrator(new HyperelasticNLFIntegrator(material_model));
+      dynamic_cast<ParNonlinearForm*>(op)->SetEssentialTrueDofs(ess_tdof_list);
+   }
+   else
+   {
+      op = new ParBilinearForm(fes);
+      dynamic_cast<ParBilinearForm*>(op)->AddDomainIntegrator(new ElasticityIntegrator(c1_cf,c2_cf));
+      K = new HypreParMatrix();
+      dynamic_cast<ParBilinearForm*>(op)->Assemble();
+      dynamic_cast<ParBilinearForm*>(op)->FormSystemMatrix(ess_tdof_list,*K);
+   }
+}
+
+void ElasticityOperator::FormLinearSystem()
+{
+   if (!formsystem) 
+   {
+      formsystem = true;
+      b->Assemble();
+      B.SetSize(ntdofs);
+      b->ParallelAssemble(B);
+      B.SetSubVector(ess_tdof_list, 0.0);
+      if (!nonlinear)
+      {
+         x.GetTrueDofs(X);
+         dynamic_cast<ParBilinearForm*>(op)->EliminateVDofsInRHS(ess_tdof_list, X, B);
+      }
+   }
+}
+
+void ElasticityOperator::UpdateRHS()
+{
+   formsystem = false;
+   delete b;
+   b = new ParLinearForm(fes);
+}
+
+void ElasticityOperator::SetNeumanPressureData(ConstantCoefficient &f, Array<int> & bdr_marker)
+{
+   pressure_cf.constant = f.constant;
+   b->AddBoundaryIntegrator(new VectorBoundaryFluxLFIntegrator(pressure_cf),bdr_marker);
+}
+
+void ElasticityOperator::SetDisplacementDirichletData(const Vector & delta, Array<int> essbdr)
+{
+   VectorConstantCoefficient delta_cf(delta);
+   x.ProjectBdrCoefficient(delta_cf,essbdr);
+} 
+
+void ElasticityOperator::ResetDisplacementDirichletData() { x = 0.0; }
+
+void ElasticityOperator::UpdateEssentialBC(Array<int> & ess_bdr_attr_, Array<int> & ess_bdr_attr_comp_)
+{
+   ess_bdr_attr = ess_bdr_attr_;
+   ess_bdr_attr_comp = ess_bdr_attr_comp_;
+   SetEssentialBC();
+}
+
+real_t ElasticityOperator::GetEnergy(const Vector & u) const
+{
+   if (nonlinear)
+   {
+      real_t energy = 0.0;
+      Vector tu(xref); tu += u;
+      ParGridFunction u_gf(fes);
+      u_gf.SetFromTrueDofs(tu);
+      energy += dynamic_cast<ParNonlinearForm*>(op)->GetEnergy(u_gf);
+      energy -= InnerProduct(comm, B, u);
+      return energy;
+   }
+   else
+   {
+      Vector ku(K->Height());
+      K->Mult(u,ku);
+      return 0.5 * InnerProduct(comm,u, ku) - InnerProduct(comm,u, B);
+   }
+}
+
+void ElasticityOperator::GetGradient(const Vector & u, Vector & gradE) const
+{
+   if (nonlinear)
+   {
+      Vector tu(xref); tu += u;
+      gradE.SetSize(op->Height());
+      dynamic_cast<ParNonlinearForm*>(op)->Mult(tu, gradE);
+   }
+   else
+   {
+      gradE.SetSize(K->Height());
+      K->Mult(u, gradE);
+   }
+   gradE.Add(-1.0, B); 
+}
+
+HypreParMatrix * ElasticityOperator::GetHessian(const Vector & u)
+{
+   if (nonlinear)
+   {
+      Vector tu(xref); tu += u;
+      return dynamic_cast<HypreParMatrix *>(&dynamic_cast<ParNonlinearForm*>(op)->GetGradient(tu));
+   }
+   else
+   {
+      return K;
+   }
+}
+
+ElasticityOperator::~ElasticityOperator()
+{
+   delete op;
+   delete b;
+   delete fes;
+   delete fec;
+   if (K) delete K;
+   if (material_model) delete material_model;
+}
+
+
+OptContactProblem::OptContactProblem(ElasticityOperator * problem_, 
+                     const std::set<int> & mortar_attrs_, 
+                     const std::set<int> & nonmortar_attrs_,
+                     ParGridFunction * coords_, bool doublepass_,
+                     const Vector & xref_, 
+                     const Vector & xrefbc_, 
+                     double tribol_ratio_,
+                     int tribol_nranks_,
+                     bool qp_,
+		     bool bound_constraints_,
+		     bool mass_weights_)
+: problem(problem_), mortar_attrs(mortar_attrs_), nonmortar_attrs(nonmortar_attrs_),
+  coords(coords_), doublepass(doublepass_), xref(xref_), xrefbc(xrefbc_), 
+   tribol_ratio(tribol_ratio_), tribol_nranks(tribol_nranks_), qp(qp_), 
+	bound_constraints(bound_constraints_), useMassWeights(mass_weights_), block_offsetsg(4)
+{
+   comm = problem->GetComm(); 
+   pmesh = problem->GetMesh();
+   vfes = problem->GetFESpace(); 
+   dim = pmesh->Dimension();
+   ComputeGapJacobian();
+
+   if (problem->IsNonlinear() && qp)
+   {
+      energy_ref = problem->GetEnergy(xrefbc);
+      problem->GetGradient(xrefbc,grad_ref);
+      Kref = problem->GetHessian(xrefbc);
+   }
+   dimU = J->Width();//problem->GetNumTDofs();
+   dimG = J->Height();
+   block_offsetsg[0] = 0;
+   block_offsetsg[1] = dimG;
+   block_offsetsg[2] = dimU;
+   block_offsetsg[3] = dimU;
+   block_offsetsg.PartialSum();
+
+   Vector diagVec(dimU); diagVec = 0.0;
+   SparseMatrix * tempSparse;
+
+   diagVec = 1.0;
+   tempSparse = new SparseMatrix(diagVec);
+   Iu = new HypreParMatrix(comm, GetGlobalNumDofs(), GetDofStarts(), tempSparse);
+   HypreStealOwnership(*Iu, *tempSparse);
+   delete tempSparse;
+
+   diagVec = -1.0;
+   tempSparse = new SparseMatrix(diagVec);
+   negIu = new HypreParMatrix(comm, GetGlobalNumDofs(), GetDofStarts(), tempSparse);
+   HypreStealOwnership(*negIu, *tempSparse);
+   delete tempSparse;
+
+
+   dl.SetSize(dimU); dl = 0.0;
+   eps.SetSize(dimU); eps = 1.e6;
+
+   if (bound_constraints)
+   {
+      dimM = dimG + 2 * dimU;
+   }
+   else
+   {
+      dimM = dimG;
+   }
+   dimC = dimM;
+
+   ml.SetSize(dimM); ml = 0.0;
+
+   MFEM_VERIFY(vfes, "space is null");
+   ParBilinearForm MassForm(vfes);
+   MassForm.AddDomainIntegrator(new VectorMassIntegrator);
+   MassForm.Assemble(); 
+
+   Array<int> empty_tdof_list;
+   Mv = new HypreParMatrix();
+   MassForm.FormSystemMatrix(empty_tdof_list,*Mv);
+
+   Vector onev(Mv->Width()); onev = 1.0;
+   Mvlump.SetSize(Mv->Height());
+   Mv->Mult(onev, Mvlump);
+}
+
+
+
+
+
+void OptContactProblem::ComputeGapJacobian()
+{
+   if (J) delete J;
+   Vector gap1;
+   const HypreParMatrix * J1 = SetupTribol(pmesh,coords,problem->GetEssentialDofs(),
+                                     mortar_attrs, nonmortar_attrs,gap1, tribol_ratio, tribol_nranks);
+   if (doublepass)
+   {
+      Vector gap2;
+      const HypreParMatrix * J2 = SetupTribol(pmesh,coords,problem->GetEssentialDofs(),
+                                       nonmortar_attrs, mortar_attrs, gap2, tribol_ratio, tribol_nranks);
+      gapv.SetSize(gap1.Size()+gap2.Size());
+      gapv.SetVector(gap1,0);
+      gapv.SetVector(gap2,gap1.Size());
+      Array2D<const HypreParMatrix *> A_array(2,1);
+      A_array(0,0) = J1;
+      A_array(1,0) = J2;
+      J = HypreParMatrixFromBlocks(A_array);
+      delete J1;
+      delete J2;
+   }
+   else
+   {
+      gapv.SetSize(gap1.Size());
+      gapv.SetVector(gap1,0);
+      J = const_cast<HypreParMatrix *>(J1);
+   }
+
+   dof_starts.SetSize(2);
+   dof_starts[0] = J->ColPart()[0];
+   dof_starts[1] = J->ColPart()[1];
+   
+   constraints_starts.SetSize(2);
+   if (bound_constraints)
+   {
+      constraints_starts[0] = J->RowPart()[0] + 2 * J->ColPart()[0];
+      constraints_starts[1] = J->RowPart()[1] + 2 * J->ColPart()[1];
+   }
+   else
+   {
+      constraints_starts[0] = J->RowPart()[0];
+      constraints_starts[1] = J->RowPart()[1];
+   }
+}
+
+
+HypreParMatrix * OptContactProblem::Duuf(const BlockVector & x)
+{
+   return DddE(x.GetBlock(0));
+}
+
+HypreParMatrix * OptContactProblem::Dumf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+HypreParMatrix * OptContactProblem::Dmuf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+HypreParMatrix * OptContactProblem::Dmmf(const BlockVector & x)
+{
+   return nullptr;
+}
+
+HypreParMatrix * OptContactProblem::Duc(const BlockVector & x)
+{
+   if (bound_constraints)
+   {
+      Array2D<const HypreParMatrix *> dcduBlockMatrix(3, 1);
+      dcduBlockMatrix(0, 0) = J;
+      dcduBlockMatrix(1, 0) = Iu;
+      dcduBlockMatrix(2, 0) = negIu;
+      if(dcdu)
+      {
+         delete dcdu;
+      }
+      dcdu = HypreParMatrixFromBlocks(dcduBlockMatrix);
+      return dcdu;
+   }
+   else
+   {
+      return J;
+   }
+}
+
+HypreParMatrix * OptContactProblem::Dmc(const BlockVector &)
+{
+   if (!NegId)
+   {
+      Vector negone(dimM); negone = -1.0;
+      SparseMatrix diag(negone);
+      NegId = new HypreParMatrix(comm, GetGlobalNumConstraints(), GetConstraintsStarts(), &diag);
+      HypreStealOwnership(*NegId, diag);
+   }
+   return NegId;
+}
+
+HypreParMatrix * OptContactProblem::lDuuc(const BlockVector &, const Vector &)
+{
+   return nullptr;
+}
+
+HypreParMatrix * OptContactProblem::GetRestrictionToInteriorDofs()
+{
+   if (!Pnc)
+   {
+      if (!Jt) 
+      {
+         Jt = J->Transpose();
+         Jt->EliminateRows(problem->GetEssentialDofs());
+      }
+
+      int hJt = Jt->Height();
+      SparseMatrix mergedJt;
+      Jt->MergeDiagAndOffd(mergedJt);
+      Array<int> zerorows;
+      for (int i = 0; i<hJt; i++)
+      {
+         if (mergedJt.RowIsEmpty(i))
+         {
+            zerorows.Append(i);
+         }
+      }
+
+      int hi = zerorows.Size();
+      SparseMatrix Pit(hi,vfes->GlobalTrueVSize());
+
+      for (int i = 0; i<hi; i++)
+      {
+         int col = zerorows[i]+vfes->GetMyTDofOffset();//prob->GetFESpace()->GetMyTDofOffset();
+         Pit.Set(i,col,1.0);
+      }
+      Pit.Finalize();
+
+      int rows_i[2];
+      int cols_i[2];
+      int nrows_i = Pit.Height();
+
+      int row_offset_i;
+      MPI_Scan(&nrows_i,&row_offset_i,1,MPI_INT,MPI_SUM,comm);
+
+      row_offset_i-=nrows_i;
+      rows_i[0] = row_offset_i;
+      rows_i[1] = row_offset_i+nrows_i;
+      for (int i = 0; i < 2; i++)
+      {
+         cols_i[i] = vfes->GetTrueDofOffsets()[i];
+      }
+      int glob_nrows_i;
+      int glob_ncols_i = vfes->GlobalTrueVSize();
+      MPI_Allreduce(&nrows_i, &glob_nrows_i,1,MPI_INT,MPI_SUM,comm);
+      HypreParMatrix * P_it = new HypreParMatrix(comm, nrows_i, glob_nrows_i,
+                                 glob_ncols_i, Pit.GetI(), Pit.GetJ(),
+                                 Pit.GetData(), rows_i,cols_i); 
+      // HypreStealOwnership(*P_it, Pit);
+      Pnc = P_it->Transpose();
+      delete P_it;
+   }
+
+   return Pnc;
+}
+
+HypreParMatrix * OptContactProblem::GetRestrictionToContactDofs()
+{
+   if (!Pc)
+   {
+      if (!Jt)
+      {
+         Jt = J->Transpose();
+         Jt->EliminateRows(problem->GetEssentialDofs());
+      }
+      int hJt = Jt->Height();
+      SparseMatrix mergedJt;
+      Jt->MergeDiagAndOffd(mergedJt);
+      Array<int> nonzerorows;
+      for (int i = 0; i<hJt; i++)
+      {
+         if (!mergedJt.RowIsEmpty(i))
+         {
+            nonzerorows.Append(i);
+         }
+      }
+      int hc = nonzerorows.Size();
+      SparseMatrix Pct(hc,vfes->GlobalTrueVSize());
+
+      for (int i = 0; i<hc; i++)
+      {
+         int col = nonzerorows[i]+vfes->GetMyTDofOffset();
+         Pct.Set(i,col,1.0);
+      }
+      Pct.Finalize();
+
+      int rows_c[2];
+      int cols_c[2];
+      int nrows_c = Pct.Height();
+
+      int row_offset_c;
+      MPI_Scan(&nrows_c,&row_offset_c,1,MPI_INT,MPI_SUM,comm);
+
+      row_offset_c-=nrows_c;
+      rows_c[0] = row_offset_c;
+      rows_c[1] = row_offset_c+nrows_c;
+      for (int i = 0; i < 2; i++)
+      {
+         cols_c[i] = vfes->GetTrueDofOffsets()[i];
+      }
+      int glob_nrows_c;
+      int glob_ncols_c = vfes->GlobalTrueVSize();
+      MPI_Allreduce(&nrows_c, &glob_nrows_c,1,MPI_INT,MPI_SUM,comm);
+
+      HypreParMatrix * P_ct = new HypreParMatrix(comm, nrows_c, glob_nrows_c,
+                                                glob_ncols_c, Pct.GetI(), Pct.GetJ(),
+                                                Pct.GetData(), rows_c,cols_c); 
+      // HypreStealOwnership(*P_ct, Pct);
+      Pc = P_ct->Transpose();
+      delete P_ct;                         
+   }
+
+   return Pc;
+}
+
+void OptContactProblem::g(const Vector & d, Vector & gd)
+{
+   Vector temp(dimU); temp = 0.0;
+   temp.Set(1.0, d);
+   temp.Add(-1.0, xref);
+   J->Mult(temp, gd);
+   gd.Add(1.0, gapv);
+}
+
+
+//           [     g1(d)      ]  
+// c(d, s) = [ eps + (d - dl) ] - s
+//           [ eps - (d - dl) ]
+
+
+void OptContactProblem::c(const BlockVector & x, Vector & y)
+{
+   const Vector disp  = x.GetBlock(0);
+   const Vector slack = x.GetBlock(1); 
+   
+   if (bound_constraints)
+   {
+      BlockVector yblock(block_offsetsg); yblock = 0.0;
+      
+      g(disp, yblock.GetBlock(0));
+      yblock.GetBlock(1).Set( 1.0, disp );
+      yblock.GetBlock(1).Add(-1.0, dl);
+      yblock.GetBlock(2).Set(-1.0, yblock.GetBlock(1));
+      yblock.GetBlock(1).Add(1.0, eps);
+      yblock.GetBlock(2).Add(1.0, eps);
+      y.Set(1.0, yblock);
+      y.Add(-1.0, slack);
+   }
+   else
+   {
+      g(disp, y);
+      y.Add(-1., slack);
+   }
+}
+
+real_t OptContactProblem::CalcObjective(const BlockVector & x, int & eval_err)
+{
+   return E(x.GetBlock(0), eval_err);
+}
+
+void OptContactProblem::CalcObjectiveGrad(const BlockVector & x, BlockVector & y)
+{
+   DdE(x.GetBlock(0), y.GetBlock(0));
+   y.GetBlock(1) = 0.0;
+}
+
+real_t OptContactProblem::E(const Vector & d, int & eval_err)
+{
+   if (problem->IsNonlinear() && qp)
+   {
+      // (d - xref)^T [ 1/2 K * (d - xref) + gradEQP] + EQP
+      double energy = 0.0;
+      Vector dx(dimU); dx = 0.0;
+      Vector temp(dimU); temp = 0.0;
+      dx.Set(1.0, d);
+      dx.Add(-1.0, xrefbc);
+      Kref->Mult(dx, temp);
+      temp *= 0.5;
+      temp.Add(1.0, grad_ref);
+      energy = InnerProduct(comm, dx, temp);
+      energy += energy_ref;
+      eval_err = 0;
+      return energy;
+   }
+   else
+   {
+      double energy = problem->GetEnergy(d);
+      if (IsFinite(energy))
+      {
+         eval_err = 0;
+      }
+      else
+      {
+         eval_err = 1;
+      }
+      if (Mpi::Root() && eval_err == 1)
+      {
+         cout << "energy = " << energy << endl;
+	 cout << "eval_err = " << eval_err << endl;
+      }
+      return energy;
+   }
+}
+
+
+
+
+void OptContactProblem::DdE(const Vector & d, Vector & gradE)
+{
+   if (problem->IsNonlinear() && qp)
+   {
+      // KQP * (d - xref) + gradEQP
+      Vector dx(dimU); dx = 0.0;
+      dx.Set(1.0, d);
+      dx.Add(-1.0, xrefbc);
+      Kref->Mult(dx, gradE);
+      gradE.Add(1.0, grad_ref);
+   }
+   else
+   {
+      return problem->GetGradient(d, gradE);
+   }
+}
+
+HypreParMatrix * OptContactProblem::DddE(const Vector & d)
+{
+   if (problem->IsNonlinear() && qp)
+   {
+      return Kref;
+   }
+   else
+   {
+      return problem->GetHessian(d);
+   }
+}
+
+void OptContactProblem::SetBoundConstraints(const Vector & dl_, const Vector & eps_)
+{
+   MFEM_VERIFY(dl_.Size() == dimU, "constraint vector dl is not of the correct size");
+   MFEM_VERIFY(eps_.Size() == dimU, "constraint vector eps is not the correct size");
+   dl.Set(1.0, dl_);
+   eps.Set(1.0, eps_);
+}
+
+OptContactProblem::~OptContactProblem()
+{
+   delete J;
+   delete Jt;
+   delete Pc;
+   delete Pnc;
+   delete NegId;
+   delete Iu;
+   delete negIu;
+   delete Mv;
+   delete Mcs;
+   if (dcdu)
+   {
+      delete dcdu;
+   }
+}
+
+HypreParMatrix * OptContactProblem::SetupTribol(ParMesh * pmesh, ParGridFunction * coords, 
+                             const Array<int> & ess_tdofs, const std::set<int> & mortar_attrs, 
+                              const std::set<int> & non_mortar_attrs, 
+                              Vector &gap, double ratio, int tribol_nranks)
+{
+   axom::slic::SimpleLogger logger;
+   axom::slic::setIsRoot(mfem::Mpi::Root());
+
+   // Initialize Tribol contact library
+   tribol::initialize(pmesh->Dimension(), pmesh->GetComm());
+
+   int coupling_scheme_id = 0;
+   int mesh1_id = 0; int mesh2_id = 1;
+   
+   tribol::registerMfemCouplingScheme(
+      coupling_scheme_id, mesh1_id, mesh2_id,
+      *pmesh, *coords, mortar_attrs, non_mortar_attrs,
+      tribol::SURFACE_TO_SURFACE,
+      tribol::NO_SLIDING,
+      tribol::SINGLE_MORTAR,
+      tribol::FRICTIONLESS,
+      tribol::LAGRANGE_MULTIPLIER,
+      tribol::BINNING_GRID
+   );
+
+   tribol::setBinningProximityScale(coupling_scheme_id, ratio);
+   tribol::CouplingSchemeManager::getInstance().findData( coupling_scheme_id )->getParameters().gap_separation_ratio = ratio;
+
+   // Access Tribol's pressure grid function (on the contact surface)
+   auto& pressure = tribol::getMfemPressure(coupling_scheme_id);
+   int vsize = pressure.ParFESpace()->GlobalTrueVSize();
+   
+   ParBilinearForm acs_form(pressure.ParFESpace());
+   acs_form.AddDomainIntegrator(new MassIntegrator);
+   acs_form.Assemble();
+   Array<int> empty_tdof_list;
+   Mcs = new HypreParMatrix();
+   acs_form.FormSystemMatrix(empty_tdof_list,*Mcs);
+   
+   Vector onecs(Mcs->Width()); onecs = 1.0;
+   Mcslumpfull.SetSize(Mcs->Height()); //Vector
+   Mcs->Mult(onecs, Mcslumpfull); 
+
+
+
+   if (mfem::Mpi::Root())
+   {
+      std::cout << "Number of pressure unknowns: " <<
+                vsize << std::endl;
+   }
+
+   // Set Tribol options for Lagrange multiplier enforcement
+   tribol::setLagrangeMultiplierOptions(
+      coupling_scheme_id,
+      tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN
+   );
+
+   // Update contact mesh decomposition
+   //tribol::updateMfemParallelDecomposition(tribol_nranks);
+   tribol::updateMfemParallelDecomposition();
+
+   // Update contact gaps, forces, and tangent stiffness
+   int cycle = 1;   // pseudo cycle
+   double t = 1.0;  // pseudo time
+   double dt = 1.0; // pseudo dt
+   tribol::update(cycle, t, dt);
+
+   // Return contact contribution to the tangent stiffness matrix
+   auto A_blk = tribol::getMfemBlockJacobian(coupling_scheme_id);
+   
+   HypreParMatrix * Mfull = (HypreParMatrix *)(&A_blk->GetBlock(1,0));
+   if (useMassWeights)
+   { 
+      Mfull->InvScaleRows(Mcslumpfull); // scaling
+   }
+   HypreParMatrix * Me = Mfull->EliminateCols(ess_tdofs);
+   delete Me;
+
+   int h = Mfull->Height();
+   SparseMatrix merged;
+   Mfull->MergeDiagAndOffd(merged);
+   Array<int> nonzero_rows;
+   
+   double max_l1_row_norm = 0.0;
+   double rel_row_norm_threshold = 1.e-5;
+   for (int i = 0; i < h; i++)
+   {
+      if (!merged.RowIsEmpty(i))
+      {
+         max_l1_row_norm = max( max_l1_row_norm, merged.GetRowNorml1(i));
+      }
+   }
+   
+   for (int i = 0; i<h; i++)
+   {
+      if (!merged.RowIsEmpty(i))
+      {
+         if (merged.GetRowNorml1(i) > rel_row_norm_threshold * max_l1_row_norm)
+         {
+            nonzero_rows.Append(i);
+         }
+      }
+   }
+
+   int hnew = nonzero_rows.Size();
+   SparseMatrix P(hnew,h);
+
+   for (int i = 0; i<hnew; i++)
+   {
+      int col = nonzero_rows[i];
+      P.Set(i,col,1.0);
+   }
+   P.Finalize();
+
+   SparseMatrix * reduced_merged = Mult(P,merged);
+
+   int rows[2];
+   int cols[2];
+   cols[0] = Mfull->ColPart()[0];
+   cols[1] = Mfull->ColPart()[1];
+   int nrows = reduced_merged->Height();
+
+   int row_offset;
+   MPI_Scan(&nrows,&row_offset,1,MPI_INT,MPI_SUM,Mfull->GetComm());
+
+   row_offset-=nrows;
+   rows[0] = row_offset;
+   rows[1] = row_offset+nrows;
+   int glob_nrows;
+   MPI_Allreduce(&nrows, &glob_nrows,1,MPI_INT,MPI_SUM,Mfull->GetComm());
+
+
+   int glob_ncols = reduced_merged->Width();
+   HypreParMatrix * M = new HypreParMatrix(Mfull->GetComm(), nrows, glob_nrows,
+                          glob_ncols, reduced_merged->GetI(), reduced_merged->GetJ(),
+                          reduced_merged->GetData(), rows,cols); 
+   // HypreStealOwnership(*M, *reduced_merged);
+   delete reduced_merged;                          
+
+   Vector gap_full;
+   tribol::getMfemGap(coupling_scheme_id, gap_full);
+
+   auto& P_submesh = *pressure.ParFESpace()->GetProlongationMatrix();
+   Vector gap_true(P_submesh.Width());
+   
+   
+   P_submesh.MultTranspose(gap_full,gap_true);
+   gap.SetSize(nrows);
+   Mcslump.SetSize(nrows);
+
+   for (int i = 0; i<nrows; i++) 
+   {
+      gap[i] = gap_true[nonzero_rows[i]];
+      Mcslump(i) = Mcslumpfull(nonzero_rows[i]);
+   }
+   if (useMassWeights)
+   {
+      gap /= Mcslump;
+   }
+   tribol::finalize();
+   return M;
+}

miniapps/contact/problems/parproblems.hpp

@@ -0,0 +1,207 @@
+
+#include "parproblems_util.hpp"
+
+class ElasticityOperator
+{
+private:
+   MPI_Comm comm;
+   bool nonlinear = false;
+   bool formsystem = false;
+   ParMesh * pmesh = nullptr;
+   Array<int> ess_bdr_attr, ess_bdr_attr_comp;
+   Array<int> ess_bdr, ess_tdof_list;
+   int order=1, ndofs, ntdofs, gndofs;
+   FiniteElementCollection * fec = nullptr;
+   ParFiniteElementSpace * fes = nullptr;
+   Operator * op = nullptr; // Bilinear or Nonlinear form
+   ParLinearForm * b = nullptr;
+   ParGridFunction x;
+   HypreParMatrix *K=nullptr; // Gradient
+   Vector B, X; // Rhs and Solution vector
+
+   ConstantCoefficient pressure_cf;
+   // linear elasticity:     
+   // c1 = λ (1ˢᵗ Lame parameter), c2 = μ (2ⁿᵈ Lame parameter or shear modulus)
+   // non linear elasticity: 
+   // c1 = G (shear modulus μ ), c2 = K (bulk modulus) 
+   Vector c1, c2;
+   PWConstCoefficient c1_cf, c2_cf;
+   NeoHookeanModel * material_model = nullptr;
+
+   Vector xref;
+   void Init();
+   void SetEssentialBC();
+   void SetUpOperator();
+
+public:
+   ElasticityOperator(ParMesh * pmesh_, Array<int> & ess_bdr_attr_, Array<int> & ess_bdr_attr_comp_,
+                       const Vector & E, const Vector & nu, bool nonlinear_ = false); 
+   void SetParameters(const Vector & E, const Vector & nu); 
+   void SetNeumanPressureData(ConstantCoefficient &f, Array<int> & bdr_marker);
+   void SetDisplacementDirichletData(const Vector & delta, Array<int> essbdr); 
+   void ResetDisplacementDirichletData();
+   void UpdateEssentialBC(Array<int> & ess_bdr_attr_, Array<int> & ess_bdr_attr_comp_);
+   void FormLinearSystem();
+   void UpdateLinearSystem();
+   void UpdateRHS();
+
+   ParMesh * GetMesh() const { return pmesh; };
+   MPI_Comm GetComm() const { return comm; };
+
+   ParFiniteElementSpace * GetFESpace() const { return fes; };
+   const FiniteElementCollection * GetFECol() const { return fec; };
+   int GetNumDofs() const { return ndofs; };
+   int GetNumTDofs() const { return ntdofs; };
+   int GetGlobalNumDofs() const { return gndofs; };
+   const HypreParMatrix * GetOperator() const { return K; };
+   const Vector & GetRHS() const { return B; };
+
+   const ParGridFunction & GetDisplacementGridFunction() const { return x; };
+   const Array<int> & GetEssentialDofs() const { return ess_tdof_list; };
+
+   real_t GetEnergy(const Vector & u) const;
+   void GetGradient(const Vector & u, Vector & gradE) const;
+   HypreParMatrix * GetHessian(const Vector & u);
+   bool IsNonlinear() { return nonlinear; }
+
+   ~ElasticityOperator();
+};
+
+class OptContactProblem
+{
+private:
+   MPI_Comm comm;
+   ElasticityOperator * problem = nullptr;
+   ParFiniteElementSpace * vfes = nullptr;
+   int dim;
+   int dimU, dimM, dimC;
+   int dimG;
+   Vector ml;
+   HypreParMatrix * NegId = nullptr;
+   HypreParMatrix * Kref=nullptr;
+   Vector grad_ref;
+   real_t energy_ref;
+
+   ParMesh * pmesh = nullptr;
+   std::set<int> mortar_attrs;
+   std::set<int> nonmortar_attrs;
+   ParGridFunction * coords = nullptr;
+   bool doublepass = false;
+
+   Vector xref;
+   Vector xrefbc;
+   bool qp;
+
+   void ComputeGapJacobian();
+   Vector gapv;
+   // Jacobian of gap
+   HypreParMatrix * J = nullptr;
+   // Transpose of the Jacobian of gap
+   HypreParMatrix * Jt = nullptr;
+   // Restriction operator to the contact dofs
+   HypreParMatrix * Pc = nullptr;
+   // Restriction operator to the non-contact dofs
+   HypreParMatrix * Pnc = nullptr;
+   Array<HYPRE_BigInt> constraints_starts;
+   Array<HYPRE_BigInt> dof_starts;
+
+
+   // with additional constraints
+   //         [ g ]
+   // g_new = [ eps + (d - dl) ]
+   //         [ eps - (d - dl) ]
+   // there are additional components to the Jacobian
+   //         [ J ] 
+   // J_new = [ I ]
+   //         [-I ]      
+   HypreParMatrix * Iu = nullptr;
+   HypreParMatrix * negIu = nullptr;
+
+   HypreParMatrix * dcdu = nullptr;
+
+   HypreParMatrix * Mv = nullptr; // mass matrix in the volume
+   HypreParMatrix * Mcs = nullptr; // mass matrix on the contact surface
+   Vector Mvlump;
+   Vector Mcslumpfull;
+   Vector Mcslump;
+   bool useMassWeights;
+   
+
+   Vector dl;
+   Vector eps;
+   Array<int> block_offsetsg;
+   bool bound_constraints;
+   int tribol_nranks;
+   double tribol_ratio;
+public:
+   OptContactProblem(ElasticityOperator * problem_, 
+                     const std::set<int> & mortar_attrs_, 
+                     const std::set<int> & nonmortar_attrs_,
+                     ParGridFunction * coords_, bool doublepass_,
+                     const Vector & xref_, 
+                     const Vector & xrefbc_, 
+                     double tribol_ratio_,
+                     int tribol_nranks_,
+                     bool qp_ = true,
+		     bool bound_constraints_=true,
+		     bool mass_weights_=false);
+   int GetDimU() {return dimU;}
+   int GetDimM() {return dimM;}
+   int GetDimC() {return dimC;}
+   Vector & Getml() {return ml;}
+   MPI_Comm GetComm() {return comm ;}
+   HYPRE_BigInt * GetConstraintsStarts() {return constraints_starts.GetData();} 
+   HYPRE_BigInt GetGlobalNumConstraints() {
+	   if (bound_constraints)
+	   {	   
+	      return J->GetGlobalNumRows() + 2 * J->GetGlobalNumCols();
+	   }
+	   else
+	   {
+	      return J->GetGlobalNumRows();
+	   }
+   }
+
+   HYPRE_BigInt * GetDofStarts() {return dof_starts.GetData();}
+   HYPRE_BigInt GetGlobalNumDofs() {return J->GetGlobalNumCols(); }
+   ElasticityOperator * GetElasticityOperator() {return problem;}
+
+   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();
+   HypreParMatrix * GetRestrictionToContactDofs(); 
+
+   void c(const BlockVector &, Vector &);
+   void g(const Vector &, Vector &);
+   double CalcObjective(const BlockVector &, int &);
+   void CalcObjectiveGrad(const BlockVector &, BlockVector &);
+
+   //double E(const Vector & d);
+   double E(const Vector & d, int & eval_err);
+   void DdE(const Vector & d, Vector & gradE);
+   HypreParMatrix * DddE(const Vector & d);
+   
+   void SetBoundConstraints(const Vector & dl_, const Vector & eps_);
+   HypreParMatrix *  SetupTribol(ParMesh * pmesh, ParGridFunction * coords,
+                              const Array<int> & ess_tdofs,
+                              const std::set<int> & mortar_attrs, 
+                              const std::set<int> & non_mortar_attrs, 
+                              Vector &gap,  double tribol_ratio, int tribol_nranks);
+   
+   void GetLumpedMassWeights(Vector & Mcslump_, Vector & Mvlump_)
+   {
+      Mcslump_.SetSize(Mcslump.Size()); Mcslump_ = 0.0;
+      Mcslump_.Set(1.0, Mcslump);
+      Mvlump_.SetSize(Mvlump.Size()); Mvlump_ = 0.0;
+      Mvlump_.Set(1.0, Mvlump);
+   };
+   ~OptContactProblem();
+};
+
+

miniapps/contact/problems/parproblems_util.cpp

@@ -0,0 +1,176 @@
+#include "parproblems_util.hpp"
+
+int get_rank(int tdof, std::vector<int> & tdof_offsets)
+{
+   int size = tdof_offsets.size();
+   if (size == 1) { return 0; }
+   std::vector<int>::iterator up;
+   up=std::upper_bound(tdof_offsets.begin(), tdof_offsets.end(),tdof); //
+   return std::distance(tdof_offsets.begin(),up)-1;
+}
+
+void ComputeTdofOffsets(const ParFiniteElementSpace * pfes,
+                        std::vector<int> & tdof_offsets)
+{
+   MPI_Comm comm = pfes->GetComm();
+   int num_procs;
+   MPI_Comm_size(comm, &num_procs);
+   tdof_offsets.resize(num_procs);
+   int mytoffset = pfes->GetMyTDofOffset();
+   MPI_Allgather(&mytoffset,1,MPI_INT,&tdof_offsets[0],1,MPI_INT,comm);
+}
+
+void ComputeTdofOffsets(MPI_Comm comm, int mytoffset, std::vector<int> & tdof_offsets)
+{
+   int num_procs;
+   MPI_Comm_size(comm,&num_procs);
+   tdof_offsets.resize(num_procs);
+   MPI_Allgather(&mytoffset,1,MPI_INT,&tdof_offsets[0],1,MPI_INT,comm);
+}
+
+void ComputeTdofs(MPI_Comm comm, int mytoffs, std::vector<int> & tdofs)
+{
+   int num_procs;
+   MPI_Comm_size(comm,&num_procs);
+   tdofs.resize(num_procs);
+   MPI_Allgather(&mytoffs,1,MPI_INT,&tdofs,1,MPI_INT,comm);
+}         
+
+
+// Performs Pᵀ * A * P for BlockOperator  P (with blocks as HypreParMatrices)
+// and A a HypreParMatrix, i.e., this handles the special case 
+// where P = [P₁ P₂ ⋅⋅⋅ Pₙ] 
+// C = Pᵀ * A * P 
+void RAP(const HypreParMatrix & A, const BlockOperator & P, 
+         BlockOperator & C)
+{
+   int nblocks = P.NumColBlocks();
+
+   const HypreParMatrix * Pi = nullptr;
+   const HypreParMatrix * Pj = nullptr;
+   HypreParMatrix * PitAPj = nullptr; 
+
+   for (int i = 0; i< nblocks; i++)
+   {
+      if (P.IsZeroBlock(0,i)) continue;
+      Pi = dynamic_cast<const HypreParMatrix*>(&P.GetBlock(0,i));
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (P.IsZeroBlock(0,j)) continue;
+         Pj = dynamic_cast<const HypreParMatrix*>(&P.GetBlock(0,j));
+         if (i == j) 
+         {
+            PitAPj = RAP(&A, Pj); 
+         }
+         else
+         {
+            PitAPj = RAP(Pi, &A, Pj); 
+         }
+         C.SetBlock(i,j,PitAPj);
+      }
+   }
+}
+
+void ParAdd(const BlockOperator & A, const BlockOperator & B, BlockOperator & C)
+{
+   int n = A.NumRowBlocks();
+   int m = A.NumColBlocks();
+   MFEM_VERIFY(B.NumRowBlocks() == n, "Inconsistent number of row blocks");
+   MFEM_VERIFY(B.NumColBlocks() == m, "Inconsistent number of column blocks");
+
+   const HypreParMatrix * a;
+   const HypreParMatrix * b;
+   for (int i = 0; i<n; i++)
+   {
+      for (int j = 0; j<m; j++)
+      {
+         a = nullptr;
+         b = nullptr;
+         if (!A.IsZeroBlock(i,j))
+         {
+            a = dynamic_cast<const HypreParMatrix*>(&A.GetBlock(i,j));
+         }
+         if (!B.IsZeroBlock(i,j))
+         {
+            b = dynamic_cast<const HypreParMatrix*>(&B.GetBlock(i,j));
+         }
+         if (a && b)
+         {
+            C.SetBlock(i,j,ParAdd(a,b));
+         }
+         else if (a)
+         {
+            C.SetBlock(i,j,new HypreParMatrix(*a));
+         }
+         else if (b)
+         {
+            C.SetBlock(i,j,new HypreParMatrix(*b));
+         }
+         else
+         {
+            // do nothing
+         }
+      }
+   }
+}
+
+
+GeneralSolutionMonitor::GeneralSolutionMonitor(ParFiniteElementSpace * fes_, HypreParMatrix * A, Vector & B, int output_rate_) : fes(fes_), output_rate(output_rate_)
+{
+   Solver * directsolver;
+#ifdef MFEM_USE_MUMPS
+   directsolver = new MUMPSSolver(MPI_COMM_WORLD);
+   auto dsolver = dynamic_cast<MUMPSSolver *>(directsolver);
+   dsolver->SetPrintLevel(0);
+   dsolver->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE);
+#else
+#ifdef MFEM_USE_MKL_CPARDISO
+   directsolver = new CPardisoSolver(MPI_COMM_WORLD);
+   auto dsolver = dynamic_cast<CPardisoSolver *>(directsolver);
+   dsolver->SetMatrixType(CPardisoSolver::MatType::REAL_NONSYMMETRIC);
+#else
+   MFEM_ABORT("This solver choice requires compiling with MUMPS or MKL");
+#endif
+#endif
+   
+   Vector X(B.Size()); X = 0.0;
+   directsolver->SetOperator(*A);
+   directsolver->Mult(B, X);
+   delete directsolver;
+   
+   true_gf = new ParGridFunction(fes);
+   error_gf = new ParGridFunction(fes);
+   true_gf->SetFromTrueDofs(X);
+   std::ostringstream paraview_file_name;
+   paraview_file_name << "QPContact-Monitor";
+   ParMesh * pmesh = fes->GetParMesh();
+   pgf = new ParGridFunction(fes);
+   paraview_dc = new ParaViewDataCollection(paraview_file_name.str(), pmesh);
+   paraview_dc->SetPrefixPath("ParaView");
+   paraview_dc->SetLevelsOfDetail(1);
+   paraview_dc->SetDataFormat(VTKFormat::BINARY);
+   paraview_dc->SetHighOrderOutput(true);
+   paraview_dc->RegisterField("u", pgf);
+   paraview_dc->RegisterField("true_u", true_gf);
+   paraview_dc->RegisterField("error", error_gf);
+   mfem::out << "GeneralSolutionMonitor" << endl;
+}
+
+void GeneralSolutionMonitor::MonitorResidual(int it, real_t norm,
+                                             const Vector &r, bool final)
+{
+
+}
+
+void GeneralSolutionMonitor::MonitorSolution(int it, real_t norm, const Vector &x, bool final)
+{
+   if (it%output_rate == 0 || final)
+   {
+      mfem::out << "GeneralSolutionMonitor::Saving iteration: " << it << endl;
+      pgf->SetFromTrueDofs(x);
+      add(1.0, *true_gf, -1.0, *pgf, *error_gf);
+      paraview_dc->SetCycle(counter++);
+      paraview_dc->SetTime(double(it));
+      paraview_dc->Save();
+   }
+}

miniapps/contact/problems/parproblems_util.hpp

@@ -0,0 +1,52 @@
+
+#include "mfem.hpp"
+
+
+using namespace std;
+using namespace mfem;
+
+#include "axom/slic.hpp"
+
+#include "tribol/interface/tribol.hpp"
+#include "tribol/interface/mfem_tribol.hpp"
+#include "tribol/mesh/CouplingScheme.hpp"
+
+int get_rank(int tdof, std::vector<int> & tdof_offsets);
+void ComputeTdofOffsets(const ParFiniteElementSpace * pfes,
+                        std::vector<int> & tdof_offsets);
+void ComputeTdofOffsets(MPI_Comm comm, int mytoffset, std::vector<int> & tdof_offsets);
+void ComputeTdofs(MPI_Comm comm, int mytoffs, std::vector<int> & tdofs);
+
+
+// Performs Pᵀ * A * P for BlockOperator  P (with blocks as HypreParMatrices)
+// and A a HypreParMatrix, i.e., this handles the special case 
+// where P = [P₁ P₂ ⋅⋅⋅ Pₙ] 
+void RAP(const HypreParMatrix & A, const BlockOperator & P, BlockOperator & C);
+void ParAdd(const BlockOperator & A, const BlockOperator & B, BlockOperator & C);
+
+class GeneralSolutionMonitor : public IterativeSolverMonitor
+{
+public:
+   GeneralSolutionMonitor(ParFiniteElementSpace * fes_, HypreParMatrix * A, Vector & B, int output_rate);
+
+   void MonitorResidual(int it, real_t norm, const Vector &r, bool final) override;
+   void MonitorSolution(int it, real_t norm, const Vector &x, bool final) override;
+
+   ~GeneralSolutionMonitor()
+   {
+      delete pgf;
+      delete error_gf;
+      delete true_gf;
+      delete paraview_dc;
+   }
+private:
+
+    ParFiniteElementSpace * fes = nullptr;
+    ParGridFunction * true_gf = nullptr;
+    ParGridFunction * error_gf = nullptr;
+    ParGridFunction * pgf = nullptr;
+    ParaViewDataCollection * paraview_dc = nullptr;
+    int counter = 0;
+    int output_rate;
+};
+