adding no contact AMG solves for every step

fem/bilinearform.cpp

@@ -289,10 +289,9 @@ void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat) const
       return;
    }
 
-   const FiniteElement &fe = *fes->GetFE(i);
-
    if (domain_integs.Size())
    {
+      const FiniteElement &fe = *fes->GetFE(i);
       ElementTransformation *eltrans = fes->GetElementTransformation(i);
       domain_integs[0]->AssembleElementMatrix(fe, *eltrans, elmat);
       for (int k = 1; k < domain_integs.Size(); k++)
@@ -303,18 +302,17 @@ void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat) const
    }
    else
    {
-      const int ndof = fe.GetDof() * fes->GetVDim();
-      elmat.SetSize(ndof);
+      fes->GetElementVDofs(i, vdofs);
+      elmat.SetSize(vdofs.Size());
       elmat = 0.0;
    }
 }
 
 void BilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat) const
 {
-   const FiniteElement &be = *fes->GetBE(i);
-
    if (boundary_integs.Size())
    {
+      const FiniteElement &be = *fes->GetBE(i);
       ElementTransformation *eltrans = fes->GetBdrElementTransformation(i);
       boundary_integs[0]->AssembleElementMatrix(be, *eltrans, elmat);
       for (int k = 1; k < boundary_integs.Size(); k++)
@@ -325,8 +323,8 @@ void BilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat) const
    }
    else
    {
-      const int ndof = be.GetDof() * fes->GetVDim();
-      elmat.SetSize(ndof);
+      fes->GetBdrElementVDofs(i, vdofs);
+      elmat.SetSize(vdofs.Size());
       elmat = 0.0;
    }
 }
@@ -1431,50 +1429,32 @@ void MixedBilinearForm::GetBlocks(Array2D<SparseMatrix *> &blocks) const
    mat->GetBlocks(blocks);
 }
 
-void MixedBilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi)
+void MixedBilinearForm::AddDomainIntegrator (BilinearFormIntegrator * bfi)
 {
-   domain_integs.Append(bfi);
+   domain_integs.Append (bfi);
    domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
 }
 
-void MixedBilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi,
-                                            Array<int> &elem_marker)
+void MixedBilinearForm::AddDomainIntegrator (BilinearFormIntegrator * bfi,
+                                             Array<int> &elem_marker)
 {
-   domain_integs.Append(bfi);
+   domain_integs.Append (bfi);
    domain_integs_marker.Append(&elem_marker);
 }
 
-void MixedBilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
+void MixedBilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi)
 {
-   boundary_integs.Append(bfi);
+   boundary_integs.Append (bfi);
    boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
 }
 
-void MixedBilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi,
-                                              Array<int> &bdr_marker)
+void MixedBilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
+                                               Array<int> &bdr_marker)
 {
-   boundary_integs.Append(bfi);
+   boundary_integs.Append (bfi);
    boundary_integs_marker.Append(&bdr_marker);
 }
 
-void MixedBilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
-{
-   interior_face_integs.Append(bfi);
-}
-
-void MixedBilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
-{
-   boundary_face_integs.Append(bfi);
-   boundary_face_integs_marker.Append(NULL); // NULL marker means apply everywhere
-}
-
-void MixedBilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
-                                             Array<int> &bdr_marker)
-{
-   boundary_face_integs.Append(bfi);
-   boundary_face_integs_marker.Append(&bdr_marker);
-}
-
 void MixedBilinearForm::AddTraceFaceIntegrator (BilinearFormIntegrator * bfi)
 {
    trace_face_integs.Append (bfi);
@@ -1607,108 +1587,6 @@ void MixedBilinearForm::Assemble(int skip_zeros)
       }
    }
 
-   if (interior_face_integs.Size())
-   {
-      FaceElementTransformations *ftr;
-      Array<int> trial_vdofs2, test_vdofs2;
-      const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
-
-      int nfaces = mesh->GetNumFaces();
-      for (int i = 0; i < nfaces; i++)
-      {
-         ftr = mesh->GetInteriorFaceTransformations(i);
-         if (ftr != NULL)
-         {
-            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
-            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
-            trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
-            test_fe1 = test_fes->GetFE(ftr->Elem1No);
-            if (ftr->Elem2No >= 0)
-            {
-               trial_fes->GetElementVDofs(ftr->Elem2No, trial_vdofs2);
-               test_fes->GetElementVDofs(ftr->Elem2No, test_vdofs2);
-               trial_vdofs.Append(trial_vdofs2);
-               test_vdofs.Append(test_vdofs2);
-               trial_fe2 = trial_fes->GetFE(ftr->Elem2No);
-               test_fe2 = test_fes->GetFE(ftr->Elem2No);
-            }
-            else
-            {
-               // The test_fe2 object is really a dummy and not used on the
-               // boundaries, but we can't dereference a NULL pointer, and we don't
-               // want to actually make a fake element.
-               trial_fe2 = trial_fe1;
-               test_fe2 = test_fe1;
-            }
-            for (int k = 0; k < interior_face_integs.Size(); k++)
-            {
-               interior_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
-                                                           *test_fe2,
-                                                           *ftr, elemmat);
-               mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
-            }
-         }
-      }
-   }
-
-   if (boundary_face_integs.Size())
-   {
-      FaceElementTransformations *ftr;
-      Array<int> tr_vdofs2, te_vdofs2;
-      const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
-
-      // Which boundary attributes need to be processed?
-      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
-                                 mesh->bdr_attributes.Max() : 0);
-      bdr_attr_marker = 0;
-      for (int k = 0; k < boundary_face_integs.Size(); k++)
-      {
-         if (boundary_face_integs_marker[k] == NULL)
-         {
-            bdr_attr_marker = 1;
-            break;
-         }
-         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
-         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
-                     "invalid boundary marker for boundary face integrator #"
-                     << k << ", counting from zero");
-         for (int i = 0; i < bdr_attr_marker.Size(); i++)
-         {
-            bdr_attr_marker[i] |= bdr_marker[i];
-         }
-      }
-
-      for (int i = 0; i < trial_fes -> GetNBE(); i++)
-      {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
-         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
-
-         ftr = mesh -> GetBdrFaceTransformations (i);
-         if (ftr != NULL)
-         {
-            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
-            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
-            trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
-            test_fe1 = test_fes->GetFE(ftr->Elem1No);
-            // The test_fe2 object is really a dummy and not used on the
-            // boundaries, but we can't dereference a NULL pointer, and we don't
-            // want to actually make a fake element.
-            trial_fe2 = trial_fe1;
-            test_fe2 = test_fe1;
-            for (int k = 0; k < boundary_face_integs.Size(); k++)
-            {
-               if (boundary_face_integs_marker[k] &&
-                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
-
-               boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
-                                                           *test_fe2,
-                                                           *ftr, elemmat);
-               mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
-            }
-         }
-      }
-   }
-
    if (trace_face_integs.Size())
    {
       FaceElementTransformations *ftr;
@@ -1889,11 +1767,10 @@ void MixedBilinearForm::ConformingAssemble()
 
 void MixedBilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat) const
 {
-   const FiniteElement &trial_fe = *trial_fes->GetFE(i);
-   const FiniteElement &test_fe = *test_fes->GetFE(i);
-
    if (domain_integs.Size())
    {
+      const FiniteElement &trial_fe = *trial_fes->GetFE(i);
+      const FiniteElement &test_fe = *test_fes->GetFE(i);
       ElementTransformation *eltrans = test_fes->GetElementTransformation(i);
       domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
                                                elmat);
@@ -1906,21 +1783,19 @@ void MixedBilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat) const
    }
    else
    {
-      const int tr_dofs = trial_fe.GetDof() * trial_fes->GetVDim();
-      const int te_dofs = test_fe.GetDof() * test_fes->GetVDim();
-
-      elmat.SetSize(te_dofs, tr_dofs);
+      trial_fes->GetElementVDofs(i, trial_vdofs);
+      test_fes->GetElementVDofs(i, test_vdofs);
+      elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
       elmat = 0.0;
    }
 }
 
 void MixedBilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat) const
 {
-   const FiniteElement &trial_be = *trial_fes->GetBE(i);
-   const FiniteElement &test_be = *test_fes->GetBE(i);
-
    if (boundary_integs.Size())
    {
+      const FiniteElement &trial_be = *trial_fes->GetBE(i);
+      const FiniteElement &test_be = *test_fes->GetBE(i);
       ElementTransformation *eltrans = test_fes->GetBdrElementTransformation(i);
       boundary_integs[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
                                                  elmat);
@@ -1933,103 +1808,9 @@ void MixedBilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat) const
    }
    else
    {
-      const int tr_dofs = trial_be.GetDof() * trial_fes->GetVDim();
-      const int te_dofs = test_be.GetDof() * test_fes->GetVDim();
-
-      elmat.SetSize(te_dofs, tr_dofs);
-      elmat = 0.0;
-   }
-}
-
-void MixedBilinearForm::ComputeFaceMatrix(int i, DenseMatrix &elmat) const
-{
-   FaceElementTransformations *ftr;
-   Mesh *mesh = test_fes -> GetMesh();
-   ftr = mesh->GetFaceElementTransformations(i);
-   MFEM_ASSERT(ftr, "No associated face transformations.");
-
-   const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
-
-   trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
-   test_fe1 = test_fes->GetFE(ftr->Elem1No);
-   if (ftr->Elem2No >= 0)
-   {
-      trial_fe2 = trial_fes->GetFE(ftr->Elem2No);
-      test_fe2 = test_fes->GetFE(ftr->Elem2No);
-   }
-   else
-   {
-      // The test_fe2 object is really a dummy and not used on the
-      // boundaries, but we can't dereference a NULL pointer, and we don't
-      // want to actually make a fake element.
-      trial_fe2 = trial_fe1;
-      test_fe2 = test_fe1;
-   }
-
-   if (interior_face_integs.Size())
-   {
-      interior_face_integs[0]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
-                                                  *test_fe2,
-                                                  *ftr, elmat);
-      for (int k = 1; k < interior_face_integs.Size(); k++)
-      {
-         interior_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
-                                                     *test_fe2,
-                                                     *ftr, elemmat);
-         elmat += elemmat;
-      }
-   }
-   else
-   {
-      int tr_dofs = trial_fe1->GetDof() * trial_fes->GetVDim();
-      int te_dofs = test_fe1->GetDof() * test_fes->GetVDim();
-      if (ftr->Elem2No >= 0)
-      {
-         tr_dofs += trial_fe2->GetDof() * trial_fes->GetVDim();
-         te_dofs += test_fe2->GetDof() * test_fes->GetVDim();
-      }
-
-      elmat.SetSize(te_dofs, tr_dofs);
-      elmat = 0.0;
-   }
-}
-
-void MixedBilinearForm::ComputeBdrFaceMatrix(int i, DenseMatrix &elmat) const
-{
-   FaceElementTransformations *ftr;
-   Mesh *mesh = test_fes -> GetMesh();
-   ftr = mesh->GetBdrFaceTransformations(i);
-   MFEM_ASSERT(ftr, "No associated boundary face.");
-
-   const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
-
-   trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
-   test_fe1 = test_fes->GetFE(ftr->Elem1No);
-   // The test_fe2 object is really a dummy and not used on the
-   // boundaries, but we can't dereference a NULL pointer, and we don't
-   // want to actually make a fake element.
-   trial_fe2 = trial_fe1;
-   test_fe2 = test_fe1;
-
-   if (boundary_face_integs.Size())
-   {
-      boundary_face_integs[0]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
-                                                  *test_fe2,
-                                                  *ftr, elmat);
-      for (int k = 1; k < boundary_face_integs.Size(); k++)
-      {
-         boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
-                                                     *test_fe2,
-                                                     *ftr, elemmat);
-         elmat += elemmat;
-      }
-   }
-   else
-   {
-      const int tr_dofs = trial_fe1->GetDof() * trial_fes->GetVDim();
-      const int te_dofs = test_fe1->GetDof() * test_fes->GetVDim();
-
-      elmat.SetSize(te_dofs, tr_dofs);
+      trial_fes->GetBdrElementVDofs(i, trial_vdofs);
+      test_fes->GetBdrElementVDofs(i, test_vdofs);
+      elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
       elmat = 0.0;
    }
 }
@@ -2333,10 +2114,6 @@ MixedBilinearForm::~MixedBilinearForm()
       for (i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
       for (i = 0; i < boundary_integs.Size(); i++)
       { delete boundary_integs[i]; }
-      for (i = 0; i < interior_face_integs.Size(); i++)
-      { delete interior_face_integs[i]; }
-      for (i = 0; i < boundary_face_integs.Size(); i++)
-      { delete boundary_face_integs[i]; }
       for (i = 0; i < trace_face_integs.Size(); i++)
       { delete trace_face_integs[i]; }
       for (i = 0; i < boundary_trace_face_integs.Size(); i++)

fem/bilinearform.hpp

@@ -772,14 +772,6 @@ protected:
    /// Entries are not owned.
    Array<Array<int>*> boundary_integs_marker;
 
-   /// Interior face integrators.
-   Array<BilinearFormIntegrator*> interior_face_integs;
-
-   /// Boundary face integrators.
-   Array<BilinearFormIntegrator*> boundary_face_integs;
-   /// Entries are not owned.
-   Array<Array<int>*> boundary_face_integs_marker;
-
    /// Trace face (skeleton) integrators.
    Array<BilinearFormIntegrator*> trace_face_integs;
 
@@ -901,16 +893,6 @@ public:
    void AddBoundaryIntegrator(BilinearFormIntegrator * bfi,
                               Array<int> &bdr_marker);
 
-   /// Adds an interior face integrator. Assumes ownership of @a bfi.
-   void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi);
-
-   /// Adds a boundary face integrator. Assumes ownership of @a bfi.
-   void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi);
-
-   /// Adds a boundary face integrator. Assumes ownership of @a bfi.
-   void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
-                             Array<int> &bdr_marker);
-
    /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
 
        This type of integrator assembles terms over all faces of the mesh using
@@ -941,16 +923,6 @@ public:
        corresponding pointer (to Array<int>) will be NULL. */
    Array<Array<int>*> *GetBBFI_Marker() { return &boundary_integs_marker; }
 
-   /// Access all integrators added with AddInteriorFaceIntegrator().
-   Array<BilinearFormIntegrator*> *GetFBFI() { return &interior_face_integs; }
-
-   /// Access all integrators added with AddBdrFaceIntegrator().
-   Array<BilinearFormIntegrator*> *GetBFBFI() { return &boundary_face_integs; }
-   /** @brief Access all boundary markers added with AddBdrFaceIntegrator().
-       If no marker was specified when the integrator was added, the
-       corresponding pointer (to Array<int>) will be NULL. */
-   Array<Array<int>*> *GetBFBFI_Marker() { return &boundary_face_integs_marker; }
-
    /// Access all integrators added with AddTraceFaceIntegrator().
    Array<BilinearFormIntegrator*> *GetTFBFI() { return &trace_face_integs; }
 
@@ -1016,13 +988,6 @@ public:
    /** @note The boundary attribute markers of the integrators are ignored. */
    void ComputeBdrTraceFaceMatrix(int i, DenseMatrix &elmat) const;
 
-   /// Compute the face matrix of the given face element
-   void ComputeFaceMatrix(int i, DenseMatrix &elmat) const;
-
-   /// Compute the boundary face matrix of the given boundary element
-   /** @note The boundary attribute markers of the integrators are ignored. */
-   void ComputeBdrFaceMatrix(int i, DenseMatrix &elmat) const;
-
    /// Assemble the given element matrix
    /** The element matrix @a elmat is assembled for the element @a i, i.e.
        added to the system matrix. The flag @a skip_zeros skips the zero

fem/bilininteg.cpp

@@ -170,16 +170,6 @@ void BilinearFormIntegrator::AssembleFaceMatrix(
               "   is not implemented for this class.");
 }
 
-void BilinearFormIntegrator::AssembleFaceMatrix(
-   const FiniteElement &trial_fe1, const FiniteElement &test_fe1,
-   const FiniteElement &trial_fe2, const FiniteElement &test_fe2,
-   FaceElementTransformations &Trans,
-   DenseMatrix &elmat)
-{
-   MFEM_ABORT("AssembleFaceMatrix (mixed form) is not implemented for this"
-              " Integrator class.");
-}
-
 void BilinearFormIntegrator::AssembleFaceMatrix(
    const FiniteElement &trial_face_fe, const FiniteElement &test_fe1,
    const FiniteElement &test_fe2, FaceElementTransformations &Trans,
@@ -233,38 +223,28 @@ void TransposeIntegrator::SetIntRule(const IntegrationRule *ir)
    bfi->SetIntRule(ir);
 }
 
-void TransposeIntegrator::AssembleElementMatrix(
+void TransposeIntegrator::AssembleElementMatrix (
    const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
 {
-   bfi->AssembleElementMatrix(el, Trans, bfi_elmat);
+   bfi -> AssembleElementMatrix (el, Trans, bfi_elmat);
    // elmat = bfi_elmat^t
    elmat.Transpose (bfi_elmat);
 }
 
-void TransposeIntegrator::AssembleElementMatrix2(
+void TransposeIntegrator::AssembleElementMatrix2 (
    const FiniteElement &trial_fe, const FiniteElement &test_fe,
    ElementTransformation &Trans, DenseMatrix &elmat)
 {
-   bfi->AssembleElementMatrix2(test_fe, trial_fe, Trans, bfi_elmat);
+   bfi -> AssembleElementMatrix2 (test_fe, trial_fe, Trans, bfi_elmat);
    // elmat = bfi_elmat^t
    elmat.Transpose (bfi_elmat);
 }
 
-void TransposeIntegrator::AssembleFaceMatrix(
+void TransposeIntegrator::AssembleFaceMatrix (
    const FiniteElement &el1, const FiniteElement &el2,
    FaceElementTransformations &Trans, DenseMatrix &elmat)
 {
-   bfi->AssembleFaceMatrix(el1, el2, Trans, bfi_elmat);
-   // elmat = bfi_elmat^t
-   elmat.Transpose (bfi_elmat);
-}
-
-void TransposeIntegrator::AssembleFaceMatrix(
-   const FiniteElement &tr_el1, const FiniteElement &te_el1,
-   const FiniteElement &tr_el2, const FiniteElement &te_el2,
-   FaceElementTransformations &Trans, DenseMatrix &elmat)
-{
-   bfi->AssembleFaceMatrix(te_el1, tr_el1, te_el2, tr_el2, Trans, bfi_elmat);
+   bfi -> AssembleFaceMatrix (el1, el2, Trans, bfi_elmat);
    // elmat = bfi_elmat^t
    elmat.Transpose (bfi_elmat);
 }
@@ -3518,150 +3498,6 @@ void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
    }
 }
 
-void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &trial_fe1,
-                                           const FiniteElement &test_fe1,
-                                           const FiniteElement &trial_fe2,
-                                           const FiniteElement &test_fe2,
-                                           FaceElementTransformations &Trans,
-                                           DenseMatrix &elmat)
-{
-   int tr_ndof1, te_ndof1, tr_ndof2, te_ndof2;
-
-   real_t un, a, b, w;
-
-   dim = test_fe1.GetDim();
-   tr_ndof1 = trial_fe1.GetDof();
-   te_ndof1 = test_fe1.GetDof();
-   Vector vu(dim), nor(dim);
-
-   if (Trans.Elem2No >= 0)
-   {
-      tr_ndof2 = trial_fe2.GetDof();
-      te_ndof2 = test_fe2.GetDof();
-   }
-   else
-   {
-      tr_ndof2 = 0;
-      te_ndof2 = 0;
-   }
-
-   tr_shape1.SetSize(tr_ndof1);
-   te_shape1.SetSize(te_ndof1);
-   tr_shape2.SetSize(tr_ndof2);
-   te_shape2.SetSize(te_ndof2);
-   elmat.SetSize(te_ndof1 + te_ndof2, tr_ndof1 + tr_ndof2);
-   elmat = 0.0;
-
-   const IntegrationRule *ir = IntRule;
-   if (ir == NULL)
-   {
-      int order;
-      // Assuming order(u)==order(mesh)
-      if (Trans.Elem2No >= 0)
-         order = (min(Trans.Elem1->OrderW(), Trans.Elem2->OrderW()) +
-                  max(trial_fe1.GetOrder(), trial_fe2.GetOrder()) +
-                  max(test_fe1.GetOrder(), test_fe2.GetOrder()));
-      else
-      {
-         order = Trans.Elem1->OrderW() + trial_fe1.GetOrder() + test_fe1.GetOrder();
-      }
-      if (trial_fe1.Space() == FunctionSpace::Pk)
-      {
-         order++;
-      }
-      ir = &IntRules.Get(Trans.FaceGeom, order);
-   }
-
-   for (int p = 0; p < ir->GetNPoints(); p++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(p);
-      IntegrationPoint eip1, eip2;
-      Trans.Loc1.Transform(ip, eip1);
-      Trans.Elem1->SetIntPoint(&eip1);
-      if (tr_ndof2 && te_ndof2)
-      {
-         Trans.Loc2.Transform(ip, eip2);
-         Trans.Elem2->SetIntPoint(&eip2);
-      }
-      trial_fe1.CalcPhysShape(*Trans.Elem1, tr_shape1);
-      test_fe1.CalcPhysShape(*Trans.Elem1, te_shape1);
-
-      Trans.Face->SetIntPoint(&ip);
-
-      u->Eval(vu, *Trans.Elem1, eip1);
-
-      if (dim == 1)
-      {
-         nor(0) = 2*eip1.x - 1.0;
-      }
-      else
-      {
-         CalcOrtho(Trans.Face->Jacobian(), nor);
-      }
-
-      un = vu * nor;
-      a = 0.5 * alpha * un;
-      b = beta * fabs(un);
-      // note: if |alpha/2|==|beta| then |a|==|b|, i.e. (a==b) or (a==-b)
-      //       and therefore two blocks in the element matrix contribution
-      //       (from the current quadrature point) are 0
-
-      if (rho)
-      {
-         real_t rho_p;
-         if (un >= 0.0 && tr_ndof2 && te_ndof2)
-         {
-            Trans.Elem2->SetIntPoint(&eip2);
-            rho_p = rho->Eval(*Trans.Elem2, eip2);
-         }
-         else
-         {
-            rho_p = rho->Eval(*Trans.Elem1, eip1);
-         }
-         a *= rho_p;
-         b *= rho_p;
-      }
-
-      w = ip.weight * (a+b);
-      if (w != 0.0)
-      {
-         for (int i = 0; i < te_ndof1; i++)
-            for (int j = 0; j < tr_ndof1; j++)
-            {
-               elmat(i, j) += w * te_shape1(i) * tr_shape1(j);
-            }
-      }
-
-      if (tr_ndof2 && te_ndof2)
-      {
-         trial_fe2.CalcPhysShape(*Trans.Elem2, tr_shape2);
-         test_fe2.CalcPhysShape(*Trans.Elem2, te_shape2);
-
-         if (w != 0.0)
-            for (int i = 0; i < te_ndof2; i++)
-               for (int j = 0; j < tr_ndof1; j++)
-               {
-                  elmat(te_ndof1+i, j) -= w * te_shape2(i) * tr_shape1(j);
-               }
-
-         w = ip.weight * (b-a);
-         if (w != 0.0)
-         {
-            for (int i = 0; i < te_ndof2; i++)
-               for (int j = 0; j < tr_ndof2; j++)
-               {
-                  elmat(te_ndof1+i, tr_ndof1+j) += w * te_shape2(i) * tr_shape2(j);
-               }
-
-            for (int i = 0; i < te_ndof1; i++)
-               for (int j = 0; j < tr_ndof2; j++)
-               {
-                  elmat(i, tr_ndof1+j) -= w * te_shape1(i) * tr_shape2(j);
-               }
-         }
-      }
-   }
-}
 
 const IntegrationRule &DGTraceIntegrator::GetRule(
    Geometry::Type geom, int order, FaceElementTransformations &T)

fem/bilininteg.hpp

@@ -159,13 +159,6 @@ public:
                                    FaceElementTransformations &Trans,
                                    DenseMatrix &elmat);
 
-   virtual void AssembleFaceMatrix(const FiniteElement &trial_fe1,
-                                   const FiniteElement &test_fe1,
-                                   const FiniteElement &trial_fe2,
-                                   const FiniteElement &test_fe2,
-                                   FaceElementTransformations &Trans,
-                                   DenseMatrix &elmat);
-
    /** Abstract method used for assembling TraceFaceIntegrators in a
        MixedBilinearForm. */
    virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe,
@@ -342,13 +335,6 @@ public:
                            FaceElementTransformations &Trans,
                            DenseMatrix &elmat) override;
 
-   void AssembleFaceMatrix(const FiniteElement &trial_fe1,
-                           const FiniteElement &test_fe1,
-                           const FiniteElement &trial_fe2,
-                           const FiniteElement &test_fe2,
-                           FaceElementTransformations &Trans,
-                           DenseMatrix &elmat) override;
-
    void AssemblePA(const FiniteElementSpace& fes) override
    {
       bfi->AssemblePA(fes);
@@ -3246,7 +3232,6 @@ protected:
 
 private:
    Vector shape1, shape2;
-   Vector tr_shape1, te_shape1, tr_shape2, te_shape2;
 
 public:
    /// Construct integrator with $\rho = 1$, $\beta = \alpha/2$.
@@ -3267,13 +3252,6 @@ public:
                            FaceElementTransformations &Trans,
                            DenseMatrix &elmat) override;
 
-   void AssembleFaceMatrix(const FiniteElement &trial_fe1,
-                           const FiniteElement &test_fe1,
-                           const FiniteElement &trial_fe2,
-                           const FiniteElement &test_fe2,
-                           FaceElementTransformations &Trans,
-                           DenseMatrix &elmat) override;
-
    void AssemblePAInteriorFaces(const FiniteElementSpace &fes) override;
 
    void AssemblePABoundaryFaces(const FiniteElementSpace &fes) override;

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

@@ -37,4 +37,3 @@ add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
 add_subdirectory(hdiv-linear-solver)
-add_subdirectory(interiorpointsolver)

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;
+};
+

miniapps/interiorpointsolver/IPsolver.cpp

@@ -1,964 +0,0 @@
-#include "mfem.hpp"
-#include "IPsolver.hpp"
-#include "Problem.hpp"
-#include <fstream>
-#include <iostream>
-#include <cstdlib>
-
-using namespace std;
-using namespace mfem;
-
-
-ParInteriorPointSolver::ParInteriorPointSolver(ParGeneralOptProblem * problem_) 
-                     : problem(problem_), 
-                       block_offsetsumlz(5), block_offsetsuml(4), block_offsetsx(3),
-                       Huu(nullptr), Hum(nullptr), Hmu(nullptr), 
-                       Hmm(nullptr), Wmm(nullptr), D(nullptr), 
-                       Ju(nullptr), Jm(nullptr), JuT(nullptr), JmT(nullptr), 
-                       saveIterates(false)
-{
-   OptTol  = 1.e-2;
-   max_iter = 20;
-   mu_k     = 1.0;
-
-   sMax     = 1.e2;
-   kSig     = 1.e10;   // control deviation from primal Hessian
-   tauMin   = 0.99;     // control rate at which iterates can approach the boundary
-   eta      = 1.e-4;   // backtracking constant
-   thetaMin = 1.e-4;   // allowed violation of the equality constraints
-
-   // constants in line-step A-5.4
-   delta    = 1.0;
-   sTheta   = 1.1;
-   sPhi     = 2.3;
-
-   // control the rate at which the penalty parameter is decreased
-   kMu     = 0.2;
-   thetaMu = 1.5;
-
-   thetaMax = 1.e6; // maximum constraint violation
-   // data for the second order correction
-   kSoc     = 0.99;
-
-   // equation (18)
-   gTheta = 1.e-5;
-   gPhi   = 1.e-5;
-
-   kEps   = 1.e1;
-
-   dimU = problem->GetDimU();
-   dimM = problem->GetDimM();
-   dimC = problem->GetDimC();
-   ckSoc.SetSize(dimC);
-  
-   block_offsetsumlz[0] = 0;
-   block_offsetsumlz[1] = dimU; // u
-   block_offsetsumlz[2] = dimM; // m
-   block_offsetsumlz[3] = dimC; // lambda
-   block_offsetsumlz[4] = dimM; // zl
-   block_offsetsumlz.PartialSum();
-   
-   MPI_Allreduce(&dimU, &dimUGlb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-   MPI_Allreduce(&dimM, &dimMGlb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-   MPI_Allreduce(&dimC, &dimCGlb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-
-
-   for(int i = 0; i < block_offsetsuml.Size(); i++)  
-   { 
-      block_offsetsuml[i] = block_offsetsumlz[i]; 
-   }
-   for(int i = 0; i < block_offsetsx.Size(); i++)    
-   { 
-      block_offsetsx[i] = block_offsetsuml[i] ; 
-   }
-
-  
-   ml = problem->Getml();
-  
-   lk.SetSize(dimC);  lk  = 0.0;
-   zlk.SetSize(dimM); zlk = 0.0;
-
-   linSolver = 0;
-   linSolveTol = 1.e-8;
-   MyRank = Mpi::WorldRank();
-   iAmRoot = MyRank == 0 ? true : false;
-}
-
-double ParInteriorPointSolver::MaxStepSize(Vector &x, Vector &xl, Vector &xhat, double tau)
-{
-   double alphaMaxloc = 1.0;
-   double alphaTmp;
-   for(int i = 0; i < x.Size(); i++)
-   {   
-      if( xhat(i) < 0. )
-      {
-         alphaTmp = -1. * tau * (x(i) - xl(i)) / xhat(i);
-         alphaMaxloc = min(alphaMaxloc, alphaTmp);
-      } 
-   }
-
-   // alphaMaxloc is the local maximum step size which is
-   // distinct on each MPI process. Need to compute
-   // the global maximum step size 
-   double alphaMaxglb;
-   MPI_Allreduce(&alphaMaxloc, &alphaMaxglb, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
-   return alphaMaxglb;
-}
-
-double ParInteriorPointSolver::MaxStepSize(Vector &x, Vector &xhat, double tau)
-{
-   Vector zero(x.Size()); zero = 0.0;
-   return MaxStepSize(x, zero, xhat, tau);
-}
-
-
-void ParInteriorPointSolver::Mult(const Vector &x0, Vector &xf)
-{
-   BlockVector x0block(block_offsetsx); x0block = 0.0;
-   x0block.GetBlock(0).Set(1.0, x0);
-   
-   
-   ParOptProblem * OptProblem = dynamic_cast<ParOptProblem *>(problem);
-   if (dimM > 0)
-   {
-      if (true)//OptProblem == nullptr)
-      { 
-	 // fixed initialization     
-         x0block.GetBlock(1) = 1.0;
-	 x0block.GetBlock(1).Add(1.0, ml);
-      }
-      else
-      {
-         // use g(d) - s = 0
-         x0block.GetBlock(1) = 0.0;
-         x0block.GetBlock(1).Add(1.0, ml);
-         Vector c0(dimC); c0 = 0.0;	 
-	 problem->c(x0block, c0);
-	 Vector dm(dimM); dm = 0.0;
-	 for (int i = 0; i < dimM; i++)
-	 {
-	    dm(i) = max(1.e0 - x0block(dimU + i), c0(i));
-	 }
-	 x0block.GetBlock(1).Add(1.0, dm);
-      }
-   }  
-   
-   BlockVector xfblock(block_offsetsx); xfblock = 0.0;
-   Mult(x0block, xfblock);
-   xf.Set(1.0, xfblock.GetBlock(0));
-}
-
-
-void ParInteriorPointSolver::Mult(const BlockVector &x0, BlockVector &xf)
-{
-   converged = false;
-   
-   BlockVector xk(block_offsetsx), xhat(block_offsetsx); xk = 0; xhat = 0.0;
-   BlockVector Xk(block_offsetsumlz), Xhat(block_offsetsumlz); Xk = 0.0; Xhat = 0.0;
-   BlockVector Xhatuml(block_offsetsuml); Xhatuml = 0.0;
-   Vector zlhat(dimM); zlhat = 0.0;
-
-   xk.GetBlock(0).Set(1.0, x0.GetBlock(0));
-   xk.GetBlock(1).Set(1.0, x0.GetBlock(1));
-   
-   // running estimate of the final values of the Lagrange multipliers
-   lk = 0.0;
-   zlk = 0.0;
-
-   for(int i = 0; i < dimM; i++)
-   {
-      zlk(i) = 1.e1 * mu_k / (xk(i+dimU) - ml(i));
-   }
-
-   Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
-   Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
-   Xk.GetBlock(2).Set(1.0, lk);
-   Xk.GetBlock(3).Set(1.0, zlk);
-
-  /* set theta0 = theta(x0)
-   *     thetaMin
-   *     thetaMax
-   * when theta(xk) < thetaMin and the switching condition holds
-   * then we ask for the Armijo sufficient decrease of the barrier
-   * objective to be satisfied, in order to accept the trial step length alphakl
-   * 
-   * thetaMax controls how the filter is initialized for each log-barrier subproblem
-   * F0 = {(th, phi) s.t. th > thetaMax}
-   * that is the filter does not allow for iterates where the constraint violation
-   * is larger than that of thetaMax
-   */
-   double theta0 = theta(xk);
-   thetaMin = 1.e-4 * max(1.0, theta0);
-   thetaMax = 1.e8  * thetaMin; // 1.e4 * max(1.0, theta0)
-
-   double Eeval, maxBarrierSolves, Eevalmu0;
-   bool printOptimalityError; // control optimality error print to console for log-barrier subproblems
-   
-   maxBarrierSolves = 10;
-   bool smallStep;
-   int numSmallSteps = 0;
-   for(jOpt = 0; jOpt < max_iter; jOpt++)
-   {
-      if(iAmRoot)
-      {
-         cout << "interior-point solve step " << jOpt << endl;
-      }
-      // A-2. Check convergence of overall optimization problem
-      printOptimalityError = false;
-      Eevalmu0 = E(xk, lk, zlk, printOptimalityError);
-      if(Eevalmu0 < OptTol)
-      {
-         converged = true;
-	 if(iAmRoot)
-         {
-            cout << "solved optimization problem :)\n";
-         }
-         break;
-      }
-      
-      if(jOpt > 0) { maxBarrierSolves = 1; }
-      
-      for(int i = 0; i < maxBarrierSolves; i++)
-      {
-         // A-3. Check convergence of the barrier subproblem
-         printOptimalityError = true;
-         Eeval = E(xk, lk, zlk, mu_k, printOptimalityError);
-         if(iAmRoot)
-         {
-            cout << "E = " << Eeval << endl;
-         }
-         if(Eeval < kEps * mu_k)
-         {
-            if(iAmRoot)
-            {
-               cout << "solved barrier subproblem :), for mu = " << mu_k << endl;
-            }
-            // A-3.1. Recompute the barrier parameter
-            mu_k  = max(OptTol / 10., min(kMu * mu_k, pow(mu_k, thetaMu)));
-            // A-3.2. Re-initialize the filter
-            F1.DeleteAll();
-            F2.DeleteAll();
-         }
-         else
-         {
-            break;
-         }
-      }
-    
-      // A-4. Compute the search direction
-      // solve for (uhat, mhat, lhat)
-      if(iAmRoot)
-      {
-         cout << "\n** A-4. IP-Newton solve **\n";
-      }
-      zlhat = 0.0; Xhatuml = 0.0;
-      IPNewtonSolve(xk, lk, zlk, zlhat, Xhatuml, smallStep, mu_k, false); 
-      if (smallStep && numSmallSteps < 5)
-      {
-         numSmallSteps += 1;
-	 xk.GetBlock(0).Add(1.0, Xhatuml.GetBlock(0));
-         xk.GetBlock(1).Add(1.0, Xhatuml.GetBlock(1));
-         lk.Add(1.0,   Xhatuml.GetBlock(2));
-         zlk.Add(1.0, zlhat);
-         continue;	 
-      }
-      numSmallSteps = 0;
-
-      // assign data stack, X = (u, m, l, zl)
-      Xk = 0.0;
-      Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
-      Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
-      Xk.GetBlock(2).Set(1.0, lk);
-      Xk.GetBlock(3).Set(1.0, zlk);
-
-      // assign data stack, Xhat = (uhat, mhat, lhat, zlhat)
-      Xhat = 0.0;
-      for(int i = 0; i < 3; i++)
-      {
-         Xhat.GetBlock(i).Set(1.0, Xhatuml.GetBlock(i));
-      }
-      Xhat.GetBlock(3).Set(1.0, zlhat);
-
-      // A-5. Backtracking line search.
-      if(iAmRoot)
-      {
-         cout << "\n** A-5. Linesearch **\n";
-         cout << "mu = " << mu_k << endl;
-      }
-      lineSearch(Xk, Xhat, mu_k);
-
-      if(lineSearchSuccess)
-      {
-         if(iAmRoot)
-         {
-            cout << "lineSearch successful :)\n";
-         }
-         if(!switchCondition || !sufficientDecrease)
-         {
-            F1.Append( (1. - gTheta) * thx0);
-            F2.Append( phx0 - gPhi * thx0);
-         }
-         // ----- A-6: Accept the trial point
-         xk.GetBlock(0).Add(alpha, Xhat.GetBlock(0));
-         xk.GetBlock(1).Add(alpha, Xhat.GetBlock(1));
-         lk.Add(alpha,   Xhat.GetBlock(2));
-         zlk.Add(alphaz, Xhat.GetBlock(3));
-         projectZ(xk, zlk, mu_k);
-      }
-      else
-      {
-         if(iAmRoot)
-         {
-            cout << "lineSearch not successful :(\n";
-            cout << "attempting feasibility restoration with theta = " << thx0 << endl;
-         }
-	 FeasibilityRestoration(xk, lk, zlk, Xk, mu_k);
-         xk.GetBlock(0).Set(1.0, Xk.GetBlock(0));
-         xk.GetBlock(1).Set(1.0, Xk.GetBlock(1));
-         lk.Set(1.0,   Xk.GetBlock(2));
-         zlk.Set(1.0,  Xk.GetBlock(3));
-      }
-      if(jOpt + 1 == max_iter && iAmRoot) 
-      {  
-         cout << "maximum optimization iterations :(\n";
-      }
-   }
-   // done with optimization routine, just reassign data to xf reference so
-   // that the application code has access to the optimal point
-   xf = 0.0;
-   xf.GetBlock(0).Set(1.0, xk.GetBlock(0));
-   xf.GetBlock(1).Set(1.0, xk.GetBlock(1));
-}
-
-void ParInteriorPointSolver::FormIPNewtonMat(BlockVector & x, Vector & l, Vector &zl, BlockOperator &Ak)
-{
-   // WARNING: Huu, Hum, Hmu, Hmm should all be Hessian terms of the Lagrangian, currently we 
-   //          them by Hessian terms of the objective function and neglect the Hessian of l^T c
-
-   Huu = problem->Duuf(x); 
-   Hum = problem->Dumf(x);
-   Hmu = problem->Dmuf(x);
-   Hmm = problem->Dmmf(x);
-
-   Vector DiagLogBar(dimM); DiagLogBar = 0.0;
-   for(int ii = 0; ii < dimM; ii++)
-   {
-      DiagLogBar(ii) = zl(ii) / (x(ii+dimU) - ml(ii));
-   }
-   
-   if (saveIterates)
-   {
-      std::ofstream diagStream;
-      char diagString[100];
-      snprintf(diagString, 100, "logBarrierHessiandata/D%d.dat", jOpt);
-      diagStream.open(diagString, ios::out | ios::trunc);
-      for(int ii = 0; ii < dimM; ii++)
-      {
-         diagStream << setprecision(30) << DiagLogBar(ii) << endl;
-      }
-      diagStream.close();
-
-      std::ofstream sStream;
-      char sString[100];
-      snprintf(sString, 100, "logBarrierHessiandata/s%d.dat", jOpt);
-      sStream.open(sString, ios::out | ios::trunc);
-      for(int ii = 0; ii < dimM; ii++)
-      {
-         sStream << setprecision(30) << x(ii+dimU) << endl;
-      }
-      sStream.close();
-      
-      std::ofstream lStream;
-      char lString[100];
-      snprintf(lString, 100, "logBarrierHessiandata/l%d.dat", jOpt);
-      lStream.open(lString, ios::out | ios::trunc);
-      for(int ii = 0; ii < dimM; ii++)
-      {
-         lStream << setprecision(30) << l(ii) << endl;
-      }
-      lStream.close();
-      
-      std::ofstream zlStream;
-      char zlString[100];
-      snprintf(zlString, 100, "logBarrierHessiandata/zl%d.dat", jOpt);
-      zlStream.open(zlString, ios::out | ios::trunc);
-      for(int ii = 0; ii < dimM; ii++)
-      {
-         zlStream << setprecision(30) << zl(ii) << endl;
-      }
-      zlStream.close();
-      
-      std::ofstream dStream;
-      char dString[100];
-      snprintf(dString, 100, "logBarrierHessiandata/d%d.dat", jOpt);
-      dStream.open(dString, ios::out | ios::trunc);
-      for(int ii = 0; ii < dimU; ii++)
-      {
-         dStream << setprecision(30) << x(ii) << endl;
-      }
-      dStream.close();
-   } 
-
-   D = GenerateHypreParMatrixFromDiagonal(problem->GetDofOffsetsM(), DiagLogBar);
-  
-   if(Hmm != nullptr)
-   {
-      Wmm = Hmm;
-      Wmm->Add(1.0, *D);
-   }
-   else
-   {
-      Wmm = D;
-   }
-
-   Ju = problem->Duc(x); JuT = Ju->Transpose();
-   Jm = problem->Dmc(x); JmT = Jm->Transpose();
-   //         IP-Newton system matrix
-   //    Ak = [[H_(u,u)  H_(u,m)   J_u^T]
-   //          [H_(m,u)  W_(m,m)   J_m^T]
-   //          [ J_u      J_m       0  ]]
-
-   Ak.SetBlock(0, 0, Huu);                         Ak.SetBlock(0, 2, JuT);
-                           Ak.SetBlock(1, 1, Wmm); Ak.SetBlock(1, 2, JmT);
-   Ak.SetBlock(2, 0,  Ju); Ak.SetBlock(2, 1,  Jm);
-
-   if(Hum != nullptr) { Ak.SetBlock(0, 1, Hum); Ak.SetBlock(1, 0, Hmu); }
-}
-
-// perturbed KKT system solve
-// determine the search direction
-void ParInteriorPointSolver::IPNewtonSolve(BlockVector &x, Vector &l, Vector &zl, Vector &zlhat, BlockVector &Xhat, bool & smallStep, double mu, bool socSolve)
-{
-   // solve A x = b, where A is the IP-Newton matrix
-   BlockOperator A(block_offsetsuml, block_offsetsuml); BlockVector b(block_offsetsuml); b = 0.0;
-   FormIPNewtonMat(x, l, zl, A);
-
-   //       [grad_u phi + Ju^T l]
-   // b = - [grad_m phi + Jm^T l]
-   //       [          c        ]
-   BlockVector gradphi(block_offsetsx); gradphi = 0.0;
-   BlockVector JTl(block_offsetsx); JTl = 0.0;
-   Dxphi(x, mu, gradphi);
-   
-   (A.GetBlock(0,2)).Mult(l, JTl.GetBlock(0));
-   (A.GetBlock(1,2)).Mult(l, JTl.GetBlock(1));
-
-   for(int ii = 0; ii < 2; ii++)
-   {
-      b.GetBlock(ii).Set(1.0, gradphi.GetBlock(ii));
-      b.GetBlock(ii).Add(1.0, JTl.GetBlock(ii));
-   }
-   if(!socSolve) 
-   {
-      problem->c(x, b.GetBlock(2));
-   }
-   else
-   {
-      b.GetBlock(2).Set(1.0, ckSoc);
-   }
-   b *= -1.0; 
-   Xhat = 0.0;
-
-
-   // Direct solver (default)
-   if(linSolver == 0)
-   {
-      Array2D<HypreParMatrix *> ABlockMatrix(3,3);
-      for(int ii = 0; ii < 3; ii++)
-      {
-      for(int jj = 0; jj < 3; jj++)
-      {
-         if(!A.IsZeroBlock(ii, jj))
-         {
-            ABlockMatrix(ii, jj) = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(ii, jj))));
-	 }
-	 else
-	 {
-	    ABlockMatrix(ii, jj) = nullptr;
-	 }
-      }
-      }
-      
-      HypreParMatrix * Ah = HypreParMatrixFromBlocks(ABlockMatrix);   
-      
-      /* direct solve of the 3x3 IP-Newton linear system */
-      #ifdef MFEM_USE_MUMPS
-        MUMPSSolver ASolver;
-        ASolver.SetPrintLevel(0);
-        ASolver.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_INDEFINITE);
-        ASolver.SetOperator(*Ah);
-        ASolver.Mult(b, Xhat);
-      #else 
-        #ifdef MFEM_USE_MKL_CPARDISO
-          CPardisoSolver ASolver(MPI_COMM_WORLD);
-          ASolver.SetOperator(*Ah);
-          ASolver.Mult(b, Xhat);
-	#else
-	  MFEM_VERIFY(false, "linSolver 0 will not work unless compiled with MUMPS or MKL");
-        #endif
-      #endif
-
-      delete Ah;
-   }
-   else if(linSolver == 1 || linSolver == 2)
-   {
-     
-      // assuming Jm = -I and Hum = 0, Hmu = 0, Hmm = 0
-      ParOptProblem * tempProblem = dynamic_cast<ParOptProblem *>(problem);
-      MFEM_VERIFY(tempProblem != nullptr, "linSolver option 1 and 2 are only applicable to ParOptProblem's");      
-      
-
-
-      // form A = Huu + Ju^T D Ju, Wmm = D for contact
-      HypreParMatrix * Huuloc = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(0, 0))));
-      HypreParMatrix * Wmmloc = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(1, 1))));
-      HypreParMatrix * Juloc  = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(2, 0))));
-      HypreParMatrix * JuTloc = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(0, 2))));
-      
-      
-      HypreParMatrix *JuTDJu   = RAP(Wmmloc, Juloc);      // Ju^T D Ju
-      HypreParMatrix *Areduced = ParAdd(Huuloc, JuTDJu);  // Huu + Ju^T D Ju
-      /* prepare the reduced rhs 
-       * breduced = bu + Ju^T (bm + Wmm bl) */
-      Vector breduced(dimU); breduced = 0.0;
-      Vector tempVec(dimM); tempVec = 0.0;
-      Wmmloc->Mult(b.GetBlock(2), tempVec);
-      tempVec.Add(1.0, b.GetBlock(1));
-      JuTloc->Mult(tempVec, breduced);
-      breduced.Add(1.0, b.GetBlock(0));
-      
-      if(linSolver == 1)
-      {
-         // setup the solver for the reduced linear system
-         #ifdef MFEM_USE_MUMPS
-	   MUMPSSolver AreducedSolver;   
-           AreducedSolver.SetPrintLevel(0);
-           AreducedSolver.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_INDEFINITE);
-	   AreducedSolver.SetOperator(*Areduced);
-	   AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-         #else 
-           #ifdef MFEM_USE_MKL_CPARDISO
-	     CPardisoSolver AreducedSolver(MPI_COMM_WORLD);
-	     AreducedSolver.SetOperator(*Areduced);
-	     AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-           #else
-	     MFEM_VERIFY(false, "linSolver 1 will not work unless compiled with MUMPS or MKL");
-           #endif
-         #endif
-      }
-      else
-      {
-         HyprePCG AreducedSolver(MPI_COMM_WORLD);
-	 AreducedSolver.SetOperator(*Areduced);
-	 HypreBoomerAMG AreducedPrec;
-         AreducedSolver.SetTol(linSolveTol);
-         AreducedSolver.SetMaxIter(500);
-         AreducedSolver.SetPreconditioner(AreducedPrec);
-         AreducedSolver.SetResidualConvergenceOptions(); // convergence criteria based on residual norm
-	 AreducedSolver.SetPrintLevel(2);
-         AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-      }
-
-      // now propagate solved uhat to obtain mhat and lhat
-      // xm = Ju xu - bl
-      Juloc->Mult(Xhat.GetBlock(0), Xhat.GetBlock(1));
-      Xhat.GetBlock(1).Add(-1.0, b.GetBlock(2));
-
-      // xl = Wmm xm - bm
-      Wmmloc->Mult(Xhat.GetBlock(1), Xhat.GetBlock(2));
-      Xhat.GetBlock(2).Add(-1.0, b.GetBlock(1));
-
-      delete JuTDJu;
-      delete Areduced;
-   }
-
-   /* backsolve to determine zlhat */
-   for(int ii = 0; ii < dimM; ii++)
-   {
-      zlhat(ii) = -1.*(zl(ii) + (zl(ii) * Xhat(ii + dimU) - mu) / (x(ii + dimU) - ml(ii)) );
-   }
-   
-   Vector smallStepCheckVec(dimU + dimM); smallStepCheckVec = 0.0;
-   for(int ii = 0; ii < dimU + dimM; ii++)
-   {
-      smallStepCheckVec(ii) = abs(Xhat(ii)) / (1. + abs(x(ii)));
-   }
-   double smallStepCheckVal = GlobalLpNorm(infinity(), smallStepCheckVec.Normlinf(), MPI_COMM_WORLD);
-   smallStep = (smallStepCheckVal < 1.e-15) ? true : false;
-   if (smallStep && iAmRoot)
-   {
-      cout << "SMALL STEP\n";
-   }
-   
-   
-   // free memory
-   delete D;
-   delete JuT;
-   delete JmT;
-   if(Hmm != nullptr)
-   {
-      delete Wmm;
-   }
-}
-
-// here Xhat, X will be BlockVectors w.r.t. the 4 partitioning X = (u, m, l, zl)
-
-void ParInteriorPointSolver::lineSearch(BlockVector& X0, BlockVector& Xhat, double mu)
-{
-   double tau  = tauMin;
-   Vector u0   = X0.GetBlock(0);
-   Vector m0   = X0.GetBlock(1);
-   Vector l0   = X0.GetBlock(2);
-   Vector z0   = X0.GetBlock(3);
-   Vector uhat = Xhat.GetBlock(0);
-   Vector mhat = Xhat.GetBlock(1);
-   Vector lhat = Xhat.GetBlock(2);
-   Vector zhat = Xhat.GetBlock(3);
-   double alphaMax  = MaxStepSize(m0, ml, mhat, tau);
-   double alphaMaxz = MaxStepSize(z0, zhat, tau);
-   alphaz = alphaMaxz;
-
-   BlockVector x0(block_offsetsx); x0 = 0.0;
-   x0.GetBlock(0).Set(1.0, u0);
-   x0.GetBlock(1).Set(1.0, m0);
-   
-   BlockVector xhat(block_offsetsx); xhat = 0.0;
-   xhat.GetBlock(0).Set(1.0, uhat);
-   xhat.GetBlock(1).Set(1.0, mhat);
-   
-   BlockVector xtrial(block_offsetsx); xtrial = 0.0;
-   BlockVector Dxphi0(block_offsetsx); Dxphi0 = 0.0;
-   int maxBacktrack = 20;
-   alpha = alphaMax;
-
-   Vector ck0(dimC); ck0 = 0.0;
-   Vector zhatsoc(dimM); zhatsoc = 0.0;
-   BlockVector Xhatumlsoc(block_offsetsuml); Xhatumlsoc = 0.0;
-   BlockVector xhatsoc(block_offsetsx); xhatsoc = 0.0;
-   Vector uhatsoc(dimU); uhatsoc = 0.0;
-   Vector mhatsoc(dimM); mhatsoc = 0.0;
-
-   Dxphi(x0, mu, Dxphi0);
-
-   Dxphi0_xhat = InnerProduct(MPI_COMM_WORLD, Dxphi0, xhat);
-   descentDirection = Dxphi0_xhat < 0. ? true : false;
-   if(descentDirection)
-   {
-      if (iAmRoot)
-      {
-         cout << "is a descent direction for the log-barrier objective\n";
-      }
-   }
-   else
-   {
-      cout << "is not a descent direction for the log-barrier objective\n";
-   }
-   thx0 = theta(x0);
-   phx0 = phi(x0, mu);
-
-   lineSearchSuccess = false;
-   for(int i = 0; i < maxBacktrack; i++)
-   {
-      if (iAmRoot)
-      {
-         cout << "\n--------- alpha = " << alpha << " ---------\n";
-      }
-      // ----- A-5.2. Compute trial point: xtrial = x0 + alpha_i xhat
-      xtrial.Set(1.0, x0);
-      xtrial.Add(alpha, xhat);
-
-      // ------ A-5.3. if not in filter region go to A.5.4 otherwise go to A-5.5.
-      thxtrial = theta(xtrial);
-      phxtrial = phi(xtrial, mu);
-      filterCheck(thxtrial, phxtrial);    
-      if (!inFilterRegion)
-      {
-         if (iAmRoot)
-         {
-            cout << "not in filter region :)\n";
-         }
-         // ------ A.5.4: Check sufficient decrease
-         if(!descentDirection)
-         {
-            switchCondition = false;
-         }
-         else
-         {
-            switchCondition = (alpha * pow(abs(Dxphi0_xhat), sPhi) > delta * pow(thx0, sTheta)) ? true : false;
-         }
-         if (iAmRoot)
-         {
-            cout << "theta(x0) = "     << thx0     << ", thetaMin = "                  << thetaMin             << endl;
-            cout << "theta(xtrial) = " << thxtrial << ", (1-gTheta) *theta(x0) = "     << (1. - gTheta) * thx0 << endl;
-            cout << "phi(xtrial) = "   << phxtrial << ", phi(x0) - gPhi *theta(x0) = " << phx0 - gPhi * thx0   << endl;
-         }      
-         // Case I      
-         if(thx0 <= thetaMin && switchCondition)
-         {
-            sufficientDecrease = (phxtrial <= phx0 + eta * alpha * Dxphi0_xhat) ? true : false;
-            if(sufficientDecrease)
-            {
-               if(iAmRoot) { cout << "Line search successful: sufficient decrease in log-barrier objective.\n"; } 
-               // accept the trial step
-               lineSearchSuccess = true;
-               break;
-            }
-         }
-         else
-         {
-            if(thxtrial <= (1. - gTheta) * thx0 || phxtrial <= phx0 - gPhi * thx0)
-            {
-               if(iAmRoot) { cout << "Line search successful: infeasibility or log-barrier objective decreased.\n"; } 
-               // accept the trial step
-               lineSearchSuccess = true;
-               break;
-            }
-         }
-      }
-      else
-      {
-         if (iAmRoot)
-         {
-            cout << "in filter region :(\n"; 
-         }
-      }
-      alpha *= 0.5;
-
-   } 
-}
-
-
-void ParInteriorPointSolver::projectZ(const Vector &x, Vector &z, double mu)
-{
-   double zi;
-   double mudivmml;
-   for(int i = 0; i < dimM; i++)
-   {
-      zi = z(i);
-      mudivmml = mu / (x(i + dimU) - ml(i));
-      z(i) = max(min(zi, kSig * mudivmml), mudivmml / kSig);
-   }
-}
-
-void ParInteriorPointSolver::filterCheck(double th, double ph)
-{
-   inFilterRegion = false;
-   if(th > thetaMax)
-   {
-      inFilterRegion = true;
-   }
-   else
-   {
-      for(int i = 0; i < F1.Size(); i++)
-      {
-         if(th >= F1[i] && ph >= F2[i])
-         {
-            inFilterRegion = true;
-            break;
-         }
-      }
-   }
-}
-
-void ParInteriorPointSolver::FeasibilityRestoration(const BlockVector & x, const Vector &l, const Vector &zl, BlockVector &X, double mu)
-{
-   ParOptProblem * OptProblem = dynamic_cast<ParOptProblem *>(problem);
-   X.GetBlock(0).Set(1.0, x.GetBlock(0));
-   X.GetBlock(1).Set(1.0, x.GetBlock(1));
-   X.GetBlock(2).Set(1.0, l);
-   X.GetBlock(3).Set(1.0, zl);
-
-   if (OptProblem != nullptr)
-   {
-      Vector g(dimC); g = 0.0;
-      OptProblem->g(x.GetBlock(0), g);
-      Ju = OptProblem->Ddg(x.GetBlock(0));
-      Ju->DropSmallEntries(1.e-16);
-      SparseMatrix JuMerged;
-      Ju->MergeDiagAndOffd(JuMerged);
-      
-      int num_loc_modified_constraints = 0;
-      int num_glb_modified_constraints = 0;
-      for (int i = 0; i < dimC; i++)
-      {
-         if(JuMerged.RowIsEmpty(i))
-	 {
-	    if (g(i) < 1.e-15)
-	    {
-	       if (iAmRoot)
-	       {
-	           cout << "WARNING: LICQ violation detected (g_i and grad(g_i) both zero)\n";
-	       }
-	       continue;
-	    }
-	    X(dimU + i) = g(i); // s_i = \gamma
-            X(dimU + dimM + i) = -1. * mu / g(i); // l_i = -mu/\gamma	 
-	    X(dimU + dimM + dimC + i) = mu / g(i);
-	    num_loc_modified_constraints += 1;
-	 }
-      }
-      MPI_Allreduce(&num_loc_modified_constraints, &num_glb_modified_constraints, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-      if (num_glb_modified_constraints == 0)
-      {
-         if(iAmRoot)
-         {
-            cout << "trying feasibility restoration with no null rows in Jacobian\n";
-            cout << "exiting\n";
-         }
-         exit(1);
-      }   
-   }
-}
-
-
-
-double ParInteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, double mu, bool printEeval)
-{
-   double E1, E2, E3;
-   double sc, sd;
-   BlockVector gradL(block_offsetsx); gradL = 0.0; // stationarity grad L = grad f + J^T l - z
-   Vector cx(dimC); cx = 0.0;     // feasibility c = c(x)
-   Vector comp(dimM); comp = 0.0; // complementarity M Z - mu 1
-
-   DxL(x, l, zl, gradL);
-   E1 = GlobalLpNorm(infinity(), gradL.Normlinf(), MPI_COMM_WORLD);
-
-   problem->c(x, cx);
-   E2 = GlobalLpNorm(infinity(), cx.Normlinf(), MPI_COMM_WORLD); 
-
-   for(int ii = 0; ii < dimM; ii++) 
-   { 
-      comp(ii) = x(dimU + ii) * zl(ii) - mu;
-   }
-   E3 = GlobalLpNorm(infinity(), comp.Normlinf(), MPI_COMM_WORLD); 
-   double ll1, zl1;
-
-   zl1 = GlobalLpNorm(1, zl.Norml1(), MPI_COMM_WORLD); 
-   ll1 = GlobalLpNorm(1, l.Norml1(), MPI_COMM_WORLD);
-   sc = max(sMax, zl1 / (double(dimMGlb)) ) / sMax;
-   sd = max(sMax, (ll1 + zl1) / (double(dimCGlb + dimMGlb))) / sMax;
-   if(iAmRoot)
-   {
-      cout << "evaluating optimality error for mu = " << mu << endl;
-      cout << "stationarity measure = "    << E1 / sd << endl;
-      cout << "feasibility measure  = "    << E2      << endl;
-      cout << "complimentarity measure = " << E3 / sc << endl;
-   }
-   return max(max(E1 / sd, E2), E3 / sc);
-}
-
-double ParInteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, bool printEeval)
-{
-  return E(x, l, zl, 0.0, printEeval);
-}
-
-double ParInteriorPointSolver::theta(const BlockVector &x)
-{
-  Vector cx(dimC); cx = 0.0;
-  problem->c(x, cx);
-  return GlobalLpNorm(2, cx.Norml2(), MPI_COMM_WORLD);
-}
-
-// log-barrier objective
-double ParInteriorPointSolver::phi(const BlockVector &x, double mu)
-{
-   double fx = problem->CalcObjective(x); 
-   double logBarrierLoc = 0.0;
-   for(int i = 0; i < dimM; i++) 
-   { 
-     logBarrierLoc += log(x(dimU+i) - ml(i));
-   }
-   double logBarrierGlb;
-   MPI_Allreduce(&logBarrierLoc, &logBarrierGlb, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-   return fx - mu * logBarrierGlb;
-}
-
-// gradient of log-barrier objective with respect to x = (u, m)
-void ParInteriorPointSolver::Dxphi(const BlockVector &x, double mu, BlockVector &y)
-{
-   problem->CalcObjectiveGrad(x, y);
-   
-   for(int i = 0; i < dimM; i++) 
-   { 
-      y(dimU + i) -= mu / (x(dimU + i));
-   } 
-}
-
-// Lagrangian function evaluation
-// L(x, l, zl) = f(x) + l^T c(x) - zl^T m
-double ParInteriorPointSolver::L(const BlockVector &x, const Vector &l, const Vector &zl)
-{
-   double fx = problem->CalcObjective(x);
-   Vector cx(dimC); problem->c(x, cx);
-   return (fx + InnerProduct(MPI_COMM_WORLD, cx, l) - InnerProduct(MPI_COMM_WORLD, x.GetBlock(1), zl));
-}
-
-void ParInteriorPointSolver::DxL(const BlockVector &x, const Vector &l, const Vector &zl, BlockVector &y)
-{
-   // evaluate the gradient of the objective with respect to the primal variables x = (u, m)
-   BlockVector gradxf(block_offsetsx); gradxf = 0.0;
-   problem->CalcObjectiveGrad(x, gradxf);
-   
-   HypreParMatrix *Jacu, *Jacm;
-   Jacu = problem->Duc(x); 
-   Jacm = problem->Dmc(x);
-   Jacu->MultTranspose(l, y.GetBlock(0));
-   Jacm->MultTranspose(l, y.GetBlock(1));   
-   
-   y.Add(1.0, gradxf);
-   (y.GetBlock(1)).Add(-1.0, zl);
-}
-
-bool ParInteriorPointSolver::GetConverged() const
-{
-   return converged;
-}
-
-void ParInteriorPointSolver::SetTol(double Tol)
-{
-   OptTol = Tol;
-}
-
-void ParInteriorPointSolver::SetMaxIter(int max_it)
-{
-   max_iter = max_it;
-}
-
-void ParInteriorPointSolver::SetBarrierParameter(double mu_0)
-{
-   mu_k = mu_0;
-}
-
-void ParInteriorPointSolver::SaveIterates(bool save)
-{
-   saveIterates = save;
-}
-
-
-void ParInteriorPointSolver::SetLinearSolver(int LinSolver)
-{
-   linSolver = LinSolver;
-}
-
-void ParInteriorPointSolver::SetLinearSolveTol(double Tol)
-{
-  linSolveTol = Tol;
-}
-
-void ParInteriorPointSolver::GetLagrangeMultiplier(Vector & y)
-{
-  y.SetSize(dimM); y = 0.;
-  y.Set(1.0, zlk);
-}
-
-
-
-ParInteriorPointSolver::~ParInteriorPointSolver() 
-{
-   F1.DeleteAll();
-   F2.DeleteAll();
-   block_offsetsx.DeleteAll();
-   block_offsetsumlz.DeleteAll();
-   block_offsetsuml.DeleteAll();
-   ml.SetSize(0);
-}

miniapps/interiorpointsolver/IPsolver.hpp

@@ -1,83 +0,0 @@
-#include "mfem.hpp"
-#include "Problem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-#ifndef PARIPSOLVER 
-#define PARIPSOLVER
-
-class ParInteriorPointSolver
-{
-protected:
-    ParGeneralOptProblem* problem;
-    double OptTol;
-    int  max_iter;
-    double mu_k; // \mu_k
-    Vector lk, zlk;
-
-    double sMax, kSig, tauMin, eta, thetaMin, delta, sTheta, sPhi, kMu, thetaMu;
-    double thetaMax, kSoc, gTheta, gPhi, kEps;
-	
-    // filter
-    Array<double> F1, F2;
-	
-    // quantities computed in lineSearch
-    double alpha, alphaz;
-    double thx0, thxtrial;
-    double phx0, phxtrial;
-    bool descentDirection, switchCondition, sufficientDecrease, lineSearchSuccess, inFilterRegion;
-    double Dxphi0_xhat;
-
-    int dimU, dimM, dimC;
-    int dimUGlb, dimMGlb, dimCGlb;
-    Array<int> block_offsetsumlz, block_offsetsuml, block_offsetsx;
-    Vector ml;
-
-    Vector ckSoc;
-    HypreParMatrix * Huu, * Hum, * Hmu, * Hmm, * Wmm, *D, * Ju, * Jm, * JuT, * JmT;
-    
-    int jOpt;
-    bool converged;
-    
-    int MyRank;
-    bool iAmRoot;
-
-    bool saveLogBarrierIterates;
-    bool saveIterates;
-    int linSolver;
-    double linSolveTol;
-public:
-    ParInteriorPointSolver(ParGeneralOptProblem*);
-    double MaxStepSize(Vector& , Vector& , Vector& , double);
-    double MaxStepSize(Vector& , Vector& , double);
-    void Mult(const BlockVector& , BlockVector&);
-    void Mult(const Vector&, Vector &); 
-    void GetLagrangeMultiplier(Vector &);
-    void FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
-    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , bool &, double, bool);
-    void lineSearch(BlockVector& , BlockVector& , double);
-    void projectZ(const Vector & , Vector &, double);
-    void filterCheck(double, double);
-    double E(const BlockVector &, const Vector &, const Vector &, double, bool);
-    double E(const BlockVector &, const Vector &, const Vector &, bool);
-    bool GetConverged() const;
-    double theta(const BlockVector &);
-    double phi(const BlockVector &, double);
-    void Dxphi(const BlockVector &, double, BlockVector &);
-    double L(const BlockVector &, const Vector &, const Vector &);
-    void DxL(const BlockVector &, const Vector &, const Vector &, BlockVector &);
-    void SetTol(double);
-    void SetMaxIter(int);
-    void SetBarrierParameter(double);    
-    void SaveIterates(bool);
-    void SetLinearSolver(int);
-    void SetLinearSolveTol(double);
-    void FeasibilityRestoration(const BlockVector &, const Vector &, const Vector &, BlockVector &, double); 
-    virtual ~ParInteriorPointSolver();
-};
-
-#endif

miniapps/interiorpointsolver/Problem.cpp

@@ -1,404 +0,0 @@
-#include "mfem.hpp"
-#include "Problem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-
-ParGeneralOptProblem::ParGeneralOptProblem() : block_offsetsx(3) {}
-
-void ParGeneralOptProblem::Init(HYPRE_BigInt * dofOffsetsU_, HYPRE_BigInt * dofOffsetsM_)
-{
-  dofOffsetsU = new HYPRE_BigInt[2];
-  dofOffsetsM = new HYPRE_BigInt[2];
-  for(int i = 0; i < 2; i++)
-  {
-    dofOffsetsU[i] = dofOffsetsU_[i];
-    dofOffsetsM[i] = dofOffsetsM_[i];
-  }
-  dimU = dofOffsetsU[1] - dofOffsetsU[0];
-  dimM = dofOffsetsM[1] - dofOffsetsM[0];
-  dimC = dimM;
-  
-  block_offsetsx[0] = 0;
-  block_offsetsx[1] = dimU;
-  block_offsetsx[2] = dimM;
-  block_offsetsx.PartialSum();
-  
-  MPI_Allreduce(&dimU, &dimUglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-  MPI_Allreduce(&dimM, &dimMglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-}
-
-
-void ParGeneralOptProblem::CalcObjectiveGrad(const BlockVector &x, BlockVector &y) const
-{
-   Duf(x, y.GetBlock(0));
-   Dmf(x, y.GetBlock(1));
-}
-
-ParGeneralOptProblem::~ParGeneralOptProblem()
-{
-   block_offsetsx.DeleteAll();
-}
-
-
-// min E(d) s.t. g(d) >= 0
-// min_(d,s) E(d) s.t. c(d,s) := g(d) - s = 0, s >= 0
-ParOptProblem::ParOptProblem() : ParGeneralOptProblem()
-{
-}
-
-void ParOptProblem::Init(HYPRE_BigInt * dofOffsetsU_, HYPRE_BigInt * dofOffsetsM_)
-{
-  dofOffsetsU = new HYPRE_BigInt[2];
-  dofOffsetsM = new HYPRE_BigInt[2];
-  for(int i = 0; i < 2; i++)
-  {
-    dofOffsetsU[i] = dofOffsetsU_[i];
-    dofOffsetsM[i] = dofOffsetsM_[i];
-  }
-
-  dimU = dofOffsetsU[1] - dofOffsetsU[0];
-  dimM = dofOffsetsM[1] - dofOffsetsM[0];
-  dimC = dimM;
-  
-  block_offsetsx[0] = 0;
-  block_offsetsx[1] = dimU;
-  block_offsetsx[2] = dimM;
-  block_offsetsx.PartialSum();
-  
-  MPI_Allreduce(&dimU, &dimUglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-  MPI_Allreduce(&dimM, &dimMglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-  
-  ml.SetSize(dimM); ml = 0.0;
-  Vector negIdentDiag(dimM);
-  negIdentDiag = -1.0;
-  Ih = GenerateHypreParMatrixFromDiagonal(dofOffsetsM, negIdentDiag);
-}
-
-
-double ParOptProblem::CalcObjective(const BlockVector &x) const { return E(x.GetBlock(0)); }
-
-void ParOptProblem::Duf(const BlockVector &x, Vector &y) const { DdE(x.GetBlock(0), y); }
-
-void ParOptProblem::Dmf(const BlockVector &x, Vector &y) const { y = 0.0; }
-
-HypreParMatrix * ParOptProblem::Duuf(const BlockVector &x) 
-{ 
-   return DddE(x.GetBlock(0)); 
-}
-
-HypreParMatrix * ParOptProblem::Dumf(const BlockVector &x) { return nullptr; }
-
-HypreParMatrix * ParOptProblem::Dmuf(const BlockVector &x) { return nullptr; }
-
-HypreParMatrix * ParOptProblem::Dmmf(const BlockVector &x) { return nullptr; }
-
-void ParOptProblem::c(const BlockVector &x, Vector &y) const // c(u,m) = g(u) - m 
-{
-   g(x.GetBlock(0), y);
-   y.Add(-1.0, x.GetBlock(1));  
-}
-
-HypreParMatrix * ParOptProblem::Duc(const BlockVector &x) 
-{ 
-   return Ddg(x.GetBlock(0)); 
-}
-
-HypreParMatrix * ParOptProblem::Dmc(const BlockVector &x) 
-{ 
-   return Ih;
-} 
-
-ParOptProblem::~ParOptProblem() 
-{
-  delete[] dofOffsetsU;
-  delete[] dofOffsetsM;
-  delete Ih;
-}
-
-
-
-
-
-// Obstacle Problem, no essential boundary conditions enforced
-// Hessian of energy term is K + M (stiffness + mass)
-ParObstacleProblem::ParObstacleProblem(ParFiniteElementSpace *Vh_, 
-                                       double (*fSource)(const Vector &),
-				       double (*obstacleSource)(const Vector &)) :
-                                       ParOptProblem(), Vh(Vh_), J(nullptr) 
-{
-   Init(Vh->GetTrueDofOffsets(), Vh->GetTrueDofOffsets());
-   cout << "dimU = " << dimU;
-   f.SetSize(dimU); f = 0.0;
-   psi.SetSize(dimU); psi = 0.0;
-
-
-   Kform = new ParBilinearForm(Vh);
-   Kform->AddDomainIntegrator(new MassIntegrator);
-   Kform->AddDomainIntegrator(new DiffusionIntegrator);
-   Kform->Assemble();
-   Kform->Finalize();
-   Kform->FormSystemMatrix(ess_tdof_list, K);
-
-   FunctionCoefficient fcoeff(fSource);
-   fform = new ParLinearForm(Vh);
-   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
-   fform->Assemble();
-   Vector F(dimU);
-   fform->ParallelAssemble(F);
-   f.SetSize(dimU);
-   f.Set(1.0, F);
-
-   Vector iDiag(dimU); iDiag = 1.0;
-   SparseMatrix * Jacg = new SparseMatrix(iDiag);
-   
-   J = new HypreParMatrix(MPI_COMM_WORLD, dimUglb, dofOffsetsU, Jacg);
-   HypreStealOwnership(*J, *Jacg);
-   delete Jacg;
-}
-
-// Obstacle Problem, essential boundary conditions enforced
-// Hessian of energy term is K (stiffness)
-ParObstacleProblem::ParObstacleProblem(ParFiniteElementSpace *Vh_, 
-				       double (*fSource)(const Vector &),
-				       double (*obstacleSource)(const Vector &),
-				       Array<int> tdof_list, Vector &xDC) : ParOptProblem(), 
-	                                                                    Vh(Vh_), J(nullptr)
-{
-   Init(Vh->GetTrueDofOffsets(), Vh->GetTrueDofOffsets());
-   f.SetSize(dimU); f = 0.0;
-   psi.SetSize(dimU); psi = 0.0;
-   // elastic energy functional terms	
-   ess_tdof_list = tdof_list;
-   Kform = new ParBilinearForm(Vh);
-   Kform->AddDomainIntegrator(new DiffusionIntegrator);
-   Kform->Assemble();
-   Kform->Finalize();
-   Kform->FormSystemMatrix(ess_tdof_list, K);
-
-   FunctionCoefficient fcoeff(fSource);
-   fform = new ParLinearForm(Vh);
-   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
-   fform->Assemble();
-   Vector F(dimU);
-   fform->ParallelAssemble(F);
-   f.SetSize(dimU);
-   f.Set(1.0, F);
-   Kform->EliminateVDofsInRHS(ess_tdof_list, xDC, f);
-   
-   // obstacle constraints --  
-   Vector iDiag(dimU); iDiag = 1.0;
-   for(int i = 0; i < ess_tdof_list.Size(); i++)
-   {
-     iDiag(ess_tdof_list[i]) = 0.0;
-   }
-   SparseMatrix * Jacg = new SparseMatrix(iDiag);
-
-   J = new HypreParMatrix(MPI_COMM_WORLD, dimUglb, dofOffsetsU, Jacg);
-   HypreStealOwnership(*J, *Jacg);
-   delete Jacg;
-
-   FunctionCoefficient psi_fc(obstacleSource);
-   ParGridFunction psi_gf(Vh);
-   psi_gf.ProjectCoefficient(psi_fc);
-   psi.Set(1.0, (*psi_gf.GetTrueDofs()));
-   for(int i = 0; i < ess_tdof_list.Size(); i++)
-   {
-     psi(ess_tdof_list[i]) = xDC(ess_tdof_list[i]) - 1.e-8;
-   }
-}
-
-
-
-double ParObstacleProblem::E(const Vector &d) const
-{
-   Vector Kd(K.Height()); Kd = 0.0;
-   MFEM_VERIFY(d.Size() == K.Width(), "ParObstacleProblem::E - Inconsistent dimensions");
-   K.Mult(d, Kd);
-   return 0.5 * InnerProduct(MPI_COMM_WORLD, d, Kd) - InnerProduct(MPI_COMM_WORLD, f, d);
-}
-
-void ParObstacleProblem::DdE(const Vector &d, Vector &gradE) const
-{
-   gradE.SetSize(K.Height());
-   MFEM_VERIFY(d.Size() == K.Width(), "ParObstacleProblem::DdE - Inconsistent dimensions");
-   K.Mult(d, gradE);
-   MFEM_VERIFY(f.Size() == K.Height(), "ParObstacleProblem::DdE - Inconsistent dimensions");
-   gradE.Add(-1.0, f);
-}
-
-HypreParMatrix * ParObstacleProblem::DddE(const Vector &d)
-{
-   return &K; 
-}
-
-// g(d) = d >= \psi
-void ParObstacleProblem::g(const Vector &d, Vector &gd) const
-{
-   MFEM_VERIFY(d.Size() == J->Width(), "ParObstacleProblem::g - Inconsistent dimensions");
-   J->Mult(d, gd);
-   MFEM_VERIFY(gd.Size() == J->Height(), "ParObstacleProblem::g - Inconsistent dimensions");
-   gd.Add(-1.0, psi);
-}
-
-HypreParMatrix * ParObstacleProblem::Ddg(const Vector &d)
-{
-   return J;
-}
-
-ParObstacleProblem::~ParObstacleProblem()
-{
-   delete Kform;
-   delete fform;
-   delete J;
-}
-
-
-
-ReducedProblem::ReducedProblem(ParOptProblem * problem_, HYPRE_Int * constraintMask)
-{
-  problem = problem_;
-  J = nullptr;
-  P = nullptr;
-  
-  int nprocs = Mpi::WorldSize();
-  int myrank = Mpi::WorldRank();
-
-  HYPRE_BigInt * dofOffsets = problem->GetDofOffsetsU();
-
-  // given a constraint mask, lets update the constraintOffsets
-  // from the original problem
-  int nLocConstraints = 0;
-  int nProblemConstraints = problem->GetDimM();
-  for (int i = 0; i < nProblemConstraints; i++)
-  {
-    if (constraintMask[i] == 1)
-    {
-      nLocConstraints += 1;
-    }
-  }
-
-  HYPRE_BigInt * constraintOffsets_reduced;
-  constraintOffsets_reduced = offsetsFromLocalSizes(nLocConstraints);
-
-
-  for (int i = 0; i < 2; i++)
-  {
-    cout << "constraintOffsetsReduced_" << i << " = " << constraintOffsets_reduced[i] << ", (rank = " << myrank << ")\n";
-  }
-
-  HYPRE_BigInt * constraintOffsets;
-  constraintOffsets = offsetsFromLocalSizes(nProblemConstraints);
-  for (int i = 0; i < 2; i++)
-  {
-    cout << "constraintOffsets_" << i << " = " << constraintOffsets[i] << ", (rank = " << myrank << ")\n";
-  }
-  
-  
-  P = GenerateProjector(constraintOffsets, constraintOffsets_reduced, constraintMask);
-
-  Init(dofOffsets, constraintOffsets_reduced);
-  delete[] constraintOffsets_reduced;
-  delete[] constraintOffsets;
-}
-
-ReducedProblem::ReducedProblem(ParOptProblem * problem_, HypreParVector & constraintMask)
-{
-  problem = problem_;
-  J = nullptr;
-  P = nullptr;
-  
-  int nprocs = Mpi::WorldSize();
-  int myrank = Mpi::WorldRank();
-
-  HYPRE_BigInt * dofOffsets = problem->GetDofOffsetsU();
-
-  // given a constraint mask, lets update the constraintOffsets
-  // from the original problem
-  int nLocConstraints = 0;
-  int nProblemConstraints = problem->GetDimM();
-  for (int i = 0; i < nProblemConstraints; i++)
-  {
-    if (constraintMask[i] == 1)
-    {
-      nLocConstraints += 1;
-    }
-  }
-  cout << "nLocConstraints = " << nLocConstraints << ", (rank = " << myrank << ")\n";
-
-  HYPRE_BigInt * constraintOffsets_reduced;
-  constraintOffsets_reduced = offsetsFromLocalSizes(nLocConstraints);
-
-
-  for (int i = 0; i < 2; i++)
-  {
-    cout << "constraintOffsetsReduced_" << i << " = " << constraintOffsets_reduced[i] << ", (rank = " << myrank << ")\n";
-  }
-
-  HYPRE_BigInt * constraintOffsets;
-  constraintOffsets = offsetsFromLocalSizes(nProblemConstraints);
-  for (int i = 0; i < 2; i++)
-  {
-    cout << "constraintOffsets_" << i << " = " << constraintOffsets[i] << ", (rank = " << myrank << ")\n";
-  }
-  
-  
-  P = GenerateProjector(constraintOffsets, constraintOffsets_reduced, constraintMask);
-
-  Init(dofOffsets, constraintOffsets_reduced);
-  delete[] constraintOffsets_reduced;
-  delete[] constraintOffsets;
-}
-
-// energy objective E(d)
-double ReducedProblem::E(const Vector &d) const
-{
-  return problem->E(d);
-}
-
-
-// gradient of energy objective
-void ReducedProblem::DdE(const Vector &d, Vector & gradE) const
-{
-  problem->DdE(d, gradE);
-}
-
-
-HypreParMatrix * ReducedProblem::DddE(const Vector &d)
-{
-  return problem->DddE(d);
-}
-
-void ReducedProblem::g(const Vector &d, Vector &gd) const
-{
-  Vector gdfull(problem->GetDimM()); gdfull = 0.0;
-  problem->g(d, gdfull);
-  P->Mult(gdfull, gd);
-}
-
-
-HypreParMatrix * ReducedProblem::Ddg(const Vector &d)
-{
-  HypreParMatrix * Jfull = problem->Ddg(d);
-  if (J != nullptr)
-  {
-    delete J; J = nullptr;
-  }
-  J = ParMult(P, Jfull, true);
-  return J;
-}
-
-ReducedProblem::~ReducedProblem()
-{
-  delete P;
-  if (J != nullptr)
-  {
-    delete J;
-  }
-}
-
-

miniapps/interiorpointsolver/Problem.hpp

@@ -1,158 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-#include "utilities.hpp"
-
-using namespace std;
-using namespace mfem;
-
-#ifndef PARPROBLEM_DEFS
-#define PARPROBLEM_DEFS
-
-// abstract ParGeneralOptProblem class
-// of the form
-// min_(u,m) f(u,m) s.t. c(u,m)=0 and m>=ml
-// the primal variable (u, m) is represented as a BlockVector
-class ParGeneralOptProblem
-{
-protected:
-    int dimU, dimM, dimC;
-    int dimUglb, dimMglb;
-    HYPRE_BigInt * dofOffsetsU;
-    HYPRE_BigInt * dofOffsetsM;
-    Array<int> block_offsetsx;
-    Vector ml;
-public:
-    ParGeneralOptProblem();
-    virtual void Init(HYPRE_BigInt * dofOffsetsU_, HYPRE_BigInt * dofOffsetsM_);
-    virtual double CalcObjective(const BlockVector &) const = 0;
-    virtual void Duf(const BlockVector &, Vector &) const = 0;
-    virtual void Dmf(const BlockVector &, Vector &) const = 0;
-    void CalcObjectiveGrad(const BlockVector &, BlockVector &) const;
-    virtual HypreParMatrix * Duuf(const BlockVector &) = 0;
-    virtual HypreParMatrix * Dumf(const BlockVector &) = 0;
-    virtual HypreParMatrix * Dmuf(const BlockVector &) = 0;
-    virtual HypreParMatrix * Dmmf(const BlockVector &) = 0;
-    virtual HypreParMatrix * Duc(const BlockVector &) = 0;
-    virtual HypreParMatrix * Dmc(const BlockVector &) = 0;
-    virtual void c(const BlockVector &, Vector &) const = 0;
-    int GetDimU() const { return dimU; };
-    int GetDimM() const { return dimM; }; 
-    int GetDimC() const { return dimC; };
-    int GetDimUGlb() const { return dimUglb; };
-    int GetDimMGlb() const { return dimMglb; };
-    HYPRE_BigInt * GetDofOffsetsU() const { return dofOffsetsU; };
-    HYPRE_BigInt * GetDofOffsetsM() const { return dofOffsetsM; }; 
-    Vector Getml() const { return ml; };
-    ~ParGeneralOptProblem();
-};
-
-
-// abstract ContactProblem class
-// of the form
-// min_d e(d) s.t. g(d) >= 0
-class ParOptProblem : public ParGeneralOptProblem
-{
-protected:
-    HypreParMatrix * Ih;
-public:
-    ParOptProblem();
-    void Init(HYPRE_BigInt *, HYPRE_BigInt *);
-    
-    // ParGeneralOptProblem methods are defined in terms of
-    // ParOptProblem specific methods: E, DdE, DddE, g, Ddg
-    double CalcObjective(const BlockVector &) const; 
-    void Duf(const BlockVector &, Vector &) const;
-    void Dmf(const BlockVector &, Vector &) const;
-    HypreParMatrix * Duuf(const BlockVector &);
-    HypreParMatrix * Dumf(const BlockVector &);
-    HypreParMatrix * Dmuf(const BlockVector &);
-    HypreParMatrix * Dmmf(const BlockVector &);
-    void c(const BlockVector &, Vector &) const;
-    HypreParMatrix * Duc(const BlockVector &);
-    HypreParMatrix * Dmc(const BlockVector &);
-    
-    // ParOptProblem specific methods:
-    
-    // energy objective function e(d)
-    // input: d an mfem::Vector
-    // output: e(d) a double
-    virtual double E(const Vector &d) const = 0;
-
-    // gradient of energy objective De / Dd
-    // input: d an mfem::Vector,
-    //        gradE an mfem::Vector, which will be the gradient of E at d
-    // output: none    
-    virtual void DdE(const Vector &d, Vector &gradE) const = 0;
-
-    // Hessian of energy objective D^2 e / Dd^2
-    // input:  d, an mfem::Vector
-    // output: The Hessian of the energy objective at d, a pointer to a HypreParMatrix
-    virtual HypreParMatrix * DddE(const Vector &d) = 0;
-
-    // Constraint function g(d) >= 0, e.g., gap function
-    // input: d, an mfem::Vector,
-    //       gd, an mfem::Vector, which upon successfully calling the g method will be
-    //                            the evaluation of the function g at d
-    // output: none
-    virtual void g(const Vector &d, Vector &gd) const = 0;
-
-    // Jacobian of constraint function Dg / Dd, e.g., gap function Jacobian
-    // input:  d, an mfem::Vector,
-    // output: The Jacobain of the constraint function g at d, a pointer to a HypreParMatrix
-    virtual HypreParMatrix * Ddg(const Vector &) = 0;
-    virtual ~ParOptProblem();
-};
-
-
-
-
-
-
-class ParObstacleProblem : public ParOptProblem
-{
-protected:
-   // data to define energy objective function e(d) = 0.5 d^T K d - f^T d, g(d) = d >= \psi
-   // stiffness matrix used to define objective
-   ParBilinearForm *Kform;
-   ParLinearForm   *fform;
-   Array<int> ess_tdof_list; // needed for calls to FormSystemMatrix
-   HypreParMatrix  K;
-   HypreParMatrix  *J;
-   ParFiniteElementSpace *Vh;
-   Vector f;
-   Vector psi;
-public :
-   ParObstacleProblem(ParFiniteElementSpace*, double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &));
-   ParObstacleProblem(ParFiniteElementSpace*, double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &), Array<int> tdof_list, Vector &);
-   double E(const Vector &) const;
-   void DdE(const Vector &, Vector &) const;
-   HypreParMatrix* DddE(const Vector &);
-   void g(const Vector &, Vector &) const;
-   HypreParMatrix* Ddg(const Vector &);
-   virtual ~ParObstacleProblem();
-};
-
-
-
-class ReducedProblem : public ParOptProblem
-{
-protected:
-  HypreParMatrix *J;
-  HypreParMatrix *P; // projector
-  ParOptProblem  *problem;
-public:
-  ReducedProblem(ParOptProblem *problem, HYPRE_Int * constraintMask);
-  ReducedProblem(ParOptProblem *problem, HypreParVector & constraintMask);
-  double E(const Vector &) const;
-  void DdE(const Vector &, Vector &) const;
-  HypreParMatrix * DddE(const Vector &);
-  void g(const Vector &, Vector &) const;
-  HypreParMatrix * Ddg(const Vector &);
-  virtual ~ReducedProblem();
-};
-
-
-
-
-#endif

miniapps/interiorpointsolver/TestProblem1.cpp

@@ -1,181 +0,0 @@
-//                         Example Problem 1
-//
-//
-// Compile with: make ParTestProblem1
-//
-// Sample runs: mpirun -np 4 ./ParTestProblem1
-//
-//
-// Description: This example code demonstrates the use of the MFEM based
-//              interior-point solver to solve the
-//              bound-constrained minimization problem
-//
-//              minimize_(x \in R^n) 1/2 x^T x subject to x - xl ≥ 0 (component-wise).
-//              
-#include "mfem.hpp"
-#include "Problem.hpp"
-#include "IPsolver.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-
-
-// solve min 1/2 x^T K x s.t. J x - xl >= 0
-// where K and J are identity matrices and xl
-// has uniform random values in [-1, 1]
-// for the Lagrangian L(x, s, l, z) = 1/2 x^T x + l^T (x - xl - s) - z^T s
-// the optimal solution is x*_i = max{0, (xl)_i}, z*_i = x*_i
-
-
-class ParEx1Problem : public ParOptProblem
-{
-protected:
-   HypreParMatrix *K;
-   HypreParMatrix *J;
-   Vector xl;
-   //HYPRE_BigInt * dofOffsets;
-public:
-   // create offsets internally only pass problem size
-   //ParEx1Problem(HYPRE_BigInt * offsets);
-   ParEx1Problem(int n);
-   double E(const Vector &) const;
-   void DdE(const Vector &, Vector &) const;
-   HypreParMatrix* DddE(const Vector &);
-   void g(const Vector &, Vector &) const;
-   HypreParMatrix* Ddg(const Vector &);
-   virtual ~ParEx1Problem();
-};
-
-
-
-void mfemIPSolve(ParGeneralOptProblem & problem, Vector &x, Vector &lambda)
-{
-   ParInteriorPointSolver IPoptimizer(&problem);
-
-   int dimX = problem.GetDimU();
-   Vector x0(dimX); x0 = 100.0;
-   x.SetSize(dimX); x = 0.0;
-   
-   double OptTol = 1.e-6;
-   double LinSolveTol = 1.e-10;
-   int linSolveStrategy = 2;
-   int MaxOptIter = 30; 
-   IPoptimizer.SetTol(OptTol);
-   IPoptimizer.SetLinearSolveTol(LinSolveTol);
-   IPoptimizer.SetLinearSolver(linSolveStrategy);
-   IPoptimizer.SetMaxIter(MaxOptIter);
-   IPoptimizer.Mult(x0, x);
-   
-   int dimM = problem.GetDimM();
-   lambda.SetSize(dimM);
-   IPoptimizer.GetLagrangeMultiplier(lambda);
-}
-
-
-
-
-
-
-int main(int argc, char *argv[])
-{
-  // Initialize MPI
-   Mpi::Init();
-   Hypre::Init();   
-
-   int n = 10;
-   OptionsParser args(argc, argv);
-   args.AddOption(&n, "-n", "--n", \
-		   "Size of the optimization problem (dimension of primal variable)");
-   args.ParseCheck();
-   
-   
-   ParEx1Problem problem(n);
-   
-   Vector xOptimal, lambdaOptimal;
-   mfemIPSolve(problem, xOptimal, lambdaOptimal);
-   for(int i = 0; i < xOptimal.Size(); i++)
-   {
-     cout << "optimal (x, z)_" << i << " = (" << xOptimal(i) << ", " << lambdaOptimal(i) << ")\n";
-   }
-
-   Mpi::Finalize();
-   return 0;
-}
-
-
-// Ex1Problem
-// min 1/2 x^T K x such that J x - xl >= 0
-// where K and J are identity matrices
-ParEx1Problem::ParEx1Problem(int n) : ParOptProblem(), K(nullptr), J(nullptr)
-{
-  // generate the parallel partition of the 
-  // variable x and the 
-  int nprocs = Mpi::WorldSize();
-  int myrank = Mpi::WorldRank();
-  
-  
-  HYPRE_BigInt * dofOffsets = new HYPRE_BigInt[2];
-  dofOffsets[0] = HYPRE_BigInt(myrank * n / nprocs);
-  dofOffsets[1] = HYPRE_BigInt((myrank + 1) * n / nprocs);
-  
-  Init(dofOffsets, dofOffsets);
-
-  Vector iDiag(dofOffsets[1] - dofOffsets[0]); iDiag = 1.0;
-  
-  K = GenerateHypreParMatrixFromDiagonal(dofOffsets, iDiag);
-
-  J = GenerateHypreParMatrixFromDiagonal(dofOffsets, iDiag);
-
-  xl.SetSize(dofOffsets[1] - dofOffsets[0]);
-  xl.Randomize(myrank);
-  xl *= 2.0;
-  xl -= 1.0;
-  delete[] dofOffsets;
-}
-
-
-
-double ParEx1Problem::E(const Vector & x) const
-{
-   Vector Kx(K->Height()); Kx = 0.0;
-   MFEM_VERIFY(x.Size() == K->Width(), "ParEx1Problem::E - Inconsistent dimensions");
-   K->Mult(x, Kx);
-   return 0.5 * InnerProduct(MPI_COMM_WORLD, x, Kx);
-}
-
-void ParEx1Problem::DdE(const Vector &x, Vector &gradE) const
-{
-   gradE.SetSize(K->Height());
-   MFEM_VERIFY(x.Size() == K->Width(), "ParEx1Problem::DdE - Inconsistent dimensions");
-   K->Mult(x, gradE);
-}
-
-HypreParMatrix * ParEx1Problem::DddE(const Vector &x)
-{
-   return K; 
-}
-
-// g(x) = x - xl >= 0
-void ParEx1Problem::g(const Vector &x, Vector &gx) const
-{
-   MFEM_VERIFY(x.Size() == J->Width(), "ParEx1Problem::g - Inconsistent dimensions");
-   J->Mult(x, gx);
-   MFEM_VERIFY(gx.Size() == J->Height(), "ParEx1Problem::g - Inconsistent dimensions");
-   gx.Add(-1.0, xl);
-}
-
-HypreParMatrix * ParEx1Problem::Ddg(const Vector &)
-{
-  return J;
-}
-
-ParEx1Problem::~ParEx1Problem()
-{
-   delete K;
-   delete J;
-}
-

miniapps/interiorpointsolver/TestProblem2.cpp

@@ -1,161 +0,0 @@
-//                       Spherical Obstacle Problem
-//
-//
-// Compile with: make ParSphericalObstacleProblem
-//
-// Sample runs: mpirun -np 4 ./ParSphericalObstacleProblem -linSolver 0
-//              mpirun -np 4 ./ParSphericalObstacleProblem -linSolver 1
-//              mpirun -np 4 ./ParSphericalObstacleProblem -linSolver 2
-//
-//
-// Description: This example code demonstrates the use of MFEM to solve the
-//              bound-constrained energy minimization problem
-//
-//                      minimize ||∇u||² subject to u ≥ ϕ in H¹₀.
-
-#include "mfem.hpp"
-#include "Problem.hpp"
-#include "IPsolver.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-double fRhs(const Vector &);
-double spherical_obstacle(const Vector &);
-double exact_solution_obstacle(const Vector &);
-
-
-int main(int argc, char *argv[])
-{
-  // Initialize MPI
-   Mpi::Init();
-   int num_procs = Mpi::WorldSize();
-   int myid = Mpi::WorldRank();
-   Hypre::Init();
-
-   int FEorder = 1; // order of the finite elements
-   int linSolver = 2;
-   int maxIPMiters = 30;
-   int ref_levels = 3;
-   OptionsParser args(argc, argv);
-   args.AddOption(&FEorder, "-o", "--order",\
-	 	  "Order of the finite elements.");
-   args.AddOption(&linSolver, "-linSolver", "--linearSolver", \
-        "IP-Newton linear system solution strategy.");
-   args.AddOption(&maxIPMiters, "-IPMiters", "--IPMiters",\
-	 	  "Maximum number of IPM iterations");
-   args.AddOption(&ref_levels, "-r", "--mesh_refinement", \
-		  "Mesh Refinement");
-  
-   args.ParseCheck();
-
-   const char *meshFile = "disk.mesh";
-   Mesh mesh(meshFile, 1, 1);
-   int dim = mesh.Dimension(); // geometric dimension of the meshed domain
-   {
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   
-   FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
-   ParFiniteElementSpace   *Vh  = new ParFiniteElementSpace(&pmesh, fec);
-   Array<int> boundary_dofs;
-   Vh->GetBoundaryTrueDofs(boundary_dofs);
-   int dimD = Vh->GetTrueVSize();
-   Vector xDC(dimD); xDC = 0.0;
-
-   ParObstacleProblem problem(Vh, &fRhs, &spherical_obstacle, boundary_dofs, xDC);
-   Vector x0(dimD); x0.Set(1.0, xDC);
-   Vector xf(dimD); xf = 0.0;
-
-   ParInteriorPointSolver optimizer(&problem); 
-   optimizer.SetTol(1.e-7);
-   optimizer.SetLinearSolveTol(1.e-9);
-   optimizer.SetLinearSolver(linSolver);
-   optimizer.SetMaxIter(maxIPMiters);
-   optimizer.Mult(x0, xf);
-
-   ParGridFunction d_gf(Vh);
-
-   d_gf.SetFromTrueDofs(xf);
-
-
-   FunctionCoefficient dtrue_fc(exact_solution_obstacle); // analytic solution
-   ParGridFunction dtrue_gf(Vh);
-   dtrue_gf.ProjectCoefficient(dtrue_fc);
-
-   double L2error = d_gf.ComputeL2Error(dtrue_fc);
-   if (myid == 0)
-   {
-     cout << "\n|| u_h - u ||_{L^2} = " << L2error << '\n' << endl;
-   } 
-   
-   ParaViewDataCollection paraview_dc("SphericalObstacleProblem", &pmesh);
-   paraview_dc.SetPrefixPath("ParaView");
-   paraview_dc.SetLevelsOfDetail(FEorder);
-   paraview_dc.SetDataFormat(VTKFormat::BINARY);
-   paraview_dc.SetHighOrderOutput(true);
-   paraview_dc.SetCycle(0);
-   paraview_dc.SetTime(0.0);
-   paraview_dc.RegisterField("u(x,y) (analytic)", &dtrue_gf);
-   paraview_dc.RegisterField("u(x,y) (numerical)", &d_gf);
-   paraview_dc.Save();
-
-   delete Vh;
-   delete fec;
-   return 0;
-}
-
-
-double fRhs(const Vector &x)
-{
-  return 0.;
-}
-
-
-double spherical_obstacle(const Vector &pt)
-{
-   double x = pt(0), y = pt(1);
-   double r = sqrt(x*x + y*y);
-   double r0 = 0.5;
-   double beta = 0.9;
-
-   double b = r0*beta;
-   double tmp = sqrt(r0*r0 - b*b);
-   double B = tmp + b*b/tmp;
-   double C = -b/tmp;
-
-   if (r > b)
-   {
-      return B + r * C;
-   }
-   else
-   {
-      return sqrt(r0*r0 - r*r);
-   }
-}
-
-double exact_solution_obstacle(const Vector &pt)
-{
-   double x = pt(0), y = pt(1);
-   double r = sqrt(x*x + y*y);
-   double r0 = 0.5;
-   double a =  0.348982574111686;
-   double A = -0.340129705945858;
-
-   if (r > a)
-   {
-      return A * log(r);
-   }
-   else
-   {
-      return sqrt(r0*r0-r*r);
-   }
-} 
-

miniapps/interiorpointsolver/disk.mesh

@@ -1,109 +0,0 @@
-MFEM NURBS mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# SEGMENT     = 1
-# SQUARE      = 3
-# CUBE        = 5
-#
-
-dimension
-2
-
-elements
-5
-1 3 4 5 6 7
-1 3 0 1 5 4
-1 3 1 2 6 5
-1 3 3 7 6 2
-1 3 0 4 7 3
-
-boundary
-4
-1 1 0 1
-1 1 2 3
-1 1 1 2
-1 1 3 0
-
-edges
-12
-0 0 1
-0 4 5
-0 7 6
-0 3 2
-1 1 2
-1 5 6
-1 4 7
-1 0 3
-2 0 4
-2 1 5
-2 2 6
-2 3 7
-
-vertices
-8
-
-knotvectors
-3
-2 3 0 0 0 1 1 1
-2 3 0 0 0 1 1 1
-2 3 0 0 0 1 1 1
-
-weights
-1
-1
-1
-1
-1
-1
-1
-1
-0.70710678118655
-1
-1
-0.70710678118655
-0.70710678118655
-1
-1
-0.70710678118655
-1
-1
-1
-1
-1
-0.85355339059327
-0.85355339059327
-0.85355339059327
-0.85355339059327
-
-FiniteElementSpace
-FiniteElementCollection: NURBS2
-VDim: 2
-Ordering: 1
-
--0.70710678118 -0.70710678118
-0.70710678118 -0.70710678118
-0.70710678118 0.70710678118
--0.70710678118 0.70710678118
--0.35355339059 -0.35355339059
-0.35355339059 -0.35355339059
-0.35355339059 0.35355339059
--0.35355339059 0.35355339059
-0 -1.41421356236
-0 -0.35355339059
-0 0.35355339059
-0 1.41421356236
-1.41421356236 0
-0.35355339059 0
--0.35355339059 0
--1.41421356236 0
--0.530330085885 -0.530330085885
-0.530330085885 -0.530330085885
-0.530330085885 0.530330085885
--0.530330085885 0.530330085885
-0 0
-0 -0.883883476475
-0.883883476475 0
-0 0.883883476475
--0.883883476475 0

miniapps/interiorpointsolver/inline-quad.mesh

@@ -1,7 +0,0 @@
-MFEM INLINE mesh v1.0
-
-type = quad
-nx = 4
-ny = 4
-sx = 1.0
-sy = 1.0

miniapps/interiorpointsolver/makefile

@@ -1,35 +0,0 @@
-# Use the MFEM build directory
-MFEM_DIR ?= ../..
-MFEM_BUILD_DIR ?= ../..
-SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/interiorpointsolver/,)
-CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
-   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
-
-MFEM_LIB_FILE = mfem_is_not_built
--include $(CONFIG_MK)
-
-# Remove built-in rule
-%: %.cpp
-%: %.c
-%: %.o
-
-%.o: %.cpp $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
-
-FRAMEWORK_SRC = Problem.cpp IPsolver.cpp utilities.cpp
-FRAMEWORK_OBJ = $(FRAMEWORK_SRC:.cpp=.o)
-
-TestProblem1: TestProblem1.o $(FRAMEWORK_OBJ) $(MFEM_LIB_FILE)
-	$(MFEM_CXX) $(MFEM_FLAGS) TestProblem1.o $(FRAMEWORK_OBJ) -o $@ $(MFEM_LIBS)
-
-TestProblem2: TestProblem2.o $(FRAMEWORK_OBJ) $(MFEM_LIB_FILE)
-	$(MFEM_CXX) $(MFEM_FLAGS) TestProblem2.o $(FRAMEWORK_OBJ) -o $@ $(MFEM_LIBS)
-
-# Generate an error message if the MFEM library is not built and exit
-$(MFEM_LIB_FILE):
-	$(error The MFEM library is not built)
-
-
-.PHONY: clean
-clean:
-	rm -f *.o TestProblem1 TestProblem2

miniapps/interiorpointsolver/utilities.cpp

@@ -1,178 +0,0 @@
-#include "mfem.hpp"
-#include "utilities.hpp"
-
-using namespace mfem;
-
-HypreParMatrix * GenerateHypreParMatrixFromSparseMatrix(HYPRE_BigInt * colOffsetsloc, HYPRE_BigInt * rowOffsetsloc, SparseMatrix * Asparse)
-{
-  int ncols_loc = colOffsetsloc[1] - colOffsetsloc[0];
-  int nrows_loc = rowOffsetsloc[1] - rowOffsetsloc[0];
-  HYPRE_BigInt ncols_glb, nrows_glb;
-  MPI_Allreduce(&nrows_loc, &nrows_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-  MPI_Allreduce(&ncols_loc, &ncols_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-
-  int * AI          = Asparse->GetI();
-  HYPRE_BigInt * AJ = Asparse->GetJ();
-  double * Adata    = Asparse->GetData();
-
-  HypreParMatrix * Ahypre = nullptr;
-  Ahypre = new HypreParMatrix(MPI_COMM_WORLD, nrows_loc, nrows_glb, ncols_glb, AI, AJ, Adata, rowOffsetsloc, colOffsetsloc);
-  return Ahypre;
-}
-
-
-HypreParMatrix * GenerateHypreParMatrixFromDiagonal(HYPRE_BigInt * offsetsloc, 
-		Vector & diag)
-{
-   int n_loc = offsetsloc[1] - offsetsloc[0];
-   int n_glb = 0;
-   MPI_Allreduce(&n_loc, &n_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-   
-   SparseMatrix * Dsparse = new SparseMatrix(n_loc, n_glb);
-   Array<int> cols;
-   Vector entries;
-   cols.SetSize(1);
-   entries.SetSize(1);
-   for(int j = 0; j < n_loc; j++)
-   {
-     cols[0] = offsetsloc[0] + j;
-     entries(0) = diag(j);
-     Dsparse->SetRow(j, cols, entries);
-   }   
-   Dsparse->Finalize();
-   HypreParMatrix * Dhypre = nullptr;
-   Dhypre = GenerateHypreParMatrixFromSparseMatrix(offsetsloc, offsetsloc, Dsparse);
-   delete Dsparse;
-   return Dhypre;   
-}
-
-HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, HYPRE_Int * mask)
-{
-  int n_cols_loc = offsets[1] - offsets[0];
-  int n_cols_glb = 0;
-  MPI_Allreduce(&n_cols_loc, &n_cols_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-
-  int n_rows_loc = reduced_offsets[1] - reduced_offsets[0];
-
-  SparseMatrix * Psparse = new SparseMatrix(n_rows_loc, n_cols_glb);
-  Array<int> cols;
-  Vector entries;
-  cols.SetSize(1);
-  entries.SetSize(1);
-
-  int row = 0;
-  for(int j = 0; j < n_cols_loc; j++)
-  {
-    if (mask[j] == 1)
-    {
-      cols[0] = offsets[0] + j;
-      entries(0) = 1.0;
-      Psparse->SetRow(row, cols, entries);
-      row += 1;
-    }
-  }
-  Psparse->Finalize();
-  HypreParMatrix * Phypre = nullptr;
-  Phypre = GenerateHypreParMatrixFromSparseMatrix(offsets, reduced_offsets, Psparse);
-  delete Psparse;
-  return Phypre;
-}
-
-HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, const HypreParVector & mask)
-{
-  int n_cols_loc = offsets[1] - offsets[0];
-  int n_cols_glb = 0;
-  MPI_Allreduce(&n_cols_loc, &n_cols_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-
-  int n_rows_loc = reduced_offsets[1] - reduced_offsets[0];
-
-  SparseMatrix * Psparse = new SparseMatrix(n_rows_loc, n_cols_glb);
-  Array<int> cols;
-  Vector entries;
-  cols.SetSize(1);
-  entries.SetSize(1);
-
-  int row = 0;
-  for(int j = 0; j < n_cols_loc; j++)
-  {
-    if (mask(j) > 0.5)
-    {
-      cols[0] = offsets[0] + j;
-      entries(0) = 1.0;
-      Psparse->SetRow(row, cols, entries);
-      row += 1;
-    }
-  }
-  Psparse->Finalize();
-  HypreParMatrix * Phypre = nullptr;
-  Phypre = GenerateHypreParMatrixFromSparseMatrix(offsets, reduced_offsets, Psparse);
-  delete Psparse;
-  return Phypre;
-}
-
-
-
-
-HYPRE_BigInt * offsetsFromLocalSizes(int n)
-{
-  HYPRE_BigInt * offsets = new HYPRE_BigInt[2];
-  
-  int nprocs = Mpi::WorldSize();
-  int myrank = Mpi::WorldRank();
-  
-  if (myrank == 0)
-  {
-    offsets[0] = 0;
-    offsets[1] = n;
-  }
-  else
-  {
-    offsets[0] = 0;
-    offsets[1] = 0;
-  }
-  
-  // receive then send
-  
-  // Receive local size info from processes with rank less than myrank 
-  // Populate that as entries of helper
-  HYPRE_BigInt * helper;
-  if (myrank > 0)
-  {
-    helper = new HYPRE_BigInt[myrank];
-  }
-  int tag;
-  for (int i = 0; i < myrank; i++)
-  {
-    tag = myrank + i * nprocs;
-    MPI_Recv (&(helper[i]), 1, MPI_INT, i, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offsets[0] += helper[i];
-  }
-
-  if (myrank > 0)
-  {
-    delete[] helper;
-  }
-  offsets[1] = offsets[0] + n;
-  
-  // Send local size info to all processes with rank greater than myrank
-  for (int i = myrank + 1; i < nprocs; i++)
-  {
-    tag = i + myrank * nprocs;
-    MPI_Send (&n, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
-  }
-  return offsets;
-}
-
-
-void HypreToMfemOffsets(HYPRE_BigInt * offsets)
-{
-  if (offsets[1] < offsets[0])
-  {
-    offsets[1] = offsets[0];
-  }
-  else
-  {
-    offsets[1] = offsets[1] + 1;
-  }
-}
-

miniapps/interiorpointsolver/utilities.hpp

@@ -1,23 +0,0 @@
-#include "mfem.hpp"
-
-using namespace mfem;
-#ifndef UTILITY_FUNCTIONS
-#define UTILITY_FUNCTIONS
-
-void HypreToMfemOffsets(HYPRE_BigInt * offsets);
-
-HypreParMatrix * GenerateHypreParMatrixFromSparseMatrix(HYPRE_BigInt * colOffsetsloc, HYPRE_BigInt * rowOffsetsloc, SparseMatrix * Asparse);
-
-HypreParMatrix * GenerateHypreParMatrixFromDiagonal(HYPRE_BigInt * offsetsloc, 
-		Vector & diag);
-
-
-HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, HYPRE_Int * mask);
-
-HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, const HypreParVector & mask);
-
-
-HYPRE_BigInt * offsetsFromLocalSizes(int n);
-
-
-#endif