Switching to single point Dirichlet BC

Add support for shared builds on Windows with MSVC using CMake

.appveyor.yml

@@ -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

.github/workflows/builds-and-tests.yml

@@ -74,6 +74,10 @@ 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
@@ -210,7 +214,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:
@@ -222,7 +226,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

.gitignore

@@ -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
@@ -131,12 +134,6 @@ examples/hiop/ex9-mesh.*
 examples/hiop/ex9-init.*
 examples/hiop/ex9-final.*
 
-examples/ipopt/exContactBlockTL
-examples/ipopt/exContactBlockTL.mesh
-examples/ipopt/exContactBlockTL-mesh.*
-examples/ipopt/exContactBlockTL-init.*
-examples/ipopt/exContactBlockTL-final.*
-
 examples/petsc/ex[1-69]p
 examples/petsc/ex1[0-1]p
 examples/petsc/mesh.*

CHANGELOG

@@ -47,6 +47,10 @@ New and updated examples and miniapps
 
 Integrations, testing and documentation
 ---------------------------------------
+- 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.
+
 - Removed the support for the Mesquite toolkit. We recommend using MFEM's TMOP
   functionality instead for mesh optimization. See the mesh-optimizer miniapp.
 

CMakeLists.txt

@@ -13,6 +13,7 @@
 # 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.")
 
@@ -24,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
@@ -183,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)
@@ -404,15 +425,6 @@ if (MFEM_USE_HIOP)
   # find_package updates HIOP_FOUND, HIOP_INCLUDE_DIRS, HIOP_LIBRARIES
 endif()
 
-# IpOpt optimizer
-if (MFEM_USE_IPOPT)
-  find_package(IPOPT REQUIRED)
-  message(
-      STATUS
-      "IPOPT_INCLUDE_DIRS=${IPOPT_INCLUDE_DIRS}, IPOPT_LIBRARIES=${IPOPT_LIBRARIES}, IPOPT_DIR=${IPOPT_DIR}")
-  # find_package updates IPOPT_FOUND, IPOPT_INCLUDE_DIRS, IPOPT_LIBRARIES
-endif()
-
 # CoDiPack package
 if (MFEM_USE_CODIPACK)
   find_package(CODIPACK REQUIRED)
@@ -508,7 +520,7 @@ find_package(Threads REQUIRED)
 #    be before SuiteSparse.
 set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
     SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
-    NETCDF MPFR PUMI HIOP IPOPT POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
+    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)
 
@@ -528,6 +540,7 @@ list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_INCLUDE_DIRS)
 # message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")
 
+message(STATUS "MFEM shared library: BUILD_SHARED_LIBS = ${BUILD_SHARED_LIBS}")
 message(STATUS "MFEM build type: CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 message(STATUS "MFEM version: v${MFEM_VERSION_STRING}")
 message(STATUS "MFEM git string: ${MFEM_GIT_STRING}")
@@ -595,6 +608,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")
@@ -766,10 +780,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

CONTRIBUTING.md

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

INSTALL

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

config/cmake/MFEMConfig.cmake.in

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

config/cmake/config.hpp.in

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

config/cmake/modules/FindIPOPT.cmake

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

config/cmake/modules/MfemCmakeUtilities.cmake

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

config/config.hpp

@@ -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.

config/config.hpp.in

@@ -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
@@ -141,9 +144,6 @@
 // Enable MFEM functionality based on the HIOP library.
 // #define MFEM_USE_HIOP
 
-// Enable MFEM functionality based on the IPOPT library.
-// #define MFEM_USE_IPOPT
-
 // Enable MFEM functionality based on the GSLIB library
 // #define MFEM_USE_GSLIB
 

config/config.mk.in

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

config/defaults.cmake

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

config/defaults.mk

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

config/sample-runs.sh

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

doc/CodeDocumentation.conf.in

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

examples/CMakeLists.txt

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

examples/contact.cpp

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

examples/ex25.cpp

@@ -469,7 +469,7 @@ int main(int argc, char *argv[])
 
       Operator *pc_r = nullptr;
       Operator *pc_i = nullptr;
-      int s = (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0;
+      double s = (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0;
       if (pa)
       {
          // Jacobi Smoother

examples/ex30.cpp

@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    int nc_limit = 1;
-   int max_elems = 1e5;
+   int max_elems = 100*1000;
    double double_max_elems = double(max_elems);
    bool visualization = true;
    double osc_threshold = 1e-3;

examples/ex33.cpp

@@ -123,7 +123,7 @@ int main(int argc, char *argv[])
 
    // 2. Compute the rational expansion coefficients that define the
    //    integer-order PDEs.
-   const int power_of_laplace = floor(alpha);
+   const int power_of_laplace = (int)floor(alpha);
    double exponent_to_approximate = alpha - power_of_laplace;
    bool integer_order = false;
    // Check if alpha is an integer or not.

examples/ipopt/CMakeLists.txt

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

examples/ipopt/README

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

examples/ipopt/block1.mesh

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

examples/ipopt/block2.mesh

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

examples/ipopt/contact.cpp

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

examples/ipopt/exContactBlockTL.hpp

@@ -1,230 +0,0 @@
-//                                Contact example
-//
-// Compile with: make exContactBlockTL
-//
-// Sample runs:  ./exContactBlockTL -m1 block1.mesh -m2 block2.mesh -at "5 6 7 8"
-// Sample runs:  ./exContactBlockTL -m1 block1_d.mesh -m2 block2_d.mesh -at "5 6 7 8"
-
-#ifndef EXCONTACTBLOCKTL_HPP
-#define EXCONTACTBLOCKTL_HPP
-
-#include "mfem.hpp"
-#include "IpTNLP.hpp"
-
-using namespace std;
-using namespace mfem;
-using namespace Ipopt;
-
-
-class ExContactBlockTL: public TNLP
-{
-public:
-   /** default constructor */
-   ExContactBlockTL(int argc, char *argv[]);
-
-   /** default destructor */
-   virtual ~ExContactBlockTL();
-
-   /**@name Overloaded from TNLP */
-   /** Method to return some info about the nlp */
-   virtual bool get_nlp_info(
-      Index&          n,
-      Index&          m,
-      Index&          nnz_jac_g,
-      Index&          nnz_h_lag,
-      IndexStyleEnum& index_style
-   );
-
-   /** Method to return the bounds for my problem */
-   virtual bool get_bounds_info(
-      Index   n,
-      Number* x_l,
-      Number* x_u,
-      Index   m,
-      Number* g_l,
-      Number* g_u
-   );
-
-   /** Method to return the starting point for the algorithm */
-   virtual bool get_starting_point(
-      Index   n,
-      bool    init_x,
-      Number* x,
-      bool    init_z,
-      Number* z_L,
-      Number* z_U,
-      Index   m,
-      bool    init_lambda,
-      Number* lambda
-   );
-
-   /** Method to return the objective value */
-   virtual bool eval_f(
-      Index         n,
-      const Number* x,
-      bool          new_x,
-      Number&       obj_value
-   );
-
-   /** Method to return the gradient of the objective */
-   virtual bool eval_grad_f(
-      Index         n,
-      const Number* x,
-      bool          new_x,
-      Number*       grad_f
-   );
-
-   /** Method to return the constraint residuals */
-   virtual bool eval_g(
-      Index         n,
-      const Number* x,
-      bool          new_x,
-      Index         m,
-      Number*       cons
-   );
-
-   /** Method to return:
-    *   1) The structure of the Jacobian (if "values" is NULL)
-    *   2) The values of the Jacobian (if "values" is not NULL)
-    */
-   virtual bool eval_jac_g(
-      Index         n,
-      const Number* x,
-      bool          new_x,
-      Index         m,
-      Index         nele_jac,
-      Index*        iRow,
-      Index*        jCol,
-      Number*       values
-   );
-
-   /** Method to return:
-    *   1) The structure of the Hessian of the Lagrangian (if "values" is NULL)
-    *   2) The values of the Hessian of the Lagrangian (if "values" is not NULL)
-    */
-   virtual bool eval_h(
-      Index         n,
-      const Number* x,
-      bool          new_x,
-      Number        obj_factor,
-      Index         m,
-      const Number* lambda,
-      bool          new_lambda,
-      Index         nele_hess,
-      Index*        iRow,
-      Index*        jCol,
-      Number*       values
-   );
-
-   /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
-   virtual void finalize_solution(
-      SolverReturn               status,
-      Index                      n,
-      const Number*              x,
-      const Number*              z_L,
-      const Number*              z_U,
-      Index                      m,
-      const Number*              g,
-      const Number*              lambda,
-      Number                     obj_value,
-      const IpoptData*           ip_data,
-      IpoptCalculatedQuantities* ip_cq
-   );
-
-private:
-   void update_g();
-   void update_jac();
-   void update_hess();
-
-private:
-   /**@name Methods to block default compiler methods.
-    *
-    * The compiler automatically generates the following three methods.
-    *  Since the default compiler implementation is generally not what
-    *  you want (for all but the most simple classes), we usually
-    *  put the declarations of these methods in the private section
-    *  and never implement them. This prevents the compiler from
-    *  implementing an incorrect "default" behavior without us
-    *  knowing. (See Scott Meyers book, "Effective C++")
-    */
-   ExContactBlockTL(
-      const ExContactBlockTL&
-   );
-
-   ExContactBlockTL& operator=(
-      const ExContactBlockTL&
-   );
-
-   Array<int> attr;
-   Array<int> m_attr;
-   Array<int> s_conn; // connectivity of the second/slave mesh
-   std::string mesh_file1;
-   std::string mesh_file2;
-   Mesh* mesh1;
-   Mesh* mesh2;
-   FiniteElementCollection* fec1;
-   FiniteElementCollection* fec2;
-   FiniteElementSpace* fespace1;
-   FiniteElementSpace* fespace2;
-   Array<int> ess_tdof_list1;
-   Array<int> ess_tdof_list2;
-   GridFunction nodes0;
-   GridFunction* nodes1;
-   GridFunction* nodes2;
-   GridFunction* x1;
-   GridFunction* x2;
-   LinearForm* b1;
-   LinearForm* b2;
-   PWConstCoefficient* lambda1_func;
-   PWConstCoefficient* lambda2_func;
-   PWConstCoefficient* mu1_func;
-   PWConstCoefficient* mu2_func;
-   BilinearForm* a1;
-   BilinearForm* a2;
-
-   mfem::Vector lambda1;
-   mfem::Vector lambda2;
-   mfem::Vector mu1;
-   mfem::Vector mu2;
-   mfem::Vector xyz;
-
-   std::set<int> bdryVerts2;
-
-   int dim;
-   // degrees of freedom of both meshes
-   int ndof_1;
-   int ndof_2;
-   int ndofs;
-   // number of nodes for each mesh
-   int nnd_1;
-   int nnd_2;
-   int nnd;
-
-   int npoints;
-
-   SparseMatrix A1;
-   mfem::Vector B1, X1;
-   SparseMatrix A2;
-   mfem::Vector B2, X2;
-
-   SparseMatrix* K;
-   mfem::Vector g;
-   mfem::Vector m_xi;
-   mfem::Vector xs;
-
-   Array<int> m_conn; // only works for linear elements that have 4 vertices!
-   DenseMatrix* coordsm;
-   SparseMatrix* M;
-
-   std::vector<SparseMatrix>* dM;
-
-   Array<int> Dirichlet_dof;
-   Array<double> Dirichlet_val;
-
-public: 
-   Mesh * GetMesh1() {return mesh1;}
-   Mesh * GetMesh2() {return mesh2;}
-
-};
-
-#endif

examples/ipopt/makefile

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

examples/ipopt/nodepair.hpp

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

examples/makefile

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

examples/nodepair.hpp

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

fem/ceed/interface/util.cpp

@@ -19,11 +19,11 @@
 
 #include <sys/types.h>
 #include <sys/stat.h>
-#ifndef _WIN32
+#if !defined(_WIN32) || !defined(_MSC_VER)
 typedef struct stat struct_stat;
 #else
 #define stat(dir, buf) _stat(dir, buf)
-#define S_ISDIR(mode) _S_IFDIR(mode)
+#define S_ISDIR(mode) (((mode) & _S_IFMT) == _S_IFDIR)
 typedef struct _stat struct_stat;
 #endif
 

fem/fe/fe_base.hpp

@@ -1154,7 +1154,7 @@ public:
    ~Poly_1D();
 };
 
-extern Poly_1D poly1d;
+extern MFEM_EXPORT Poly_1D poly1d;
 
 /// An element defined as an ND tensor product of 1D elements on a segment,
 /// square, or cube

fem/geom.hpp

@@ -29,7 +29,7 @@ namespace mfem
     Geometry::PRISM - w/ vert. (0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,0,1),(0,1,1)
     Geometry::PYRAMID - w/ vert. (0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1)
 */
-class Geometry
+class MFEM_EXPORT Geometry
 {
 public:
    enum Type
@@ -126,7 +126,11 @@ public:
    int NumBdr(int GeomType) { return NumBdrArray[GeomType]; }
 };
 
-template <> struct Geometry::Constants<Geometry::POINT>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::POINT>
 {
    static const int Dimension = 0;
    static const int NumVert = 1;
@@ -136,7 +140,11 @@ template <> struct Geometry::Constants<Geometry::POINT>
    static const int InvOrient[NumOrient];
 };
 
-template <> struct Geometry::Constants<Geometry::SEGMENT>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::SEGMENT>
 {
    static const int Dimension = 1;
    static const int NumVert = 2;
@@ -148,7 +156,11 @@ template <> struct Geometry::Constants<Geometry::SEGMENT>
    static const int InvOrient[NumOrient];
 };
 
-template <> struct Geometry::Constants<Geometry::TRIANGLE>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::TRIANGLE>
 {
    static const int Dimension = 2;
    static const int NumVert = 3;
@@ -174,7 +186,11 @@ template <> struct Geometry::Constants<Geometry::TRIANGLE>
    static const int InvOrient[NumOrient];
 };
 
-template <> struct Geometry::Constants<Geometry::SQUARE>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::SQUARE>
 {
    static const int Dimension = 2;
    static const int NumVert = 4;
@@ -194,7 +210,11 @@ template <> struct Geometry::Constants<Geometry::SQUARE>
    static const int InvOrient[NumOrient];
 };
 
-template <> struct Geometry::Constants<Geometry::TETRAHEDRON>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::TETRAHEDRON>
 {
    static const int Dimension = 3;
    static const int NumVert = 4;
@@ -216,7 +236,11 @@ template <> struct Geometry::Constants<Geometry::TETRAHEDRON>
    static const int InvOrient[NumOrient];
 };
 
-template <> struct Geometry::Constants<Geometry::CUBE>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::CUBE>
 {
    static const int Dimension = 3;
    static const int NumVert = 8;
@@ -234,7 +258,11 @@ template <> struct Geometry::Constants<Geometry::CUBE>
    };
 };
 
-template <> struct Geometry::Constants<Geometry::PRISM>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::PRISM>
 {
    static const int Dimension = 3;
    static const int NumVert = 6;
@@ -252,7 +280,11 @@ template <> struct Geometry::Constants<Geometry::PRISM>
    };
 };
 
-template <> struct Geometry::Constants<Geometry::PYRAMID>
+template <> struct
+/// @cond Suppress_Doxygen_warnings
+   MFEM_EXPORT
+/// @endcond
+   Geometry::Constants<Geometry::PYRAMID>
 {
    static const int Dimension = 3;
    static const int NumVert = 5;
@@ -272,7 +304,7 @@ template <> struct Geometry::Constants<Geometry::PYRAMID>
 
 // Defined in fe.cpp to ensure construction after 'mfem::TriangleFE' and
 // `mfem::TetrahedronFE`.
-extern Geometry Geometries;
+extern MFEM_EXPORT Geometry Geometries;
 
 
 class RefinedGeometry
@@ -321,7 +353,7 @@ public:
    ~GeometryRefiner();
 };
 
-extern GeometryRefiner GlobGeometryRefiner;
+extern MFEM_EXPORT GeometryRefiner GlobGeometryRefiner;
 
 }
 

fem/intrules.hpp

@@ -376,10 +376,10 @@ public:
 };
 
 /// A global object with all integration rules (defined in intrules.cpp)
-extern IntegrationRules IntRules;
+extern MFEM_EXPORT IntegrationRules IntRules;
 
 /// A global object with all refined integration rules
-extern IntegrationRules RefinedIntRules;
+extern MFEM_EXPORT IntegrationRules RefinedIntRules;
 
 }
 

general/device.hpp

@@ -125,7 +125,7 @@ private:
    enum MODES {SEQUENTIAL, ACCELERATED};
 
    static bool device_env, mem_host_env, mem_device_env, mem_types_set;
-   static Device device_singleton;
+   static MFEM_EXPORT Device device_singleton;
 
    MODES mode = Device::SEQUENTIAL;
    int dev = 0;   ///< Device ID of the configured device.

general/globals.hpp

@@ -63,12 +63,12 @@ public:
 /** @brief Global stream used by the library for standard output. Initially it
     uses the same std::streambuf as std::cout, however that can be changed.
     @sa OutStream. */
-extern OutStream out;
+extern MFEM_EXPORT OutStream out;
 /** @brief Global stream used by the library for standard error output.
     Initially it uses the same std::streambuf as std::cerr, however that can be
     changed.
     @sa OutStream. */
-extern OutStream err;
+extern MFEM_EXPORT OutStream err;
 
 
 /** @brief Construct a string of the form "<prefix><myid><suffix>" where the

general/isockstream.cpp

@@ -25,11 +25,11 @@
 #include <winsock.h>
 #ifdef _MSC_VER
 typedef int ssize_t;
+// Link with ws2_32.lib
+#pragma comment(lib, "ws2_32.lib")
 #endif
 typedef int socklen_t;
 #define close closesocket
-// Link with ws2_32.lib
-#pragma comment(lib, "ws2_32.lib")
 #endif
 
 using namespace std;

general/mem_manager.hpp

@@ -64,7 +64,7 @@ constexpr int DeviceMemoryType = static_cast<int>(MemoryType::MANAGED);
 constexpr int DeviceMemoryTypeSize = MemoryTypeSize - DeviceMemoryType;
 
 /// Memory type names, used during Device:: configuration.
-extern const char *MemoryTypeName[MemoryTypeSize];
+extern MFEM_EXPORT const char *MemoryTypeName[MemoryTypeSize];
 
 /// Memory classes identify sets of memory types.
 /** This type is used by kernels that can work with multiple MemoryType%s.
@@ -588,7 +588,7 @@ private:
 /** The MFEM memory manager class. Host-side pointers are inserted into this
     manager which keeps track of the associated device pointer, and where the
     data currently resides. */
-class MemoryManager
+class MFEM_EXPORT MemoryManager
 {
 private:
 
@@ -1281,7 +1281,7 @@ inline int Memory<T>::CompareHostAndDevice(int size) const
 
 
 /// The (single) global memory manager object
-extern MemoryManager mm;
+extern MFEM_EXPORT MemoryManager mm;
 
 } // namespace mfem
 

general/socketstream.cpp

@@ -30,10 +30,10 @@
 #include <winsock.h>
 #ifdef _MSC_VER
 typedef int ssize_t;
-#endif
 // Link with ws2_32.lib
 #pragma comment(lib, "ws2_32.lib")
 #endif
+#endif
 
 #ifdef MFEM_USE_GNUTLS
 #include <cstdlib>  // getenv

general/tic_toc.hpp

@@ -66,7 +66,7 @@ public:
 };
 
 
-extern StopWatch tic_toc;
+extern MFEM_EXPORT StopWatch tic_toc;
 
 /// Start the tic_toc timer
 extern void tic();

general/tinyxml2.cpp

@@ -21,6 +21,8 @@ must not be misrepresented as being the original software.
 distribution.
 */
 
+// This file was modified by the MFEM team.
+
 #include "tinyxml2.h"
 
 #include <new>		// yes, this one new style header, is in the Android SDK.
@@ -106,7 +108,7 @@ distribution.
 #elif defined(__APPLE__) || (__FreeBSD__)
 	#define TIXML_FSEEK fseeko
 	#define TIXML_FTELL ftello
-#elif defined(__unix__) && defined(__x86_64__)
+#elif defined(__unix__) && defined(__x86_64__) && !defined(__CYGWIN__)
 	#define TIXML_FSEEK fseeko64
 	#define TIXML_FTELL ftello64
 #else

general/version.cpp

@@ -97,9 +97,6 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_HIOP
       "MFEM_USE_HIOP\n"
 #endif
-#ifdef MFEM_USE_IPOPT
-      "MFEM_USE_IPOPT\n"
-#endif
 #ifdef MFEM_USE_HIP
       "MFEM_USE_HIP\n"
 #endif

linalg/ode.hpp

@@ -204,7 +204,7 @@ public:
 class RK6Solver : public ExplicitRKSolver
 {
 private:
-   static const double a[28], b[8], c[7];
+   static MFEM_EXPORT const double a[28], b[8], c[7];
 
 public:
    RK6Solver() : ExplicitRKSolver(8, a, b, c) { }
@@ -258,7 +258,7 @@ public:
 class AB1Solver : public AdamsBashforthSolver
 {
 private:
-   static const double a[1];
+   static MFEM_EXPORT const double a[1];
 
 public:
    AB1Solver() : AdamsBashforthSolver(1, a) { }
@@ -268,7 +268,7 @@ public:
 class AB2Solver : public AdamsBashforthSolver
 {
 private:
-   static const double a[2];
+   static MFEM_EXPORT const double a[2];
 
 public:
    AB2Solver() : AdamsBashforthSolver(2, a) { }
@@ -278,7 +278,7 @@ public:
 class AB3Solver : public AdamsBashforthSolver
 {
 private:
-   static const double a[3];
+   static MFEM_EXPORT const double a[3];
 
 public:
    AB3Solver() : AdamsBashforthSolver(3, a) { }
@@ -288,7 +288,7 @@ public:
 class AB4Solver : public AdamsBashforthSolver
 {
 private:
-   static const double a[4];
+   static MFEM_EXPORT const double a[4];
 
 public:
    AB4Solver() : AdamsBashforthSolver(4, a) { }
@@ -298,7 +298,7 @@ public:
 class AB5Solver : public AdamsBashforthSolver
 {
 private:
-   static const double a[5];
+   static MFEM_EXPORT const double a[5];
 
 public:
    AB5Solver() : AdamsBashforthSolver(5, a) { }
@@ -339,7 +339,7 @@ public:
 class AM0Solver : public AdamsMoultonSolver
 {
 private:
-   static const double a[1];
+   static MFEM_EXPORT const double a[1];
 
 public:
    AM0Solver() : AdamsMoultonSolver(0, a) { }
@@ -350,7 +350,7 @@ public:
 class AM1Solver : public AdamsMoultonSolver
 {
 private:
-   static const double a[2];
+   static MFEM_EXPORT const double a[2];
 
 public:
    AM1Solver() : AdamsMoultonSolver(1, a) { }
@@ -360,7 +360,7 @@ public:
 class AM2Solver : public AdamsMoultonSolver
 {
 private:
-   static const double a[3];
+   static MFEM_EXPORT const double a[3];
 
 public:
    AM2Solver() : AdamsMoultonSolver(2, a) { }
@@ -370,7 +370,7 @@ public:
 class AM3Solver : public AdamsMoultonSolver
 {
 private:
-   static const double a[4];
+   static MFEM_EXPORT const double a[4];
 
 public:
    AM3Solver() : AdamsMoultonSolver(3, a) { }
@@ -380,7 +380,7 @@ public:
 class AM4Solver : public AdamsMoultonSolver
 {
 private:
-   static const double a[5];
+   static MFEM_EXPORT const double a[5];
 
 public:
    AM4Solver() : AdamsMoultonSolver(4, a) { }

makefile

@@ -119,7 +119,7 @@ $(if $(word 2,$(SRC)),$(error Spaces in SRC = "$(SRC)" are not supported))
 MFEM_GIT_STRING = $(shell [ -d $(MFEM_DIR)/.git ] && git -C $(MFEM_DIR) \
    describe --all --long --abbrev=40 --dirty --always 2> /dev/null)
 
-EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop ipopt petsc pumi sundials superlu moonolith
+EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
@@ -198,6 +198,7 @@ ALL_LIBS =
 # Building static and/or shared libraries:
 MFEM_STATIC ?= $(STATIC)
 MFEM_SHARED ?= $(SHARED)
+MFEM_SHARED_BUILD = $(MFEM_SHARED)
 
 # Internal shortcuts
 override static = $(if $(MFEM_STATIC:YES=),,YES)
@@ -275,7 +276,7 @@ endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = ENZYME SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
- SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP IPOPT\
+ SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER PARELAG BENCHMARK\
  MOONOLITH ALGOIM
 
@@ -341,12 +342,12 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
- MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_IPOPT MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
+ MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX\
  MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
  MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_ALGOIM MFEM_USE_ENZYME\
- MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -589,7 +590,8 @@ install: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))
 	      > $(PREFIX_INC)/$$hdr && chmod 640 $(PREFIX_INC)/$$hdr; done
 # install config include
 	mkdir -p $(PREFIX_INC)/mfem/config
-	$(INSTALL) -m 640 $(BLD)config/_config.hpp $(PREFIX_INC)/mfem/config/config.hpp
+	$(INSTALL) -m 640 $(BLD)config/_config.hpp $(PREFIX_INC)/mfem/config/_config.hpp
+	$(INSTALL) -m 640 $(SRC)config/config.hpp $(PREFIX_INC)/mfem/config/config.hpp
 	$(INSTALL) -m 640 $(SRC)config/tconfig.hpp $(PREFIX_INC)/mfem/config
 # install remaining includes in each subdirectory
 	for dir in $(DIRS); do \
@@ -690,7 +692,6 @@ status info:
 	$(info MFEM_USE_CONDUIT       = $(MFEM_USE_CONDUIT))
 	$(info MFEM_USE_PUMI          = $(MFEM_USE_PUMI))
 	$(info MFEM_USE_HIOP          = $(MFEM_USE_HIOP))
-	$(info MFEM_USE_IPOPT          = $(MFEM_USE_IPOPT))
 	$(info MFEM_USE_GSLIB         = $(MFEM_USE_GSLIB))
 	$(info MFEM_USE_CUDA          = $(MFEM_USE_CUDA))
 	$(info MFEM_USE_HIP           = $(MFEM_USE_HIP))

mesh/hexahedron.hpp

@@ -68,7 +68,7 @@ public:
    virtual ~Hexahedron() { }
 };
 
-extern class TriLinear3DFiniteElement HexahedronFE;
+extern MFEM_EXPORT class TriLinear3DFiniteElement HexahedronFE;
 
 }
 

mesh/point.hpp

@@ -63,7 +63,7 @@ public:
 };
 
 class PointFiniteElement;
-extern PointFiniteElement PointFE;
+extern MFEM_EXPORT PointFiniteElement PointFE;
 
 }
 

mesh/quadrilateral.hpp

@@ -69,7 +69,7 @@ public:
    virtual ~Quadrilateral() { }
 };
 
-extern class BiLinear2DFiniteElement QuadrilateralFE;
+extern MFEM_EXPORT class BiLinear2DFiniteElement QuadrilateralFE;
 
 }
 

mesh/segment.hpp

@@ -69,7 +69,7 @@ public:
 };
 
 class Linear1DFiniteElement;
-extern Linear1DFiniteElement SegmentFE;
+extern MFEM_EXPORT Linear1DFiniteElement SegmentFE;
 
 }
 

mesh/tetrahedron.hpp

@@ -117,7 +117,7 @@ public:
 };
 
 // Defined in fe.cpp to ensure construction before 'mfem::Geometries'.
-extern class Linear3DFiniteElement TetrahedronFE;
+extern MFEM_EXPORT class Linear3DFiniteElement TetrahedronFE;
 
 }
 

mesh/triangle.hpp

@@ -98,7 +98,7 @@ public:
 };
 
 // Defined in fe.cpp to ensure construction before 'mfem::Geometries'.
-extern Linear2DFiniteElement TriangleFE;
+extern MFEM_EXPORT Linear2DFiniteElement TriangleFE;
 
 }
 

mesh/wedge.hpp

@@ -71,7 +71,7 @@ public:
    virtual ~Wedge() { }
 };
 
-extern class LinearWedgeFiniteElement WedgeFE;
+extern MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE;
 
 }
 

miniapps/common/pfem_extras.cpp

@@ -115,9 +115,15 @@ IrrotationalProjector
      ownsWeakDiv_(weakDiv == NULL),
      ownsGrad_(grad == NULL)
 {
+   int myid = H1FESpace_->GetMyRank();
+
    ess_bdr_.SetSize(H1FESpace_->GetParMesh()->bdr_attributes.Max());
-   ess_bdr_ = 1;
-   H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
+   ess_bdr_ = 0;
+   ess_bdr_tdofs_.SetSize((myid == 0) ? 1 : 0);
+   if (myid == 0)
+   {
+      ess_bdr_tdofs_[0] = 0;
+   }
 
    int geom = H1FESpace_->GetFE(0)->GetGeomType();
    const IntegrationRule * ir = &IntRules.Get(geom, irOrder);