serial/parallel builds can run serial example codes, serial build can run serial examples

miniapps/contact/ObstacleProblem.cpp

@@ -0,0 +1,119 @@
+//                         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 "ipsolver/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 = 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
+   {
+     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();
+      }
+   }
+
+   FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+   FiniteElementSpace   *Vh  = new FiniteElementSpace(&mesh, fec);
+
+   ObstacleProblem problem(Vh,&fRhs, &obstacle);
+  
+   int dimD = problem.GetDimU();
+   Vector x0(dimD); x0 = 100.0;
+   Vector xf(dimD); xf = 0.0;
+
+   InteriorPointSolver optimizer(&problem); 
+   optimizer.SetTol(1.e-8);
+   optimizer.SetLinearSolveTol(1.e-10);
+   optimizer.SetLinearSolver(linSolver);
+   optimizer.SetMaxIter(maxIPMiters);
+   optimizer.Mult(x0, xf);
+
+   GridFunction d_gf(Vh);
+   d_gf.Set(1.0, xf);
+
+
+
+   FunctionCoefficient dm_fc(dmanufacturedFun); // 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;
+   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;
+}
+
+double obstacle(const Vector &x)
+{
+  return 0.0;
+}

miniapps/contact/contact_driver.cpp

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

miniapps/contact/ipsolver/IPsolver.hpp

@@ -0,0 +1,88 @@
+#include "mfem.hpp"
+#include "../problems/problems.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+#ifndef IP_SOLVER 
+#define IP_SOLVER
+
+class InteriorPointSolver
+{
+protected:
+   GeneralOptProblem* 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;
+
+   SparseMatrix * Huu, * Hum, * Hmu, * Hmm, * Wmm, * D, * Ju, * Jm, * JuT, * JmT;
+#ifdef MFEM_USE_MPI
+   HypreParMatrix * Huuh, * Humh, * Hmuh, * Hmmh, * Wmmh, * Dh, * Juh, * Jmh, * JuTh, * JmTh;
+#endif
+   
+   int jOpt;
+   bool converged;
+   
+   int MyRank;
+   bool iAmRoot;
+   bool parallel;
+
+   bool saveLogBarrierIterates;
+
+   int linSolver;
+   double linSolveTol;
+public:
+   InteriorPointSolver(GeneralOptProblem*);
+   double MaxStepSize(Vector& , Vector& , Vector& , double);
+   double MaxStepSize(Vector& , Vector& , double);
+   void Mult(const BlockVector& , BlockVector&);
+   void Mult(const Vector&, Vector &); 
+   void GetLagrangeMultiplier(Vector &);
+   void FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
+   void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , 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 ~InteriorPointSolver();
+};
+
+#endif

miniapps/contact/makefile

@@ -0,0 +1,126 @@
+# 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 = ./
+
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/contact/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+CONTACT_SEQ_SRC = problems/problems.cpp problems/problems_util.cpp util/util.cpp util/mpicomm.cpp ipsolver/IPsolver.cpp
+CONTACT_SEQ_OBJ = $(CONTACT_PAR_SRC:.cpp=.o)
+CONTACT_PAR_SRC = $(CONTACT_SEQ_SRC)
+CONTACT_PAR_OBJ = $(CONTACT_PAR_SRC:.cpp=.o)
+
+OBSTACLE_SEQ_SRC = $(CONTACT_SEQ_SRC)
+OBSTACLE_SEQ_OBJ = $(OBSTACLE_SEQ_SRC:.cpp=.o)
+
+OBSTACLE_PAR_SRC = $(OBSTACLE_SEQ_SRC)
+OBSTACLE_PAR_OBJ = $(OBSTACLE_PAR_SRC:.cpp=.o)
+
+OBSTACLE_SRC = ObstacleProblem.cpp $(OBSTACLE_SEQ_SRC)
+OBSTACLE_OBJ = $(OBSTACLE_SRC:.cpp=.o)
+
+POBSTACLE_SRC = pObstacleProblem.cpp $(OBSTACLE_PAR_SRC)
+POBSTACLE_OBJ = $(POBSTACLE_SRC:.cpp=.o)
+
+CONTACT_SRC = contact_driver.cpp $(CONTACT_SEQ_SRC)
+CONTACT_OBJ = $(CONTACT_SRC:.cpp=.o)
+
+PCONTACT_SRC = pcontact_driver.cpp $(CONTACT_PAR_SRC)
+PCONTACT_OBJ = $(PCONTACT_SRC:.cpp=.o)
+
+SEQ_MINIAPPS = contact_driver obstacle_problem
+PAR_MINIAPPS = pcontact_driver pobstacle_problem
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
+
+# If MFEM_SHARED is set, add the ../common rpath
+COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all lib-common clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+%.o: %.cpp
+
+%.o: $(SRC)%.cpp $(wildcard $(SRC)%.hpp) $(MFEM_LIB_FILE)\
+ $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+problems/%.o: $(SRC)problems/%.cpp $(wildcard $(SRC)problems/%.hpp) $(MFEM_LIB_FILE)\
+ $(CONFIG_MK) | lib-common
+	mkdir -p $(@D)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+
+all: $(MINIAPPS)
+
+obstacle_problem: $(OBSTACLE_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(OBSTACLE_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+pobstacle_problem: $(POBSTACLE_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(POBSTACLE_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+contact_driver: $(CONTACT_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(CONTACT_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+pcontact_driver: $(PCONTACT_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(PCONTACT_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+
+
+# Rule for building lib-common
+lib-common:
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Specific execution options
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+contact-test-seq: diffusion
+	@$(call mfem-test,$<,, contact miniapp,)
+pcontact-test-par: pcontact
+	@$(call mfem-test,$<, $(RUN_MPI), pcontact miniapp,)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+	rm -f $(CONTACT_OBJ) $(PCONTACT_OBJ)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf ParaView 

miniapps/contact/meshes/block1.mesh

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

miniapps/contact/meshes/block2.mesh

@@ -0,0 +1,68 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+3
+
+# 1 nothing
+elements
+4
+1   5     0     1     3     2     6     7    9     8
+1   5     2     3     5     4     8     9    11    10
+1   5     6     7     9     8    12    13    15    14
+1   5     8     9    11    10    14    15    17    16
+
+# 0 nothing
+# 1 dirichlet bc
+# 2 contact
+boundary
+16
+1   3     1     0     2     3
+1   3     3     2     4     5
+1   3    12    13    15    14
+1   3    14    15    17    16
+3   3     2     0     6     8
+3   3     4     2     8    10
+3   3     8     6    12    14
+3   3    10     8    14    16
+2   3     1     3     9     7
+2   3     3     5    11    9
+2   3     7    9    15    13
+2   3     9    11    17    15
+1   3     6     0     1     7
+1   3    12     6     7    13
+1   3     4    10    11     5
+1   3    10    16    17    11
+
+vertices
+18
+3
+         0    0.2464    0.2464
+    0.5071    0.2464    0.2464
+         0    0.5000    0.2464
+    0.5071    0.5000    0.2464
+         0    0.7536    0.2464
+    0.5071    0.7536    0.2464
+         0    0.2464    0.5000
+    0.5071    0.2464    0.5000
+         0    0.5000    0.5000
+    0.5071    0.5000    0.5000
+         0    0.7536    0.5000
+    0.5071    0.7536    0.5000
+         0    0.2464    0.7536
+    0.5071    0.2464    0.7536
+         0    0.5000    0.7536
+    0.5071    0.5000    0.7536
+         0    0.7536    0.7536
+    0.5071    0.7536    0.7536

miniapps/contact/meshes/rotatedblock2.mesh

@@ -0,0 +1,70 @@
+
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+3
+
+# 1 nothing
+elements
+4
+1   5     0     1     3     2     6     7    9     8
+1   5     2     3     5     4     8     9    11    10
+1   5     6     7     9     8    12    13    15    14
+1   5     8     9    11    10    14    15    17    16
+
+# 0 nothing
+# 1 dirichlet bc
+# 2 contact
+boundary
+16
+1   3     1     0     2     3
+1   3     3     2     4     5
+1   3    12    13    15    14
+1   3    14    15    17    16
+3   3     2     0     6     8
+3   3     4     2     8    10
+3   3     8     6    12    14
+3   3    10     8    14    16
+2   3     1     3     9     7
+2   3     3     5    11    9
+2   3     7    9    15    13
+2   3     9    11    17    15
+1   3     6     0     1     7
+1   3    12     6     7    13
+1   3     4    10    11     5
+1   3    10    16    17    11
+
+vertices
+18
+3
+
+0.000 0.362753 0.168656
+0.507100 0.362753 0.168656
+0.000 0.597049 0.265704
+0.507100 0.597049 0.265704
+0.0000 0.831344 0.362753
+0.507100 0.831344 0.362753
+0.0000 0.265704 0.402951
+0.507100 0.265704 0.402951
+0.0000 0.500000 0.500000
+0.507100 0.500000 0.500000
+0.0000 0.734296 0.597049
+0.507100 0.734296 0.597049
+0.0000 0.168656 0.637247
+0.507100 0.168656 0.637247
+0.0000 0.402951 0.734296
+0.507100 0.402951 0.734296
+0.0000 0.637247 0.831344
+0.507100 0.637247 0.831344

miniapps/contact/pObstacleProblem.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 "ipsolver/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[])
+{
+  // 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);
+
+   ObstacleProblem problem(Vh,&fRhs, &obstacle);
+  
+   int dimD = problem.GetDimU();
+   Vector x0(dimD); x0 = 100.0;
+   Vector xf(dimD); xf = 0.0;
+
+   InteriorPointSolver 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);
+   ParaViewDataCollection paraview_dc("BarrierProblemSolution", &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("d(x) (numerical)", &d_gf);
+   paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
+   paraview_dc.Save();
+
+   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;
+}
+
+double obstacle(const Vector &x)
+{
+  return 0.1;
+}

miniapps/contact/pcontact_driver.cpp

@@ -0,0 +1,344 @@
+//                               Parallel contact example
+//
+// Compile with: make pcontact_driver
+// sample run
+// mpirun -np 6 ./pcontact_driver -sr 2 -pr 2
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "ipsolver/IPsolver.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
+   Hypre::Init();
+   // 1. Parse command-line options.
+   const char *mesh_file1 = "meshes/block1.mesh";
+   const char *mesh_file2 = "meshes/rotatedblock2.mesh";
+   int order = 1;
+   int sref = 0;
+   int pref = 0;
+   Array<int> attr;
+   Array<int> m_attr;
+   bool visualization = true;
+   bool paraview = false;
+   double linsolvertol = 1e-12;
+   int relax_type = 8;
+   double optimizer_tol = 1e-10;
+   bool elasticity_options = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file1, "-m1", "--mesh1",
+                  "First mesh file to use.");
+   args.AddOption(&mesh_file2, "-m2", "--mesh2",
+                  "Second mesh file to use.");
+   args.AddOption(&attr, "-at", "--attributes-surf",
+                  "Attributes of boundary faces on contact surface for mesh 2.");
+   args.AddOption(&sref, "-sr", "--serial-refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&pref, "-pr", "--parallel-refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&linsolvertol, "-stol", "--solver-tol",
+                  "Linear Solver Tolerance.");
+   args.AddOption(&optimizer_tol, "-otol", "--optimizer-tol",
+                  "Interior Point Solver Tolerance.");
+   args.AddOption(&relax_type, "-rt", "--relax-type",
+                  "Selection of Smoother for AMG");
+   args.AddOption(&elasticity_options, "-elast", "--elasticity-options",
+                  "-no-elast",
+                  "--no-elasticity-options",
+                  "Enable or disable Elasticity options for the AMG preconditioner.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
+                  "--no-paraview",
+                  "Enable or disable ParaView visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   ElasticityProblem * prob1 = nullptr;
+   ElasticityProblem * prob2 = nullptr;
+
+   ParMesh * pmesh1 = nullptr;
+   ParMesh * pmesh2 = nullptr;
+   ParMesh * pmesh  = nullptr;
+
+   bool own_mesh = true;
+   if (elasticity_options) { own_mesh = true; }
+
+   if (own_mesh)
+   {
+      Mesh * mesh1 = new Mesh(mesh_file1,1);
+      Mesh * mesh2 = new Mesh(mesh_file2,1);
+
+      for (int i = 0; i<sref; i++)
+      {
+         mesh1->UniformRefinement();
+         mesh2->UniformRefinement();
+      }
+
+      int * part1 = mesh1->GeneratePartitioning(num_procs);
+      int * part2 = mesh2->GeneratePartitioning(num_procs);
+
+      Array<int> part_array1(part1, mesh1->GetNE());
+      Array<int> part_array2(part2, mesh2->GetNE());
+
+      for (int i = 0; i<mesh1->GetNE(); i++)
+      {
+         mesh1->SetAttribute(i,1);
+      }
+      mesh1->SetAttributes();
+      for (int i = 0; i<mesh2->GetNE(); i++)
+      {
+         mesh2->SetAttribute(i,2);
+      }
+      mesh2->SetAttributes();
+
+      Mesh * mesh_array[2];
+      mesh_array[0] = mesh1;
+      mesh_array[1] = mesh2;
+      Mesh * mesh = new Mesh(mesh_array, 2);
+
+
+
+      Array<int> part_array;
+      part_array.Append(part_array1);
+      part_array.Append(part_array2);
+
+      pmesh = new ParMesh(MPI_COMM_WORLD, *mesh, part_array.GetData());
+      MFEM_VERIFY(pmesh->GetNE(), "Empty partition");
+      for (int i = 0; i<pref; i++)
+      {
+         pmesh->UniformRefinement();
+      }
+
+      Array<int> attr1(1); attr1 = 1;
+      Array<int> attr2(1); attr2 = 2;
+
+
+      pmesh1 = new ParMesh(MPI_COMM_WORLD,*mesh1, part1);
+      pmesh2 = new ParMesh(MPI_COMM_WORLD,*mesh2, part2);
+
+      for (int i = 0; i<pref; i++)
+      {
+         pmesh1->UniformRefinement();
+         pmesh2->UniformRefinement();
+      }
+
+      MFEM_VERIFY(pmesh1->GetNE(), "Empty partition mesh1");
+      MFEM_VERIFY(pmesh2->GetNE(), "Empty partition mesh2");
+
+      prob1 = new ElasticityProblem(pmesh1,order);
+      prob2 = new ElasticityProblem(pmesh2,order);
+   }
+   else
+   {
+      prob1 = new ElasticityProblem(MPI_COMM_WORLD, mesh_file1,sref,pref,order);
+      prob2 = new ElasticityProblem(MPI_COMM_WORLD, mesh_file2,sref,pref,order);
+   }
+   cout << "constructed ElasticityProblem\n";
+
+   Vector lambda1(prob1->GetMesh()->attributes.Max()); lambda1 = 57.6923076923;
+   Vector mu1(prob1->GetMesh()->attributes.Max()); mu1 = 38.4615384615;
+   Vector lambda2(prob2->GetMesh()->attributes.Max()); lambda2 = 57.6923076923;
+   Vector mu2(prob2->GetMesh()->attributes.Max()); mu2 = 38.4615384615;
+
+   prob1->SetLambda(lambda1); prob1->SetMu(mu1);
+   prob2->SetLambda(lambda2); prob2->SetMu(mu2);
+
+   ContactProblem contact(prob1,prob2);
+   InteriorPointSolver optimizer(&contact);
+   ParFiniteElementSpace *pfes = nullptr;
+   if (pmesh)
+   {
+      pfes = new ParFiniteElementSpace(pmesh,prob1->GetFECol(),3,Ordering::byVDIM);
+      pmesh->SetNodalFESpace(pfes);
+      //if (elasticity_options)
+      //{
+      //   optimizer.SetFiniteElementSpace(pfes);
+      //}
+   }
+
+   optimizer.SetTol(optimizer_tol);
+   optimizer.SetMaxIter(30);
+
+   //int linsolver = 2;
+   //optimizer.SetLinearSolver(linsolver);
+   //optimizer.SetLinearSolveTol(linsolvertol);
+   //optimizer.SetLinearSolveRelaxType(relax_type);
+
+   ParGridFunction x1 = prob1->GetDisplacementParGridFunction();
+   ParGridFunction x2 = prob2->GetDisplacementParGridFunction();
+
+   int ndofs1 = prob1->GetNumTDofs();
+   int ndofs2 = prob2->GetNumTDofs();
+   int gndofs1 = prob1->GetGlobalNumDofs();
+   int gndofs2 = prob2->GetGlobalNumDofs();
+   int ndofs = ndofs1 + ndofs2;
+
+   Vector X1 = x1.GetTrueVector();
+   Vector X2 = x2.GetTrueVector();
+
+   Vector x0(ndofs); x0 = 0.0;
+   x0.SetVector(X1,0);
+   x0.SetVector(X2,X1.Size());
+
+   Vector xf(ndofs); xf = 0.0;
+   optimizer.SetLinearSolveTol(1.e-14);
+   optimizer.Mult(x0, xf);
+
+   double Einitial = contact.E(x0);
+   double Efinal = contact.E(xf);
+   //Array<int> & CGiterations = optimizer.GetCGIterNumbers();
+   //if (Mpi::Root())
+   //{
+   //   mfem::out << endl;
+   //   mfem::out << " Initial Energy objective     = " << Einitial << endl;
+   //   mfem::out << " Final Energy objective       = " << Efinal << endl;
+   //   mfem::out << " Global number of dofs        = " << gndofs1 + gndofs2 << endl;
+   //   mfem::out << " Global number of constraints = " << numconstr << endl;
+   //   mfem::out << " CG iteration numbers         = " ;
+   //   CGiterations.Print(mfem::out, CGiterations.Size());
+   //}
+
+   MFEM_VERIFY(optimizer.GetConverged(),
+               "Interior point solver did not converge.");
+
+   if (visualization || paraview)
+   {
+      ParFiniteElementSpace * fes1 = dynamic_cast<ParFiniteElementSpace *>(prob1->GetFESpace());
+      ParFiniteElementSpace * fes2 = dynamic_cast<ParFiniteElementSpace *>(prob2->GetFESpace());
+
+      ParMesh * mesh1 = fes1->GetParMesh();
+      ParMesh * mesh2 = fes2->GetParMesh();
+
+      Vector X1_new(xf.GetData(),fes1->GetTrueVSize());
+      Vector X2_new(&xf.GetData()[fes1->GetTrueVSize()],fes2->GetTrueVSize());
+
+      //contact.ComputeGapFunctionAndDerivatives(X1_new, X2_new);
+      //Vector gxf;
+      //contact.g(xf, gxf);
+      //for(int i = 0; i < gxf.Size(); i++)
+      //{
+      //  mfem::out << "g(xf)_i = " << gxf(i) << endl;
+      //}
+
+      ParGridFunction x1_gf(fes1);
+      ParGridFunction x2_gf(fes2);
+
+      x1_gf.SetFromTrueDofs(X1_new);
+      x2_gf.SetFromTrueDofs(X2_new);
+
+      mesh1->MoveNodes(x1_gf);
+      mesh2->MoveNodes(x2_gf);
+
+
+      ParGridFunction * xgf = nullptr;
+      if (pmesh)
+      {
+         xgf = new ParGridFunction(pfes);
+         *xgf = 0.0;
+         // auto map1 = ParSubMesh::CreateTransferMap(x1_gf, *xgf);
+         // auto map2 = ParSubMesh::CreateTransferMap(x2_gf, *xgf);
+         // map1.Transfer(x1_gf,*xgf);
+         // map2.Transfer(x2_gf,*xgf);
+         xgf->SetVector(x1_gf,0);
+         xgf->SetVector(x2_gf,x1_gf.Size());
+         ParFiniteElementSpace * pfes1 = x1_gf.ParFESpace();
+         ParFiniteElementSpace * pfes2 = x2_gf.ParFESpace();
+         pmesh->MoveNodes(*xgf);
+         
+      }
+      if (paraview)
+      {
+         ParaViewDataCollection paraview_dc1("QPContactBody1", mesh1);
+         paraview_dc1.SetPrefixPath("ParaView");
+         paraview_dc1.SetLevelsOfDetail(1);
+         paraview_dc1.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc1.SetHighOrderOutput(true);
+         paraview_dc1.SetCycle(0);
+         paraview_dc1.SetTime(0.0);
+         paraview_dc1.RegisterField("Body1", &x1_gf);
+         paraview_dc1.Save();
+
+         ParaViewDataCollection paraview_dc2("QPContactBody2", mesh2);
+         paraview_dc2.SetPrefixPath("ParaView");
+         paraview_dc2.SetLevelsOfDetail(1);
+         paraview_dc2.SetDataFormat(VTKFormat::BINARY);
+         paraview_dc2.SetHighOrderOutput(true);
+         paraview_dc2.SetCycle(0);
+         paraview_dc2.SetTime(0.0);
+         paraview_dc2.RegisterField("Body2", &x2_gf);
+         paraview_dc2.Save();
+
+         if (pmesh)
+         {
+            ParaViewDataCollection paraview_dc("QPContact2Bodies", pmesh);
+            paraview_dc.SetPrefixPath("ParaView");
+            paraview_dc.SetLevelsOfDetail(1);
+            paraview_dc.SetDataFormat(VTKFormat::BINARY);
+            paraview_dc.SetHighOrderOutput(true);
+            paraview_dc.SetCycle(0);
+            paraview_dc.SetTime(0.0);
+            paraview_dc.RegisterField("BothBodies", xgf);
+            paraview_dc.Save();
+         }
+      }
+
+
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int visport = 19916;
+
+         {
+            socketstream sol_sock(vishost, visport);
+            sol_sock.precision(8);
+            sol_sock << "parallel " << 2*num_procs << " " << myid << "\n"
+                     << "solution\n" << *mesh1 << x1_gf << flush;
+         }
+         {
+            socketstream sol_sock(vishost, visport);
+            sol_sock.precision(8);
+            sol_sock << "parallel " << 2*num_procs << " " << myid+num_procs << "\n"
+                     << "solution\n" << *mesh2 << x2_gf << flush;
+         }
+         if (pmesh)
+         {
+            socketstream sol_sock1(vishost, visport);
+            sol_sock1.precision(8);
+            sol_sock1 << "parallel " << num_procs << " " << myid << "\n"
+                      << "solution\n" << *pmesh << *xgf << flush;
+         }
+      }
+      delete xgf;
+   }
+
+   delete pfes;
+   delete prob2;
+   delete prob1;
+   delete pmesh2;
+   delete pmesh1;
+
+   return 0;
+}

miniapps/contact/problems/problems.hpp

@@ -0,0 +1,296 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "problems_util.hpp"
+
+using namespace std;
+using namespace mfem;
+
+#ifndef PROBLEM_DEFS
+#define PROBLEM_DEFS
+
+// abstract GeneralOptProblem class
+// of the form
+// min_(u,m) f(u,m) s.t. c(u,m)=0 and m>=ml
+// the primal variable (u, m) is represented as a BlockVector
+class GeneralOptProblem
+{
+protected:
+    int dimU, dimM, dimC;
+    int dimUGlb, dimMGlb, dimCGlb;
+    #ifdef MFEM_USE_MPI
+      HYPRE_BigInt * dofOffsetsU;
+      HYPRE_BigInt * dofOffsetsM;
+    #endif
+    Array<int> block_offsetsx;
+    Vector ml;
+    bool parallel;
+public:
+    GeneralOptProblem();
+    virtual double CalcObjective(const BlockVector &) = 0;
+    virtual void Duf(const BlockVector &, Vector &) = 0;
+    virtual void Dmf(const BlockVector &, Vector &) = 0;
+    void CalcObjectiveGrad(const BlockVector &, BlockVector &);
+    #ifdef MFEM_USE_MPI
+       void InitGeneral(HYPRE_BigInt * dofOffsetsU_, HYPRE_BigInt * dofOffsetsM_);
+       HYPRE_BigInt * GetDofOffsetsU() const { return dofOffsetsU; };
+       HYPRE_BigInt * GetDofOffsetsM() const { return dofOffsetsM; };
+    #endif
+    virtual void InitGeneral(int dimU, int dimM);
+    virtual Operator * Duuf(const BlockVector &) = 0;
+    virtual Operator * Dumf(const BlockVector &) = 0;
+    virtual Operator * Dmuf(const BlockVector &) = 0;
+    virtual Operator * Dmmf(const BlockVector &) = 0;
+    virtual Operator * Duc(const BlockVector &) = 0;
+    virtual Operator * Dmc(const BlockVector &) = 0;
+    virtual void c(const BlockVector &, Vector &) = 0;
+    int GetDimU() const { return dimU; };
+    int GetDimM() const { return dimM; }; 
+    int GetDimC() const { return dimC; };
+    int GetDimUGlb() const { return dimUGlb; };
+    int GetDimMGlb() const { return dimMGlb; };
+    int GetDimCGlb() const { return dimCGlb; };
+    bool IsParallel() const { return parallel; };
+    Vector Getml() const { return ml; };
+    ~GeneralOptProblem();
+};
+
+
+// Specialized optimization problem
+// of the form
+// min_d e(d) s.t. g(d) >= 0
+// suited for contact mechanics problems that can be formualted
+// as an optimization problem
+class OptProblem : public GeneralOptProblem
+{
+protected:
+#ifdef MFEM_USE_MPI
+   HypreParMatrix * Ih;
+#endif
+   SparseMatrix * Isparse;
+public:
+   OptProblem();
+    
+   // GeneralOptProblem methods are defined in terms of
+   // OptProblem specific methods: E, DdE, DddE, g, Ddg
+   double CalcObjective(const BlockVector &); 
+   void Duf(const BlockVector &, Vector &);
+   void Dmf(const BlockVector &, Vector &);
+#ifdef MFEM_USE_MPI
+   void Init(HYPRE_BigInt *, HYPRE_BigInt *);
+#endif
+    void Init(int, int);
+    Operator * Duuf(const BlockVector &);
+    Operator * Dumf(const BlockVector &);
+    Operator * Dmuf(const BlockVector &);
+    Operator * Dmmf(const BlockVector &);
+    Operator * Duc(const BlockVector &);
+    Operator * Dmc(const BlockVector &);
+
+    void c(const BlockVector &, Vector &);
+    
+    // ParOptProblem specific methods:
+    
+    // energy objective function e(d)
+    // input: d an mfem::Vector
+    // output: e(d) a double
+    virtual double E(const Vector &d) = 0;
+
+    // gradient of energy objective De / Dd
+    // input: d an mfem::Vector,
+    //        gradE an mfem::Vector, which will be the gradient of E at d
+    // output: none    
+    virtual void DdE(const Vector &d, Vector &gradE) = 0;
+
+    // Hessian of energy objective D^2 e / Dd^2
+    // input:  d, an mfem::Vector
+    // output: The Hessian of the energy objective at d, a pointer to a (HyprePar or Sparse) Matrix
+    virtual Operator * DddE(const Vector &d) = 0;
+
+    // Constraint function g(d) >= 0, e.g., gap function
+    // input: d, an mfem::Vector,
+    //       gd, an mfem::Vector, which upon successfully calling the g method will be
+    //                            the evaluation of the function g at d
+    // output: none
+    virtual void g(const Vector &d, Vector &gd) = 0;
+
+    // Jacobian of constraint function Dg / Dd, e.g., gap function Jacobian
+    // input:  d, an mfem::Vector,
+    // output: The Jacobain of the constraint function g at d, a pointer to a (HyprePar or Sparse) Matrix
+    virtual Operator * Ddg(const Vector &) = 0;
+    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
+   // stiffness matrix used to define objective
+#ifdef MFEM_USE_MPI
+   ParFiniteElementSpace * Vhp;
+   ParBilinearForm * Kformp;
+   ParLinearForm   * fformp;
+   HypreParMatrix    Kh;
+   HypreParMatrix  * Jh;
+#endif
+   FiniteElementSpace * Vh;
+   BilinearForm * Kform;
+   LinearForm   * fform;
+   SparseMatrix K;
+   SparseMatrix * J;
+   Array<int> ess_tdof_list; 
+   Vector f;
+   Vector psi;
+public :
+   double E(const Vector & d);
+   void DdE(const Vector & d, Vector & dE);
+   void   g(const Vector & d, Vector & gd);
+   Operator * DddE(const Vector & d);
+   Operator * Ddg (const Vector & d);
+   ObstacleProblem(FiniteElementSpace*, double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &));
+   virtual ~ObstacleProblem();
+};
+
+
+class ElasticityProblem
+{
+   private:
+      bool formsystem = false;
+      bool own_mesh;
+      bool parallel;
+      Mesh * mesh = nullptr;
+      int order;
+      int ndofs;
+      int ntdofs;
+      int gndofs;
+      FiniteElementCollection * fec = nullptr;
+      FiniteElementSpace      * fes = nullptr;
+      Vector lambda, mu;
+      PWConstCoefficient lambda_cf, mu_cf;
+      Array<int> ess_bdr, ess_tdof_list;
+      BilinearForm *a = nullptr;
+      LinearForm    b;
+      GridFunction  x;
+      SparseMatrix A;
+      Vector X, B;
+      
+#ifdef MFEM_USE_MPI 
+      MPI_Comm comm;
+      ParMesh * pmesh = nullptr;
+      ParFiniteElementSpace * fesp = nullptr;
+      ParBilinearForm * ap = nullptr;
+      ParLinearForm bp;
+      ParGridFunction xp;
+      HypreParMatrix Ap;
+#endif
+      void Init();
+      
+   public:
+      ElasticityProblem(const char * mesh_file, int ref, int order_);
+#ifdef MFEM_USE_MPI
+      ElasticityProblem(MPI_Comm comm_, const char *mesh_file , int sref, int pref, int order_);
+      ElasticityProblem(ParMesh * pmesh_, int order_);
+#endif
+      Mesh * GetMesh();
+      FiniteElementCollection * GetFECol() { return fec; };
+      FiniteElementSpace * GetFESpace();
+      int GetNumDofs() { return ndofs; };
+      int GetNumTDofs() { return ntdofs; };
+      int GetGlobalNumDofs() { return gndofs; };
+      Operator & GetOperator();
+      Vector & GetRHS();
+      void SetLambda(const Vector & lambda_);
+      void SetMu(const Vector & mu_);
+      void FormLinearSystem();
+      void UpdateLinearSystem();
+      void SetDisplacementDirichletData(const Vector & delta);
+      GridFunction & GetDisplacementGridFunction();
+#ifdef MFEM_USE_MPI
+      ParGridFunction & GetDisplacementParGridFunction();
+#endif
+      Array<int> & GetEssentialDofs() { return ess_tdof_list; };
+      bool IsParallel() const { return parallel; };
+      ~ElasticityProblem();
+};
+
+
+class ContactProblem : public OptProblem
+{
+private:
+   int numprocs;
+   int myid;
+   ElasticityProblem * prob1 = nullptr;
+   ElasticityProblem * prob2 = nullptr;
+   FiniteElementSpace * vfes1 = nullptr;
+   FiniteElementSpace * vfes2 = nullptr;
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+   ParFiniteElementSpace * vfes1p = nullptr;
+   ParFiniteElementSpace * vfes2p = nullptr;
+#endif
+   int dim;
+   GridFunction nodes0;
+   GridFunction *nodes1 = nullptr;
+   std::set<int> contact_vertices;
+   bool recompute = true;
+   bool compute_hessians = true;
+   std::vector<int> dof_offsets;
+   std::vector<int> vertex_offsets; 
+   std::vector<int> constraints_offsets; 
+   Array<int> tdof_offsets;
+   Array<int> constraints_starts;
+   Array<int> globalvertices1;
+   Array<int> globalvertices2;
+   Array<int> vertices2;
+   Array<int> vertices1;
+   
+   protected:
+      bool parallel;
+      int npoints=0;
+      int gnpoints=0;
+      int nv, gnv;
+      SparseMatrix * K = nullptr;
+      BlockVector *B = nullptr;
+      Vector gapv;
+      SparseMatrix * M = nullptr;
+      Array<SparseMatrix *> dM;
+#ifdef MFEM_USE_MPI
+      HypreParMatrix * Kp = nullptr;
+      HypreParMatrix * Mp = nullptr;
+      Array<HypreParMatrix*> dMp;
+      void ParComputeContactVertices();
+#endif
+      void ComputeContactVertices();
+   
+   public:
+      ContactProblem(ElasticityProblem * prob1_, ElasticityProblem * prob2_);
+#ifdef MFEM_USE_MPI
+      MPI_Comm GetComm() { return comm; };
+      void ParComputeGapFunctionAndDerivatives(const Vector & displ1, const Vector & displ2);
+#endif
+      void ComputeGapFunctionAndDerivatives(const Vector & displ1, const Vector & displ2);
+      int GetNumDofs();
+      int GetGlobalNumDofs();
+      int GetNumConstraints() { return npoints; };
+      int GetGlobalNumConstraints() { return gnpoints; };
+      std::vector<int> & GetDofOffets() { return dof_offsets; }
+      std::vector<int> & GetVertexOffsets() { return vertex_offsets; }
+      std::vector<int> & GetConstraintsOffsets() { return constraints_offsets; }
+      Array<int> & GetConstraintsStarts() { return constraints_starts; }
+      Vector & GetGapFunction() { return gapv; }
+      Operator * GetJacobian();
+      Array<Operator *> GetHessian();
+      double E(const Vector & d);
+      void DdE(const Vector &d, Vector &gradE);
+      Operator * DddE(const Vector &d);
+      void g(const Vector &d, Vector &gd); // todo : add in cpp file
+      Operator * Ddg(const Vector &d);
+};
+
+
+#endif

miniapps/contact/problems/problems_util.hpp

@@ -0,0 +1,88 @@
+#include "mfem.hpp"
+#include "../util/mpicomm.hpp"
+
+using namespace std;
+using namespace mfem;
+
+void BasisEval(const Vector xi, Vector &N, DenseMatrix &dNdxi); // dNdxi is 2*4
+void BasisEvalDerivs(const Vector xi, Vector& N, DenseMatrix& dNdxi,
+                     DenseMatrix& dN2dxi);
+// returns the vector and matrix form of the shape functions and its derivative
+void BasisVectorDerivs(const Vector xi, DenseMatrix& N, DenseMatrix& dNdxi,
+                       DenseMatrix& ddNdxi);
+void cross(const Vector a, const Vector b, Vector& c);
+// a outer b
+void outer(const Vector a, const Vector b, DenseMatrix& c);
+// dphidxi 2*4
+// coords 4*3
+void ComputeNormal(const DenseMatrix& dphidxi, const DenseMatrix& coords,
+                   Vector& normal, double& nnorm);
+void SlaveToMaster(const DenseMatrix& m_coords, const Vector& s_x, Vector& xi);
+
+// m_coords is expected to be 4 * 3
+void  ComputeGapJacobian(const Vector x_s, const Vector xi,
+                         const DenseMatrix m_coords,
+                         double& gap, Vector& normal, Vector& dgdxm, Vector& dgdxs);
+
+void ComputeGapHessian(const Vector x_s, const Vector xi,
+                       const DenseMatrix m_coords,
+                       DenseMatrix& dg2dx);
+void NodeSegConPairs(const Vector x1, const Vector xi2,
+                     const DenseMatrix coords2,
+                     double& node_g, Vector& node_dg, DenseMatrix& node_dg2);
+// coordsm : (npoints*4, 3) use what class?
+// m_conn: (npoints*4)
+void Assemble_Contact(const Vector x_s,
+                      const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector& g, SparseMatrix& M,
+                      Array<SparseMatrix *> & dM);
+
+void Assemble_Contact(const Vector x_s,
+                      const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector & g, SparseMatrix & M1, SparseMatrix & M2, 
+                      Array<SparseMatrix *> & dM11,
+                      Array<SparseMatrix *> & dM12,
+                      Array<SparseMatrix *> & dM21,
+                      Array<SparseMatrix *> & dM22);
+void Assemble_Contact(const Vector x_s,
+                      const Vector xi, const DenseMatrix coordsm, const Array<int> s_conn,
+                      const Array<int> m_conn, Vector & g, SparseMatrix & M1, SparseMatrix & M2,const Array<int> & points_map); 
+
+void FindSurfaceToProject(Mesh& mesh, const int elem, int& cbdrface);
+
+Vector GetNormalVector(Mesh & mesh, const int elem, const double *ref,
+                       int & refFace, int & refNormal, bool & interior);
+int GetHexVertex(int cdim, int c, int fa, int fb, Vector & refCrd);
+
+// Coordinates in xyz are assumed to be ordered as [X, Y, Z]
+// where X is the list of x-coordinates for all points and so on.
+// conn: connectivity of the target surface elements
+// xi: surface reference cooridnates for the cloest point, involves a linear transformation from [0,1] to [-1,1]
+void FindPointsInMesh(Mesh & mesh, Vector const& xyz, Array<int>& conn, Vector& xi);
+
+#ifdef MFEM_USE_MPI
+// Coordinates in xyz are assumed to be ordered as [X, Y, Z]
+// where X is the list of x-coordinates for all points and so on.
+// conn: connectivity of the target surface elements
+// xi: surface reference cooridnates for the cloest point, involves a linear transformation from [0,1] to [-1,1]
+void FindPointsInMesh(Mesh & mesh, const Array<int> & gvert, const Vector & xyz, const Array<int> & s_conn, Array<int>& conn,
+                      Vector & xyz2, Array<int> & s_conn2, Vector& xi, DenseMatrix & coords);
+
+// somewhat simplified version of the above
+void FindPointsInMesh(Mesh & mesh, const Array<int> & gvert, Array<int> & s_conn, const Vector &x1, Vector & xyz, Array<int>& conn,
+                      Vector& xi, DenseMatrix & coords);                   
+
+int get_rank(int tdof, std::vector<int> & tdof_offsets);
+void ComputeTdofOffsets(const ParFiniteElementSpace * pfes,
+                        std::vector<int> & tdof_offsets);
+void ComputeTdofOffsets(MPI_Comm comm, int mytoffset, std::vector<int> & tdof_offsets);
+void ComputeTdofs(MPI_Comm comm, int mytoffs, std::vector<int> & tdofs);
+
+
+// Performs Pᵀ * A * P for BlockOperator  P (with blocks as HypreParMatrices)
+// and A a HypreParMatrix, i.e., this handles the special case 
+// where P = [P₁ P₂ ⋅⋅⋅ Pₙ] 
+void RAP(const HypreParMatrix & A, const BlockOperator & P, BlockOperator & C);
+void ParAdd(const BlockOperator & A, const BlockOperator & B, BlockOperator & C);
+
+#endif

miniapps/contact/util/mpicomm.cpp

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

miniapps/contact/util/mpicomm.hpp

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

miniapps/contact/util/util.cpp

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

miniapps/contact/util/util.hpp

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