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

removing bug in the evaluation of the optimality measures wherein all processors now see the same value of the optimality error
clean up... using a GeneralOptimizationProblem class and an OptimizationProblem class and no longer using the specific ContactProblem class as the generality of the Parent class does not have anything per say to do with contact
2024-03-01 14:58:16 -08:00 · 2023-11-10 10:43:20 -08:00 · 2023-10-10 11:33:29 -07:00 · 2023-08-07 10:07:55 -07:00 · 2023-08-07 10:06:35 -07:00 · 2023-07-25 19:20:04 -07:00
869 changed files with 18416 additions and 5589 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
+- cmake --build build --target metis -j 4
 - 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
+- cmake --build build -j 4
 - cmake --build build --target install
 - cd ../..

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

 # MFEM
 before_build:
- 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
+- 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

 build_script:
- 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
+- 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

 after_build:
- cd build_serial
+# - cd build_serial
+# - ctest -C Release --output-on-failure
+# - cd ..
+- cd build_shared_serial
 - ctest -C Release --output-on-failure
@@ -62,7 +62,7 @@ jobs:

      - name: GHCR Login
        if: (github.event_name != 'pull_request')
-        uses: docker/login-action@v1
+        uses: docker/login-action@v2
        with:
          registry: ghcr.io
          username: ${{ github.actor }}
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -74,12 +74,19 @@ 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
@@ -94,7 +101,7 @@ jobs:
    # This external action allows to interrupt a workflow already running on
    # the same branch to save resource
    - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.9.0
+      uses: styfle/cancel-workflow-action@0.11.0
      with:
        access_token: ${{ github.token }}

@@ -102,7 +109,7 @@ jobs:
    # /home/runner/work/mfem/mfem/mfem
    # Note: Done now to access "install-hypre" and "install-metis" actions.
    - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
      with:
        path: ${{ env.MFEM_TOP_DIR }}
        # Fetch the complete history for codecov to access commits ID
@@ -115,25 +122,25 @@ jobs:
      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-latest'
      run: |
        sudo apt-get install mpich libmpich-dev
-        export MAKE_CXX_FLAG="MPICXX=mpic++"

    - name: get lcov (Linux)
      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-latest'
      run: |
        sudo apt-get install lcov

-    - name: Set up Homebrew
-      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
-      uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3
+    # Keep the following section in case we need it again in the future,
+    # see: https://github.com/mfem/mfem/pull/3385#discussion_r1058013032
+    # - name: Set up Homebrew
+    #   if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
+    #   uses: Homebrew/actions/setup-homebrew@master

    - name: get MPI (MacOS)
      if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
      run: |
        export HOMEBREW_NO_INSTALL_CLEANUP=1
        brew install openmpi
-        export MAKE_CXX_FLAG="MPICXX=mpic++"

-    - name: get MPI (MacOS)
+    - name: get lcov (MacOS)
      if: matrix.codecov == 'YES' && matrix.os == 'macos-latest'
      run: |
        export HOMEBREW_NO_INSTALL_CLEANUP=1
@@ -141,14 +148,14 @@ jobs:

    - name: get MPI (Windows)
      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
-      uses: mpi4py/setup-mpi@v1.0.3
+      uses: mpi4py/setup-mpi@v1.1.4

    # Get Hypre through cache, or build it.
    # Install will only run on cache miss.
    - name: cache hypre
      id: hypre-cache
      if: matrix.mpi == 'par'
-      uses: actions/cache@v2
+      uses: actions/cache@v3
      with:
        path: ${{ env.HYPRE_TOP_DIR }}
        key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2
@@ -176,7 +183,7 @@ jobs:
    - name: cache metis
      id: metis-cache
      if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
-      uses: actions/cache@v2
+      uses: actions/cache@v3
      with:
        path: ${{ env.METIS_TOP_DIR }}
        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
@@ -190,12 +197,13 @@ jobs:

    - name: cache vcpkg (Windows)
      id: vcpkg-cache
+      if: matrix.os == 'windows-latest'
      uses: actions/cache@v3
      with:
        path: vcpkg_cache
        key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1

-    - name: prepare binary cache location
+    - name: prepare vcpkg binary cache location (Windows)
      if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
      run: |
        mkdir -p vcpkg_cache
@@ -209,7 +217,7 @@ jobs:

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

    # Run checks (and only checks) on debug targets
@@ -230,11 +238,7 @@ jobs:
      run: |
        cd ${{ env.MFEM_TOP_DIR }}  && make check

-    - name: unit tests
-      if: matrix.build-system == 'make' && matrix.target == 'opt'
-      run: |
-        cd ${{ env.MFEM_TOP_DIR }} && make unittest
-
+    # Note: 'tests' include the unit tests
    - name: tests
      if: matrix.build-system == 'make' && matrix.target == 'opt'
      run: |
@@ -247,7 +251,7 @@ jobs:
        cd ${{ env.MFEM_TOP_DIR }} && cmake --build build --target check --config ${CTEST_CONFIG}
      shell: bash

-    - name: cmake unit tests (Ubuntu 20.04)
+    - name: cmake unit tests (Ubuntu)
      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
      run: |
        CTEST_CONFIG="Release"
@@ -265,7 +269,7 @@ jobs:
    # Code coverage (process and upload reports)
    - name: codecov
      if: matrix.codecov == 'YES'
-      uses: mfem/github-actions/upload-coverage@v2.0
+      uses: mfem/github-actions/upload-coverage@v2.2
      with:
        name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
        project_dir: ${{ env.MFEM_TOP_DIR }}
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -35,23 +35,22 @@ jobs:

    steps:
    - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.9.0
+      uses: styfle/cancel-workflow-action@0.11.0
      with:
        access_token: ${{ github.token }}

    - name: checkout MFEM
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
      with:
        path: mfem

    - name: Get MPI (Linux)
      run: |
        sudo apt-get install mpich libmpich-dev
-        export MAKE_CXX_FLAG="MPICXX=mpic++"

    - name: Cache Hypre Install
      id: hypre-cache
-      uses: actions/cache@v2
+      uses: actions/cache@v3
      with:
        path: ${{ env.HYPRE_TOP_DIR }}
        key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2.2
@@ -66,7 +65,7 @@ jobs:

    - name: Cache Metis Install
      id: metis-cache
-      uses: actions/cache@v2
+      uses: actions/cache@v3
      with:
        path: ${{ env.METIS_TOP_DIR }}
        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -34,12 +34,12 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.9.0
+      uses: styfle/cancel-workflow-action@0.11.0
      with:
        access_token: ${{ github.token }}

    - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3

    - name: copyright check
      id: copyright
@@ -84,7 +84,7 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3

    - name: get astyle
      run: |
@@ -101,11 +101,14 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3

    - name: get doxygen and graphviz
      run: |
        sudo apt-get install doxygen graphviz
+
+    - name: update doxygen config file
+      run: |
        cd doc
        doxygen -u CodeDocumentation.conf.in

@@ -123,7 +126,7 @@ jobs:
    runs-on: ubuntu-latest
    steps:
    - name: checkout mfem
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
      with:
        fetch-depth: 0

@@ -21,6 +21,9 @@ CMakeFiles/
 # Default install location
 /mfem/

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

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

 examples/hiop/ex9
 examples/hiop/ex9p
@@ -148,8 +149,8 @@ examples/petsc/velocity.*
 examples/petsc/elastic_energy.*
 examples/petsc/mode_*

-examples/pumi/ex1
-examples/pumi/ex[126]p
+examples/pumi/ex[12]
+examples/pumi/ex[16]p
 examples/pumi/refined.mesh
 examples/pumi/sol.gf
 examples/pumi/mesh.*
@@ -194,6 +195,9 @@ 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
@@ -203,8 +207,10 @@ miniapps/meshing/mesh-explorer
 miniapps/meshing/shaper
 miniapps/meshing/extruder
 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
@@ -214,9 +220,12 @@ miniapps/meshing/toroid-*.mesh
 miniapps/meshing/twist-*.mesh
 miniapps/meshing/mesh-explorer.mesh
 miniapps/meshing/partitioning.txt
+miniapps/meshing/mesh-explorer-visit*
+miniapps/meshing/mesh-explorer-paraview/
 miniapps/meshing/shaper.mesh
 miniapps/meshing/extruder.mesh
 miniapps/meshing/trimmer.mesh
+miniapps/meshing/reflected.mesh
 miniapps/meshing/optimized*
 miniapps/meshing/perturbed*
 miniapps/meshing/polar-nc.mesh
@@ -241,6 +250,7 @@ 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

@@ -276,6 +286,7 @@ miniapps/tools/convert-dc
 miniapps/tools/lor-transfer
 miniapps/tools/get-values
 miniapps/tools/check-tmop-metric
+miniapps/tools/tmop-metric-magnitude

 miniapps/toys/automata
 miniapps/toys/life
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -93,7 +93,7 @@ report_baseline:
          git pull && \
          git add ${rundir} && \
          git commit -m "${msg}" && \
-          git push origin master
+          ${CI_PROJECT_DIR}/.gitlab/scripts/git_try_to_push
        else
          for file in ${rundir}/*; do
            echo "------------------------------"
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -27,11 +27,6 @@ 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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -0,0 +1,41 @@
+#!/bin/bash
+
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Try to push to the remote 5 times. If the push fails, and the local and remote
+# have diverged, then pull from the remote to merge changes, and try pushing
+# again. If some other failure happens
+for i in {1..5}; do
+   git push origin master && exit 0
+   # Wait for 20 seconds in case someone else is pushing to the remote
+   # concurrently
+   sleep 20
+   # Fetch any updates from the remote
+   git remote update
+   # Get the latest commit on the local branch
+   LOCAL=$(git rev-parse @)
+   # Get the latest commit on the remote
+   REMOTE=$(git rev-parse @{u})
+   # Get the common ancestor
+   BASE=$(git merge-base @ @{u})
+   # Have the local and remote diverged?
+   if [[ $LOCAL != $REMOTE && $LOCAL != $BASE && $REMOTE != $BASE ]]; then
+      git pull
+      if [[ $? == 0 ]]; then
+         continue
+      else
+         exit 1 # Something else went wrong trying to pull
+      fi
+   fi
+done
+
+exit 1 # Did not succeed in 5 attempts
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -32,7 +32,7 @@ if [[ "$AUTOTEST_COMMIT" != "NO" ]]; then
  git pull && \
  git add ${rundir} && \
  git commit -m "${msg}" && \
-  git push origin master
+  ${CI_PROJECT_DIR}/.gitlab/scripts/git_try_to_push
 else
  for file in ${rundir}/*; do
    echo "------------------------------"
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -29,7 +29,7 @@ if [[ "$AUTOTEST_COMMIT" != "NO" ]]; then
  git pull && \
  git add ${rundir} && \
  git commit -m "${msg}" && \
-  git push origin master
+  ${CI_PROJECT_DIR}/.gitlab/scripts/git_try_to_push
 else
  for file in ${rundir}/*; do
    echo "------------------------------"
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -8,6 +8,84 @@
                               https://mfem.org


+Version 4.5.3 (development)
+===========================
+
+- 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
+=========================================
+
+- 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.
+
+- 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 support for GridFunction::GetGradients() and GetVectorGradient() on
+  face-neighbor elements.
+
+- Added support for pyramids in Gmsh meshes.
+
+- 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.
+
+- VisItDataCollection now correctly handles data collection names containing
+  underscores.
+
+- 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.
+
+- The following integrations have updated minimum version requirements:
+  * RAJA >= 2022.10.3
+
+API changes
+-----------
+- The implicit cast methods of class Vector to 'double *' and 'const double *'
+  have been deprecated and generate deprecation warnings if used. They will be
+  removed in a future release.
+
+- The methods Mesh::GetFaceBaseGeometry and Mesh::GetFaceGeometryType have been
+  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
 =========================================

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -10,7 +10,10 @@
 # CONTRIBUTING.md for details.

 # The variable CMAKE_CXX_STANDARD and related were introduced in CMake v3.1
-cmake_minimum_required(VERSION 3.1)
+# 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.")

@@ -22,8 +25,11 @@ 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.
-message(STATUS "(optional) USER_CONFIG = ${USER_CONFIG}")
-include("${USER_CONFIG}" OPTIONAL)
+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})")
 include("${CMAKE_CURRENT_SOURCE_DIR}/config/defaults.cmake")

 # Allow overwriting of the compiler by setting CXX/MPICXX on the command line or
@@ -51,7 +57,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.0)
+set(${PROJECT_NAME}_VERSION 4.5.3)

 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -66,8 +72,7 @@ set(XSDK_ENABLE_C OFF)
 set(XSDK_ENABLE_Fortran OFF)

 # Check if we need to enable C or Fortran.
-if (CMAKE_VERSION VERSION_LESS 3.2 OR
-    MFEM_USE_CONDUIT OR
+if (MFEM_USE_CONDUIT OR
    MFEM_USE_SIDRE OR
    MFEM_USE_PETSC)
   # This seems to be needed by:
@@ -81,11 +86,13 @@ if (MFEM_USE_STRUMPACK)
  # Just needed to find the MPI_Fortran libraries to link with
  set(XSDK_ENABLE_Fortran ON)
 endif()
-# SUNDIALS >= 6.4.0 requires C++14:
-if (MFEM_USE_SUNDIALS AND ("${CMAKE_CXX_STANDARD}" LESS "14"))
-  set(CMAKE_CXX_STANDARD 14)
-endif()
-if (MFEM_USE_GINKGO AND ("${CMAKE_CXX_STANDARD}" LESS "14"))
+# SUNDIALS, STRUMPACK, Ginkgo, RAJA and Umpire require C++14:
+if ((MFEM_USE_SUNDIALS OR
+     MFEM_USE_STRUMPACK OR
+     MFEM_USE_GINKGO OR
+     MFEM_USE_RAJA OR
+     MFEM_USE_UMPIRE) AND
+    ("${CMAKE_CXX_STANDARD}" LESS "14"))
  set(CMAKE_CXX_STANDARD 14)
 endif()

@@ -103,17 +110,12 @@ if (MFEM_USE_CUDA)
   if (MFEM_USE_HIP)
      message(FATAL_ERROR " *** MFEM_USE_HIP cannot be combined with MFEM_USE_CUDA.")
   endif()
-   # MFEM_USE_CUDA requires CMake 3.8 or newer (for direct CUDA support)
-   cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
   # Use ${CMAKE_CXX_COMPILER} as the cuda host compiler.
   if (NOT CMAKE_CUDA_HOST_COMPILER)
      set(CMAKE_CUDA_HOST_COMPILER ${CMAKE_CXX_COMPILER})
   endif()
   enable_language(CUDA)
-   set(CMAKE_CUDA_STANDARD 11)
-   if (MFEM_USE_GINKGO)
-     set(CMAKE_CUDA_STANDARD 14)
-   endif()
+   set(CMAKE_CUDA_STANDARD ${CMAKE_CXX_STANDARD})
   set(CMAKE_CUDA_STANDARD_REQUIRED ON)
   set(CMAKE_CUDA_EXTENSIONS OFF)
   set(CUDA_FLAGS "--expt-extended-lambda")
@@ -136,8 +138,7 @@ if (MFEM_USE_CUDA)
      set(CUDA_FLAGS "-ccbin=${CMAKE_CXX_COMPILER} ${CUDA_FLAGS}")
      set(CMAKE_CUDA_HOST_LINK_LAUNCHER ${CMAKE_CXX_COMPILER})
   endif()
-   set(CMAKE_CUDA_FLAGS "${CUDA_FLAGS}" CACHE STRING
-       "CUDA flags set for MFEM" FORCE)
+   set(CMAKE_CUDA_FLAGS ${CMAKE_CUDA_FLAGS} ${CUDA_FLAGS})
   set(CUSPARSE_FOUND TRUE)
   set(CUSPARSE_LIBRARIES "cusparse")
   set(CUBLAS_FOUND TRUE)
@@ -186,6 +187,23 @@ 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)
@@ -199,6 +217,26 @@ if (MFEM_USE_HIP)
  find_package(HIPSPARSE REQUIRED)
 endif()

+# OpenMP
+if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
+  if (NOT MFEM_THREAD_SAFE AND MFEM_USE_LEGACY_OPENMP)
+    message(FATAL_ERROR " *** MFEM_USE_LEGACY_OPENMP requires MFEM_THREAD_SAFE=ON.")
+  endif()
+  find_package(OpenMP REQUIRED)
+  set(OPENMP_LIBRARIES ${OpenMP_CXX_LIBRARIES})
+  if(APPLE)
+    # On macOS, the compiler needs additional help to find the <omp.h> header.
+    # See issue #2642 for more information.
+    set(OPENMP_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS})
+  endif(APPLE)
+  if (OPENMP_FOUND)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+    if (MFEM_USE_CUDA)
+      set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -Xcompiler=${OpenMP_CXX_FLAGS}")
+    endif()
+  endif()
+endif()
+
 # MPI -> hypre; PETSc (optional)
 if (MFEM_USE_MPI)
  find_package(MPI REQUIRED)
@@ -264,20 +302,6 @@ if (MFEM_USE_LAPACK)
  find_package(LAPACK REQUIRED)
 endif()

-# OpenMP
-if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
-  if (NOT MFEM_THREAD_SAFE AND MFEM_USE_LEGACY_OPENMP)
-    message(FATAL_ERROR " *** MFEM_USE_LEGACY_OPENMP requires MFEM_THREAD_SAFE=ON.")
-  endif()
-  find_package(OpenMP REQUIRED)
-  set(OPENMP_LIBRARIES ${OpenMP_CXX_LIBRARIES})
-  if(APPLE)
-    # On macOS, the compiler needs additional help to find the <omp.h> header.
-    # See issue #2642 for more information.
-    set(OPENMP_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS})
-  endif(APPLE)
-endif()
-
 # SuiteSparse (before SUNDIALS which may depend on KLU)
 if (MFEM_USE_SUITESPARSE)
  find_package(SuiteSparse REQUIRED
@@ -296,11 +320,6 @@ if (MFEM_USE_SUNDIALS)
  find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS})
 endif()

-# Mesquite
-if (MFEM_USE_MESQUITE)
-  find_package(Mesquite REQUIRED)
-endif()
-
 # SuperLU_DIST can only be enabled in parallel
 if (MFEM_USE_SUPERLU)
  if (MFEM_USE_MPI)
@@ -488,6 +507,10 @@ endif()
 #   an ALIAS target is missing?
 # Call Stack (most recent call first):
 #   CMakeLists.txt:474 (mfem_add_library)
+#
+# NOTE: We need to figure out which TPL library adds the dependency on
+#       "Threads::Threads" and call the next line only when that TPL library is
+#       enabled.  -V. Dobrev
 find_package(Threads REQUIRED)

 # List all possible libraries in order of dependencies.
@@ -495,8 +518,8 @@ find_package(Threads REQUIRED)
 #    With newer versions of SuiteSparse which include METIS header using 64-bit
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
-set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist METIS SuiteSparse SUNDIALS
-    PETSC SLEPC MESQUITE MUMPS STRUMPACK AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
+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
    ADIOS2 CUBLAS CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
    MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)
@@ -511,18 +534,13 @@ foreach(TPL IN LISTS MFEM_TPLS)
    list(APPEND TPL_INCLUDE_DIRS ${${TPL}_INCLUDE_DIRS})
  endif()
 endforeach(TPL)
+list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
+list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_INCLUDE_DIRS)
 # message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")

-if (OPENMP_FOUND)
-  message(STATUS "MFEM: using package OpenMP")
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-  if (MFEM_USE_CUDA)
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -Xcompiler=${OpenMP_CXX_FLAGS}")
-  endif()
-endif()
-
+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}")
@@ -578,7 +596,19 @@ 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}")

@@ -590,6 +620,7 @@ 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")
@@ -701,13 +732,15 @@ add_subdirectory(doc)
 message(STATUS "CMAKE_INSTALL_PREFIX = ${CMAKE_INSTALL_PREFIX}")
 set(INSTALL_INCLUDE_DIR include
  CACHE PATH "Relative path for installing header files.")
+set(INSTALL_BIN_DIR bin
+  CACHE PATH "Relative path for installing the binaries.")
 set(INSTALL_LIB_DIR lib
  CACHE PATH "Relative path for installing the library.")
 # other options: "share/mfem/cmake", "lib/mfem/cmake"
 set(INSTALL_CMAKE_DIR lib/cmake/mfem
  CACHE PATH "Relative path for installing cmake config files.")

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

@@ -719,7 +752,9 @@ set(CMAKE_INSTALL_DEFAULT_COMPONENT_NAME Development)
 # Install the library
 install(TARGETS ${PROJECT_NAME}
  EXPORT ${PROJECT_NAME_UC}Targets
-  DESTINATION ${INSTALL_LIB_DIR})
+  RUNTIME DESTINATION ${INSTALL_BIN_DIR}
+  LIBRARY DESTINATION ${INSTALL_LIB_DIR}
+  ARCHIVE DESTINATION ${INSTALL_LIB_DIR})

 # Install the master headers
 foreach(Header mfem.hpp mfem-performance.hpp)
@@ -757,10 +792,8 @@ endif()

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

 # Package the whole thing up nicely
@@ -337,10 +337,6 @@ MFEM_USE_SUNDIALS = YES/NO
   library. When enabled, this option uses the SUNDIALS_* library options,
   see below.

-MFEM_USE_MESQUITE = YES/NO
-   Enable MFEM functionality based on the Mesquite library. When enabled, this
-   option uses the MESQUITE_* library options, see below.
-
 MFEM_USE_SUITESPARSE = YES/NO
   Enable MFEM functionality based on the SuiteSparse library. Currently, this
   option adds the classes UMFPackSolver and KLUSolver (both sparse serial
@@ -636,11 +632,6 @@ The specific libraries and their options are:
  Options: SUNDIALS_OPT, SUNDIALS_LIB.
  Versions: SUNDIALS >= 5.0.0, SUNDIALS >= 5.4.0 for CUDA support.

- Mesquite (optional), used when MFEM_USE_MESQUITE = YES.
-  URL: http://trilinos.org/oldsite/packages/mesquite
-  Options: MESQUITE_OPT, MESQUITE_LIB.
-  The Mesquite support is deprecated and will be removed in the future.
-
 - SuiteSparse (optional), used when MFEM_USE_SUITESPARSE = YES.
  URL: http://faculty.cse.tamu.edu/davis/suitesparse.html
  Options: SUITESPARSE_OPT, SUITESPARSE_LIB.
@@ -802,10 +793,10 @@ The specific libraries and their options are:
  Versions: libCEED >= 0.10.

 - RAJA (optional), used when MFEM_USE_RAJA = YES.
-  Beginning with MFEM v4.3, only RAJA v0.14.0+ is supported.
+  Beginning with MFEM v4.5.1, only RAJA v2022.10.3+ is supported.
  URL: https://github.com/LLNL/RAJA
  Options: RAJA_DIR, RAJA_OPT, RAJA_LIB.
-  Versions: RAJA >= 0.14.0.
+  Versions: RAJA >= 2022.10.3.

 - Moonolith (optional), use when MFEM_USE_MOONOLITH = YES.
  URL: https://bitbucket.org/zulianp/par_moonolith
@@ -969,7 +960,6 @@ MFEM_USE_LEGACY_OPENMP
 MFEM_USE_OPENMP
 MFEM_USE_MEMALLOC
 MFEM_TIMER_TYPE - Set automatically, can be overwritten.
-MFEM_USE_MESQUITE
 MFEM_USE_SUITESPARSE
 MFEM_USE_SUPERLU
 MFEM_USE_MUMPS
@@ -1034,7 +1024,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
 - HYPRE
 - METIS - The option MFEM_USE_METIS_5 is auto-detected.
 - ParMETIS
- - MESQUITE
 - SuiteSparse
 - SuperLUDist, STRUMPACK
 - Ginkgo
@@ -1,6 +1,6 @@
 BSD 3-Clause License

-Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC
+Copyright (c) 2010-2023, 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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -212,10 +212,6 @@ IF (DEFINED TPL_ENABLE_SUNDIALS)
  SET(MFEM_USE_SUNDIALS ${TPL_ENABLE_SUNDIALS} CACHE BOOL "Enable SUNDIALS usage" FORCE)
 ENDIF()

-IF (DEFINED TPL_ENABLE_MESQUITE)
-  SET(MFEM_USE_MESQUITE ${TPL_ENABLE_MESQUITE} CACHE BOOL "Enable MESQUITE usage" FORCE)
-ENDIF()
-
 IF (DEFINED TPL_ENABLE_SUITESPARSE)
  SET(MFEM_USE_SUITESPARSE ${TPL_ENABLE_SUITESPARSE} CACHE BOOL "Enable SuiteSparse usage" FORCE)
 ENDIF()
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -29,7 +29,6 @@ set(MFEM_USE_LEGACY_OPENMP @MFEM_USE_LEGACY_OPENMP@)
 set(MFEM_USE_MEMALLOC @MFEM_USE_MEMALLOC@)
 set(MFEM_TIMER_TYPE @MFEM_TIMER_TYPE@)
 set(MFEM_USE_SUNDIALS @MFEM_USE_SUNDIALS@)
-set(MFEM_USE_MESQUITE @MFEM_USE_MESQUITE@)
 set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)
 set(MFEM_USE_MUMPS @MFEM_USE_MUMPS@)
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -39,6 +39,9 @@
 // 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
@@ -77,9 +80,6 @@
 // Internal MFEM option: enable group/batch allocation for some small objects.
 #cmakedefine MFEM_USE_MEMALLOC

-// Enable MFEM functionality based on the Mesquite library.
-#cmakedefine MFEM_USE_MESQUITE
-
 // Enable MFEM functionality based on the SuiteSparse library.
 #cmakedefine MFEM_USE_SUITESPARSE

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -16,7 +16,7 @@

 include(MfemCmakeUtilities)
 mfem_find_package(Algoim ALGOIM ALGOIM_DIR
-       	"include" "algoim_quad.hpp"
+        "include;src" "algoim_quad.hpp"
       	"" ""
-	"Paths to headers required by Algoim."
-	"Libraries required by Algoim.")
+        "Paths to headers required by Algoim."
+        "Libraries required by Algoim.")
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -16,7 +16,22 @@

 include(MfemCmakeUtilities)
 mfem_find_package(Caliper CALIPER CALIPER_DIR
-       	"include" "caliper/cali.h"
-       	"lib" "caliper"
-  	"Paths to headers required by Caliper."
-       	"Libraries required by Caliper.")
+      "include" "caliper/cali.h"
+      "lib" "caliper"
+      "Paths to headers required by Caliper."
+      "Libraries required by Caliper.")
+
+# Append adiak path/lib if the user provided ADIAK_DIR
+if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
+    find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})
+    list(APPEND CALIPER_INCLUDE_DIRS ${adiak_INCLUDE_DIRS})
+    list(APPEND CALIPER_LIBRARIES ${adiak_LIBRARIES})
+endif()
+
+# Append gotcha path/lib if the user provided GOTCHA_DIR
+if(GOTCHA_DIR AND EXISTS ${GOTCHA_DIR})
+    find_package(gotcha NO_DEFAULT_PATH REQUIRED PATHS ${GOTCHA_DIR}/lib/cmake/gotcha ${GOTCHA_DIR})
+    list(APPEND CALIPER_INCLUDE_DIRS ${gotcha_INCLUDE_DIRS})
+    list(APPEND CALIPER_LIBRARIES ${gotcha_LIBRARIES})
+endif()
+
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -21,8 +21,15 @@ 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.")
-  else()
-    message(STATUS "Conduit Relay HDF5 Support is DISABLED")
+  # 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")
 endif()

 include(MfemCmakeUtilities)
@@ -32,4 +39,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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -14,9 +14,21 @@
 #   - HDF5_LIBRARIES - The HDF5 libraries
 #   - HDF5_INCLUDE_DIRS - The HDF5 include directories

+# NOTE: Using this FindHDF5.cmake instead of the CMake provided version may lead
+#       to issues with some TPL libraries that depend (or may depend) on HDF5.
+#       For this reason, we should consider removing this file, or at least
+#       making it use the CMake provided version by default and apply the logic
+#       below only when specifically requested by a user.  -V. Dobrev
+
 # First Check for HDF5_DIR
 if(NOT HDF5_DIR)
-    MESSAGE(FATAL_ERROR "Could not find HDF5. HDF5 support needs explicit HDF5_DIR")
+  message(FATAL_ERROR
+    "Could not find HDF5. HDF5 support needs explicit HDF5_DIR")
+endif()
+
+if (NOT HDF5_FIND_QUIETLY)
+  message(STATUS "Looking for HDF5 ...")
+  message(STATUS "   in HDF5_DIR = ${HDF5_DIR}")
 endif()

 # Find includes
@@ -50,8 +62,9 @@ include(FindPackageHandleStandardArgs)

 # Handle the QUIETLY and REQUIRED arguments and set HDF5_FOUND to TRUE if all
 # listed variables are TRUE
-find_package_handle_standard_args(HDF5  DEFAULT_MSG
+find_package_handle_standard_args(HDF5
+                                  " *** HDF5 not found. Please set HDF5_DIR."
+                                  HDF5_LIBRARIES
                                  HDF5_INCLUDE_DIRS
                                  __HDF5_LIBRARY
-                                  __HDF5_HL_LIBRARY
-                                  HDF5_LIBRARIES )
+                                  __HDF5_HL_LIBRARY)
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,20 +0,0 @@
-# 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.
-
-# Defines the following variables:
-#   - MESQUITE_FOUND
-#   - MESQUITE_LIBRARIES
-#   - MESQUITE_INCLUDE_DIRS
-
-include(MfemCmakeUtilities)
-mfem_find_package(Mesquite MESQUITE MESQUITE_DIR
-  "include" "Mesquite_all_headers.hpp" "lib" "mesquite"
-  "Paths to headers required by Mesquite." "Libraries required by Mesquite.")
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -14,21 +14,38 @@
 #   - 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 "")
+
 # 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
 # as already having been found), so we reset the variable to force FindHDF5.cmake
 # to be called for a second time
 set(HDF5_FOUND OFF)
 enable_language(C) # FindHDF5.cmake uses the C compiler
+
 mfem_find_package(NetCDF NETCDF NETCDF_DIR "include" netcdf.h "lib" netcdf
  "Paths to headers required by NetCDF." "Libraries required by NetCDF.")
+
 # The HL extension libraries are in a separate variable and must precede
 # the "regular" hdf5 library, as hdf5_hl depends on hdf5
 # 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
-list(INSERT NETCDF_LIBRARIES 1 ${HDF5_C_LIBRARY_hdf5_hl})
+# (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()
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -14,17 +14,27 @@
 #   - RAJA_LIBRARIES
 #   - RAJA_INCLUDE_DIRS

-include(MfemCmakeUtilities)
-mfem_find_package(RAJA RAJA RAJA_DIR "include" "RAJA/RAJA.hpp" "lib" "RAJA"
-  "Paths to headers required by RAJA." "Libraries required by RAJA.")
-
-if (NOT RAJA_CONFIG_CMAKE)
-   set(RAJA_CONFIG_CMAKE "${RAJA_DIR}/share/raja/cmake/raja-config.cmake")
+if (RAJA_FOUND)
+   return()
 endif()
-if (EXISTS "${RAJA_CONFIG_CMAKE}")
-   include("${RAJA_CONFIG_CMAKE}")
-   if (ENABLE_CUDA AND NOT MFEM_USE_CUDA)
-      message(FATAL_ERROR
-         "RAJA is built with CUDA: MFEM_USE_CUDA=YES is required")
+message(STATUS "Looking for RAJA ...")
+if (RAJA_DIR)
+   message(STATUS "   in RAJA_DIR = ${RAJA_DIR}")
+   find_package(RAJA CONFIG NO_DEFAULT_PATH PATHS "${RAJA_DIR}")
+endif()
+if (NOT RAJA_FOUND)
+   message(STATUS "   in standard CMake locations")
+   find_package(RAJA CONFIG)
+endif()
+if (RAJA_FOUND)
+   set(RAJA_LIBRARIES "RAJA" CACHE STRING "RAJA imported target." FORCE)
+   set(RAJA_INCLUDE_DIRS "" CACHE STRING "RAJA include dirs (not used)" FORCE)
+   message(STATUS
+      "Found RAJA target: ${RAJA_LIBRARIES} (version: ${RAJA_VERSION})")
+else()
+   set(msg STATUS)
+   if (RAJA_FIND_REQUIRED)
+      set(msg FATAL_ERROR)
   endif()
+   message(${msg} "RAJA not found. Please set RAJA_DIR to the RAJA prefix.")
 endif()
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -14,6 +14,23 @@
 #   - UMPIRE_LIBRARIES
 #   - UMPIRE_INCLUDE_DIRS

-find_package(umpire REQUIRED CONFIG)
+if (NOT umpire_DIR AND UMPIRE_DIR)
+  set(umpire_DIR ${UMPIRE_DIR}/lib/cmake/umpire)
+endif()
+message(STATUS "Looking for UMPIRE ...")
+message(STATUS "   in UMPIRE_DIR = ${UMPIRE_DIR}")
+message(STATUS "      umpire_DIR = ${umpire_DIR}")
+find_package(umpire CONFIG)
 set(UMPIRE_FOUND ${umpire_FOUND})
 set(UMPIRE_LIBRARIES "umpire")
+if (UMPIRE_FOUND)
+  message(STATUS
+    "Found UMPIRE target: ${UMPIRE_LIBRARIES} (version: ${umpire_VERSION})")
+else()
+  set(msg STATUS)
+  if (UMPIRE_FIND_REQUIRED)
+    set(msg FATAL_ERROR)
+  endif()
+  message(${msg}
+    "UMPIRE not found. Please set UMPIRE_DIR to the install prefix.")
+endif()
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -46,6 +46,10 @@ endfunction()
 # Wrapper for add_executable
 macro(mfem_add_executable NAME)
  add_executable(${NAME} ${ARGN})
+  if (MFEM_USE_CUDA)
+    set_target_properties(${NAME} PROPERTIES
+      CUDA_RESOLVE_DEVICE_SYMBOLS ON)
+  endif()
 endmacro()

 # Wrapper for add_library
@@ -99,25 +103,6 @@ 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()

@@ -158,53 +143,16 @@ macro(add_mfem_miniapp MFEM_EXE_NAME)

  # Append the additional libraries and options
  if (LIBRARIES_LIST)
-    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
-      target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${LIBRARIES_LIST})
-    else()
-      target_link_libraries(${MFEM_EXE_NAME} ${LIBRARIES_LIST})
-    endif()
+    target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${LIBRARIES_LIST})
  endif()
  if (EXTRA_OPTIONS_LIST)
    string(REPLACE ";" " " EXTRA_OPTIONS_STRING "${EXTRA_OPTIONS_LIST}")
    message(STATUS "${MFEM_EXE_NAME}: add flags \"${EXTRA_OPTIONS_STRING}\"")
-    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
-      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${EXTRA_OPTIONS_LIST})
-    else()
-      get_target_property(THIS_COMPILE_FLAGS ${MFEM_EXE_NAME} COMPILE_FLAGS)
-      if (THIS_COMPILE_FLAGS)
-        set(THIS_COMPILE_FLAGS "${THIS_COMPILE_FLAGS} ${EXTRA_OPTIONS_STRING}")
-      else()
-        set(THIS_COMPILE_FLAGS "${EXTRA_OPTIONS_STRING}")
-      endif()
-      set_target_properties(${MFEM_EXE_NAME}
-        PROPERTIES COMPILE_FLAGS ${THIS_COMPILE_FLAGS})
-    endif()
+    target_compile_options(${MFEM_EXE_NAME} PRIVATE ${EXTRA_OPTIONS_LIST})
  endif()
  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.
-    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
-      target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_LIBRARIES})
-    else()
-      target_link_libraries(${MFEM_EXE_NAME} ${MPI_CXX_LIBRARIES})
-    endif()
-
-    if (MPI_CXX_INCLUDE_PATH)
-      target_include_directories(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_INCLUDE_PATH})
-    endif()
-    if (MPI_CXX_COMPILE_FLAGS)
-      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_COMPILE_FLAGS})
-    endif()
-
-    if (MPI_CXX_LINK_FLAGS)
-      set_target_properties(${MFEM_EXE_NAME} PROPERTIES
-        LINK_FLAGS "${MPI_CXX_LINK_FLAGS}")
-    endif()
-  endif()
 endmacro()


@@ -449,9 +397,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()
@@ -517,14 +465,18 @@ 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}")
-            get_target_property(ImpConfigs ${TargetName} IMPORTED_CONFIGURATIONS)
-            list(FIND ImpConfigs ${ImportConfig} _Index)
-            if ((_Index EQUAL -1) OR ("${ImportConfig}" STREQUAL "NO_CONFIG"))
+            if ((NOT ImportConfig) OR ("${ImportConfig}" STREQUAL "NO_CONFIG"))
              set(ImportConfig "NO_CONFIG")
              set(ImportConfigSuffix "")
              # message(FATAL_ERROR " *** ${ReqPack}: configuration "
@@ -542,7 +494,8 @@ 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})
+                IMPORTED_LINK_INTERFACE_LIBRARIES${ImportConfigSuffix}
+                INTERFACE_LINK_LIBRARIES)
              get_target_property(_value ${TargetName} ${_prop})
              if (_value)
                list(APPEND _Pack_LIBS ${_value})
@@ -793,6 +746,7 @@ 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)
@@ -814,18 +768,7 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
      message(STATUS " *** Warning: [${tgt}] LOCATION not defined!")
    endif()
  elseif ("${type}" STREQUAL "INTERFACE_LIBRARY")
-    get_target_property(Libs ${tgt} INTERFACE_LINK_LIBRARIES)
-    if (Libs)
-      foreach(Lib ${Libs})
-        if (NOT (TARGET ${Lib}))
-          list(APPEND LinkOpts "${Lib}")
-        else()
-          mfem_get_target_options(${Lib} COpts LOpts)
-          list(APPEND CompileOpts ${COpts})
-          list(APPEND LinkOpts ${LOpts})
-        endif()
-      endforeach()
-    endif()
+    # The INTERFACE_LINK_LIBRARIES property is handled below.
    # Other properties we may need to handle:
    # INTERFACE_LINK_DEPENDS
    # INTERFACE_LINK_DIRECTORIES
@@ -833,6 +776,26 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
  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()
+        else()
+          mfem_get_target_options(${Lib} COpts LOpts)
+          list(APPEND CompileOpts ${COpts})
+          list(APPEND LinkOpts ${LOpts})
+        endif()
+      endforeach()
+    endif()
  endif()

  # Other potentially relevant properties:
@@ -856,6 +819,9 @@ 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)


@@ -878,7 +844,7 @@ function(mfem_export_mk_files)
  set(CONFIG_MK_BOOL_VARS MFEM_USE_MPI MFEM_USE_METIS MFEM_USE_METIS_5
      MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB MFEM_USE_LIBUNWIND
      MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_LEGACY_OPENMP MFEM_USE_OPENMP
-      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
+      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_SUITESPARSE
      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
@@ -907,7 +873,7 @@ function(mfem_export_mk_files)
         "${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
         MFEM_CXXFLAGS)
  set(MFEM_TPLFLAGS "")
-  foreach(dir ${MFEM_TPL_INCLUDE_DIRS})
+  foreach(dir ${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-2022, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -30,11 +30,30 @@
 #endif

 // Windows specific options
-#ifdef _WIN32
-#ifndef _USE_MATH_DEFINES
+#if defined(_WIN32) && !defined(_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-2022, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -39,6 +39,9 @@
 // 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
@@ -85,9 +88,6 @@
 // Enable MFEM functionality based on the SUNDIALS libraries.
 // #define MFEM_USE_SUNDIALS

-// Enable MFEM functionality based on the Mesquite library.
-// #define MFEM_USE_MESQUITE
-
 // Enable MFEM functionality based on the SuiteSparse library.
 // #define MFEM_USE_SUITESPARSE

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -29,7 +29,6 @@ MFEM_USE_OPENMP        = @MFEM_USE_OPENMP@
 MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
-MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -19,6 +19,7 @@ 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)
@@ -30,7 +31,6 @@ option(MFEM_USE_OPENMP "Enable the OpenMP backend" OFF)
 option(MFEM_USE_LEGACY_OPENMP "Enable legacy OpenMP usage" OFF)
 option(MFEM_USE_MEMALLOC "Enable the internal MEMALLOC option." ON)
 option(MFEM_USE_SUNDIALS "Enable SUNDIALS usage" OFF)
-option(MFEM_USE_MESQUITE "Enable MESQUITE usage" OFF)
 option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" OFF)
 option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
 option(MFEM_USE_SUPERLU5 "Use the old SuperLU_DIST 5.1 version" OFF)
@@ -124,9 +124,6 @@ set(SUNDIALS_DIR "${MFEM_DIR}/../sundials-5.0.0/instdir" CACHE PATH
 # set(SUNDIALS_REQUIRED_PACKAGES "SuiteSparse/KLU/AMD/BTF/COLAMD/config"
 #     CACHE STRING "Additional packages required by SUNDIALS.")

-set(MESQUITE_DIR "${MFEM_DIR}/../mesquite-2.99" CACHE PATH
-    "Path to the Mesquite library.")
-
 set(SuiteSparse_DIR "${MFEM_DIR}/../SuiteSparse" CACHE PATH
    "Path to the SuiteSparse library.")
 set(SuiteSparse_REQUIRED_PACKAGES "BLAS" "METIS"
@@ -188,6 +185,7 @@ set(GNUTLS_DIR "" CACHE PATH "Path to the GnuTLS library.")

 set(GSLIB_DIR "" CACHE PATH "Path to the GSLIB library.")

+set(HDF5_DIR "/usr" CACHE PATH "Path to the HDF5 library.")
 set(NETCDF_DIR "" CACHE PATH "Path to the NetCDF library.")
 set(NetCDF_REQUIRED_PACKAGES "HDF5/C/HL" CACHE STRING
    "Additional packages required by NetCDF.")
@@ -230,17 +228,12 @@ set(MKL_LIBRARY_DIR "" CACHE STRING "Custom library subdirectory")

 set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
 set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
-# If RAJA is built with external CAMP:
-# set(RAJA_REQUIRED_PACKAGES "camp"
-#     CACHE STRING "Packages that RAJA depends on.")
-# set(camp_DIR "${MFEM_DIR}/../camp/lib/cmake/camp"
-#     CACHE PATH "Path to CAMP CMake files.")
 set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
 set(UMPIRE_DIR "${MFEM_DIR}/../umpire" CACHE PATH "Path to Umpire")
 set(CALIPER_DIR "${MFEM_DIR}/../caliper" CACHE PATH "Path to Caliper")
 set(BLITZ_DIR "${MFEM_DIR}/../blitz" CACHE PATH "Path to Blitz")
 set(ALGOIM_DIR "${MFEM_DIR}/../algoim" CACHE PATH "Path to Algoim")
-set(ALGOIM_REQUIRED_PACKAGES "BLITZ" CACHE STRING
+set(Algoim_REQUIRED_PACKAGES "Blitz" CACHE STRING
    "Packages that ALGOIM depends on.")

 set(BENCHMARK_DIR "${MFEM_DIR}/../google-benchmark" CACHE PATH
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -58,7 +58,7 @@ HIP_CXX = hipcc
 # example: gfx600 (tahiti), gfx700 (kaveri), gfx701 (hawaii), gfx801 (carrizo),
 # gfx900, gfx1010, etc.
 HIP_ARCH = gfx900
-HIP_FLAGS = --amdgpu-target=$(HIP_ARCH)
+HIP_FLAGS = --offload-arch=$(HIP_ARCH)
 HIP_XCOMPILER =
 HIP_XLINKER   = -Wl,

@@ -131,7 +131,6 @@ MFEM_USE_LEGACY_OPENMP = NO
 MFEM_USE_MEMALLOC      = YES
 MFEM_TIMER_TYPE        = $(if $(NOTMAC),2,4)
 MFEM_USE_SUNDIALS      = NO
-MFEM_USE_MESQUITE      = NO
 MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
@@ -270,11 +269,6 @@ endif
 # If SUNDIALS was built with KLU:
 # MFEM_USE_SUITESPARSE = YES

-# MESQUITE library configuration
-MESQUITE_DIR = @MFEM_DIR@/../mesquite-2.99
-MESQUITE_OPT = -I$(MESQUITE_DIR)/include
-MESQUITE_LIB = -L$(MESQUITE_DIR)/lib -lmesquite
-
 # SuiteSparse library configuration
 LIB_RT = $(if $(NOTMAC),-lrt,)
 SUITESPARSE_DIR = @MFEM_DIR@/../SuiteSparse
@@ -324,6 +318,9 @@ MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps\

 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build
+ifeq ($(MFEM_USE_STRUMPACK),YES)
+   BASE_FLAGS = -std=c++14
+endif
 STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)
 # If STRUMPACK was build with OpenMP support, the following may be need:
 # STRUMPACK_OPT += $(OPENMP_OPT)
@@ -476,7 +473,17 @@ OCCA_LIB = $(XLINKER)-rpath,$(OCCA_DIR)/lib -L$(OCCA_DIR)/lib -locca
 # CALIPER library configuration
 CALIPER_DIR = @MFEM_DIR@/../caliper
 CALIPER_OPT = -I$(CALIPER_DIR)/include
-CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib64 -lcaliper
+CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 $(XLINKER)-rpath,$(CALIPER_DIR)/lib -L$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib -lcaliper
+
+ifdef ADIAK_DIR
+   CALIPER_OPT += -I$(ADIAK_DIR)/include
+   CALIPER_LIB += $(XLINKER)-rpath,$(ADIAK_DIR)/lib64 $(XLINKER)-rpath,$(ADIAK_DIR)/lib -L$(ADIAK_DIR)/lib64 -L$(ADIAK_DIR)/lib -ladiak
+endif
+ifdef GOTCHA_DIR
+   CALIPER_OPT += -I$(GOTCHA_DIR)/include
+   CALIPER_LIB += $(XLINKER)-rpath,$(GOTCHA_DIR)/lib64 $(XLINKER)-rpath,$(GOTCHA_DIR)/lib -L$(GOTCHA_DIR)/lib64 -L$(GOTCHA_DIR)/lib -lgotcha
+endif
+

 # BLITZ library configuration
 BLITZ_DIR = @MFEM_DIR@/../blitz
@@ -499,20 +506,29 @@ CEED_OPT = -I$(CEED_DIR)/include
 CEED_LIB = $(XLINKER)-rpath,$(CEED_DIR)/lib -L$(CEED_DIR)/lib -lceed

 # RAJA library configuration
+ifeq ($(MFEM_USE_RAJA),YES)
+   BASE_FLAGS = -std=c++14
+endif
 RAJA_DIR = @MFEM_DIR@/../raja
 RAJA_OPT = -I$(RAJA_DIR)/include
 ifdef CUB_DIR
   RAJA_OPT += -I$(CUB_DIR)
 endif
+
+CAMP_LIB = -lcamp
 ifdef CAMP_DIR
   RAJA_OPT += -I$(CAMP_DIR)/include
+   CAMP_LIB = $(XLINKER)-rpath,$(CAMP_DIR)/lib -L$(CAMP_DIR)/lib -lcamp
 endif
-RAJA_LIB = $(XLINKER)-rpath,$(RAJA_DIR)/lib -L$(RAJA_DIR)/lib -lRAJA
+RAJA_LIB = $(XLINKER)-rpath,$(RAJA_DIR)/lib -L$(RAJA_DIR)/lib -lRAJA $(CAMP_LIB)

 # UMPIRE library configuration
+ifeq ($(MFEM_USE_UMPIRE),YES)
+   BASE_FLAGS = -std=c++14
+endif
 UMPIRE_DIR = @MFEM_DIR@/../umpire
 UMPIRE_OPT = -I$(UMPIRE_DIR)/include $(if $(CAMP_DIR), -I$(CAMP_DIR)/include)
-UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -lumpire
+UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -lumpire $(CAMP_LIB)

 # MKL CPardiso library configuration
 MKL_CPARDISO_DIR ?=
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -35,8 +35,8 @@ cd $(git rev-parse --show-toplevel)
 # copyright check
 copyright=true
 if [[ "${option}" == "--copyright" || "${option}" == "" ]]; then
-  if git grep -l "^\(#\|//\).*\(\-2021\|\ 2010,\)" > matches.txt; then
-    echo "Please update the following files to Copyright (c) 2010-2022:"
+  if git grep -n "^\(#\|//\).*Copyright.*2010-20\(2[^3]\|[^2].\)" > matches.txt; then
+    echo "Please update the following files to Copyright (c) 2010-2023:"
    cat matches.txt
    copyright=false
  fi
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,6 +1,6 @@
 #!/bin/bash

-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -138,7 +138,7 @@ groups_parallel=(
 '"meshing"
   "Meshing miniapps:"
   "miniapps/meshing"
-   "pmesh-optimizer.cpp pminimal-surface.cpp"'
+   "pmesh-optimizer.cpp pmesh-fitting.cpp pminimal-surface.cpp"'
 '"electromagnetics"
   "Electromagnetics miniapps:"
   "miniapps/electromagnetics"
@@ -227,7 +227,7 @@ groups_all=(
   "Meshing miniapps:"
   "miniapps/meshing"
   "mobius-strip.cpp klein-bottle.cpp extruder.cpp toroid.cpp
-    {,p}mesh-optimizer.cpp {,p}minimal-surface.cpp"'
+    {,p}mesh-optimizer.cpp pmesh-fitting.cpp {,p}minimal-surface.cpp"'
 '"electromagnetics"
   "Electromagnetics miniapps:"
   "miniapps/electromagnetics"
@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -0,0 +1,109 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+5
+1 3 4 5 6 7
+1 3 0 1 5 4
+1 3 1 2 6 5
+1 3 3 7 6 2
+1 3 0 4 7 3
+
+boundary
+4
+1 1 0 1
+1 1 2 3
+1 1 1 2
+1 1 3 0
+
+edges
+12
+0 0 1
+0 4 5
+0 7 6
+0 3 2
+1 1 2
+1 5 6
+1 4 7
+1 0 3
+2 0 4
+2 1 5
+2 2 6
+2 3 7
+
+vertices
+8
+
+knotvectors
+3
+2 3 0 0 0 1 1 1
+2 3 0 0 0 1 1 1
+2 3 0 0 0 1 1 1
+
+weights
+1
+1
+1
+1
+1
+1
+1
+1
+0.70710678118655
+1
+1
+0.70710678118655
+0.70710678118655
+1
+1
+0.70710678118655
+1
+1
+1
+1
+1
+0.85355339059327
+0.85355339059327
+0.85355339059327
+0.85355339059327
+
+FiniteElementSpace
+FiniteElementCollection: NURBS2
+VDim: 2
+Ordering: 1
+
+-0.70710678118 -0.70710678118
+0.70710678118 -0.70710678118
+0.70710678118 0.70710678118
+-0.70710678118 0.70710678118
+-0.35355339059 -0.35355339059
+0.35355339059 -0.35355339059
+0.35355339059 0.35355339059
+-0.35355339059 0.35355339059
+0 -1.41421356236
+0 -0.35355339059
+0 0.35355339059
+0 1.41421356236
+1.41421356236 0
+0.35355339059 0
+-0.35355339059 0
+-1.41421356236 0
+-0.530330085885 -0.530330085885
+0.530330085885 -0.530330085885
+0.530330085885 0.530330085885
+-0.530330085885 0.530330085885
+0 0
+0 -0.883883476475
+0.883883476475 0
+0 0.883883476475
+-0.883883476475 0
@@ -0,0 +1,46 @@
+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-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -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.0
+PROJECT_NUMBER         = v4.5.3

 # 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
@@ -763,36 +763,54 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                         @MFEM_SOURCE_DIR@/config \
                         @MFEM_SOURCE_DIR@/general \
                         @MFEM_SOURCE_DIR@/linalg \
+                         @MFEM_SOURCE_DIR@/linalg/simd \
                         @MFEM_SOURCE_DIR@/mesh \
+                         @MFEM_SOURCE_DIR@/mesh/submesh \
                         @MFEM_SOURCE_DIR@/fem \
-                         @MFEM_SOURCE_DIR@/fem/moonolith \
+                         @MFEM_SOURCE_DIR@/fem/ceed \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/convection \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/diffusion \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/mass \
+                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/nlconvection \
+                         @MFEM_SOURCE_DIR@/fem/ceed/interface \
+                         @MFEM_SOURCE_DIR@/fem/ceed/solvers \
                         @MFEM_SOURCE_DIR@/fem/fe \
                         @MFEM_SOURCE_DIR@/fem/lor \
+                         @MFEM_SOURCE_DIR@/fem/moonolith \
+                         @MFEM_SOURCE_DIR@/fem/qinterp \
+                         @MFEM_SOURCE_DIR@/fem/tmop \
                         @MFEM_SOURCE_DIR@/examples \
-                         @MFEM_SOURCE_DIR@/examples/caliper \
                         @MFEM_SOURCE_DIR@/examples/amgx \
+                         @MFEM_SOURCE_DIR@/examples/caliper \
                         @MFEM_SOURCE_DIR@/examples/ginkgo \
-                         @MFEM_SOURCE_DIR@/examples/moonolith \
                         @MFEM_SOURCE_DIR@/examples/hiop \
+                         @MFEM_SOURCE_DIR@/examples/moonolith \
                         @MFEM_SOURCE_DIR@/examples/petsc \
                         @MFEM_SOURCE_DIR@/examples/pumi \
                         @MFEM_SOURCE_DIR@/examples/sundials \
                         @MFEM_SOURCE_DIR@/examples/superlu \
                         @MFEM_SOURCE_DIR@/miniapps/adjoint \
+                         @MFEM_SOURCE_DIR@/miniapps/autodiff \
                         @MFEM_SOURCE_DIR@/miniapps/common \
                         @MFEM_SOURCE_DIR@/miniapps/electromagnetics \
                         @MFEM_SOURCE_DIR@/miniapps/gslib \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/kernels \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/materials \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/operators \
+                         @MFEM_SOURCE_DIR@/miniapps/hooke/preconditioners \
                         @MFEM_SOURCE_DIR@/miniapps/meshing \
                         @MFEM_SOURCE_DIR@/miniapps/mtop \
-                         @MFEM_SOURCE_DIR@/miniapps/autodiff \
+                         @MFEM_SOURCE_DIR@/miniapps/multidomain \
                         @MFEM_SOURCE_DIR@/miniapps/navier \
                         @MFEM_SOURCE_DIR@/miniapps/nurbs \
+                         @MFEM_SOURCE_DIR@/miniapps/parelag \
                         @MFEM_SOURCE_DIR@/miniapps/performance \
                         @MFEM_SOURCE_DIR@/miniapps/shifted \
                         @MFEM_SOURCE_DIR@/miniapps/solvers \
                         @MFEM_SOURCE_DIR@/miniapps/tools \
-                         @MFEM_SOURCE_DIR@/miniapps/toys \
-                         @MFEM_SOURCE_DIR@/miniapps/parelag
+                         @MFEM_SOURCE_DIR@/miniapps/toys

 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -33,13 +33,13 @@ add_mfem_examples(CALIPER_EXE_SRCS ${PREFIX})
 if (MFEM_ENABLE_TESTING)
  foreach(SRC_FILE ${CALIPER_EXE_SRCS})
    get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
-    string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
+    string(REPLACE ".cpp" "" TEST_NAME ${PREFIX}${SRC_FILENAME})

    set(THIS_TEST_OPTIONS "-no-vis")

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