including an obstacle problem variant wherein the essential dofs are not given to the optimizer as free variables that it must choose the values for.

data/disk.mesh

@@ -0,0 +1,109 @@
+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

examples/contact/DirichletObstacleProblem.cpp

@@ -0,0 +1,141 @@
+#include "mfem.hpp"
+#include "Problems.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+
+double fRhs(const Vector &pt);
+double obstacle(const Vector &pt);
+double dmanufacturedFun(const Vector &pt);
+
+int main(int argc, char *argv[])
+{
+  int FEorder = 1; // order of the finite elements
+  int linSolver = 0;
+  int maxIPMiters = 30;
+  bool iAmRoot = true;
+  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.Parse();
+  if(!args.Good())
+  {
+    args.PrintUsage(cout);
+    return 1;
+  }
+  else
+  {
+    if( iAmRoot )
+    {
+      args.PrintOptions(cout);
+    }
+  }
+
+  const char *meshFile = "../../data/inline-quad.mesh";
+  Mesh *mesh = new Mesh(meshFile, 1, 1);
+  int dim = mesh->Dimension(); // geometric dimension of the domain
+  for (int l = 0; l < ref_levels; l++)
+  {
+      mesh->UniformRefinement();
+  }
+
+  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
+  Array<int> ess_tdof_list;
+  if (mesh->bdr_attributes.Size())
+  {
+     Array<int> ess_bdr(mesh->bdr_attributes.Max());
+     ess_bdr = 1;
+     Vh->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+  }
+  
+  double DC_val = 0.0;
+  int dimD = Vh->GetTrueVSize();
+  Vector x0(dimD); x0 = DC_val;
+  Vector xf(dimD); xf = 0.0;
+
+  ObstacleProblem problem(Vh, x0, &fRhs, &obstacle, ess_tdof_list);
+  
+
+  InteriorPointSolver optimizer(&problem); 
+  optimizer.SetTol(1.e-7);
+  optimizer.SetLinearSolver(linSolver);
+  optimizer.SetMaxIter(maxIPMiters);
+  optimizer.Mult(x0, xf);
+  
+  double Einitial = problem.E(x0);
+  double Efinal = problem.E(xf);
+  cout << "Energy objective at initial point = " << Einitial << endl;
+  cout << "Energy objective at optimizer = " << Efinal << endl;
+
+  GridFunction d_gf(Vh);
+
+  d_gf = xf;
+  
+  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
+  GridFunction dm_gf(Vh);
+  dm_gf.ProjectCoefficient(dm_fc);
+
+  
+  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
+  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("d(x) (numerical)", &d_gf);
+  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
+  paraview_dc.Save();
+  
+
+  delete Vh;
+  delete fec;
+  delete mesh;
+  return 0;
+}
+
+
+
+double dmanufacturedFun(const Vector &pt)
+{
+  double alpha = 16.5;
+  return sin(M_PI * pt(1)) * (sin(M_PI * pt(0)) - alpha * pow(pt(0) * (1. - pt(0)), 2));
+}
+
+
+// f(x) forcing term... which enters the objective energy functional
+// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
+// of f(x). f(x) is such that in the absence of bound-constraints then
+// the solution of the optimization problem satisfies the PDE
+// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
+// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
+
+double fRhs(const Vector &pt)
+{
+  double alpha = 16.5;
+  double fx;
+  fx = pow(M_PI, 2) * sin(M_PI * pt(0));
+  fx += alpha * (2. * pow(pt(0), 2) + 2. * pow(1.-pt(0), 2) - 8. * pt(0) * (1.-pt(0)));
+  fx += pow(M_PI, 2) * sin(M_PI * pt(0)) * dmanufacturedFun(pt);
+  fx *= sin(M_PI * pt(1));
+  return fx;
+}
+
+double obstacle(const Vector &pt)
+{
+  return 0.0;
+}

examples/contact/DirichletObstacleProblemVariant.cpp

@@ -0,0 +1,156 @@
+#include "mfem.hpp"
+#include "Problems.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+
+double fRhs(const Vector &pt);
+double obstacle(const Vector &pt);
+double dmanufacturedFun(const Vector &pt);
+
+int main(int argc, char *argv[])
+{
+  int FEorder = 1; // order of the finite elements
+  int linSolver = 0;
+  int maxIPMiters = 30;
+  bool iAmRoot = true;
+  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.Parse();
+  if(!args.Good())
+  {
+    args.PrintUsage(cout);
+    return 1;
+  }
+  else
+  {
+    if( iAmRoot )
+    {
+      args.PrintOptions(cout);
+    }
+  }
+
+  const char *meshFile = "../../data/inline-quad.mesh";
+  Mesh *mesh = new Mesh(meshFile, 1, 1);
+  int dim = mesh->Dimension(); // geometric dimension of the domain
+  for (int l = 0; l < ref_levels; l++)
+  {
+      mesh->UniformRefinement();
+  }
+
+  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
+  Array<int> ess_tdof_list;
+  if (mesh->bdr_attributes.Size())
+  {
+     Array<int> ess_bdr(mesh->bdr_attributes.Max());
+     ess_bdr = 1;
+     Vh->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+  }
+  
+  double DC_val = 0.06;
+  Vector x0DC(Vh->GetTrueVSize()); x0DC = DC_val;
+  int dimD = Vh->GetTrueVSize() - ess_tdof_list.Size();
+  Vector x0(dimD); x0 = 0.0;
+  Vector xf(dimD); xf = 0.0;
+
+  ObstacleProblemVariant problem(Vh, x0DC, &fRhs, &obstacle, ess_tdof_list);
+
+  InteriorPointSolver optimizer(&problem); 
+  optimizer.SetTol(1.e-7);
+  optimizer.SetLinearSolver(linSolver);
+  optimizer.SetMaxIter(maxIPMiters);
+  optimizer.Mult(x0, xf);
+  
+  double Einitial = problem.E(x0);
+  double Efinal = problem.E(xf);
+  cout << "Energy objective at initial point = " << Einitial << endl;
+  cout << "Energy objective at optimizer = " << Efinal << endl;
+
+  Array<int> noness_tdof_list;
+  noness_tdof_list.SetSize(dimD);
+  int i = 0;
+  for(int j = 0; j < Vh->GetTrueVSize(); j++)
+  {
+     if(ess_tdof_list.Find(j) == -1)
+     {
+        noness_tdof_list[i] = j;
+        i += 1;
+     }
+  }
+
+
+  GridFunction d_gf(Vh);
+  d_gf.Set(1.0, x0DC);
+  d_gf.SetSubVector(noness_tdof_list, xf);
+
+
+
+  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
+  GridFunction dm_gf(Vh);
+  dm_gf.ProjectCoefficient(dm_fc);
+
+  
+  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
+  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("d(x) (numerical)", &d_gf);
+  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
+  paraview_dc.Save();
+ 
+
+  delete Vh;
+  delete fec;
+  delete mesh;
+  return 0;
+}
+
+
+
+double dmanufacturedFun(const Vector &pt)
+{
+  double alpha = 16.5;
+  return sin(M_PI * pt(1)) * (sin(M_PI * pt(0)) - alpha * pow(pt(0) * (1. - pt(0)), 2));
+}
+
+
+// f(x) forcing term... which enters the objective energy functional
+// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
+// of f(x). f(x) is such that in the absence of bound-constraints then
+// the solution of the optimization problem satisfies the PDE
+// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
+// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
+
+double fRhs(const Vector &pt)
+{
+  double alpha = 16.5;
+  double fx;
+  fx = pow(M_PI, 2) * sin(M_PI * pt(0));
+  fx += alpha * (2. * pow(pt(0), 2) + 2. * pow(1.-pt(0), 2) - 8. * pt(0) * (1.-pt(0)));
+  fx += pow(M_PI, 2) * sin(M_PI * pt(0)) * dmanufacturedFun(pt);
+  fx *= sin(M_PI * pt(1));
+  return fx;
+}
+
+double obstacle(const Vector &pt)
+{
+  return 0.0;
+}

examples/contact/IPsolver.cpp

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

examples/contact/IPsolver.hpp

@@ -0,0 +1,80 @@
+#include "mfem.hpp"
+#include "Problems.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+#ifndef IPSOLVER 
+#define IPSOLVER
+
+class InteriorPointSolver
+{
+protected:
+    GeneralOptProblem* optProblem;
+    double rel_tol;
+    int  max_iter;
+    double mu_k; // \mu_k
+    Vector lk, zlk;
+
+    double sMax, kSig, tauMin, eta, thetaMin, delta, sTheta, sPhi, kMu, thetaMu;
+    double thetaMax, kSoc, gTheta, gPhi, kEps;
+	
+    // filter
+    Array<double> F1, F2;
+	
+    // quantities computed in lineSearch
+    double alpha, alphaz;
+    double thx0, thxtrial;
+    double phx0, phxtrial;
+    bool descentDirection, switchCondition, sufficientDecrease, lineSearchSuccess, inFilterRegion;
+    double Dxphi0_xhat;
+
+    int dimU, dimM, dimC;
+    Array<int> block_offsetsumlz, block_offsetsuml, block_offsetsx;
+    Vector ml;
+
+    Vector ckSoc;
+    SparseMatrix * Huu, * Hum, * Hmu, * Hmm, * Wmm, *D, * Ju, * Jm, * JuT, * JmT;
+    SparseMatrix * Huucl, *HLuu;    
+    int jOpt;
+    bool converged;
+    
+    int MyRank;
+    bool iAmRoot;
+
+    bool saveLogBarrierIterates;
+
+    int linSolver;
+
+public:
+    InteriorPointSolver(GeneralOptProblem*);
+    void Mult(const BlockVector& , BlockVector&); // used when the user wants to be aware of bound-constrained variable m >= ml
+    void Mult(const Vector&, Vector &); // useful when the user doesn't need to know about bound-constrained variable   m >= ml, e.g., when m is a slack variable
+    double MaxStepSize(Vector& , Vector& , Vector& , double);
+    double MaxStepSize(Vector& , Vector& , double);
+    void FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
+    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , double, bool);
+    void lineSearch(BlockVector& , BlockVector& , double);
+    void projectZ(const Vector & , Vector &, double);
+    void filterCheck(double, double);
+    double E(const BlockVector &, const Vector &, const Vector &, double, bool);
+    double E(const BlockVector &, const Vector &, const Vector &, bool);
+    bool GetConverged() const;
+    // TO DO: include Hessian of Lagrangian
+    double theta(const BlockVector &);
+    double phi(const BlockVector &, double);
+    void Dxphi(const BlockVector &, double, BlockVector &);
+    double L(const BlockVector &, const Vector &, const Vector &);
+    void DxL(const BlockVector &, const Vector &, const Vector &, BlockVector &);
+    void SetTol(double);
+    void SetMaxIter(int);
+    void SetBarrierParameter(double);    
+    void SaveLogBarrierHessianIterates(bool);
+    void SetLinearSolver(int);
+    virtual ~InteriorPointSolver();
+};
+
+#endif

examples/contact/ObstacleProblem.cpp

@@ -0,0 +1,125 @@
+#include "mfem.hpp"
+#include "Problems.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+
+double dmanufacturedFun(const Vector &);
+double fRhs(const Vector &);
+double obstacle(const Vector &);
+
+int main(int argc, char *argv[])
+{
+  int FEorder = 1; // order of the finite elements
+  int linSolver = 0;
+  int maxIPMiters = 30;
+  bool iAmRoot = true;
+  
+  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.Parse();
+  if(!args.Good())
+  {
+    args.PrintUsage(cout);
+    return 1;
+  }
+  else
+  {
+    if( iAmRoot )
+    {
+      args.PrintOptions(cout);
+    }
+  }
+
+  const char *meshFile = "../../data/inline-quad.mesh";
+  Mesh *mesh = new Mesh(meshFile, 1, 1);
+  int dim = mesh->Dimension(); // geometric dimension of the domain
+  {
+     int ref_levels = 3;
+     for (int l = 0; l < ref_levels; l++)
+     {
+        mesh->UniformRefinement();
+     }
+  }
+
+  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
+  ObstacleProblem problem(Vh, &fRhs, &obstacle);
+  
+  int dimD = problem.GetDimD();
+  Vector x0(dimD); x0 = 0.0;
+  Vector xf(dimD); xf = 0.0;
+
+  InteriorPointSolver optimizer(&problem); 
+  optimizer.SetTol(1.e-7);
+  optimizer.SetLinearSolver(linSolver);
+  optimizer.SetMaxIter(maxIPMiters);
+  optimizer.Mult(x0, xf);
+
+  double Einitial = problem->E(x0);
+  double Efinal = problem->E(xf);
+  cout << "Energy objective at initial point = " << Einitial << endl;
+  cout << "Energy objective at QP optimizer = " << Efinal << endl;
+
+
+  GridFunction d_gf(Vh);
+
+  d_gf = xf;
+
+  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
+  GridFunction dm_gf(Vh);
+  dm_gf.ProjectCoefficient(dm_fc);
+  
+  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
+  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("d(x) (numerical)", &d_gf);
+  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
+  paraview_dc.Save();
+  
+  
+  delete Vh;
+  delete fec;
+  delete mesh;
+  return 0;
+}
+
+
+double dmanufacturedFun(const Vector &x)
+{
+  return cos(2*M_PI*x(0)) + 0.2 - 2.0*(pow(x(0),3) - 1.5*pow(x(0),2));
+}
+
+
+// f(x) forcing term... which enters the objective energy functional
+// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
+// of f(x). f(x) is such that in the absence of bound-constraints then
+// the solution of the optimization problem satisfies the PDE
+// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
+// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
+double fRhs(const Vector &x)
+{
+  double fx = 0.;
+  fx = 0.2 - 2.0 * (pow(x(0),3)- 1.5*pow(x(0),2.) - 6 * x(0) + 3.) + (1. + pow(2.*M_PI,2))*cos(2.*M_PI*x(0));
+  return fx;
+}
+
+double obstacle(const Vector &x)
+{
+  return 0.0;
+}

examples/contact/ParIPsolver.cpp

@@ -0,0 +1,834 @@
+#include "mfem.hpp"
+#include "ParIPsolver.hpp"
+#include "ParProblems.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), 
+                       saveLogBarrierIterates(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.8;     // control rate at which iterates can approach the boundary
+   eta      = 1.e-4;   // backtracking constant
+   thetaMin = 1.e-4;   // allowed violation of the equality constraints
+
+   // constants in line-step A-5.4
+   delta    = 1.0;
+   sTheta   = 1.1;
+   sPhi     = 2.3;
+
+   // control the rate at which the penalty parameter is decreased
+   kMu     = 0.2;
+   thetaMu = 1.5;
+
+   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();
+   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);
+   ckSoc.SetSize(dimC);
+  
+   block_offsetsumlz[0] = 0;
+   block_offsetsumlz[1] = dimU; // u
+   block_offsetsumlz[2] = dimM; // m
+   block_offsetsumlz[3] = dimC; // lambda
+   block_offsetsumlz[4] = dimM; // zl
+   block_offsetsumlz.PartialSum();
+
+   for(int i = 0; i < block_offsetsuml.Size(); i++)  
+   { 
+      block_offsetsuml[i] = block_offsetsumlz[i]; 
+   }
+   for(int i = 0; i < block_offsetsx.Size(); i++)    
+   { 
+      block_offsetsx[i] = block_offsetsuml[i] ; 
+   }
+
+  
+   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);
+   x0block.GetBlock(1) = 100.;
+   x0block.GetBlock(1).Add(1.0, ml);
+   BlockVector xfblock(block_offsetsx); xfblock = 0.0;
+   Mult(x0block, xfblock);
+   xf.Set(1.0, xfblock.GetBlock(0));
+}
+
+
+void 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;
+
+   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;
+      // why do we have Xhatuml ....???
+      // TO DO: remove Xhatuml in favor of passing Xhat
+      IPNewtonSolve(xk, lk, zlk, zlhat, Xhatuml, mu_k, false); 
+
+      // assign data stack, X = (u, m, l, zl)
+      Xk = 0.0;
+      Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
+      Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
+      Xk.GetBlock(2).Set(1.0, lk);
+      Xk.GetBlock(3).Set(1.0, zlk);
+
+      // assign data stack, Xhat = (uhat, mhat, lhat, zlhat)
+      Xhat = 0.0;
+      for(int i = 0; i < 3; i++)
+      {
+         Xhat.GetBlock(i).Set(1.0, Xhatuml.GetBlock(i));
+      }
+      Xhat.GetBlock(3).Set(1.0, zlhat);
+
+      // A-5. Backtracking line search.
+      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
+         // print info regarding zl...
+         xk.GetBlock(0).Add(alpha, Xhat.GetBlock(0));
+         xk.GetBlock(1).Add(alpha, Xhat.GetBlock(1));
+         lk.Add(alpha,   Xhat.GetBlock(2));
+         zlk.Add(alphaz, Xhat.GetBlock(3));
+         projectZ(xk, zlk, mu_k);
+      }
+      else
+      {
+         if(iAmRoot)
+         {
+            cout << "lineSearch not successful :(\n";
+            cout << "attempting feasibility restoration with theta = " << thx0 << endl;
+            cout << "no feasibility restoration implemented, exiting now \n";
+         }
+         break;
+      }
+      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(saveLogBarrierIterates)
+   {
+      std::ofstream diagStream;
+      char diagString[100];
+      snprintf(diagString, 100, "logBarrierHessiandata/D%d.dat", jOpt);
+      diagStream.open(diagString, ios::out | ios::trunc);
+      for(int ii = 0; ii < dimM; ii++)
+      {
+         diagStream << setprecision(30) << DiagLogBar(ii) << endl;
+      }
+      diagStream.close();
+   } 
+
+   SparseMatrix * Ds = new SparseMatrix(DiagLogBar);
+   ParFiniteElementSpace * fes = problem->GetfesM();
+   D = new HypreParMatrix(fes->GetComm(), fes->GlobalTrueVSize(), fes->GetTrueDofOffsets(), Ds);
+   HypreStealOwnership(*D,*Ds);
+   delete Ds;
+  
+   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, 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 *>(&(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)
+   {
+      // form A = Huu + Ju^T D Ju, Wmm = D for contact
+      HypreParMatrix * Huuloc = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(0, 0)));
+      HypreParMatrix * Wmmloc = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(1, 1)));
+      HypreParMatrix * Juloc  = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(2, 0)));
+      HypreParMatrix * JuTloc = dynamic_cast<HypreParMatrix *>(&(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)) );
+   }
+
+   // 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  = max(tauMin, 1.0 - mu);
+   Vector u0   = X0.GetBlock(0);
+   Vector m0   = X0.GetBlock(1);
+   Vector l0   = X0.GetBlock(2);
+   Vector z0   = X0.GetBlock(3);
+   Vector uhat = Xhat.GetBlock(0);
+   Vector mhat = Xhat.GetBlock(1);
+   Vector lhat = Xhat.GetBlock(2);
+   Vector zhat = Xhat.GetBlock(3);
+   double alphaMax  = MaxStepSize(m0, ml, mhat, tau);
+   double alphaMaxz = MaxStepSize(z0, zhat, tau);
+   alphaz = alphaMaxz;
+
+   BlockVector x0(block_offsetsx); x0 = 0.0;
+   x0.GetBlock(0).Set(1.0, u0);
+   x0.GetBlock(1).Set(1.0, m0);
+   
+   BlockVector xhat(block_offsetsx); xhat = 0.0;
+   xhat.GetBlock(0).Set(1.0, uhat);
+   xhat.GetBlock(1).Set(1.0, mhat);
+   
+   BlockVector xtrial(block_offsetsx); xtrial = 0.0;
+   BlockVector Dxphi0(block_offsetsx); Dxphi0 = 0.0;
+   int maxBacktrack = 20;
+   alpha = alphaMax;
+
+   Vector ck0(dimC); ck0 = 0.0;
+   Vector zhatsoc(dimM); zhatsoc = 0.0;
+   BlockVector Xhatumlsoc(block_offsetsuml); Xhatumlsoc = 0.0;
+   BlockVector xhatsoc(block_offsetsx); xhatsoc = 0.0;
+   Vector uhatsoc(dimU); uhatsoc = 0.0;
+   Vector mhatsoc(dimM); mhatsoc = 0.0;
+
+   Dxphi(x0, mu, Dxphi0);
+
+   Dxphi0_xhat = InnerProduct(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;
+            }
+         }
+         // A-5.5: Initialize the second-order correction
+         if((!(thx0 < thxtrial)) && i == 0)
+         {
+            if (iAmRoot)
+            {
+               cout << "second order correction\n";
+            }
+            problem->c(xtrial, ckSoc);
+            problem->c(x0, ck0);
+            ckSoc.Add(alphaMax, ck0);
+            // A-5.6 Compute the second-order correction.
+            IPNewtonSolve(x0, l0, z0, zhatsoc, Xhatumlsoc, mu, true);
+            mhatsoc.Set(1.0, Xhatumlsoc.GetBlock(1));
+            //WARNING: not complete but currently solver isn't entering this region
+         }
+      }
+      else
+      {
+         if (iAmRoot)
+         {
+            cout << "in filter region :(\n"; 
+         }
+      }
+      // include more if needed
+      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;
+         }
+      }
+   }
+}
+
+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)/ double(dimCglb + dimMglb); 
+   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 && printEeval)
+   {
+      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 sqrt(InnerProduct(MPI_COMM_WORLD,cx, cx));
+}
+
+// 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, *JacuT, *JacmT;
+   Jacu = problem->Duc(x); 
+   Jacm = problem->Dmc(x);
+   JacuT = Jacu->Transpose();
+   JacmT = Jacm->Transpose();
+   
+   JacuT->Mult(l, y.GetBlock(0));
+   JacmT->Mult(l, y.GetBlock(1));
+   
+   delete JacuT;
+   delete JacmT;
+   
+   y.Add(1.0, gradxf);
+   (y.GetBlock(1)).Add(-1.0, zl);
+}
+
+bool 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::SaveLogBarrierHessianIterates(bool save)
+{
+   MFEM_ASSERT(MyRank == 0 || save == false, "currently can only save logbarrier hessian in serial codes");
+   saveLogBarrierIterates = save;
+}
+
+void ParInteriorPointSolver::SetLinearSolver(int LinSolver)
+{
+   linSolver = LinSolver;
+}
+
+void ParInteriorPointSolver::SetLinearSolveTol(double Tol)
+{
+  linSolveTol = Tol;
+}
+
+
+ParInteriorPointSolver::~ParInteriorPointSolver() 
+{
+   F1.DeleteAll();
+   F2.DeleteAll();
+   block_offsetsx.DeleteAll();
+   block_offsetsumlz.DeleteAll();
+   block_offsetsuml.DeleteAll();
+   ml.SetSize(0);
+}

examples/contact/ParIPsolver.hpp

@@ -0,0 +1,82 @@
+#include "mfem.hpp"
+#include "ParProblems.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;
+
+    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 FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
+    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , double, bool);
+    void lineSearch(BlockVector& , BlockVector& , double);
+    void projectZ(const Vector & , Vector &, double);
+    void filterCheck(double, double);
+    double E(const BlockVector &, const Vector &, const Vector &, double, bool);
+    double E(const BlockVector &, const Vector &, const Vector &, bool);
+    bool GetConverged() const;
+    // TO DO: include Hessian of Lagrangian
+    double theta(const BlockVector &);
+    double phi(const BlockVector &, double);
+    void Dxphi(const BlockVector &, double, BlockVector &);
+    double L(const BlockVector &, const Vector &, const Vector &);
+    void DxL(const BlockVector &, const Vector &, const Vector &, BlockVector &);
+    void SetTol(double);
+    void SetMaxIter(int);
+    void SetBarrierParameter(double);    
+    void SaveLogBarrierHessianIterates(bool);
+    void SetLinearSolver(int);
+    void SetLinearSolveTol(double);
+    virtual ~ParInteriorPointSolver();
+};
+
+#endif

examples/contact/ParObstacleProblem.cpp

@@ -0,0 +1,129 @@
+//                         Obstacle Problem
+//
+//
+// Compile with: make ParObstacleProblem
+//
+// Sample runs: mpirun -np 4 ./ParObstacleProblem
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve the
+//              bound-constrained energy minimization problem
+//
+//                      minimize (||∇u||² + ||u||²) subject to u ≥ ϕ in H¹.
+
+#include "mfem.hpp"
+#include "ParProblems.hpp"
+#include "ParIPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double dmanufacturedFun(const Vector &);
+double fRhs(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.Parse();
+   if(!args.Good())
+   {
+     args.PrintUsage(cout);
+     return 1;
+   }
+   else
+   {
+     if(Mpi::Root())
+     {  
+       args.PrintOptions(cout);
+     }
+   }
+
+   const char *meshFile = "../../data/inline-quad.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);
+
+   ParObstacleProblem problem(Vh,Vh,&fRhs);
+  
+   int dimD = problem.GetDimD();
+   Vector x0(dimD); x0 = 100.0;
+   Vector xf(dimD); xf = 0.0;
+
+   ParInteriorPointSolver optimizer(&problem); 
+   optimizer.SetTol(1.e-8);
+   optimizer.SetLinearSolveTol(1.e-10);
+   optimizer.SetLinearSolver(linSolver);
+   optimizer.SetMaxIter(maxIPMiters);
+   optimizer.Mult(x0, xf);
+
+   ParGridFunction d_gf(Vh);
+
+   d_gf.SetFromTrueDofs(xf);
+
+
+   FunctionCoefficient dm_fc(dmanufacturedFun); // manufactured solution
+   ParGridFunction dm_gf(Vh);
+   dm_gf.ProjectCoefficient(dm_fc);
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream exact_sock(vishost, visport);
+   exact_sock.precision(8);
+   exact_sock << "parallel " << num_procs << " " << myid << "\n";
+   exact_sock << "solution\n" << pmesh << dm_gf
+              << "window_title 'Manufactured solution'" << flush;
+
+   socketstream sol_sock(vishost, visport);
+   sol_sock.precision(8);
+   sol_sock << "parallel " << num_procs << " " << myid << "\n";
+   sol_sock << "solution\n" << pmesh << d_gf 
+            << "window_title 'Numerical solution'" << flush;
+
+   delete Vh;
+   delete fec;
+   return 0;
+}
+
+
+double dmanufacturedFun(const Vector &x)
+{
+  return cos(2*M_PI*x(0)) + 0.2 - 2.0*(pow(x(0),3) - 1.5*pow(x(0),2));
+}
+
+double fRhs(const Vector &x)
+{
+  double fx = 0.;
+  fx = 0.2 - 2.0 * (pow(x(0),3)- 1.5*pow(x(0),2.) - 6 * x(0) + 3.) + (1. + pow(2.*M_PI,2))*cos(2.*M_PI*x(0));
+  return fx;
+}

examples/contact/ParProblems.cpp

@@ -0,0 +1,212 @@
+#include "mfem.hpp"
+#include "ParProblems.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+ParGeneralOptProblem::ParGeneralOptProblem(ParFiniteElementSpace * fesU_, ParFiniteElementSpace * fesM_) 
+: fesU(fesU_), fesM(fesM_) 
+{ 
+   dimU = fesU->GetTrueVSize();
+   dimM = fesM->GetTrueVSize();
+   dimC = fesM->GetTrueVSize();
+}
+
+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(ParFiniteElementSpace * fesU_, 
+                                     ParFiniteElementSpace * fesM_)
+ : ParGeneralOptProblem(fesU_, fesM_), block_offsetsx(3)
+{
+   block_offsetsx[0] = 0;
+   block_offsetsx[1] = dimU;
+   block_offsetsx[2] = dimM;
+   block_offsetsx.PartialSum();
+   ml.SetSize(dimM); ml = 0.0;
+   Vector negIdentDiag(dimM);
+   negIdentDiag = -1.0;
+   SparseMatrix * diag = new SparseMatrix(negIdentDiag);
+   Ih = new HypreParMatrix(fesM->GetComm(), fesM->GlobalTrueVSize(),
+                                             fesM->GetTrueDofOffsets(), diag);
+   HypreStealOwnership(*Ih, *diag);
+   delete diag;
+}
+
+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 Ih;
+}
+
+
+// Obstacle Problem, no essential boundary conditions enforced
+// Hessian of energy term is K + M (stiffness + mass)
+ParObstacleProblem::ParObstacleProblem(ParFiniteElementSpace *fesU_, 
+                                       ParFiniteElementSpace *fesM_, 
+                                       double (*fSource)(const Vector &)) : 
+                                       ParOptProblem(fesU_,fesM_), f(dimU), psi(dimU), J(nullptr)
+{
+   Kform = new ParBilinearForm(fesU);
+   Kform->AddDomainIntegrator(new MassIntegrator);
+   Kform->AddDomainIntegrator(new DiffusionIntegrator);
+   Kform->Assemble();
+   Kform->Finalize();
+   Kform->FormSystemMatrix(ess_tdof_list, K);
+
+   FunctionCoefficient fcoeff(fSource);
+   fform = new ParLinearForm(fesU);
+   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
+   fform->Assemble();
+   Vector F(dimU);
+   fform->ParallelAssemble(F);
+   f.SetSize(dimU);
+   f.Set(1.0, F);
+
+   psi = 0.0;
+   
+   Vector iDiag(dimU); iDiag = 1.0;
+   SparseMatrix * Jacg = new SparseMatrix(iDiag);
+   
+   J = new HypreParMatrix(fesU->GetComm(),fesU->GlobalTrueVSize(),fesU->GetTrueDofOffsets(),Jacg);
+   HypreStealOwnership(*J, *Jacg);
+   delete Jacg;
+}
+
+// Obstacle Problem, essential boundary conditions enforced
+// Hessian of energy term is K (stiffness)
+ParObstacleProblem::ParObstacleProblem(ParFiniteElementSpace *fesU_, 
+                                       ParFiniteElementSpace *fesM_, 
+				       double (*fSource)(const Vector &),
+				       double (*obstacleSource)(const Vector &),
+				       Array<int> tdof_list, Vector &xDC) : ParOptProblem(fesU_,fesM_), f(dimU), psi(dimU), J(nullptr)
+{
+   // elastic energy functional terms	
+   ess_tdof_list = tdof_list;
+   Kform = new ParBilinearForm(fesU);
+   Kform->AddDomainIntegrator(new DiffusionIntegrator);
+   Kform->Assemble();
+   Kform->Finalize();
+   Kform->FormSystemMatrix(ess_tdof_list, K);
+
+   FunctionCoefficient fcoeff(fSource);
+   fform = new ParLinearForm(fesU);
+   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(fesU->GetComm(),fesU->GlobalTrueVSize(),fesU->GetTrueDofOffsets(),Jacg);
+   HypreStealOwnership(*J, *Jacg);
+   delete Jacg;
+
+   FunctionCoefficient psi_fc(obstacleSource);
+   ParGridFunction psi_gf(fesU);
+   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]) -= 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;
+}

examples/contact/ParProblems.hpp

@@ -0,0 +1,107 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+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
+// think about supporting general lower and upper bounds (see HiOP user manual) 
+
+class ParGeneralOptProblem
+{
+protected:
+    int dimU, dimM, dimC;
+    ParFiniteElementSpace * fesU = nullptr;
+    ParFiniteElementSpace * fesM = nullptr;
+    Array<int> block_offsetsx;
+    Vector ml;
+public:
+    ParGeneralOptProblem(ParFiniteElementSpace * fesU_, ParFiniteElementSpace * fesM_); // constructor
+    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;
+    // TO DO: include Hessian terms of constraint c
+    virtual void c(const BlockVector &, Vector &) const = 0;
+    int GetDimU() const { return dimU; };
+    int GetDimM() const { return dimM; }; 
+    int GetDimC() const { return dimC; };
+    ParFiniteElementSpace * GetfesU() {return fesU;}
+    ParFiniteElementSpace * GetfesM() {return fesM;}
+    Vector Getml() const { return ml; };
+    ~ParGeneralOptProblem(); // destructor
+};
+
+
+// abstract ContactProblem class
+// of the form
+// min_d e(d) s.t. g(d) >= 0
+class ParOptProblem : public ParGeneralOptProblem
+{
+protected:
+    Array<int> block_offsetsx;
+    HypreParMatrix * Ih;
+public:
+    ParOptProblem(ParFiniteElementSpace * fesU_, ParFiniteElementSpace * fesM_);        // constructor
+    double CalcObjective(const BlockVector &) const; // objective e
+    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 &);
+    HypreParMatrix * Duc(const BlockVector &);
+    HypreParMatrix * Dmc(const BlockVector &);
+
+    void c(const BlockVector &, Vector &) const;
+    virtual double E(const Vector &) const = 0;               // objective e(d) (energy function)
+    virtual void DdE(const Vector &, Vector &) const = 0;      // gradient of objective De / Dd
+    virtual HypreParMatrix * DddE(const Vector &) = 0;
+    // Hessian of objective D^2 e / D d^2
+    virtual HypreParMatrix * Ddg(const Vector &) = 0;
+    // Jacobian of inequality constraint Dg / Dd
+    virtual void g(const Vector &, Vector &) const = 0;       // inequality constraint g(d) >= 0 (gap function)
+    int GetDimD() const { return fesU->GetTrueVSize(); };
+    int GetDimS() const { return fesM->GetTrueVSize(); };
+    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*, ParFiniteElementSpace*, double (*fSource)(const Vector &));
+   ParObstacleProblem(ParFiniteElementSpace*, 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();
+};
+
+#endif

examples/contact/ParSphericalObstacleProblem.cpp

@@ -0,0 +1,173 @@
+//                       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 "ParProblems.hpp"
+#include "ParIPsolver.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.Parse();
+   if(!args.Good())
+   {
+     args.PrintUsage(cout);
+     return 1;
+   }
+   else
+   {
+     if(myid == 0)
+     {  
+       args.PrintOptions(cout);
+     }
+   }
+
+   const char *meshFile = "../../data/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, 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-10);
+   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);
+   }
+} 
+

examples/contact/Problems.hpp

@@ -0,0 +1,300 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <set>
+
+using namespace std;
+using namespace mfem;
+
+#ifndef PROBLEM_DEFS
+#define PROBLEM_DEFS
+
+
+
+// abstract GeneralOptProblem class
+// for the problem
+// min_(u,m) f(u,m) 
+// such that c(u,m)=0 and m >= ml
+class GeneralOptProblem
+{
+protected:
+    int dimU, dimM, dimC;
+    Array<int> block_offsetsx;
+    Vector ml;
+public:
+    GeneralOptProblem();
+    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 SparseMatrix* Duuf(const BlockVector &) = 0;
+    virtual SparseMatrix* Dumf(const BlockVector &) = 0;
+    virtual SparseMatrix* Dmuf(const BlockVector &) = 0;
+    virtual SparseMatrix* Dmmf(const BlockVector &) = 0;
+    virtual void c(const BlockVector &, Vector &) const = 0;
+    virtual SparseMatrix* Duc(const BlockVector &) = 0;
+    virtual SparseMatrix* Dmc(const BlockVector &) = 0;
+    virtual SparseMatrix* lDuuc(const BlockVector &, const Vector &) = 0;
+    virtual SparseMatrix* lDumc(const BlockVector &, const Vector &) = 0;
+    virtual SparseMatrix* lDmuc(const BlockVector &, const Vector &) = 0;
+    virtual SparseMatrix* lDmmc(const BlockVector &, const Vector &) = 0;
+    // TO DO: include log-barrier lumped-mass and pass that
+    // to the optimizer
+    //virtual SparseMatrix* GetLogBarrierLumpedMass() = 0;
+    int GetDimU() const { return dimU; };
+    int GetDimM() const { return dimM; }; 
+    int GetDimC() const { return dimC; };
+    Vector Getml() const { return ml; };
+    ~GeneralOptProblem();
+};
+
+
+// abstract OptProblem class
+// of the form
+// min_d e(d) s.t. g(d) >= 0
+class OptProblem : public GeneralOptProblem
+{
+protected:
+    int dimD;
+    int dimS;
+    Array<int> block_offsetsx;
+    SparseMatrix * negIdentity;
+    SparseMatrix * zeroMatum;
+    SparseMatrix * zeroMatmu;
+    SparseMatrix * zeroMatmm;
+public:
+    //OptProblem(int, int);        // constructor
+    OptProblem();
+    void InitializeParentData(int, int);
+    double CalcObjective(const BlockVector &) const; // objective e
+    void Duf(const BlockVector &, Vector &) const;
+    void Dmf(const BlockVector &, Vector &) const;
+    SparseMatrix* Duuf(const BlockVector &);
+    SparseMatrix* Dumf(const BlockVector &);
+    SparseMatrix* Dmuf(const BlockVector &);
+    SparseMatrix* Dmmf(const BlockVector &);
+    void c(const BlockVector &, Vector &) const;
+    SparseMatrix* Duc(const BlockVector &);
+    SparseMatrix* Dmc(const BlockVector &);
+    SparseMatrix* lDuuc(const BlockVector &, const Vector &);
+    SparseMatrix* lDumc(const BlockVector &, const Vector &);
+    SparseMatrix* lDmuc(const BlockVector &, const Vector &);
+    SparseMatrix* lDmmc(const BlockVector &, const Vector &);
+    virtual double E(const Vector &) const = 0;           // objective e(d) (energy function)
+    virtual void DdE(const Vector &, Vector &) const = 0; // gradient of objective De / Dd
+    virtual SparseMatrix* DddE(const Vector &) = 0;       // Hessian of objective D^2 e / D d^2
+    virtual void g(const Vector &, Vector &) const = 0;   // inequality constraint g(d) >= 0 (gap function)
+    virtual SparseMatrix* Ddg(const Vector &) = 0;        // Jacobian of inequality constraint Dg / Dd
+    virtual SparseMatrix* lDddg(const Vector &, const Vector &) = 0;
+    int GetDimD() const { return dimD; };
+    int GetDimS() const { return dimS; };
+    virtual ~OptProblem();
+};
+
+
+class ObstacleProblem : public OptProblem
+{
+protected:
+   // data to define energy objective function e(d) = 0.5 d^T K d - f^T d, g(d) = d + \psi >= 0
+   // stiffness matrix used to define objective
+   BilinearForm *Kform;
+   LinearForm   *fform;
+   Array<int> ess_tdof_list;
+   SparseMatrix  *K;
+   SparseMatrix  *J;
+   SparseMatrix  *Hcl;
+   FiniteElementSpace *Vh;
+   Vector f;
+   Vector psil;
+   Vector psiu;
+   bool twoBounds;
+   Vector xDC;
+   double Ce;
+public : 
+   ObstacleProblem(FiniteElementSpace*, double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &));
+   ObstacleProblem(FiniteElementSpace*, Vector&,  double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &), Array<int> tdof_list);
+   ObstacleProblem(FiniteElementSpace*, Vector &, double (*fSource)(const Vector &), double (*obstacleSourcel)(const Vector &), double (*obstacleSourceu)(const Vector &), Array<int> tdof_list);
+   double E(const Vector &) const;
+   void DdE(const Vector &, Vector &) const;
+   SparseMatrix* DddE(const Vector &);
+   void g(const Vector &, Vector &) const;
+   SparseMatrix* Ddg(const Vector &);
+   SparseMatrix * lDddg(const Vector &, const Vector &);
+   virtual ~ObstacleProblem();
+};
+
+
+SparseMatrix * GenerateProjector(int n, Array<int> ess_tdof_list);
+
+
+class ObstacleProblemVariant : public OptProblem
+{
+protected:
+   // data to define energy objective function e(d) = 0.5 d^T K d - f^T d, g(d) = d + \psi >= 0
+   // stiffness matrix used to define objective
+   BilinearForm *Kform;
+   LinearForm   *fform;
+   Array<int> ess_tdof_list;
+   Array<int> noness_tdof_list;
+   SparseMatrix  *K;
+   SparseMatrix  *RKP; // R K P = R K R^T
+   SparseMatrix  *J;
+   SparseMatrix  *Hcl;
+   SparseMatrix *R;
+   FiniteElementSpace *Vh;
+   Vector f;
+   Vector psil;
+   Vector xDC;
+   double Ce;
+public : 
+   ObstacleProblemVariant(FiniteElementSpace*, Vector&,  double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &), Array<int> tdof_list);
+   double E(const Vector &) const;
+   void DdE(const Vector &, Vector &) const;
+   SparseMatrix* DddE(const Vector &);
+   void g(const Vector &, Vector &) const;
+   SparseMatrix* Ddg(const Vector &);
+   SparseMatrix * lDddg(const Vector &, const Vector &);
+   virtual ~ObstacleProblemVariant();
+};
+
+
+
+
+class QPOptProblem : public OptProblem
+{
+protected:
+  SparseMatrix *K;
+  SparseMatrix *J;
+  SparseMatrix *zeroMatdd;
+  Vector f;
+  Vector g0;
+public:
+  QPOptProblem(const SparseMatrix, const SparseMatrix, const Vector, const Vector);
+  double E(const Vector &) const;
+  void DdE(const Vector &, Vector &) const;
+  SparseMatrix* DddE(const Vector &);
+  void g(const Vector &, Vector &) const;
+  SparseMatrix* Ddg(const Vector &);
+  SparseMatrix * lDddg(const Vector &, const Vector &);
+  virtual ~QPOptProblem();
+};
+
+
+class ExContactBlockTL : public OptProblem
+{
+public:
+   double E(const Vector &) const;
+   void DdE(const Vector &, Vector &) const;
+   SparseMatrix* DddE(const Vector &);
+   void g(const Vector &, Vector &) const;
+   SparseMatrix* Ddg(const Vector &);
+   SparseMatrix * lDddg(const Vector &, const Vector &);
+   FiniteElementSpace GetVh1();
+   FiniteElementSpace GetVh2();
+   SparseMatrix *zeroMatdd;
+public:
+   /** default constructor */
+   ExContactBlockTL(Mesh *, Mesh *, int);
+   
+
+   /** default destructor */
+   virtual ~ExContactBlockTL();
+
+private:
+   void update_g() const;
+
+private:
+   /**@name Methods to block default compiler methods.
+    *
+    * The compiler automatically generates the following three methods.
+    *  Since the default compiler implementation is generally not what
+    *  you want (for all but the most simple classes), we usually
+    *  put the declarations of these methods in the private section
+    *  and never implement them. This prevents the compiler from
+    *  implementing an incorrect "default" behavior without us
+    *  knowing. (See Scott Meyers book, "Effective C++")
+    */
+   ExContactBlockTL(
+      const ExContactBlockTL&
+   );
+
+   ExContactBlockTL& operator=(
+      const ExContactBlockTL&
+   );
+
+   Array<int> attr;
+   Array<int> m_attr;
+   Array<int> s_conn; // connectivity of the second/slave mesh
+   std::string mesh_file1;
+   std::string mesh_file2;
+   Mesh* mesh1;
+   Mesh* mesh2;
+   FiniteElementCollection* fec1;
+   FiniteElementCollection* fec2;
+   FiniteElementSpace* fespace1;
+   FiniteElementSpace* fespace2;
+   Array<int> ess_tdof_list1;
+   Array<int> ess_tdof_list2;
+   GridFunction nodes0;
+   GridFunction* nodes1;
+   GridFunction* nodes2;
+   mutable GridFunction* x1;
+   mutable GridFunction* x2;
+   PWConstCoefficient* lambda1_func;
+   PWConstCoefficient* lambda2_func;
+   PWConstCoefficient* mu1_func;
+   PWConstCoefficient* mu2_func;
+   BilinearForm* a1;
+   BilinearForm* a2;
+
+   mfem::Vector lambda1;
+   mfem::Vector lambda2;
+   mfem::Vector mu1;
+   mfem::Vector mu2;
+   mutable mfem::Vector xyz;
+
+   std::set<int> bdryVerts2;
+
+   int dim;
+   // degrees of freedom of both meshes
+   int ndof_1;
+   int ndof_2;
+   int ndofs;
+   // number of nodes for each mesh
+   int nnd_1;
+   int nnd_2;
+   int nnd;
+
+   int npoints;
+
+   SparseMatrix A1;
+   mfem::Vector B1, X1;
+   SparseMatrix A2;
+   mfem::Vector B2, X2;
+   BlockVector *B;
+   SparseMatrix* K;
+   mutable mfem::Vector gapv;
+   mutable mfem::Vector m_xi;
+   mutable mfem::Vector xs;
+
+   mutable Array<int> m_conn; // only works for linear elements that have 4 vertices!
+   mutable DenseMatrix* coordsm;
+   mutable SparseMatrix* M;
+
+   mutable std::vector<SparseMatrix>* dM;
+
+   Array<int> Dirichlet_dof;
+   Array<double> Dirichlet_val;
+   Array<int> block_offsets;
+public: 
+   Mesh * GetMesh1() {return mesh1;}
+   Mesh * GetMesh2() {return mesh2;}
+   GridFunction & GetMesh1GridFunction() {return *x1;}
+   GridFunction & GetMesh2GridFunction() {return *x2;}
+   Array<int> & GetMesh1DirichletDofs() {return ess_tdof_list1;}
+   Array<int> & GetMesh2DirichletDofs() {return ess_tdof_list2;}
+};
+
+
+
+#endif

examples/contact/SphericalObstacleProblem.cpp

@@ -0,0 +1,175 @@
+//                          Spherical Obstacle Problem
+//
+//
+// Compile with: make SphericalobstacleProblem
+//
+// Sample runs: ./SphericalobstacleProblem
+//
+//
+// 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 "Problems.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[])
+{
+  int FEorder = 1;     // finite element order
+  int linSolver = 0;   // linear solver 0 (direct), 1 (iterative) or 2 (iterative)
+  int maxIPMiters = 30;
+  bool iAmRoot = true;
+  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.Parse();
+  if(!args.Good())
+  {
+    args.PrintUsage(cout);
+    return 1;
+  }
+  else
+  {
+    if( iAmRoot )
+    {
+      args.PrintOptions(cout);
+    }
+  }
+
+  const char *meshFile = "../../data/disk.mesh";
+  Mesh *mesh = new Mesh(meshFile, 1, 1);
+  int dim = mesh->Dimension(); // geometric dimension of the domain
+  {
+     for (int l = 0; l < ref_levels; l++)
+     {
+        mesh->UniformRefinement();
+     }
+  }
+  double h_min, h_max, kappa_min, kappa_max;
+  mesh->GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
+
+  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
+  Array<int> ess_tdof_list;
+  if (mesh->bdr_attributes.Size())
+  {
+     Array<int> ess_bdr(mesh->bdr_attributes.Max());
+     ess_bdr = 1;
+     Vh->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+  }
+  
+  int dimD = Vh->GetTrueVSize();
+  Vector x0(dimD); x0 = 0.0;
+  Vector xf(dimD); xf = 0.0;
+  
+  ObstacleProblem problem(Vh, x0, &fRhs, &spherical_obstacle, ess_tdof_list);
+  InteriorPointSolver optimizer(&problem); 
+  
+  optimizer.SetTol(1.e-7);
+  optimizer.SetLinearSolver(linSolver);
+  optimizer.SetMaxIter(maxIPMiters);
+  optimizer.Mult(x0, xf);
+  
+  
+  double Einitial = problem.E(x0);
+  double Efinal = problem.E(xf);
+  cout << "Energy objective at initial point = " << Einitial << endl;
+  cout << "Energy objective at optimizer = " << Efinal << endl;
+
+  GridFunction d_gf(Vh);
+  d_gf = xf;
+
+  FunctionCoefficient dtrue_fc(exact_solution_obstacle); // exact solution
+  GridFunction dtrue_gf(Vh);
+  dtrue_gf.ProjectCoefficient(dtrue_fc);
+  
+  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
+  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("d(x) (numerical)", &d_gf);
+  paraview_dc.RegisterField("d(x) (true)", &dtrue_gf);
+  paraview_dc.Save();
+  
+  FunctionCoefficient exact_coef(exact_solution_obstacle);
+  double L2_error = d_gf.ComputeL2Error(exact_coef);
+  cout << "||u - u_true||_L^2(Omega) = " << L2_error << ", hmax = " << h_max << ", hmin = " << h_min << endl;
+  
+  delete Vh;
+  delete fec;
+  delete mesh;
+  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);
+   }
+} 
+

examples/contact/TwoSidedDirichletObstacleProblem.cpp

@@ -0,0 +1,143 @@
+#include "mfem.hpp"
+#include "Problems.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+
+double fRhs(const Vector &pt);
+double obstaclel(const Vector &pt);
+double obstacleu(const Vector &pt);
+double dmanufacturedFun(const Vector &pt);
+
+int main(int argc, char *argv[])
+{
+  int FEorder = 1; // order of the finite elements
+  int linSolver = 0;
+  int maxIPMiters = 30;
+  bool iAmRoot = true;
+  int ref_levels = 1; 
+  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.Parse();
+  if(!args.Good())
+  {
+    args.PrintUsage(cout);
+    return 1;
+  }
+  else
+  {
+    if( iAmRoot )
+    {
+      args.PrintOptions(cout);
+    }
+  }
+
+  const char *meshFile = "../../data/inline-quad.mesh";
+  Mesh *mesh = new Mesh(meshFile, 1, 1);
+  int dim = mesh->Dimension(); // geometric dimension of the domain
+  for (int l = 0; l < ref_levels; l++)
+  {
+      mesh->UniformRefinement();
+  }
+
+  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
+  Array<int> ess_tdof_list;
+  if (mesh->bdr_attributes.Size())
+  {
+     Array<int> ess_bdr(mesh->bdr_attributes.Max());
+     ess_bdr = 1;
+     Vh->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+  }
+  
+  double DC_val = 0.0;
+  int dimD = Vh->GetTrueVSize();
+
+  Vector x0(dimD); x0 = DC_val;
+  Vector xf(dimD); xf = 0.0;
+
+  ObstacleProblem problem(Vh, x0, &fRhs, &obstaclel, &obstacleu, ess_tdof_list);
+
+  InteriorPointSolver optimizer(&problem); 
+  optimizer.SetTol(1.e-7);
+  optimizer.SetLinearSolver(linSolver);
+  optimizer.SetMaxIter(maxIPMiters);
+  optimizer.Mult(x0, xf);
+  
+
+  GridFunction d_gf(Vh);
+
+  d_gf = xf;
+  
+  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
+  GridFunction dm_gf(Vh);
+  dm_gf.ProjectCoefficient(dm_fc);
+
+  
+  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
+  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("d(x) (numerical)", &d_gf);
+  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
+  paraview_dc.Save();
+  
+
+  delete Vh;
+  delete fec;
+  delete mesh;
+  return 0;
+}
+
+
+
+double dmanufacturedFun(const Vector &pt)
+{
+  double alpha = 16.5;
+  return sin(M_PI * pt(1)) * (sin(M_PI * pt(0)) - alpha * pow(pt(0) * (1. - pt(0)), 2));
+}
+
+
+// f(x) forcing term... which enters the objective energy functional
+// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
+// of f(x). f(x) is such that in the absence of bound-constraints then
+// the solution of the optimization problem satisfies the PDE
+// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
+// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
+
+double fRhs(const Vector &pt)
+{
+  double alpha = 16.5;
+  double fx;
+  fx = pow(M_PI, 2) * sin(M_PI * pt(0));
+  fx += alpha * (2. * pow(pt(0), 2) + 2. * pow(1.-pt(0), 2) - 8. * pt(0) * (1.-pt(0)));
+  fx += pow(M_PI, 2) * sin(M_PI * pt(0)) * dmanufacturedFun(pt);
+  fx *= sin(M_PI * pt(1));
+  return fx;
+}
+
+double obstaclel(const Vector &pt)
+{
+  return 0.0;
+}
+
+double obstacleu(const Vector &pt)
+{
+  return 0.08;
+}

examples/contact/block1.mesh

@@ -0,0 +1,103 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+3
+
+elements
+9
+1   5     0      1     3     2     8    9     11    10
+1   5     2      3     5     4    10    11    13    12
+1   5     4      5     7     6    12    13    15    14
+1   5     8      9    11    10    16    17    19    18
+1   5     10    11    13    12    18    19    21    20
+1   5     12    13    15    14    20    21    23    22
+1   5     16    17    19    18    24    25    27    26
+1   5     18    19    21    20    26    27    29    28
+1   5     20    21    23    22    28    29    31    30
+
+
+
+# 0 nothing
+# 1 dirichlet bc
+# 2 contact
+boundary
+30
+1   3    1     0     2     3
+1   3    3     2     4     5
+1   3    5     4     6     7
+1   3    24    25    27    26
+1   3    26    27    29    28
+1   3    28    29    31    30
+2   3    2     0     8    10
+2   3    4     2    10    12
+2   3    6     4    12    14
+2  3    10    8    16    18
+2  3    12    10    18    20
+2  3    14    12    20    22
+2  3    18    16    24    26
+2  3    20    18    26    28
+2  3    22    20    28    30
+3  3    1     3    11    9
+3  3     3     5    13    11
+3  3     5     7    15    13
+3  3    9    11    19    17
+3  3    11    13    21    19
+3  3    13    15    23    21
+3  3    17    19    27    25
+3  3    19    21    29    27
+3  3    21    23    31    29
+1  3     8     0     1    9
+1  3    16     8    9    17
+1  3    24    16    17    25
+1  3     6    14    15     7
+1  3    14    22    23    15
+1  3    22    30    31    23
+
+
+vertices
+32
+3
+   -1.0000         0         0
+         0         0         0
+   -1.0000    0.3000         0
+         0    0.3000         0
+   -1.0000    0.6500         0
+         0    0.6500         0
+   -1.0000    1.0000         0
+         0    1.0000         0
+   -1.0000         0    0.3000
+         0         0    0.3000
+   -1.0000    0.3000    0.3500
+         0    0.3000    0.3500
+   -1.0000    0.6500    0.3000
+         0    0.6500    0.3000
+   -1.0000    1.0000    0.3000
+         0    1.0000    0.3000
+   -1.0000         0    0.6500
+         0         0    0.6500
+   -1.0000    0.3000    0.6500
+         0    0.3000    0.6500
+   -1.0000    0.6500    0.6500
+         0    0.6500    0.6500
+   -1.0000    1.0000    0.6500
+         0    1.0000    0.6500
+   -1.0000         0    1.0000
+         0         0    1.0000
+   -1.0000    0.3000    1.0000
+         0    0.3000    1.0000
+   -1.0000    0.6500    1.0000
+         0    0.6500    1.0000
+   -1.0000    1.0000    1.0000
+         0    1.0000    1.0000

examples/contact/exQPContactBlockTL.cpp

@@ -0,0 +1,246 @@
+//                                Quadratic-Programming (QP) Contact example
+//
+// Compile with: make exQPContactBlockTL
+//
+// Sample runs:  ./exQPContactBlockTL
+
+
+#include <fstream>
+#include <iostream>
+#include <array>
+
+#include "mfem.hpp"
+#include "Problems.hpp"
+#include "IPsolver.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+  int linSolver = 0;
+  int maxIPMiters = 30;
+  bool iAmRoot = true;
+  int ref_levels = 0;
+
+  OptionsParser args(argc, argv);
+  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.Parse();
+  if(!args.Good())
+  {
+    args.PrintUsage(cout);
+    return 1;
+  }
+  else
+  {
+    if( iAmRoot )
+    {
+      args.PrintOptions(cout);
+    }
+  }
+	
+  Mesh * mesh1 = new Mesh("block1.mesh", 1, 1);
+  Mesh * mesh2 = new Mesh("rotatedblock2.mesh", 1, 1);
+  for(int i = 0; i < ref_levels; i++)
+  {
+     mesh1->UniformRefinement();
+     mesh2->UniformRefinement(); 
+  }
+
+  // Create an instance of the nlp
+  ExContactBlockTL * contact = new ExContactBlockTL(mesh1, mesh2, 1);
+  int ndofs = contact->GetDimD();
+  int nconstraints = contact->GetDimS();
+  
+  // set up a QP-problem
+  // E(d) = 1 / 2 d^T K d + f^T d
+  // g(d) = J d + g0
+  // where K, J, f and g0 are evaluated at d0 (a valid configuration)
+  
+  // to do: seems more appropriate to evaluate at a valid configuration...
+  //        that is one where the Dirichlet conditions hold... need to pull
+  //        this data from contactBlockTL...
+  Vector d0(ndofs); d0 = 0.0;
+  Array<int> ess_tdofs1 = contact->GetMesh1DirichletDofs();
+  Array<int> ess_tdofs2 = contact->GetMesh2DirichletDofs();
+  int sz1 = ess_tdofs1.Size();
+  int sz2 = ess_tdofs2.Size();
+  Array<int> DirichletDofs(sz1+sz2);
+  for (int i = 0; i<sz1; i++)
+  {
+    DirichletDofs[i] = ess_tdofs1[i];
+  }
+  for (int i = 0; i<sz2; i++)
+  {
+    DirichletDofs[i+sz1] = ess_tdofs2[i]+contact->GetVh1().GetTrueVSize();
+  }
+  GridFunction x1 = contact->GetMesh1GridFunction();
+  GridFunction x2 = contact->GetMesh2GridFunction();
+
+  SparseMatrix *K;
+  Vector f(ndofs); f = 0.0;
+  contact->DdE(d0, f); K = contact->DddE(d0);
+  d0.SetVector(x1,0);
+  d0.SetVector(x2,x1.Size());
+  SparseMatrix *J;
+  Vector g0(nconstraints); g0 = 0.0;
+  contact->g(d0, g0); J = contact->Ddg(d0);
+  Vector temp(nconstraints);
+  J->Mult(d0, temp);
+  g0.Add(-1.0, temp);
+
+  // check which rows of the Jacobian are zero!
+  Vector ei(nconstraints); ei = 0.0;
+  Vector JTei(ndofs); JTei = 0.0;
+
+  double normJTei;
+
+  Array<int> nonZeroRows;
+  for(int i = 0; i < nconstraints; i++)
+  {
+    Array<int> col_tmp;
+    Vector v_tmp; v_tmp = 0.0;
+    J->GetRow(i, col_tmp, v_tmp);
+    normJTei = v_tmp.Norml2();
+    if (normJTei > 1.e-12)
+    {
+      nonZeroRows.Append(i);
+    }
+  }
+  mfem::out << J->Height() << " linearized constraints\n";
+  mfem::out << nonZeroRows.Size() << " (reduced) linearized constraints\n"; 
+  
+  // remove zero rows of the gap function Jacobian and corresponding gap function entries
+  SparseMatrix * Jreduced = new SparseMatrix(nonZeroRows.Size(), ndofs);
+  Vector g0reduced(nonZeroRows.Size()); g0reduced = 0.0; 
+  
+  
+  for(int i = 0; i < nonZeroRows.Size(); i++)
+  {
+    Array<int> col_tmp;
+    Vector v_tmp; v_tmp = 0.0;
+    J->GetRow(nonZeroRows[i], col_tmp, v_tmp);
+    
+    /* obtain subset of columns of the given nonZero Jacobian row that are not Dirichlet constrained */
+    bool freeDof;
+    Array<int> free_col_indicies;
+    for(int j = 0; j < col_tmp.Size(); j++)
+    {
+      freeDof = true;
+      for(int k = 0; k < DirichletDofs.Size(); k++)
+      {
+        if(col_tmp[j] == DirichletDofs[k])
+        {
+          freeDof = false;
+        }
+      }
+      if(freeDof)
+      {
+        free_col_indicies.Append(j);
+      }
+    }
+
+    Array<int> col_tmp_reduced(free_col_indicies.Size());
+    Vector v_tmp_reduced(free_col_indicies.Size());
+    for(int j = 0; j < free_col_indicies.Size(); j++)
+    {
+      col_tmp_reduced[j] = col_tmp[free_col_indicies[j]];
+      v_tmp_reduced(j)   = v_tmp(free_col_indicies[j]);
+    }
+
+    Jreduced->SetRow(i, col_tmp_reduced, v_tmp_reduced);
+    g0reduced(i) = g0(nonZeroRows[i]);
+  }
+
+  QPOptProblem *QPContact = new QPOptProblem(*K, *Jreduced, f, g0reduced);
+  
+  InteriorPointSolver * QPContactOptimizer = new InteriorPointSolver(QPContact);
+  QPContactOptimizer->SetTol(1.e-6);
+  QPContactOptimizer->SetLinearSolver(linSolver);
+  Vector x0(ndofs); x0 = 0.0;
+  x0.SetVector(x1,0);
+  x0.SetVector(x2,x1.Size());
+
+  Vector xf(ndofs); xf = 0.0;
+  QPContactOptimizer->Mult(x0, xf);
+  
+  MFEM_VERIFY(QPContactOptimizer->GetConverged(), "Interior point solver did not converge.");
+  double Einitial = QPContact->E(x0);
+  double Efinal = QPContact->E(xf);
+  cout << "Energy objective at initial point = " << Einitial << endl;
+  cout << "Energy objective at QP optimizer = " << Efinal << endl;
+ 
+  
+  
+  int gdim = mesh1->Dimension();
+  FiniteElementCollection * fec = new H1_FECollection(1, gdim);
+  FiniteElementSpace * fespace1 = new FiniteElementSpace(mesh1, fec, gdim, Ordering::byVDIM);
+  FiniteElementSpace * fespace2 = new FiniteElementSpace(mesh2, fec, gdim, Ordering::byVDIM);
+   
+  GridFunction x1_gf(fespace1);
+  GridFunction x2_gf(fespace2);
+
+  int ndof1 = fespace1->GetTrueVSize();
+  int ndof2 = fespace2->GetTrueVSize();
+  int ndof  = ndof1 + ndof2;
+  for(int i = 0; i < ndof1; i++)
+  {
+    x1_gf(i) = xf(i);
+  }
+  for(int i = ndof1; i < ndof; i++)
+  {
+    x2_gf(i - ndof1) = xf(i);
+  }
+
+  mesh1->SetNodalFESpace(fespace1);
+  mesh2->SetNodalFESpace(fespace2);
+  GridFunction *nodes1 = mesh1->GetNodes();
+  GridFunction *nodes2 = mesh2->GetNodes();
+
+  {
+    *nodes1 += x1_gf;
+    *nodes2 += x2_gf;
+  }
+  
+
+  ParaViewDataCollection paraview_dc1("QPContactBody1", mesh1);
+  paraview_dc1.SetPrefixPath("ParaView");
+  paraview_dc1.SetLevelsOfDetail(1);
+  paraview_dc1.SetDataFormat(VTKFormat::BINARY);
+  paraview_dc1.SetHighOrderOutput(true);
+  paraview_dc1.SetCycle(0);
+  paraview_dc1.SetTime(0.0);
+  paraview_dc1.RegisterField("Body1", &x1_gf);
+  paraview_dc1.Save();
+  
+  ParaViewDataCollection paraview_dc2("QPContactBody2", mesh2);
+  paraview_dc2.SetPrefixPath("ParaView");
+  paraview_dc2.SetLevelsOfDetail(1);
+  paraview_dc2.SetDataFormat(VTKFormat::BINARY);
+  paraview_dc2.SetHighOrderOutput(true);
+  paraview_dc2.SetCycle(0);
+  paraview_dc2.SetTime(0.0);
+  paraview_dc2.RegisterField("Body2", &x2_gf);
+  paraview_dc2.Save();
+
+  delete fespace1;
+  delete fespace2;
+  delete fec;
+  delete mesh1;
+  delete mesh2;
+  
+  delete QPContact;
+  delete QPContactOptimizer;
+  
+  delete Jreduced;
+  delete contact;
+  return 0;
+}

examples/contact/makefile

@@ -0,0 +1,125 @@
+# 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)/examples/contact/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ObstacleProblem SphericalObstacleProblem DirichletObstacleProblem exQPContactBlockTL
+PAR_EXAMPLES = ParObstacleProblem
+EXAMPLES = $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+
+ifeq ($(MFEM_USE_SUITESPARSE),NO)
+$(SEQ_EXAMPLES):
+	$(error MFEM is not configured with SUITESPARSE)
+endif
+
+ifeq ($(MFEM_USE_MUMPS),NO)
+ifeq ($(MFEM_USE_MKL_CPARDISO), NO)
+$(PAR_EXAMPLES):
+	$(error MFEM is not configured with MUMPS or CPARDISO)
+endif
+endif
+
+all: $(EXAMPLES)
+
+ObstacleProblem: ObstacleProblem.o Problems.o IPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) ObstacleProblem.o Problems.o IPsolver.o -o $@ $(MFEM_LIBS)
+
+SphericalObstacleProblem: SphericalObstacleProblem.o Problems.o IPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) SphericalObstacleProblem.o Problems.o IPsolver.o -o $@ $(MFEM_LIBS)
+
+DirichletObstacleProblem: DirichletObstacleProblem.o Problems.o IPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) DirichletObstacleProblem.o Problems.o IPsolver.o -o $@ $(MFEM_LIBS)
+
+DirichletObstacleProblemVariant: DirichletObstacleProblemVariant.o Problems.o IPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) DirichletObstacleProblemVariant.o Problems.o IPsolver.o -o $@ $(MFEM_LIBS)
+
+TwoSidedDirichletObstacleProblem: TwoSidedDirichletObstacleProblem.o Problems.o IPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) TwoSidedDirichletObstacleProblem.o Problems.o IPsolver.o -o $@ $(MFEM_LIBS)
+
+
+exQPContactBlockTL: exQPContactBlockTL.o Problems.o IPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) exQPContactBlockTL.o Problems.o IPsolver.o -o $@ $(MFEM_LIBS)
+
+ParTest: ParTest.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) ParTest.o -o $@ $(MFEM_LIBS)
+
+
+ObstacleProblem.o: $(SRC)ObstacleProblem.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< 
+
+SphericalObstacleProblem.o: $(SRC)SphericalObstacleProblem.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< 
+
+DirichletObstacleProblem.o: $(SRC)DirichletObstacleProblem.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+
+DirichletObstacleProblemVariant.o: $(SRC)DirichletObstacleProblemVariant.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+
+TwoSidedDirichletObstacleProblem.o: $(SRC)TwoSidedDirichletObstacleProblem.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+
+exQPContactBlockTL.o: $(SRC)exQPContactBlockTL.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+
+Problems.o: $(SRC)Problems.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< 
+
+IPsolver.o: $(SRC)IPsolver.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+
+ParObstacleProblem: ParObstacleProblem.o ParProblems.o ParIPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) ParObstacleProblem.o ParProblems.o ParIPsolver.o -o $@ $(MFEM_LIBS)
+
+ParSphericalObstacleProblem: ParSphericalObstacleProblem.o ParProblems.o ParIPsolver.o $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) ParSphericalObstacleProblem.o ParProblems.o ParIPsolver.o -o $@ $(MFEM_LIBS)
+
+ParObstacleProblem.o: $(SRC)ParObstacleProblem.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< 
+
+ParSphericalObstacleProblem.o: $(SRC)ParSphericalObstacleProblem.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< 
+
+ParProblems.o: $(SRC)ParProblems.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< 
+
+ParIPsolver.o: $(SRC)ParIPsolver.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+
+ParTest.o: $(SRC)ParTest.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
+
+clean: clean-build 
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+# For out-of-source builds, link the data files from the source tree:
+ifneq ($(SRC),)
+DATA_FILES = block1.mesh rotatedblock2.mesh
+$(DATA_FILES): %: $(SRC)%
+	ln -sf $(<) .
+copy-data: | $(DATA_FILES)
+# For out-of-source builds, the test and sample runs for 'field-interp' need
+# data from the meshing miniapps directory:
+exQPContactBlockTL: | mesh-data
+.PHONY: mesh-data
+mesh-data:
+	$(MAKE) -C ./ copy-data
+endif

examples/contact/nodepair.hpp

@@ -0,0 +1,896 @@
+
+using namespace std;
+using namespace mfem;
+
+void BasisEval(const Vector xi, Vector &N, DenseMatrix &dNdxi) // dNdxi is 2*4
+{
+   N[0] = 0.25*(1-xi[0])*(1-xi[1]);
+   N[1] = 0.25*(1+xi[0])*(1-xi[1]);
+   N[2] = 0.25*(1+xi[0])*(1+xi[1]);
+   N[3] = 0.25*(1-xi[0])*(1+xi[1]);
+
+   dNdxi(0,0) = 0.25*(-1+xi[1]);
+   dNdxi(0,1) = 0.25*(1-xi[1]);
+   dNdxi(0,2) = 0.25*(1+xi[1]);
+   dNdxi(0,3) = 0.25*(-1-xi[1]);
+   dNdxi(1,0) = 0.25*(-1+xi[0]);
+   dNdxi(1,1) = 0.25*(-1-xi[0]);
+   dNdxi(1,2) = 0.25*(1+xi[0]);
+   dNdxi(1,3) = 0.25*(1-xi[0]);
+}
+
+
+void BasisEvalDerivs(const Vector xi, Vector& N, DenseMatrix& dNdxi,
+                     DenseMatrix& dN2dxi)
+{
+   N[0] = 0.25*(1-xi[0])*(1-xi[1]);
+   N[1] = 0.25*(1+xi[0])*(1-xi[1]);
+   N[2] = 0.25*(1+xi[0])*(1+xi[1]);
+   N[3] = 0.25*(1-xi[0])*(1+xi[1]);
+
+   dNdxi.SetSize(2,4); dNdxi = 0.0;
+   dN2dxi.SetSize(3,4);
+   dN2dxi = 0.0;  // first row dxi2, second detadxi, third deta2
+
+   dNdxi(0,0) = 0.25*(-1+xi[1]); dNdxi(0,1) = 0.25*(1-xi[1]);
+   dNdxi(0,2) = 0.25*(1+xi[1]);  dNdxi(0,3) = 0.25*(-1-xi[1]);
+   dNdxi(1,0) = 0.25*(-1+xi[0]); dNdxi(1,1) = 0.25*(-1-xi[0]);
+   dNdxi(1,2) = 0.25*(1+xi[0]);  dNdxi(1,3) = 0.25*(1-xi[0]);
+
+   dN2dxi(1,0) = 0.25; dN2dxi(1,1) = -0.25; dN2dxi(1,2) = 0.25;
+   dN2dxi(1,3) = -0.25;
+}
+
+// returns the vector and matrix form of the shape functions and its derivative
+void BasisVectorDerivs(const Vector xi, DenseMatrix& N, DenseMatrix& dNdxi,
+                       DenseMatrix& ddNdxi)
+{
+   N.SetSize(3,12); N = 0.0;
+   N(0,0) = 0.25*(1-xi[0])*(1-xi[1]); N(0,3) = 0.25*(1+xi[0])*(1-xi[1]);
+   N(0,6) = 0.25*(1+xi[0])*(1+xi[1]); N(0,9) = 0.25*(1-xi[0])*(1+xi[1]);
+
+   N(1,1) = 0.25*(1-xi[0])*(1-xi[1]); N(1,4) = 0.25*(1+xi[0])*(1-xi[1]);
+   N(1,7) = 0.25*(1+xi[0])*(1+xi[1]); N(1,10) = 0.25*(1-xi[0])*(1+xi[1]);
+
+   N(2,2) = 0.25*(1-xi[0])*(1-xi[1]); N(2,5) = 0.25*(1+xi[0])*(1-xi[1]);
+   N(2,8) = 0.25*(1+xi[0])*(1+xi[1]); N(2,11) = 0.25*(1-xi[0])*(1+xi[1]);
+
+   dNdxi.SetSize(3*2, 3*4);  dNdxi = 0.0;
+   dNdxi(0,0) = 0.25*(-1+xi[1]);  dNdxi(0,3) = 0.25*(1-xi[1]);
+   dNdxi(0,6) = 0.25*(1+xi[1]);  dNdxi(0,9) = 0.25*(-1-xi[1]);
+   dNdxi(1,1) = 0.25*(-1+xi[1]);  dNdxi(1,4) = 0.25*(1-xi[1]);
+   dNdxi(1,7) = 0.25*(1+xi[1]);  dNdxi(1,10) = 0.25*(-1-xi[1]);
+   dNdxi(2,2) = 0.25*(-1+xi[1]);  dNdxi(2,5) = 0.25*(1-xi[1]);
+   dNdxi(2,8) = 0.25*(1+xi[1]);  dNdxi(2,11) = 0.25*(-1-xi[1]);
+
+   dNdxi(3,0) = 0.25*(-1+xi[0]);  dNdxi(3,3) = 0.25*(-1-xi[0]);
+   dNdxi(3,6) = 0.25*(1+xi[0]);  dNdxi(3,9) = 0.25*(1-xi[0]);
+   dNdxi(4,1) = 0.25*(-1+xi[0]);  dNdxi(4,4) = 0.25*(-1-xi[0]);
+   dNdxi(4,7) = 0.25*(1+xi[0]);  dNdxi(4,10) = 0.25*(1-xi[0]);
+   dNdxi(5,2) = 0.25*(-1+xi[0]);  dNdxi(5,5) = 0.25*(-1-xi[0]);
+   dNdxi(5,8) = 0.25*(1+xi[0]);  dNdxi(5,11) = 0.25*(1-xi[0]);
+
+   ddNdxi.SetSize(3*4, 3*4); ddNdxi = 0.0;
+   ddNdxi(3,0) = 0.25; ddNdxi(3,3) = -0.25;
+   ddNdxi(3,6) = 0.25; ddNdxi(3,9) = -0.25;
+   ddNdxi(4,1) = 0.25; ddNdxi(4,4) = -0.25;
+   ddNdxi(4,7) = 0.25; ddNdxi(4,10) = -0.25;
+   ddNdxi(5,2) = 0.25; ddNdxi(5,5) = -0.25;
+   ddNdxi(5,8) = 0.25; ddNdxi(5,11) = -0.25;
+
+   ddNdxi(6,0) = 0.25; ddNdxi(6,3) = -0.25;
+   ddNdxi(6,6) = 0.25; ddNdxi(6,9) = -0.25;
+   ddNdxi(7,1) = 0.25; ddNdxi(7,4) = -0.25;
+   ddNdxi(7,7) = 0.25; ddNdxi(7,10) = -0.25;
+   ddNdxi(8,2) = 0.25; ddNdxi(8,5) = -0.25;
+   ddNdxi(8,8) = 0.25; ddNdxi(8,11) = -0.25;
+}
+
+
+void cross(const Vector a, const Vector b, Vector& c)
+{
+   assert(a.Size()==3);
+   c.SetSize(3);
+   c[0] = a[1]*b[2] - a[2]*b[1];
+   c[1] = -a[0]*b[2] + b[0]*a[2];
+   c[2] = a[0]*b[1] - a[1]*b[0];
+
+}
+// a outer b
+void outer(const Vector a, const Vector b, DenseMatrix& c)
+{
+   int m = a.Size();
+   int n = b.Size();
+   assert(c.Height()==m);
+   assert(c.Width() ==n);
+   for (int i=0; i<m; i++)
+   {
+      for (int j=0; j<n; j++)
+      {
+         c(i,j) = a[i]*b[j];
+      }
+   }
+}
+// dphidxi 2*4
+// coords 4*3
+void ComputeNormal(const DenseMatrix& dphidxi, const DenseMatrix& coords,
+                   Vector& normal, double& nnorm)
+{
+
+   DenseMatrix dxdxi(2,3);
+   Mult(dphidxi, coords, dxdxi);
+   Vector dxdxi1(3);
+   Vector dxdxi2(3);
+
+   dxdxi.GetRow(0,dxdxi1);
+   dxdxi.GetRow(1,dxdxi2);
+
+   cross(dxdxi1, dxdxi2, normal);  // is there a cross product?  no
+   // VectorCrossProductCoefficient::Eval has hard-coded cross product
+   nnorm = normal.Norml2( );
+   normal /=  nnorm;
+}
+
+void SlaveToMaster(const DenseMatrix& m_coords, const Vector& s_x, Vector& xi)
+{
+   bool converged = false;
+   bool pt_on_elem = false;
+   int dim = 3;
+   xi.SetSize(dim-1);
+   xi = 0.0;
+   int max_iter = 15;
+   double off_el_xi = 1e-2;
+   double proj_newton_tol = 1e-13;
+   double proj_max_gap = 0.5;
+   Vector gap_v(dim);
+   // warm start from linear solution
+
+   for (int it=0; it<max_iter; it++)
+   {
+      //cout<<it<<endl;
+      Vector m_N(4);
+      m_N = 0.;
+      DenseMatrix m_dN(2,4);
+      m_dN = 0.;
+      DenseMatrix m_dN2(3,4);
+      m_dN2 = 0.;
+      BasisEvalDerivs(xi, m_N, m_dN, m_dN2);
+
+      Vector x_c(dim);
+      m_coords.MultTranspose(m_N, x_c);
+
+      gap_v = s_x;
+      gap_v -= x_c;
+
+      DenseMatrix m_dx(2,3);
+      m_dx = 0.;
+      Mult(m_dN, m_coords, m_dx);
+
+      Vector r(dim-1);
+      r = 0.0;
+      m_dx.Mult(gap_v, r);
+
+      if (r.Normlinf() < proj_newton_tol)
+      {
+         converged = true;
+         break;
+      }
+
+      DenseMatrix drdxi(dim-1,dim-1);
+      drdxi = 0.;
+      MultABt(m_dx, m_dx, drdxi); // m_dx * m_dx.T
+      drdxi *= -1.0;
+
+      DenseMatrix m_dx2(3,3); m_dx2 = 0.0;
+      Mult(m_dN2,m_coords, m_dx2);
+
+      //m_d2x = m_dN(:,:,2) * m_elem_coords(1:4,:);  //m_dN(:,:,2) is 3*4
+      for (int d=0; d<3; d++)
+      {
+         DenseMatrix Mtemp(2,2); Mtemp = 0.0;
+         Mtemp(0,0) = m_dx2(0,d); Mtemp(0,1) = m_dx2(1,d);
+         Mtemp(1,0) = m_dx2(1,d); Mtemp(1,1) = m_dx2(2,d);
+
+         drdxi.Add(gap_v[d], Mtemp);
+      }
+
+      //cond_num = rcond(drdxi);  condition number?
+      //drdxi.TestInversion();
+      DenseMatrixInverse drdxi_inv(drdxi);
+      Vector xi_tmp(dim-1);
+
+      drdxi_inv.Mult(r,xi_tmp);
+      xi -= xi_tmp;
+   }
+   if (!converged)
+   {
+      xi = 0.0;
+   }
+   off_el_xi += 1 ; // tolerance of offset of xi outside [-1,1]
+
+   //cout<<gap_v.Norml2()<<" " <<xi.Normlinf()<<endl;
+   //
+   // Discuss with Frank... what is happening here
+   if (gap_v.Norml2() < proj_max_gap && xi.Normlinf() <= off_el_xi)
+   {
+      pt_on_elem = true;
+   }
+
+   if (pt_on_elem)
+   {
+     //cout << "convergence of node to segment projection? " << converged << endl;
+     //for(int i = 0; i < 2; i++)
+     //{
+     //  cout << "xi_" << i << " = " << xi(i) << endl;
+     //}
+   }
+   MFEM_VERIFY(pt_on_elem == true, "xi went out of bounds");
+   MFEM_VERIFY(converged == true, "projection didn't converge");
+}
+
+
+
+// m_coords is expected to be 4 * 3
+void  ComputeGapJacobian(const Vector x_s, const Vector xi,
+                         const DenseMatrix m_coords,
+                         double& gap, Vector& normal, Vector& dgdxm, Vector& dgdxs)
+{
+   Vector m_N(4);
+   DenseMatrix m_dN(2,4);
+   DenseMatrix m_dN2(3,4);
+   BasisEvalDerivs(xi, m_N, m_dN, m_dN2);
+
+   Vector x_c(3);
+   m_coords.MultTranspose(m_N, x_c);
+
+   Vector gap_v(3); gap_v = 0.0;
+   gap_v  = x_s;
+   gap_v -= x_c;
+
+   DenseMatrix m_dx(2,3);
+   Mult(m_dN, m_coords, m_dx);
+
+   double nnorm = 0;
+   ComputeNormal(m_dN, m_coords, normal, nnorm);
+
+   gap = gap_v * normal; // gap function value, dot product between vectors
+
+   //dr_dx = zeros(2,4,3); % nsegment, nodes in quad, ndim
+
+   DenseMatrix dr_dx_res1(4,3); dr_dx_res1 = 0.;
+   DenseMatrix dr_dx_res2(4,3); dr_dx_res2 = 0.;
+
+   Vector m_dxrow1(3);
+   m_dx.GetRow(0, m_dxrow1);
+   outer(m_N, m_dxrow1, dr_dx_res1);// 4*1 times 1*3
+   dr_dx_res1 *= -1.0;
+
+   Vector m_dxrow2(3);
+   m_dx.GetRow(1, m_dxrow2);
+   outer(m_N, m_dxrow2, dr_dx_res2);// 4*1 times 1*3
+   dr_dx_res2 *= -1.0;
+
+   Vector m_dNrow1(4); m_dN.GetRow(0, m_dNrow1);
+   Vector m_dNrow2(4); m_dN.GetRow(1, m_dNrow2);
+
+   DenseMatrix dr_dx_res1_tmp(4,3); dr_dx_res1_tmp = 0.;
+   DenseMatrix dr_dx_res2_tmp(4,3); dr_dx_res2_tmp = 0.;
+   outer(m_dNrow1, gap_v, dr_dx_res1_tmp);// 4*1 times 1*3
+   outer(m_dNrow2, gap_v, dr_dx_res2_tmp);// 4*1 times 1*3
+
+   dr_dx_res1 += dr_dx_res1_tmp; // outer product in vector?
+   dr_dx_res2 += dr_dx_res2_tmp;
+
+
+   DenseMatrix K_dxidx1(2,2); // 2*2
+   K_dxidx1 = 0.;
+   MultABt(m_dx, m_dx, K_dxidx1); // m_dx * m_dx.T
+
+   Vector v_dxidx2(4);
+   m_coords.Mult(gap_v, v_dxidx2); // m_coords * gap_v;  // 4*3 * 3 = 4
+
+   DenseMatrix K_dxidx2(2,2); K_dxidx2 = 0.0;
+
+   Vector m_dN2row1(4); m_dN2.GetRow(0, m_dN2row1);
+   Vector m_dN2row2(4); m_dN2.GetRow(1, m_dN2row2);
+   Vector m_dN2row3(4); m_dN2.GetRow(2, m_dN2row3);
+   // how to get 2nd order? multidimensional matrix?
+   K_dxidx2(0,0) = m_dN2row1 * v_dxidx2; // how would 4*1 * 1*4 be computed?
+   K_dxidx2(0,1) = m_dN2row2 * v_dxidx2;
+   K_dxidx2(1,0) = m_dN2row2 * v_dxidx2;
+   K_dxidx2(1,1) = m_dN2row3 * v_dxidx2;
+
+   DenseMatrix K_dxidx(2,2);
+   K_dxidx -= K_dxidx1;
+   K_dxidx += K_dxidx2;
+
+   // resize the vectors and matrices
+   Vector dxidx(24); dxidx = 0.0;
+   Vector drdx_r(24); drdx_r = 0.0;
+
+   for (int i=0; i<4; i++)
+   {
+      for (int j=0; j<3; j++)
+      {
+         drdx_r[4*j+i] = dr_dx_res1(i,j);
+         drdx_r[4*j+i+12] = dr_dx_res2(i,j);
+
+      }
+   }
+   //drdx_r(1:4*3,1) = reshape(dr_dx_res(:,:,1),4*3,1);
+   //drdx_r(4*3+1:2*4*3,1) = reshape(dr_dx_res(:,:,2),4*3,1);
+   DenseMatrix drdx_K(24,24); drdx_K = 0.;
+   for (int i =0; i<12; i++)
+   {
+      drdx_K(i,i) = K_dxidx(0,0);
+      drdx_K(i,12+i) = K_dxidx(0,1);
+      drdx_K(12+i,i) = K_dxidx(1,0);
+      drdx_K(12+i,12+i) = K_dxidx(1,1);
+   }
+
+   DenseMatrixInverse drdxK_inv(drdx_K);
+   drdxK_inv.Mult(drdx_r,dxidx);
+   // LinearSolve (drdx_K,drdx_r, dxidx) ; //???
+   dxidx *= -1.0;
+
+
+
+   Vector drdxs_r(6);
+   drdxs_r[0] = m_dx(0,0); drdxs_r[1] = m_dx(0,1); drdxs_r[2] = m_dx(0,2);
+   drdxs_r[3] = m_dx(1,0); drdxs_r[4] = m_dx(1,1); drdxs_r[5] = m_dx(1,2);
+
+   DenseMatrix drdxs_K(6,6); drdxs_K = 0.;
+   for (int i=0; i<3; i++)
+   {
+      drdxs_K(i,i) = K_dxidx(0,0);
+      drdxs_K(i,3+i) = K_dxidx(0,1);
+      drdxs_K(i+3,i) = K_dxidx(1,0);
+      drdxs_K(i+3,i+3) = K_dxidx(1,1);
+   }
+
+   Vector dxidxs(6); dxidxs = 0.0;
+   DenseMatrixInverse drdxsK_inv(drdxs_K);
+   drdxsK_inv.Mult(drdxs_r,dxidxs);
+   dxidxs *= -1.0;
+   //dxidxs = -drdxs_K\drdxs_r;
+
+   //dxidx = reshape(dxidx, 4,3,2); dxidxs = reshape(dxidxs, 1,3,2);
+
+   dgdxm.SetSize(12); dgdxm = 0.;
+   DenseMatrix dgdxm_tmp(4,3);
+   outer(m_N, normal,dgdxm_tmp);
+   for (int i=0; i<4; i++)
+   {
+      for (int j=0; j<3; j++)
+      {
+         dgdxm[3*i+j] = -dgdxm_tmp(i,j);
+      }
+   }
+   //dxidx_M = -m_dN(1:2,:,1) * (m_coords(1:4,:)*normal'); % this turns out to be 0
+
+   dgdxs.SetSize(3);
+   dgdxs += normal;
+   //dgdxs = dgdxs + dxidx_M(1) * dxidxs(:,:,1) + dxidx_M(2) * dxidxs(:,:,2);
+};
+
+void ComputeGapHessian(const Vector x_s, const Vector xi,
+                       const DenseMatrix m_coords,
+                       DenseMatrix& dg2dx)
+{
+   Vector m_N(4);
+   DenseMatrix m_dN(2,4);
+   DenseMatrix m_dN2(3,4);
+   BasisEvalDerivs(xi, m_N, m_dN, m_dN2);
+
+   int dim = 3;
+   int num_dofs1 = dim;
+   int num_dofs2 = 4*dim;
+   int num_dofs = num_dofs1 + num_dofs2;
+   dg2dx.SetSize(num_dofs,num_dofs); dg2dx = 0.0;
+
+   Vector x_c(3);
+   m_coords.MultTranspose(m_N,x_c);
+
+   Vector gap_v(3); gap_v = 0.0;
+   gap_v  = x_s;
+   gap_v -= x_c;
+
+   DenseMatrix m_dx(2,3);
+   Mult(m_dN, m_coords, m_dx);
+
+   DenseMatrix m_dx2(3,3); m_dx2 = 0.0;
+   Mult(m_dN2,m_coords, m_dx2);
+   double nnorm = 0.0;
+   Vector normal(3); normal = 0.0;
+   ComputeNormal(m_dN, m_coords, normal, nnorm);
+
+   double gap = gap_v * normal; // gap function value, dot product between vectors
+
+   DenseMatrix M(2,2); M = 0.0;
+   MultABt(m_dx, m_dx, M);
+
+   DenseMatrix f(2, num_dofs2); f = 0.0;
+
+   for (int d=0; d<3; d++)
+   {
+      DenseMatrix Mtemp(2,2); Mtemp = 0.0;
+      Mtemp(0,0) = m_dx2(0,d); Mtemp(0,1) = m_dx2(1,d);
+      Mtemp(1,0) = m_dx2(1,d); Mtemp(1,1) = m_dx2(2,d);
+
+      M.Add(-gap_v[d], Mtemp);
+
+      Vector m_dxcol(2); m_dx.GetColumn(d, m_dxcol);
+      DenseMatrix ftmp(2,4);
+      outer(m_dxcol, m_N, ftmp);
+      ftmp *= -1;
+      ftmp.Add( gap_v[d], m_dN);   // 2*4
+
+      for (int j=0; j<4; j++)
+      {
+         assert(d+3*j<num_dofs2);
+         f(0,d+j*3) = ftmp(0,j);
+         f(1,d+j*3) = ftmp(1,j);
+      }
+   }
+   //fprintf('hess dxidxm\n');
+   DenseMatrixInverse Minv(M);
+   DenseMatrix dxidxm(2,num_dofs2); dxidxm = 0.0;
+   Minv.Mult(f, dxidxm);
+   //LinearSolve??
+   //dxidxm = M\f;
+
+   DenseMatrix nde2(2,2); nde2 = 0.0;
+   DenseMatrix Nndx2(2,num_dofs2); Nndx2 = 0.0;
+
+   for (int d=0; d<3; d++)
+   {
+      DenseMatrix ndetmp(2,2); ndetmp = 0.0;
+      ndetmp(0,0) = normal(d)*m_dx2(0,d); ndetmp(0,1) = normal(d)*m_dx2(1,d);
+      ndetmp(1,0) = normal(d)*m_dx2(1,d); ndetmp(1,1) = normal(d)*m_dx2(2,d);
+
+      nde2 += ndetmp;
+
+      for (int j=0; j<4; j++)
+      {
+         assert(d+3*j<num_dofs2);
+         Nndx2(0,d+j*3) = normal[d]*m_dN(0,j);
+         Nndx2(1,d+j*3) = normal[d]*m_dN(1,j);
+      }
+
+   }
+
+   DenseMatrix Ndn(2,num_dofs2); Ndn = 0.0;
+   Ndn += Nndx2;
+   AddMult(nde2, dxidxm, Ndn);
+
+
+   DenseMatrix M2(2,2); M2 = 0.0;
+   MultABt(m_dx, m_dx, M2);
+   DenseMatrixInverse M2inv(M2);
+   DenseMatrix diag2(2,2); diag2(0,0) = 1.0; diag2(1,1) = 1.0;
+   DenseMatrix m_con(2,2); m_con = 0.0;
+
+   M2inv.Mult(diag2, m_con);
+
+   DenseMatrix dg2dxm(num_dofs2, num_dofs2); dg2dxm = 0.0;
+
+   DenseMatrix dg2dxm_tmp(num_dofs2,2); dg2dxm_tmp = 0.0;
+   MultAtB(Ndn, m_con, dg2dxm_tmp);
+   Mult(dg2dxm_tmp, Ndn, dg2dxm);
+   dg2dxm *= gap;
+
+   DenseMatrix dg2dxm_tmp2(num_dofs2,num_dofs2); dg2dxm_tmp2 = 0.0;
+   MultAtB(Nndx2, dxidxm, dg2dxm_tmp2);
+   dg2dxm.Add(-1.0, dg2dxm_tmp2);
+
+   dg2dxm_tmp = 0.0;
+   MultAtB(dxidxm, nde2, dg2dxm_tmp);
+
+   AddMult_a(-1.0, dg2dxm_tmp, dxidxm, dg2dxm);
+
+   dg2dxm_tmp2 = 0.0;
+   MultAtB(dxidxm, Nndx2, dg2dxm_tmp2);
+   dg2dxm.Add(-1.0, dg2dxm_tmp2);
+
+   Vector v_dxidx2(4);
+   m_coords.Mult(gap_v, v_dxidx2); // m_coords * gap_v;  // 4*3 * 3 = 4
+
+   DenseMatrix K_dxidx2(2,2); K_dxidx2 = 0.0;
+
+   Vector m_dN2row1(4); m_dN2.GetRow(0, m_dN2row1);
+   Vector m_dN2row2(4); m_dN2.GetRow(1, m_dN2row2);
+   Vector m_dN2row3(4); m_dN2.GetRow(2, m_dN2row3);
+   K_dxidx2(0,0) = m_dN2row1 * v_dxidx2; // how would 4*1 * 1*4 be computed?
+   K_dxidx2(0,1) = m_dN2row2 * v_dxidx2;
+   K_dxidx2(1,0) = m_dN2row2 * v_dxidx2;
+   K_dxidx2(1,1) = m_dN2row3 * v_dxidx2;
+
+   DenseMatrix K_dxidx(2,2);
+   K_dxidx -= M2;
+   K_dxidx += K_dxidx2;
+
+   Vector drdxs_r(6);
+   drdxs_r[0] = m_dx(0,0); drdxs_r[1] = m_dx(0,1); drdxs_r[2] = m_dx(0,2);
+   drdxs_r[3] = m_dx(1,0); drdxs_r[4] = m_dx(1,1); drdxs_r[5] = m_dx(1,2);
+
+   DenseMatrix drdxs_K(6,6); drdxs_K = 0.;
+   for (int i=0; i<3; i++)
+   {
+      drdxs_K(i,i) = K_dxidx(0,0);
+      drdxs_K(i,3+i) = K_dxidx(0,1);
+      drdxs_K(i+3,i) = K_dxidx(1,0);
+      drdxs_K(i+3,i+3) = K_dxidx(1,1);
+   }
+   Vector dxidxs(6);
+
+   DenseMatrixInverse drdxsK_inv(drdxs_K);
+   drdxsK_inv.Mult(drdxs_r,dxidxs);
+   dxidxs *= -1.0;
+   //dxidxs = -drdxs_K\drdxs_r;
+
+   DenseMatrix dxidxs_m(2,3); dxidxs_m = 0.0;
+   dxidxs_m(0,0) = dxidxs[0]; dxidxs_m(0,1) = dxidxs[1]; dxidxs_m(0,2) = dxidxs[2];
+   dxidxs_m(1,0) = dxidxs[3]; dxidxs_m(1,1) = dxidxs[4]; dxidxs_m(1,2) = dxidxs[5];
+
+   DenseMatrix dtao1dxs(3,3); dtao1dxs = 0.0;
+   DenseMatrix dtao2dxs(3,3); dtao2dxs = 0.0;
+
+   Vector dxidxs_row1(3); dxidxs_row1 = 0.0; Vector dxidxs_row2(3);
+   dxidxs_row2 = 0.0;
+   Vector mdx2_row1(3); mdx2_row1 = 0.0; Vector mdx2_row2(3); mdx2_row2 = 0.0;
+   Vector mdx2_row3(3); mdx2_row3 = 0.0;
+   dxidxs_m.GetRow(0,dxidxs_row1);
+   dxidxs_m.GetRow(1,dxidxs_row2);
+   m_dx2.GetRow(0,mdx2_row1);
+   m_dx2.GetRow(1,mdx2_row2);
+   m_dx2.GetRow(2,mdx2_row3);
+
+   DenseMatrix dtaotmp(3,3); dtaotmp = 0.0;
+   outer(mdx2_row1, dxidxs_row1,dtaotmp);
+   dtao1dxs += dtaotmp; dtaotmp = 0.0;
+   outer(mdx2_row2, dxidxs_row1,dtaotmp);
+   dtao1dxs += dtaotmp; dtaotmp = 0.0;
+
+   outer(mdx2_row2, dxidxs_row2, dtaotmp);
+   dtao2dxs += dtaotmp; dtaotmp = 0.0;
+   outer(mdx2_row3, dxidxs_row2, dtaotmp);
+   dtao2dxs += dtaotmp; dtaotmp = 0.0;
+
+   DenseMatrix dtaodxs(3,3); dtaodxs = 0.0; //tao = tao1 cross tao2
+
+   for (int d=0; d<3; d++)
+   {
+      Vector dtao1dxs_tmp(3); dtao1dxs_tmp = 0.0;
+      dtao1dxs.GetColumn(d,dtao1dxs_tmp);
+      Vector m_dxrow(3); m_dx.GetRow(1, m_dxrow);
+
+      Vector dtaodxs_tmp(3); dtaodxs_tmp = 0.0;
+      cross(dtao1dxs_tmp, m_dxrow, dtaodxs_tmp);
+
+      Vector dtaodxs_tmp2(3); dtaodxs_tmp2 = 0.0;
+      m_dx.GetRow(0, m_dxrow);
+      dtao1dxs_tmp = 0.0;  // reuse the same vector for dtao2
+      dtao2dxs.GetColumn(d,dtao1dxs_tmp);
+      cross(m_dxrow, dtao1dxs_tmp, dtaodxs_tmp2);
+
+      dtaodxs_tmp2 += dtaodxs_tmp;
+      dtaodxs.SetCol(d, dtaodxs_tmp2);
+   }
+
+   DenseMatrix dndxs(3,3); dndxs = 0.0; dndxs += dtaodxs; dndxs *= 1.0/nnorm;
+   DenseMatrix dndxs_tmp(3,3); dndxs_tmp = 0.0;
+   outer(normal, normal, dndxs_tmp);
+   AddMult_a(-1/nnorm, dndxs_tmp, dtaodxs, dndxs);
+
+   DenseMatrix dgvdxs(3,3); dgvdxs = 0.0;
+   MultAtB(m_dx, dxidxs_m, dgvdxs);
+   dgvdxs *= -1;
+   for (int d=0; d<3; d++)
+   {
+      dgvdxs(d,d) += 1.0;
+   }
+   //dxidxs: 2*3
+
+   DenseMatrix dg2dxs(3,3); dg2dxs = 0.0;
+   DenseMatrix dg2dxs_tmp(3,2); dg2dxs_tmp = 0.0;
+   MultAtB(dxidxs_m, nde2, dg2dxs_tmp);
+   AddMult_a(-1.0, dg2dxs_tmp, dxidxs_m, dg2dxs);
+   DenseMatrix dg2dxs_tmp2(3,3); dg2dxs_tmp2 = 0.0;
+   MultAtB(dgvdxs, dndxs, dg2dxs_tmp2);
+   dg2dxs += dg2dxs_tmp2;
+   dg2dxs_tmp2 = 0.0;
+   MultAtB(dndxs, dndxs_tmp, dg2dxs_tmp2);
+   AddMult(dg2dxs_tmp2, dgvdxs, dg2dxs);
+
+   DenseMatrix Ne(3,12), Be(6,12), dBe(12,12);
+   BasisVectorDerivs(xi, Ne, Be, dBe);
+
+   DenseMatrix dtao1dxm(3,12); dtao1dxm.CopyRows(Be, 0, 2);
+   DenseMatrix dtao2dxm(3,12); dtao2dxm.CopyRows(Be, 3, 5);
+
+   Vector m_coords_v(12);
+   for (int i=0; i<4; i++)
+   {
+      for (int j=0; j<3; j++)
+      {
+         m_coords_v[i*3+j] = m_coords(i,j);
+      }
+   }
+
+   for (int i=0; i<2; i++)
+   {
+      Vector dxidxm_tmp(num_dofs2); dxidxm_tmp = 0.0;
+      dxidxm.GetRow(i,dxidxm_tmp);
+
+      DenseMatrix dBe_tmp(3,12);
+      dBe_tmp.CopyRows(dBe,i*3,(i+1)*3-1);
+
+      DenseMatrix dtaodxm_tmp(12,12); dtaodxm_tmp = 0.0;
+      outer(m_coords_v, dxidxm_tmp, dtaodxm_tmp);
+      AddMult(dBe_tmp, dtaodxm_tmp, dtao1dxm);
+
+      //dtao1dxm += dBe(:,:,i)*reshape(m_coords(1:4,:)',12,1)*reshape(dxidxm(i,:),1,12); % 3*12
+      dBe_tmp = 0.0;
+      dBe_tmp.CopyRows(dBe,(i+2)*3,(i+3)*3-1);
+      AddMult(dBe_tmp, dtaodxm_tmp, dtao2dxm);
+
+   }
+
+   DenseMatrix dtaodxm(3,12); dtaodxm = 0.0;//tao = tao1 cross tao2
+
+   for (int d=0; d<12; d++)
+   {
+      Vector dtaodxm_tmp(3); dtaodxm_tmp = 0.0;
+      Vector dtaodxm_tmp2(3); dtaodxm_tmp2 = 0.0;
+      Vector tmp1(3); tmp1 = 0.0;  dtao1dxm.GetColumn(d,tmp1);
+      Vector m_dxrow2(3); m_dx.GetRow(1, m_dxrow2);
+      Vector m_dxrow1(3); m_dx.GetRow(0, m_dxrow1);
+      Vector tmp2(3); tmp2 = 0.0;  dtao2dxm.GetColumn(d,tmp2);
+
+      cross(tmp1, m_dxrow2, dtaodxm_tmp);
+      cross(m_dxrow1,tmp2, dtaodxm_tmp2);
+      dtaodxm_tmp += dtaodxm_tmp2;
+
+      dtaodxm.SetCol(d, dtaodxm_tmp);
+   }
+
+   DenseMatrix dndxm(3,12); dndxm = 0.0;
+   dndxm += dtaodxm;
+   dndxm *= 1.0/nnorm;
+   AddMult_a(-1/nnorm, dndxs_tmp, dtaodxm, dndxm); //dndxs_tmp = normal'*normal
+
+   DenseMatrix dgvdxm(3,12); dgvdxm = 0.0;
+   dgvdxm -= Ne;
+
+   for (int i=0; i<2; i++)
+   {
+      Vector dxidxm_tmp(num_dofs2); dxidxm_tmp = 0.0;
+      dxidxm.GetRow(i,dxidxm_tmp);
+
+      DenseMatrix Be_tmp(3,12);
+      Be_tmp.CopyRows(Be,i*3,(i+1)*3-1);
+
+      DenseMatrix dgvdxm_tmp(12,12); dgvdxm_tmp = 0.0;
+      outer(m_coords_v, dxidxm_tmp, dgvdxm_tmp);
+      AddMult_a(-1.0, Be_tmp, dgvdxm_tmp, dgvdxm);
+
+   }
+
+   DenseMatrix dg2dxsxm(3,12); dg2dxsxm = 0.0;
+   DenseMatrix dg2dxsxm_tmp(3,3); dg2dxsxm_tmp = 0.0;
+   MultAtB(dgvdxs, dndxm, dg2dxsxm);
+
+   MultAtB(dndxs, dndxs_tmp, dg2dxsxm_tmp);
+   AddMult(dg2dxsxm_tmp, dgvdxm, dg2dxsxm); // += dndxs'*normal'*normal*dgvdxm;
+
+   DenseMatrix dgvdxsxmn(3,12); dgvdxsxmn = 0.0;
+   DenseMatrix dgvdxsxmn_tmp(3,2); dgvdxsxmn_tmp = 0.0;
+   MultAtB(dxidxs_m, nde2, dgvdxsxmn_tmp);   //dxidxs_m: 2*3
+
+   AddMult_a(-1.0, dgvdxsxmn_tmp, dxidxm, dgvdxsxmn);
+
+
+   for (int i =0; i<2; i++)
+   {
+      DenseMatrix Be_tmp(3,12);
+      Be_tmp.CopyRows(Be,i*3,(i+1)*3-1);
+
+      Vector dxidxs_row(3); dxidxs_row = 0.0; dxidxs_m.GetRow(i,dxidxs_row);
+      DenseMatrix dgvdxsxmn_tmp2(3,3); dgvdxsxmn_tmp2 = 0.0;
+      outer(dxidxs_row, normal, dgvdxsxmn_tmp2);
+      AddMult_a(-1.0, dgvdxsxmn_tmp2, Be_tmp, dgvdxsxmn);
+   }
+
+   dg2dxsxm += dgvdxsxmn;
+
+   DenseMatrix dg2dxmxs(12,3); dg2dxmxs = 0.0;
+   DenseMatrix dg2dxmxs_tmp(12,3); dg2dxmxs_tmp = 0.0;
+   MultAtB(dgvdxm, dndxs, dg2dxmxs);
+   MultAtB(dndxm, dndxs_tmp, dg2dxmxs_tmp);
+   AddMult(dg2dxmxs_tmp, dgvdxs, dg2dxmxs);
+
+   DenseMatrix dgvdxmxsn(12,3); dgvdxmxsn = 0.0;
+   DenseMatrix dgvdxmxsn_tmp(12,2); dgvdxmxsn_tmp = 0.0;
+
+   MultAtB(dxidxm, nde2, dgvdxmxsn_tmp);
+   dgvdxmxsn_tmp *= -1.0;
+   AddMult(dgvdxmxsn_tmp, dxidxs_m, dgvdxmxsn);
+
+   for (int i =0; i<2; i++)
+   {
+      DenseMatrix Be_tmp(3,12);
+      Be_tmp.CopyRows(Be,i*3,(i+1)*3-1);
+      Be_tmp.Transpose(); // Be is now 12*3
+
+      Vector dxidxs_row(3); dxidxs_row = 0.0; dxidxs_m.GetRow(i,dxidxs_row);
+      DenseMatrix dgvdxmxsn_tmp2(3,3); dgvdxmxsn_tmp2 = 0.0;
+      outer(normal, dxidxs_row, dgvdxmxsn_tmp2);
+      AddMult_a(-1.0, Be_tmp, dgvdxmxsn_tmp2, dgvdxmxsn);
+
+   }
+
+   dg2dxmxs += dgvdxmxsn;
+
+   dg2dx.CopyMN(dg2dxs, 0, 0);
+   dg2dx.CopyMN(dg2dxm, 3, 3);
+   dg2dx.CopyMN(dg2dxsxm, 0, 3);
+   dg2dx.CopyMN(dg2dxmxs, 3, 0);
+
+};
+
+
+
+void NodeSegConPairs(const Vector x1, const Vector xi2,
+                     const DenseMatrix coords2,
+                     double& node_g, Vector& node_dg, DenseMatrix& node_dg2)
+{
+   double gap = 0.0;
+   Vector normal(3); normal = 0.0;
+   Vector dgdxm(12); dgdxm = 0.0;
+   Vector dgdxs(3); dgdxs = 0.0;
+
+   ComputeGapJacobian(x1, xi2, coords2, gap, normal, dgdxm, dgdxs);
+   node_g = gap;
+
+   node_dg.SetSize(12+3);
+   for (int i=0; i<3; i++) { node_dg[i] = dgdxs[i]; }
+   for (int i=0; i<12; i++) { node_dg[i+3] = dgdxm[i]; }
+
+   DenseMatrix dg2dx(15,15); dg2dx = 0.0;
+   DenseMatrix dgvdxmxsn(12,3); dgvdxmxsn = 0.0;
+   ComputeGapHessian(x1, xi2, coords2, dg2dx);
+
+   node_dg2.SetSize(15,15);
+   node_dg2 = dg2dx;
+
+   /*
+   if(obj.space1.conns{e1}(i)==150)  % for debugging purpose
+
+   v1 = 1:3;
+   v2 = 1:12;
+   %v1 = ones(1,3)
+   %v2 = ones(1,12)
+   v2 = reshape(v2,4,3);
+   x1n1 = x1 + 0.01*v1;
+   coords2n1 = coords2 + 0.001*v2;
+      [xi2n1, gapv1, ~, ~] = SlaveToMaster(obj, coords2n1, x1n1);
+      [gapn1, n1,dgdxmn1, dgdxsn1] = ComputeGapJacobian(obj, x1n1, xi2n1, coords2n1);
+   x1n2 = x1 - 0.01*v1;
+   coords2n2 = coords2 - 0.001*v2;
+      [xi2n2, gapv2, ~, ~] = SlaveToMaster(obj, coords2n2, x1n2);
+      [gapn2, n2,dgdxmn2, dgdxsn2] = ComputeGapJacobian(obj, x1n2, xi2n2, coords2n2);
+   fprintf('fd\n');
+      %gapv1-gapv2
+   [dgdxsn1(:)',dgdxmn1(:)'] - [dgdxsn2(:)',dgdxmn2(:)']
+
+   %dgdxsn1-dgdxsn2
+   fprintf('code\n');
+   v2n = v2';
+   %dg2dx(1:3,1:3)*0.04*ones(3,1)
+   temp = zeros(12,3);
+   for i = 1:4
+   temp1 = dg2dx(3+(i-1)*3+1:3+i*3,1:3);
+   temp((i-1)*3+1:i*3,:) = temp1';
+   end
+   temp2 = zeros(3,12);
+   for i = 1:4
+   temp3 = dg2dx(1:3,3+(i-1)*3+1:3+i*3);
+   temp2(:,(i-1)*3+1:i*3) = temp3';
+   end
+   %dg2dx
+   %dg2dx(4:end,1:3) = temp;
+   %dg2dx(1:3,4:end) = temp2;
+   %dgvdxm * 0.002*v2n(:)
+   (dg2dx*[0.02*v1(:)',0.002*v2n(:)']')'
+   %dg2dx(4:end,1:3)
+   end*/
+
+};
+
+
+// coordsm : (npoints*4, 3) use what class?
+// m_conn: (npoints*4)
+void Assemble_Contact(const int m, const int npoints, const int ndofs,
+                      const Vector x_s,
+                      const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector& g, SparseMatrix& M,
+                      std::vector<SparseMatrix>& dM)
+{
+   int ndim = 3;
+
+   g.SetSize(m);
+   g = 0.0;
+
+   //SparseMatrix M(m, n);  // M needs to be the correct size
+
+   //dM.resize(m); // needs to clear?
+
+   double g_tmp = 0.;
+   Vector dg(4*ndim+ndim);
+   dg = 0.;
+   DenseMatrix dg2(4*ndim+ndim,4*ndim+ndim);
+   dg2 = 0.;
+
+   for (int i=0; i<npoints; i++)
+   {
+      Vector x1(ndim);
+      x1[0] = x_s[i*ndim];
+      x1[1] = x_s[i*ndim+1];
+      x1[2] = x_s[i*ndim+2];
+
+      Vector xi2(ndim-1);
+      xi2[0] = xi[i*(ndim-1)];
+      xi2[1] = xi[i*(ndim-1)+1];
+
+      DenseMatrix coords2(4,3);
+      coords2.CopyRows(coordsm, i*4,(i+1)*4-1);
+
+      //how to get coords2?
+      dg = 0.0;
+      dg2 = 0.;
+      NodeSegConPairs(x1, xi2, coords2, g_tmp, dg, dg2);
+      //x1.Print();
+      //xi2.Print();
+      //coords2.Print();
+      g[s_conn[i]] = g_tmp; // should be unique
+      Array<int> m_conn_i(4);
+      m_conn.GetSubArray(4*i, 4, m_conn_i);
+
+      Array<int> node_conn(5);
+      node_conn[0] = s_conn[i];
+      for (int j=0; j<4; j++)
+      {
+         node_conn[j+1] = m_conn_i[j];
+      }
+
+      Array<int>  M_i_tmp(1);
+      M_i_tmp[0] = s_conn[i];
+
+      //j_idx = (node_conn-1)*obj.disp_field.num_components +repmat((1:obj.disp_field.num_components)', 1, length(node_conn{i}));
+      Array<int> j_idx(5*ndim); j_idx = 0;
+      for (int j=0; j< 5; j++)
+      {
+         for (int k=0; k<ndim; k++)
+         {
+            j_idx[j*ndim+k] = node_conn[j]*ndim+k;
+         }
+      }
+      DenseMatrix M_v_tmp(1, ndim*(4+1));   // SetData now?
+      M_v_tmp.SetRow(0, dg);
+
+      M.AddSubMatrix(M_i_tmp, j_idx, M_v_tmp);
+
+      Array<int> dM_i(ndim*(4+1));
+      Array<int> dM_j(ndim*(4+1));
+
+      for (int j=0; j< ndim*(4+1); j++)
+      {
+         dM_i[j] = j_idx[j];
+         dM_j[j] = j_idx[j];
+      }
+      //dg2.Print();
+      //dM[s_conn[i]].Print();
+      dM[s_conn[i]].AddSubMatrix(dM_i,dM_j, dg2);
+   }
+};
+

examples/contact/obstacleProblem.cpp

@@ -0,0 +1,119 @@
+#include "mfem.hpp"
+#include "Problems.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+
+double dmanufacturedFun(const Vector &);
+double fRhs(const Vector &);
+double obstacle(const Vector &);
+
+int main(int argc, char *argv[])
+{
+  int FEorder = 1; // order of the finite elements
+  int linSolver = 0;
+  int maxIPMiters = 30;
+  bool iAmRoot = true;
+  
+  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.Parse();
+  if(!args.Good())
+  {
+    args.PrintUsage(cout);
+    return 1;
+  }
+  else
+  {
+    if( iAmRoot )
+    {
+      args.PrintOptions(cout);
+    }
+  }
+
+  const char *meshFile = "../../data/inline-quad.mesh";
+  Mesh *mesh = new Mesh(meshFile, 1, 1);
+  int dim = mesh->Dimension(); // geometric dimension of the domain
+  {
+     int ref_levels = 3;
+     for (int l = 0; l < ref_levels; l++)
+     {
+        mesh->UniformRefinement();
+     }
+  }
+
+  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
+  ObstacleProblem problem(Vh, &fRhs, &obstacle);
+  
+  int dimD = problem.GetDimD();
+  Vector x0(dimD); x0 = 0.0;
+  Vector xf(dimD); xf = 0.0;
+
+  InteriorPointSolver optimizer(&problem); 
+  optimizer.SetTol(1.e-7);
+  optimizer.SetLinearSolver(linSolver);
+  optimizer.SetMaxIter(maxIPMiters);
+  optimizer.Mult(x0, xf);
+
+  GridFunction d_gf(Vh);
+
+  d_gf = xf;
+
+  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
+  GridFunction dm_gf(Vh);
+  dm_gf.ProjectCoefficient(dm_fc);
+  
+  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
+  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("d(x) (numerical)", &d_gf);
+  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
+  paraview_dc.Save();
+  
+  
+  delete Vh;
+  delete fec;
+  delete mesh;
+  return 0;
+}
+
+
+double dmanufacturedFun(const Vector &x)
+{
+  return cos(2*M_PI*x(0)) + 0.2 - 2.0*(pow(x(0),3) - 1.5*pow(x(0),2));
+}
+
+
+// f(x) forcing term... which enters the objective energy functional
+// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
+// of f(x). f(x) is such that in the absence of bound-constraints then
+// the solution of the optimization problem satisfies the PDE
+// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
+// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
+double fRhs(const Vector &x)
+{
+  double fx = 0.;
+  fx = 0.2 - 2.0 * (pow(x(0),3)- 1.5*pow(x(0),2.) - 6 * x(0) + 3.) + (1. + pow(2.*M_PI,2))*cos(2.*M_PI*x(0));
+  return fx;
+}
+
+double obstacle(const Vector &x)
+{
+  return 0.0;
+}

examples/contact/rotatedblock2.mesh

@@ -0,0 +1,70 @@
+
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+3
+
+# 1 nothing
+elements
+4
+1   5     0     1     3     2     6     7    9     8
+1   5     2     3     5     4     8     9    11    10
+1   5     6     7     9     8    12    13    15    14
+1   5     8     9    11    10    14    15    17    16
+
+# 0 nothing
+# 1 dirichlet bc
+# 2 contact
+boundary
+16
+1   3     1     0     2     3
+1   3     3     2     4     5
+1   3    12    13    15    14
+1   3    14    15    17    16
+3   3     2     0     6     8
+3   3     4     2     8    10
+3   3     8     6    12    14
+3   3    10     8    14    16
+2   3     1     3     9     7
+2   3     3     5    11    9
+2   3     7    9    15    13
+2   3     9    11    17    15
+1   3     6     0     1     7
+1   3    12     6     7    13
+1   3     4    10    11     5
+1   3    10    16    17    11
+
+vertices
+18
+3
+
+0.000000000000 0.145770950245 0.443895630208
+0.507100000000 0.145770950245 0.443895630208
+0.000000000000 0.350937660019 0.294833290227
+0.507100000000 0.350937660019 0.294833290227
+0.000000000000 0.556104369792 0.145770950245
+0.507100000000 0.556104369792 0.145770950245
+0.000000000000 0.294833290227 0.649062339981
+0.507100000000 0.294833290227 0.649062339981
+0.000000000000 0.500000000000 0.500000000000
+0.507100000000 0.500000000000 0.500000000000
+0.000000000000 0.705166709773 0.350937660019
+0.507100000000 0.705166709773 0.350937660019
+0.000000000000 0.443895630208 0.854229049755
+0.507100000000 0.443895630208 0.854229049755
+0.000000000000 0.649062339981 0.705166709773
+0.507100000000 0.649062339981 0.705166709773
+0.000000000000 0.854229049755 0.556104369792
+0.507100000000 0.854229049755 0.556104369792

examples/makefile

@@ -36,7 +36,7 @@ ifeq ($(MFEM_USE_MPI),NO)
 else
    EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
 endif
-SUBDIRS =
+SUBDIRS = contact 
 ifeq ($(MFEM_USE_AMGX),YES)
    SUBDIRS += amgx
 endif

makefile

@@ -119,7 +119,7 @@ $(if $(word 2,$(SRC)),$(error Spaces in SRC = "$(SRC)" are not supported))
 MFEM_GIT_STRING = $(shell [ -d $(MFEM_DIR)/.git ] && git -C $(MFEM_DIR) \
    describe --all --long --abbrev=40 --dirty --always 2> /dev/null)
 
-EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
+EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith contact
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples