minor bug

remove comments
fixing bugs
2023-05-31 11:00:58 -07:00 · 2023-04-07 18:43:48 -07:00 · 2023-04-07 18:43:09 -07:00 · 2023-04-04 16:40:28 -07:00 · 2023-04-04 12:14:57 -07:00 · 2023-04-04 00:09:53 -07:00
868 changed files with 6434 additions and 11494 deletions
@@ -31,7 +31,7 @@ install:
 - ps: ( get-content "GKlib\gk_arch.h") | % { If ($_.ReadCount -ge 52) {$_ -replace "#ifdef __MSC__","#ifdef DISABLE_THIS_ANCIENT_MSC_CHECK"} Else {$_} } | set-content "GKlib\gk_arch.h"
 - cmake -H. -Bbuild
 # -DCMAKE_BUILD_TYPE=Release
- cmake --build build --target metis -j 4
+- cmake --build build
 - set METIS_PATH=%cd%
 - cd ..

@@ -40,7 +40,7 @@ install:
 - 7z x v2.19.0.tar.gz -so | 7z x -si -ttar > nul
 - cd hypre-2.19.0/src
 - cmake -H. -Bbuild
- cmake --build build -j 4
+- cmake --build build
 - cmake --build build --target install
 - cd ../..

@@ -49,20 +49,14 @@ install:

 # MFEM
 before_build:
- cmake . -Bbuild_parallel -DMFEM_USE_MPI=TRUE -DMFEM_USE_METIS_5=TRUE -DHYPRE_DIR=%cd%\..\hypre-2.19.0\src\hypre -DMETIS_LIBRARIES=%METIS_PATH%\build\libmetis\Debug\metis.lib -DMETIS_INCLUDE_DIRS=%METIS_PATH%\include
-# - cmake . -Bbuild_serial -DMFEM_USE_MPI=FALSE
- cmake . -Bbuild_shared_serial -DMFEM_USE_MPI=FALSE -DBUILD_SHARED_LIBS=ON
+- cmake -H. -DCMAKE_INSTALL_PREFIX=install -Bbuild_parallel -DMFEM_USE_MPI=TRUE -DMFEM_USE_METIS_5=TRUE -DHYPRE_DIR=%cd%\..\hypre-2.19.0\src\hypre -DMETIS_LIBRARIES=%METIS_PATH%\build\libmetis\Debug\metis.lib -DMETIS_INCLUDE_DIRS=%METIS_PATH%\include
+- cmake -H. -DCMAKE_INSTALL_PREFIX=install -Bbuild_serial -DMFEM_USE_MPI=FALSE

 build_script:
- cmake --build build_parallel --config Debug -j 4
-# - cmake --build build_serial --config Release -j 4
-# - cmake --build build_serial --target exec --config Release -j 4
- cmake --build build_shared_serial --config Release -j 4
- cmake --build build_shared_serial --target exec --config Release -j 4
+- cmake --build build_parallel --config Release -j 4
+- cmake --build build_serial --config Release -j 4
+- cmake --build build_serial --target exec --config Release -j 4

 after_build:
-# - cd build_serial
-# - ctest -C Release --output-on-failure
-# - cd ..
- cd build_shared_serial
+- cd build_serial
 - ctest -C Release --output-on-failure
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -74,19 +74,12 @@ jobs:
            codecov: YES
          - os: windows-latest
            codecov: NO
-          - os: windows-latest
-            target: opt
-            mpi: par
-            config-opts: '-DBUILD_SHARED_LIBS=ON'
          - os: ubuntu-latest
            target: opt
            codecov: NO
            mpi: par
            build-system: cmake
            hypre-target: int32
-            # This option can be set to pass additional configuration options to
-            # the MFEM configuration command.
-            # config-opts: '-DCMAKE_VERBOSE_MAKEFILE=ON'
          - os: ubuntu-latest
            target: opt
            codecov: NO
@@ -217,7 +210,7 @@ jobs:

    # MFEM build and test
    - name: build
-      uses: mfem/github-actions/build-mfem@v2.3
+      uses: mfem/github-actions/build-mfem@v2.2
      env:
        VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
      with:
@@ -229,7 +222,7 @@ jobs:
        hypre-dir: ${{ env.HYPRE_TOP_DIR }}
        metis-dir: ${{ env.METIS_TOP_DIR }}
        mfem-dir: ${{ env.MFEM_TOP_DIR }}
-        config-options: ${{ matrix.config-opts }}
+        config-options: ${{ env.MFEM_EXTRA_CONFIG }}
        library-only: ${{ matrix.target == 'dbg' }}

    # Run checks (and only checks) on debug targets
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -21,9 +21,6 @@ CMakeFiles/
 # Default install location
 /mfem/

-# Typical build directory
-/build/
-
 # Generated files in main directory, config/ and docs/
 /deps.mk
 config/_config.hpp
@@ -121,9 +118,11 @@ examples/amgx/sol.gf
 examples/amgx/mesh.*
 examples/amgx/sol.*

-examples/ginkgo/ex1
-examples/ginkgo/refined.mesh
-examples/ginkgo/sol.gf
+examples/gingko/ex1
+examples/gingko/refined.mesh
+examples/gingko/sol.gf
+examples/gingko/mesh.*
+examples/gingko/sol.*

 examples/hiop/ex9
 examples/hiop/ex9p
@@ -132,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.*
@@ -149,8 +154,8 @@ examples/petsc/velocity.*
 examples/petsc/elastic_energy.*
 examples/petsc/mode_*

-examples/pumi/ex[12]
-examples/pumi/ex[16]p
+examples/pumi/ex1
+examples/pumi/ex[126]p
 examples/pumi/refined.mesh
 examples/pumi/sol.gf
 examples/pumi/mesh.*
@@ -195,9 +200,6 @@ miniapps/gslib/field-diff
 miniapps/gslib/field-interp
 miniapps/gslib/findpts
 miniapps/gslib/pfindpts
-miniapps/gslib/schwarz_ex1
-miniapps/gslib/schwarz_ex1p
-miniapps/gslib/interpolated.gf

 miniapps/meshing/mobius-strip
 miniapps/meshing/klein-bottle
@@ -210,7 +212,6 @@ miniapps/meshing/trimmer
 miniapps/meshing/reflector
 miniapps/meshing/mesh-optimizer
 miniapps/meshing/pmesh-optimizer
-miniapps/meshing/pmesh-fitting
 miniapps/meshing/minimal-surface
 miniapps/meshing/pminimal-surface
 miniapps/meshing/polar-nc
@@ -250,7 +251,6 @@ miniapps/navier/navier_tgv
 miniapps/navier/navier_shear
 miniapps/navier/navier_3dfoc
 miniapps/navier/navier_turbchan
-miniapps/navier/navier_cht
 miniapps/navier/tgv_out*.txt
 miniapps/navier/*_output

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -27,6 +27,11 @@ allocate_resource:
  timeout: 6h

 # GitLab jobs for the Quartz machine at LLNL
+debug_ser_gcc_4_9_3:
+  variables:
+    SPEC: "%gcc@4.9.3 +debug~mpi"
+  extends: .build_and_test_on_quartz
+
 debug_ser_gcc_6_1_0:
  variables:
    SPEC: "%gcc@6.1.0 +debug~mpi"
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -8,69 +8,61 @@
                               https://mfem.org


-Version 4.5.3 (development)
+Version 4.5.1 (development)
 ===========================
+- When using discontinuous (L2) spaces, use local (element-wise) L2 projection
+  as the coarsening operator for non-conforming AMR meshes.

- Added new methods in the Mesh class to set and get attributes on NURBS patches
-  and patch boundaries.
-
-New and updated examples and miniapps
-------------------------------------
- Added a miniapp pmesh-fitting in miniapps/meshing for interface and boundary fitting to implicit domains defined using level-set functions.
-
- Moved the distance solver methods from miniapps/shifted to miniapps/common. 
-
-Discretization improvements
---------------------------
- Face restriction operators for Nedelec and Raviart-Thomas finite element
-  spaces are now supported through the ConformingFaceRestriction class.
-
- VectorFEBoundaryFluxLFIntegrator is now supported on device/GPU.
-
-
-Version 4.5.2, released on March 23, 2023
-=========================================
-
+Meshing improvements
+--------------------
 - Added support for pyramids in non-conforming meshes. Currently only isotropic
  refinement is supported in this case.

- Removed the support for the Mesquite toolkit. We recommend using MFEM's TMOP
-  functionality instead for mesh optimization. See the mesh-optimizer miniapp.
-
- Added a fast normalization-based distance solver, see the Distance miniapp
-  in the miniapps/shifted/ directory.
-
- Added a new meshing miniapp, Reflector, which reflects a high-order or NURBS
-  hexahedral mesh about a plane.
-
 - Updated logic in FindPointsGSLIB to ignore points found near (but outside) the
  domain boundary.

- Added an option to auto-balance compound TMOP metrics.
+- Added support for pyramids in Gmsh meshes.

 - Fixed a bug in TMOP metric 301.

- When using discontinuous (L2) spaces, use local (element-wise) L2 projection
-  as the coarsening operator for non-conforming AMR meshes.
+- Added an option to auto-balance compound TMOP metrics.

- Added support for GridFunction::GetGradients() and GetVectorGradient() on
-  face-neighbor elements.
+Discretization improvements
+---------------------------
+- TBD

- Added support for pyramids in Gmsh meshes.
+Linear and nonlinear solvers
+----------------------------
+- Added a fast normalization-based distance solver, see the Distance miniapp
+  in the miniapps/shifted/ directory.

- The Mesh Explorer miniapp can now save mesh files in the VisIt or ParaView
-  formats using the corresponding DataCollection. See option 'D' in the menu.
+New and updated examples and miniapps
+-------------------------------------
+- Added a new meshing miniapp, reflector, which reflects a high-order or NURBS
+  hexahedral mesh about a plane.

- VisItDataCollection now correctly handles data collection names containing
-  underscores.
+- The mesh-explorer miniapp can now save mesh files in the VisIt or ParaView
+  formats using the corresponding DataCollection objects. See option 'D' in the
+  main menu.

- Added support for shared Windows builds with MSVC through CMake.
-  Developers note: this enhancement is facilitated by the addition of the macro
-  MFEM_EXPORT, see config.hpp for more details on its usage.
+Integrations, testing and documentation
+---------------------------------------
+- Removed the support for the Mesquite toolkit. We recommend using MFEM's TMOP
+  functionality instead for mesh optimization. See the mesh-optimizer miniapp.

 - The following integrations have updated minimum version requirements:
  * RAJA >= 2022.10.3

+Miscellaneous
+-------------
+- VisItDataCollection now correctly handles collection names containing
+  underscores.
+
+- VisItDataCollection::SetPadDigits() no longer alters the number of digits
+  used to represent the MPI rank because VisIt seems to require 6 digits.
+  This parameter can still be explicitly overridden with
+  VisItDataCollection::SetPadDigitsRank().
+
 API changes
 -----------
 - The implicit cast methods of class Vector to 'double *' and 'const double *'
@@ -81,11 +73,6 @@ API changes
  deprecated, and Mesh::GetFaceGeometry (which provides identical functionality)
  should be used instead.

- VisItDataCollection::SetPadDigits() no longer alters the number of digits
-  used to represent the MPI rank because VisIt seems to require 6 digits.
-  This parameter can still be explicitly overridden with SetPadDigitsRank().
-
-
 Version 4.5, released on October 22, 2022
 =========================================

@@ -129,6 +116,9 @@ Discretization improvements
 - Added a class CoefficientVector for efficient access of variable coefficient
  values at quadrature points (in particular for GPU/device kernels).

+- Added support for GridFunction::GetGradients() and
+  GriFunction::GetVectorGradient() on face-neighbor elements.
+
 - Added WhiteGaussianNoiseDomainLFIntegrator: a LinearFormIntegrator class for
  spatial Gaussian white noise.

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -13,7 +13,6 @@
 # Version 3.8 fixes the handling of CMAKE_CXX_STANDARD for try_compile.
 # Version 3.8 or newer is required for direct CUDA support.
 cmake_minimum_required(VERSION 3.8)
-message(STATUS "CMake version: ${CMAKE_VERSION}")
 set(USER_CONFIG "${CMAKE_CURRENT_SOURCE_DIR}/config/user.cmake" CACHE PATH
  "Path to optional user configuration file.")

@@ -25,11 +24,8 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 # Load user settings before the defaults - this way the defaults will not
 # overwrite the user set options. If the user has not set all options, we still
 # have the defaults.
-include("${USER_CONFIG}" OPTIONAL RESULT_VARIABLE USER_CONFIG_LOADED)
-if (USER_CONFIG_LOADED)
-  set(USER_CONFIG_LOADED "LOADED")
-endif()
-message(STATUS "Loading USER_CONFIG = ${USER_CONFIG} (${USER_CONFIG_LOADED})")
+message(STATUS "(optional) USER_CONFIG = ${USER_CONFIG}")
+include("${USER_CONFIG}" OPTIONAL)
 include("${CMAKE_CURRENT_SOURCE_DIR}/config/defaults.cmake")

 # Allow overwriting of the compiler by setting CXX/MPICXX on the command line or
@@ -57,7 +53,7 @@ project(mfem NONE)
 # Current version of MFEM, see also `makefile`.
 #   mfem_VERSION = (string)
 #   MFEM_VERSION = (int)   [automatically derived from mfem_VERSION]
-set(${PROJECT_NAME}_VERSION 4.5.3)
+set(${PROJECT_NAME}_VERSION 4.5.1)

 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -187,23 +183,6 @@ else()
  set(MFEM_DEBUG OFF)
 endif()

-# Shared build on Windows
-if (WIN32 AND BUILD_SHARED_LIBS)
-  # CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS works only with MSVC?
-  set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON CACHE BOOL
-      "Automatically export symbols for shared Windows DLL build")
-  # Place all .dll and .exe files in the same sub-directory, 'bin/<config>', of
-  # the build directory, so that the .exe files can find the .dll files at
-  # runtime:
-  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin")
-  # Note: CMake tests are run from the build-tree sub-directory where the test
-  # is defined, e.g. <build-dir>/tests/unit, so hard coded paths to meshes still
-  # work fine.
-  # Not strictly necessary:
-  # set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib")
-  # set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib")
-endif()
-
 # AMD HIP
 if (MFEM_USE_HIP)
  if (HIP_ARCH)
@@ -425,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)
@@ -520,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)

@@ -540,7 +528,6 @@ list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_INCLUDE_DIRS)
 # message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")

-message(STATUS "MFEM shared library: BUILD_SHARED_LIBS = ${BUILD_SHARED_LIBS}")
 message(STATUS "MFEM build type: CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 message(STATUS "MFEM version: v${MFEM_VERSION_STRING}")
 message(STATUS "MFEM git string: ${MFEM_GIT_STRING}")
@@ -596,19 +583,7 @@ target_include_directories(mfem
  PUBLIC
    $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
-    # The INSTALL_INTERFACE is handled below
    ${TPL_INCLUDE_DIRS})
-if (MFEM_USE_MPI)
-  if (MPI_CXX_COMPILE_FLAGS)
-    separate_arguments(MPI_CXX_COMPILE_ARGS UNIX_COMMAND
-      "${MPI_CXX_COMPILE_FLAGS}")
-    target_compile_options(mfem PUBLIC ${MPI_CXX_COMPILE_ARGS})
-  endif()
-  if (MPI_CXX_LINK_FLAGS)
-    target_link_libraries(mfem PUBLIC ${MPI_CXX_LINK_FLAGS})
-  endif()
-endif()
-
 set_target_properties(mfem PROPERTIES VERSION "${mfem_VERSION}")
 set_target_properties(mfem PROPERTIES SOVERSION "${mfem_VERSION}")

@@ -620,7 +595,6 @@ if (NOT ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}"))
 endif()

 # Generate configuration file in the build directory: config/_config.hpp.
-set(MFEM_SHARED_BUILD ${BUILD_SHARED_LIBS})
 configure_file(
  "${PROJECT_SOURCE_DIR}/config/cmake/config.hpp.in"
  "${PROJECT_BINARY_DIR}/config/_config.hpp")
@@ -740,7 +714,7 @@ set(INSTALL_LIB_DIR lib
 set(INSTALL_CMAKE_DIR lib/cmake/mfem
  CACHE PATH "Relative path for installing cmake config files.")

-target_include_directories(mfem BEFORE
+target_include_directories(mfem
  PUBLIC
    $<INSTALL_INTERFACE:${INSTALL_INCLUDE_DIR}>)

@@ -792,8 +766,10 @@ endif()

 # Install the configuration header files
 install(FILES ${PROJECT_BINARY_DIR}/config/_config.hpp
-  ${PROJECT_SOURCE_DIR}/config/config.hpp
-  ${PROJECT_SOURCE_DIR}/config/tconfig.hpp
+  DESTINATION ${INSTALL_INCLUDE_DIR}/mfem/config
+  RENAME config.hpp)
+
+install(FILES ${PROJECT_SOURCE_DIR}/config/tconfig.hpp
  DESTINATION ${INSTALL_INCLUDE_DIR}/mfem/config)

 # Package the whole thing up nicely
@@ -112,6 +112,7 @@ The MFEM source code has the following structure:
  │   ├── caliper
  │   ├── ginkgo
  │   ├── hiop
+  │   ├── ipopt
  │   ├── jupyter
  │   ├── moonolith
  │   ├── petsc
@@ -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
@@ -1,6 +1,6 @@
 BSD 3-Clause License

-Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC
+Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC
 All rights reserved.

 Redistribution and use in source and binary forms, with or without
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -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@)
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// 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.
 //
@@ -39,9 +39,6 @@
 // Description of the git commit used to build MFEM.
 #cmakedefine MFEM_GIT_STRING "@MFEM_GIT_STRING@"

-// Enable shared library build of MFEM.
-#cmakedefine MFEM_SHARED_BUILD
-
 // Build the parallel MFEM library.
 // Requires an MPI compiler, and the libraries HYPRE and METIS.
 #cmakedefine MFEM_USE_MPI
@@ -134,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
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -21,15 +21,8 @@ if(EXISTS ${CONDUIT_DIR}/include/conduit/conduit_relay_hdf5.hpp)
  # we only need  HDF5 if Conduit was built with HDF5 support
  set(Conduit_REQUIRED_PACKAGES "HDF5" CACHE STRING
      "Additional packages required by Conduit.")
-  # Suppress warning about HDF5_ROOT being set
-  if (POLICY CMP0074)
-    cmake_policy(SET CMP0074 NEW)
-  endif()
-  # HDF5_ROOT is needed in some cases, e.g. when HDF5_TARGET_NAMES is set and
-  # MFEM's FindHDF5.cmake is not used.
-  set(HDF5_ROOT ${HDF5_DIR} CACHE PATH "")
-else()
-  message(STATUS "Conduit Relay HDF5 Support is DISABLED")
+  else()
+    message(STATUS "Conduit Relay HDF5 Support is DISABLED")
 endif()

 include(MfemCmakeUtilities)
@@ -39,4 +32,4 @@ mfem_find_package(Conduit CONDUIT CONDUIT_DIR
  ADD_COMPONENT relay
    "include;include/conduit" conduit_relay.hpp "lib" conduit_relay
  ADD_COMPONENT blueprint
-    "include;include/conduit" conduit_blueprint.hpp "lib" conduit_blueprint)
+      "include;include/conduit" conduit_blueprint.hpp "lib" conduit_blueprint)
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -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.")
+
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -14,17 +14,11 @@
 #   - NETCDF_LIBRARIES
 #   - NETCDF_INCLUDE_DIRS

-# Suppress warning about HDF5_ROOT being set
-if (POLICY CMP0074)
-  cmake_policy(SET CMP0074 NEW)
-endif()
-
 include(MfemCmakeUtilities)

 # FindHDF5.cmake uses HDF5_ROOT, so we "translate" from the MFEM convention
-# (Needed in some cases, e.g. when HDF5_TARGET_NAMES is set and MFEM's
-# FindHDF5.cmake is not used.)
-set(HDF5_ROOT ${HDF5_DIR} CACHE PATH "")
+# (MFEM's FindHDF5.cmake does not need HDF5_ROOT)
+# set(HDF5_ROOT ${HDF5_DIR} CACHE PATH "")

 # We need to guard against the case where HDF5 was already found but without
 # the HL extensions (in which case mfem_find_package will treat the package
@@ -41,11 +35,6 @@ mfem_find_package(NetCDF NETCDF NETCDF_DIR "include" netcdf.h "lib" netcdf
 # The netcdf library will always be the first element of NETCDF_LIBRARIES
 # and we need to insert after that library but before the hdf5 library, so
 # position 1 is used
-# (Needed in some cases, e.g. when HDF5_TARGET_NAMES is set and MFEM's
-# FindHDF5.cmake is not used.)
-if (HDF5_C_LIBRARY_hdf5_hl)
-  list(INSERT NETCDF_LIBRARIES 1 ${HDF5_C_LIBRARY_hdf5_hl})
-  if (NOT NetCDF_FIND_QUIETLY)
-    message(STATUS "Updated NetCDF libraries: ${NETCDF_LIBRARIES}")
-  endif()
-endif()
+# (MFEM's FindHDF5.cmake does not set HDF5_C_LIBRARY_hdf5_hl and the HL library
+# is already added to NETCDF_LIBRARIES)
+# list(INSERT NETCDF_LIBRARIES 1 ${HDF5_C_LIBRARY_hdf5_hl})
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -103,6 +103,25 @@ macro(add_mfem_examples EXE_SRCS)
      ${MFEM_EXEC_PREREQUISITES_TARGET_NAME} ${EXE_PREREQUISITE})

    target_link_libraries(${EXE_NAME} mfem)
+    if (MFEM_USE_MPI)
+      # Not needed: (mfem already links with MPI_CXX_LIBRARIES)
+      # target_link_libraries(${EXE_NAME} ${MPI_CXX_LIBRARIES})
+
+      # Language-specific include directories:
+      if (MPI_CXX_INCLUDE_PATH)
+        target_include_directories(${EXE_NAME} PRIVATE "${MPI_CXX_INCLUDE_PATH}")
+      endif()
+      if (MPI_CXX_COMPILE_FLAGS)
+        separate_arguments(MPI_CXX_COMPILE_ARGS UNIX_COMMAND
+          "${MPI_CXX_COMPILE_FLAGS}")
+        target_compile_options(${EXE_NAME} PRIVATE ${MPI_CXX_COMPILE_ARGS})
+      endif()
+
+      if (MPI_CXX_LINK_FLAGS)
+        set_target_properties(${EXE_NAME} PROPERTIES
+          LINK_FLAGS "${MPI_CXX_LINK_FLAGS}")
+      endif()
+    endif()
  endforeach(SRC_FILE)
 endmacro()

@@ -153,6 +172,26 @@ macro(add_mfem_miniapp MFEM_EXE_NAME)
  if (EXTRA_DEFINES_LIST)
    target_compile_definitions(${MFEM_EXE_NAME} PRIVATE ${EXTRA_DEFINES_LIST})
  endif()
+
+  # Handle the MPI separately
+  if (MFEM_USE_MPI)
+    # Add MPI_CXX_LIBRARIES, in case this target does not link with mfem.
+    target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_LIBRARIES})
+
+    if (MPI_CXX_INCLUDE_PATH)
+      target_include_directories(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_INCLUDE_PATH})
+    endif()
+    if (MPI_CXX_COMPILE_FLAGS)
+      separate_arguments(MPI_CXX_COMPILE_ARGS UNIX_COMMAND
+        "${MPI_CXX_COMPILE_FLAGS}")
+      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_COMPILE_ARGS})
+    endif()
+
+    if (MPI_CXX_LINK_FLAGS)
+      set_target_properties(${MFEM_EXE_NAME} PROPERTIES
+        LINK_FLAGS "${MPI_CXX_LINK_FLAGS}")
+    endif()
+  endif()
 endmacro()


@@ -397,9 +436,9 @@ function(mfem_find_package Name Prefix DirVar IncSuffixes Header LibSuffixes
              break()
            endif()
          endforeach()
-          if (Required AND NOT ${ReqPack}_FOUND)
+          if (${Required} AND NOT ${ReqPack}_FOUND)
            message(FATAL_ERROR " *** Required package ${ReqPack} not found."
-              " Checked target names: ${ReqPack} ${${ReqPack}_TARGET_NAMES}")
+              "Checked target names: ${ReqPack} ${${ReqPack}_TARGET_NAMES}")
          endif()
        endif()
      endif()
@@ -465,18 +504,14 @@ function(mfem_find_package Name Prefix DirVar IncSuffixes Header LibSuffixes
          endif()
          get_target_property(IsImported ${TargetName} IMPORTED)
          if (IsImported)
-            get_target_property(ImpConfigs ${TargetName} IMPORTED_CONFIGURATIONS)
-            # message(STATUS "${ReqPack} imported configs: ${ImpConfigs}")
            set(ImportConfig ${${ReqPack}_IMPORT_CONFIG})
            if (NOT ImportConfig)
              set(ImportConfig RELEASE)
-              list(FIND ImpConfigs ${ImportConfig} _Index)
-              if (_Index EQUAL -1)
-                list(GET ImpConfigs 0 ImportConfig)
-              endif()
            endif()
            set(ImportConfigSuffix "_${ImportConfig}")
-            if ((NOT ImportConfig) OR ("${ImportConfig}" STREQUAL "NO_CONFIG"))
+            get_target_property(ImpConfigs ${TargetName} IMPORTED_CONFIGURATIONS)
+            list(FIND ImpConfigs ${ImportConfig} _Index)
+            if ((_Index EQUAL -1) OR ("${ImportConfig}" STREQUAL "NO_CONFIG"))
              set(ImportConfig "NO_CONFIG")
              set(ImportConfigSuffix "")
              # message(FATAL_ERROR " *** ${ReqPack}: configuration "
@@ -494,8 +529,7 @@ function(mfem_find_package Name Prefix DirVar IncSuffixes Header LibSuffixes
          else()
            # Set _Pack_LIBS from the target properties for ImportConfig
            foreach (_prop IMPORTED_LOCATION${ImportConfigSuffix}
-                IMPORTED_LINK_INTERFACE_LIBRARIES${ImportConfigSuffix}
-                INTERFACE_LINK_LIBRARIES)
+                IMPORTED_LINK_INTERFACE_LIBRARIES${ImportConfigSuffix})
              get_target_property(_value ${TargetName} ${_prop})
              if (_value)
                list(APPEND _Pack_LIBS ${_value})
@@ -746,7 +780,6 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
      list(APPEND CompileOpts "-isystem \"${SysDir}\"")
    endforeach()
  endif()
-  set(AddInterfaceLinkLibs TRUE)
  if ("${type}" STREQUAL "STATIC_LIBRARY")
    get_target_property(Location ${tgt} LOCATION)
    if (Location)
@@ -768,27 +801,11 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
      message(STATUS " *** Warning: [${tgt}] LOCATION not defined!")
    endif()
  elseif ("${type}" STREQUAL "INTERFACE_LIBRARY")
-    # The INTERFACE_LINK_LIBRARIES property is handled below.
-    # Other properties we may need to handle:
-    # INTERFACE_LINK_DEPENDS
-    # INTERFACE_LINK_DIRECTORIES
-    # INTERFACE_LINK_OPTIONS
-  else()
-    message(STATUS " *** Warning: [${tgt}] uses target type '${type}'"
-      " which is not supported!")
-    set(AddInterfaceLinkLibs FALSE)
-  endif()
-  if (AddInterfaceLinkLibs)
    get_target_property(Libs ${tgt} INTERFACE_LINK_LIBRARIES)
-    # message(STATUS "${tgt}[INTERFACE_LINK_LIBRARIES]: ${Libs}")
    if (Libs)
      foreach(Lib ${Libs})
        if (NOT (TARGET ${Lib}))
-          # message(STATUS "Lib = ${Lib}")
-          # Filter-out generator expressions
-          if (NOT ("${Lib}" MATCHES "^\\$"))
-            list(APPEND LinkOpts "${Lib}")
-          endif()
+          list(APPEND LinkOpts "${Lib}")
        else()
          mfem_get_target_options(${Lib} COpts LOpts)
          list(APPEND CompileOpts ${COpts})
@@ -796,6 +813,13 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
        endif()
      endforeach()
    endif()
+    # Other properties we may need to handle:
+    # INTERFACE_LINK_DEPENDS
+    # INTERFACE_LINK_DIRECTORIES
+    # INTERFACE_LINK_OPTIONS
+  else()
+    message(STATUS " *** Warning: [${tgt}] uses target type '${type}'"
+      " which is not supported!")
  endif()

  # Other potentially relevant properties:
@@ -819,9 +843,6 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
  set(${CompileOptsVar} "${CompileOpts}" PARENT_SCOPE)
  set(${LinkOptsVar} "${LinkOpts}" PARENT_SCOPE)

-  # message(STATUS "${tgt}[CompileOpts]: ${CompileOpts}")
-  # message(STATUS "${tgt}[LinkOpts]: ${LinkOpts}")
-
 endfunction(mfem_get_target_options)


@@ -848,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
@@ -873,7 +894,7 @@ function(mfem_export_mk_files)
         "${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
         MFEM_CXXFLAGS)
  set(MFEM_TPLFLAGS "")
-  foreach(dir ${TPL_INCLUDE_DIRS})
+  foreach(dir ${MFEM_TPL_INCLUDE_DIRS})
    set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
  endforeach()
  # TODO: MFEM_TPLFLAGS: add other TPL flags, in addition to the -I flags.
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// 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.
 //
@@ -30,30 +30,11 @@
 #endif

 // Windows specific options
-#if defined(_WIN32) && !defined(_USE_MATH_DEFINES)
+#ifdef _WIN32
+#ifndef _USE_MATH_DEFINES
 // Macro needed to get defines like M_PI from <cmath>. (Visual Studio C++ only?)
 #define _USE_MATH_DEFINES
 #endif
-// Macro MFEM_EXPORT: this macro is used when declaring exported global
-// variables and static class variables in public header files, e.g.:
-//    extern MFEM_EXPORT Geometry Geometries;
-//    static MFEM_EXPORT Device device_singleton;
-// In cases where a class contains multiple static variables, instead of marking
-// all such variables with MFEM_EXPORT, one can mark the class with MFEM_EXPORT,
-// e.g.:
-//    class MFEM_EXPORT MemoryManager ...
-// Note: MFEM's GitHub CI includes a shared MSVC build that will fail if a
-// variable that needs MFEM_EXPORT does not have it. However, builds with
-// optional external libraries are not tested and may require separate checks to
-// determine the necessity of MFEM_EXPORT.
-#if defined(_MSC_VER) && defined(MFEM_SHARED_BUILD)
-#ifdef mfem_EXPORTS
-#define MFEM_EXPORT __declspec(dllexport)
-#else
-#define MFEM_EXPORT __declspec(dllimport)
-#endif
-#else
-#define MFEM_EXPORT
 #endif
 // On Cygwin the option -std=c++11 prevents the definition of M_PI. Defining
 // the following macro allows us to get M_PI and some needed functions, e.g.
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// 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.
 //
@@ -39,9 +39,6 @@
 // Description of the git commit used to build MFEM.
 // #define MFEM_GIT_STRING "@MFEM_GIT_STRING@"

-// Enable shared library build of MFEM.
-// #define MFEM_SHARED_BUILD
-
 // Build the parallel MFEM library.
 // Requires an MPI compiler, and the libraries HYPRE and METIS.
 // #define MFEM_USE_MPI
@@ -144,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

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -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@
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -19,7 +19,6 @@ if (NOT CMAKE_BUILD_TYPE)
 endif()

 # MFEM options. Set to mimic the default "defaults.mk" file.
-option(BUILD_SHARED_LIBS "Enable shared library build of MFEM" OFF)
 option(MFEM_USE_MPI "Enable MPI parallel build" OFF)
 option(MFEM_USE_METIS "Enable METIS usage" ${MFEM_USE_MPI})
 option(MFEM_USE_EXCEPTIONS "Enable the use of exceptions" OFF)
@@ -49,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)
@@ -221,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")
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -58,7 +58,7 @@ HIP_CXX = hipcc
 # example: gfx600 (tahiti), gfx700 (kaveri), gfx701 (hawaii), gfx801 (carrizo),
 # gfx900, gfx1010, etc.
 HIP_ARCH = gfx900
-HIP_FLAGS = --offload-arch=$(HIP_ARCH)
+HIP_FLAGS = --amdgpu-target=$(HIP_ARCH)
 HIP_XCOMPILER =
 HIP_XLINKER   = -Wl,

@@ -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)
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// 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.
 //
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// 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.
 //
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -35,8 +35,8 @@ cd $(git rev-parse --show-toplevel)
 # copyright check
 copyright=true
 if [[ "${option}" == "--copyright" || "${option}" == "" ]]; then
-  if git grep -n "^\(#\|//\).*Copyright.*2010-20\(2[^3]\|[^2].\)" > matches.txt; then
-    echo "Please update the following files to Copyright (c) 2010-2023:"
+  if git grep -l "^\(#\|//\).*\(\-2021\|\ 2010,\)" > matches.txt; then
+    echo "Please update the following files to Copyright (c) 2010-2022:"
    cat matches.txt
    copyright=false
  fi
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -58,6 +58,10 @@ groups_serial=(
   "HiOp examples:"
   "examples/hiop"
   "ex9.cpp"'
+'"ipopt"
+   "IpOpt examples:"
+   "examples/ipopt"
+   "ex10.cpp"'
 '"pumi"
   "PUMI examples:"
   "examples/pumi"
@@ -138,7 +142,7 @@ groups_parallel=(
 '"meshing"
   "Meshing miniapps:"
   "miniapps/meshing"
-   "pmesh-optimizer.cpp pmesh-fitting.cpp pminimal-surface.cpp"'
+   "pmesh-optimizer.cpp pminimal-surface.cpp"'
 '"electromagnetics"
   "Electromagnetics miniapps:"
   "miniapps/electromagnetics"
@@ -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"
@@ -227,7 +235,7 @@ groups_all=(
   "Meshing miniapps:"
   "miniapps/meshing"
   "mobius-strip.cpp klein-bottle.cpp extruder.cpp toroid.cpp
-    {,p}mesh-optimizer.cpp pmesh-fitting.cpp {,p}minimal-surface.cpp"'
+    {,p}mesh-optimizer.cpp {,p}minimal-surface.cpp"'
 '"electromagnetics"
   "Electromagnetics miniapps:"
   "miniapps/electromagnetics"
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+// 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.
 //
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,109 +0,0 @@
-MFEM NURBS mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# SEGMENT     = 1
-# SQUARE      = 3
-# CUBE        = 5
-#
-
-dimension
-2
-
-elements
-5
-1 3 4 5 6 7
-1 3 0 1 5 4
-1 3 1 2 6 5
-1 3 3 7 6 2
-1 3 0 4 7 3
-
-boundary
-4
-1 1 0 1
-1 1 2 3
-1 1 1 2
-1 1 3 0
-
-edges
-12
-0 0 1
-0 4 5
-0 7 6
-0 3 2
-1 1 2
-1 5 6
-1 4 7
-1 0 3
-2 0 4
-2 1 5
-2 2 6
-2 3 7
-
-vertices
-8
-
-knotvectors
-3
-2 3 0 0 0 1 1 1
-2 3 0 0 0 1 1 1
-2 3 0 0 0 1 1 1
-
-weights
-1
-1
-1
-1
-1
-1
-1
-1
-0.70710678118655
-1
-1
-0.70710678118655
-0.70710678118655
-1
-1
-0.70710678118655
-1
-1
-1
-1
-1
-0.85355339059327
-0.85355339059327
-0.85355339059327
-0.85355339059327
-
-FiniteElementSpace
-FiniteElementCollection: NURBS2
-VDim: 2
-Ordering: 1
-
-0.70710678118 -0.70710678118
-0.70710678118 -0.70710678118
-0.70710678118 0.70710678118
-0.70710678118 0.70710678118
-0.35355339059 -0.35355339059
-0.35355339059 -0.35355339059
-0.35355339059 0.35355339059
-0.35355339059 0.35355339059
-0 -1.41421356236
-0 -0.35355339059
-0 0.35355339059
-0 1.41421356236
-1.41421356236 0
-0.35355339059 0
-0.35355339059 0
-1.41421356236 0
-0.530330085885 -0.530330085885
-0.530330085885 -0.530330085885
-0.530330085885 0.530330085885
-0.530330085885 0.530330085885
-0 0
-0 -0.883883476475
-0.883883476475 0
-0 0.883883476475
-0.883883476475 0
@@ -1,46 +0,0 @@
-MFEM NURBS mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# SEGMENT     = 1
-# SQUARE      = 3
-# CUBE        = 5
-#
-
-dimension
-1
-
-elements
-1
-1 1 0 1
-
-boundary
-2
-1 0 0
-2 0 1
-
-edges
-1
-0 0 1
-
-vertices
-2
-
-knotvectors
-1
-1 2 0 0 1 1
-
-weights
-1
-1
-
-FiniteElementSpace
-FiniteElementCollection: NURBS1
-VDim: 1
-Ordering: 1
-
-0
-1
-
-
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -38,7 +38,7 @@ PROJECT_NAME           = "MFEM"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.

-PROJECT_NUMBER         = v4.5.3
+PROJECT_NUMBER         = v4.5.1

 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
@@ -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 \
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -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)
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# 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.
 #
@@ -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;
+}
@@ -1,141 +0,0 @@
-#include "mfem.hpp"
-#include "Problems.hpp"
-#include "IPsolver.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-
-double fRhs(const Vector &pt);
-double obstacle(const Vector &pt);
-double dmanufacturedFun(const Vector &pt);
-
-int main(int argc, char *argv[])
-{
-  int FEorder = 1; // order of the finite elements
-  int linSolver = 0;
-  int maxIPMiters = 30;
-  bool iAmRoot = true;
-  int ref_levels = 3; 
-  OptionsParser args(argc, argv);
-  args.AddOption(&FEorder, "-o", "--order",\
-		  "Order of the finite elements.");
-  args.AddOption(&linSolver, "-linSolver", "--linearSolver", \
-       "IP-Newton linear system solution strategy.");
-  args.AddOption(&maxIPMiters, "-IPMiters", "--IPMiters",\
-		  "Maximum number of IPM iterations");
-  args.AddOption(&ref_levels, "-r", "--mesh_refinement", \
-		  "Mesh Refinement");
-  
-  args.Parse();
-  if(!args.Good())
-  {
-    args.PrintUsage(cout);
-    return 1;
-  }
-  else
-  {
-    if( iAmRoot )
-    {
-      args.PrintOptions(cout);
-    }
-  }
-
-  const char *meshFile = "../../data/inline-quad.mesh";
-  Mesh *mesh = new Mesh(meshFile, 1, 1);
-  int dim = mesh->Dimension(); // geometric dimension of the domain
-  for (int l = 0; l < ref_levels; l++)
-  {
-      mesh->UniformRefinement();
-  }
-
-  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
-  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
-  Array<int> ess_tdof_list;
-  if (mesh->bdr_attributes.Size())
-  {
-     Array<int> ess_bdr(mesh->bdr_attributes.Max());
-     ess_bdr = 1;
-     Vh->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-  }
-  
-  double DC_val = 0.0;
-  int dimD = Vh->GetTrueVSize();
-  Vector x0(dimD); x0 = DC_val;
-  Vector xf(dimD); xf = 0.0;
-
-  ObstacleProblem problem(Vh, x0, &fRhs, &obstacle, ess_tdof_list);
-  
-
-  InteriorPointSolver optimizer(&problem); 
-  optimizer.SetTol(1.e-7);
-  optimizer.SetLinearSolver(linSolver);
-  optimizer.SetMaxIter(maxIPMiters);
-  optimizer.Mult(x0, xf);
-  
-  double Einitial = problem.E(x0);
-  double Efinal = problem.E(xf);
-  cout << "Energy objective at initial point = " << Einitial << endl;
-  cout << "Energy objective at optimizer = " << Efinal << endl;
-
-  GridFunction d_gf(Vh);
-
-  d_gf = xf;
-  
-  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
-  GridFunction dm_gf(Vh);
-  dm_gf.ProjectCoefficient(dm_fc);
-
-  
-  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
-  paraview_dc.SetPrefixPath("ParaView");
-  paraview_dc.SetLevelsOfDetail(FEorder);
-  paraview_dc.SetDataFormat(VTKFormat::BINARY);
-  paraview_dc.SetHighOrderOutput(true);
-  paraview_dc.SetCycle(0);
-  paraview_dc.SetTime(0.0);
-  paraview_dc.RegisterField("d(x) (numerical)", &d_gf);
-  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
-  paraview_dc.Save();
-  
-
-  delete Vh;
-  delete fec;
-  delete mesh;
-  return 0;
-}
-
-
-
-double dmanufacturedFun(const Vector &pt)
-{
-  double alpha = 16.5;
-  return sin(M_PI * pt(1)) * (sin(M_PI * pt(0)) - alpha * pow(pt(0) * (1. - pt(0)), 2));
-}
-
-
-// f(x) forcing term... which enters the objective energy functional
-// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
-// of f(x). f(x) is such that in the absence of bound-constraints then
-// the solution of the optimization problem satisfies the PDE
-// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
-// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
-
-double fRhs(const Vector &pt)
-{
-  double alpha = 16.5;
-  double fx;
-  fx = pow(M_PI, 2) * sin(M_PI * pt(0));
-  fx += alpha * (2. * pow(pt(0), 2) + 2. * pow(1.-pt(0), 2) - 8. * pt(0) * (1.-pt(0)));
-  fx += pow(M_PI, 2) * sin(M_PI * pt(0)) * dmanufacturedFun(pt);
-  fx *= sin(M_PI * pt(1));
-  return fx;
-}
-
-double obstacle(const Vector &pt)
-{
-  return 0.0;
-}
@@ -1,156 +0,0 @@
-#include "mfem.hpp"
-#include "Problems.hpp"
-#include "IPsolver.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-
-double fRhs(const Vector &pt);
-double obstacle(const Vector &pt);
-double dmanufacturedFun(const Vector &pt);
-
-int main(int argc, char *argv[])
-{
-  int FEorder = 1; // order of the finite elements
-  int linSolver = 0;
-  int maxIPMiters = 30;
-  bool iAmRoot = true;
-  int ref_levels = 3; 
-  OptionsParser args(argc, argv);
-  args.AddOption(&FEorder, "-o", "--order",\
-		  "Order of the finite elements.");
-  args.AddOption(&linSolver, "-linSolver", "--linearSolver", \
-       "IP-Newton linear system solution strategy.");
-  args.AddOption(&maxIPMiters, "-IPMiters", "--IPMiters",\
-		  "Maximum number of IPM iterations");
-  args.AddOption(&ref_levels, "-r", "--mesh_refinement", \
-		  "Mesh Refinement");
-  
-  args.Parse();
-  if(!args.Good())
-  {
-    args.PrintUsage(cout);
-    return 1;
-  }
-  else
-  {
-    if( iAmRoot )
-    {
-      args.PrintOptions(cout);
-    }
-  }
-
-  const char *meshFile = "../../data/inline-quad.mesh";
-  Mesh *mesh = new Mesh(meshFile, 1, 1);
-  int dim = mesh->Dimension(); // geometric dimension of the domain
-  for (int l = 0; l < ref_levels; l++)
-  {
-      mesh->UniformRefinement();
-  }
-
-  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
-  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
-  Array<int> ess_tdof_list;
-  if (mesh->bdr_attributes.Size())
-  {
-     Array<int> ess_bdr(mesh->bdr_attributes.Max());
-     ess_bdr = 1;
-     Vh->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-  }
-  
-  double DC_val = 0.06;
-  Vector x0DC(Vh->GetTrueVSize()); x0DC = DC_val;
-  int dimD = Vh->GetTrueVSize() - ess_tdof_list.Size();
-  Vector x0(dimD); x0 = 0.0;
-  Vector xf(dimD); xf = 0.0;
-
-  ObstacleProblemVariant problem(Vh, x0DC, &fRhs, &obstacle, ess_tdof_list);
-
-  InteriorPointSolver optimizer(&problem); 
-  optimizer.SetTol(1.e-7);
-  optimizer.SetLinearSolver(linSolver);
-  optimizer.SetMaxIter(maxIPMiters);
-  optimizer.Mult(x0, xf);
-  
-  double Einitial = problem.E(x0);
-  double Efinal = problem.E(xf);
-  cout << "Energy objective at initial point = " << Einitial << endl;
-  cout << "Energy objective at optimizer = " << Efinal << endl;
-
-  Array<int> noness_tdof_list;
-  noness_tdof_list.SetSize(dimD);
-  int i = 0;
-  for(int j = 0; j < Vh->GetTrueVSize(); j++)
-  {
-     if(ess_tdof_list.Find(j) == -1)
-     {
-        noness_tdof_list[i] = j;
-        i += 1;
-     }
-  }
-
-
-  GridFunction d_gf(Vh);
-  d_gf.Set(1.0, x0DC);
-  d_gf.SetSubVector(noness_tdof_list, xf);
-
-
-
-  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
-  GridFunction dm_gf(Vh);
-  dm_gf.ProjectCoefficient(dm_fc);
-
-  
-  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
-  paraview_dc.SetPrefixPath("ParaView");
-  paraview_dc.SetLevelsOfDetail(FEorder);
-  paraview_dc.SetDataFormat(VTKFormat::BINARY);
-  paraview_dc.SetHighOrderOutput(true);
-  paraview_dc.SetCycle(0);
-  paraview_dc.SetTime(0.0);
-  paraview_dc.RegisterField("d(x) (numerical)", &d_gf);
-  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
-  paraview_dc.Save();
- 
-
-  delete Vh;
-  delete fec;
-  delete mesh;
-  return 0;
-}
-
-
-
-double dmanufacturedFun(const Vector &pt)
-{
-  double alpha = 16.5;
-  return sin(M_PI * pt(1)) * (sin(M_PI * pt(0)) - alpha * pow(pt(0) * (1. - pt(0)), 2));
-}
-
-
-// f(x) forcing term... which enters the objective energy functional
-// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
-// of f(x). f(x) is such that in the absence of bound-constraints then
-// the solution of the optimization problem satisfies the PDE
-// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
-// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
-
-double fRhs(const Vector &pt)
-{
-  double alpha = 16.5;
-  double fx;
-  fx = pow(M_PI, 2) * sin(M_PI * pt(0));
-  fx += alpha * (2. * pow(pt(0), 2) + 2. * pow(1.-pt(0), 2) - 8. * pt(0) * (1.-pt(0)));
-  fx += pow(M_PI, 2) * sin(M_PI * pt(0)) * dmanufacturedFun(pt);
-  fx *= sin(M_PI * pt(1));
-  return fx;
-}
-
-double obstacle(const Vector &pt)
-{
-  return 0.0;
-}
@@ -1,827 +0,0 @@
-#include "mfem.hpp"
-#include "IPsolver.hpp"
-#include "Problems.hpp"
-#include <fstream>
-#include <iostream>
-#include <cstdlib>
-
-using namespace std;
-using namespace mfem;
-
-
-
-InteriorPointSolver::InteriorPointSolver(GeneralOptProblem * Problem) : optProblem(Problem), block_offsetsumlz(5), block_offsetsuml(4), block_offsetsx(3),
-Huu(nullptr), Hum(nullptr), Hmu(nullptr), Hmm(nullptr), Wmm(nullptr), D(nullptr), Ju(nullptr), Jm(nullptr), JuT(nullptr), JmT(nullptr), Huucl(nullptr), HLuu(nullptr), saveLogBarrierIterates(false)
-{
-  rel_tol  = 1.e-2;
-  max_iter = 20;
-  mu_k     = 1.0;
-
-  sMax     = 1.e2;
-  kSig     = 1.e10;   // control deviation from primal Hessian
-  tauMin   = 0.8;     // control rate at which iterates can approach the boundary
-  eta      = 1.e-4;   // backtracking constant
-  thetaMin = 1.e-4;   // allowed violation of the equality constraints
-
-  // constants in line-step A-5.4
-  delta    = 1.0;
-  sTheta   = 1.1;
-  sPhi     = 2.3;
-
-  // control the rate at which the penalty parameter is decreased
-  kMu     = 0.2;
-  thetaMu = 1.5;
-
-  // TO DO -- include the filter
-
-  thetaMax = 1.e6; // maximum constraint violation
-  // data for the second order correction
-  kSoc     = 0.99;
-
-  // equation (18)
-  gTheta = 1.e-5;
-  gPhi   = 1.e-5;
-
-  kEps   = 1.e1;
-
-  dimU = optProblem->GetDimU();
-  dimM = optProblem->GetDimM();
-  dimC = optProblem->GetDimC();
-  ckSoc.SetSize(dimC);
-  
-  block_offsetsumlz[0] = 0;
-  block_offsetsumlz[1] = dimU; // u
-  block_offsetsumlz[2] = dimM; // m
-  block_offsetsumlz[3] = dimC; // lambda
-  block_offsetsumlz[4] = dimM; // zl
-  block_offsetsumlz.PartialSum();
-
-  for(int i = 0; i < block_offsetsuml.Size(); i++)  { block_offsetsuml[i] = block_offsetsumlz[i]; }
-  for(int i = 0; i < block_offsetsx.Size(); i++)    { block_offsetsx[i]   = block_offsetsuml[i] ; }
-
-  // lower-bound for the inequality constraint m >= ml
-  ml = optProblem->Getml();
-  
-  lk.SetSize(dimC);  lk  = 0.0;
-  zlk.SetSize(dimM); zlk = 0.0;
-
-  linSolver = 0;
-  MyRank = 0;
-  iAmRoot = MyRank == 0 ? true : false;
-}
-
-double InteriorPointSolver::MaxStepSize(Vector &x, Vector &xl, Vector &xhat, double tau)
-{
-  double alphaMaxloc = 1.0;
-  double alphaTmp;
-  for(int i = 0; i < x.Size(); i++)
-  {
-    if( xhat(i) < 0. )
-    {
-      alphaTmp = -1. * tau * (x(i) - xl(i)) / xhat(i);
-      alphaMaxloc = min(alphaMaxloc, alphaTmp);
-    } 
-  }
-
-  // alphaMaxloc is the local maximum step size which is
-  // distinct on each MPI process. Need to compute
-  // the global maximum step size 
-  double alphaMaxglb;
-  alphaMaxglb = alphaMaxloc;
-  return alphaMaxglb;
-}
-
-double InteriorPointSolver::MaxStepSize(Vector &x, Vector &xhat, double tau)
-{
-  Vector zero(x.Size()); zero = 0.0;
-  return MaxStepSize(x, zero, xhat, tau);
-}
-
-
-void InteriorPointSolver::Mult(const Vector &x0, Vector &xf)
-{
-  BlockVector x0block(block_offsetsx); x0block = 0.0;
-  x0block.GetBlock(0).Set(1.0, x0);
-  // To do: give options for user specificiation of initialization m0
-  x0block.GetBlock(1) = 100.;
-  x0block.GetBlock(1).Add(1.0, ml);
-  BlockVector xfblock(block_offsetsx); xfblock = 0.0;
-  Mult(x0block, xfblock);
-  xf.Set(1.0, xfblock.GetBlock(0));
-}
-
-void InteriorPointSolver::Mult(const BlockVector &x0, BlockVector &xf)
-{
-  converged = false;
-  
-  BlockVector xk(block_offsetsx), xhat(block_offsetsx); xk = 0; xhat = 0.0;
-  BlockVector Xk(block_offsetsumlz), Xhat(block_offsetsumlz); Xk = 0.0; Xhat = 0.0;
-  BlockVector Xhatuml(block_offsetsuml); Xhatuml = 0.0;
-  Vector zlhat(dimM); zlhat = 0.0;
-
-  xk.GetBlock(0).Set(1.0, x0.GetBlock(0));
-  xk.GetBlock(1).Set(1.0, x0.GetBlock(1));
-  // running estimate of the final values of the Lagrange multipliers
-  lk  = 0.0;
-  zlk = 0.0;
-
-  for(int i = 0; i < dimM; i++)
-  {
-    zlk(i) = 1.e1 * mu_k / (xk(i+dimU) - ml(i));
-  }
-
-  Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
-  Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
-  Xk.GetBlock(2).Set(1.0, lk);
-  Xk.GetBlock(3).Set(1.0, zlk);
-
-  /* set theta0 = theta(x0)
-   *     thetaMin
-   *     thetaMax
-   * when theta(xk) < thetaMin and the switching condition holds
-   * then we ask for the Armijo sufficient decrease of the barrier
-   * objective to be satisfied, in order to accept the trial step length alphakl
-   * 
-   * thetaMax controls how the filter is initialized for each log-barrier subproblem
-   * F0 = {(th, phi) s.t. th > thetaMax}
-   * that is the filter does not allow for iterates where the constraint violation
-   * is larger than that of thetaMax
-   */
-  double theta0 = theta(xk);
-  thetaMin = 1.e-4 * max(1.0, theta0);
-  thetaMax = 1.e8  * thetaMin;
-
-  double Eeval, maxBarrierSolves, Eevalmu0;
-  bool printOptimalityError; // control optimality error print to console for log-barrier subproblems
-  
-  maxBarrierSolves = 10;
-
-  for(jOpt = 0; jOpt < max_iter; jOpt++)
-  {
-    mfem::out << "interior-point solve step " << jOpt << endl;
-    // A-2. Check convergence of overall optimization problem
-    printOptimalityError = false;
-    Eevalmu0 = E(xk, lk, zlk, printOptimalityError);
-    if(Eevalmu0 < rel_tol)
-    {
-      converged = true;
-      mfem::out << "solved optimization problem :)\n";
-      break;
-    }
-    
-    if(jOpt > 0) { maxBarrierSolves = 1; }
-    
-    for(int i = 0; i < maxBarrierSolves; i++)
-    {
-      // A-3. Check convergence of the barrier subproblem
-      printOptimalityError = true;
-      Eeval = E(xk, lk, zlk, mu_k, printOptimalityError);
-      if(Eeval < kEps * mu_k)
-      {
-        mfem::out << "solved barrier subproblem, for mu = " << mu_k << endl;
-        // A-3.1. Recompute the barrier parameter
-        mu_k  = max(rel_tol / 10., min(kMu * mu_k, pow(mu_k, thetaMu)));
-        // A-3.2. Re-initialize the filter
-        F1.DeleteAll();
-        F2.DeleteAll();
-      }
-      else
-      {
-        break;
-      }
-    }
-    
-    // A-4. Compute the search direction
-    // solve for (uhat, mhat, lhat)
-    mfem::out << "\n** A-4. IP-Newton solve **\n";
-    zlhat = 0.0; Xhatuml = 0.0;
-    // why do we have Xhatuml ....???
-    // TO DO: remove Xhatuml in favor of passing Xhat
-    IPNewtonSolve(xk, lk, zlk, zlhat, Xhatuml, mu_k, false); 
-
-
-    // assign data stack, X = (u, m, l, zl)
-    Xk = 0.0;
-    Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
-    Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
-    Xk.GetBlock(2).Set(1.0, lk);
-    Xk.GetBlock(3).Set(1.0, zlk);
-
-    // assign data stack, Xhat = (uhat, mhat, lhat, zlhat)
-    Xhat = 0.0;
-    for(int i = 0; i < 3; i++)
-    {
-      Xhat.GetBlock(i).Set(1.0, Xhatuml.GetBlock(i));
-    }
-    Xhat.GetBlock(3).Set(1.0, zlhat);
-
-
-    // A-5. Backtracking line search.
-    mfem::out << "\n** A-5. Linesearch **\n";
-    mfem::out << "mu = " << mu_k << endl;
-    
-    lineSearch(Xk, Xhat, mu_k);
-    
-    if(lineSearchSuccess)
-    {
-      if(!switchCondition || !sufficientDecrease)
-      {
-        F1.Append( (1. - gTheta) * thx0);
-        F2.Append( phx0 - gPhi * thx0);
-      }
-      // ----- A-6: Accept the trial point
-      // print info regarding zl...
-      xk.GetBlock(0).Add(alpha, Xhat.GetBlock(0));
-      xk.GetBlock(1).Add(alpha, Xhat.GetBlock(1));
-      lk.Add(alpha,   Xhat.GetBlock(2));
-      zlk.Add(alphaz, Xhat.GetBlock(3));
-      projectZ(xk, zlk, mu_k);
-    }
-    else
-    {
-      mfem::out << "lineSearch not successful :(\n";
-      mfem::out << "attempting feasibility restoration with theta = " << thx0 << endl;
-      mfem::out << "no feasibility restoration implemented, exiting now \n";
-      break;
-    }
-    //
-    if(jOpt + 1 == max_iter) 
-    { 
-      mfem::out << "maximum optimization iterations :(\n";
-    }
-  }
-  // done with optimization routine, just reassign data to xf reference so
-  // that the application code has access to the optimal point
-  xf = 0.0;
-  xf.GetBlock(0).Set(1.0, xk.GetBlock(0));
-  xf.GetBlock(1).Set(1.0, xk.GetBlock(1));
-}
-
-void InteriorPointSolver::FormIPNewtonMat(BlockVector & x, Vector & l, Vector &zl, BlockOperator &Ak)
-{
-  // WARNING: Huu, Hum, Hmu, Hmm should all be Hessian terms of the Lagrangian, currently we 
-  //          them by Hessian terms of the objective function and neglect the Hessian of l^T c
-
-  Huu = optProblem->Duuf(x); Hum = optProblem->Dumf(x);
-  Hmu = optProblem->Dmuf(x); Hmm = optProblem->Dmmf(x);
-
-  Vector DiagLogBar(dimM); DiagLogBar = 0.0;
-  for(int ii = 0; ii < dimM; ii++)
-  {
-    DiagLogBar(ii) = zl(ii) / (x(ii+dimU) - ml(ii));
-  }
-  if(saveLogBarrierIterates)
-  {
-    std::ofstream diagStream;
-    char diagString[100];
-    snprintf(diagString, 100, "logBarrierHessiandata/D%d.dat", jOpt);
-    diagStream.open(diagString, ios::out | ios::trunc);
-    for(int ii = 0; ii < dimM; ii++)
-    {
-      diagStream << setprecision(30) << DiagLogBar(ii) << endl;
-    }
-    diagStream.close();
-  } 
-
-
-  D = new SparseMatrix(DiagLogBar);
-  
-  if(Hmm != nullptr)
-  {
-    Wmm = new SparseMatrix(*Hmm);
-    Wmm->Add(1.0, *D);
-  }
-  else
-  {
-    Wmm = D;
-  }
-
-  Ju = optProblem->Duc(x); JuT = Transpose(*Ju);
-  Jm = optProblem->Dmc(x); JmT = Transpose(*Jm);
- 
-  Huucl = optProblem->lDuuc(x, l);
-  if(Huucl != nullptr)
-  {
-    HLuu = Add(*Huucl, *Huu);
-    Ak.SetBlock(0, 0, HLuu);
-  } 
-  else
-  {
-    Ak.SetBlock(0, 0, Huu);
-  }
-
-  //         IP-Newton system matrix
-  //    Ak = [[H_(u,u)  H_(u,m)   J_u^T]
-  //          [H_(m,u)  W_(m,m)   J_m^T]
-  //          [ J_u      J_m       0  ]]
-
-  Ak.SetBlock(0, 2, JuT);
-  Ak.SetBlock(1, 1, Wmm); Ak.SetBlock(1, 2, JmT);
-  Ak.SetBlock(2, 0,  Ju); Ak.SetBlock(2, 1,  Jm);
-
-  if(Hum != nullptr) { Ak.SetBlock(0, 1, Hum); Ak.SetBlock(1, 0, Hmu); }
-}
-
-
-// perturbed KKT system solve
-// determine the search direction
-void InteriorPointSolver::IPNewtonSolve(BlockVector &x, Vector &l, Vector &zl, Vector &zlhat, BlockVector &Xhat, double mu, bool socSolve)
-{
-  // solve A x = b, where A is the IP-Newton matrix
-  BlockOperator A(block_offsetsuml, block_offsetsuml); BlockVector b(block_offsetsuml); b = 0.0;
-  FormIPNewtonMat(x, l, zl, A);
-
-  //       [grad_u phi + Ju^T l]
-  // b = - [grad_m phi + Jm^T l]
-  //       [          c        ]
-  BlockVector gradphi(block_offsetsx); gradphi = 0.0;
-  BlockVector JTl(block_offsetsx); JTl = 0.0;
-  Dxphi(x, mu, gradphi);
-  
-  (A.GetBlock(0,2)).Mult(l, JTl.GetBlock(0));
-  (A.GetBlock(1,2)).Mult(l, JTl.GetBlock(1));
-
-  for(int ii = 0; ii < 2; ii++)
-  {
-    b.GetBlock(ii).Set(1.0, gradphi.GetBlock(ii));
-    b.GetBlock(ii).Add(1.0, JTl.GetBlock(ii));
-  }
-  if(!socSolve) 
-  {
-    optProblem->c(x, b.GetBlock(2));
-  }
-  else
-  {
-    b.GetBlock(2).Set(1.0, ckSoc);
-  }
-  b *= -1.0; 
-  Xhat = 0.0;
-
-
-  #ifdef MFEM_USE_SUITESPARSE
-    // Direct solve for IP-Newton saddle-point system
-    // A = [ [ Huu   0    Ju^T]
-    //       [  0    D     -I ]
-    //       [ Ju   -I      0 ]]
-    if(linSolver == 0)
-    {
-      BlockMatrix ABlockMatrix(block_offsetsuml, block_offsetsuml);
-      for(int ii = 0; ii < 3; ii++)
-      {
-        for(int jj = 0; jj < 3; jj++)
-        {
-          if(!A.IsZeroBlock(ii, jj))
-          {
-            ABlockMatrix.SetBlock(ii, jj, dynamic_cast<SparseMatrix *>(&(A.GetBlock(ii, jj))));
-          }
-        }
-      }
-      /* direct solve of the 3x3 IP-Newton linear system */
-      UMFPackSolver ASolver;
-      SparseMatrix *ASparse = ABlockMatrix.CreateMonolithic();
-      ASolver.SetOperator(*ASparse);
-      ASolver.Mult(b, Xhat);
-
-      Vector residual(Xhat.Size());
-      ASparse->Mult(Xhat, residual);
-      residual.Add(-1.0, b);
-      delete ASparse;
-    }
-    else if(linSolver == 1)
-    {
-      // Direct solve for 0,0 Schur complement of IP-Newton system, Huu + Ju^T Wmm Ju,
-      // where Wmm = D for contact problems
-      SparseMatrix * Huuloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 0)));
-      SparseMatrix * Wmmloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(1, 1)));
-      SparseMatrix * Juloc  = dynamic_cast<SparseMatrix *>(&(A.GetBlock(2, 0)));
-      SparseMatrix * JuTloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 2)));
-      Vector DVec(dimM); DVec = 0.0;
-      Vector one(dimM); one = 1.0;
-      D->Mult(one, DVec);
-      SparseMatrix *JuTDJu = Mult_AtDA(*Juloc, DVec);     // Ju^T D Ju
-      SparseMatrix *Areduced = Add(*Huuloc, *JuTDJu);  // Huu + Ju^T D Ju
-
-      
-      /* prepare the reduced rhs */
-      // breduced = bu + Ju^T (bm + Wmm bl)
-      Vector breduced(dimU); breduced = 0.0;
-      Vector tempVec(dimM); tempVec = 0.0;
-      Wmmloc->Mult(b.GetBlock(2), tempVec);
-      tempVec.Add(1.0, b.GetBlock(1));
-      JuTloc->Mult(tempVec, breduced);
-      breduced.Add(1.0, b.GetBlock(0));
-
-      // solve the reduced linear system
-      UMFPackSolver AreducedSolver;
-      AreducedSolver.SetOperator(*Areduced);
-      AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-
-      // now propagate solved uhat to obtain mhat and lhat
-      // xm = Ju xu - bl
-      Juloc->Mult(Xhat.GetBlock(0), Xhat.GetBlock(1));
-      Xhat.GetBlock(1).Add(-1.0, b.GetBlock(2));
-
-      // xl = Wmm xm - bm
-      Wmmloc->Mult(Xhat.GetBlock(1), Xhat.GetBlock(2));
-      Xhat.GetBlock(2).Add(-1.0, b.GetBlock(1));
-
-      delete JuTDJu;
-      delete Areduced;
-    }
-  #else
-    MFEM_VERIFY(linSolver > 1, "linSolver = 0, 1 require MFEM_USE_SUITESPARSE=YES");
-  #endif
-  if(linSolver > 1)
-  { 
-    // Iterative solve for 0,0 Schur complement of IP-Newton system, Huu + Ju^T Wmm Ju,
-    // where Wmm = D for contact problems
-    // here the iterative solver is a Jacobi-preconditioned CG-solve
-    SparseMatrix * Huuloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 0)));
-    SparseMatrix * Wmmloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(1, 1)));
-    SparseMatrix * Juloc  = dynamic_cast<SparseMatrix *>(&(A.GetBlock(2, 0)));
-    SparseMatrix * JuTloc = dynamic_cast<SparseMatrix *>(&(A.GetBlock(0, 2)));
-    Vector DVec(dimM); DVec = 0.0;
-    Vector one(dimM); one = 1.0;
-    D->Mult(one, DVec);
-    SparseMatrix *JuTDJu   = Mult_AtDA(*Juloc, DVec);     // Ju^T D Ju
-    SparseMatrix *Areduced = Add(*Huuloc, *JuTDJu);  // Huu + Ju^T D Ju
-
-    /* prepare the reduced rhs */
-    // breduced = bu + Ju^T (bm + Wmm bl)
-    Vector breduced(dimU); breduced = 0.0;
-    Vector tempVec(dimM); tempVec = 0.0;
-    Wmmloc->Mult(b.GetBlock(2), tempVec);
-    tempVec.Add(1.0, b.GetBlock(1));
-    JuTloc->Mult(tempVec, breduced);
-    breduced.Add(1.0, b.GetBlock(0));
-    
-    if (linSolver == 2)
-    {
-      /* Jacobi preconditioned conjugate-gradient solve */
-      DSmoother AreducedPrec((SparseMatrix &)(*Areduced));
-      CGSolver AreducedSolver;
-      AreducedSolver.SetOperator(*Areduced);
-      AreducedSolver.SetAbsTol(1.e-12);
-      AreducedSolver.SetRelTol(1.e-8);
-      AreducedSolver.SetMaxIter(500);
-      AreducedSolver.SetPreconditioner(AreducedPrec);
-      AreducedSolver.SetPrintLevel(1);
-      AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-    }
-    else
-    {
-      /* Gauss-Seidel preconditioned GMRES solve */
-      GSSmoother AreducedPrec((SparseMatrix &)(*Areduced));
-      GMRESSolver AreducedSolver;
-      AreducedSolver.SetOperator(*Areduced);
-      AreducedSolver.SetAbsTol(1.e-12);
-      AreducedSolver.SetRelTol(1.e-8);
-      AreducedSolver.SetMaxIter(500);
-      AreducedSolver.SetPreconditioner(AreducedPrec);
-      AreducedSolver.SetPrintLevel(1);
-      AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-    }
-    
-    // now propagate solved uhat to obtain mhat and lhat
-    // xm = Ju xu - bl
-    Juloc->Mult(Xhat.GetBlock(0), Xhat.GetBlock(1));
-    Xhat.GetBlock(1).Add(-1.0, b.GetBlock(2));
-
-    // xl = Wmm xm - bm
-    Wmmloc->Mult(Xhat.GetBlock(1), Xhat.GetBlock(2));
-    Xhat.GetBlock(2).Add(-1.0, b.GetBlock(1));
-
-    delete JuTDJu;
-    delete Areduced;
-  }
-    
-  
-  /* backsolve to determine zlhat */
-  for(int ii = 0; ii < dimM; ii++)
-  {
-    zlhat(ii) = -1.*(zl(ii) + (zl(ii) * Xhat(ii + dimU) - mu) / (x(ii + dimU) - ml(ii)) );
-  }
-
-  // free memory
-  if(Hmm != nullptr)
-  {
-    delete Wmm;
-  }
-  if( Huucl != nullptr)
-  {
-    delete HLuu; HLuu = nullptr;
-  }
-  delete D;
-  delete JuT;
-  delete JmT;
-}
-
-// here Xhat, X will be BlockVectors w.r.t. the 4 partitioning X = (u, m, l, zl)
-
-void InteriorPointSolver::lineSearch(BlockVector& X0, BlockVector& Xhat, double mu)
-{
-  double tau  = max(tauMin, 1.0 - mu);
-  Vector u0   = X0.GetBlock(0);
-  Vector m0   = X0.GetBlock(1);
-  Vector l0   = X0.GetBlock(2);
-  Vector z0   = X0.GetBlock(3);
-  Vector uhat = Xhat.GetBlock(0);
-  Vector mhat = Xhat.GetBlock(1);
-  Vector lhat = Xhat.GetBlock(2);
-  Vector zhat = Xhat.GetBlock(3);
-  double alphaMax  = MaxStepSize(m0, ml, mhat, tau);
-  double alphaMaxz = MaxStepSize(z0, zhat, tau);
-  alphaz = alphaMaxz;
-
-
-  BlockVector x0(block_offsetsx); x0 = 0.0;
-  x0.GetBlock(0).Set(1.0, u0);
-  x0.GetBlock(1).Set(1.0, m0);
-  
-  BlockVector xhat(block_offsetsx); xhat = 0.0;
-  xhat.GetBlock(0).Set(1.0, uhat);
-  xhat.GetBlock(1).Set(1.0, mhat);
-  
-  BlockVector xtrial(block_offsetsx); xtrial = 0.0;
-  BlockVector Dxphi0(block_offsetsx); Dxphi0 = 0.0;
-  int maxBacktrack = 20;
-  alpha = alphaMax;
-
-
-  Vector ck0(dimC); ck0 = 0.0;
-  Vector zhatsoc(dimM); zhatsoc = 0.0;
-  BlockVector Xhatumlsoc(block_offsetsuml); Xhatumlsoc = 0.0;
-  BlockVector xhatsoc(block_offsetsx); xhatsoc = 0.0;
-  Vector uhatsoc(dimU); uhatsoc = 0.0;
-  Vector mhatsoc(dimM); mhatsoc = 0.0;
-
-  Dxphi(x0, mu, Dxphi0);
-  Dxphi0_xhat = InnerProduct(Dxphi0, xhat);
-  descentDirection = Dxphi0_xhat < 0. ? true : false;
-  if(descentDirection)
-  {
-    mfem::out << "is a descent direction for the log-barrier objective\n";
-  }
-  else
-  {
-    mfem::out << "is not a descent direction for the log-barrier objective\n";
-  }
-  mfem::out << "Dxphi^T xhat / (|| Dxphi||_2 * || xhat ||_2) = " << Dxphi0_xhat / (xhat.Norml2() * Dxphi0.Norml2()) << endl;
-  thx0 = theta(x0);
-  phx0 = phi(x0, mu);
-
-  lineSearchSuccess = false;
-  for(int i = 0; i < maxBacktrack; i++)
-  {
-    mfem::out << "\n--------- alpha = " << alpha << " ---------\n";
-
-    // ----- A-5.2. Compute trial point: xtrial = x0 + alpha_i xhat
-    xtrial.Set(1.0, x0);
-    xtrial.Add(alpha, xhat);
-
-    // ------ A-5.3. if not in filter region go to A.5.4 otherwise go to A-5.5.
-    thxtrial = theta(xtrial);
-    phxtrial = phi(xtrial, mu);
-
-    filterCheck(thxtrial, phxtrial);    
-    if(!inFilterRegion)
-    {
-      mfem::out << "not in filter region :)\n";
-      // ------ A.5.4: Check sufficient decrease
-      if(!descentDirection)
-      {
-        switchCondition = false;
-      }
-      else
-      {
-        switchCondition = (alpha * pow(abs(Dxphi0_xhat), sPhi) > delta * pow(thx0, sTheta)) ? true : false;
-      }
-      mfem::out << "theta(x0) = "     << thx0     << ", thetaMin = "                  << thetaMin             << endl;
-      mfem::out << "theta(xtrial) = " << thxtrial << ", (1-gTheta) *theta(x0) = "     << (1. - gTheta) * thx0 << endl;
-      mfem::out << "phi(xtrial) = "   << phxtrial << ", phi(x0) - gPhi *theta(x0) = " << phx0 - gPhi * thx0   << endl;
-      
-      // Case I      
-      if(thx0 <= thetaMin && switchCondition)
-      {
-        sufficientDecrease = phxtrial <= phx0 + eta * alpha * Dxphi0_xhat ? true : false;
-        if(sufficientDecrease)
-        {
-          mfem::out << "Accepted step length -- sufficient decrease in log-barrier objective.\n"; 
-          // accept the trial step
-          lineSearchSuccess = true;
-          break;
-        }
-      }
-      else
-      {
-        if(thxtrial <= (1. - gTheta) * thx0 || phxtrial <= phx0 - gPhi * thx0)
-        {
-          mfem::out << "Accepted step length -- decrease in either constraint violation or log-barrier objective.\n"; 
-          // accept the trial step
-          lineSearchSuccess = true;
-          break;
-        }
-      }
-      // A-5.5: Initialize the second-order correction
-      if((!(thx0 < thxtrial)) && i == 0)
-      {
-        mfem::out << "second order correction\n";
-        optProblem->c(xtrial, ckSoc);
-        optProblem->c(x0, ck0);
-        ckSoc.Add(alphaMax, ck0);
-        // A-5.6 Compute the second-order correction.
-        IPNewtonSolve(x0, l0, z0, zhatsoc, Xhatumlsoc, mu, true);
-        mhatsoc.Set(1.0, Xhatumlsoc.GetBlock(1));
-        // alphasoc = MaxStepSize(m0, ml, mhatsoc, tau);
-        //WARNING: not complete but currently solver isn't entering this region
-      }
-    }
-    else
-    {
-      mfem::out << "in filter region\n";
-    }
-
-    // include more if needed
-    alpha *= 0.5;
-
-  } 
-}
-
-
-void InteriorPointSolver::projectZ(const Vector &x, Vector &z, double mu)
-{
-  double zi;
-  double mudivmml;
-  for(int i = 0; i < dimM; i++)
-  {
-    zi = z(i);
-    mudivmml = mu / (x(i + dimU) - ml(i));
-    z(i) = max(min(zi, kSig * mudivmml), mudivmml / kSig);
-  }
-}
-
-void InteriorPointSolver::filterCheck(double th, double ph)
-{
-  inFilterRegion = false;
-  if(th > thetaMax)
-  {
-    inFilterRegion = true;
-  }
-  else
-  {
-    for(int i = 0; i < F1.Size(); i++)
-    {
-      if(th >= F1[i] && ph >= F2[i])
-      {
-        inFilterRegion = true;
-        break;
-      }
-    }
-  }
-}
-
-double InteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, double mu, bool print)
-{
-  double E1, E2, E3;
-  double sc, sd;
-  BlockVector gradL(block_offsetsx); gradL = 0.0; // stationarity grad L = grad f + J^T l - z
-  Vector cx(dimC); cx = 0.0;     // feasibility c = c(x)
-  Vector comp(dimM); comp = 0.0; // complementarity M Z - mu 1
-
-  DxL(x, l, zl, gradL);
-  E1 = gradL.Normlinf();
-
-  optProblem->c(x, cx);
-  E2 = cx.Normlinf();
-
-  for(int ii = 0; ii < dimM; ii++) 
-  { 
-    comp(ii) = x(dimU + ii) * zl(ii) - mu;
-  }
-  E3 = comp.Normlinf();
-
-  double ll1, zl1;
-  zl1 = zl.Norml1() / double(dimC + dimM);
-  ll1 =  l.Norml1();
-  sc = max(sMax, zl1 / (double(dimM)) ) / sMax;
-  sd = max(sMax, (ll1 + zl1) / (double(dimC + dimM))) / sMax;
-  if(print)
-  {
-    mfem::out << "evaluating optimality error for mu = " << mu << endl;
-    mfem::out << "stationarity measure = "    << E1 / sd << endl;
-    mfem::out << "feasibility measure  = "    << E2      << endl;
-    mfem::out << "complimentarity measure = " << E3 / sc << endl;
-  }
-  return max(max(E1 / sd, E2), E3 / sc);
-}
-
-double InteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, bool print)
-{
-  return E(x, l, zl, 0.0, print);
-}
-
-double InteriorPointSolver::theta(const BlockVector &x)
-{
-  Vector cx(dimC); cx = 0.0;
-  optProblem->c(x, cx);
-  return cx.Norml2();
-}
-
-// log-barrier objective
-double InteriorPointSolver::phi(const BlockVector &x, double mu)
-{
-  double fx = optProblem->CalcObjective(x); 
-  double logBarrierLoc = 0.0;
-  for(int i = 0; i < dimM; i++) 
-  { 
-    logBarrierLoc += log(x(dimU+i)-ml(i));
-  }
-  double logBarrierGlb = 0.0;
-  logBarrierGlb = logBarrierLoc;
-  return fx - mu * logBarrierGlb;
-}
-
-
-
-
-// gradient of log-barrier objective with respect to x = (u, m)
-void InteriorPointSolver::Dxphi(const BlockVector &x, double mu, BlockVector &y)
-{
-  optProblem->CalcObjectiveGrad(x, y);
-  for(int i = 0; i < dimM; i++) 
-  { 
-    y(dimU + i) -= mu / (x(dimU + i));
-  } 
-}
-
-
-// Lagrangian function evaluation
-// L(x, l, zl) = f(x) + l^T c(x) - zl^T m
-double InteriorPointSolver::L(const BlockVector &x, const Vector &l, const Vector &zl)
-{
-  double fx = optProblem->CalcObjective(x);
-  Vector cx(dimC); optProblem->c(x, cx);
-  return (fx + InnerProduct(cx, l) - InnerProduct(x.GetBlock(1), zl));
-}
-
-void InteriorPointSolver::DxL(const BlockVector &x, const Vector &l, const Vector &zl, BlockVector &y)
-{
-  // evaluate the gradient of the objective with respect to the primal variables x = (u, m)
-  BlockVector gradxf(block_offsetsx); gradxf = 0.0;
-  optProblem->CalcObjectiveGrad(x, gradxf);
-  
-  SparseMatrix *Jacu, *Jacm, *JacuT, *JacmT;
-  Jacu = optProblem->Duc(x); Jacm = optProblem->Dmc(x);
-  JacuT = Transpose(*Jacu);
-  JacmT = Transpose(*Jacm);
-  JacuT->Mult(l, y.GetBlock(0));
-  JacmT->Mult(l, y.GetBlock(1));
-  delete JacuT; 
-  delete JacmT;
-  y.Add(1.0, gradxf);
-  (y.GetBlock(1)).Add(-1.0, zl);
-}
-
-
-bool InteriorPointSolver::GetConverged() const
-{
-  return converged;
-}
-
-void InteriorPointSolver::SetTol(double Tol)
-{
-  rel_tol = Tol;
-}
-
-void InteriorPointSolver::SetMaxIter(int max_it)
-{
-  max_iter = max_it;
-}
-
-void InteriorPointSolver::SetBarrierParameter(double mu_0)
-{
-  mu_k = mu_0;
-}
-
-void InteriorPointSolver::SaveLogBarrierHessianIterates(bool save)
-{
-  MFEM_ASSERT(MyRank == 0 || save == false, "currently can only save logbarrier hessian in serial codes");
-  saveLogBarrierIterates = save;
-}
-
-void InteriorPointSolver::SetLinearSolver(int LinSolver)
-{
-  linSolver = LinSolver;
-}
-
-
-
-InteriorPointSolver::~InteriorPointSolver() 
-{
-  F1.DeleteAll();
-  F2.DeleteAll();
-  block_offsetsx.DeleteAll();
-  block_offsetsumlz.DeleteAll();
-  block_offsetsuml.DeleteAll();
-  ml.SetSize(0);
-}
@@ -1,80 +0,0 @@
-#include "mfem.hpp"
-#include "Problems.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-#ifndef IPSOLVER 
-#define IPSOLVER
-
-class InteriorPointSolver
-{
-protected:
-    GeneralOptProblem* optProblem;
-    double rel_tol;
-    int  max_iter;
-    double mu_k; // \mu_k
-    Vector lk, zlk;
-
-    double sMax, kSig, tauMin, eta, thetaMin, delta, sTheta, sPhi, kMu, thetaMu;
-    double thetaMax, kSoc, gTheta, gPhi, kEps;
-	
-    // filter
-    Array<double> F1, F2;
-	
-    // quantities computed in lineSearch
-    double alpha, alphaz;
-    double thx0, thxtrial;
-    double phx0, phxtrial;
-    bool descentDirection, switchCondition, sufficientDecrease, lineSearchSuccess, inFilterRegion;
-    double Dxphi0_xhat;
-
-    int dimU, dimM, dimC;
-    Array<int> block_offsetsumlz, block_offsetsuml, block_offsetsx;
-    Vector ml;
-
-    Vector ckSoc;
-    SparseMatrix * Huu, * Hum, * Hmu, * Hmm, * Wmm, *D, * Ju, * Jm, * JuT, * JmT;
-    SparseMatrix * Huucl, *HLuu;    
-    int jOpt;
-    bool converged;
-    
-    int MyRank;
-    bool iAmRoot;
-
-    bool saveLogBarrierIterates;
-
-    int linSolver;
-
-public:
-    InteriorPointSolver(GeneralOptProblem*);
-    void Mult(const BlockVector& , BlockVector&); // used when the user wants to be aware of bound-constrained variable m >= ml
-    void Mult(const Vector&, Vector &); // useful when the user doesn't need to know about bound-constrained variable   m >= ml, e.g., when m is a slack variable
-    double MaxStepSize(Vector& , Vector& , Vector& , double);
-    double MaxStepSize(Vector& , Vector& , double);
-    void FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
-    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , double, bool);
-    void lineSearch(BlockVector& , BlockVector& , double);
-    void projectZ(const Vector & , Vector &, double);
-    void filterCheck(double, double);
-    double E(const BlockVector &, const Vector &, const Vector &, double, bool);
-    double E(const BlockVector &, const Vector &, const Vector &, bool);
-    bool GetConverged() const;
-    // TO DO: include Hessian of Lagrangian
-    double theta(const BlockVector &);
-    double phi(const BlockVector &, double);
-    void Dxphi(const BlockVector &, double, BlockVector &);
-    double L(const BlockVector &, const Vector &, const Vector &);
-    void DxL(const BlockVector &, const Vector &, const Vector &, BlockVector &);
-    void SetTol(double);
-    void SetMaxIter(int);
-    void SetBarrierParameter(double);    
-    void SaveLogBarrierHessianIterates(bool);
-    void SetLinearSolver(int);
-    virtual ~InteriorPointSolver();
-};
-
-#endif
@@ -1,125 +0,0 @@
-#include "mfem.hpp"
-#include "Problems.hpp"
-#include "IPsolver.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-
-double dmanufacturedFun(const Vector &);
-double fRhs(const Vector &);
-double obstacle(const Vector &);
-
-int main(int argc, char *argv[])
-{
-  int FEorder = 1; // order of the finite elements
-  int linSolver = 0;
-  int maxIPMiters = 30;
-  bool iAmRoot = true;
-  
-  OptionsParser args(argc, argv);
-  args.AddOption(&FEorder, "-o", "--order",\
-		  "Order of the finite elements.");
-  args.AddOption(&linSolver, "-linSolver", "--linearSolver", \
-       "IP-Newton linear system solution strategy.");
-  args.AddOption(&maxIPMiters, "-IPMiters", "--IPMiters",\
-		  "Maximum number of IPM iterations");
-  
-  args.Parse();
-  if(!args.Good())
-  {
-    args.PrintUsage(cout);
-    return 1;
-  }
-  else
-  {
-    if( iAmRoot )
-    {
-      args.PrintOptions(cout);
-    }
-  }
-
-  const char *meshFile = "../../data/inline-quad.mesh";
-  Mesh *mesh = new Mesh(meshFile, 1, 1);
-  int dim = mesh->Dimension(); // geometric dimension of the domain
-  {
-     int ref_levels = 3;
-     for (int l = 0; l < ref_levels; l++)
-     {
-        mesh->UniformRefinement();
-     }
-  }
-
-  FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
-  FiniteElementSpace      *Vh  = new FiniteElementSpace(mesh, fec);
-  ObstacleProblem problem(Vh, &fRhs, &obstacle);
-  
-  int dimD = problem.GetDimD();
-  Vector x0(dimD); x0 = 0.0;
-  Vector xf(dimD); xf = 0.0;
-
-  InteriorPointSolver optimizer(&problem); 
-  optimizer.SetTol(1.e-7);
-  optimizer.SetLinearSolver(linSolver);
-  optimizer.SetMaxIter(maxIPMiters);
-  optimizer.Mult(x0, xf);
-
-  double Einitial = problem->E(x0);
-  double Efinal = problem->E(xf);
-  cout << "Energy objective at initial point = " << Einitial << endl;
-  cout << "Energy objective at QP optimizer = " << Efinal << endl;
-
-
-  GridFunction d_gf(Vh);
-
-  d_gf = xf;
-
-  FunctionCoefficient dm_fc(dmanufacturedFun); // pseudo-manufactured solution
-  GridFunction dm_gf(Vh);
-  dm_gf.ProjectCoefficient(dm_fc);
-  
-  ParaViewDataCollection paraview_dc("BarrierProblemSolution", mesh);
-  paraview_dc.SetPrefixPath("ParaView");
-  paraview_dc.SetLevelsOfDetail(FEorder);
-  paraview_dc.SetDataFormat(VTKFormat::BINARY);
-  paraview_dc.SetHighOrderOutput(true);
-  paraview_dc.SetCycle(0);
-  paraview_dc.SetTime(0.0);
-  paraview_dc.RegisterField("d(x) (numerical)", &d_gf);
-  paraview_dc.RegisterField("d(x) (pseudo-manufactured)", &dm_gf);
-  paraview_dc.Save();
-  
-  
-  delete Vh;
-  delete fec;
-  delete mesh;
-  return 0;
-}
-
-
-double dmanufacturedFun(const Vector &x)
-{
-  return cos(2*M_PI*x(0)) + 0.2 - 2.0*(pow(x(0),3) - 1.5*pow(x(0),2));
-}
-
-
-// f(x) forcing term... which enters the objective energy functional
-// E(d) = 0.5 d^T K d - f^T d, where f is a discrete vector representation
-// of f(x). f(x) is such that in the absence of bound-constraints then
-// the solution of the optimization problem satisfies the PDE
-// -div(grad(d)) + d = f + homogeneous Neumann conditions on the unit interval,
-// for d(x) = cos(2 \pi x) + a0 + a3 (x^3 - 1.5 x^2), a2 = 0.2, a3 = -2
-double fRhs(const Vector &x)
-{
-  double fx = 0.;
-  fx = 0.2 - 2.0 * (pow(x(0),3)- 1.5*pow(x(0),2.) - 6 * x(0) + 3.) + (1. + pow(2.*M_PI,2))*cos(2.*M_PI*x(0));
-  return fx;
-}
-
-double obstacle(const Vector &x)
-{
-  return 0.0;
-}
@@ -1,834 +0,0 @@
-#include "mfem.hpp"
-#include "ParIPsolver.hpp"
-#include "ParProblems.hpp"
-#include <fstream>
-#include <iostream>
-#include <cstdlib>
-
-using namespace std;
-using namespace mfem;
-
-
-ParInteriorPointSolver::ParInteriorPointSolver(ParGeneralOptProblem * problem_) 
-                     : problem(problem_), 
-                       block_offsetsumlz(5), block_offsetsuml(4), block_offsetsx(3),
-                       Huu(nullptr), Hum(nullptr), Hmu(nullptr), 
-                       Hmm(nullptr), Wmm(nullptr), D(nullptr), 
-                       Ju(nullptr), Jm(nullptr), JuT(nullptr), JmT(nullptr), 
-                       saveLogBarrierIterates(false)
-{
-   OptTol  = 1.e-2;
-   max_iter = 20;
-   mu_k     = 1.0;
-
-   sMax     = 1.e2;
-   kSig     = 1.e10;   // control deviation from primal Hessian
-   tauMin   = 0.8;     // control rate at which iterates can approach the boundary
-   eta      = 1.e-4;   // backtracking constant
-   thetaMin = 1.e-4;   // allowed violation of the equality constraints
-
-   // constants in line-step A-5.4
-   delta    = 1.0;
-   sTheta   = 1.1;
-   sPhi     = 2.3;
-
-   // control the rate at which the penalty parameter is decreased
-   kMu     = 0.2;
-   thetaMu = 1.5;
-
-   thetaMax = 1.e6; // maximum constraint violation
-   // data for the second order correction
-   kSoc     = 0.99;
-
-   // equation (18)
-   gTheta = 1.e-5;
-   gPhi   = 1.e-5;
-
-   kEps   = 1.e1;
-
-   dimU = problem->GetDimU();
-   dimM = problem->GetDimM();
-   dimC = problem->GetDimC();
-   MPI_Allreduce(&dimU, &dimUglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-   MPI_Allreduce(&dimM, &dimMglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-   MPI_Allreduce(&dimC, &dimCglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-   ckSoc.SetSize(dimC);
-  
-   block_offsetsumlz[0] = 0;
-   block_offsetsumlz[1] = dimU; // u
-   block_offsetsumlz[2] = dimM; // m
-   block_offsetsumlz[3] = dimC; // lambda
-   block_offsetsumlz[4] = dimM; // zl
-   block_offsetsumlz.PartialSum();
-
-   for(int i = 0; i < block_offsetsuml.Size(); i++)  
-   { 
-      block_offsetsuml[i] = block_offsetsumlz[i]; 
-   }
-   for(int i = 0; i < block_offsetsx.Size(); i++)    
-   { 
-      block_offsetsx[i] = block_offsetsuml[i] ; 
-   }
-
-  
-   ml = problem->Getml();
-  
-   lk.SetSize(dimC);  lk  = 0.0;
-   zlk.SetSize(dimM); zlk = 0.0;
-
-   linSolver = 0;
-   linSolveTol = 1.e-8;
-   MyRank = Mpi::WorldRank();
-   iAmRoot = MyRank == 0 ? true : false;
-}
-
-double ParInteriorPointSolver::MaxStepSize(Vector &x, Vector &xl, Vector &xhat, double tau)
-{
-   double alphaMaxloc = 1.0;
-   double alphaTmp;
-   for(int i = 0; i < x.Size(); i++)
-   {   
-      if( xhat(i) < 0. )
-      {
-         alphaTmp = -1. * tau * (x(i) - xl(i)) / xhat(i);
-         alphaMaxloc = min(alphaMaxloc, alphaTmp);
-      } 
-   }
-
-   // alphaMaxloc is the local maximum step size which is
-   // distinct on each MPI process. Need to compute
-   // the global maximum step size 
-   double alphaMaxglb;
-   MPI_Allreduce(&alphaMaxloc, &alphaMaxglb, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
-   return alphaMaxglb;
-}
-
-double ParInteriorPointSolver::MaxStepSize(Vector &x, Vector &xhat, double tau)
-{
-   Vector zero(x.Size()); zero = 0.0;
-   return MaxStepSize(x, zero, xhat, tau);
-}
-
-
-void ParInteriorPointSolver::Mult(const Vector &x0, Vector &xf)
-{
-   BlockVector x0block(block_offsetsx); x0block = 0.0;
-   x0block.GetBlock(0).Set(1.0, x0);
-   x0block.GetBlock(1) = 100.;
-   x0block.GetBlock(1).Add(1.0, ml);
-   BlockVector xfblock(block_offsetsx); xfblock = 0.0;
-   Mult(x0block, xfblock);
-   xf.Set(1.0, xfblock.GetBlock(0));
-}
-
-
-void ParInteriorPointSolver::Mult(const BlockVector &x0, BlockVector &xf)
-{
-   converged = false;
-   
-   BlockVector xk(block_offsetsx), xhat(block_offsetsx); xk = 0; xhat = 0.0;
-   BlockVector Xk(block_offsetsumlz), Xhat(block_offsetsumlz); Xk = 0.0; Xhat = 0.0;
-   BlockVector Xhatuml(block_offsetsuml); Xhatuml = 0.0;
-   Vector zlhat(dimM); zlhat = 0.0;
-
-   xk.GetBlock(0).Set(1.0, x0.GetBlock(0));
-   xk.GetBlock(1).Set(1.0, x0.GetBlock(1));
-   // running estimate of the final values of the Lagrange multipliers
-   lk  = 0.0;
-   zlk = 0.0;
-
-   for(int i = 0; i < dimM; i++)
-   {
-      zlk(i) = 1.e1 * mu_k / (xk(i+dimU) - ml(i));
-   }
-
-   Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
-   Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
-   Xk.GetBlock(2).Set(1.0, lk);
-   Xk.GetBlock(3).Set(1.0, zlk);
-
-  /* set theta0 = theta(x0)
-   *     thetaMin
-   *     thetaMax
-   * when theta(xk) < thetaMin and the switching condition holds
-   * then we ask for the Armijo sufficient decrease of the barrier
-   * objective to be satisfied, in order to accept the trial step length alphakl
-   * 
-   * thetaMax controls how the filter is initialized for each log-barrier subproblem
-   * F0 = {(th, phi) s.t. th > thetaMax}
-   * that is the filter does not allow for iterates where the constraint violation
-   * is larger than that of thetaMax
-   */
-   double theta0 = theta(xk);
-   thetaMin = 1.e-4 * max(1.0, theta0);
-   thetaMax = 1.e8  * thetaMin; // 1.e4 * max(1.0, theta0)
-
-   double Eeval, maxBarrierSolves, Eevalmu0;
-   bool printOptimalityError; // control optimality error print to console for log-barrier subproblems
-   
-   maxBarrierSolves = 10;
-
-   for(jOpt = 0; jOpt < max_iter; jOpt++)
-   {
-      if(iAmRoot)
-      {
-         cout << "interior-point solve step " << jOpt << endl;
-      }
-      // A-2. Check convergence of overall optimization problem
-      printOptimalityError = false;
-      Eevalmu0 = E(xk, lk, zlk, printOptimalityError);
-      if(Eevalmu0 < OptTol)
-      {
-         converged = true;
-         if(iAmRoot)
-         {
-            cout << "solved optimization problem :)\n";
-         }
-         break;
-      }
-      
-      if(jOpt > 0) { maxBarrierSolves = 1; }
-      
-      for(int i = 0; i < maxBarrierSolves; i++)
-      {
-         // A-3. Check convergence of the barrier subproblem
-         printOptimalityError = true;
-         Eeval = E(xk, lk, zlk, mu_k, printOptimalityError);
-         if(iAmRoot)
-         {
-            cout << "E = " << Eeval << endl;
-         }
-         if(Eeval < kEps * mu_k)
-         {
-            if(iAmRoot)
-            {
-               cout << "solved barrier subproblem :), for mu = " << mu_k << endl;
-            }
-            // A-3.1. Recompute the barrier parameter
-            mu_k  = max(OptTol / 10., min(kMu * mu_k, pow(mu_k, thetaMu)));
-            // A-3.2. Re-initialize the filter
-            F1.DeleteAll();
-            F2.DeleteAll();
-         }
-         else
-         {
-            break;
-         }
-      }
-    
-      // A-4. Compute the search direction
-      // solve for (uhat, mhat, lhat)
-      if(iAmRoot)
-      {
-         cout << "\n** A-4. IP-Newton solve **\n";
-      }
-      zlhat = 0.0; Xhatuml = 0.0;
-      // why do we have Xhatuml ....???
-      // TO DO: remove Xhatuml in favor of passing Xhat
-      IPNewtonSolve(xk, lk, zlk, zlhat, Xhatuml, mu_k, false); 
-
-      // assign data stack, X = (u, m, l, zl)
-      Xk = 0.0;
-      Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
-      Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
-      Xk.GetBlock(2).Set(1.0, lk);
-      Xk.GetBlock(3).Set(1.0, zlk);
-
-      // assign data stack, Xhat = (uhat, mhat, lhat, zlhat)
-      Xhat = 0.0;
-      for(int i = 0; i < 3; i++)
-      {
-         Xhat.GetBlock(i).Set(1.0, Xhatuml.GetBlock(i));
-      }
-      Xhat.GetBlock(3).Set(1.0, zlhat);
-
-      // A-5. Backtracking line search.
-      if(iAmRoot)
-      {
-         cout << "\n** A-5. Linesearch **\n";
-         cout << "mu = " << mu_k << endl;
-      }
-      lineSearch(Xk, Xhat, mu_k);
-
-      if(lineSearchSuccess)
-      {
-         if(iAmRoot)
-         {
-            cout << "lineSearch successful :)\n";
-         }
-         if(!switchCondition || !sufficientDecrease)
-         {
-            F1.Append( (1. - gTheta) * thx0);
-            F2.Append( phx0 - gPhi * thx0);
-         }
-         // ----- A-6: Accept the trial point
-         // print info regarding zl...
-         xk.GetBlock(0).Add(alpha, Xhat.GetBlock(0));
-         xk.GetBlock(1).Add(alpha, Xhat.GetBlock(1));
-         lk.Add(alpha,   Xhat.GetBlock(2));
-         zlk.Add(alphaz, Xhat.GetBlock(3));
-         projectZ(xk, zlk, mu_k);
-      }
-      else
-      {
-         if(iAmRoot)
-         {
-            cout << "lineSearch not successful :(\n";
-            cout << "attempting feasibility restoration with theta = " << thx0 << endl;
-            cout << "no feasibility restoration implemented, exiting now \n";
-         }
-         break;
-      }
-      if(jOpt + 1 == max_iter && iAmRoot) 
-      {  
-         cout << "maximum optimization iterations :(\n";
-      }
-   }
-   // done with optimization routine, just reassign data to xf reference so
-   // that the application code has access to the optimal point
-   xf = 0.0;
-   xf.GetBlock(0).Set(1.0, xk.GetBlock(0));
-   xf.GetBlock(1).Set(1.0, xk.GetBlock(1));
-}
-
-void ParInteriorPointSolver::FormIPNewtonMat(BlockVector & x, Vector & l, Vector &zl, BlockOperator &Ak)
-{
-   // WARNING: Huu, Hum, Hmu, Hmm should all be Hessian terms of the Lagrangian, currently we 
-   //          them by Hessian terms of the objective function and neglect the Hessian of l^T c
-
-   Huu = problem->Duuf(x); 
-   Hum = problem->Dumf(x);
-   Hmu = problem->Dmuf(x);
-   Hmm = problem->Dmmf(x);
-
-   Vector DiagLogBar(dimM); DiagLogBar = 0.0;
-   for(int ii = 0; ii < dimM; ii++)
-   {
-      DiagLogBar(ii) = zl(ii) / (x(ii+dimU) - ml(ii));
-   }
-   if(saveLogBarrierIterates)
-   {
-      std::ofstream diagStream;
-      char diagString[100];
-      snprintf(diagString, 100, "logBarrierHessiandata/D%d.dat", jOpt);
-      diagStream.open(diagString, ios::out | ios::trunc);
-      for(int ii = 0; ii < dimM; ii++)
-      {
-         diagStream << setprecision(30) << DiagLogBar(ii) << endl;
-      }
-      diagStream.close();
-   } 
-
-   SparseMatrix * Ds = new SparseMatrix(DiagLogBar);
-   ParFiniteElementSpace * fes = problem->GetfesM();
-   D = new HypreParMatrix(fes->GetComm(), fes->GlobalTrueVSize(), fes->GetTrueDofOffsets(), Ds);
-   HypreStealOwnership(*D,*Ds);
-   delete Ds;
-  
-   if(Hmm != nullptr)
-   {
-      Wmm = Hmm;
-      Wmm->Add(1.0, *D);
-   }
-   else
-   {
-      Wmm = D;
-   }
-
-   Ju = problem->Duc(x); JuT = Ju->Transpose();
-   Jm = problem->Dmc(x); JmT = Jm->Transpose();
-  
-   //         IP-Newton system matrix
-   //    Ak = [[H_(u,u)  H_(u,m)   J_u^T]
-   //          [H_(m,u)  W_(m,m)   J_m^T]
-   //          [ J_u      J_m       0  ]]
-
-   Ak.SetBlock(0, 0, Huu);                         Ak.SetBlock(0, 2, JuT);
-                           Ak.SetBlock(1, 1, Wmm); Ak.SetBlock(1, 2, JmT);
-   Ak.SetBlock(2, 0,  Ju); Ak.SetBlock(2, 1,  Jm);
-
-   if(Hum != nullptr) { Ak.SetBlock(0, 1, Hum); Ak.SetBlock(1, 0, Hmu); }
-}
-
-// perturbed KKT system solve
-// determine the search direction
-void ParInteriorPointSolver::IPNewtonSolve(BlockVector &x, Vector &l, Vector &zl, Vector &zlhat, BlockVector &Xhat, double mu, bool socSolve)
-{
-   // solve A x = b, where A is the IP-Newton matrix
-   BlockOperator A(block_offsetsuml, block_offsetsuml); BlockVector b(block_offsetsuml); b = 0.0;
-   FormIPNewtonMat(x, l, zl, A);
-
-   //       [grad_u phi + Ju^T l]
-   // b = - [grad_m phi + Jm^T l]
-   //       [          c        ]
-   BlockVector gradphi(block_offsetsx); gradphi = 0.0;
-   BlockVector JTl(block_offsetsx); JTl = 0.0;
-   Dxphi(x, mu, gradphi);
-   
-   (A.GetBlock(0,2)).Mult(l, JTl.GetBlock(0));
-   (A.GetBlock(1,2)).Mult(l, JTl.GetBlock(1));
-
-   for(int ii = 0; ii < 2; ii++)
-   {
-      b.GetBlock(ii).Set(1.0, gradphi.GetBlock(ii));
-      b.GetBlock(ii).Add(1.0, JTl.GetBlock(ii));
-   }
-   if(!socSolve) 
-   {
-      problem->c(x, b.GetBlock(2));
-   }
-   else
-   {
-      b.GetBlock(2).Set(1.0, ckSoc);
-   }
-   b *= -1.0; 
-   Xhat = 0.0;
-
-
-   // Direct solver (default)
-   if(linSolver == 0)
-   {
-      Array2D<HypreParMatrix *> ABlockMatrix(3,3);
-      for(int ii = 0; ii < 3; ii++)
-      {
-      for(int jj = 0; jj < 3; jj++)
-      {
-         if(!A.IsZeroBlock(ii, jj))
-         {
-            ABlockMatrix(ii, jj) = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(ii, jj)));
-	 }
-	 else
-	 {
-	    ABlockMatrix(ii, jj) = nullptr;
-	 }
-      }
-      }
-      
-      HypreParMatrix * Ah = HypreParMatrixFromBlocks(ABlockMatrix);   
-      
-      /* direct solve of the 3x3 IP-Newton linear system */
-      #ifdef MFEM_USE_MUMPS
-        MUMPSSolver ASolver;
-        ASolver.SetPrintLevel(0);
-        ASolver.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_INDEFINITE);
-        ASolver.SetOperator(*Ah);
-        ASolver.Mult(b, Xhat);
-      #else 
-        #ifdef MFEM_USE_MKL_CPARDISO
-          CPardisoSolver ASolver(MPI_COMM_WORLD);
-          ASolver.SetOperator(*Ah);
-          ASolver.Mult(b, Xhat);
-	#else
-	  MFEM_VERIFY(false, "linSolver 0 will not work unless compiled with MUMPS or MKL");
-        #endif
-      #endif
-
-      delete Ah;
-   }
-   else if(linSolver == 1 || linSolver == 2)
-   {
-      // form A = Huu + Ju^T D Ju, Wmm = D for contact
-      HypreParMatrix * Huuloc = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(0, 0)));
-      HypreParMatrix * Wmmloc = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(1, 1)));
-      HypreParMatrix * Juloc  = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(2, 0)));
-      HypreParMatrix * JuTloc = dynamic_cast<HypreParMatrix *>(&(A.GetBlock(0, 2)));
-      
-      
-      HypreParMatrix *JuTDJu   = RAP(Wmmloc, Juloc);     // Ju^T D Ju
-      HypreParMatrix *Areduced = ParAdd(Huuloc, JuTDJu);  // Huu + Ju^T D Ju
-      /* prepare the reduced rhs */
-      // breduced = bu + Ju^T (bm + Wmm bl)
-      Vector breduced(dimU); breduced = 0.0;
-      Vector tempVec(dimM); tempVec = 0.0;
-      Wmmloc->Mult(b.GetBlock(2), tempVec);
-      tempVec.Add(1.0, b.GetBlock(1));
-      JuTloc->Mult(tempVec, breduced);
-      breduced.Add(1.0, b.GetBlock(0));
-      
-      if(linSolver == 1)
-      {
-         // setup the solver for the reduced linear system
-         #ifdef MFEM_USE_MUMPS
-	   MUMPSSolver AreducedSolver;   
-           AreducedSolver.SetPrintLevel(0);
-           AreducedSolver.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_INDEFINITE);
-	   AreducedSolver.SetOperator(*Areduced);
-	   AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-         #else 
-           #ifdef MFEM_USE_MKL_CPARDISO
-	     CPardisoSolver AreducedSolver(MPI_COMM_WORLD);
-	     AreducedSolver.SetOperator(*Areduced);
-	     AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-           #else
-	     MFEM_VERIFY(false, "linSolver 1 will not work unless compiled with MUMPS or MKL");
-           #endif
-         #endif
-      }
-      else
-      {
-         HyprePCG AreducedSolver(MPI_COMM_WORLD);
-	 AreducedSolver.SetOperator(*Areduced);
-	 HypreBoomerAMG AreducedPrec;
-         AreducedSolver.SetTol(linSolveTol);
-         AreducedSolver.SetMaxIter(500);
-         AreducedSolver.SetPreconditioner(AreducedPrec);
-         AreducedSolver.SetResidualConvergenceOptions(); // convergence criteria based on residual norm
-	 AreducedSolver.SetPrintLevel(2);
-         AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
-      }
-
-      // now propagate solved uhat to obtain mhat and lhat
-      // xm = Ju xu - bl
-      Juloc->Mult(Xhat.GetBlock(0), Xhat.GetBlock(1));
-      Xhat.GetBlock(1).Add(-1.0, b.GetBlock(2));
-
-      // xl = Wmm xm - bm
-      Wmmloc->Mult(Xhat.GetBlock(1), Xhat.GetBlock(2));
-      Xhat.GetBlock(2).Add(-1.0, b.GetBlock(1));
-
-      delete JuTDJu;
-      delete Areduced;
-   }
-
-   /* backsolve to determine zlhat */
-   for(int ii = 0; ii < dimM; ii++)
-   {
-      zlhat(ii) = -1.*(zl(ii) + (zl(ii) * Xhat(ii + dimU) - mu) / (x(ii + dimU) - ml(ii)) );
-   }
-
-   // free memory
-   delete D;
-   delete JuT;
-   delete JmT;
-   if(Hmm != nullptr)
-   {
-      delete Wmm;
-   }
-}
-
-// here Xhat, X will be BlockVectors w.r.t. the 4 partitioning X = (u, m, l, zl)
-
-void ParInteriorPointSolver::lineSearch(BlockVector& X0, BlockVector& Xhat, double mu)
-{
-   double tau  = max(tauMin, 1.0 - mu);
-   Vector u0   = X0.GetBlock(0);
-   Vector m0   = X0.GetBlock(1);
-   Vector l0   = X0.GetBlock(2);
-   Vector z0   = X0.GetBlock(3);
-   Vector uhat = Xhat.GetBlock(0);
-   Vector mhat = Xhat.GetBlock(1);
-   Vector lhat = Xhat.GetBlock(2);
-   Vector zhat = Xhat.GetBlock(3);
-   double alphaMax  = MaxStepSize(m0, ml, mhat, tau);
-   double alphaMaxz = MaxStepSize(z0, zhat, tau);
-   alphaz = alphaMaxz;
-
-   BlockVector x0(block_offsetsx); x0 = 0.0;
-   x0.GetBlock(0).Set(1.0, u0);
-   x0.GetBlock(1).Set(1.0, m0);
-   
-   BlockVector xhat(block_offsetsx); xhat = 0.0;
-   xhat.GetBlock(0).Set(1.0, uhat);
-   xhat.GetBlock(1).Set(1.0, mhat);
-   
-   BlockVector xtrial(block_offsetsx); xtrial = 0.0;
-   BlockVector Dxphi0(block_offsetsx); Dxphi0 = 0.0;
-   int maxBacktrack = 20;
-   alpha = alphaMax;
-
-   Vector ck0(dimC); ck0 = 0.0;
-   Vector zhatsoc(dimM); zhatsoc = 0.0;
-   BlockVector Xhatumlsoc(block_offsetsuml); Xhatumlsoc = 0.0;
-   BlockVector xhatsoc(block_offsetsx); xhatsoc = 0.0;
-   Vector uhatsoc(dimU); uhatsoc = 0.0;
-   Vector mhatsoc(dimM); mhatsoc = 0.0;
-
-   Dxphi(x0, mu, Dxphi0);
-
-   Dxphi0_xhat = InnerProduct(MPI_COMM_WORLD, Dxphi0, xhat);
-   descentDirection = Dxphi0_xhat < 0. ? true : false;
-   if(descentDirection)
-   {
-      if (iAmRoot)
-      {
-         cout << "is a descent direction for the log-barrier objective\n";
-      }
-   }
-   else
-   {
-      cout << "is not a descent direction for the log-barrier objective\n";
-   }
-   thx0 = theta(x0);
-   phx0 = phi(x0, mu);
-
-   lineSearchSuccess = false;
-   for(int i = 0; i < maxBacktrack; i++)
-   {
-      if (iAmRoot)
-      {
-         cout << "\n--------- alpha = " << alpha << " ---------\n";
-      }
-      // ----- A-5.2. Compute trial point: xtrial = x0 + alpha_i xhat
-      xtrial.Set(1.0, x0);
-      xtrial.Add(alpha, xhat);
-
-      // ------ A-5.3. if not in filter region go to A.5.4 otherwise go to A-5.5.
-      thxtrial = theta(xtrial);
-      phxtrial = phi(xtrial, mu);
-      filterCheck(thxtrial, phxtrial);    
-      if(!inFilterRegion)
-      {
-         if (iAmRoot)
-         {
-            cout << "not in filter region :)\n";
-         }
-         // ------ A.5.4: Check sufficient decrease
-         if(!descentDirection)
-         {
-            switchCondition = false;
-         }
-         else
-         {
-            switchCondition = (alpha * pow(abs(Dxphi0_xhat), sPhi) > delta * pow(thx0, sTheta)) ? true : false;
-         }
-         if (iAmRoot)
-         {
-            cout << "theta(x0) = "     << thx0     << ", thetaMin = "                  << thetaMin             << endl;
-            cout << "theta(xtrial) = " << thxtrial << ", (1-gTheta) *theta(x0) = "     << (1. - gTheta) * thx0 << endl;
-            cout << "phi(xtrial) = "   << phxtrial << ", phi(x0) - gPhi *theta(x0) = " << phx0 - gPhi * thx0   << endl;
-         }      
-         // Case I      
-         if(thx0 <= thetaMin && switchCondition)
-         {
-            sufficientDecrease = (phxtrial <= phx0 + eta * alpha * Dxphi0_xhat) ? true : false;
-            if(sufficientDecrease)
-            {
-               if(iAmRoot) { cout << "Line search successful: sufficient decrease in log-barrier objective.\n"; } 
-               // accept the trial step
-               lineSearchSuccess = true;
-               break;
-            }
-         }
-         else
-         {
-            if(thxtrial <= (1. - gTheta) * thx0 || phxtrial <= phx0 - gPhi * thx0)
-            {
-               if(iAmRoot) { cout << "Line search successful: infeasibility or log-barrier objective decreased.\n"; } 
-               // accept the trial step
-               lineSearchSuccess = true;
-               break;
-            }
-         }
-         // A-5.5: Initialize the second-order correction
-         if((!(thx0 < thxtrial)) && i == 0)
-         {
-            if (iAmRoot)
-            {
-               cout << "second order correction\n";
-            }
-            problem->c(xtrial, ckSoc);
-            problem->c(x0, ck0);
-            ckSoc.Add(alphaMax, ck0);
-            // A-5.6 Compute the second-order correction.
-            IPNewtonSolve(x0, l0, z0, zhatsoc, Xhatumlsoc, mu, true);
-            mhatsoc.Set(1.0, Xhatumlsoc.GetBlock(1));
-            //WARNING: not complete but currently solver isn't entering this region
-         }
-      }
-      else
-      {
-         if (iAmRoot)
-         {
-            cout << "in filter region :(\n"; 
-         }
-      }
-      // include more if needed
-      alpha *= 0.5;
-
-   } 
-}
-
-
-void ParInteriorPointSolver::projectZ(const Vector &x, Vector &z, double mu)
-{
-   double zi;
-   double mudivmml;
-   for(int i = 0; i < dimM; i++)
-   {
-      zi = z(i);
-      mudivmml = mu / (x(i + dimU) - ml(i));
-      z(i) = max(min(zi, kSig * mudivmml), mudivmml / kSig);
-   }
-}
-
-void ParInteriorPointSolver::filterCheck(double th, double ph)
-{
-   inFilterRegion = false;
-   if(th > thetaMax)
-   {
-      inFilterRegion = true;
-   }
-   else
-   {
-      for(int i = 0; i < F1.Size(); i++)
-      {
-         if(th >= F1[i] && ph >= F2[i])
-         {
-            inFilterRegion = true;
-            break;
-         }
-      }
-   }
-}
-
-double ParInteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, double mu, bool printEeval)
-{
-   double E1, E2, E3;
-   double sc, sd;
-   BlockVector gradL(block_offsetsx); gradL = 0.0; // stationarity grad L = grad f + J^T l - z
-   Vector cx(dimC); cx = 0.0;     // feasibility c = c(x)
-   Vector comp(dimM); comp = 0.0; // complementarity M Z - mu 1
-
-   DxL(x, l, zl, gradL);
-   E1 = GlobalLpNorm(infinity(), gradL.Normlinf(), MPI_COMM_WORLD); 
-
-   problem->c(x, cx);
-   E2 = GlobalLpNorm(infinity(), cx.Normlinf(), MPI_COMM_WORLD); 
-
-
-   for(int ii = 0; ii < dimM; ii++) 
-   { 
-      comp(ii) = x(dimU + ii) * zl(ii) - mu;
-   }
-   E3 = GlobalLpNorm(infinity(), comp.Normlinf(), MPI_COMM_WORLD); 
-
-   double ll1, zl1;
-
-   zl1 = GlobalLpNorm(1, zl.Norml1(), MPI_COMM_WORLD)/ double(dimCglb + dimMglb); 
-   ll1 = GlobalLpNorm(1, l.Norml1(), MPI_COMM_WORLD);
-   sc = max(sMax, zl1 / (double(dimMglb)) ) / sMax;
-   sd = max(sMax, (ll1 + zl1) / (double(dimCglb + dimMglb))) / sMax;
-   if(iAmRoot && printEeval)
-   {
-      cout << "evaluating optimality error for mu = " << mu << endl;
-      cout << "stationarity measure = "    << E1 / sd << endl;
-      cout << "feasibility measure  = "    << E2      << endl;
-      cout << "complimentarity measure = " << E3 / sc << endl;
-   }
-   return max(max(E1 / sd, E2), E3 / sc);
-}
-
-double ParInteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, bool printEeval)
-{
-  return E(x, l, zl, 0.0, printEeval);
-}
-
-double ParInteriorPointSolver::theta(const BlockVector &x)
-{
-  Vector cx(dimC); cx = 0.0;
-  problem->c(x, cx);
-  return sqrt(InnerProduct(MPI_COMM_WORLD,cx, cx));
-}
-
-// log-barrier objective
-double ParInteriorPointSolver::phi(const BlockVector &x, double mu)
-{
-   double fx = problem->CalcObjective(x); 
-   double logBarrierLoc = 0.0;
-   for(int i = 0; i < dimM; i++) 
-   { 
-     logBarrierLoc += log(x(dimU+i)-ml(i));
-   }
-   double logBarrierGlb;
-   MPI_Allreduce(&logBarrierLoc, &logBarrierGlb, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-   return fx - mu * logBarrierGlb;
-}
-
-// gradient of log-barrier objective with respect to x = (u, m)
-void ParInteriorPointSolver::Dxphi(const BlockVector &x, double mu, BlockVector &y)
-{
-   problem->CalcObjectiveGrad(x, y);
-   
-   for(int i = 0; i < dimM; i++) 
-   { 
-      y(dimU + i) -= mu / (x(dimU + i));
-   } 
-}
-
-// Lagrangian function evaluation
-// L(x, l, zl) = f(x) + l^T c(x) - zl^T m
-double ParInteriorPointSolver::L(const BlockVector &x, const Vector &l, const Vector &zl)
-{
-   double fx = problem->CalcObjective(x);
-   Vector cx(dimC); problem->c(x, cx);
-   return (fx + InnerProduct(MPI_COMM_WORLD,cx, l) - InnerProduct(MPI_COMM_WORLD, x.GetBlock(1), zl));
-}
-
-void ParInteriorPointSolver::DxL(const BlockVector &x, const Vector &l, const Vector &zl, BlockVector &y)
-{
-   // evaluate the gradient of the objective with respect to the primal variables x = (u, m)
-   BlockVector gradxf(block_offsetsx); gradxf = 0.0;
-   problem->CalcObjectiveGrad(x, gradxf);
-   
-   HypreParMatrix *Jacu, *Jacm, *JacuT, *JacmT;
-   Jacu = problem->Duc(x); 
-   Jacm = problem->Dmc(x);
-   JacuT = Jacu->Transpose();
-   JacmT = Jacm->Transpose();
-   
-   JacuT->Mult(l, y.GetBlock(0));
-   JacmT->Mult(l, y.GetBlock(1));
-   
-   delete JacuT;
-   delete JacmT;
-   
-   y.Add(1.0, gradxf);
-   (y.GetBlock(1)).Add(-1.0, zl);
-}
-
-bool ParInteriorPointSolver::GetConverged() const
-{
-   return converged;
-}
-
-void ParInteriorPointSolver::SetTol(double Tol)
-{
-   OptTol = Tol;
-}
-
-void ParInteriorPointSolver::SetMaxIter(int max_it)
-{
-   max_iter = max_it;
-}
-
-void ParInteriorPointSolver::SetBarrierParameter(double mu_0)
-{
-   mu_k = mu_0;
-}
-
-void ParInteriorPointSolver::SaveLogBarrierHessianIterates(bool save)
-{
-   MFEM_ASSERT(MyRank == 0 || save == false, "currently can only save logbarrier hessian in serial codes");
-   saveLogBarrierIterates = save;
-}
-
-void ParInteriorPointSolver::SetLinearSolver(int LinSolver)
-{
-   linSolver = LinSolver;
-}
-
-void ParInteriorPointSolver::SetLinearSolveTol(double Tol)
-{
-  linSolveTol = Tol;
-}
-
-
-ParInteriorPointSolver::~ParInteriorPointSolver() 
-{
-   F1.DeleteAll();
-   F2.DeleteAll();
-   block_offsetsx.DeleteAll();
-   block_offsetsumlz.DeleteAll();
-   block_offsetsuml.DeleteAll();
-   ml.SetSize(0);
-}
@@ -1,82 +0,0 @@
-#include "mfem.hpp"
-#include "ParProblems.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-
-#ifndef PARIPSOLVER 
-#define PARIPSOLVER
-
-class ParInteriorPointSolver
-{
-protected:
-    ParGeneralOptProblem* problem;
-    double OptTol;
-    int  max_iter;
-    double mu_k; // \mu_k
-    Vector lk, zlk;
-
-    double sMax, kSig, tauMin, eta, thetaMin, delta, sTheta, sPhi, kMu, thetaMu;
-    double thetaMax, kSoc, gTheta, gPhi, kEps;
-	
-    // filter
-    Array<double> F1, F2;
-	
-    // quantities computed in lineSearch
-    double alpha, alphaz;
-    double thx0, thxtrial;
-    double phx0, phxtrial;
-    bool descentDirection, switchCondition, sufficientDecrease, lineSearchSuccess, inFilterRegion;
-    double Dxphi0_xhat;
-
-    int dimU, dimM, dimC;
-    int dimUglb, dimMglb, dimCglb;
-    Array<int> block_offsetsumlz, block_offsetsuml, block_offsetsx;
-    Vector ml;
-
-    Vector ckSoc;
-    HypreParMatrix * Huu, * Hum, * Hmu, * Hmm, * Wmm, *D, * Ju, * Jm, * JuT, * JmT;
-    
-    int jOpt;
-    bool converged;
-    
-    int MyRank;
-    bool iAmRoot;
-
-    bool saveLogBarrierIterates;
-
-    int linSolver;
-    double linSolveTol;
-public:
-    ParInteriorPointSolver(ParGeneralOptProblem*);
-    double MaxStepSize(Vector& , Vector& , Vector& , double);
-    double MaxStepSize(Vector& , Vector& , double);
-    void Mult(const BlockVector& , BlockVector&);
-    void Mult(const Vector&, Vector &); 
-    void FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
-    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , double, bool);
-    void lineSearch(BlockVector& , BlockVector& , double);
-    void projectZ(const Vector & , Vector &, double);
-    void filterCheck(double, double);
-    double E(const BlockVector &, const Vector &, const Vector &, double, bool);
-    double E(const BlockVector &, const Vector &, const Vector &, bool);
-    bool GetConverged() const;
-    // TO DO: include Hessian of Lagrangian
-    double theta(const BlockVector &);
-    double phi(const BlockVector &, double);
-    void Dxphi(const BlockVector &, double, BlockVector &);
-    double L(const BlockVector &, const Vector &, const Vector &);
-    void DxL(const BlockVector &, const Vector &, const Vector &, BlockVector &);
-    void SetTol(double);
-    void SetMaxIter(int);
-    void SetBarrierParameter(double);    
-    void SaveLogBarrierHessianIterates(bool);
-    void SetLinearSolver(int);
-    void SetLinearSolveTol(double);
-    virtual ~ParInteriorPointSolver();
-};
-
-#endif
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Socratis Petrides	bc39660801	minor bug	2023-05-31 11:00:58 -07:00
Frank Wang	a72cd6b7a8	remove comments	2023-04-07 18:43:48 -07:00
Frank Wang	bcc308b263	fixing bugs	2023-04-07 18:43:09 -07:00
nychiang	23b5e402e3	clean code	2023-04-04 16:40:28 -07:00
nychiang	ee0dd42064	code style	2023-04-04 12:14:57 -07:00
nychiang	39a78ff3df	IPOPT in mfem	2023-04-04 00:09:53 -07:00
Frank Wang	20460cafda	add bc	2023-03-31 10:57:16 -07:00
Frank Wang	2f188e9fcf	fix bug	2023-03-29 10:29:51 -07:00
Frank Wang	e772d71ecb	update	2023-03-24 21:43:41 -07:00
Frank Wang	d9c4cad220	latest update	2023-03-21 12:35:20 -07:00
Frank Wang	0c572d4779	comment out nodepair for now	2023-03-16 13:44:16 -07:00
Frank Wang	4320ad2272	adding things needed for Jacobian computation	2023-03-16 13:42:17 -07:00
Dylan Copeland	5c09d2b4c9	Added computation of face reference coordinates, as well as the global vertex indices corresponding to the corners of the face.	2023-03-15 20:16:36 -07:00
Dylan Copeland	65e75e271b	Enabling an example with points outside domain 1.	2023-03-15 17:35:55 -07:00
Frank Wang	f8e5fe77e6	update contact	2023-03-12 19:07:25 -07:00
Frank Wang	4ff880f0cf	remove binary	2023-03-10 11:04:44 -08:00
Frank Wang	a1f57a6375	adding x field	2023-03-10 10:28:01 -08:00
Dylan Copeland	5ffb605333	Fixing contact example.	2023-03-08 18:52:04 -08:00
Frank Wang	f8ea695e13	add contactcpp	2023-03-08 12:02:01 -08:00