minor bug

.gitignore

@@ -131,6 +131,12 @@ examples/hiop/ex9-mesh.*
 examples/hiop/ex9-init.*
 examples/hiop/ex9-final.*
 
+examples/ipopt/exContactBlockTL
+examples/ipopt/exContactBlockTL.mesh
+examples/ipopt/exContactBlockTL-mesh.*
+examples/ipopt/exContactBlockTL-init.*
+examples/ipopt/exContactBlockTL-final.*
+
 examples/petsc/ex[1-69]p
 examples/petsc/ex1[0-1]p
 examples/petsc/mesh.*

CMakeLists.txt

@@ -404,6 +404,15 @@ if (MFEM_USE_HIOP)
   # find_package updates HIOP_FOUND, HIOP_INCLUDE_DIRS, HIOP_LIBRARIES
 endif()
 
+# IpOpt optimizer
+if (MFEM_USE_IPOPT)
+  find_package(IPOPT REQUIRED)
+  message(
+      STATUS
+      "IPOPT_INCLUDE_DIRS=${IPOPT_INCLUDE_DIRS}, IPOPT_LIBRARIES=${IPOPT_LIBRARIES}, IPOPT_DIR=${IPOPT_DIR}")
+  # find_package updates IPOPT_FOUND, IPOPT_INCLUDE_DIRS, IPOPT_LIBRARIES
+endif()
+
 # CoDiPack package
 if (MFEM_USE_CODIPACK)
   find_package(CODIPACK REQUIRED)
@@ -499,7 +508,7 @@ find_package(Threads REQUIRED)
 #    be before SuiteSparse.
 set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
     SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
-    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
+    NETCDF MPFR PUMI HIOP IPOPT POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
     ADIOS2 CUBLAS CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
     MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)
 

CONTRIBUTING.md

@@ -112,6 +112,7 @@ The MFEM source code has the following structure:
   │   ├── caliper
   │   ├── ginkgo
   │   ├── hiop
+  │   ├── ipopt
   │   ├── jupyter
   │   ├── moonolith
   │   ├── petsc

INSTALL

@@ -471,6 +471,9 @@ MFEM_USE_HIOP = YES/NO
    Enable the usage of HiOp (https://github.com/LLNL/hiop) in MFEM. HiOp is an
    HPC solver for nonlinear optimization problems.
 
+MFEM_USE_IPOPT = YES/NO
+   Enable the usage of Ipopt in MFEM.
+
 MFEM_USE_CODIPACK = YES/NO
   Enable automatic differentiation using the CoDiPack library.
   www.scicomp.uni-kl.de/codi/
@@ -738,6 +741,11 @@ The specific libraries and their options are:
   Options: HIOP_OPT, HIOP_LIB.
   Versions: HIOP >= 0.4.6.
 
+- Ipopt (optional), used when MFEM_USE_IPOPT = YES.
+  URL: https://github.com/coin-or/Ipopt
+  Options: IPOPT_OPT, IPOPT_LIB.
+  Versions: IPOPT >= 3.14
+
 - CoDiPack (optional), used with MFEM_USE_CODIPACK = YES
   URL: https://www.scicomp.uni-kl.de/codi/
   Options: CODIPACK_OPT
@@ -972,6 +980,7 @@ MFEM_USE_MPFR
 MFEM_USE_ZLIB
 MFEM_USE_PUMI
 MFEM_USE_HIOP
+MFEM_USE_IPOPT
 MFEM_USE_CODIPACK
 MFEM_USE_ADFORWARD
 MFEM_USE_CUDA
@@ -1035,6 +1044,7 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - POSIXCLOCKS
  - PUMI
  - HIOP
+ - IPOPT
  - CoDiPack
  - OCCA
  - RAJA

config/cmake/MFEMConfig.cmake.in

@@ -36,6 +36,7 @@ set(MFEM_USE_STRUMPACK @MFEM_USE_STRUMPACK@)
 set(MFEM_USE_GINKGO @MFEM_USE_GINKGO@)
 set(MFEM_USE_AMGX @MFEM_USE_AMGX@)
 set(MFEM_USE_HIOP @MFEM_USE_HIOP@)
+set(MFEM_USE_IPOPT @MFEM_USE_IPOPT@)
 set(MFEM_USE_GNUTLS @MFEM_USE_GNUTLS@)
 set(MFEM_USE_GSLIB @MFEM_USE_GSLIB@)
 set(MFEM_USE_NETCDF @MFEM_USE_NETCDF@)

config/cmake/config.hpp.in

@@ -131,6 +131,9 @@
 // Enable MFEM functionality based on the HiOp library
 #cmakedefine MFEM_USE_HIOP
 
+// Enable MFEM functionality based on the Ipopt library
+#cmakedefine MFEM_USE_IPOPT
+
 // Build the GPU/CUDA-enabled version of the MFEM library.
 // Requires a CUDA compiler (nvcc).
 #cmakedefine MFEM_USE_CUDA

config/cmake/modules/FindIPOPT.cmake

@@ -0,0 +1,23 @@
+# Copyright (c) 2010-2022, 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.
+
+# Sets the following variables:
+#   - IPOPT_FOUND
+#   - IPOPT_INCLUDE_DIRS
+#   - IPOPT_LIBRARIES
+
+include(MfemCmakeUtilities)
+mfem_find_package(IPOPT IPOPT IPOPT_DIR
+  "include" "IpTNLP.hpp"
+  "lib" "ipopt"
+  "Paths to headers required by IPOPT."
+  "Libraries required by IPOPT.")
+

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -869,7 +869,7 @@ function(mfem_export_mk_files)
       MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS MFEM_USE_STRUMPACK
       MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_GNUTLS MFEM_USE_NETCDF
       MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS
-      MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA
+      MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_IPOPT MFEM_USE_GSLIB MFEM_USE_CUDA
       MFEM_USE_HIP MFEM_USE_RAJA MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER
       MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO
       MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG

config/config.hpp.in

@@ -141,6 +141,9 @@
 // Enable MFEM functionality based on the HIOP library.
 // #define MFEM_USE_HIOP
 
+// Enable MFEM functionality based on the IPOPT library.
+// #define MFEM_USE_IPOPT
+
 // Enable MFEM functionality based on the GSLIB library
 // #define MFEM_USE_GSLIB
 

config/config.mk.in

@@ -46,6 +46,7 @@ MFEM_USE_FMS           = @MFEM_USE_FMS@
 MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@
 MFEM_USE_HIOP          = @MFEM_USE_HIOP@
+MFEM_USE_IPOPT          = @MFEM_USE_IPOPT@
 MFEM_USE_GSLIB         = @MFEM_USE_GSLIB@
 MFEM_USE_CUDA          = @MFEM_USE_CUDA@
 MFEM_USE_HIP           = @MFEM_USE_HIP@

config/defaults.cmake

@@ -48,6 +48,7 @@ option(MFEM_USE_FMS "Enable FMS usage" OFF)
 option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
 option(MFEM_USE_PUMI "Enable PUMI" OFF)
 option(MFEM_USE_HIOP "Enable HiOp" OFF)
+option(MFEM_USE_IPOPT "Enable Ipopt" OFF)
 option(MFEM_USE_CUDA "Enable CUDA" OFF)
 option(MFEM_USE_HIP "Enable HIP" OFF)
 option(MFEM_USE_OCCA "Enable OCCA" OFF)
@@ -220,6 +221,10 @@ set(HIOP_DIR "${MFEM_DIR}/../hiop/install" CACHE STRING
     "Directory where HiOp is installed")
 set(HIOP_REQUIRED_PACKAGES "BLAS" "LAPACK" CACHE STRING
     "Packages that HiOp depends on.")
+set(IPOPT_DIR "${MFEM_DIR}/../ipopt/install" CACHE STRING
+    "Directory where IpOpt is installed")
+set(IPOPT_REQUIRED_PACKAGES "BLAS" "LAPACK" CACHE STRING
+    "Packages that IpOpt depends on.")
 
 set(MKL_CPARDISO_DIR "" CACHE STRING "MKL installation path.")
 set(MKL_MPI_WRAPPER_LIB "mkl_blacs_mpich_lp64" CACHE STRING "MKL MPI wrapper library")

config/defaults.mk

@@ -148,6 +148,7 @@ MFEM_USE_FMS           = NO
 MFEM_USE_CONDUIT       = NO
 MFEM_USE_PUMI          = NO
 MFEM_USE_HIOP          = NO
+MFEM_USE_IPOPT         = NO
 MFEM_USE_GSLIB         = NO
 MFEM_USE_CUDA          = NO
 MFEM_USE_HIP           = NO
@@ -447,6 +448,11 @@ HIOP_DIR = @MFEM_DIR@/../hiop/install
 HIOP_OPT = -I$(HIOP_DIR)/include
 HIOP_LIB = -L$(HIOP_DIR)/lib -lhiop $(LAPACK_LIB)
 
+# IPOPT
+IPOPT_DIR = @MFEM_DIR@/../ipopt/install
+IPOPT_OPT = -I$(IPOPT_DIR)/include
+IPOPT_LIB = -L$(IPOPT_DIR)/lib -lipopt $(LAPACK_LIB)
+
 # CoDiPack
 CODIPACK_DIR = @MFEM_DIR@/../CoDiPack
 CODIPACK_OPT = -I$(CODIPACK_DIR)

config/sample-runs.sh

@@ -58,6 +58,10 @@ groups_serial=(
    "HiOp examples:"
    "examples/hiop"
    "ex9.cpp"'
+'"ipopt"
+   "IpOpt examples:"
+   "examples/ipopt"
+   "ex10.cpp"'
 '"pumi"
    "PUMI examples:"
    "examples/pumi"
@@ -215,6 +219,10 @@ groups_all=(
    "HiOp examples:"
    "examples/hiop"
    "ex9.cpp ex9p.cpp"'
+'"ipopt"
+   "IpOpt examples:"
+   "examples/ipopt"
+   "ex10.cpp"'
 '"pumi"
    "PUMI examples:"
    "examples/pumi"

doc/CodeDocumentation.conf.in

@@ -785,6 +785,7 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                          @MFEM_SOURCE_DIR@/examples/caliper \
                          @MFEM_SOURCE_DIR@/examples/ginkgo \
                          @MFEM_SOURCE_DIR@/examples/hiop \
+                         @MFEM_SOURCE_DIR@/examples/ipopt \
                          @MFEM_SOURCE_DIR@/examples/moonolith \
                          @MFEM_SOURCE_DIR@/examples/petsc \
                          @MFEM_SOURCE_DIR@/examples/pumi \

examples/CMakeLists.txt

@@ -178,6 +178,11 @@ if (MFEM_USE_HIOP)
   add_subdirectory(hiop)
 endif()
 
+# Include the examples/ipopt directory if IpOpt is enabled
+if (MFEM_USE_IPOPT)
+  add_subdirectory(ipopt)
+endif()
+
 # Include the examples/petsc directory if PETSc is enabled.
 if (MFEM_USE_PETSC)
   add_subdirectory(petsc)

examples/contact.cpp

@@ -0,0 +1,810 @@
+//                                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 "nodepair.hpp"
+
+using namespace std;
+using namespace mfem;
+
+bool ifequalarray(const Array<int> a1, const Array<int> a2)
+{
+   if (a1.Size()!=a2.Size())
+   {
+      return false;
+   }
+   for (int i=0; i<a1.Size(); i++)
+   {
+      if (a1[i] != a2[i])
+      {
+         return false;
+      }
+   }
+   return true;
+}
+
+void FindSurfaceToProject(Mesh& mesh, const int elem, int& cbdrface)
+{
+   Array<int> attr;
+   attr.Append(2);
+   Array<int> faces;
+   Array<int> ori;
+   std::vector<Array<int> > facesVertices;
+   std::vector<int > faceid;
+   mesh.GetElementFaces(elem, faces, ori);
+   int face = -1;
+   for (int i=0; i<faces.Size(); i++)
+   {
+      face = faces[i];
+      Array<int> faceVert;
+      if (!mesh.FaceIsInterior(face)) // if on the boundary
+      {
+         mesh.GetFaceVertices(face, faceVert);
+         faceVert.Sort();
+         facesVertices.push_back(faceVert);
+         faceid.push_back(face);
+      }
+   }
+   int bdrface = facesVertices.size();
+
+   Array<int> bdryFaces;
+   // This shoulnd't need to be rebuilt
+   std::vector<Array<int> > bdryVerts;
+   for (int b=0; b<mesh.GetNBE(); ++b)
+   {
+      if (attr.FindSorted(mesh.GetBdrAttribute(b)) >= 0)  // found the contact surface
+      {
+         bdryFaces.Append(b);
+         Array<int> vert;
+         mesh.GetBdrElementVertices(b, vert);
+         vert.Sort();
+         bdryVerts.push_back(vert);
+      }
+   }
+
+   int bdrvert = bdryVerts.size();
+   cbdrface = -1;  // the face number of the contact surface element
+   int count_cbdrface = 0;  // the number of matching surfaces, used for checks
+
+   for (int i=0; i<bdrface; i++)
+   {
+      for (int j=0; j<bdrvert; j++)
+      {
+         if (ifequalarray(facesVertices[i], bdryVerts[j]))
+         {
+            cbdrface = faceid[i];
+            count_cbdrface += 1;
+         }
+      }
+   }
+   MFEM_VERIFY(count_cbdrface == 1,"projection surface not found");
+
+};
+
+Vector GetNormalVector(Mesh & mesh, const int elem, const double *ref,
+                       int & refFace, int & refNormal, bool & interior)
+{
+   ElementTransformation *trans = mesh.GetElementTransformation(elem);
+   const int dim = mesh.Dimension();
+   const int spaceDim = trans->GetSpaceDim();
+
+   MFEM_VERIFY(spaceDim == 3, "");
+
+   Vector n(spaceDim);
+
+   IntegrationPoint ip;
+   ip.Set(ref, dim);
+
+   trans->SetIntPoint(&ip);
+   //CalcOrtho(trans->Jacobian(), n);  // Works only for face transformations
+   const DenseMatrix jac = trans->Jacobian();
+
+   int dimNormal = -1;
+   int normalSide = -1;
+
+   const double tol = 1.0e-8;
+   for (int i=0; i<dim; ++i)
+   {
+      const double d0 = std::abs(ref[i]);
+      const double d1 = std::abs(ref[i] - 1.0);
+
+      const double d = std::min(d0, d1);
+      // TODO: this works only for hexahedral meshes!
+
+      if (d < tol)
+      {
+         MFEM_VERIFY(dimNormal == -1, "");
+         dimNormal = i;
+
+         if (d0 < tol)
+         {
+            normalSide = 0;
+         }
+         else
+         {
+            normalSide = 1;
+         }
+      }
+   }
+   // closest point on the boundary
+   if (dimNormal < 0 || normalSide < 0) // node is inside the element
+   {
+      interior = 1;
+      Vector n(3);
+      n = 0.0;
+      return n;
+   }
+
+   MFEM_VERIFY(dimNormal >= 0 && normalSide >= 0, "");
+   refNormal = dimNormal;
+
+   MFEM_VERIFY(dim == 3, "");
+
+   {
+      // Find the reference face
+      if (dimNormal == 0)
+      {
+         refFace = (normalSide == 1) ? 2 : 4;
+      }
+      else if (dimNormal == 1)
+      {
+         refFace = (normalSide == 1) ? 3 : 1;
+      }
+      else
+      {
+         refFace = (normalSide == 1) ? 5 : 0;
+      }
+   }
+
+   std::vector<Vector> tang(2);
+
+   int tangDir[2] = {-1, -1};
+   {
+      int t = 0;
+      for (int i=0; i<dim; ++i)
+      {
+         if (i != dimNormal)
+         {
+            tangDir[t] = i;
+            t++;
+         }
+      }
+
+      MFEM_VERIFY(t == 2, "");
+   }
+
+   for (int i=0; i<2; ++i)
+   {
+      tang[i].SetSize(3);
+
+      Vector tangRef(3);
+      tangRef = 0.0;
+      tangRef[tangDir[i]] = 1.0;
+
+      jac.Mult(tangRef, tang[i]);
+   }
+
+   Vector c(3);  // Cross product
+
+   c[0] = (tang[0][1] * tang[1][2]) - (tang[0][2] * tang[1][1]);
+   c[1] = (tang[0][2] * tang[1][0]) - (tang[0][0] * tang[1][2]);
+   c[2] = (tang[0][0] * tang[1][1]) - (tang[0][1] * tang[1][0]);
+
+   c /= c.Norml2();
+
+   Vector nref(3);
+   nref = 0.0;
+   nref[dimNormal] = 1.0;
+
+   Vector ndir(3);
+   jac.Mult(nref, ndir);
+
+   ndir /= ndir.Norml2();
+
+   const double dp = ndir * c;
+
+   // TODO: eliminate c?
+   n = c;
+   if (dp < 0.0)
+   {
+      n *= -1.0;
+   }
+   interior = 0;
+   return n;
+}
+
+// WARNING: global variable, just for this little example.
+std::array<std::array<int, 3>, 8> HEX_VERT =
+{
+   {  {0,0,0},
+      {1,0,0},
+      {1,1,0},
+      {0,1,0},
+      {0,0,1},
+      {1,0,1},
+      {1,1,1},
+      {0,1,1}
+   }
+};
+
+int GetHexVertex(int cdim, int c, int fa, int fb, Vector & refCrd)
+{
+   int ref[3];
+   ref[cdim] = c;
+   ref[cdim == 0 ? 1 : 0] = fa;
+   ref[cdim == 2 ? 1 : 2] = fb;
+
+   for (int i=0; i<3; ++i) { refCrd[i] = ref[i]; }
+
+   int refv = -1;
+
+   for (int i=0; i<8; ++i)
+   {
+      bool match = true;
+      for (int j=0; j<3; ++j)
+      {
+         if (ref[j] != HEX_VERT[i][j]) { match = false; }
+      }
+
+      if (match) { refv = i; }
+   }
+
+   MFEM_VERIFY(refv >= 0, "");
+
+   return refv;
+}
+
+// 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)
+{
+   const int dim = mesh.Dimension();
+   const int np = xyz.Size() / dim;
+
+   MFEM_VERIFY(np * dim == xyz.Size(), "");
+
+   mesh.EnsureNodes();
+
+   //FindPointsGSLIB finder(MPI_COMM_WORLD);
+   FindPointsGSLIB finder;
+
+   finder.SetDistanceToleranceForPointsFoundOnBoundary(0.5);
+
+   const double bb_t = 0.5;
+   finder.Setup(mesh, bb_t);
+
+   finder.FindPoints(xyz);
+
+   /// Return code for each point searched by FindPoints: inside element (0), on
+   /// element boundary (1), or not found (2).
+   Array<unsigned int> codes = finder.GetCode();
+
+   /// Return element number for each point found by FindPoints.
+   Array<unsigned int> elems = finder.GetElem();
+
+   /// Return reference coordinates for each point found by FindPoints.
+   Vector refcrd = finder.GetReferencePosition();
+
+   /// Return distance between the sought and the found point in physical space,
+   /// for each point found by FindPoints.
+   Vector dist = finder.GetDist();
+
+   MFEM_VERIFY(dist.Size() == np, "");
+   MFEM_VERIFY(refcrd.Size() == np * dim, "");
+   MFEM_VERIFY(elems.Size() == np, "");
+   MFEM_VERIFY(codes.Size() == np, "");
+
+   bool allfound = true;
+   for (auto code : codes)
+      if (code == 2) { allfound = false; }
+
+   MFEM_VERIFY(allfound, "A point was not found");
+
+   cout << "Maximum distance of projected points: " << dist.Max() << endl;
+
+   // extract information
+   for (int i=0; i<np; ++i)
+   {
+      /*cout << "Point " << i << ": (";
+      for (int j=0; j<dim; ++j)
+      {
+         cout << xyz[i + (j*np)];
+         if (j == dim-1) {cout << ")" << endl;}
+         else{cout << ", ";}
+      }*/
+      //cout << "  element: " << elems[i] << endl;
+      //cout << "  element " << elems[i] << " vertices:" << endl;
+      //Array<int> vert;
+      //mesh.GetElementVertices(elems[i], vert);
+      //for (auto v : vert)
+      //{
+      //   cout << "    " << v << endl;
+      //}
+
+      /*cout << "  reference coordinates: (";
+      for (int j=0; j<dim; ++j)
+      {
+         cout << refcrd[(i*dim) + j];
+         if (j == dim-1)
+         {
+            cout << ")" << endl;
+         }
+         else
+         {
+            cout << ", ";
+         }
+      }*/
+
+      int refFace, refNormal, refNormalSide;
+      bool is_interior = -1;
+      Vector normal = GetNormalVector(mesh, elems[i], refcrd.GetData() + (i*dim),
+                                      refFace, refNormal, is_interior);
+      int phyFace;
+      if (is_interior)
+      {
+         phyFace = -1; // the id of the face that has the closest point
+         FindSurfaceToProject(mesh, elems[i], phyFace);
+
+         Array<int> cbdrVert;
+         mesh.GetFaceVertices(phyFace, cbdrVert);
+         Vector xs(dim);
+         xs[0] = xyz[i + 0*np];
+         xs[1] = xyz[i + 1*np];
+         xs[2] = xyz[i + 2*np];
+         Vector xi_tmp(dim-1);
+         // get nodes!
+
+         GridFunction *nodes = mesh.GetNodes();
+         DenseMatrix coords(4,3);
+         for (int i=0; i<4; i++)
+         {
+            for (int j=0; j<3; j++)
+            {
+               coords(i,j) = (*nodes)[cbdrVert[i]*3+j];
+            }
+         }
+         SlaveToMaster(coords, xs, xi_tmp);
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi[i*(dim-1)+j] = xi_tmp[j];
+         }
+         // now get get the projection to the surface
+      }
+      else
+      {
+         Vector faceRefCrd(dim-1);
+         {
+            int fd = 0;
+            for (int j=0; j<dim; ++j)
+            {
+               if (j == refNormal)
+               {
+                  refNormalSide = (refcrd[(i*dim) + j] > 0.5);
+               }
+               else
+               {
+                  faceRefCrd[fd] = refcrd[(i*dim) + j];
+                  fd++;
+               }
+            }
+
+            MFEM_VERIFY(fd == dim-1, "");
+         }
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi[i*(dim-1)+j] = faceRefCrd[j]*2.0 - 1.0;
+         }
+         //cout << "  face reference coordinates: (";
+         /*for (int j=0; j<dim-1; ++j)
+         {
+            cout << faceRefCrd[j];
+            if (j == dim-2){cout << ")" << endl;}
+            else{cout << ", ";}
+         }*/
+      }
+      //cout << "  normal vector: ";
+      //normal.Print();
+
+      // ask, does this do anything?
+      /*
+      IntegrationPoint ip;
+      ip.Set(refcrd.GetData() + (i*dim), dim);
+      ElementTransformation *trans = mesh.GetElementTransformation(elems[i]);
+      Vector phys(trans->GetSpaceDim());
+      trans->Transform(ip, phys);
+      cout << "  physical coordinates: ";
+      phys.Print();
+      */
+
+      // Get the element face
+      Array<int> faces;
+      Array<int> ori;
+      int face;
+
+      if (is_interior)
+      {
+         face = phyFace;
+      }
+      else
+      {
+         mesh.GetElementFaces(elems[i], faces, ori);
+         face = faces[refFace];
+      }
+
+      Array<int> faceVert;
+      mesh.GetFaceVertices(face, faceVert);
+
+      //cout << "  face " << face << " vertices:" << endl;
+      //for (auto v : faceVert){ cout << "    " << v << endl;}
+
+      for (int p=0; p<4; p++)
+      {
+         conn[4*i+p] = faceVert[p];
+      }
+      /*
+      Vector ref(dim);
+
+      for (int p=0; p<2; ++p)
+         for (int q=0; q<2; ++q)
+         {
+            const int refv = GetHexVertex(refNormal, refNormalSide, p, q, ref);
+            cout << "  face reference vertex (" << p << "," << q
+                 << ") is global vertex " << vert[refv] << endl;
+
+            {
+               // Sanity check
+               ip.Set(ref.GetData(), dim);
+               trans->Transform(ip, phys);
+               for (int j=0; j<dim; ++j)
+               {
+                  phys[j] -= mesh.GetVertex(vert[refv])[j];
+               }
+          phys.Print();
+          cout<<vert[refv]<<endl;
+               cout<<mesh.GetVertex(vert[refv])[0]<<endl;
+               cout<<mesh.GetVertex(vert[refv])[1]<<endl;
+               cout<<mesh.GetVertex(vert[refv])[2]<<endl;
+               MFEM_VERIFY(phys.Norml2() < 1.0e-12, "Sanity check failed");
+            }
+         }*/
+   }
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file1 = "block1.mesh";
+   const char *mesh_file2 = "block2.mesh";
+
+   Array<int> attr;
+   Array<int> m_attr;
+
+   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.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh1(mesh_file1, 1, 1);
+   Mesh mesh2(mesh_file2, 1, 1);
+
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mesh1a_sock(vishost, visport);
+      mesh1a_sock.precision(8);
+      mesh1a_sock << "mesh\n" << mesh1 << flush;
+      socketstream mesh2a_sock(vishost, visport);
+      mesh2a_sock.precision(8);
+      mesh2a_sock << "mesh\n" << mesh2 << flush;
+   }
+
+   const int dim = mesh1.Dimension();
+   MFEM_VERIFY(dim == mesh2.Dimension(), "");
+
+   // boundary attribute 2 is the potential contact surface of nodes
+   attr.Append(2);
+   // boundary attribute 2 is the potential contact surface for master surface
+   m_attr.Append(2);
+
+   //  Define a finite element space on the mesh. Here we use vector finite
+   //  elements, i.e. dim copies of a scalar finite element space. The vector
+   //  dimension is specified by the last argument of the FiniteElementSpace
+   //  constructor.
+   FiniteElementCollection *fec1;
+   FiniteElementSpace *fespace1;
+   fec1 = new H1_FECollection(1, dim);
+   fespace1 = new FiniteElementSpace(&mesh1, fec1, dim, Ordering::byVDIM);
+   cout << "Number of finite element unknowns for mesh1: "
+        << fespace1->GetTrueVSize() << endl;
+   mesh1.SetNodalFESpace(fespace1);
+   GridFunction nodes0 = *mesh1.GetNodes(); // undeformed mesh1 nodal grid function
+   GridFunction *nodes1 = mesh1.GetNodes();
+
+   FiniteElementCollection *fec2 = new H1_FECollection(1, dim);
+   FiniteElementSpace *fespace2 = new FiniteElementSpace(&mesh2, fec2, dim,
+                                                         Ordering::byVDIM);
+   cout << "Number of finite element unknowns for mesh2: "
+        << fespace2->GetTrueVSize() << endl;
+
+   // degrees of freedom of both meshes
+   int ndof_1 = fespace1->GetTrueVSize();
+   int ndof_2 = fespace2->GetTrueVSize();
+   int ndofs = ndof_1 + ndof_2;
+   // number of nodes for each mesh
+   int nnd_1 = mesh1.GetNV();
+   int nnd_2 = mesh2.GetNV();
+   int nnd = nnd_1 + nnd_2;
+   // Determine the list of true (i.e. conforming) essential boundary dofs.
+   // In this example, the boundary conditions are defined by marking only
+   // boundary attribute 1 from the mesh as essential and converting it to a
+   // list of true dofs.
+   Array<int> ess_tdof_list1, ess_bdr1(mesh1.bdr_attributes.Max());
+   ess_bdr1 = 0;
+   //ess_bdr1[0] = 1;
+   // Not ready to be passed on yet
+   // fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+   Array<int> ess_tdof_list2, ess_bdr2(mesh2.bdr_attributes.Max());
+   ess_bdr2 = 0;
+   //ess_bdr2[0] = 1;
+
+   // Define the displacement vector x as a finite element grid function
+   // corresponding to fespace. GridFunction is a derived class of Vector.
+   GridFunction x1(fespace1);
+   x1 = 0.0;
+   GridFunction x2(fespace2);
+   x2 = 0.0;
+
+   // Generate force
+   LinearForm *b1 = new LinearForm(fespace1);
+   b1->Assemble();
+
+   LinearForm *b2 = new LinearForm(fespace2);
+   b2->Assemble();
+
+   // Set up the bilinear form a(.,.) on the finite element space
+   //  corresponding to the linear elasticity integrator with piece-wise
+   //  constants coefficient lambda and mu.
+   Vector lambda1(mesh1.attributes.Max());
+   lambda1 = 57.6923076923;
+   PWConstCoefficient lambda1_func(lambda1);
+   Vector mu1(mesh1.attributes.Max());
+   mu1 = 38.4615384615;
+   PWConstCoefficient mu1_func(mu1);
+
+   BilinearForm *a1 = new BilinearForm(fespace1);
+   a1->AddDomainIntegrator(new ElasticityIntegrator(lambda1_func,mu1_func));
+
+   Vector lambda2(mesh2.attributes.Max());
+   lambda2 = 57.6923076923;
+   PWConstCoefficient lambda2_func(lambda2);
+   Vector mu2(mesh2.attributes.Max());
+   mu2 = 38.4615384615;
+   PWConstCoefficient mu2_func(mu2);
+
+   BilinearForm *a2 = new BilinearForm(fespace2);
+   a2->AddDomainIntegrator(new ElasticityIntegrator(lambda2_func,mu2_func));
+
+   a1->Assemble();
+   SparseMatrix A1;
+   Vector B1, X1;
+   a1->FormLinearSystem(ess_tdof_list1, x1, *b1, A1, X1, B1);
+
+   a2->Assemble();
+   SparseMatrix A2;
+   Vector B2, X2;
+   a2->FormLinearSystem(ess_tdof_list2, x2, *b2, A2, X2, B2);
+
+   // Combine elasticity operator for two meshes into one.
+   // Block Matrix
+   SparseMatrix K(ndofs,ndofs);
+   for (int i=0; i<A1.Height(); i++)
+   {
+      Array<int> col_tmp;
+      Vector v_tmp;
+      col_tmp = 0;
+      v_tmp = 0.0;
+      A1.GetRow(i, col_tmp, v_tmp);
+      K.SetRow(i, col_tmp, v_tmp);
+   }
+   for (int i=0; i<A2.Height(); i++)
+   {
+      Array<int> col_tmp;
+      Vector v_tmp;
+      col_tmp = 0;
+      v_tmp = 0.0;
+      A2.GetRow(i, col_tmp, v_tmp);
+      for (int j=0; j<col_tmp.Size(); j++)
+      {
+         col_tmp[j] += ndof_1;
+      }
+      K.SetRow(i+ndof_1, col_tmp, v_tmp);  // mesh1 top left corner
+   }
+
+   // Construct node to segment contact constraint.
+
+   attr.Sort();
+   cout << "Boundary attributes for contact surface faces in mesh 2" << endl;
+   for (auto a : attr) { cout << a << endl; }
+
+   Array<int> bdryFaces2;  // TODO: remove this?
+
+   std::set<int> bdryVerts2;
+   for (int b=0; b<mesh2.GetNBE(); ++b)
+   {
+      if (attr.FindSorted(mesh2.GetBdrAttribute(b)) >= 0)
+      {
+         bdryFaces2.Append(b);
+         Array<int> vert;
+         mesh2.GetBdrElementVertices(b, vert);
+         for (auto v : vert)
+         {
+            bdryVerts2.insert(v);
+         }
+      }
+   }
+
+   int npoints = bdryVerts2.size();
+   Array<int> s_conn(npoints); // connectivity of the second/slave mesh
+   Vector xyz(dim * npoints);
+   xyz = 0.0;
+
+   cout << "Boundary vertices for contact surface vertices in mesh 2" << endl;
+
+   // construct the nodal coordinates on mesh2 to be projected, including displacement
+   int count = 0;
+   for (auto v : bdryVerts2)
+   {
+      cout << v << ": " << mesh2.GetVertex(v)[0] << ", "
+           << mesh2.GetVertex(v)[1] << ", "
+           << mesh2.GetVertex(v)[2] << endl;
+
+      for (int i=0; i<dim; ++i)
+      {
+         xyz[count + (i * npoints)] = mesh2.GetVertex(v)[i] + x2[v*dim+i];
+      }
+
+      s_conn[count] = v + nnd_1; // dof1 is the master
+      count++;
+   }
+
+   MFEM_VERIFY(count == npoints, "");
+
+   // gap function
+   Vector g(npoints*dim);
+   g = -1.0;
+   // segment reference coordinates of the closest point
+   Vector m_xi(npoints*(dim-1));
+   m_xi = -1.0;
+   Vector xs(dim*npoints);
+   xs = 0.0;
+   for (int i=0; i<npoints; i++)
+   {
+      for (int j=0; j<dim; j++)
+      {
+         xs[i*dim+j] = xyz[i + (j*npoints)];
+      }
+   }
+
+   Array<int> m_conn(
+      npoints*4); // only works for linear elements that have 4 vertices!
+   DenseMatrix coordsm(npoints*4, dim);
+
+   // adding displacement to mesh1 using a fixed grid function from mesh1
+   x1 = 1e-4; // x1 order: [xyz xyz... xyz]
+   add(nodes0, x1, *nodes1);
+
+   FindPointsInMesh(mesh1, xyz, m_conn, m_xi);
+
+   for (int i=0; i<npoints; i++)
+   {
+      for (int j=0; j<4; j++)
+      {
+         for (int k=0; k<dim; k++)
+         {
+            coordsm(i*4+j,k) = mesh1.GetVertex(m_conn[i*4+j])[k]+x1[dim*m_conn[i*4+j]+k];
+         }
+      }
+   }
+   //coordsm.Print();
+   SparseMatrix M(nnd,ndofs);
+   std::vector<SparseMatrix> dM(nnd, SparseMatrix(ndofs,ndofs));
+
+   Assemble_Contact(nnd, npoints, ndofs, xs, m_xi, coordsm,
+                    s_conn, m_conn, g, M, dM);
+
+   std::set<int> dirbdryv2;
+   for (int b=0; b<mesh2.GetNBE(); ++b)
+   {
+      if (mesh2.GetBdrAttribute(b) == 1)
+      {
+         Array<int> vert;
+         mesh2.GetBdrElementVertices(b, vert);
+         for (auto v : vert)
+         {
+            dirbdryv2.insert(v);
+         }
+      }
+   }
+   std::set<int> dirbdryv1;
+   for (int b=0; b<mesh1.GetNBE(); ++b)
+   {
+      if (mesh1.GetBdrAttribute(b) == 1)
+      {
+         Array<int> vert;
+         mesh1.GetBdrElementVertices(b, vert);
+         for (auto v : vert)
+         {
+            dirbdryv1.insert(v);
+         }
+      }
+   }
+
+   Array<int> Dirichlet_dof;
+   Array<double> Dirichlet_val;
+
+   for (auto v : dirbdryv2)
+   {
+      for (int i=0; i<dim; ++i)
+      {
+         Dirichlet_dof.Append(v*dim + i + ndof_1);
+         Dirichlet_val.Append(0.);
+      }
+   }
+   double delta = 0.1;
+   for (auto v : dirbdryv1)
+   {
+      Dirichlet_dof.Append(v*dim + 0);
+      Dirichlet_val.Append(delta);
+      Dirichlet_dof.Append(v*dim + 1);
+      Dirichlet_val.Append(0.);
+      Dirichlet_dof.Append(v*dim + 2);
+      Dirichlet_val.Append(0.);
+   }
+
+
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mesh1_sock(vishost, visport);
+      mesh1_sock.precision(8);
+      mesh1_sock << "mesh\n" << mesh1 << flush;
+      socketstream mesh2_sock(vishost, visport);
+      mesh2_sock.precision(8);
+      mesh2_sock << "mesh\n" << mesh2 << flush;
+   }
+
+   //M.Print();
+   /*Vector eps(ndofs);
+   Vector sol(ndofs); sol = 0.;
+   for(int i=0;i<ndofs;i++) eps[i] = 1e-5 * i ;
+   for(int i=0;i<9;i++)
+   {
+     cout<<i<<endl;
+     dM[s_conn[i]].Mult(eps,sol);
+     sol.Print();
+   }
+   */
+   return 0;
+}

examples/ipopt/CMakeLists.txt

@@ -0,0 +1,60 @@
+# Copyright (c) 2010-2022, 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.
+
+set(IPOPT_EXAMPLES_SRCS)
+list(APPEND IPOPT_EXAMPLES_SRCS exContactBlockTL.cpp)
+
+# Include the source directory where mfem.hpp and mfem-performance.hpp are.
+include_directories(BEFORE ${PROJECT_BINARY_DIR})
+
+# Add "test_ipopt" target, see below.
+add_custom_target(test_ipopt
+  ${CMAKE_CTEST_COMMAND} -R ipopt USES_TERMINAL)
+
+# Add one executable per cpp file, adding "ipopt_" as prefix. Sets
+# "test_ipopt" as a target that depends on the given examples.
+set(PFX ipopt_)
+add_mfem_examples(IPOPT_EXAMPLES_SRCS ${PFX} "" test_ipopt)
+
+# Testing.
+# The IPOPT tests can be run separately using the target "test_ipopt"
+# which builds the examples and runs:
+#   ctest -R ipopt
+
+if (MFEM_ENABLE_TESTING)
+  # Command line options for the tests.
+  # Example 9:
+  set(EXCONTACTBTL_COMMON_OPTS -m ../../data/periodic-segment.mesh -p 0 -dt 0.005)
+  set(EXCONTACTBTL_TEST_OPTS ${EXCONTACTBTL_COMMON_OPTS}  -r  2 )
+
+  # Add the tests: one test per source file.
+  foreach(SRC_FILE ${IPOPT_EXAMPLES_SRCS})
+    get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
+    string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
+    string(TOUPPER ${TEST_NAME} UP_TEST_NAME)
+    set(TEST_NAME ${PFX}${TEST_NAME})
+
+    set(THIS_TEST_OPTIONS "-no-vis")
+    list(APPEND THIS_TEST_OPTIONS ${${UP_TEST_NAME}_TEST_OPTS})
+    # message(STATUS "Test ${TEST_NAME} options: ${THIS_TEST_OPTIONS}")
+
+    if (NOT (${TEST_NAME} MATCHES ".*p$"))
+      add_test(NAME ${TEST_NAME}_ser
+        COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
+    else()
+      add_test(NAME ${TEST_NAME}_np=4
+        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4
+        ${MPIEXEC_PREFLAGS}
+        $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS}
+        ${MPIEXEC_POSTFLAGS})
+    endif()
+  endforeach()
+endif()

examples/ipopt/README

@@ -0,0 +1,19 @@
+                    Finite Element Discretization Library
+                                   __
+                       _ __ ___   / _|  ___  _ __ ___
+                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
+                      | | | | | ||  _||  __/| | | | | |
+                      |_| |_| |_||_|   \___||_| |_| |_|
+
+                               https://mfem.org
+
+This directory contains modifications of the example codes that illustrate the
+use of MFEM for solving nonlinear constrained optimization problems, including
+features based on the IpOpt, a lightweight HPC solver for nonlinear optimization
+problems.
+
+To use the Ipopt features, make sure that MFEM is configured with the option
+"MFEM_USE_IPOPT = YES", see the top-level INSTALL file for details.
+
+We recommend comparing the original example codes with the corresponding files
+in the current directory.

examples/ipopt/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
+0   3    1     0     2     3
+0   3    3     2     4     5
+0   3    5     4     6     7
+0   3    24    25    27    26
+0   3    26    27    29    28
+0   3    28    29    31    30
+1   3    2     0     8    10
+1   3    4     2    10    12
+1   3    6     4    12    14
+1  3    10    8    16    18
+1  3    12    10    18    20
+1  3    14    12    20    22
+1  3    18    16    24    26
+1  3    20    18    26    28
+1  3    22    20    28    30
+2  3    1     3    11    9
+2  3     3     5    13    11
+2  3     5     7    15    13
+2  3    9    11    19    17
+2  3    11    13    21    19
+2  3    13    15    23    21
+2  3    17    19    27    25
+2  3    19    21    29    27
+2  3    21    23    31    29
+0  3     8     0     1    9
+0  3    16     8    9    17
+0  3    24    16    17    25
+0  3     6    14    15     7
+0  3    14    22    23    15
+0  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.6667         0
+         0    0.6667         0
+   -1.0000    1.0000         0
+         0    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.6667    0.3333
+         0    0.6667    0.3333
+   -1.0000    1.0000    0.3333
+         0    1.0000    0.3333
+   -1.0000         0    0.6667
+         0         0    0.6667
+   -1.0000    0.3333    0.6667
+         0    0.3333    0.6667
+   -1.0000    0.6667    0.6667
+         0    0.6667    0.6667
+   -1.0000    1.0000    0.6667
+         0    1.0000    0.6667
+   -1.0000         0    1.0000
+         0         0    1.0000
+   -1.0000    0.3333    1.0000
+         0    0.3333    1.0000
+   -1.0000    0.6667    1.0000
+         0    0.6667    1.0000
+   -1.0000    1.0000    1.0000
+         0    1.0000    1.0000

examples/ipopt/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
+0   3     1     0     2     3
+0   3     3     2     4     5
+0   3    12    13    15    14
+0   3    14    15    17    16
+2   3     2     0     6     8
+2   3     4     2     8    10
+2   3     8     6    12    14
+2   3    10     8    14    16
+1   3     1     3     9     7
+1   3     3     5    11    9
+1   3     7    9    15    13
+1   3     9    11    17    15
+0   3     6     0     1     7
+0   3    12     6     7    13
+0   3     4    10    11     5
+0   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

examples/ipopt/contact.cpp

@@ -0,0 +1,742 @@
+//                                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 "nodepair.hpp"
+
+using namespace std;
+using namespace mfem;
+
+bool ifequalarray(const Array<int> a1, const Array<int> a2)
+{
+   if (a1.Size()!=a2.Size())
+   {
+      return false;
+   }
+   for (int i=0; i<a1.Size(); i++)
+   {
+      if (a1[i] != a2[i])
+      {
+         return false;
+      }
+   }
+   return true;
+}
+
+void FindSurfaceToProject(Mesh& mesh, const int elem, int& cbdrface)
+{
+   Array<int> attr;
+   attr.Append(2);
+   Array<int> faces;
+   Array<int> ori;
+   std::vector<Array<int> > facesVertices;
+   std::vector<int > faceid;
+   mesh.GetElementFaces(elem, faces, ori);
+   int face = -1;
+   for (int i=0; i<faces.Size(); i++)
+   {
+      face = faces[i];
+      Array<int> faceVert;
+      if (!mesh.FaceIsInterior(face)) // if on the boundary
+      {
+         mesh.GetFaceVertices(face, faceVert);
+         faceVert.Sort();
+         facesVertices.push_back(faceVert);
+         faceid.push_back(face);
+      }
+   }
+   int bdrface = facesVertices.size();
+
+   Array<int> bdryFaces;
+   // This shoulnd't need to be rebuilt
+   std::vector<Array<int> > bdryVerts;
+   for (int b=0; b<mesh.GetNBE(); ++b)
+   {
+      if (attr.FindSorted(mesh.GetBdrAttribute(b)) >= 0)  // found the contact surface
+      {
+         bdryFaces.Append(b);
+         Array<int> vert;
+         mesh.GetBdrElementVertices(b, vert);
+         vert.Sort();
+         bdryVerts.push_back(vert);
+      }
+   }
+
+   int bdrvert = bdryVerts.size();
+   cbdrface = -1;  // the face number of the contact surface element
+   int count_cbdrface = 0;  // the number of matching surfaces, used for checks
+
+   for (int i=0; i<bdrface; i++)
+   {
+      for (int j=0; j<bdrvert; j++)
+      {
+         if (ifequalarray(facesVertices[i], bdryVerts[j]))
+         {
+            cbdrface = faceid[i];
+            count_cbdrface += 1;
+         }
+      }
+   }
+   MFEM_VERIFY(count_cbdrface == 1,"projection surface not found");
+
+};
+
+Vector GetNormalVector(Mesh & mesh, const int elem, const double *ref,
+                       int & refFace, int & refNormal, bool & interior)
+{
+   ElementTransformation *trans = mesh.GetElementTransformation(elem);
+   const int dim = mesh.Dimension();
+   const int spaceDim = trans->GetSpaceDim();
+
+   MFEM_VERIFY(spaceDim == 3, "");
+
+   Vector n(spaceDim);
+
+   IntegrationPoint ip;
+   ip.Set(ref, dim);
+
+   trans->SetIntPoint(&ip);
+   //CalcOrtho(trans->Jacobian(), n);  // Works only for face transformations
+   const DenseMatrix jac = trans->Jacobian();
+
+   int dimNormal = -1;
+   int normalSide = -1;
+
+   const double tol = 1.0e-8;
+   for (int i=0; i<dim; ++i)
+   {
+      const double d0 = std::abs(ref[i]);
+      const double d1 = std::abs(ref[i] - 1.0);
+
+      const double d = std::min(d0, d1);
+      // TODO: this works only for hexahedral meshes!
+
+      if (d < tol)
+      {
+         MFEM_VERIFY(dimNormal == -1, "");
+         dimNormal = i;
+
+         if (d0 < tol)
+         {
+            normalSide = 0;
+         }
+         else
+         {
+            normalSide = 1;
+         }
+      }
+   }
+   // closest point on the boundary
+   if (dimNormal < 0 || normalSide < 0) // node is inside the element
+   {
+      interior = 1;
+      Vector n(3);
+      n = 0.0;
+      return n;
+   }
+
+   MFEM_VERIFY(dimNormal >= 0 && normalSide >= 0, "");
+   refNormal = dimNormal;
+
+   MFEM_VERIFY(dim == 3, "");
+
+   {
+      // Find the reference face
+      if (dimNormal == 0)
+      {
+         refFace = (normalSide == 1) ? 2 : 4;
+      }
+      else if (dimNormal == 1)
+      {
+         refFace = (normalSide == 1) ? 3 : 1;
+      }
+      else
+      {
+         refFace = (normalSide == 1) ? 5 : 0;
+      }
+   }
+
+   std::vector<Vector> tang(2);
+
+   int tangDir[2] = {-1, -1};
+   {
+      int t = 0;
+      for (int i=0; i<dim; ++i)
+      {
+         if (i != dimNormal)
+         {
+            tangDir[t] = i;
+            t++;
+         }
+      }
+
+      MFEM_VERIFY(t == 2, "");
+   }
+
+   for (int i=0; i<2; ++i)
+   {
+      tang[i].SetSize(3);
+
+      Vector tangRef(3);
+      tangRef = 0.0;
+      tangRef[tangDir[i]] = 1.0;
+
+      jac.Mult(tangRef, tang[i]);
+   }
+
+   Vector c(3);  // Cross product
+
+   c[0] = (tang[0][1] * tang[1][2]) - (tang[0][2] * tang[1][1]);
+   c[1] = (tang[0][2] * tang[1][0]) - (tang[0][0] * tang[1][2]);
+   c[2] = (tang[0][0] * tang[1][1]) - (tang[0][1] * tang[1][0]);
+
+   c /= c.Norml2();
+
+   Vector nref(3);
+   nref = 0.0;
+   nref[dimNormal] = 1.0;
+
+   Vector ndir(3);
+   jac.Mult(nref, ndir);
+
+   ndir /= ndir.Norml2();
+
+   const double dp = ndir * c;
+
+   // TODO: eliminate c?
+   n = c;
+   if (dp < 0.0)
+   {
+      n *= -1.0;
+   }
+   interior = 0;
+   return n;
+}
+
+// WARNING: global variable, just for this little example.
+std::array<std::array<int, 3>, 8> HEX_VERT =
+{
+   {  {0,0,0},
+      {1,0,0},
+      {1,1,0},
+      {0,1,0},
+      {0,0,1},
+      {1,0,1},
+      {1,1,1},
+      {0,1,1}
+   }
+};
+
+int GetHexVertex(int cdim, int c, int fa, int fb, Vector & refCrd)
+{
+   int ref[3];
+   ref[cdim] = c;
+   ref[cdim == 0 ? 1 : 0] = fa;
+   ref[cdim == 2 ? 1 : 2] = fb;
+
+   for (int i=0; i<3; ++i) { refCrd[i] = ref[i]; }
+
+   int refv = -1;
+
+   for (int i=0; i<8; ++i)
+   {
+      bool match = true;
+      for (int j=0; j<3; ++j)
+      {
+         if (ref[j] != HEX_VERT[i][j]) { match = false; }
+      }
+
+      if (match) { refv = i; }
+   }
+
+   MFEM_VERIFY(refv >= 0, "");
+
+   return refv;
+}
+
+// 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)
+{
+   const int dim = mesh.Dimension();
+   const int np = xyz.Size() / dim;
+
+   MFEM_VERIFY(np * dim == xyz.Size(), "");
+
+   mesh.EnsureNodes();
+
+   //FindPointsGSLIB finder(MPI_COMM_WORLD);
+   FindPointsGSLIB finder;
+
+   finder.SetDistanceToleranceForPointsFoundOnBoundary(0.5);
+
+   const double bb_t = 0.5;
+   finder.Setup(mesh, bb_t);
+
+   finder.FindPoints(xyz);
+
+   /// Return code for each point searched by FindPoints: inside element (0), on
+   /// element boundary (1), or not found (2).
+   Array<unsigned int> codes = finder.GetCode();
+
+   /// Return element number for each point found by FindPoints.
+   Array<unsigned int> elems = finder.GetElem();
+
+   /// Return reference coordinates for each point found by FindPoints.
+   Vector refcrd = finder.GetReferencePosition();
+
+   /// Return distance between the sought and the found point in physical space,
+   /// for each point found by FindPoints.
+   Vector dist = finder.GetDist();
+
+   MFEM_VERIFY(dist.Size() == np, "");
+   MFEM_VERIFY(refcrd.Size() == np * dim, "");
+   MFEM_VERIFY(elems.Size() == np, "");
+   MFEM_VERIFY(codes.Size() == np, "");
+
+   bool allfound = true;
+   for (auto code : codes)
+      if (code == 2) { allfound = false; }
+
+   MFEM_VERIFY(allfound, "A point was not found");
+
+   cout << "Maximum distance of projected points: " << dist.Max() << endl;
+
+   // extract information
+   for (int i=0; i<np; ++i)
+   {
+      cout << "Point " << i << ": (";
+      for (int j=0; j<dim; ++j)
+      {
+         cout << xyz[i + (j*np)];
+         if (j == dim-1) {cout << ")" << endl;}
+         else {cout << ", ";}
+      }
+      //cout << "  element: " << elems[i] << endl;
+      //cout << "  element " << elems[i] << " vertices:" << endl;
+      //Array<int> vert;
+      //mesh.GetElementVertices(elems[i], vert);
+      //for (auto v : vert)
+      //{
+      //   cout << "    " << v << endl;
+      //}
+
+      /*cout << "  reference coordinates: (";
+      for (int j=0; j<dim; ++j)
+      {
+         cout << refcrd[(i*dim) + j];
+         if (j == dim-1)
+         {
+            cout << ")" << endl;
+         }
+         else
+         {
+            cout << ", ";
+         }
+      }*/
+
+      int refFace, refNormal, refNormalSide;
+      bool is_interior = -1;
+      Vector normal = GetNormalVector(mesh, elems[i], refcrd.GetData() + (i*dim),
+                                      refFace, refNormal, is_interior);
+      int phyFace;
+      if (is_interior)
+      {
+         phyFace = -1; // the id of the face that has the closest point
+         FindSurfaceToProject(mesh, elems[i], phyFace);
+
+         Array<int> cbdrVert;
+         mesh.GetFaceVertices(phyFace, cbdrVert);
+         Vector xs(dim);
+         xs[0] = xyz[i + 0*np];
+         xs[1] = xyz[i + 1*np];
+         xs[2] = xyz[i + 2*np];
+         Vector xi_tmp(dim-1);
+         // get nodes!
+
+         GridFunction *nodes = mesh.GetNodes();
+         DenseMatrix coords(4,3);
+         for (int i=0; i<4; i++)
+         {
+            for (int j=0; j<3; j++)
+            {
+               coords(i,j) = (*nodes)[cbdrVert[i]*3+j];
+            }
+         }
+         SlaveToMaster(coords, xs, xi_tmp);
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi[i*(dim-1)+j] = xi_tmp[j];
+         }
+         // now get get the projection to the surface
+      }
+      else
+      {
+         Vector faceRefCrd(dim-1);
+         {
+            int fd = 0;
+            for (int j=0; j<dim; ++j)
+            {
+               if (j == refNormal)
+               {
+                  refNormalSide = (refcrd[(i*dim) + j] > 0.5);
+               }
+               else
+               {
+                  faceRefCrd[fd] = refcrd[(i*dim) + j];
+                  fd++;
+               }
+            }
+
+            MFEM_VERIFY(fd == dim-1, "");
+         }
+
+         for (int j=0; j<dim-1; ++j)
+         {
+            xi[i*(dim-1)+j] = faceRefCrd[j]*2.0 - 1.0;
+         }
+         //cout << "  face reference coordinates: (";
+         for (int j=0; j<dim-1; ++j)
+         {
+            cout << faceRefCrd[j];
+            if (j == dim-2) {cout << ")" << endl;}
+            else {cout << ", ";}
+         }
+      }
+      //cout << "  normal vector: ";
+      //normal.Print();
+
+      // ask, does this do anything?
+      /*
+      IntegrationPoint ip;
+      ip.Set(refcrd.GetData() + (i*dim), dim);
+      ElementTransformation *trans = mesh.GetElementTransformation(elems[i]);
+      Vector phys(trans->GetSpaceDim());
+      trans->Transform(ip, phys);
+      cout << "  physical coordinates: ";
+      phys.Print();
+      */
+
+      // Get the element face
+      Array<int> faces;
+      Array<int> ori;
+      int face;
+
+      if (is_interior)
+      {
+         face = phyFace;
+      }
+      else
+      {
+         mesh.GetElementFaces(elems[i], faces, ori);
+         face = faces[refFace];
+      }
+
+      Array<int> faceVert;
+      mesh.GetFaceVertices(face, faceVert);
+
+      //cout << "  face " << face << " vertices:" << endl;
+      //for (auto v : faceVert){ cout << "    " << v << endl;}
+
+      for (int p=0; p<4; p++)
+      {
+         conn[4*i+p] = faceVert[p];
+      }
+      /*
+      Vector ref(dim);
+
+      for (int p=0; p<2; ++p)
+         for (int q=0; q<2; ++q)
+         {
+            const int refv = GetHexVertex(refNormal, refNormalSide, p, q, ref);
+            cout << "  face reference vertex (" << p << "," << q
+                 << ") is global vertex " << vert[refv] << endl;
+
+            {
+               // Sanity check
+               ip.Set(ref.GetData(), dim);
+               trans->Transform(ip, phys);
+               for (int j=0; j<dim; ++j)
+               {
+                  phys[j] -= mesh.GetVertex(vert[refv])[j];
+               }
+          phys.Print();
+          cout<<vert[refv]<<endl;
+               cout<<mesh.GetVertex(vert[refv])[0]<<endl;
+               cout<<mesh.GetVertex(vert[refv])[1]<<endl;
+               cout<<mesh.GetVertex(vert[refv])[2]<<endl;
+               MFEM_VERIFY(phys.Norml2() < 1.0e-12, "Sanity check failed");
+            }
+         }*/
+   }
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file1 = "block1.mesh";
+   const char *mesh_file2 = "block2.mesh";
+
+   Array<int> attr;
+   Array<int> m_attr;
+
+   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.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh1(mesh_file1, 1, 1);
+   Mesh mesh2(mesh_file2, 1, 1);
+
+   const int dim = mesh1.Dimension();
+   MFEM_VERIFY(dim == mesh2.Dimension(), "");
+
+   // boundary attribute 2 is the potential contact surface of nodes
+   attr.Append(2);
+   // boundary attribute 2 is the potential contact surface for master surface
+   m_attr.Append(2);
+
+   //  Define a finite element space on the mesh. Here we use vector finite
+   //  elements, i.e. dim copies of a scalar finite element space. The vector
+   //  dimension is specified by the last argument of the FiniteElementSpace
+   //  constructor.
+   FiniteElementCollection *fec1;
+   FiniteElementSpace *fespace1;
+   fec1 = new H1_FECollection(1, dim);
+   fespace1 = new FiniteElementSpace(&mesh1, fec1, dim, Ordering::byVDIM);
+   cout << "Number of finite element unknowns for mesh1: "
+        << fespace1->GetTrueVSize() << endl;
+   mesh1.SetNodalFESpace(fespace1);
+   GridFunction nodes0 = *mesh1.GetNodes(); // undeformed mesh1 nodal grid function
+   GridFunction *nodes1 = mesh1.GetNodes();
+
+   FiniteElementCollection *fec2 = new H1_FECollection(1, dim);
+   FiniteElementSpace *fespace2 = new FiniteElementSpace(&mesh2, fec2, dim,
+                                                         Ordering::byVDIM);
+   cout << "Number of finite element unknowns for mesh2: "
+        << fespace2->GetTrueVSize() << endl;
+
+   // degrees of freedom of both meshes
+   int ndof_1 = fespace1->GetTrueVSize();
+   int ndof_2 = fespace2->GetTrueVSize();
+   int ndofs = ndof_1 + ndof_2;
+   // number of nodes for each mesh
+   int nnd_1 = mesh1.GetNV();
+   int nnd_2 = mesh2.GetNV();
+   int nnd = nnd_1 + nnd_2;
+   // Determine the list of true (i.e. conforming) essential boundary dofs.
+   // In this example, the boundary conditions are defined by marking only
+   // boundary attribute 1 from the mesh as essential and converting it to a
+   // list of true dofs.
+   Array<int> ess_tdof_list1, ess_bdr1(mesh1.bdr_attributes.Max());
+   cout<<mesh1.bdr_attributes.Max()<<endl;
+   ess_bdr1 = 0;
+   //ess_bdr1[0] = 1;
+   // Not ready to be passed on yet
+   // fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+   Array<int> ess_tdof_list2, ess_bdr2(mesh2.bdr_attributes.Max());
+   ess_bdr2 = 0;
+   //ess_bdr2[0] = 1;
+
+   // Define the displacement vector x as a finite element grid function
+   // corresponding to fespace. GridFunction is a derived class of Vector.
+   GridFunction x1(fespace1);
+   x1 = 0.0;
+   GridFunction x2(fespace2);
+   x2 = 0.0;
+
+   // Generate force
+   LinearForm *b1 = new LinearForm(fespace1);
+   b1->Assemble();
+
+   LinearForm *b2 = new LinearForm(fespace2);
+   b2->Assemble();
+
+   // Set up the bilinear form a(.,.) on the finite element space
+   //  corresponding to the linear elasticity integrator with piece-wise
+   //  constants coefficient lambda and mu.
+   Vector lambda1(mesh1.attributes.Max());
+   lambda1 = 57.6923076923;
+   PWConstCoefficient lambda1_func(lambda1);
+   Vector mu1(mesh1.attributes.Max());
+   mu1 = 38.4615384615;
+   PWConstCoefficient mu1_func(mu1);
+
+   BilinearForm *a1 = new BilinearForm(fespace1);
+   a1->AddDomainIntegrator(new ElasticityIntegrator(lambda1_func,mu1_func));
+
+   Vector lambda2(mesh2.attributes.Max());
+   lambda2 = 57.6923076923;
+   PWConstCoefficient lambda2_func(lambda2);
+   Vector mu2(mesh2.attributes.Max());
+   mu2 = 38.4615384615;
+   PWConstCoefficient mu2_func(mu2);
+
+   BilinearForm *a2 = new BilinearForm(fespace2);
+   a2->AddDomainIntegrator(new ElasticityIntegrator(lambda2_func,mu2_func));
+
+   a1->Assemble();
+   SparseMatrix A1;
+   Vector B1, X1;
+   a1->FormLinearSystem(ess_tdof_list1, x1, *b1, A1, X1, B1);
+
+   a2->Assemble();
+   SparseMatrix A2;
+   Vector B2, X2;
+   a2->FormLinearSystem(ess_tdof_list2, x2, *b2, A2, X2, B2);
+
+   // Combine elasticity operator for two meshes into one.
+   // Block Matrix
+   SparseMatrix K(ndofs,ndofs);
+   for (int i=0; i<A1.Height(); i++)
+   {
+      Array<int> col_tmp;
+      Vector v_tmp;
+      col_tmp = 0;
+      v_tmp = 0.0;
+      A1.GetRow(i, col_tmp, v_tmp);
+      K.SetRow(i, col_tmp, v_tmp);
+   }
+   for (int i=0; i<A2.Height(); i++)
+   {
+      Array<int> col_tmp;
+      Vector v_tmp;
+      col_tmp = 0;
+      v_tmp = 0.0;
+      A2.GetRow(i, col_tmp, v_tmp);
+      for (int j=0; j<col_tmp.Size(); j++)
+      {
+         col_tmp[j] += ndof_1;
+      }
+      K.SetRow(i+ndof_1, col_tmp, v_tmp);  // mesh1 top left corner
+   }
+
+   // Construct node to segment contact constraint.
+
+   attr.Sort();
+   cout << "Boundary attributes for contact surface faces in mesh 2" << endl;
+   for (auto a : attr)
+   {
+      cout << a << endl;
+   }
+
+   Array<int> bdryFaces2;  // TODO: remove this?
+
+   std::set<int> bdryVerts2;
+   for (int b=0; b<mesh2.GetNBE(); ++b)
+   {
+      if (attr.FindSorted(mesh2.GetBdrAttribute(b)) >= 0)
+      {
+         bdryFaces2.Append(b);
+         Array<int> vert;
+         mesh2.GetBdrElementVertices(b, vert);
+         for (auto v : vert)
+         {
+            bdryVerts2.insert(v);
+         }
+      }
+   }
+
+   int npoints = bdryVerts2.size();
+   Array<int> s_conn(npoints); // connectivity of the second/slave mesh
+   Vector xyz(dim * npoints);
+   xyz = 0.0;
+
+   cout << "Boundary vertices for contact surface vertices in mesh 2" << endl;
+
+   // construct the nodal coordinates on mesh2 to be projected, including displacement
+   int count = 0;
+   for (auto v : bdryVerts2)
+   {
+      cout << v << ": " << mesh2.GetVertex(v)[0] << ", "
+           << mesh2.GetVertex(v)[1] << ", "
+           << mesh2.GetVertex(v)[2] << endl;
+
+      for (int i=0; i<dim; ++i)
+      {
+         xyz[count + (i * npoints)] = mesh2.GetVertex(v)[i] + x2[v*dim+i];
+      }
+
+      s_conn[count] = v + nnd_1; // dof1 is the master
+      count++;
+   }
+
+   MFEM_VERIFY(count == npoints, "");
+
+   // gap function
+   Vector g(npoints*dim);
+   g = -1.0;
+   // segment reference coordinates of the closest point
+   Vector m_xi(npoints*(dim-1));
+   m_xi = -1.0;
+   Vector xs(dim*npoints);
+   xs = 0.0;
+   for (int i=0; i<npoints; i++)
+   {
+      for (int j=0; j<dim; j++)
+      {
+         xs[i*dim+j] = xyz[i + (j*npoints)];
+      }
+   }
+
+   Array<int> m_conn(
+      npoints*4); // only works for linear elements that have 4 vertices!
+   DenseMatrix coordsm(npoints*4, dim);
+
+   // adding displacement to mesh1 using a fixed grid function from mesh1
+   x1 = 1e-4; // x1 order: [xyz xyz... xyz]
+   add(nodes0, x1, *nodes1);
+
+   FindPointsInMesh(mesh1, xyz, m_conn, m_xi);
+
+   for (int i=0; i<npoints; i++)
+   {
+      for (int j=0; j<4; j++)
+      {
+         for (int k=0; k<dim; k++)
+         {
+            coordsm(i*4+j,k) = mesh1.GetVertex(m_conn[i*4+j])[k]+x1[dim*m_conn[i*4+j]+k];
+         }
+      }
+   }
+   //coordsm.Print();
+   SparseMatrix M(nnd,ndofs);
+   std::vector<SparseMatrix> dM(nnd, SparseMatrix(ndofs,ndofs));
+
+   Assemble_Contact(nnd, npoints, ndofs, xs, m_xi, coordsm,
+                    s_conn, m_conn, g, M, dM);
+
+   //M.Print();
+   /*Vector eps(ndofs);
+   Vector sol(ndofs); sol = 0.;
+   for(int i=0;i<ndofs;i++) eps[i] = 1e-5 * i ;
+   for(int i=0;i<9;i++)
+   {
+     cout<<i<<endl;
+     dM[s_conn[i]].Mult(eps,sol);
+     sol.Print();
+   }
+   */
+   return 0;
+}

examples/ipopt/exContactBlockTL.hpp

@@ -0,0 +1,230 @@
+//                                Contact example
+//
+// Compile with: make exContactBlockTL
+//
+// Sample runs:  ./exContactBlockTL -m1 block1.mesh -m2 block2.mesh -at "5 6 7 8"
+// Sample runs:  ./exContactBlockTL -m1 block1_d.mesh -m2 block2_d.mesh -at "5 6 7 8"
+
+#ifndef EXCONTACTBLOCKTL_HPP
+#define EXCONTACTBLOCKTL_HPP
+
+#include "mfem.hpp"
+#include "IpTNLP.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace Ipopt;
+
+
+class ExContactBlockTL: public TNLP
+{
+public:
+   /** default constructor */
+   ExContactBlockTL(int argc, char *argv[]);
+
+   /** default destructor */
+   virtual ~ExContactBlockTL();
+
+   /**@name Overloaded from TNLP */
+   /** Method to return some info about the nlp */
+   virtual bool get_nlp_info(
+      Index&          n,
+      Index&          m,
+      Index&          nnz_jac_g,
+      Index&          nnz_h_lag,
+      IndexStyleEnum& index_style
+   );
+
+   /** Method to return the bounds for my problem */
+   virtual bool get_bounds_info(
+      Index   n,
+      Number* x_l,
+      Number* x_u,
+      Index   m,
+      Number* g_l,
+      Number* g_u
+   );
+
+   /** Method to return the starting point for the algorithm */
+   virtual bool get_starting_point(
+      Index   n,
+      bool    init_x,
+      Number* x,
+      bool    init_z,
+      Number* z_L,
+      Number* z_U,
+      Index   m,
+      bool    init_lambda,
+      Number* lambda
+   );
+
+   /** Method to return the objective value */
+   virtual bool eval_f(
+      Index         n,
+      const Number* x,
+      bool          new_x,
+      Number&       obj_value
+   );
+
+   /** Method to return the gradient of the objective */
+   virtual bool eval_grad_f(
+      Index         n,
+      const Number* x,
+      bool          new_x,
+      Number*       grad_f
+   );
+
+   /** Method to return the constraint residuals */
+   virtual bool eval_g(
+      Index         n,
+      const Number* x,
+      bool          new_x,
+      Index         m,
+      Number*       cons
+   );
+
+   /** Method to return:
+    *   1) The structure of the Jacobian (if "values" is NULL)
+    *   2) The values of the Jacobian (if "values" is not NULL)
+    */
+   virtual bool eval_jac_g(
+      Index         n,
+      const Number* x,
+      bool          new_x,
+      Index         m,
+      Index         nele_jac,
+      Index*        iRow,
+      Index*        jCol,
+      Number*       values
+   );
+
+   /** Method to return:
+    *   1) The structure of the Hessian of the Lagrangian (if "values" is NULL)
+    *   2) The values of the Hessian of the Lagrangian (if "values" is not NULL)
+    */
+   virtual bool eval_h(
+      Index         n,
+      const Number* x,
+      bool          new_x,
+      Number        obj_factor,
+      Index         m,
+      const Number* lambda,
+      bool          new_lambda,
+      Index         nele_hess,
+      Index*        iRow,
+      Index*        jCol,
+      Number*       values
+   );
+
+   /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
+   virtual void finalize_solution(
+      SolverReturn               status,
+      Index                      n,
+      const Number*              x,
+      const Number*              z_L,
+      const Number*              z_U,
+      Index                      m,
+      const Number*              g,
+      const Number*              lambda,
+      Number                     obj_value,
+      const IpoptData*           ip_data,
+      IpoptCalculatedQuantities* ip_cq
+   );
+
+private:
+   void update_g();
+   void update_jac();
+   void update_hess();
+
+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;
+   GridFunction* x1;
+   GridFunction* x2;
+   LinearForm* b1;
+   LinearForm* b2;
+   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;
+   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;
+
+   SparseMatrix* K;
+   mfem::Vector g;
+   mfem::Vector m_xi;
+   mfem::Vector xs;
+
+   Array<int> m_conn; // only works for linear elements that have 4 vertices!
+   DenseMatrix* coordsm;
+   SparseMatrix* M;
+
+   std::vector<SparseMatrix>* dM;
+
+   Array<int> Dirichlet_dof;
+   Array<double> Dirichlet_val;
+
+public: 
+   Mesh * GetMesh1() {return mesh1;}
+   Mesh * GetMesh2() {return mesh2;}
+
+};
+
+#endif

examples/ipopt/makefile

@@ -0,0 +1,68 @@
+# Copyright (c) 2010-2022, 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/ipopt/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = exContactBlockTL
+EXAMPLES = $(SEQ_EXAMPLES)
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+ifeq ($(MFEM_USE_IPOPT),NO)
+$(EXAMPLES):
+	$(error MFEM is not configured with IPOPT)
+endif
+
+MFEM_TESTS = EXAMPLES
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel example)
+%-test-seq: %
+	@$(call mfem-test,$<,, Serial example)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# 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 *~ $(SEQ_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -f exContactBlockTL.mesh exContactBlockTL-mesh.* exContactBlockTL-init.* exContactBlockTL-final.* ExampleContactBlockTL*

examples/ipopt/nodepair.hpp

@@ -0,0 +1,888 @@
+
+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;
+   double r = 1e10;
+   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;
+   if (gap_v.Norml2() < proj_max_gap && xi.Normlinf() <= off_el_xi)
+   {
+      pt_on_elem = true;
+   }
+
+   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 n = ndofs;
+   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/makefile

@@ -46,6 +46,9 @@ endif
 ifeq ($(MFEM_USE_HIOP),YES)
    SUBDIRS += hiop
 endif
+ifeq ($(MFEM_USE_IPOPT),YES)
+   SUBDIRS += ipopt
+endif
 ifeq ($(MFEM_USE_PETSC),YES)
    SUBDIRS += petsc
 endif

examples/nodepair.hpp

@@ -0,0 +1,888 @@
+
+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;
+   double r = 1e10;
+   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;
+   if (gap_v.Norml2() < proj_max_gap && xi.Normlinf() <= off_el_xi)
+   {
+      pt_on_elem = true;
+   }
+
+   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 n = ndofs;
+   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);
+   }
+};
+

general/version.cpp

@@ -97,6 +97,9 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_HIOP
       "MFEM_USE_HIOP\n"
 #endif
+#ifdef MFEM_USE_IPOPT
+      "MFEM_USE_IPOPT\n"
+#endif
 #ifdef MFEM_USE_HIP
       "MFEM_USE_HIP\n"
 #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 ipopt petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
@@ -275,7 +275,7 @@ endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = ENZYME SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
- SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP\
+ SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP IPOPT\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER PARELAG BENCHMARK\
  MOONOLITH ALGOIM
 
@@ -341,7 +341,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
- MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
+ MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_IPOPT MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX\
  MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
@@ -690,6 +690,7 @@ status info:
 	$(info MFEM_USE_CONDUIT       = $(MFEM_USE_CONDUIT))
 	$(info MFEM_USE_PUMI          = $(MFEM_USE_PUMI))
 	$(info MFEM_USE_HIOP          = $(MFEM_USE_HIOP))
+	$(info MFEM_USE_IPOPT          = $(MFEM_USE_IPOPT))
 	$(info MFEM_USE_GSLIB         = $(MFEM_USE_GSLIB))
 	$(info MFEM_USE_CUDA          = $(MFEM_USE_CUDA))
 	$(info MFEM_USE_HIP           = $(MFEM_USE_HIP))