diff --git a/integrated-tests/Test-scripts/test_fierro_cuda.py b/integrated-tests/Test-scripts/test_fierro_cuda.py
new file mode 100644
index 000000000..471ac648f
--- /dev/null
+++ b/integrated-tests/Test-scripts/test_fierro_cuda.py
@@ -0,0 +1,152 @@
+#!/usr/bin/python3
+
+import os
+import os.path
+import sys
+import math
+
+
+solvers = ["fierro-parallel-explicit"]
+
+executables = []
+tests = []
+
+# Add paths to all tested executables
+for i in range(len(solvers)):
+    executables.append("./../../build-fierro-cuda/bin/"+solvers[i])
+
+# Check that each executable exists
+for i in range(len(solvers)):
+    if not os.path.exists(executables[i]):
+        raise ValueError(solvers[i]+" executable not found in build-fierro-openmp directory")
+\
+# Add names of each test
+parallel_explicit_tests = ["Noh", "Sedov", "Sod"]
+# parallel_implicit_tests = ["Beam"]
+
+inputs = []
+standard_results = []
+
+tests.append(parallel_explicit_tests)
+# tests.append(parallel_implicit_tests)
+
+position_keyword = "POINTS"
+
+for i in range(len(solvers)):
+    tmp1 = []
+    tmp2 = []
+    for j in range(len(tests[i])):
+        tmp1.append("Solver-Inputs/"+solvers[i]+"/"+tests[i][j]+".yaml")
+        tmp2.append("standard-results/"+solvers[i]+"/"+tests[i][j]+"/vtk/data/VTK0.vtk")
+    inputs.append(tmp1)
+    standard_results.append(tmp2)
+
+# Extract vector valued data from vtk output file
+def extract_vector_data(filename, keyword):
+    data = []
+    found_keyword = False
+
+    with open(filename, 'r') as file:
+        for line in file:
+            if found_keyword:
+                if line.strip():  # If the line is not blank
+                    # Split each line into fields.
+                    fields = line.split()
+                    xx = float(fields[0])
+                    yy = float(fields[1])
+                    zz = float(fields[2])
+                    # Group and append to results.
+                    vec = [xx, yy, zz]
+                    data.append(vec)
+                else:
+                    break  # Exit the loop when encountering a blank line
+            elif keyword in line:
+                found_keyword = True
+
+    return data
+
+# Extract scalar valued data from vtk output file
+def extract_scalar_data(filename, keyword):
+    data = []
+    found_keyword = False
+
+    with open(filename, 'r') as file:
+        for line in file:
+            if found_keyword:
+                if line.strip():  # If the line is not blank
+                    fields = line.split()
+                    if fields[0] != "LOOKUP_TABLE":
+                        data.append(float(fields[0]))
+                else:
+                    break  # Exit the loop when encountering a blank line
+            elif keyword in line:
+                found_keyword = True
+                # file.next() # skip first line "LOOKUP_TABLE"
+ 
+    return data
+
+# Calculate the percent difference between two arrays of scalars
+def percent_difference_scalars(array1, array2):
+    if len(array1) != len(array2):
+        raise ValueError("Arrays must have the same length")
+
+    percent_diff = []
+    for i in range(len(array1)):
+        diff = array2[i] - array1[i]
+        percent_diff.append((diff / array1[i]) * 100 if array1[i] != 0 else 0)
+
+    return percent_diff
+
+# Calculate the percent difference between two arrays of vectors
+def percent_difference_vectors(array1, array2):
+    if len(array1) != len(array2):
+        raise ValueError("Arrays must have the same length")
+
+    percent_diff = []
+    for i in range(len(array1)):
+        if len(array1[i]) != 3 or len(array2[i]) != 3:
+            raise ValueError("Subarrays must have length 3")
+        
+        diff = [array2[i][j] - array1[i][j] for j in range(3)]
+        percent_diff.append([(diff[j] / array1[i][j]) * 100 if array1[i][j] != 0 else 0 for j in range(3)])
+
+    percent_diff_mag = []
+    for i in range(len(array1)):
+        percent_diff_mag.append(magnitude(percent_diff[i]))
+
+    return percent_diff_mag
+
+# Calculate the magnitude of a vector
+def magnitude(array):
+    mag = math.sqrt(sum(x**2 for x in array))
+    return mag
+
+# Run each test
+for i in range(len(solvers)):
+    for j in range(len(tests[i])):
+        
+        # Run simulation
+        print("Running "+tests[i][j])
+        os.system(executables[i] + ' ' + inputs[i][j])
+
+        GT_positions = extract_vector_data(standard_results[i][j], position_keyword)
+
+        # Read simulation results
+        results_filename = "vtk/data/VTK0.vtk"
+
+        if os.path.exists(results_filename):
+            print("Simulation Finished")
+        else:
+            print("Simulation did not finish")
+            raise ValueError("Simulation did not finish")
+
+        results_positions = extract_vector_data(results_filename, position_keyword)
+        position_diff = percent_difference_vectors(GT_positions, results_positions)
+
+
+        for k in range(len(position_diff)):
+            if position_diff[k] >= 1.0e-6:
+                raise ValueError(" ****************** ERROR: Position difference out of range for "+tests[i][j]+" problem ****************** ")
+
+        print("Removing simulation outputs")
+        os.system('rm -rf  vtk' )
\ No newline at end of file
diff --git a/lib/Elements b/lib/Elements
index 19b5cad11..38bc48fa8 160000
--- a/lib/Elements
+++ b/lib/Elements
@@ -1 +1 @@
-Subproject commit 19b5cad11d411a914cf188814afb2bd6180c85b7
+Subproject commit 38bc48fa84882a975e5a80a8781af9f51cdd73d7
diff --git a/scripts/build-fierro.sh b/scripts/build-fierro.sh
index 83a555bf8..e6dbff417 100755
--- a/scripts/build-fierro.sh
+++ b/scripts/build-fierro.sh
@@ -47,6 +47,11 @@ show_help() {
     echo "          mac                         A Mac computer. This option does not allow for cuda and hip builds, and build_cores will be set to 1"
     echo "          msu                         A linux computer managed by the HPCC group at Mississippi State University"
     echo " "
+    echo "      --intel_mkl                     Decides whether to build Trilinos using the Intel MKL library"
+    echo " "
+    echo "          enabled                     Links and builds Trilinos with the Intel MKL library"
+    echo "          disabled                    Links and builds Trilinos using LAPACK and BLAS"
+    echo " "
     echo "      --heffte_build_type             The build type for the heffte installation. The default is 'fftw'"
     echo " "
     echo "          fftw                        General heffte run type"
diff --git a/scripts/trilinos-install.sh b/scripts/trilinos-install.sh
index 7bc7e6768..a057bc527 100644
--- a/scripts/trilinos-install.sh
+++ b/scripts/trilinos-install.sh
@@ -26,6 +26,13 @@ then
     mkdir -p ${TRILINOS_BUILD_DIR} 
 fi
 
+if [ "$kokkos_build_type" = "cuda" ]; then
+    export OMPI_CXX=${TRILINOS_SOURCE_DIR}/packages/kokkos/bin/nvcc_wrapper
+    export CUDA_LAUNCH_BLOCKING=1
+elif [ "$kokkos_build_type" = *"hip"* ]; then
+    export OMPI_CXX=hipcc
+fi
+
 #check if Trilinos library files were installed, install them otherwise.
 [ -d "${TRILINOS_BUILD_DIR}/lib" ] && echo "Directory ${TRILINOS_BUILD_DIR}/lib exists, assuming successful installation; delete build folder and run build script again if there was an environment error that has been corrected."
 
@@ -50,6 +57,7 @@ CUDA_ADDITIONS=(
 -D Tpetra_ENABLE_CUDA=ON
 -D Xpetra_ENABLE_Kokkos_Refactor=ON
 -D MueLu_ENABLE_Kokkos_Refactor=ON
+-D Tpetra_ASSUME_GPU_AWARE_MPI:BOOL=FALSE
 )
 
 # Kokkos flags for Hip
@@ -64,6 +72,7 @@ export OMPI_CXX=hipcc
 -D KokkosKernels_ENABLE_TPL_CUSPARSE=OFF
 -D Tpetra_INST_HIP=ON
 -D Xpetra_ENABLE_Kokkos_Refactor=ON
+-D Tpetra_ASSUME_GPU_AWARE_MPI:BOOL=FALSE
 )
 
 # Kokkos flags for OpenMP
@@ -141,13 +150,10 @@ if [ "$kokkos_build_type" = "openmp" ]; then
         ${OPENMP_ADDITIONS[@]}
     )
 elif [ "$kokkos_build_type" = "cuda" ]; then
-    export OMPI_CXX=${TRILINOS_SOURCE_DIR}/packages/kokkos/bin/nvcc_wrapper
-    export CUDA_LAUNCH_BLOCKING=1
     cmake_options+=(
         ${CUDA_ADDITIONS[@]}
     )
 elif [ "$kokkos_build_type" = "hip" ]; then
-    export OMPI_CXX=hipcc
     cmake_options+=(
         ${HIP_ADDITIONS[@]}
     )
diff --git a/single-node-refactor/CMakeLists.txt b/single-node-refactor/CMakeLists.txt
index c8f008e1f..18c3e7bdd 100755
--- a/single-node-refactor/CMakeLists.txt
+++ b/single-node-refactor/CMakeLists.txt
@@ -40,6 +40,44 @@ add_subdirectory(../lib/Elements/matar cbin)
 # include_directories(Mesh-Builder)
 # add_subdirectory(Mesh-Builder)
 
+if (FIERRO_ENABLE_TRILINOS)
+  find_package(Trilinos REQUIRED) #new
+  # Assume if the CXX compiler exists, the rest do too.
+  if (EXISTS ${Trilinos_CXX_COMPILER})
+    set(CMAKE_CXX_COMPILER ${Trilinos_CXX_COMPILER})
+    set(CMAKE_C_COMPILER ${Trilinos_C_COMPILER})
+    set(CMAKE_Fortran_COMPILER ${Trilinos_Fortran_COMPILER})
+  endif()
+  if(NOT DISTRIBUTION)
+    # Make sure to use same compilers and flags as Trilinos
+    set(CMAKE_CXX_FLAGS  "${Trilinos_CXX_COMPILER_FLAGS} ${CMAKE_CXX_FLAGS}")
+    set(CMAKE_C_FLAGS  "${Trilinos_C_COMPILER_FLAGS} ${CMAKE_C_FLAGS}")
+    set(CMAKE_Fortran_FLAGS  "${Trilinos_Fortran_COMPILER_FLAGS} ${CMAKE_Fortran_FLAGS}")
+  endif()
+
+  message("\nFound Trilinos!  Here are the details: ")
+  message("   Trilinos_DIR = ${Trilinos_DIR}")
+  message("   Trilinos_VERSION = ${Trilinos_VERSION}")
+  message("   Trilinos_PACKAGE_LIST = ${Trilinos_PACKAGE_LIST}")
+  message("   Trilinos_LIBRARIES = ${Trilinos_LIBRARIES}")
+  message("   Trilinos_INCLUDE_DIRS = ${Trilinos_INCLUDE_DIRS}")
+  message("   Trilinos_LIBRARY_DIRS = ${Trilinos_LIBRARY_DIRS}")
+  message("   Trilinos_TPL_LIST = ${Trilinos_TPL_LIST}")
+  message("   Trilinos_TPL_INCLUDE_DIRS = ${Trilinos_TPL_INCLUDE_DIRS}")
+  message("   Trilinos_TPL_LIBRARIES = ${Trilinos_TPL_LIBRARIES}")
+  message("   Trilinos_TPL_LIBRARY_DIRS = ${Trilinos_TPL_LIBRARY_DIRS}")
+  message("   Trilinos_BUILD_SHARED_LIBS = ${Trilinos_BUILD_SHARED_LIBS}")
+  message("End of Trilinos details\n")
+
+  include_directories(${Trilinos_INCLUDE_DIRS} ${Trilinos_TPL_INCLUDE_DIRS})
+  list(APPEND LINKING_LIBRARIES Trilinos::all_selected_libs)
+  add_definitions(-DTRILINOS_INTERFACE=1 -DHAVE_MPI=1)
+else()
+  find_package(Kokkos REQUIRED)
+  list(APPEND LINKING_LIBRARIES Kokkos::kokkos)
+endif()
+find_package(Matar REQUIRED)
+
 include_directories(src/material_models/artificial_viscosity)
 include_directories(src/material_models/eos)
 include_directories(src/material_models/erosion)
diff --git a/single-node-refactor/scripts/build-fierro.sh b/single-node-refactor/scripts/build-fierro.sh
index fa0b2ebe6..9dcaa18b7 100755
--- a/single-node-refactor/scripts/build-fierro.sh
+++ b/single-node-refactor/scripts/build-fierro.sh
@@ -7,6 +7,8 @@ show_help() {
     echo "  --build_action=<full-app|set-env|install-kokkos|fierro>. Default is 'full-app'"
     echo "  --machine=<darwin|chicoma|linux|mac>. Default is 'linux'"
     echo "  --build_cores=<Integers greater than 0>. Default is set 1"
+    echo "  --trilinos=<enabled|disabled>. Default is 'disabled'"
+    echo "  --intel_mkl=<enabled|disabled>. Default is 'disabled'"
     echo "  --help: Display this help message"
     echo " "
     echo " "
@@ -36,6 +38,16 @@ show_help() {
     echo "          linux                       A general linux machine (that does not use modules)"
     echo "          mac                         A Mac computer. This option does not allow for cuda and hip builds, and build_cores will be set to 1"
     echo " "
+    echo "      --trilinos                      Decides if Trilinos is available for certain MATAR functionality"
+    echo " "
+    echo "          disabled                    Trilinos is not being used"
+    echo "          enabled                     Trilinos will be linked with MATAR to enable relevant functionality"
+    echo " "
+    echo "      --intel_mkl                     Decides whether to build Trilinos using the Intel MKL library"
+    echo " "
+    echo "          enabled                     Links and builds Trilinos with the Intel MKL library"
+    echo "          disabled                    Links and builds Trilinos using LAPACK and BLAS"
+    echo " "
     echo "      --build_cores                   The number of build cores to be used by make and make install commands. The default is 1" 
     echo " "
     echo "      --debug                         Build with debug. Default is false." 
@@ -48,6 +60,8 @@ solver="SGH"
 machine="linux"
 kokkos_build_type="openmp"
 build_cores="1"
+trilinos="disabled"
+intel_mkl="disabled"
 debug="false"
 
 # Define arrays of valid options
@@ -55,6 +69,8 @@ valid_build_action=("full-app" "set-env" "install-kokkos" "fierro")
 valid_solver=("SGH")
 valid_kokkos_build_types=("serial" "openmp" "pthreads" "cuda" "hip")
 valid_machines=("darwin" "chicoma" "linux" "mac")
+valid_trilinos=("disabled" "enabled")
+valid_intel_mkl=("disabled" "enabled")
 valid_debug=("true" "false")
 
 # Parse command line arguments
@@ -110,6 +126,26 @@ for arg in "$@"; do
                 return 1
             fi
             ;;
+        --trilinos=*)
+            option="${arg#*=}"
+            if [[ " ${valid_trilinos[*]} " == *" $option "* ]]; then
+                trilinos="$option"
+            else
+                echo "Error: Invalid --kokkos_build_type specified."
+                show_help
+                return 1
+            fi
+            ;;
+        --intel_mkl=*)
+            option="${arg#*=}"
+            if [[ " ${valid_intel_mkl[*]} " == *" $option "* ]]; then
+                intel_mkl="$option"
+            else
+                echo "Error: Invalid --intel_mkl specified."
+                show_help
+                return 1
+            fi
+            ;;
         --debug=*)
             option="${arg#*=}"
             if [[ " ${valid_debug[*]} " == *" $option "* ]]; then
@@ -155,6 +191,8 @@ echo "Building based on these argument options:"
 echo "Build action - ${build_action}"
 echo "Solver - ${solver}"
 echo "Kokkos backend - ${kokkos_build_type}"
+echo "Trilinos - ${trilinos}"
+echo "Intel MKL library - ${intel_mkl}"
 echo "make -j ${build_cores}"
 
 cd "$( dirname "${BASH_SOURCE[0]}" )"
@@ -164,13 +202,17 @@ source setup-env.sh ${machine} ${kokkos_build_type} ${build_cores}
 
 # Next, do action based on args
 if [ "$build_action" = "full-app" ]; then
-    source kokkos-install.sh ${kokkos_build_type} ${debug}
-    source matar-install.sh ${kokkos_build_type} ${debug}
-    source cmake_build.sh ${solver} ${debug}
+    if [ "$trilinos" = "disabled" ]; then    
+        source kokkos-install.sh ${kokkos_build_type} ${debug} 
+    elif [ "$trilinos" = "enabled" ]; then    
+        source trilinos-install.sh ${kokkos_build_type}  ${intel_mkl} ${debug}
+    fi
+    source matar-install.sh ${kokkos_build_type} ${debug} ${trilinos}
+    source cmake_build.sh ${solver} ${debug} ${trilinos}
 elif [ "$build_action" = "install-kokkos" ]; then
     source kokkos-install.sh ${kokkos_build_type}
 elif [ "$build_action" = "fierro" ]; then
-    source cmake_build.sh ${solver}
+    source cmake_build.sh ${solver} ${debug} ${trilinos}
 else
     echo "No build action, only setup the environment."
 fi
diff --git a/single-node-refactor/scripts/cmake_build.sh b/single-node-refactor/scripts/cmake_build.sh
index 27f7f5608..96061c8bb 100755
--- a/single-node-refactor/scripts/cmake_build.sh
+++ b/single-node-refactor/scripts/cmake_build.sh
@@ -2,16 +2,31 @@
 
 solver="${1}"
 debug="${2}"
+trilinos="${3}"
 
 echo "Removing old Kokkos build and installation directory"
 rm -rf ${SGH_BUILD_DIR}
 mkdir -p ${SGH_BUILD_DIR}
 
-cmake_options=(
--D BUILD_EXPLICIT_SOLVER=OFF
--D CMAKE_PREFIX_PATH="${MATAR_INSTALL_DIR};${KOKKOS_INSTALL_DIR}"
-#-D CMAKE_CXX_FLAGS="-I${matardir}/src"
-)
+
+if [ "$trilinos" = "enabled" ]; then
+    if [ ! -d "${TRILINOS_INSTALL_DIR}/lib" ]; then
+        Trilinos_DIR=${TRILINOS_INSTALL_DIR}/lib64/cmake/Trilinos
+    else
+        Trilinos_DIR=${TRILINOS_INSTALL_DIR}/lib/cmake/Trilinos
+    fi
+    cmake_options+=(
+        -D CMAKE_PREFIX_PATH="${MATAR_INSTALL_DIR}"
+        -D Trilinos_DIR="$Trilinos_DIR"
+        -D FIERRO_ENABLE_TRILINOS=ON
+    )
+else
+    cmake_options=(
+    -D BUILD_EXPLICIT_SOLVER=OFF
+    -D CMAKE_PREFIX_PATH="${MATAR_INSTALL_DIR};${KOKKOS_INSTALL_DIR}"
+    #-D CMAKE_CXX_FLAGS="-I${matardir}/src"
+    )
+fi
 
 if [ "$debug" = "true" ]; then
     echo "Setting debug to true for CMAKE build type"
diff --git a/single-node-refactor/scripts/matar-install.sh b/single-node-refactor/scripts/matar-install.sh
index 525182ce4..4b4afde71 100755
--- a/single-node-refactor/scripts/matar-install.sh
+++ b/single-node-refactor/scripts/matar-install.sh
@@ -8,7 +8,6 @@ mkdir -p ${MATAR_BUILD_DIR}
 
 cmake_options=(
     -D CMAKE_INSTALL_PREFIX="${MATAR_INSTALL_DIR}"
-    -D CMAKE_PREFIX_PATH="${KOKKOS_INSTALL_DIR}"
 )
 
 if [ "$debug" = "true" ]; then
@@ -19,13 +18,37 @@ if [ "$debug" = "true" ]; then
     )
 fi
 
+
 if [ "$kokkos_build_type" = "none" ]; then
     cmake_options+=(
         -D Matar_ENABLE_KOKKOS=OFF
     )
+elif [ "$trilinos" = "enabled" ]; then
+    if [ ! -d "${TRILINOS_INSTALL_DIR}/lib" ]; then
+        Trilinos_DIR=${TRILINOS_INSTALL_DIR}/lib64/cmake/Trilinos
+    else
+        Trilinos_DIR=${TRILINOS_INSTALL_DIR}/lib/cmake/Trilinos
+    fi
+    cmake_options+=(
+        -D Trilinos_DIR="$Trilinos_DIR"
+        -D Matar_ENABLE_TRILINOS=ON
+        -D Matar_ENABLE_KOKKOS=ON
+    )
 else
     cmake_options+=(
         -D Matar_ENABLE_KOKKOS=ON
+        -D CMAKE_PREFIX_PATH="${KOKKOS_INSTALL_DIR}"
+    )
+    if [ "$kokkos_build_type" = "cuda" ]; then
+        cmake_options+=(
+            -D Matar_CUDA_BUILD=ON
+        )
+    fi
+fi
+
+if [[ "$kokkos_build_type" = *"mpi"* ]] || [ "$trilinos" = "enabled" ]; then
+    cmake_options+=(
+    -D Matar_ENABLE_MPI=ON
     )
 fi
 
diff --git a/single-node-refactor/scripts/setup-env.sh b/single-node-refactor/scripts/setup-env.sh
index 866625b52..461429719 100755
--- a/single-node-refactor/scripts/setup-env.sh
+++ b/single-node-refactor/scripts/setup-env.sh
@@ -29,6 +29,10 @@ export KOKKOS_SOURCE_DIR=${MATAR_SOURCE_DIR}/src/Kokkos/kokkos
 export KOKKOS_BUILD_DIR=${builddir}/kokkos
 export KOKKOS_INSTALL_DIR=${installdir}/kokkos-${kokkos_build_type}
 
+export TRILINOS_SOURCE_DIR=${libdir}/Trilinos
+export TRILINOS_BUILD_DIR=${TRILINOS_SOURCE_DIR}/build-${kokkos_build_type}
+export TRILINOS_INSTALL_DIR=${TRILINOS_BUILD_DIR}
+
 export FIERRO_BUILD_CORES=$build_cores
 
 cd $scriptdir
diff --git a/single-node-refactor/scripts/trilinos-install.sh b/single-node-refactor/scripts/trilinos-install.sh
new file mode 100644
index 000000000..fb32c9005
--- /dev/null
+++ b/single-node-refactor/scripts/trilinos-install.sh
@@ -0,0 +1,177 @@
+#!/bin/bash -e
+
+kokkos_build_type="${1}"
+intel_mkl="${2}"
+debug="${3}"
+
+# If all arguments are valid, you can use them in your script as needed
+echo "Trilinos Kokkos Build Type: $kokkos_build_type"
+
+#check if Trilinos directory exists, git clone Trilinos if it doesn't
+[ -d "${TRILINOS_SOURCE_DIR}" ] && echo "Directory Trilinos exists, skipping Trilinos download"
+
+if [ ! -d "${TRILINOS_SOURCE_DIR}" ]
+then
+  echo "Directory Trilinos does not exist, downloading Trilinos...."
+  git clone --depth 1 https://github.com/trilinos/Trilinos.git ${TRILINOS_SOURCE_DIR}
+fi
+
+#check if Trilinos build directory exists, create Trilinos/build if it doesn't
+[ -d "${TRILINOS_BUILD_DIR}" ] && echo "Directory ${TRILINOS_BUILD_DIR} exists, moving on"
+
+if [ ! -d "${TRILINOS_BUILD_DIR}" ]
+then
+  echo "Directory ${TRILINOS_BUILD_DIR} does not exist, creating it...."
+    rm -rf ${TRILINOS_BUILD_DIR} ${TRILINOS_INSTALL_DIR}
+    mkdir -p ${TRILINOS_BUILD_DIR} 
+fi
+
+
+if [ "$kokkos_build_type" = "cuda" ] || [ "$kokkos_build_type" = "cuda_mpi" ]; then
+    export OMPI_CXX=${TRILINOS_SOURCE_DIR}/packages/kokkos/bin/nvcc_wrapper
+    export CUDA_LAUNCH_BLOCKING=1
+elif [ "$kokkos_build_type" = *"hip"* ] || [ "$kokkos_build_type" = *"hip_mpi"* ]; then
+    export OMPI_CXX=hipcc
+fi
+
+#check if Trilinos library files were installed, install them otherwise.
+[ -d "${TRILINOS_BUILD_DIR}/lib" ] && echo "Directory ${TRILINOS_BUILD_DIR}/lib exists, assuming successful installation; delete build folder and run build script again if there was an environment error that has been corrected."
+
+[ -d "${TRILINOS_BUILD_DIR}/lib64" ] && echo "Directory ${TRILINOS_BUILD_DIR}/lib64 exists, assuming successful installation; delete build folder and run build script again if there was an environment error that has been corrected."
+
+if [ ! -d "${TRILINOS_BUILD_DIR}/lib" ] && [ ! -d "${TRILINOS_BUILD_DIR}/lib64" ]
+then
+  echo "Directory Trilinos/build/lib does not exist, compiling Trilinos (this might take a while)...."
+
+CUDA_ADDITIONS=(
+-D TPL_ENABLE_CUDA=ON
+-D TPL_ENABLE_CUBLAS=ON
+-D TPL_ENABLE_CUSPARSE=ON
+-D Kokkos_ENABLE_CUDA=ON
+-D Kokkos_ENABLE_CUDA_LAMBDA=ON
+-D Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE=ON
+-D Kokkos_ENABLE_DEPRECATED_CODE=OFF
+-D Kokkos_ENABLE_CUDA_UVM=OFF
+-D Trilinos_ENABLE_KokkosKernels=ON
+-D KokkosKernels_ENABLE_TPL_CUBLAS=ON
+-D KokkosKernels_ENABLE_TPL_CUSPARSE=ON
+-D Tpetra_ENABLE_CUDA=ON
+-D MueLu_ENABLE_Kokkos_Refactor=OFF
+-D Tpetra_ASSUME_GPU_AWARE_MPI:BOOL=FALSE
+)
+
+# Kokkos flags for Hip
+HIP_ADDITIONS=(
+export OMPI_CXX=hipcc
+-D Kokkos_ENABLE_HIP=ON
+-D Kokkos_ENABLE_HIP_RELOCATABLE_DEVICE_CODE=ON
+-D Kokkos_ENABLE_DEPRECATED_CODE=OFF
+-D Kokkos_ARCH_VEGA90A=ON
+-D Trilinos_ENABLE_KokkosKernels=ON
+-D KokkosKernels_ENABLE_TPL_CUBLAS=OFF
+-D KokkosKernels_ENABLE_TPL_CUSPARSE=OFF
+-D Tpetra_INST_HIP=ON
+-D Tpetra_ASSUME_GPU_AWARE_MPI:BOOL=FALSE
+)
+
+# Kokkos flags for OpenMP
+OPENMP_ADDITIONS=(
+-D Trilinos_ENABLE_OpenMP=ON
+)
+
+# Flags for building with MKL, which is supported at MSU HPCC
+MSU_ADDITIONS=(
+-D BLAS_LIBRARY_NAMES="libmkl_rt.so"
+-D BLAS_LIBRARY_DIRS="/apps/spack-managed/gcc-11.3.1/intel-oneapi-mkl-2022.2.1-7l7jlsd56x2kljiskrcvsoenmq4y3cu7/mkl/2022.2.1/lib/intel64"
+-D LAPACK_LIBRARY_NAMES="libmkl_rt.so"
+-D LAPACK_LIBRARY_DIRS="/apps/spack-managed/gcc-11.3.1/intel-oneapi-mkl-2022.2.1-7l7jlsd56x2kljiskrcvsoenmq4y3cu7/mkl/2022.2.1/lib/intel64"
+-D TPL_ENABLE_MKL:BOOL=ON
+-D MKL_LIBRARY_DIRS:FILEPATH="/apps/spack-managed/gcc-11.3.1/intel-oneapi-mkl-2022.2.1-7l7jlsd56x2kljiskrcvsoenmq4y3cu7/mkl/2022.2.1/lib/intel64"
+-D MKL_LIBRARY_NAMES:STRING="mkl_rt"
+-D MKL_INCLUDE_DIRS:FILEPATH="/apps/spack-managed/gcc-11.3.1/intel-oneapi-mkl-2022.2.1-7l7jlsd56x2kljiskrcvsoenmq4y3cu7/mkl/2022.2.1/include"
+)
+
+# Configure kokkos using CMake
+cmake_options=(
+-D CMAKE_BUILD_TYPE=Release
+-D Trilinos_MUST_FIND_ALL_TPL_LIBS=TRUE
+-D CMAKE_CXX_STANDARD=17
+-D TPL_ENABLE_MPI=ON
+)
+
+echo "**** Machine = ${machine} ****"
+if [ "$machine" = "msu" ]; then
+    echo "**** WARNING: Verify MKL path in trilinos-install.sh ****"
+    cmake_options+=(
+        ${MSU_ADDITIONS[@]}
+    )
+fi
+
+cmake_options+=(
+-D Trilinos_ENABLE_Kokkos=ON
+${ADDITIONS[@]}
+-D Trilinos_ENABLE_Amesos2=OFF
+-D Trilinos_ENABLE_Belos=OFF
+-D Trilinos_ENABLE_MueLu=OFF 
+-D Trilinos_ENABLE_ROL=OFF 
+-D Trilinos_ENABLE_Ifpack2=OFF
+-D Trilinos_ENABLE_Zoltan2=ON 
+-D Trilinos_ENABLE_Anasazi=OFF 
+-D MueLu_ENABLE_TESTS=OFF 
+-D Trilinos_ENABLE_ALL_PACKAGES=OFF 
+-D Trilinos_ENABLE_ALL_OPTIONAL_PACKAGES=OFF 
+-D Trilinos_ENABLE_TESTS=OFF 
+-D CMAKE_INSTALL_PREFIX=${TRILINOS_INSTALL_DIR}
+-D Xpetra_ENABLE_Kokkos_Refactor=ON 
+)
+
+# Flags for building with Intel MKL library
+INTEL_MKL_ADDITIONS=(
+-D TPL_ENABLE_MKL=ON
+-D BLAS_LIBRARY_NAMES="libmkl_rt.so"
+-D BLAS_LIBRARY_DIRS="$MKLROOT/lib/intel64"
+-D LAPACK_LIBRARY_NAMES="libmkl_rt.so"
+-D LAPACK_LIBRARY_DIRS="$MKLROOT/lib/intel64"
+-D MKL_LIBRARY_DIRS="$MKLROOT/lib/intel64"
+-D MKL_LIBRARY_NAMES="mkl_rt"
+-D MKL_INCLUDE_DIRS="$MKLROOT/include"
+)
+
+echo "**** Intel MKL = ${intel_mkl} ****"
+if [ "$intel_mkl" = "enabled" ]; then
+    echo "**** assuming MKL installation at $MKLROOT ****"
+    cmake_options+=(
+        ${INTEL_MKL_ADDITIONS[@]}
+    )
+fi
+
+if [ "$kokkos_build_type" = "openmp" ] || [ "$kokkos_build_type" = "openmp_mpi" ]; then
+    cmake_options+=(
+        ${OPENMP_ADDITIONS[@]}
+    )
+elif [ "$kokkos_build_type" = "cuda" ] || [ "$kokkos_build_type" = "cuda_mpi" ]; then
+    cmake_options+=(
+        ${CUDA_ADDITIONS[@]}
+    )
+elif [ "$kokkos_build_type" = *"hip"* ] || [ "$kokkos_build_type" = *"hip_mpi"* ]; then
+    cmake_options+=(
+        ${HIP_ADDITIONS[@]}
+    )
+fi
+
+# Print CMake options for reference
+echo "CMake Options: ${cmake_options[@]}"
+
+# Configure Trilinos
+cmake "${cmake_options[@]}" -B "${TRILINOS_BUILD_DIR}" -S "${TRILINOS_SOURCE_DIR}"
+
+# Build Trilinos
+echo "Building Trilinos..."
+make -C "${TRILINOS_BUILD_DIR}" -j${MATAR_BUILD_CORES}
+
+# Install Trilinos
+echo "Installing Trilinos..."
+make -C "${TRILINOS_BUILD_DIR}" install all
+
+echo "Trilinos installation complete."
+fi
diff --git a/single-node-refactor/src/CMakeLists.txt b/single-node-refactor/src/CMakeLists.txt
index 8d06a7f19..8448379b5 100755
--- a/single-node-refactor/src/CMakeLists.txt
+++ b/single-node-refactor/src/CMakeLists.txt
@@ -35,13 +35,6 @@ set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} ${VECTOR_Fortran
 set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} ${VECTOR_C_FLAGS}")
 set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${VECTOR_CXX_FLAGS}")
 
-find_package(Kokkos REQUIRED)
-find_package(Matar REQUIRED)
-
-
-
-
-
 
 message("\n ****** ADDING FIERRO EXECUTABLE ******** \n ")
 
diff --git a/single-node-refactor/src/Solvers/SGH_solver_3D/CMakeLists.txt b/single-node-refactor/src/Solvers/SGH_solver_3D/CMakeLists.txt
index 5a17c0336..053340d61 100755
--- a/single-node-refactor/src/Solvers/SGH_solver_3D/CMakeLists.txt
+++ b/single-node-refactor/src/Solvers/SGH_solver_3D/CMakeLists.txt
@@ -1,8 +1,5 @@
 cmake_minimum_required(VERSION 3.1.3)
 
-find_package(Matar REQUIRED)
-find_package(Kokkos REQUIRED)
-
 add_definitions(-DHAVE_KOKKOS=1)
 
 if (CUDA)
diff --git a/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp b/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp
index 378604381..f6e4a04d6 100644
--- a/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp
@@ -561,7 +561,7 @@ double sum_domain_internal_energy(const DCArrayKokkos<double>& MaterialPoints_ma
     double IE_loc_sum;
 
     // loop over the material points and tally IE
-    REDUCE_SUM(matpt_lid, 0, num_mat_points, IE_loc_sum, {
+    FOR_REDUCE_SUM(matpt_lid, 0, num_mat_points, IE_loc_sum, {
         IE_loc_sum += MaterialPoints_mass(matpt_lid) * MaterialPoints_sie(1, matpt_lid);
     }, IE_sum);
     Kokkos::fence();
@@ -594,7 +594,7 @@ double sum_domain_kinetic_energy(const Mesh_t& mesh,
     double KE_sum = 0.0;
     double KE_loc_sum;
 
-    REDUCE_SUM(node_gid, 0, mesh.num_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM(node_gid, 0, mesh.num_nodes, KE_loc_sum, {
         double ke = 0;
 
         for (size_t dim = 0; dim < mesh.num_dims; dim++) {
@@ -616,7 +616,7 @@ double sum_domain_material_mass(const DCArrayKokkos<double>& MaterialPoints_mass
     double mass_domain = 0.0;
     double mass_loc_domain;
 
-    REDUCE_SUM(matpt_lid, 0, num_mat_points, mass_loc_domain, {
+    FOR_REDUCE_SUM(matpt_lid, 0, num_mat_points, mass_loc_domain, {
         mass_loc_domain += MaterialPoints_mass(matpt_lid);
     }, mass_domain);
     Kokkos::fence();
@@ -647,7 +647,7 @@ double sum_domain_node_mass(const Mesh_t& mesh,
     double mass_domain = 0.0;
     double mass_loc_domain;
 
-    REDUCE_SUM(node_gid, 0, mesh.num_nodes, mass_loc_domain, {
+    FOR_REDUCE_SUM(node_gid, 0, mesh.num_nodes, mass_loc_domain, {
         if (mesh.num_dims == 2) {
             mass_loc_domain += node_mass(node_gid) * node_coords(1, node_gid, 1);
         }
diff --git a/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_setup.cpp b/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_setup.cpp
index 0bdb70893..efa993498 100644
--- a/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_setup.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver_3D/src/sgh_setup.cpp
@@ -310,7 +310,7 @@ void SGH3D::setup(SimulationParameters_t& SimulationParamaters,
         size_t sum_local;
         size_t sum_total;
 
-        REDUCE_SUM(elem_gid, 0, num_elems, sum_local, {
+        FOR_REDUCE_SUM(elem_gid, 0, num_elems, sum_local, {
             if (elem_mat_id(elem_gid) == mat_id) {
                 // increment the number of elements the materials live in
                 sum_local++;
diff --git a/single-node-refactor/src/Solvers/SGH_solver_3D/src/time_integration.cpp b/single-node-refactor/src/Solvers/SGH_solver_3D/src/time_integration.cpp
index ce00b2c30..4dff167ed 100644
--- a/single-node-refactor/src/Solvers/SGH_solver_3D/src/time_integration.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver_3D/src/time_integration.cpp
@@ -123,7 +123,7 @@ void SGH3D::get_timestep(Mesh_t& mesh,
 
     double dt_lcl;
     double min_dt_calc;
-    REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, {
+    FOR_REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, {
         size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid);
 
         double coords0[24];  // element coords
diff --git a/single-node-refactor/src/Solvers/SGH_solver_rz/CMakeLists.txt b/single-node-refactor/src/Solvers/SGH_solver_rz/CMakeLists.txt
index f83899dd6..940c9b1e1 100755
--- a/single-node-refactor/src/Solvers/SGH_solver_rz/CMakeLists.txt
+++ b/single-node-refactor/src/Solvers/SGH_solver_rz/CMakeLists.txt
@@ -1,8 +1,5 @@
 cmake_minimum_required(VERSION 3.1.3)
 
-find_package(Matar REQUIRED)
-find_package(Kokkos REQUIRED)
-
 add_definitions(-DHAVE_KOKKOS=1)
 
 if (CUDA)
diff --git a/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp b/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp
index 7e652b1d7..ce692ae27 100644
--- a/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp
@@ -518,7 +518,7 @@ double sum_domain_internal_energy_rz(const DCArrayKokkos<double>& MaterialPoints
     double IE_loc_sum;
 
     // loop over the material points and tally IE
-    REDUCE_SUM(matpt_lid, 0, num_mat_points, IE_loc_sum, {
+    FOR_REDUCE_SUM(matpt_lid, 0, num_mat_points, IE_loc_sum, {
         IE_loc_sum += MaterialPoints_mass(matpt_lid) * MaterialPoints_sie(1,matpt_lid);
     }, IE_sum);
     Kokkos::fence();
@@ -535,7 +535,7 @@ double sum_domain_kinetic_energy_rz(const Mesh_t& mesh,
     double KE_sum = 0.0;
     double KE_loc_sum;
 
-    REDUCE_SUM(node_gid, 0, mesh.num_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM(node_gid, 0, mesh.num_nodes, KE_loc_sum, {
         double ke = 0;
 
         for (size_t dim = 0; dim < mesh.num_dims; dim++) {
@@ -560,7 +560,7 @@ double sum_domain_material_mass_rz(const DCArrayKokkos<double>& MaterialPoints_m
     double mass_domain = 0.0;
     double mass_loc_domain;
 
-    REDUCE_SUM(matpt_lid, 0, num_mat_points, mass_loc_domain, {
+    FOR_REDUCE_SUM(matpt_lid, 0, num_mat_points, mass_loc_domain, {
 
             mass_loc_domain += MaterialPoints_mass(matpt_lid);
         
@@ -577,7 +577,7 @@ double sum_domain_node_mass_rz(const CArrayKokkos<double>& node_extensive_mass,
     double mass_domain = 0.0;
     double mass_loc_domain;
 
-    REDUCE_SUM(node_gid, 0, num_nodes, mass_loc_domain, {
+    FOR_REDUCE_SUM(node_gid, 0, num_nodes, mass_loc_domain, {
 
             mass_loc_domain += node_extensive_mass(node_gid);
         
diff --git a/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_setup_rz.cpp b/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_setup_rz.cpp
index 0d4c0398a..c5929db5b 100644
--- a/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_setup_rz.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver_rz/src/sgh_setup_rz.cpp
@@ -344,7 +344,7 @@ void SGHRZ::setup(SimulationParameters_t& SimulationParamaters,
         size_t sum_local;
         size_t sum_total;
 
-        REDUCE_SUM(elem_gid, 0, num_elems, sum_local,{
+        FOR_REDUCE_SUM(elem_gid, 0, num_elems, sum_local,{
 
             if(elem_mat_id(elem_gid) == mat_id){
                 // increment the number of elements the materials live in
diff --git a/single-node-refactor/src/Solvers/SGH_solver_rz/src/time_integration_rz.cpp b/single-node-refactor/src/Solvers/SGH_solver_rz/src/time_integration_rz.cpp
index 28fdfc603..c2f14bc6c 100644
--- a/single-node-refactor/src/Solvers/SGH_solver_rz/src/time_integration_rz.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver_rz/src/time_integration_rz.cpp
@@ -126,7 +126,7 @@ void SGHRZ::get_timestep_rz(Mesh_t& mesh,
 
     double dt_lcl;
     double min_dt_calc;
-    REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, {
+    FOR_REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, {
 
         size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); 
 
diff --git a/single-node-refactor/src/Solvers/SGTM_solver_3D/CMakeLists.txt b/single-node-refactor/src/Solvers/SGTM_solver_3D/CMakeLists.txt
index ec7e1adba..7f8e6adac 100755
--- a/single-node-refactor/src/Solvers/SGTM_solver_3D/CMakeLists.txt
+++ b/single-node-refactor/src/Solvers/SGTM_solver_3D/CMakeLists.txt
@@ -1,8 +1,5 @@
 cmake_minimum_required(VERSION 3.1.3)
 
-find_package(Matar REQUIRED)
-find_package(Kokkos REQUIRED)
-
 add_definitions(-DHAVE_KOKKOS=1)
 
 if (CUDA)
diff --git a/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp b/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp
index 0c942f7f6..1b7597545 100644
--- a/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp
+++ b/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp
@@ -409,7 +409,7 @@ double SGTM3D::sum_domain_internal_energy(const DCArrayKokkos<double>& MaterialP
     double IE_loc_sum;
 
     // loop over the material points and tally IE
-    REDUCE_SUM(matpt_lid, 0, num_mat_points, IE_loc_sum, {
+    FOR_REDUCE_SUM(matpt_lid, 0, num_mat_points, IE_loc_sum, {
         IE_loc_sum += MaterialPoints_mass(matpt_lid) * MaterialPoints_sie(1, matpt_lid);
     }, IE_sum);
     Kokkos::fence();
@@ -442,7 +442,7 @@ double SGTM3D::sum_domain_kinetic_energy(const Mesh_t& mesh,
     double KE_sum = 0.0;
     double KE_loc_sum;
 
-    REDUCE_SUM(node_gid, 0, mesh.num_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM(node_gid, 0, mesh.num_nodes, KE_loc_sum, {
         double ke = 0;
 
         for (size_t dim = 0; dim < mesh.num_dims; dim++) {
@@ -464,7 +464,7 @@ double SGTM3D::sum_domain_material_mass(const DCArrayKokkos<double>& MaterialPoi
     double mass_domain = 0.0;
     double mass_loc_domain;
 
-    REDUCE_SUM(matpt_lid, 0, num_mat_points, mass_loc_domain, {
+    FOR_REDUCE_SUM(matpt_lid, 0, num_mat_points, mass_loc_domain, {
         mass_loc_domain += MaterialPoints_mass(matpt_lid);
     }, mass_domain);
     Kokkos::fence();
@@ -495,7 +495,7 @@ double SGTM3D::sum_domain_node_mass(const Mesh_t& mesh,
     double mass_domain = 0.0;
     double mass_loc_domain;
 
-    REDUCE_SUM(node_gid, 0, mesh.num_nodes, mass_loc_domain, {
+    FOR_REDUCE_SUM(node_gid, 0, mesh.num_nodes, mass_loc_domain, {
         if (mesh.num_dims == 2) {
             mass_loc_domain += node_mass(node_gid) * node_coords(1, node_gid, 1);
         }
diff --git a/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_setup.cpp b/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_setup.cpp
index a442b73e8..230c286dd 100644
--- a/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_setup.cpp
+++ b/single-node-refactor/src/Solvers/SGTM_solver_3D/src/sgtm_setup.cpp
@@ -305,7 +305,7 @@ void SGTM3D::setup(SimulationParameters_t& SimulationParamaters,
         size_t sum_local;
         size_t sum_total;
 
-        REDUCE_SUM(elem_gid, 0, num_elems, sum_local, {
+        FOR_REDUCE_SUM(elem_gid, 0, num_elems, sum_local, {
             if (elem_mat_id(elem_gid) == mat_id) {
                 // increment the number of elements the materials live in
                 sum_local++;
diff --git a/single-node-refactor/src/Solvers/SGTM_solver_3D/src/time_integration.cpp b/single-node-refactor/src/Solvers/SGTM_solver_3D/src/time_integration.cpp
index a17c34355..a6d03df39 100644
--- a/single-node-refactor/src/Solvers/SGTM_solver_3D/src/time_integration.cpp
+++ b/single-node-refactor/src/Solvers/SGTM_solver_3D/src/time_integration.cpp
@@ -125,7 +125,7 @@ void SGTM3D::get_timestep(Mesh_t& mesh,
 
     double dt_lcl;
     double min_dt_calc;
-    REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, {
+    FOR_REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, {
         size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid);
 
         double coords0[24];  // element coords
diff --git a/single-node-refactor/src/common/CMakeLists.txt b/single-node-refactor/src/common/CMakeLists.txt
index c8824f6e1..527169c4b 100644
--- a/single-node-refactor/src/common/CMakeLists.txt
+++ b/single-node-refactor/src/common/CMakeLists.txt
@@ -1,8 +1,5 @@
 cmake_minimum_required(VERSION 3.1.3)
 
-find_package(Matar REQUIRED)
-find_package(Kokkos REQUIRED)
-
 add_definitions(-DHAVE_KOKKOS=1)
 
 if (CUDA)
diff --git a/single-node-refactor/src/common/include/mesh.h b/single-node-refactor/src/common/include/mesh.h
index 8a84ba624..f92272a52 100644
--- a/single-node-refactor/src/common/include/mesh.h
+++ b/single-node-refactor/src/common/include/mesh.h
@@ -429,7 +429,7 @@ struct Mesh_t
         // find the max number of elems around a node
         size_t max_num_elems_in_node;
         size_t max_num_lcl;
-        REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
+        FOR_REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
             // num_corners_in_node = num_elems_in_node
             size_t max_num = num_corners_in_node(node_gid);
 
@@ -1299,7 +1299,7 @@ struct Mesh_t
         // find the max number of elems around a node
         size_t max_num_elems_in_node;
         size_t max_num_lcl;
-        REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
+        FOR_REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
             // num_corners_in_node = num_elems_in_node
             size_t max_num = num_corners_in_node(node_gid);
 
diff --git a/single-node-refactor/src/input/CMakeLists.txt b/single-node-refactor/src/input/CMakeLists.txt
index a3015a0e3..174196873 100644
--- a/single-node-refactor/src/input/CMakeLists.txt
+++ b/single-node-refactor/src/input/CMakeLists.txt
@@ -1,9 +1,6 @@
 cmake_minimum_required(VERSION 3.1.3)
 
 
-find_package(Matar REQUIRED)
-find_package(Kokkos REQUIRED)
-
 add_definitions(-DHAVE_KOKKOS=1)
 
 if (CUDA)
diff --git a/src/EVPFFT/src/evpfft.cpp b/src/EVPFFT/src/evpfft.cpp
index 9a782251e..7e3c1ba41 100644
--- a/src/EVPFFT/src/evpfft.cpp
+++ b/src/EVPFFT/src/evpfft.cpp
@@ -146,6 +146,11 @@ void EVPFFT::set_some_voxels_arrays_to_zero()
 void EVPFFT::init_after_reading_input_data()
 {
 
+    init_xk_gb();
+    init_disgradmacro();
+    init_ept();
+    init_evm();
+
 #ifndef ABSOLUTE_NO_OUTPUT
     if (iwfields == 1) {
       int imicro = 0;
@@ -154,11 +159,6 @@ void EVPFFT::init_after_reading_input_data()
     }
 #endif
 
-    init_xk_gb();
-    init_disgradmacro();
-    init_ept();
-    init_evm();
-
 // the variables initialized in the funcitons below are reduced into
 // and should be done once, hence the need for #if guard since the variables
 // needs to be initialized after udot and dt are know from fierro
diff --git a/src/EVPFFT/src/write_micro_state.cpp b/src/EVPFFT/src/write_micro_state.cpp
index d3820f5a0..4120b5702 100644
--- a/src/EVPFFT/src/write_micro_state.cpp
+++ b/src/EVPFFT/src/write_micro_state.cpp
@@ -474,30 +474,40 @@ void EVPFFT::write_micro_state_pvtu()
   edotp.update_host();
 
   // Calculate point positions
-  MatrixTypeRealDevice xtmp(3);
-  MatrixTypeRealDevice defgradavg(3,3);
+  MatrixTypeRealDual defgradavg(3,3);
   MatrixTypeRealDual xintp(3,npts1+1,npts2+1,npts3+1);
   MatrixTypeIntHost pid(npts1+1,npts2+1,npts3+1);
 
-  FOR_ALL(
+  for (int ii = 1; ii <= 3; ii++) {
+    for (int jj = 1; jj <= 3; jj++) {
+      if (imicro==0) {
+        defgradavg.host(ii,jj) = 0.0;
+      } else {
+        defgradavg.host(ii,jj) = disgradmacroactual(ii,jj);
+      }
+      if (ii==jj) {defgradavg.host(ii,jj) += 1.0;}
+    }
+  }
+  defgradavg.update_device();
+
+  FOR_ALL_CLASS(
           kz, 1, npts3+2,
           ky, 1, npts2+2,
           kx, 1, npts1+2, {
 
-        xtmp(1) = double(kx+local_start1)-0.5;
-        xtmp(2) = double(ky+local_start2)-0.5;
-        xtmp(3) = double(kz+local_start3)-0.5;
+        double xtmp[3];
+        xtmp[0] = (double)kx+(double)local_start1-0.5;
+        xtmp[1] = (double)ky+(double)local_start2-0.5;
+        xtmp[2] = (double)kz+(double)local_start3-0.5;
 
         for (int ii = 1; ii <= 3; ii++) {
           real_t dum = 0.0;
           for (int jj = 1; jj <= 3; jj++) {
-            defgradavg(ii,jj) = disgradmacroactual(ii,jj);
-            if (ii==jj) {defgradavg(ii,jj) += 1.0;}
-            dum += defgradavg(ii,jj)*xtmp(jj);
+            dum += defgradavg(ii,jj)*xtmp[jj-1];
           }
-          xintp(ii,kx,ky,kz) = dum;
+          xintp(ii,kx,ky,kz) = dum*delt(ii);
         }
-  }); // end FOR_ALL
+  }); // end FOR_ALL_CLASS
   Kokkos::fence();
   xintp.update_host();
 
@@ -699,7 +709,6 @@ void EVPFFT::write_micro_state_pvtu()
   out.open(filename,std::ofstream::binary);
 
   byte_offset = 0;
-  double crap = 0.0;
   
   //  Write Header
   //tmp_str = "<?xml version=\"1.0\"?>\n";
@@ -938,7 +947,7 @@ void EVPFFT::write_micro_state_pvtu()
       for (int kx = 1; kx <= npts1+1; kx++){
         ic += 1;
         for (int dim = 1; dim <= num_dims; dim++){
-          coord_tmp = xintp(dim,kx,ky,kz);
+          coord_tmp = xintp.host(dim,kx,ky,kz);
           out.write((char *) &coord_tmp,sizeof(coord_tmp));
         }
         pid(kx,ky,kz) = ic;
@@ -1006,7 +1015,6 @@ void EVPFFT::write_micro_state_pvtu()
   // for (int node_gid = 0; node_gid < num_points; node_gid++){
   //   for (int dim = 0; dim < num_dims; dim++){
   //     //out.write((char *) &node_vel.host(1, node_gid, dim),sizeof(double));
-  //     out.write((char *) &crap,sizeof(double));
   //   }
   // }
 
@@ -1169,7 +1177,6 @@ void EVPFFT::write_micro_state_pvtu()
   // for (int elem_gid = 0; elem_gid < num_cells; elem_gid++){
   //   for (int ii = 0; ii < 9; ii++){
   //     //out.write((char *) &elem_stress.host(1,elem_gid,ii),sizeof(double));
-  //     out.write((char *) &crap,sizeof(double));
   //   }
   // }
   // //  Density
diff --git a/src/LS-EVPFFT/src/evpfft.cpp b/src/LS-EVPFFT/src/evpfft.cpp
index 70d3accad..7b08c4e0c 100644
--- a/src/LS-EVPFFT/src/evpfft.cpp
+++ b/src/LS-EVPFFT/src/evpfft.cpp
@@ -150,6 +150,13 @@ void EVPFFT::set_some_voxels_arrays_to_zero()
 void EVPFFT::init_after_reading_input_data()
 {
 
+    init_xk_gb();
+    init_disgradmacro_velgradmacro();
+    init_ept();
+    init_disgrad();
+    init_evm();
+    init_defgrad();
+
 #ifndef ABSOLUTE_NO_OUTPUT
     if (iwfields == 1) {
       int imicro = 0;
@@ -158,13 +165,6 @@ void EVPFFT::init_after_reading_input_data()
     }
 #endif
 
-    init_xk_gb();
-    init_disgradmacro_velgradmacro();
-    init_ept();
-    init_disgrad();
-    init_evm();
-    init_defgrad();
-
 // the variables initialized in the funcitons below are reduced into
 // and should be done once, hence the need for #if guard since the variables
 // needs to be initialized after udot and dt are know from fierro
@@ -424,7 +424,7 @@ void EVPFFT::init_defgrad() {
       for (int ii = 1; ii <= 3; ii++) {
         defgradavg(jj,ii) = 0.0;
       }
-      defgradavg(jj,jj) = delt(jj);
+      defgradavg(jj,jj) = 1.0;
     } 
   }); // end FOR_ALL_CLASS
 
diff --git a/src/LS-EVPFFT/src/write_micro_state.cpp b/src/LS-EVPFFT/src/write_micro_state.cpp
index d4d02a25a..b86b70ab2 100644
--- a/src/LS-EVPFFT/src/write_micro_state.cpp
+++ b/src/LS-EVPFFT/src/write_micro_state.cpp
@@ -498,34 +498,42 @@ void EVPFFT::write_micro_state_pvtu()
   edotp.update_host();
 
   // Calculate point positions
-  MatrixTypeRealDevice xtmp(3);
+  MatrixTypeRealDual defgradavg_dual(3,3);
   MatrixTypeRealDual uf(3,npts1_g,npts2_g,npts3_g);
   MatrixTypeRealDual ufintp(3,npts1+1,npts2+1,npts3+1);
   MatrixTypeRealDual xintp(3,npts1+1,npts2+1,npts3+1);
   MatrixTypeIntHost pid(npts1+1,npts2+1,npts3+1);
 
-  FOR_ALL(
+  for (int ii = 1; ii <= 3; ii++) {
+    for (int jj = 1; jj <= 3; jj++) {
+      defgradavg_dual.host(ii,jj) = defgradavg(ii,jj);
+    }
+  }
+  defgradavg_dual.update_device();
+  
+  FOR_ALL_CLASS(
           kz, 1, npts3+1,
           ky, 1, npts2+1,
           kx, 1, npts1+1, {
 
-        xtmp(1) = double(kx);
-        xtmp(2) = double(ky);
-        xtmp(3) = double(kz);
+        double xtmp[3];
+        xtmp[0] = double(kx);
+        xtmp[1] = double(ky);
+        xtmp[2] = double(kz);
         for (int ii = 1; ii <= 3; ii++) {
           real_t dum = 0.0;
           for (int jj = 1; jj <= 3; jj++) {
-            dum += defgradavg(ii,jj)*xtmp(jj);
+            dum += defgradavg_dual(ii,jj)*xtmp[jj-1];
           }
           uf(ii,kx+local_start1,ky+local_start2,kz+local_start3) = xnode(ii,kx,ky,kz) - dum;
         }
-  }); // end FOR_ALL
+  }); // end FOR_ALL_CLASS
   Kokkos::fence();
   uf.update_host();
 
   MPI_Allreduce(MPI_IN_PLACE, uf.host_pointer(), n, MPI_REAL_T, MPI_SUM, mpi_comm);
 
-  FOR_ALL(
+  FOR_ALL_CLASS(
           kz, 1, npts3+2,
           ky, 1, npts2+2,
           kx, 1, npts1+2, {
@@ -572,27 +580,28 @@ void EVPFFT::write_micro_state_pvtu()
           ufintp(ii,kx,ky,kz) = 0.125*(uf(ii,kxn1,kyn1,kzn1)+uf(ii,kxn2,kyn1,kzn1)+uf(ii,kxn1,kyn2,kzn1)+uf(ii,kxn2,kyn2,kzn1)+
           uf(ii,kxn1,kyn1,kzn2)+uf(ii,kxn2,kyn1,kzn2)+uf(ii,kxn1,kyn2,kzn2)+uf(ii,kxn2,kyn2,kzn2) );
         }
-  }); // end FOR_ALL
+  }); // end FOR_ALL_CLASS
   Kokkos::fence();
   ufintp.update_host();
 
-  FOR_ALL(
+  FOR_ALL_CLASS(
           kz, 1, npts3+2,
           ky, 1, npts2+2,
           kx, 1, npts1+2, {
 
-        xtmp(1) = double(kx+local_start1)-0.5;
-        xtmp(2) = double(ky+local_start2)-0.5;
-        xtmp(3) = double(kz+local_start3)-0.5;
+        double xtmp[3];
+        xtmp[0] = double(kx+local_start1)-0.5;
+        xtmp[1] = double(ky+local_start2)-0.5;
+        xtmp[2] = double(kz+local_start3)-0.5;
 
         for (int ii = 1; ii <= 3; ii++) {
           real_t dum = 0.0;
           for (int jj = 1; jj <= 3; jj++) {
-            dum += defgradavg(ii,jj)*xtmp(jj);
+            dum += defgradavg_dual(ii,jj)*xtmp[jj-1];
           }
-          xintp(ii,kx,ky,kz) = dum + ufintp(ii,kx,ky,kz);
+          xintp(ii,kx,ky,kz) = (dum + ufintp(ii,kx,ky,kz))*delt(ii);
         }
-  }); // end FOR_ALL
+  }); // end FOR_ALL_CLASS
   Kokkos::fence();
   xintp.update_host();
 
@@ -794,7 +803,6 @@ void EVPFFT::write_micro_state_pvtu()
   out.open(filename,std::ofstream::binary);
 
   byte_offset = 0;
-  double crap = 0.0;
   
   //  Write Header
   //tmp_str = "<?xml version=\"1.0\"?>\n";
@@ -1033,7 +1041,7 @@ void EVPFFT::write_micro_state_pvtu()
       for (int kx = 1; kx <= npts1+1; kx++){
         ic += 1;
         for (int dim = 1; dim <= num_dims; dim++){
-          coord_tmp = xintp(dim,kx,ky,kz);
+          coord_tmp = xintp.host(dim,kx,ky,kz);
           out.write((char *) &coord_tmp,sizeof(coord_tmp));
         }
         pid(kx,ky,kz) = ic;
@@ -1101,7 +1109,6 @@ void EVPFFT::write_micro_state_pvtu()
   // for (int node_gid = 0; node_gid < num_points; node_gid++){
   //   for (int dim = 0; dim < num_dims; dim++){
   //     //out.write((char *) &node_vel.host(1, node_gid, dim),sizeof(double));
-  //     out.write((char *) &crap,sizeof(double));
   //   }
   // }
 
@@ -1264,7 +1271,6 @@ void EVPFFT::write_micro_state_pvtu()
   // for (int elem_gid = 0; elem_gid < num_cells; elem_gid++){
   //   for (int ii = 0; ii < 9; ii++){
   //     //out.write((char *) &elem_stress.host(1,elem_gid,ii),sizeof(double));
-  //     out.write((char *) &crap,sizeof(double));
   //   }
   // }
   // //  Density
diff --git a/src/Parallel-Solvers/FEA_Module.cpp b/src/Parallel-Solvers/FEA_Module.cpp
index 59f939db8..a1e610b73 100644
--- a/src/Parallel-Solvers/FEA_Module.cpp
+++ b/src/Parallel-Solvers/FEA_Module.cpp
@@ -99,9 +99,11 @@ FEA_Module::FEA_Module(Solver* Solver_Pointer)
     initial_node_coords_distributed         = Solver_Pointer->initial_node_coords_distributed;
     all_initial_node_coords_distributed     = Solver_Pointer->all_initial_node_coords_distributed;
     design_node_densities_distributed       = Solver_Pointer->design_node_densities_distributed;
+    design_node_coords_distributed          = Solver_Pointer->design_node_coords_distributed;
     filtered_node_densities_distributed     = Solver_Pointer->filtered_node_densities_distributed;
     test_node_densities_distributed         = Solver_Pointer->test_node_densities_distributed;
     all_node_densities_distributed          = Solver_Pointer->all_node_densities_distributed;
+    all_design_node_coords_distributed      = Solver_Pointer->all_design_node_coords_distributed;
     all_filtered_node_densities_distributed = Solver_Pointer->all_filtered_node_densities_distributed;
     Global_Element_Densities                = Solver_Pointer->Global_Element_Densities;
     Element_Types = Solver_Pointer->Element_Types;
diff --git a/src/Parallel-Solvers/FEA_Module.h b/src/Parallel-Solvers/FEA_Module.h
index defd5608a..95d453f8a 100644
--- a/src/Parallel-Solvers/FEA_Module.h
+++ b/src/Parallel-Solvers/FEA_Module.h
@@ -258,9 +258,11 @@ class FEA_Module
     Teuchos::RCP<MV>       initial_node_coords_distributed;
     Teuchos::RCP<MV>       all_initial_node_coords_distributed;
     Teuchos::RCP<MV>       design_node_densities_distributed;
+    Teuchos::RCP<MV>       design_node_coords_distributed;
     Teuchos::RCP<MV>       filtered_node_densities_distributed;
     Teuchos::RCP<const MV> test_node_densities_distributed;
     Teuchos::RCP<MV>       all_node_densities_distributed;
+    Teuchos::RCP<MV>       all_design_node_coords_distributed;
     Teuchos::RCP<MV>       all_filtered_node_densities_distributed;
     Teuchos::RCP<MV>       Global_Element_Densities;
 
diff --git a/src/Parallel-Solvers/FEA_Module_Inertial.cpp b/src/Parallel-Solvers/FEA_Module_Inertial.cpp
index 83e739da3..257929218 100644
--- a/src/Parallel-Solvers/FEA_Module_Inertial.cpp
+++ b/src/Parallel-Solvers/FEA_Module_Inertial.cpp
@@ -120,7 +120,269 @@ FEA_Module_Inertial::~FEA_Module_Inertial()
    Compute the mass of each element; estimated with quadrature
 ------------------------------------------------------------------------- */
 
-void FEA_Module_Inertial::compute_element_masses(const_host_vec_array design_densities, bool max_flag, bool use_initial_coords)
+void FEA_Module_Inertial::compute_element_masses(const_host_vec_array design_variables, bool max_flag, bool use_initial_coords)
+{
+    if(simparam->topology_optimization_on){
+        compute_element_masses_TO(design_variables, max_flag, use_initial_coords);
+    }
+    else if(simparam->shape_optimization_on){
+        compute_element_masses_SO(design_variables, max_flag, use_initial_coords);
+    }
+
+}
+
+/* ------------------------------------------------------------------------------------------
+   Compute the mass of each element when using shape optimization; estimated with quadrature
+--------------------------------------------------------------------------------------------- */
+
+void FEA_Module_Inertial::compute_element_masses_SO(const_host_vec_array design_coords, bool max_flag, bool use_initial_coords)
+{
+    // local number of uniquely assigned elements
+    size_t nonoverlap_nelements = element_map->getLocalNumElements();
+    // initialize memory for volume storage
+    host_vec_array Element_Masses = Global_Element_Masses->getLocalView<HostSpace>(Tpetra::Access::ReadWrite);
+    // local variable for host view in the dual view
+    const_host_vec_array all_node_coords;
+    if (use_initial_coords)
+    {
+        all_node_coords = all_initial_node_coords_distributed->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+    }
+    else
+    {
+        all_node_coords = all_design_node_coords_distributed->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+    }
+    const_host_vec_array all_design_densities;
+    if (nodal_density_flag)
+    {
+        all_design_densities = all_node_densities_distributed->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+    }
+    const_host_elem_conn_array nodes_in_elem  = global_nodes_in_elem_distributed->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+    int                        nodes_per_elem = elem->num_basis();
+    int                        z_quad, y_quad, x_quad, direct_product_count;
+    size_t                     local_node_id;
+    LO                         ielem;
+    GO                         global_element_index;
+
+    real_t Jacobian, current_density, weight_multiply;
+    // CArrayKokkos<real_t, array_layout, device_type, memory_traits> legendre_nodes_1D(num_gauss_points);
+    // CArrayKokkos<real_t, array_layout, device_type, memory_traits> legendre_weights_1D(num_gauss_points);
+    CArray<real_t>     legendre_nodes_1D(num_gauss_points);
+    CArray<real_t>     legendre_weights_1D(num_gauss_points);
+    
+    real_t pointer_quad_coordinate[num_dim];
+    real_t pointer_quad_coordinate_weight[num_dim];
+    real_t pointer_interpolated_point[num_dim];
+    real_t pointer_JT_row1[num_dim];
+    real_t pointer_JT_row2[num_dim];
+    real_t pointer_JT_row3[num_dim];
+    
+    ViewCArray<real_t> quad_coordinate(pointer_quad_coordinate, num_dim);
+    ViewCArray<real_t> quad_coordinate_weight(pointer_quad_coordinate_weight, num_dim);
+    ViewCArray<real_t> interpolated_point(pointer_interpolated_point, num_dim);
+    ViewCArray<real_t> JT_row1(pointer_JT_row1, num_dim);
+    ViewCArray<real_t> JT_row2(pointer_JT_row2, num_dim);
+    ViewCArray<real_t> JT_row3(pointer_JT_row3, num_dim);
+
+    real_t pointer_basis_values[elem->num_basis()];
+    real_t pointer_basis_derivative_s1[elem->num_basis()];
+    real_t pointer_basis_derivative_s2[elem->num_basis()];
+    real_t pointer_basis_derivative_s3[elem->num_basis()];
+    
+    ViewCArray<real_t> basis_values(pointer_basis_values, elem->num_basis());
+    ViewCArray<real_t> basis_derivative_s1(pointer_basis_derivative_s1, elem->num_basis());
+    ViewCArray<real_t> basis_derivative_s2(pointer_basis_derivative_s2, elem->num_basis());
+    ViewCArray<real_t> basis_derivative_s3(pointer_basis_derivative_s3, elem->num_basis());
+    
+    CArrayKokkos<real_t, array_layout, device_type, memory_traits> nodal_positions(elem->num_basis(), num_dim);
+    CArrayKokkos<real_t, array_layout, device_type, memory_traits> nodal_density(elem->num_basis());
+
+    // initialize weights
+    elements::legendre_nodes_1D(legendre_nodes_1D, num_gauss_points);
+    elements::legendre_weights_1D(legendre_weights_1D, num_gauss_points);
+
+    Solver::node_ordering_convention active_node_ordering_convention = Solver_Pointer_->active_node_ordering_convention;
+    
+    CArrayKokkos<size_t, array_layout, HostSpace, memory_traits> convert_node_order(max_nodes_per_element);
+    
+    if ((active_node_ordering_convention == Solver::ENSIGHT && num_dim == 3) || (active_node_ordering_convention == Solver::IJK && num_dim == 2))
+    {
+        convert_node_order(0) = 0;
+        convert_node_order(1) = 1;
+        convert_node_order(2) = 3;
+        convert_node_order(3) = 2;
+        if (num_dim == 3)
+        {
+            convert_node_order(4) = 4;
+            convert_node_order(5) = 5;
+            convert_node_order(6) = 7;
+            convert_node_order(7) = 6;
+        }
+    }
+    else
+    {
+        convert_node_order(0) = 0;
+        convert_node_order(1) = 1;
+        convert_node_order(2) = 2;
+        convert_node_order(3) = 3;
+        if (num_dim == 3)
+        {
+            convert_node_order(4) = 4;
+            convert_node_order(5) = 5;
+            convert_node_order(6) = 6;
+            convert_node_order(7) = 7;
+        }
+    }
+
+    // loop over elements and use quadrature rule to compute volume from Jacobian determinant
+    for (int nonoverlapping_ielem = 0; nonoverlapping_ielem < nonoverlap_nelements; nonoverlapping_ielem++)
+    {
+        global_element_index = element_map->getGlobalElement(nonoverlapping_ielem);
+        ielem = all_element_map->getLocalElement(global_element_index);
+        // acquire set of nodes for this local element
+        for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+        {
+            local_node_id = all_node_map->getLocalElement(nodes_in_elem(ielem, convert_node_order(node_loop)));
+            nodal_positions(node_loop, 0) = all_node_coords(local_node_id, 0);
+            nodal_positions(node_loop, 1) = all_node_coords(local_node_id, 1);
+            nodal_positions(node_loop, 2) = all_node_coords(local_node_id, 2);
+            if (nodal_density_flag)
+            {
+                nodal_density(node_loop) = all_design_densities(local_node_id, 0);
+            }
+            /*
+            if(myrank==1&&nodal_positions(node_loop,2)>10000000){
+                std::cout << " LOCAL MATRIX DEBUG ON TASK " << myrank << std::endl;
+                std::cout << node_loop+1 <<" " << local_node_id <<" "<< nodes_in_elem(ielem, node_loop) << " "<< nodal_positions(node_loop,2) << std::endl;
+                std::fflush(stdout);
+            }
+            */
+            // std::cout << local_node_id << " " << nodes_in_elem(ielem, node_loop) << " "
+            // << nodal_positions(node_loop,0) << " " << nodal_positions(node_loop,1) << " "<< nodal_positions(node_loop,2) << " " << nodal_density(node_loop) <<std::endl;
+        }
+
+        // debug print of index
+        // std::cout << "nonoverlap element id on TASK " << myrank << " is " << nonoverlapping_ielem << std::endl;
+        // std::fflush(stdout);
+
+        // initialize element mass
+        Element_Masses(nonoverlapping_ielem, 0) = 0;
+
+        direct_product_count = std::pow(num_gauss_points, num_dim);
+
+        // loop over quadrature points
+        for (int iquad = 0; iquad < direct_product_count; iquad++)
+        {
+            // set current quadrature point
+            if (num_dim == 3)
+            {
+                z_quad = iquad / (num_gauss_points * num_gauss_points);
+            }
+            y_quad = (iquad % (num_gauss_points * num_gauss_points)) / num_gauss_points;
+            x_quad = iquad % num_gauss_points;
+            quad_coordinate(0) = legendre_nodes_1D(x_quad);
+            quad_coordinate(1) = legendre_nodes_1D(y_quad);
+            if (num_dim == 3)
+            {
+                quad_coordinate(2) = legendre_nodes_1D(z_quad);
+            }
+
+            // set current quadrature weight
+            quad_coordinate_weight(0) = legendre_weights_1D(x_quad);
+            quad_coordinate_weight(1) = legendre_weights_1D(y_quad);
+            if (num_dim == 3)
+            {
+                quad_coordinate_weight(2) = legendre_weights_1D(z_quad);
+            }
+            else
+            {
+                quad_coordinate_weight(2) = 1;
+            }
+            weight_multiply = quad_coordinate_weight(0) * quad_coordinate_weight(1) * quad_coordinate_weight(2);
+
+            // compute shape functions at this point for the element type
+            elem->basis(basis_values, quad_coordinate);
+
+            // compute all the necessary coordinates and derivatives at this point
+
+            // compute shape function derivatives
+            elem->partial_xi_basis(basis_derivative_s1, quad_coordinate);
+            elem->partial_eta_basis(basis_derivative_s2, quad_coordinate);
+            elem->partial_mu_basis(basis_derivative_s3, quad_coordinate);
+
+            // compute derivatives of x,y,z w.r.t the s,t,w isoparametric space needed by JT (Transpose of the Jacobian)
+            // derivative of x,y,z w.r.t s
+            JT_row1(0) = 0;
+            JT_row1(1) = 0;
+            JT_row1(2) = 0;
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                JT_row1(0) += nodal_positions(node_loop, 0) * basis_derivative_s1(node_loop);
+                JT_row1(1) += nodal_positions(node_loop, 1) * basis_derivative_s1(node_loop);
+                JT_row1(2) += nodal_positions(node_loop, 2) * basis_derivative_s1(node_loop);
+            }
+
+            // derivative of x,y,z w.r.t t
+            JT_row2(0) = 0;
+            JT_row2(1) = 0;
+            JT_row2(2) = 0;
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                JT_row2(0) += nodal_positions(node_loop, 0) * basis_derivative_s2(node_loop);
+                JT_row2(1) += nodal_positions(node_loop, 1) * basis_derivative_s2(node_loop);
+                JT_row2(2) += nodal_positions(node_loop, 2) * basis_derivative_s2(node_loop);
+            }
+
+            // derivative of x,y,z w.r.t w
+            JT_row3(0) = 0;
+            JT_row3(1) = 0;
+            JT_row3(2) = 0;
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                JT_row3(0) += nodal_positions(node_loop, 0) * basis_derivative_s3(node_loop);
+                JT_row3(1) += nodal_positions(node_loop, 1) * basis_derivative_s3(node_loop);
+                JT_row3(2) += nodal_positions(node_loop, 2) * basis_derivative_s3(node_loop);
+                // debug print
+                /*if(myrank==1&&nodal_positions(node_loop,2)*basis_derivative_s3(node_loop)<-10000000){
+                std::cout << " ELEMENT VOLUME JACOBIAN DEBUG ON TASK " << myrank << std::endl;
+                std::cout << node_loop+1 << " " << JT_row3(2) << " "<< nodal_positions(node_loop,2) <<" "<< basis_derivative_s3(node_loop) << std::endl;
+                std::fflush(stdout);
+                }*/
+            }
+
+            // compute the determinant of the Jacobian
+            Jacobian = JT_row1(0) * (JT_row2(1) * JT_row3(2) - JT_row3(1) * JT_row2(2)) -
+                        JT_row1(1) * (JT_row2(0) * JT_row3(2) - JT_row3(0) * JT_row2(2)) +
+                        JT_row1(2) * (JT_row2(0) * JT_row3(1) - JT_row3(0) * JT_row2(1));
+            if (Jacobian < 0)
+            {
+                Jacobian = -Jacobian;
+            }
+
+            // compute density
+            current_density = 0;
+            if (max_flag)
+            {
+                current_density = 1;
+            }
+            else
+            {
+                for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+                {
+                    current_density += nodal_density(node_loop) * basis_values(node_loop);
+                }
+            }
+
+            Element_Masses(nonoverlapping_ielem, 0) += current_density * weight_multiply * Jacobian;
+        }
+    
+    }
+}
+
+/* ----------------------------------------------------------------------
+   Compute the mass of each element; estimated with quadrature
+------------------------------------------------------------------------- */
+
+void FEA_Module_Inertial::compute_element_masses_TO(const_host_vec_array design_densities, bool max_flag, bool use_initial_coords)
 {
     // local number of uniquely assigned elements
     size_t nonoverlap_nelements = element_map->getLocalNumElements();
@@ -385,7 +647,30 @@ void FEA_Module_Inertial::compute_element_masses(const_host_vec_array design_den
    Compute the gradients of mass function with respect to nodal densities
 ------------------------------------------------------------------------- */
 
-void FEA_Module_Inertial::compute_nodal_gradients(const_host_vec_array design_variables, host_vec_array design_gradients, bool use_initial_coords)
+void FEA_Module_Inertial::compute_nodal_gradients(const_host_vec_array design_variables, host_vec_array design_gradients, bool use_initial_coords){
+    if(simparam->topology_optimization_on){
+        compute_TO_gradients(design_variables, design_gradients, use_initial_coords);
+    }
+    else if(simparam->shape_optimization_on){
+        compute_shape_gradients(design_variables, design_gradients, use_initial_coords);
+    }
+
+}
+
+/* --------------------------------------------------------------------------------
+   Compute the gradients of mass function with respect to nodal design coordinates
+----------------------------------------------------------------------------------- */
+
+void FEA_Module_Inertial::compute_shape_gradients(const_host_vec_array design_coords, host_vec_array design_gradients, bool use_initial_coords)
+{
+
+}
+
+/* ----------------------------------------------------------------------
+   Compute the gradients of mass function with respect to nodal densities
+------------------------------------------------------------------------- */
+
+void FEA_Module_Inertial::compute_TO_gradients(const_host_vec_array design_variables, host_vec_array design_gradients, bool use_initial_coords)
 {
     // local number of uniquely assigned elements
     size_t nonoverlap_nelements = element_map->getLocalNumElements();
diff --git a/src/Parallel-Solvers/FEA_Module_Inertial.h b/src/Parallel-Solvers/FEA_Module_Inertial.h
index 252de74a5..666475449 100644
--- a/src/Parallel-Solvers/FEA_Module_Inertial.h
+++ b/src/Parallel-Solvers/FEA_Module_Inertial.h
@@ -49,13 +49,21 @@ class FEA_Module_Inertial : public FEA_Module
 
     void compute_element_volumes();
 
-    void compute_element_masses(const_host_vec_array design_densities, bool max_flag, bool use_initial_coords = false);
+    void compute_element_masses(const_host_vec_array design_variables, bool max_flag, bool use_initial_coords = false);
+    
+    void compute_element_masses_TO(const_host_vec_array design_densities, bool max_flag, bool use_initial_coords = false);
+    
+    void compute_element_masses_SO(const_host_vec_array design_coords, bool max_flag, bool use_initial_coords = false);
 
     void compute_element_moments(const_host_vec_array design_densities, bool max_flag, int moment_component, bool use_initial_coords = false);
 
     void compute_element_moments_of_inertia(const_host_vec_array design_densities, bool max_flag, int inertia_component, bool use_initial_coords = false);
 
-    void compute_nodal_gradients(const_host_vec_array design_densities, host_vec_array gradients, bool use_initial_coords = false);
+    void compute_nodal_gradients(const_host_vec_array design_variables, host_vec_array gradients, bool use_initial_coords = false);
+
+    void compute_TO_gradients(const_host_vec_array design_densities, host_vec_array gradients, bool use_initial_coords = false);
+
+    void compute_shape_gradients(const_host_vec_array design_coords, host_vec_array gradients, bool use_initial_coords = false);
 
     void compute_moment_gradients(const_host_vec_array design_densities, host_vec_array gradients, int moment_component, bool use_initial_coords = false);
 
diff --git a/src/Parallel-Solvers/Implicit-Lagrange/Topology_Optimization_Function_Headers.h b/src/Parallel-Solvers/Implicit-Lagrange/Implicit_Optimization_Function_Headers.h
similarity index 92%
rename from src/Parallel-Solvers/Implicit-Lagrange/Topology_Optimization_Function_Headers.h
rename to src/Parallel-Solvers/Implicit-Lagrange/Implicit_Optimization_Function_Headers.h
index 828628c1b..b285371be 100644
--- a/src/Parallel-Solvers/Implicit-Lagrange/Topology_Optimization_Function_Headers.h
+++ b/src/Parallel-Solvers/Implicit-Lagrange/Implicit_Optimization_Function_Headers.h
@@ -35,13 +35,14 @@
  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  **********************************************************************************************/
  
-#ifndef TOPOLOGY_OPTIMIZATION_FUNCTION_FILES_H
-#define TOPOLOGY_OPTIMIZATION_FUNCTION_FILES_H  
+#ifndef IMPLICIT_OPTIMIZATION_FUNCTION_FILES_H
+#define IMPLICIT_OPTIMIZATION_FUNCTION_FILES_H  
 
 #include "Mass_Objective.h"
-#include "Mass_Constraint.h"
+#include "Mass_Constraint_Top_Opt.h"
+#include "Mass_Constraint_Shape_Opt.h"
 #include "Center_of_Mass_Constraint.h"
-#include "Moment_of_Inertia_Constraint.h"
+#include "Moment_of_Inertia_Constraint_Top_Opt.h"
 #include "Bounded_Strain_Constraint.h"
 #include "Strain_Energy_Constraint.h"
 #include "Displacement_Constraint.h"
diff --git a/src/Parallel-Solvers/Implicit-Lagrange/Implicit_Solver.cpp b/src/Parallel-Solvers/Implicit-Lagrange/Implicit_Solver.cpp
index 72d7cad72..edfe5fc9b 100644
--- a/src/Parallel-Solvers/Implicit-Lagrange/Implicit_Solver.cpp
+++ b/src/Parallel-Solvers/Implicit-Lagrange/Implicit_Solver.cpp
@@ -82,7 +82,7 @@
 #include <ROL_TpetraMultiVector.hpp>
 
 //Objective Functions and Constraint Functions
-#include "Topology_Optimization_Function_Headers.h"
+#include "Implicit_Optimization_Function_Headers.h"
 
 
 #define BUFFER_LINES 20000
@@ -146,7 +146,7 @@ void Implicit_Solver::run(){
     MPI_Comm_size(world,&nranks);
     
     if(myrank == 0){
-      std::cout << "Running TO Solver" << std::endl;
+      std::cout << "Running Implicit Solver" << std::endl;
     }
     //initialize Trilinos communicator class
     comm = Tpetra::getDefaultComm();
@@ -219,7 +219,7 @@ void Implicit_Solver::run(){
     //set boundary conditions
     generate_tcs();
 
-    //initialize TO design variable storage
+    //initialize optmization design variable storage
     init_design();
 
     //process process list of requested FEA modules to construct list of objects
@@ -1135,14 +1135,14 @@ void Implicit_Solver::FEA_module_setup(){
 void Implicit_Solver::setup_optimization_problem(){
   int num_dim = simparam.num_dims;
   bool nodal_density_flag = simparam.nodal_density_flag;
-  int nTO_modules = simparam.TO_Module_List.size();
+  int nOpt_modules = simparam.Optimization_Module_List.size();
   int nmulti_objective_modules = simparam.Multi_Objective_Modules.size();
-  std::vector<TO_MODULE_TYPE> TO_Module_List = simparam.TO_Module_List;
-  std::vector<int> TO_Module_My_FEA_Module = simparam.TO_Module_My_FEA_Module;
+  std::vector<OPTIMIZATION_MODULE_TYPE> Optimization_Module_List = simparam.Optimization_Module_List;
+  std::vector<int> Optimization_Module_My_FEA_Module = simparam.Optimization_Module_My_FEA_Module;
   std::vector<int> Multi_Objective_Modules = simparam.Multi_Objective_Modules;
   std::vector<real_t> Multi_Objective_Weights = simparam.Multi_Objective_Weights;
   std::vector<std::vector<real_t>> Function_Arguments = simparam.Function_Arguments;
-  std::vector<FUNCTION_TYPE> TO_Function_Type = simparam.TO_Function_Type;
+  std::vector<FUNCTION_TYPE> Optimization_Function_Type = simparam.Optimization_Function_Type;
   std::vector<ROL::Ptr<ROL::Objective<real_t>>> Multi_Objective_Terms;
 
   std::string constraint_base, constraint_name;
@@ -1374,72 +1374,72 @@ void Implicit_Solver::setup_optimization_problem(){
   //Instantiate (the one) objective function for the problem
   ROL::Ptr<ROL::Objective<real_t>> obj, sub_obj;
   bool objective_declared = false;
-  for(int imodule = 0; imodule < nTO_modules; imodule++){
-    if(TO_Function_Type[imodule] == FUNCTION_TYPE::OBJECTIVE){
+  for(int imodule = 0; imodule < nOpt_modules; imodule++){
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::OBJECTIVE){
       //check if previous module already defined an objective, there must be one objective module
       if(objective_declared){
         *fos << "PROGRAM IS ENDING DUE TO ERROR; ANOTHER OBJECTIVE FUNCTION WITH NAME \"" 
-              << to_string(TO_Module_List[imodule]) <<"\" ATTEMPTED TO REPLACE A PREVIOUS OBJECTIVE; THERE MUST BE ONE OBJECTIVE." << std::endl;
+              << to_string(Optimization_Module_List[imodule]) <<"\" ATTEMPTED TO REPLACE A PREVIOUS OBJECTIVE; THERE MUST BE ONE OBJECTIVE." << std::endl;
           exit_solver(0);
       }
-      if(TO_Module_List[imodule] == TO_MODULE_TYPE::Strain_Energy_Minimize){
+      if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Strain_Energy_Minimize_TopOpt){
         //debug print
-        *fos << " STRAIN ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
+        *fos << " STRAIN ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
         if(simparam.optimization_options.method_of_moving_asymptotes){
-          sub_obj = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag);
+          sub_obj = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag);
           obj = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);
         }
         else{
-          obj = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag);
+          obj = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag);
         }
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Heat_Capacity_Potential_Minimize){
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Minimize_TopOpt){
         //debug print
-        *fos << " HEAT CAPACITY POTENTIAL OBJECTIVE EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
+        *fos << " HEAT CAPACITY POTENTIAL OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
         if(simparam.optimization_options.method_of_moving_asymptotes){
-          sub_obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag);
+          sub_obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag);
           obj = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);
         }
         else{
-          obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag);
+          obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag);
         }
       }
       //Multi-Objective case
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Multi_Objective){
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Multi_Objective){
         //allocate vector of Objective Functions to pass
         Multi_Objective_Terms = std::vector<ROL::Ptr<ROL::Objective<real_t>>>(nmulti_objective_modules);
         for(int imulti = 0; imulti < nmulti_objective_modules; imulti++){
           //get module index for objective term
           module_id = Multi_Objective_Modules[imulti];
-          if(TO_Module_List[module_id] == TO_MODULE_TYPE::Strain_Energy_Minimize){
+          if(Optimization_Module_List[module_id] == OPTIMIZATION_MODULE_TYPE::Strain_Energy_Minimize_TopOpt){
             //debug print
-            *fos << " STRAIN ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[module_id] << std::endl;
+            *fos << " STRAIN ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[module_id] << std::endl;
             if(simparam.optimization_options.method_of_moving_asymptotes){
-              sub_obj = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[module_id]], nodal_density_flag);
+              sub_obj = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[module_id]], nodal_density_flag);
               Multi_Objective_Terms[imulti] = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);
             }
             else{
-              Multi_Objective_Terms[imulti] = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[module_id]], nodal_density_flag);
+              Multi_Objective_Terms[imulti] = ROL::makePtr<StrainEnergyMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[module_id]], nodal_density_flag);
             }
           }
-          else if(TO_Module_List[module_id] == TO_MODULE_TYPE::Heat_Capacity_Potential_Minimize){
+          else if(Optimization_Module_List[module_id] == OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Minimize_TopOpt){
             //debug print
-            *fos << " HEAT CAPACITY POTENTIAL OBJECTIVE EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[module_id] << std::endl;
+            *fos << " HEAT CAPACITY POTENTIAL OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[module_id] << std::endl;
             if(simparam.optimization_options.method_of_moving_asymptotes){
-              sub_obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[module_id]], nodal_density_flag);
+              sub_obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[module_id]], nodal_density_flag);
               Multi_Objective_Terms[imulti] = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);}
             else{
-              Multi_Objective_Terms[imulti] = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[module_id]], nodal_density_flag);
+              Multi_Objective_Terms[imulti] = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[module_id]], nodal_density_flag);
             }
           }
-          else if(TO_Module_List[module_id] == TO_MODULE_TYPE::Thermo_Elastic_Strain_Energy_Minimize){
+          else if(Optimization_Module_List[module_id] == OPTIMIZATION_MODULE_TYPE::Thermo_Elastic_Strain_Energy_Minimize_TopOpt){
             //debug print
-            *fos << " HEAT CAPACITY POTENTIAL OBJECTIVE EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[module_id] << std::endl;
+            *fos << " HEAT CAPACITY POTENTIAL OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[module_id] << std::endl;
             if(simparam.optimization_options.method_of_moving_asymptotes){
-              sub_obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[module_id]], nodal_density_flag);
+              sub_obj = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[module_id]], nodal_density_flag);
               Multi_Objective_Terms[imulti] = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);}
             else{
-              Multi_Objective_Terms[imulti] = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[TO_Module_My_FEA_Module[module_id]], nodal_density_flag);
+              Multi_Objective_Terms[imulti] = ROL::makePtr<HeatCapacityPotentialMinimize_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[module_id]], nodal_density_flag);
             }
           }
         }
@@ -1448,7 +1448,7 @@ void Implicit_Solver::setup_optimization_problem(){
       }
       else{
         *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED OBJECTIVE FUNCTION REQUESTED WITH NAME \"" 
-              << TO_Module_List[imodule] <<"\"" << std::endl;
+              << Optimization_Module_List[imodule] <<"\"" << std::endl;
         exit_solver(0);
       }
       objective_declared = true;
@@ -1479,21 +1479,21 @@ void Implicit_Solver::setup_optimization_problem(){
   //problem->addConstraint("equality Constraint 3",eq_constraint3,constraint_mul3);
   //problem->addLinearConstraint("Equality Constraint",eq_constraint,constraint_mul);
   
-  for(int imodule = 0; imodule < nTO_modules; imodule++){
+  for(int imodule = 0; imodule < nOpt_modules; imodule++){
     number_union.str(constraint_base);
     number_union << imodule + 1;
     constraint_name = number_union.str();
     ROL::Ptr<std::vector<real_t> > li_ptr = ROL::makePtr<std::vector<real_t>>(1,0.0);
     ROL::Ptr<ROL::Vector<real_t> > constraint_mul = ROL::makePtr<ROL::StdVector<real_t>>(li_ptr);
-    if(TO_Function_Type[imodule] == FUNCTION_TYPE::EQUALITY_CONSTRAINT){
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::EQUALITY_CONSTRAINT){
       //pointers are reference counting
       ROL::Ptr<ROL::Constraint<real_t>> eq_constraint;
-      if(TO_Module_List[imodule] == TO_MODULE_TYPE::Mass_Constraint){
+      if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Mass_Constraint_TopOpt){
         
-        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false);
+        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Displacement_Constraint){
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Displacement_Constraint_TopOpt){
         //simple test of assignment to the vector for constraint dofs
         Teuchos::RCP<MV> target_displacements = Teuchos::rcp(new MV(local_dof_map, 1));
         Teuchos::RCP<Tpetra::MultiVector<int,LO,GO>> active_dofs = Teuchos::rcp(new Tpetra::MultiVector<int,LO,GO>(local_dof_map, 1));
@@ -1525,8 +1525,8 @@ void Implicit_Solver::setup_optimization_problem(){
           }
           disp_in->close(); // Close the file
         }
-        *fos << " DISPLACEMENT CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<DisplacementConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, target_displacements, active_dofs, 0, false);
+        *fos << " DISPLACEMENT CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<DisplacementConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, target_displacements, active_dofs, 0, false);
 
         //ROL::Ptr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>> rol_x =
         //ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>>(design_node_densities_distributed);
@@ -1543,31 +1543,31 @@ void Implicit_Solver::setup_optimization_problem(){
         //eq_constraint->checkApplyJacobian(*rol_x, *rol_d, *constraint_buf);
         //eq_constraint->checkApplyAdjointHessian(*rol_x, *constraint_buf, *rol_d, *rol_x);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Moment_of_Inertia_Constraint){
-        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], false);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt){
+        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], false);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Center_of_Mass_Constraint){
-        *fos << " CENTER OF MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<CenterOfMassConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], false);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Center_of_Mass_Constraint_TopOpt){
+        *fos << " CENTER OF MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<CenterOfMassConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], false);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Strain_Energy_Constraint){    
-        *fos << " STRAIN ENERGY CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<StrainEnergyConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Strain_Energy_Constraint_TopOpt){    
+        *fos << " STRAIN ENERGY CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<StrainEnergyConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Heat_Capacity_Potential_Constraint){
-        *fos << " HEAT CAPACITY POTENTIAL CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<HeatCapacityPotentialConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Constraint_TopOpt){
+        *fos << " HEAT CAPACITY POTENTIAL CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<HeatCapacityPotentialConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false);
       }
       else{
-        *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED EQUALITY CONSTRAINT FUNCTION REQUESTED WITH NAME \"" <<TO_Module_List[imodule] <<"\"" << std::endl;
+        *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED EQUALITY CONSTRAINT FUNCTION REQUESTED WITH NAME \"" <<Optimization_Module_List[imodule] <<"\"" << std::endl;
         exit_solver(0);
       }
       *fos << " ADDING CONSTRAINT " << constraint_name << std::endl;
       problem->addConstraint(constraint_name, eq_constraint, constraint_mul);
     }
 
-    if(TO_Function_Type[imodule] == FUNCTION_TYPE::INEQUALITY_CONSTRAINT){
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::INEQUALITY_CONSTRAINT){
       //pointers are reference counting
       ROL::Ptr<ROL::Constraint<real_t>> ineq_constraint;
       ROL::Ptr<std::vector<real_t> > ll_ptr = ROL::makePtr<std::vector<real_t>>(1,Function_Arguments[imodule][0]);
@@ -1575,39 +1575,39 @@ void Implicit_Solver::setup_optimization_problem(){
       ROL::Ptr<ROL::Vector<real_t> > ll = ROL::makePtr<ROL::StdVector<real_t>>(ll_ptr);
       ROL::Ptr<ROL::Vector<real_t> > lu = ROL::makePtr<ROL::StdVector<real_t>>(lu_ptr);
       ROL::Ptr<ROL::BoundConstraint<real_t>> constraint_bnd = ROL::makePtr<ROL::Bounds<real_t>>(ll,lu);
-      if(TO_Module_List[imodule] == TO_MODULE_TYPE::Mass_Constraint){
-        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag);
+      if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Mass_Constraint_TopOpt){
+        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Displacement_Constraint){
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Displacement_Constraint_TopOpt){
         //simple test of assignment to the vector for constraint dofs
         Teuchos::RCP<MV> target_displacements = Teuchos::rcp(new MV(local_dof_map, 1));
         Teuchos::RCP<Tpetra::MultiVector<int,LO,GO>> active_dofs = Teuchos::rcp(new Tpetra::MultiVector<int,LO,GO>(local_dof_map, 1));
         active_dofs->putScalar(0);
         host_vec_array target_displacements_view = target_displacements->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
         host_ivec_array active_dofs_view = active_dofs->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
-        *fos << " DISPLACEMENT CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<DisplacementConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, target_displacements, active_dofs);
+        *fos << " DISPLACEMENT CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<DisplacementConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, target_displacements, active_dofs);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Moment_of_Inertia_Constraint){
-        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0]);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt){
+        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0]);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Center_of_Mass_Constraint){
-        *fos << " CENTER OF MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<CenterOfMassConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0]);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Center_of_Mass_Constraint_TopOpt){
+        *fos << " CENTER OF MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<CenterOfMassConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0]);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Strain_Energy_Constraint){
-        *fos << " STRAIN ENERGY CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<StrainEnergyConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][2]);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Strain_Energy_Constraint_TopOpt){
+        *fos << " STRAIN ENERGY CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<StrainEnergyConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][2]);
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Heat_Capacity_Potential_Constraint){
-        *fos << " HEAT CAPACITY POTENTIAL CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<HeatCapacityPotentialConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][2]);
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Constraint_TopOpt){
+        *fos << " HEAT CAPACITY POTENTIAL CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<HeatCapacityPotentialConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][2]);
       }
       else{
         *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED INEQUALITY CONSTRAINT FUNCTION REQUESTED WITH NAME \"" 
-              << to_string(TO_Module_List[imodule]) <<"\"" << std::endl;
+              << to_string(Optimization_Module_List[imodule]) <<"\"" << std::endl;
         exit_solver(0);
       }
       *fos << " ADDING CONSTRAINT " << constraint_name << std::endl;
diff --git a/src/Parallel-Solvers/Optimization_Modules/Mass_Constraint_Shape_Opt.h b/src/Parallel-Solvers/Optimization_Modules/Mass_Constraint_Shape_Opt.h
new file mode 100644
index 000000000..349214af4
--- /dev/null
+++ b/src/Parallel-Solvers/Optimization_Modules/Mass_Constraint_Shape_Opt.h
@@ -0,0 +1,323 @@
+/**********************************************************************************************
+ © 2020. Triad National Security, LLC. All rights reserved.
+ This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
+ National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
+ Department of Energy/National Nuclear Security Administration. All rights in the program are
+ reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear
+ Security Administration. The Government is granted for itself and others acting on its behalf a
+ nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare
+ derivative works, distribute copies to the public, perform publicly and display publicly, and
+ to permit others to do so.
+ This program is open source under the BSD-3 License.
+ Redistribution and use in source and binary forms, with or without modification, are permitted
+ provided that the following conditions are met:
+ 
+ 1.  Redistributions of source code must retain the above copyright notice, this list of
+ conditions and the following disclaimer.
+ 
+ 2.  Redistributions in binary form must reproduce the above copyright notice, this list of
+ conditions and the following disclaimer in the documentation and/or other materials
+ provided with the distribution.
+ 
+ 3.  Neither the name of the copyright holder nor the names of its contributors may be used
+ to endorse or promote products derived from this software without specific prior
+ written permission.
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
+ CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+ OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+ OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ **********************************************************************************************/
+ 
+#ifndef MASS_CONSTRAINT_SHAPEOPT_H
+#define MASS_CONSTRAINT_SHAPEOPT_H
+
+#include "matar.h"
+#include "elements.h"
+#include <Teuchos_ScalarTraits.hpp>
+#include <Teuchos_RCP.hpp>
+#include <Teuchos_oblackholestream.hpp>
+#include <Teuchos_Tuple.hpp>
+#include <Teuchos_VerboseObject.hpp>
+
+#include <Tpetra_Core.hpp>
+#include <Tpetra_Map.hpp>
+#include <Tpetra_MultiVector.hpp>
+#include "Tpetra_Details_makeColMap.hpp"
+#include "Tpetra_Details_DefaultTypes.hpp"
+
+#include "ROL_Types.hpp"
+#include <ROL_TpetraMultiVector.hpp>
+#include "ROL_Constraint.hpp"
+#include "ROL_Elementwise_Reduce.hpp"
+#include "FEA_Module_Inertial.h"
+
+class MassConstraint_ShapeOpt : public ROL::Constraint<real_t> {
+  
+  typedef Tpetra::Map<>::local_ordinal_type LO;
+  typedef Tpetra::Map<>::global_ordinal_type GO;
+  typedef Tpetra::Map<>::node_type Node;
+  typedef Tpetra::Map<LO, GO, Node> Map;
+  typedef Tpetra::MultiVector<real_t, LO, GO, Node> MV;
+  typedef ROL::Vector<real_t> V;
+  typedef ROL::TpetraMultiVector<real_t,LO,GO,Node> ROL_MV;
+  
+  using traits = Kokkos::ViewTraits<LO*, Kokkos::LayoutLeft, void, void>;
+  using array_layout    = typename traits::array_layout;
+  using execution_space = typename traits::execution_space;
+  using device_type     = typename traits::device_type;
+  using memory_traits   = typename traits::memory_traits;
+  using global_size_t = Tpetra::global_size_t;
+
+  typedef Kokkos::View<real_t*, Kokkos::LayoutRight, device_type, memory_traits> values_array;
+  typedef Kokkos::View<GO*, array_layout, device_type, memory_traits> global_indices_array;
+  typedef Kokkos::View<LO*, array_layout, device_type, memory_traits> indices_array;
+  
+  //typedef Kokkos::DualView<real_t**, Kokkos::LayoutLeft, device_type>::t_dev vec_array;
+  typedef MV::dual_view_type::t_dev vec_array;
+  typedef MV::dual_view_type::t_host host_vec_array;
+  typedef Kokkos::View<const real_t**, array_layout, HostSpace, memory_traits> const_host_vec_array;
+  typedef MV::dual_view_type dual_vec_array;
+
+private:
+
+  FEA_Module_Inertial *FEM_;
+  ROL::Ptr<ROL_MV> ROL_Element_Masses;
+  ROL::Ptr<ROL_MV> ROL_Gradients;
+  Teuchos::RCP<MV> constraint_gradients_distributed;
+  real_t initial_mass;
+  real_t current_mass;
+  bool inequality_flag_, use_initial_coords_;
+  real_t constraint_value_;
+
+  ROL::Ptr<const MV> getVector( const V& x ) {
+    return dynamic_cast<const ROL_MV&>(x).getVector();
+  }
+
+  ROL::Ptr<MV> getVector( V& x ) {
+    return dynamic_cast<ROL_MV&>(x).getVector();
+  }
+
+public:
+  bool nodal_density_flag_;
+  int last_comm_step, current_step, last_solve_step;
+  std::string my_fea_module = "Inertial";
+
+  MassConstraint_ShapeOpt(FEA_Module *FEM, real_t constraint_value=0, bool inequality_flag=true, bool use_initial_coords=false) 
+  {
+    FEM_ = dynamic_cast<FEA_Module_Inertial*>(FEM);
+    use_initial_coords_ = use_initial_coords;
+    last_comm_step = last_solve_step = -1;
+    current_step = 0;
+    inequality_flag_ = inequality_flag;
+    constraint_value_ = constraint_value;
+    ROL_Element_Masses = ROL::makePtr<ROL_MV>(FEM_->Global_Element_Masses);
+    const_host_vec_array design_coords = FEM_->design_node_coords_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadOnly);
+    
+    FEM_->compute_element_masses(design_coords,true,use_initial_coords_);
+    FEM_->mass_init = true;
+    
+    //sum per element results across all MPI ranks
+    ROL::Elementwise::ReductionSum<real_t> sumreduc;
+    FEM_->mass = initial_mass = ROL_Element_Masses->reduce(sumreduc);
+    //debug print
+    if(FEM_->myrank==0)
+      std::cout << "INITIAL MASS: " << initial_mass << std::endl;
+    if(FEM_->mass_gradients_distributed.is_null())
+      FEM_->mass_gradients_distributed = Teuchos::rcp(new MV(FEM_->map, 1));
+    constraint_gradients_distributed = FEM_->mass_gradients_distributed;
+  }
+
+  /* --------------------------------------------------------------------------------------
+   Update solver state variables to synchronize with the current design variable vector, z
+  ----------------------------------------------------------------------------------------- */
+
+  void update(const ROL::Vector<real_t> &z, ROL::UpdateType type, int iter = -1 ) {
+    current_step++;
+
+    ROL::Ptr<const MV> zp = getVector(z);
+    const_host_vec_array design_coords = zp->getLocalView<HostSpace> (Tpetra::Access::ReadOnly);
+
+    if (type == ROL::UpdateType::Initial)  {
+
+      //initial design density data was already communicated for ghost nodes in init_design()
+    }
+    else if (type == ROL::UpdateType::Accept) {
+      // z was accepted as the new iterate
+    }
+    else if (type == ROL::UpdateType::Revert) {
+      // z has been rejected as the new iterate
+      // Revert to cached value
+      FEM_->comm_variables(zp);
+    }
+    else if (type == ROL::UpdateType::Trial) {
+      // This is a new value of x
+      FEM_->comm_variables(zp);
+    }
+    else { // ROL::UpdateType::Temp
+      // This is a new value of x used for,
+      // e.g., finite-difference checks
+      FEM_->comm_variables(zp);
+    }
+  }
+
+  /* --------------------------------------------------------------------------------------
+   Update constraint value (c) with the current design variable vector, z
+  ----------------------------------------------------------------------------------------- */
+
+  void value(ROL::Vector<real_t> &c, const ROL::Vector<real_t> &z, real_t &tol ) override {
+    //std::cout << "Started constraint value on task " <<FEM_->myrank <<std::endl;
+    ROL::Ptr<const MV> zp = getVector(z);
+    ROL::Ptr<std::vector<real_t>> cp = dynamic_cast<ROL::StdVector<real_t>&>(c).getVector();
+    const_host_vec_array design_coords = zp->getLocalView<HostSpace> (Tpetra::Access::ReadOnly);
+
+    //communicate ghosts and solve for nodal degrees of freedom as a function of the current design variables
+    /*
+    if(last_comm_step!=current_step){
+      FEM_->comm_variables(zp);
+      last_comm_step = current_step;
+    }
+    */
+    FEM_->compute_element_masses(design_coords,false,use_initial_coords_);
+    
+    //sum per element results across all MPI ranks
+    ROL::Elementwise::ReductionSum<real_t> sumreduc;
+    current_mass = ROL_Element_Masses->reduce(sumreduc);
+    FEM_->mass = current_mass;
+    FEM_->mass_update = current_step;
+    
+    //debug print
+    if(FEM_->myrank==0)
+    std::cout << "SYSTEM MASS RATIO: " << current_mass/initial_mass << std::endl;
+    
+    if(inequality_flag_)
+      (*cp)[0] = current_mass/initial_mass;
+    else
+      (*cp)[0] = current_mass/initial_mass - constraint_value_;
+
+    //std::cout << "Ended constraint value on task " <<FEM_->myrank <<std::endl;
+  }
+
+  /* ----------------------------------------------------------------------------------------
+    Update adjoint of the constraint jacobian (ajv), involves gradient vector
+    and constraint adjoint (v) for the current design variable vector, z
+  ------------------------------------------------------------------------------------------*/
+
+  void applyAdjointJacobian(ROL::Vector<real_t> &ajv, const ROL::Vector<real_t> &v, const ROL::Vector<real_t> &z, real_t &tol) override {
+    //std::cout << "Started constraint adjoint grad on task " <<FEM_->myrank << std::endl;
+     //get Tpetra multivector pointer from the ROL vector
+    ROL::Ptr<const MV> zp = getVector(z);
+    ROL::Ptr<const std::vector<real_t>> vp = dynamic_cast<const ROL::StdVector<real_t>&>(v).getVector();
+    ROL::Ptr<MV> ajvp = getVector(ajv);
+    int num_dim = FEM_->num_dim;
+    
+    //ROL::Ptr<ROL_MV> ROL_Element_Volumes;
+
+    //get local view of the data
+    const_host_vec_array design_coords = zp->getLocalView<HostSpace> (Tpetra::Access::ReadOnly);
+    //host_vec_array constraint_gradients = constraint_gradients_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+    host_vec_array constraint_gradients = ajvp->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+    //host_vec_array dual_constraint_vector = vp->getLocalView<HostSpace> (Tpetra::Access::ReadOnly);
+
+    //communicate ghosts
+    /*
+    if(last_comm_step!=current_step){
+      FEM_->comm_variables(zp);
+      last_comm_step = current_step;
+    }
+    */
+    int rnum_elem = FEM_->rnum_elem;
+
+    FEM_->compute_nodal_gradients(design_coords, constraint_gradients, use_initial_coords_);
+      //debug print of gradient
+      //std::ostream &out = std::cout;
+      //Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
+      //if(FEM_->myrank==0)
+      //*fos << "Gradient data :" << std::endl;
+      //ajvp->describe(*fos,Teuchos::VERB_EXTREME);
+      //*fos << std::endl;
+      //std::fflush(stdout);
+    for(int i = 0; i < FEM_->nlocal_nodes; i++){
+      constraint_gradients(i,0) *= (*vp)[0]/initial_mass;
+      constraint_gradients(i,1) *= (*vp)[0]/initial_mass;
+      if(num_dim==3){
+        constraint_gradients(i,2) *= (*vp)[0]/initial_mass;
+      }
+    }
+    
+    
+    //std::cout << "Ended constraint adjoint grad on task " <<FEM_->myrank  << std::endl;
+    //debug print
+    //std::cout << "Constraint Gradient value " << std::endl;
+  }
+  
+  /* ----------------------------------------------------------------------------------------
+    Update the constraint jacobian (jv), involves gradient vector
+    and design vector differential (v) for the current design variable vector, x
+  ------------------------------------------------------------------------------------------*/
+
+  void applyJacobian(ROL::Vector<real_t> &jv, const ROL::Vector<real_t> &v, const ROL::Vector<real_t> &x, real_t &tol) override {
+    //std::cout << "Started constraint grad on task " <<FEM_->myrank  << std::endl;
+     //get Tpetra multivector pointer from the ROL vector
+    ROL::Ptr<const MV> zp = getVector(x);
+    ROL::Ptr<std::vector<real_t>> jvp = dynamic_cast<ROL::StdVector<real_t>&>(jv).getVector();
+    int num_dim = FEM_->num_dim;
+    
+    //ROL::Ptr<ROL_MV> ROL_Element_Volumes;
+
+    //get local view of the data
+    const_host_vec_array design_coords = zp->getLocalView<HostSpace> (Tpetra::Access::ReadOnly);
+    host_vec_array constraint_gradients = constraint_gradients_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+
+    //communicate ghosts and solve for nodal degrees of freedom as a function of the current design variables
+    //communicate ghosts
+    /*
+    if(last_comm_step!=current_step){
+      FEM_->comm_variables(zp);
+      last_comm_step = current_step;
+    }
+    */
+
+    FEM_->compute_nodal_gradients(design_coords, constraint_gradients);
+    for(int i = 0; i < FEM_->nlocal_nodes; i++){
+      constraint_gradients(i,0) /= initial_mass;
+      constraint_gradients(i,1) /= initial_mass;
+      if(num_dim==3){
+        constraint_gradients(i,2) /= initial_mass;
+      }
+    }
+    
+
+    ROL_Gradients = ROL::makePtr<ROL_MV>(constraint_gradients_distributed);
+    real_t gradient_dot_v = ROL_Gradients->dot(v);
+    //debug print
+    //std::cout << "Constraint Gradient value " << gradient_dot_v << std::endl;
+
+    (*jvp)[0] = gradient_dot_v;
+    //std::cout << "Ended constraint grad on task " <<FEM_->myrank  << std::endl;
+  }
+  
+  /* ----------------------------------------------------------------------------------------
+    Update adjoint of the constraint Hessian vector product (ahuv), product of hessian and
+    differential design vector (v), and constraint adjoint (u) for the 
+    current design variable vector, z
+  ------------------------------------------------------------------------------------------*/
+
+  void applyAdjointHessian(ROL::Vector<real_t> &ahuv, const ROL::Vector<real_t> &u, const ROL::Vector<real_t> &v, const ROL::Vector<real_t> &z, real_t &tol) {
+    
+    // Unwrap hv
+    ROL::Ptr<MV> ahuvp = getVector(ahuv);
+
+    ahuvp->putScalar(0);
+    
+  }
+
+};
+
+#endif // end header guard
diff --git a/src/Parallel-Solvers/Optimization_Modules/Mass_Constraint.h b/src/Parallel-Solvers/Optimization_Modules/Mass_Constraint_Top_Opt.h
similarity index 90%
rename from src/Parallel-Solvers/Optimization_Modules/Mass_Constraint.h
rename to src/Parallel-Solvers/Optimization_Modules/Mass_Constraint_Top_Opt.h
index 4145e58b7..8ccf14513 100644
--- a/src/Parallel-Solvers/Optimization_Modules/Mass_Constraint.h
+++ b/src/Parallel-Solvers/Optimization_Modules/Mass_Constraint_Top_Opt.h
@@ -329,63 +329,7 @@ class MassConstraint_TopOpt : public ROL::Constraint<real_t> {
     ahuvp->putScalar(0);
     
   }
-  /*
-  void hessVec_21( ROL::Vector<real_t> &hv, const ROL::Vector<real_t> &v, 
-                   const ROL::Vector<real_t> &u, const ROL::Vector<real_t> &z, real_t &tol ) {
-                     
-    // Unwrap g
-    ROL::Ptr<MV> hvp = getVector(hv);
-
-    // Unwrap v
-    ROL::Ptr<const MV> vp = getVector(v);
-
-    // Unwrap x
-    ROL::Ptr<const MV> up = getVector(u);
-    ROL::Ptr<const MV> zp = getVector(z);
- 
-    // Apply Jacobian
-    hv.zero();
-    if ( !useLC_ ) {
-      std::MV<real_t> U;
-      U.assign(up->begin(),up->end());
-      FEM_->set_boundary_conditions(U);
-      std::MV<real_t> V;
-      V.assign(vp->begin(),vp->end());
-      FEM_->set_boundary_conditions(V);
-      FEM_->apply_adjoint_jacobian(*hvp,U,*zp,V);
-      for (size_t i=0; i<hvp->size(); i++) {
-        (*hvp)[i] *= 2.0;
-      }
-    }
-    
-  }
-
-  void hessVec_22( ROL::Vector<real_t> &hv, const ROL::Vector<real_t> &v, 
-                   const ROL::Vector<real_t> &u, const ROL::Vector<real_t> &z, real_t &tol ) {
-                     
-    ROL::Ptr<MV> hvp = getVector(hv);
 
-    // Unwrap v
-    ROL::Ptr<const MV> vp = getVector(v);
-
-    // Unwrap x
-    ROL::Ptr<const MV> up = getVector(u);
-    ROL::Ptr<const MV> zp = getVector(z);
-    
-    // Apply Jacobian
-    hv.zero();
-    if ( !useLC_ ) {
-      MV U;
-      U.assign(up->begin(),up->end());
-      FEM_->set_boundary_conditions(U);
-      MV V;
-      V.assign(vp->begin(),vp->end());
-      FEM_->set_boundary_conditions(V);
-      FEM_->apply_adjoint_jacobian(*hvp,U,*zp,*vp,U);
-    }
-    
-  }
-  */
 };
 
 #endif // end header guard
diff --git a/src/Parallel-Solvers/Optimization_Modules/Moment_of_Inertia_Constraint.h b/src/Parallel-Solvers/Optimization_Modules/Moment_of_Inertia_Constraint_Top_Opt.h
similarity index 100%
rename from src/Parallel-Solvers/Optimization_Modules/Moment_of_Inertia_Constraint.h
rename to src/Parallel-Solvers/Optimization_Modules/Moment_of_Inertia_Constraint_Top_Opt.h
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/FEA_Module_Dynamic_Elasticity.cpp b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/FEA_Module_Dynamic_Elasticity.cpp
index ea7b1bcaa..fc76ca3cf 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/FEA_Module_Dynamic_Elasticity.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/FEA_Module_Dynamic_Elasticity.cpp
@@ -1738,29 +1738,16 @@ void FEA_Module_Dynamic_Elasticity::elastic_solve()
 
     const DCArrayKokkos<boundary_t> boundary = module_params->boundary;
     const DCArrayKokkos<material_t> material = simparam->material;
-
-    int nTO_modules;
     int old_max_forward_buffer;
 
     unsigned long cycle;
 
     real_t objective_accumulation, global_objective_accumulation;
 
-    std::vector<std::vector<int>> FEA_Module_My_TO_Modules = simparam->FEA_Module_My_TO_Modules;
     problem = Explicit_Solver_Pointer_->problem; // Pointer to ROL optimization problem object
     ROL::Ptr<ROL::Objective<real_t>> obj_pointer;
 
     // reset time accumulating objective and constraints
-    /*
-    for(int imodule = 0 ; imodule < FEA_Module_My_TO_Modules[my_fea_module_index_].size(); imodule++){
-    current_module_index = FEA_Module_My_TO_Modules[my_fea_module_index_][imodule];
-    //test if module needs reset
-    if(){
-
-    }
-    }
-    */
-    // simple setup to just request KE for now; above loop to be expanded and used later for scanning modules
     if (simparam->topology_optimization_on) {
         obj_pointer = problem->getObjective();
         KineticEnergyMinimize_TopOpt& kinetic_energy_minimize_function = dynamic_cast<KineticEnergyMinimize_TopOpt&>(*obj_pointer);
@@ -1784,10 +1771,6 @@ void FEA_Module_Dynamic_Elasticity::elastic_solve()
         }
     }
 
-    if (simparam->topology_optimization_on) {
-        nTO_modules = simparam->TO_Module_List.size();
-    }
-
     int myrank = Explicit_Solver_Pointer_->myrank;
     if (simparam->output_options.output_file_format == OUTPUT_FORMAT::vtk && simparam->output_options.write_initial) {
         if (myrank == 0) {
@@ -1830,7 +1813,7 @@ void FEA_Module_Dynamic_Elasticity::elastic_solve()
     int nlocal_elem_non_overlapping = Explicit_Solver_Pointer_->nlocal_elem_non_overlapping;
 
     // extensive IE
-    REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
         IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     }, IE_sum);
     IE_t0 = IE_sum;
@@ -1838,7 +1821,7 @@ void FEA_Module_Dynamic_Elasticity::elastic_solve()
     MPI_Allreduce(&IE_t0, &global_IE_t0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // extensive KE
-    REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
         double ke = 0;
         for (size_t dim = 0; dim < num_dim; dim++) {
             ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
@@ -2332,7 +2315,7 @@ void FEA_Module_Dynamic_Elasticity::elastic_solve()
                 KE_loc_sum = 0.0;
                 KE_sum     = 0.0;
                 // extensive KE
-                REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+                FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
                     double ke = 0;
                     for (size_t dim = 0; dim < num_dim; dim++) {
                         // midpoint integration approximation
@@ -2429,7 +2412,7 @@ void FEA_Module_Dynamic_Elasticity::elastic_solve()
     KE_sum     = 0.0;
 
     // extensive IE
-    REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
         IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     }, IE_sum);
     IE_tend = IE_sum;
@@ -2438,7 +2421,7 @@ void FEA_Module_Dynamic_Elasticity::elastic_solve()
     MPI_Allreduce(&IE_tend, &global_IE_tend, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // extensive KE
-    REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
         double ke = 0;
         for (size_t dim = 0; dim < num_dim; dim++) {
             ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/elastic_optimization.cpp b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/elastic_optimization.cpp
index 6f45205dd..cffa586e4 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/elastic_optimization.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/elastic_optimization.cpp
@@ -1,5 +1,5 @@
 /**********************************************************************************************
- © 2020. Triad National Security, LLC. All rights reserved.
+ � 2020. Triad National Security, LLC. All rights reserved.
  This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
  National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
  Department of Energy/National Nuclear Security Administration. All rights in the program are
@@ -106,7 +106,6 @@ void FEA_Module_Dynamic_Elasticity::update_forward_solve(Teuchos::RCP<const MV>
     const DCArrayKokkos<material_t> material = simparam->material;
     CArray<double> current_element_nodal_densities = CArray<double>(num_nodes_in_elem);
 
-    std::vector<std::vector<int>> FEA_Module_My_TO_Modules = simparam->FEA_Module_My_TO_Modules;
     problem = Explicit_Solver_Pointer_->problem; // Pointer to ROL optimization problem object
     ROL::Ptr<ROL::Objective<real_t>> obj_pointer;
 
@@ -159,16 +158,6 @@ void FEA_Module_Dynamic_Elasticity::update_forward_solve(Teuchos::RCP<const MV>
     node_velocities_distributed->assign(*initial_node_velocities_distributed);
 
     // reset time accumulating objective and constraints
-    /*
-    for(int imodule = 0 ; imodule < FEA_Module_My_TO_Modules[my_fea_module_index_].size(); imodule++){
-      current_module_index = FEA_Module_My_TO_Modules[my_fea_module_index_][imodule];
-      //test if module needs reset
-      if(){
-
-      }
-    }
-    */
-    // simple setup to just request KE for now; above loop to be expanded and used later for scanning modules
     obj_pointer = problem->getObjective();
     KineticEnergyMinimize_TopOpt& kinetic_energy_minimize_function = dynamic_cast<KineticEnergyMinimize_TopOpt&>(*obj_pointer);
     kinetic_energy_minimize_function.objective_accumulation = 0;
@@ -1362,7 +1351,7 @@ void FEA_Module_Dynamic_Elasticity::init_assembly()
 
     // compute maximum stride
     size_t update = 0;
-    REDUCE_MAX_CLASS(inode, 0, nlocal_nodes, update, {
+    FOR_REDUCE_MAX_CLASS(inode, 0, nlocal_nodes, update, {
         if (update < Graph_Matrix_Strides_initial(inode)) {
             update = Graph_Matrix_Strides_initial(inode);
         }
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/mesh.h b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/mesh.h
index 0465f374e..a9ab1214e 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/mesh.h
+++ b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/mesh.h
@@ -1,5 +1,5 @@
 /**********************************************************************************************
- © 2020. Triad National Security, LLC. All rights reserved.
+ � 2020. Triad National Security, LLC. All rights reserved.
  This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
  National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
  Department of Energy/National Nuclear Security Administration. All rights in the program are
@@ -322,7 +322,7 @@ struct mesh_t
         // find the max number of elems around a node
         size_t max_num_elems_in_node;
         size_t max_num_lcl;
-        REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
+        FOR_REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
             // num_corners_in_node = num_elems_in_node
             size_t max_num = num_corners_in_node(node_gid);
 
@@ -681,7 +681,7 @@ struct mesh_t
         // find the max number of elems around a node
         size_t max_num_elems_in_node;
         size_t max_num_lcl;
-        REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
+        FOR_REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
             // num_corners_in_node = num_elems_in_node
             size_t max_num = num_corners_in_node(node_gid);
 
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/time_integration.cpp b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/time_integration.cpp
index 01f4ff8f4..2e9597ef5 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/time_integration.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/time_integration.cpp
@@ -1,5 +1,5 @@
 /**********************************************************************************************
- © 2020. Triad National Security, LLC. All rights reserved.
+ � 2020. Triad National Security, LLC. All rights reserved.
  This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
  National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
  Department of Energy/National Nuclear Security Administration. All rights in the program are
@@ -112,7 +112,7 @@ void FEA_Module_Dynamic_Elasticity::get_timestep(mesh_t& mesh,
     int    num_dims = num_dim;
     double dt_lcl;
     double min_dt_calc;
-    REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
+    FOR_REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
         double coords0[24];  // element coords
         ViewCArrayKokkos<double> coords(coords0, 8, 3);
 
@@ -218,7 +218,7 @@ void FEA_Module_Dynamic_Elasticity::get_timestep2D(mesh_t& mesh,
     double dt_lcl;
     double min_dt_calc;
 
-    REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
+    FOR_REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
         double coords0[8];  // element coords
         ViewCArrayKokkos<double> coords(coords0, 4, 2);
 
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/FEA_Module_Eulerian.cpp b/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/FEA_Module_Eulerian.cpp
index 164cd923f..832b3e9b5 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/FEA_Module_Eulerian.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/FEA_Module_Eulerian.cpp
@@ -833,7 +833,7 @@ void FEA_Module_Eulerian::euler_solve()
     int nlocal_elem_non_overlapping = Explicit_Solver_Pointer_->nlocal_elem_non_overlapping;
 
     // extensive IE
-    REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
         IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     }, IE_sum);
     IE_t0 = IE_sum;
@@ -841,7 +841,7 @@ void FEA_Module_Eulerian::euler_solve()
     MPI_Allreduce(&IE_t0, &global_IE_t0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // extensive KE
-    REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
         double ke = 0;
         for (size_t dim = 0; dim < num_dim; dim++) {
             ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
@@ -1248,7 +1248,7 @@ void FEA_Module_Eulerian::euler_solve()
                 KE_loc_sum = 0.0;
                 KE_sum     = 0.0;
                 // extensive KE
-                REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+                FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
                     double ke = 0;
                     for (size_t dim = 0; dim < num_dim; dim++) {
                         // midpoint integration approximation
@@ -1365,7 +1365,7 @@ void FEA_Module_Eulerian::euler_solve()
     KE_sum     = 0.0;
 
     // extensive IE
-    REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
         IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     }, IE_sum);
     IE_tend = IE_sum;
@@ -1374,7 +1374,7 @@ void FEA_Module_Eulerian::euler_solve()
     MPI_Allreduce(&IE_tend, &global_IE_tend, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // extensive KE
-    REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
         double ke = 0;
         for (size_t dim = 0; dim < num_dim; dim++) {
             ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/eulerian_time_integration.cpp b/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/eulerian_time_integration.cpp
index d6bdb6e4c..20b91cbcf 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/eulerian_time_integration.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/Eulerian_Solver/eulerian_time_integration.cpp
@@ -1,5 +1,5 @@
 /**********************************************************************************************
- © 2020. Triad National Security, LLC. All rights reserved.
+ � 2020. Triad National Security, LLC. All rights reserved.
  This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
  National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
  Department of Energy/National Nuclear Security Administration. All rights in the program are
@@ -113,7 +113,7 @@ void FEA_Module_Eulerian::get_timestep(mesh_t& mesh,
     int    num_dims = num_dim;
     double dt_lcl;
     double min_dt_calc;
-    REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
+    FOR_REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
         double coords0[24];  // element coords
         ViewCArrayKokkos<double> coords(coords0, 8, 3);
 
@@ -221,7 +221,7 @@ void FEA_Module_Eulerian::get_timestep2D(mesh_t& mesh,
     double dt_lcl;
     double min_dt_calc;
 
-    REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
+    FOR_REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
         double coords0[8];  // element coords
         ViewCArrayKokkos<double> coords(coords0, 4, 2);
 
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Explicit_Optimization_Function_Headers.h b/src/Parallel-Solvers/Parallel-Explicit/Explicit_Optimization_Function_Headers.h
new file mode 100644
index 000000000..1b8f342cc
--- /dev/null
+++ b/src/Parallel-Solvers/Parallel-Explicit/Explicit_Optimization_Function_Headers.h
@@ -0,0 +1,50 @@
+/**********************************************************************************************
+ © 2020. Triad National Security, LLC. All rights reserved.
+ This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
+ National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
+ Department of Energy/National Nuclear Security Administration. All rights in the program are
+ reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear
+ Security Administration. The Government is granted for itself and others acting on its behalf a
+ nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare
+ derivative works, distribute copies to the public, perform publicly and display publicly, and
+ to permit others to do so.
+ This program is open source under the BSD-3 License.
+ Redistribution and use in source and binary forms, with or without modification, are permitted
+ provided that the following conditions are met:
+ 
+ 1.  Redistributions of source code must retain the above copyright notice, this list of
+ conditions and the following disclaimer.
+ 
+ 2.  Redistributions in binary form must reproduce the above copyright notice, this list of
+ conditions and the following disclaimer in the documentation and/or other materials
+ provided with the distribution.
+ 
+ 3.  Neither the name of the copyright holder nor the names of its contributors may be used
+ to endorse or promote products derived from this software without specific prior
+ written permission.
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
+ CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+ OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+ OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ **********************************************************************************************/
+ 
+#ifndef EXPLICIT_OPTIMIZATION_FUNCTION_FILES_H
+#define EXPLICIT_OPTIMIZATION_FUNCTION_FILES_H 
+
+#include "Mass_Constraint_Top_Opt.h"
+#include "Mass_Constraint_Shape_Opt.h"
+#include "Moment_of_Inertia_Constraint_Top_Opt.h"
+#include "Kinetic_Energy_Minimize.h"
+#include "Internal_Energy_Minimize.h"
+#include "Internal_Energy_Minimize_Shape_Opt.h"
+#include "MMA_Objective.hpp"
+#include "Area_Normals.h"
+
+#endif // end HEADER_H
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.cpp b/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.cpp
index 771b98ab0..49f1d2b34 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.cpp
@@ -86,13 +86,7 @@
 #include <ROL_TpetraMultiVector.hpp>
 
 //Objective Functions and Constraint Functions
-//#include "Topology_Optimization_Function_Headers.h"
-#include "Mass_Constraint.h"
-#include "Moment_of_Inertia_Constraint.h"
-#include "Kinetic_Energy_Minimize.h"
-#include "Internal_Energy_Minimize.h"
-#include "MMA_Objective.hpp"
-#include "Area_Normals.h"
+#include "Explicit_Optimization_Function_Headers.h"
 
 #define BUFFER_LINES 20000
 #define MAX_WORD 30
@@ -323,10 +317,11 @@ void Explicit_Solver::run() {
   //set initial saved velocities
   initial_node_velocities_distributed->assign(*node_velocities_distributed);
     
-  if(simparam.topology_optimization_on || simparam.shape_optimization_on){
-      //design_node_densities_distributed->randomize(1,1);
-      setup_optimization_problem();
-      //problem = ROL::makePtr<ROL::Problem<real_t>>(obj,x);
+  if(simparam.topology_optimization_on){
+    setup_topology_optimization_problem();
+  }
+  else if(simparam.shape_optimization_on){
+    setup_shape_optimization_problem();
   }
   else{
     // ---------------------------------------------------------------------
@@ -1115,17 +1110,17 @@ void Explicit_Solver::FEA_module_setup(){
    Setup Optimization Problem Object, Relevant Objective, and Constraints
 ------------------------------------------------------------------------- */
 
-void Explicit_Solver::setup_optimization_problem(){
+void Explicit_Solver::setup_topology_optimization_problem(){
   int num_dim = simparam.num_dims;
   bool nodal_density_flag = simparam.nodal_density_flag;
-  int nTO_modules = simparam.TO_Module_List.size();
+  int nOpt_modules = simparam.Optimization_Module_List.size();
   //int nmulti_objective_modules = simparam->nmulti_objective_modules;
-  std::vector<TO_MODULE_TYPE> TO_Module_List = simparam.TO_Module_List;
-  std::vector<int> TO_Module_My_FEA_Module = simparam.TO_Module_My_FEA_Module;
+  std::vector<OPTIMIZATION_MODULE_TYPE> Optimization_Module_List = simparam.Optimization_Module_List;
+  std::vector<int> Optimization_Module_My_FEA_Module = simparam.Optimization_Module_My_FEA_Module;
   //std::vector<int> Multi_Objective_Modules = simparam->Multi_Objective_Modules;
   //std::vector<real_t> Multi_Objective_Weights = simparam->Multi_Objective_Weights;
   std::vector<std::vector<real_t>> Function_Arguments = simparam.Function_Arguments;
-  std::vector<FUNCTION_TYPE> TO_Function_Type = simparam.TO_Function_Type;
+  std::vector<FUNCTION_TYPE> Optimization_Function_Type = simparam.Optimization_Function_Type;
   std::vector<ROL::Ptr<ROL::Objective<real_t>>> Multi_Objective_Terms;
 
   std::string constraint_base, constraint_name;
@@ -1197,15 +1192,13 @@ void Explicit_Solver::setup_optimization_problem(){
     }
 
     if(simparam.optimization_options.method_of_moving_asymptotes){
-      Teuchos::RCP<MV> mma_upper_bound_node_densities_distributed = Teuchos::rcp(new MV(map, 1));
-      Teuchos::RCP<MV> mma_lower_bound_node_densities_distributed = Teuchos::rcp(new MV(map, 1));
-      //mma_lower_bound_node_densities_distributed->assign(*lower_bound_node_densities_distributed);
-      //mma_upper_bound_node_densities_distributed->assign(*upper_bound_node_densities_distributed);
+      Teuchos::RCP<MV> mma_upper_bound_distributed = Teuchos::rcp(new MV(map, 1));
+      Teuchos::RCP<MV> mma_lower_bound_distributed = Teuchos::rcp(new MV(map, 1));
       
-      mma_lower_bound_node_densities_distributed->putScalar(-0.1);
-      mma_upper_bound_node_densities_distributed->putScalar(0.1);
-      mma_lower_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(mma_lower_bound_node_densities_distributed);
-      mma_upper_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(mma_upper_bound_node_densities_distributed);
+      mma_lower_bound_distributed->putScalar(-0.1);
+      mma_upper_bound_distributed->putScalar(0.1);
+      mma_lower_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(mma_lower_bound_distributed);
+      mma_upper_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(mma_upper_bound_distributed);
       
       mma_bnd = ROL::makePtr<ROL::Bounds<real_t>>(mma_lower_bounds, mma_upper_bounds);
     }
@@ -1357,18 +1350,18 @@ void Explicit_Solver::setup_optimization_problem(){
   //Instantiate (the one) objective function for the problem
   ROL::Ptr<ROL::Objective<real_t>> obj, sub_obj;
   bool objective_declared = false;
-  for(int imodule = 0; imodule < nTO_modules; imodule++){
-    if(TO_Function_Type[imodule] == FUNCTION_TYPE::OBJECTIVE){
+  for(int imodule = 0; imodule < nOpt_modules; imodule++){
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::OBJECTIVE){
       //check if previous module already defined an objective, there must be one objective module
       if(objective_declared){
         // TODO: Put this validation earlier.
         *fos << "PROGRAM IS ENDING DUE TO ERROR; ANOTHER OBJECTIVE FUNCTION WITH NAME \"" 
-              << TO_Module_List[imodule] <<"\" ATTEMPTED TO REPLACE A PREVIOUS OBJECTIVE; THERE MUST BE ONE OBJECTIVE." << std::endl;
+              << Optimization_Module_List[imodule] <<"\" ATTEMPTED TO REPLACE A PREVIOUS OBJECTIVE; THERE MUST BE ONE OBJECTIVE." << std::endl;
           exit_solver(0);
       }
-      if(TO_Module_List[imodule] == TO_MODULE_TYPE::Kinetic_Energy_Minimize){
+      if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Kinetic_Energy_Minimize_TopOpt){
         //debug print
-        *fos << " KINETIC ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
+        *fos << " KINETIC ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
         if(simparam.optimization_options.method_of_moving_asymptotes){
           sub_obj = ROL::makePtr<KineticEnergyMinimize_TopOpt>(this, nodal_density_flag);
           obj = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);
@@ -1377,9 +1370,9 @@ void Explicit_Solver::setup_optimization_problem(){
           obj = ROL::makePtr<KineticEnergyMinimize_TopOpt>(this, nodal_density_flag);
         }
       }
-      else if(TO_Module_List[imodule] == TO_MODULE_TYPE::Internal_Energy_Minimize){
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_TopOpt){
         //debug print
-        *fos << " KINETIC ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
+        *fos << " KINETIC ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
         if(simparam.optimization_options.method_of_moving_asymptotes){
           sub_obj = ROL::makePtr<InternalEnergyMinimize_TopOpt>(this, nodal_density_flag);
           obj = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);
@@ -1391,7 +1384,7 @@ void Explicit_Solver::setup_optimization_problem(){
       else{
         // TODO: Put validation earlier
         *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED OBJECTIVE FUNCTION REQUESTED WITH NAME \"" 
-              << TO_Module_List[imodule] << "\"" << std::endl;
+              << Optimization_Module_List[imodule] << "\"" << std::endl;
         exit_solver(0);
       }
       objective_declared = true;
@@ -1422,35 +1415,35 @@ void Explicit_Solver::setup_optimization_problem(){
   //problem->addConstraint("equality Constraint 3",eq_constraint3,constraint_mul3);
   //problem->addLinearConstraint("Equality Constraint",eq_constraint,constraint_mul);
   
-  for(int imodule = 0; imodule < nTO_modules; imodule++){
+  for(int imodule = 0; imodule < nOpt_modules; imodule++){
     number_union.str(constraint_base);
     number_union << imodule + 1;
     constraint_name = number_union.str();
     ROL::Ptr<std::vector<real_t> > li_ptr = ROL::makePtr<std::vector<real_t>>(1,0.0);
     ROL::Ptr<ROL::Vector<real_t> > constraint_mul = ROL::makePtr<ROL::StdVector<real_t>>(li_ptr);
-    if(TO_Function_Type[imodule] == FUNCTION_TYPE::EQUALITY_CONSTRAINT){
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::EQUALITY_CONSTRAINT){
       //pointers are reference counting
       ROL::Ptr<ROL::Constraint<real_t>> eq_constraint;
-      if(TO_Module_List[imodule]==TO_MODULE_TYPE::Mass_Constraint){
+      if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Mass_Constraint_TopOpt){
         
-        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false, true);
+        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][0], false, true);
       }
-      else if(TO_Module_List[imodule]==TO_MODULE_TYPE::Moment_of_Inertia_Constraint){
-        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        eq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], false, true);
+      else if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt){
+        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], false, true);
       }
       else{
         // TODO: Put validation earlier
         *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED EQUALITY CONSTRAINT FUNCTION REQUESTED WITH NAME \"" 
-              << TO_Module_List[imodule] <<"\"" << std::endl;
+              << Optimization_Module_List[imodule] <<"\"" << std::endl;
         exit_solver(0);
       }
       *fos << " ADDING CONSTRAINT " << constraint_name << std::endl;
       problem->addConstraint(constraint_name, eq_constraint, constraint_mul);
     }
 
-    if(TO_Function_Type[imodule] == FUNCTION_TYPE::INEQUALITY_CONSTRAINT){
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::INEQUALITY_CONSTRAINT){
       //pointers are reference counting
       ROL::Ptr<ROL::Constraint<real_t>> ineq_constraint;
       ROL::Ptr<std::vector<real_t> > ll_ptr = ROL::makePtr<std::vector<real_t>>(1,Function_Arguments[imodule][0]);
@@ -1458,18 +1451,18 @@ void Explicit_Solver::setup_optimization_problem(){
       ROL::Ptr<ROL::Vector<real_t> > ll = ROL::makePtr<ROL::StdVector<real_t>>(ll_ptr);
       ROL::Ptr<ROL::Vector<real_t> > lu = ROL::makePtr<ROL::StdVector<real_t>>(lu_ptr);
       ROL::Ptr<ROL::BoundConstraint<real_t>> constraint_bnd = ROL::makePtr<ROL::Bounds<real_t>>(ll,lu);
-      if(TO_Module_List[imodule]==TO_MODULE_TYPE::Mass_Constraint){
-        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, true, true);
+      if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Mass_Constraint_TopOpt){
+        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<MassConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, 0, true, true);
       }
-      else if(TO_Module_List[imodule]==TO_MODULE_TYPE::Moment_of_Inertia_Constraint){
-        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<TO_Module_My_FEA_Module[imodule] << std::endl;
-        ineq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[TO_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], true, true);
+      else if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt){
+        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], true, true);
       }
       else{
         // TODO: Put this validation earlier
         *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED INEQUALITY CONSTRAINT FUNCTION REQUESTED WITH NAME \"" 
-              << TO_Module_List[imodule] << "\"" << std::endl;
+              << Optimization_Module_List[imodule] << "\"" << std::endl;
         exit_solver(0);
       }
       *fos << " ADDING CONSTRAINT " << constraint_name << std::endl;
@@ -1567,6 +1560,255 @@ void Explicit_Solver::setup_optimization_problem(){
     //std::cout << "Final Mass Constraint is " << final_mass/initial_mass << std::endl;
 }
 
+/* ----------------------------------------------------------------------
+   Setup Optimization Problem Object, Relevant Objective, and Constraints
+------------------------------------------------------------------------- */
+
+void Explicit_Solver::setup_shape_optimization_problem(){
+  int num_dim = simparam.num_dims;
+  int nOpt_modules = simparam.Optimization_Module_List.size();
+  //int nmulti_objective_modules = simparam->nmulti_objective_modules;
+  std::vector<OPTIMIZATION_MODULE_TYPE> Optimization_Module_List = simparam.Optimization_Module_List;
+  std::vector<int> Optimization_Module_My_FEA_Module = simparam.Optimization_Module_My_FEA_Module;
+  //std::vector<int> Multi_Objective_Modules = simparam->Multi_Objective_Modules;
+  //std::vector<real_t> Multi_Objective_Weights = simparam->Multi_Objective_Weights;
+  std::vector<std::vector<real_t>> Function_Arguments = simparam.Function_Arguments;
+  std::vector<FUNCTION_TYPE> Optimization_Function_Type = simparam.Optimization_Function_Type;
+  std::vector<ROL::Ptr<ROL::Objective<real_t>>> Multi_Objective_Terms;
+
+  std::string constraint_base, constraint_name;
+  std::stringstream number_union;
+  CArray<GO> Surface_Nodes;
+  GO current_node_index, current_element_index;
+  LO local_node_index, local_element_id;
+  int num_bdy_patches_in_set;
+  size_t node_id, patch_id, module_id;
+  int num_boundary_sets;
+  int local_surface_id;
+  const_host_vec_array design_coordinates;
+  typedef ROL::TpetraMultiVector<real_t,LO,GO,node_type> ROL_MV;
+  const_host_elem_conn_array nodes_in_elem = global_nodes_in_elem_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadOnly);
+  bool nodal_density_flag = true;
+  
+  // fill parameter list with desired algorithmic options or leave as default
+  // Read optimization input parameter list.
+
+  Teuchos::RCP<Teuchos::ParameterList> parlist;
+  if(simparam.optimization_options.optimization_parameters_xml_file){
+    std::string xmlFileName = simparam.optimization_options.xml_parameters_file_name;
+    
+    //check if parameter file exists
+    std::ifstream param_file_check(xmlFileName);
+    if (!param_file_check.is_open()) {
+        *fos << "Unable to find xml parameter file required for optimization with the ROL library"  << std::endl;
+        exit_solver(0);
+    }
+    param_file_check.close();
+
+    parlist = ROL::getParametersFromXmlFile(xmlFileName);
+  }
+  else{
+    parlist = Teuchos::rcp(new Teuchos::ParameterList("Inputs"));
+    set_rol_params(parlist);
+  }
+  //ROL::ParameterList parlist;
+
+  ROL::Ptr<ROL::BoundConstraint<real_t> > bnd, mma_bnd;
+  // Bound constraint defining the possible range of design coordinates variables
+  ROL::Ptr<ROL::Vector<real_t> > lower_bounds, mma_lower_bounds;
+  ROL::Ptr<ROL::Vector<real_t> > upper_bounds, mma_upper_bounds;
+
+  //set bounds on design variables
+  //allocate global vector information
+  upper_bound_node_coordinates_distributed = Teuchos::rcp(new MV(map, num_dim));
+  //all_lower_bound_node_coordinates_distributed = Teuchos::rcp(new MV(all_node_map, 1));
+  //lower_bound_node_coordinates_distributed = Teuchos::rcp(new MV(*all_lower_bound_node_coordinates_distributed, map));
+  lower_bound_node_coordinates_distributed = Teuchos::rcp(new MV(map, num_dim));
+  host_vec_array node_coordinates_upper_bound = upper_bound_node_coordinates_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+  host_vec_array node_coordinates_lower_bound = lower_bound_node_coordinates_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+
+  //begin by assigning values of coordinate vectors to bound vectors
+  lower_bound_node_coordinates_distributed->assign(*all_node_coords_distributed);
+  upper_bound_node_coordinates_distributed->assign(*all_node_coords_distributed);
+
+  //initialize densities to 1 for now; in the future there might be an option to read in an initial condition for each node
+  for(int inode = 0; inode < nlocal_nodes; inode++){
+    for(int idim = 0; idim < num_dim; idim++){
+      node_coordinates_upper_bound(inode,idim) += simparam.optimization_options.max_coord_move_length;
+      node_coordinates_lower_bound(inode,idim) -= simparam.optimization_options.max_coord_move_length;
+    }
+  }
+
+  if(simparam.optimization_options.method_of_moving_asymptotes){
+    Teuchos::RCP<MV> mma_upper_bound_distributed = Teuchos::rcp(new MV(map, num_dim));
+    Teuchos::RCP<MV> mma_lower_bound_distributed = Teuchos::rcp(new MV(map, num_dim));
+    
+    mma_lower_bound_distributed->putScalar(-0.1);
+    mma_upper_bound_distributed->putScalar(0.1);
+    mma_lower_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(mma_lower_bound_distributed);
+    mma_upper_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(mma_upper_bound_distributed);
+    
+    mma_bnd = ROL::makePtr<ROL::Bounds<real_t>>(mma_lower_bounds, mma_upper_bounds);
+  }
+  
+
+  //Design variables to optimize
+  ROL::Ptr<ROL::Vector<real_t>> x;
+  x = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(design_node_coords_distributed);
+  
+  //Instantiate (the one) objective function for the problem
+  ROL::Ptr<ROL::Objective<real_t>> obj, sub_obj;
+  bool objective_declared = false;
+  for(int imodule = 0; imodule < nOpt_modules; imodule++){
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::OBJECTIVE){
+      //check if previous module already defined an objective, there must be one objective module
+      if(objective_declared){
+        // TODO: Put this validation earlier.
+        *fos << "PROGRAM IS ENDING DUE TO ERROR; ANOTHER OBJECTIVE FUNCTION WITH NAME \"" 
+              << Optimization_Module_List[imodule] <<"\" ATTEMPTED TO REPLACE A PREVIOUS OBJECTIVE; THERE MUST BE ONE OBJECTIVE." << std::endl;
+          exit_solver(0);
+      }
+      if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Kinetic_Energy_Minimize_TopOpt){
+        //debug print
+        *fos << " KINETIC ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        if(simparam.optimization_options.method_of_moving_asymptotes){
+          sub_obj = ROL::makePtr<KineticEnergyMinimize_TopOpt>(this, nodal_density_flag);
+          obj = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);
+        }
+        else{
+          obj = ROL::makePtr<KineticEnergyMinimize_TopOpt>(this, nodal_density_flag);
+        }
+      }
+      else if(Optimization_Module_List[imodule] == OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_ShapeOpt){
+        //debug print
+        *fos << " KINETIC ENERGY OBJECTIVE EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        if(simparam.optimization_options.method_of_moving_asymptotes){
+          sub_obj = ROL::makePtr<InternalEnergyMinimize_ShapeOpt>(this);
+          obj = ROL::makePtr<ObjectiveMMA>(sub_obj, mma_bnd, x);
+        }
+        else{
+          obj = ROL::makePtr<InternalEnergyMinimize_ShapeOpt>(this);
+        }
+      }
+      else{
+        // TODO: Put validation earlier
+        *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED OBJECTIVE FUNCTION REQUESTED WITH NAME \"" 
+              << Optimization_Module_List[imodule] << "\"" << std::endl;
+        exit_solver(0);
+      }
+      objective_declared = true;
+    }
+  }
+  
+  //optimization problem interface that can have constraints added to it before passing to solver object
+  problem = ROL::makePtr<ROL::Problem<real_t>>(obj,x);
+  
+  for(int imodule = 0; imodule < nOpt_modules; imodule++){
+    number_union.str(constraint_base);
+    number_union << imodule + 1;
+    constraint_name = number_union.str();
+    ROL::Ptr<std::vector<real_t> > li_ptr = ROL::makePtr<std::vector<real_t>>(1,0.0);
+    ROL::Ptr<ROL::Vector<real_t> > constraint_mul = ROL::makePtr<ROL::StdVector<real_t>>(li_ptr);
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::EQUALITY_CONSTRAINT){
+      //pointers are reference counting
+      ROL::Ptr<ROL::Constraint<real_t>> eq_constraint;
+      if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Mass_Constraint_ShapeOpt){
+        
+        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<MassConstraint_ShapeOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], Function_Arguments[imodule][0], false, true);
+      }
+      else if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt){
+        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        eq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], nodal_density_flag, Function_Arguments[imodule][1], Function_Arguments[imodule][0], false, true);
+      }
+      else{
+        // TODO: Put validation earlier
+        *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED EQUALITY CONSTRAINT FUNCTION REQUESTED WITH NAME \"" 
+              << Optimization_Module_List[imodule] <<"\"" << std::endl;
+        exit_solver(0);
+      }
+      *fos << " ADDING CONSTRAINT " << constraint_name << std::endl;
+      problem->addConstraint(constraint_name, eq_constraint, constraint_mul);
+    }
+
+    if(Optimization_Function_Type[imodule] == FUNCTION_TYPE::INEQUALITY_CONSTRAINT){
+      //pointers are reference counting
+      ROL::Ptr<ROL::Constraint<real_t>> ineq_constraint;
+      ROL::Ptr<std::vector<real_t> > ll_ptr = ROL::makePtr<std::vector<real_t>>(1,Function_Arguments[imodule][0]);
+      ROL::Ptr<std::vector<real_t> > lu_ptr = ROL::makePtr<std::vector<real_t>>(1,Function_Arguments[imodule][1]);    
+      ROL::Ptr<ROL::Vector<real_t> > ll = ROL::makePtr<ROL::StdVector<real_t>>(ll_ptr);
+      ROL::Ptr<ROL::Vector<real_t> > lu = ROL::makePtr<ROL::StdVector<real_t>>(lu_ptr);
+      ROL::Ptr<ROL::BoundConstraint<real_t>> constraint_bnd = ROL::makePtr<ROL::Bounds<real_t>>(ll,lu);
+      if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Mass_Constraint_ShapeOpt){
+        *fos << " MASS CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<MassConstraint_ShapeOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], 0, true, true);
+      }
+      else if(Optimization_Module_List[imodule]==OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt){
+        *fos << " MOMENT OF INERTIA CONSTRAINT EXPECTS FEA MODULE INDEX " <<Optimization_Module_My_FEA_Module[imodule] << std::endl;
+        ineq_constraint = ROL::makePtr<MomentOfInertiaConstraint_TopOpt>(fea_modules[Optimization_Module_My_FEA_Module[imodule]], Function_Arguments[imodule][1], Function_Arguments[imodule][0], true, true);
+      }
+      else{
+        // TODO: Put this validation earlier
+        *fos << "PROGRAM IS ENDING DUE TO ERROR; UNDEFINED INEQUALITY CONSTRAINT FUNCTION REQUESTED WITH NAME \"" 
+              << Optimization_Module_List[imodule] << "\"" << std::endl;
+        exit_solver(0);
+      }
+      *fos << " ADDING CONSTRAINT " << constraint_name << std::endl;
+      problem->addConstraint(constraint_name, ineq_constraint, constraint_mul, constraint_bnd);
+    }
+    number_union.clear();
+  }
+
+  lower_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(lower_bound_node_coordinates_distributed);
+  upper_bounds = ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO>>(upper_bound_node_coordinates_distributed);
+  
+  bnd = ROL::makePtr<ROL::Bounds<real_t>>(lower_bounds, upper_bounds);
+  problem->addBoundConstraint(bnd);
+  
+  //real_t initial_mass = ROL_Element_Masses->reduce(sumreduc);
+
+  problem->setProjectionAlgorithm(*parlist);
+  //finalize problem
+  problem->finalize(false,true,*fos);
+  //problem->check(true,std::cout);
+
+  //debug checks
+  ROL::Ptr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>> rol_x =
+   ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>>(design_node_coords_distributed);
+  //construct direction vector for check
+  Teuchos::RCP<MV> directions_distributed = Teuchos::rcp(new MV(map, num_dim));
+  directions_distributed->putScalar(-0.1);
+  directions_distributed->randomize(-0.8,1);
+  Kokkos::View <real_t*, array_layout, HostSpace, memory_traits> direction_norm("gradient norm",1);
+  directions_distributed->norm2(direction_norm);
+  directions_distributed->scale(1/direction_norm(0));
+  //set all but first component to 0 for debug
+  host_vec_array directions = directions_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+  //for(int init = 1; init < nlocal_nodes; init++)
+  //directions(4,0) = -0.3;
+  ROL::Ptr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>> rol_d =
+  ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>>(directions_distributed);
+  //obj->checkGradient(*rol_x, *rol_d);
+  //obj->checkHessVec(*rol_x, *rol_d);
+  //directions_distributed->putScalar(-0.000001);
+  //obj->checkGradient(*rol_x, *rol_d);
+  //directions_distributed->putScalar(-0.0000001);
+  //obj->checkGradient(*rol_x, *rol_d);
+  
+
+  // Instantiate Solver.
+  ROL::Solver<real_t> solver(problem,*parlist);
+    
+  // Solve optimization problem.
+  //std::ostream outStream;
+  solver.solve(*fos);
+
+  //print final constraint satisfaction
+  //fea_elasticity->compute_element_masses(design_densities,false);
+  //real_t final_mass = ROL_Element_Masses->reduce(sumreduc);
+  //if(myrank==0)
+    //std::cout << "Final Mass Constraint is " << final_mass/initial_mass << std::endl;
+}
 
 /* ----------------------------------------------------------------------------
    Initialize sets of element boundary surfaces and arrays for input conditions
@@ -3363,15 +3605,27 @@ void Explicit_Solver::ensight_writer(){
 
 void Explicit_Solver::init_design(){
   int num_dim = simparam.num_dims;
-  bool nodal_density_flag = simparam.nodal_density_flag;
+  bool topology_optimization_on = simparam.topology_optimization_on;
+  bool shape_optimization_on = simparam.shape_optimization_on;
   Input_Options input_options;
   if(simparam.input_options.has_value()){
     input_options = simparam.input_options.value();
   }
-  //set densities
-  if(nodal_density_flag){
-    if(simparam.input_options.has_value()){
-      if (!input_options.topology_optimization_restart) {
+  if(topology_optimization_on){
+    bool nodal_density_flag = simparam.nodal_density_flag;
+    //set densities
+    if(nodal_density_flag){
+      if(simparam.input_options.has_value()){
+        if (!input_options.topology_optimization_restart) {
+          design_node_densities_distributed = Teuchos::rcp(new MV(map, 1));
+          host_vec_array node_densities = design_node_densities_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+        
+          for(int inode = 0; inode < nlocal_nodes; inode++){
+            node_densities(inode,0) = 1;
+          }
+        }
+      }
+      else{
         design_node_densities_distributed = Teuchos::rcp(new MV(map, 1));
         host_vec_array node_densities = design_node_densities_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
       
@@ -3379,40 +3633,56 @@ void Explicit_Solver::init_design(){
           node_densities(inode,0) = 1;
         }
       }
+        //allocate global vector information
+      all_node_densities_distributed = Teuchos::rcp(new MV(all_node_map, 1));
+
+      //communicate ghost information to the all vector
+      //create import object using local node indices map and all indices map
+      //Tpetra::Import<LO, GO> importer(map, all_node_map);
+      
+      //design_node_densities_distributed->randomize(0.1,1);
+      //comms to get ghosts
+      all_node_densities_distributed->doImport(*design_node_densities_distributed, *importer, Tpetra::INSERT);
+      
     }
     else{
-      design_node_densities_distributed = Teuchos::rcp(new MV(map, 1));
-      host_vec_array node_densities = design_node_densities_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
-    
-      for(int inode = 0; inode < nlocal_nodes; inode++){
-        node_densities(inode,0) = 1;
-      }
-    }
-      //allocate global vector information
-    all_node_densities_distributed = Teuchos::rcp(new MV(all_node_map, 1));
+      //initialize memory for volume storage
+      vec_array Element_Densities("Element Densities", rnum_elem, 1);
+      for(int ielem = 0; ielem < rnum_elem; ielem++)
+        Element_Densities(ielem,0) = 1;
 
-    //communicate ghost information to the all vector
-    //create import object using local node indices map and all indices map
-    Tpetra::Import<LO, GO> importer(map, all_node_map);
-    
-    //design_node_densities_distributed->randomize(0.1,1);
-    //comms to get ghosts
-    all_node_densities_distributed->doImport(*design_node_densities_distributed, importer, Tpetra::INSERT);
-    
-  }
-  else{
-    //initialize memory for volume storage
-    vec_array Element_Densities("Element Densities", rnum_elem, 1);
-    for(int ielem = 0; ielem < rnum_elem; ielem++)
-      Element_Densities(ielem,0) = 1;
+      //create global vector
+      Global_Element_Densities = Teuchos::rcp(new MV(all_element_map, Element_Densities));
 
-    //create global vector
-    Global_Element_Densities = Teuchos::rcp(new MV(all_element_map, Element_Densities));
+      //if(myrank==0)
+      //*fos << "Global Element Densities:" << std::endl;
+      //Global_Element_Densities->describe(*fos,Teuchos::VERB_EXTREME);
+      //*fos << std::endl;
+    }
+  }
+  else if(shape_optimization_on){
+    //set initial design coordinates (not equal to initial coordinates before optimization when restarting)
+      if(simparam.input_options.has_value()){
+        if (!input_options.shape_optimization_restart) {
+          design_node_coords_distributed = Teuchos::rcp(new MV(map, num_dim));
+          design_node_coords_distributed->assign(*node_coords_distributed);
+        }
+      }
+      else{
+        design_node_coords_distributed = Teuchos::rcp(new MV(map, num_dim));
+        design_node_coords_distributed->assign(*node_coords_distributed);
+      }
+        //allocate global vector information
+      all_design_node_coords_distributed = Teuchos::rcp(new MV(all_node_map, num_dim));
 
-    //if(myrank==0)
-    //*fos << "Global Element Densities:" << std::endl;
-    //Global_Element_Densities->describe(*fos,Teuchos::VERB_EXTREME);
-    //*fos << std::endl;
+      //communicate ghost information to the all vector
+      //create import object using local node indices map and all indices map
+      //Tpetra::Import<LO, GO> importer(map, all_node_map);
+      
+      //design_node_densities_distributed->randomize(0.1,1);
+      //comms to get ghosts
+      all_design_node_coords_distributed->doImport(*design_node_coords_distributed, *importer, Tpetra::INSERT);   
+    
   }
 
 }
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.h b/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.h
index cf432e0df..db2c26a4f 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.h
+++ b/src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.h
@@ -88,7 +88,9 @@ class Explicit_Solver : public Solver
     // process input to decide TO problem and FEA modules
     void FEA_module_setup();
 
-    void setup_optimization_problem();
+    void setup_topology_optimization_problem();
+
+    void setup_shape_optimization_problem();
 
     // initialize data for boundaries of the model and storage for boundary conditions and applied loads
     void init_boundaries();
diff --git a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/FEA_Module_SGH.h b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/FEA_Module_SGH.h
index ba7ccac63..2734eda9a 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/FEA_Module_SGH.h
+++ b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/FEA_Module_SGH.h
@@ -508,6 +508,16 @@ class FEA_Module_SGH : public FEA_Module
     void compute_topology_optimization_gradient_IVP(Teuchos::RCP<const MV> design_densities_distributed, Teuchos::RCP<MV> design_gradients_distributed,
                                                       unsigned long cycle, real_t global_dt);
 
+    void compute_shape_optimization_adjoint_full(Teuchos::RCP<const MV> design_coordinates_distributed){} // Force depends on node coords, velocity, and sie
+
+    void compute_shape_optimization_gradient_full(Teuchos::RCP<const MV> design_coordinates_distributed, Teuchos::RCP<MV> design_gradients_distributed){}
+
+    void compute_shape_optimization_gradient_tally(Teuchos::RCP<const MV> design_coordinates_distributed, Teuchos::RCP<MV> design_gradients_distributed,
+                                                      unsigned long cycle, real_t global_dt){}
+
+    void compute_shape_optimization_gradient_IVP(Teuchos::RCP<const MV> design_coordinates_distributed, Teuchos::RCP<MV> design_gradients_distributed,
+                                                      unsigned long cycle, real_t global_dt){}
+
     void boundary_adjoint(const mesh_t& mesh,
                           const DCArrayKokkos<boundary_t>& boundary,
                           vec_array& node_adjoint,
diff --git a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/mesh.h b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/mesh.h
index 883f4f156..667abf381 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/mesh.h
+++ b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/include/mesh.h
@@ -1,5 +1,5 @@
 /**********************************************************************************************
- © 2020. Triad National Security, LLC. All rights reserved.
+ � 2020. Triad National Security, LLC. All rights reserved.
  This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
  National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
  Department of Energy/National Nuclear Security Administration. All rights in the program are
@@ -318,7 +318,7 @@ struct mesh_t
         // find the max number of elems around a node
         size_t max_num_elems_in_node;
         size_t max_num_lcl;
-        REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
+        FOR_REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
             // num_corners_in_node = num_elems_in_node
             size_t max_num = num_corners_in_node(node_gid);
 
@@ -677,7 +677,7 @@ struct mesh_t
         // find the max number of elems around a node
         size_t max_num_elems_in_node;
         size_t max_num_lcl;
-        REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
+        FOR_REDUCE_MAX_CLASS(node_gid, 0, num_nodes, max_num_lcl, {
             // num_corners_in_node = num_elems_in_node
             size_t max_num = num_corners_in_node(node_gid);
 
diff --git a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/FEA_Module_SGH.cpp b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/FEA_Module_SGH.cpp
index 76694700b..974da50d9 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/FEA_Module_SGH.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/FEA_Module_SGH.cpp
@@ -945,11 +945,8 @@ void FEA_Module_SGH::sgh_solve()
     size_t num_bdy_nodes = mesh->num_bdy_nodes;
     size_t cycle;
     real_t objective_accumulation, global_objective_accumulation;
-
-    int nTO_modules;
     int old_max_forward_buffer;
 
-    std::vector<std::vector<int>> FEA_Module_My_TO_Modules = simparam->FEA_Module_My_TO_Modules;
     problem = Explicit_Solver_Pointer_->problem; // Pointer to ROL optimization problem object
     ROL::Ptr<ROL::Objective<real_t>> obj_pointer;
     
@@ -964,6 +961,8 @@ void FEA_Module_SGH::sgh_solve()
     num_active_checkpoints = 0;
     bool time_accumulation;
 
+    TpetraDFArray<double> test_nodes(num_nodes, num_dim);
+
     
     if(simparam->optimization_options.disable_forward_solve_output){
         //sets the value large enough to not write during the sgh loop
@@ -971,16 +970,6 @@ void FEA_Module_SGH::sgh_solve()
     }
 
     // reset time accumulating objective and constraints
-    /*
-    for(int imodule = 0 ; imodule < FEA_Module_My_TO_Modules[my_fea_module_index_].size(); imodule++){
-    current_module_index = FEA_Module_My_TO_Modules[my_fea_module_index_][imodule];
-    //test if module needs reset
-    if(){
-
-    }
-    }
-    */
-    // simple setup to just request KE for now; above loop to be expanded and used later for scanning modules
     if (topology_optimization_on||shape_optimization_on) {
         obj_pointer = problem->getObjective();
         objective_function = dynamic_cast<FierroOptimizationObjective*>(obj_pointer.getRawPtr());
@@ -1008,10 +997,6 @@ void FEA_Module_SGH::sgh_solve()
         }
     }
 
-    if (topology_optimization_on) {
-        nTO_modules = simparam->TO_Module_List.size();
-    }
-
     int myrank = Explicit_Solver_Pointer_->myrank;
     if (simparam->output_options.write_initial&&!simparam->optimization_options.disable_forward_solve_output) {
         if (myrank == 0) {
@@ -1054,7 +1039,7 @@ void FEA_Module_SGH::sgh_solve()
     int nlocal_elem_non_overlapping = Explicit_Solver_Pointer_->nlocal_elem_non_overlapping;
 
     // extensive IE
-    REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
         IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     }, IE_sum);
     IE_t0 = IE_sum;
@@ -1062,7 +1047,7 @@ void FEA_Module_SGH::sgh_solve()
     MPI_Allreduce(&IE_t0, &global_IE_t0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // extensive KE
-    REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
         double ke = 0;
         for (size_t dim = 0; dim < num_dim; dim++) {
             ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
@@ -1855,7 +1840,7 @@ void FEA_Module_SGH::sgh_solve()
     KE_sum     = 0.0;
 
     // extensive IE
-    REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
         IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     }, IE_sum);
     IE_tend = IE_sum;
@@ -1864,7 +1849,7 @@ void FEA_Module_SGH::sgh_solve()
     MPI_Allreduce(&IE_tend, &global_IE_tend, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // extensive KE
-    REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
         double ke = 0;
         for (size_t dim = 0; dim < num_dim; dim++) {
             ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
diff --git a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/sgh_optimization.cpp b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/sgh_optimization.cpp
index 3aba79eb2..0c2429852 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/sgh_optimization.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/sgh_optimization.cpp
@@ -105,8 +105,6 @@ void FEA_Module_SGH::update_forward_solve(Teuchos::RCP<const MV> zp, bool print_
     const DCArrayKokkos<boundary_t> boundary = module_params->boundary;
     const DCArrayKokkos<material_t> material = simparam->material;
     CArray<double> current_element_nodal_densities = CArray<double>(num_nodes_in_elem);
-
-    std::vector<std::vector<int>> FEA_Module_My_TO_Modules = simparam->FEA_Module_My_TO_Modules;
     problem = Explicit_Solver_Pointer_->problem; // Pointer to ROL optimization problem object
     ROL::Ptr<ROL::Objective<real_t>> obj_pointer;
     
@@ -162,16 +160,6 @@ void FEA_Module_SGH::update_forward_solve(Teuchos::RCP<const MV> zp, bool print_
     node_velocities_distributed->assign(*initial_node_velocities_distributed);
 
     // reset time accumulating objective and constraints
-    /*
-    for(int imodule = 0 ; imodule < FEA_Module_My_TO_Modules[my_fea_module_index_].size(); imodule++){
-      current_module_index = FEA_Module_My_TO_Modules[my_fea_module_index_][imodule];
-      //test if module needs reset
-      if(){
-
-      }
-    }
-    */
-    // simple setup to just request KE for now; above loop to be expanded and used later for scanning modules
     obj_pointer = problem->getObjective();
     objective_function = dynamic_cast<FierroOptimizationObjective*>(obj_pointer.getRawPtr());
     objective_function->objective_accumulation = 0;
@@ -2438,7 +2426,7 @@ void FEA_Module_SGH::init_assembly()
 
     // compute maximum stride
     size_t update = 0;
-    REDUCE_MAX_CLASS(inode, 0, nlocal_nodes, update, {
+    FOR_REDUCE_MAX_CLASS(inode, 0, nlocal_nodes, update, {
         if (update < Graph_Matrix_Strides_initial(inode)) {
             update = Graph_Matrix_Strides_initial(inode);
         }
@@ -2968,7 +2956,7 @@ void FEA_Module_SGH::checkpoint_solve(std::set<Dynamic_Checkpoint>::iterator sta
     // int nlocal_elem_non_overlapping = Explicit_Solver_Pointer_->nlocal_elem_non_overlapping;
 
     // // extensive IE
-    // REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    // FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
     //     IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     // }, IE_sum);
     // IE_t0 = IE_sum;
@@ -2976,7 +2964,7 @@ void FEA_Module_SGH::checkpoint_solve(std::set<Dynamic_Checkpoint>::iterator sta
     // MPI_Allreduce(&IE_t0, &global_IE_t0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // // extensive KE
-    // REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    // FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
     //     double ke = 0;
     //     for (size_t dim = 0; dim < num_dim; dim++) {
     //         ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
@@ -3556,7 +3544,7 @@ void FEA_Module_SGH::checkpoint_solve(std::set<Dynamic_Checkpoint>::iterator sta
     // KE_sum     = 0.0;
 
     // // extensive IE
-    // REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+    // FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
     //     IE_loc_sum += elem_mass(elem_gid) * elem_sie(rk_level, elem_gid);
     // }, IE_sum);
     // IE_tend = IE_sum;
@@ -3565,7 +3553,7 @@ void FEA_Module_SGH::checkpoint_solve(std::set<Dynamic_Checkpoint>::iterator sta
     // MPI_Allreduce(&IE_tend, &global_IE_tend, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
     // // extensive KE
-    // REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+    // FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
     //     double ke = 0;
     //     for (size_t dim = 0; dim < num_dim; dim++) {
     //         ke += node_vel(rk_level, node_gid, dim) * node_vel(rk_level, node_gid, dim); // 1/2 at end
diff --git a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/time_integration.cpp b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/time_integration.cpp
index c07a05eb1..f7906babd 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/time_integration.cpp
+++ b/src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/src/time_integration.cpp
@@ -1,5 +1,5 @@
 /**********************************************************************************************
- © 2020. Triad National Security, LLC. All rights reserved.
+ � 2020. Triad National Security, LLC. All rights reserved.
  This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
  National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
  Department of Energy/National Nuclear Security Administration. All rights in the program are
@@ -111,7 +111,7 @@ void FEA_Module_SGH::get_timestep(mesh_t& mesh,
     int    num_dims = num_dim;
     double dt_lcl;
     double min_dt_calc;
-    REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
+    FOR_REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
         double coords0[24];  // element coords
         ViewCArrayKokkos<double> coords(coords0, 8, 3);
 
@@ -217,7 +217,7 @@ void FEA_Module_SGH::get_timestep2D(mesh_t& mesh,
     double dt_lcl;
     double min_dt_calc;
 
-    REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
+    FOR_REDUCE_MIN_CLASS(elem_gid, 0, rnum_elem, dt_lcl, {
         double coords0[8];  // element coords
         ViewCArrayKokkos<double> coords(coords0, 4, 2);
 
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Internal_Energy_Minimize.h b/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Internal_Energy_Minimize.h
index b41516140..885468e47 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Internal_Energy_Minimize.h
+++ b/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Internal_Energy_Minimize.h
@@ -315,7 +315,7 @@ typedef MV::dual_view_type dual_vec_array;
             int nobj_volumes = FEM_SGH_->simparam->optimization_options.optimization_objective_regions.size();
             auto optimization_objective_regions = FEM_SGH_->simparam->optimization_options.optimization_objective_regions;
             const_vec_array all_initial_node_coords = FEM_SGH_->all_initial_node_coords_distributed->getLocalView<device_type>(Tpetra::Access::ReadOnly);
-            REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+            FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
                 int node_id;
                 double current_elem_coords[3];
                 bool contained = false;
@@ -341,7 +341,7 @@ typedef MV::dual_view_type dual_vec_array;
             }, IE_sum);
         }
         else{
-            REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+            FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
                 IE_loc_sum += elem_mass(elem_gid) * 0.5 * (current_elem_sie(elem_gid,0)+previous_elem_sie(elem_gid,0));
             }, IE_sum);
         }
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Internal_Energy_Minimize_Shape_Opt.h b/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Internal_Energy_Minimize_Shape_Opt.h
new file mode 100644
index 000000000..732faa87d
--- /dev/null
+++ b/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Internal_Energy_Minimize_Shape_Opt.h
@@ -0,0 +1,471 @@
+/**********************************************************************************************
+ © 2020. Triad National Security, LLC. All rights reserved.
+ This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos
+ National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S.
+ Department of Energy/National Nuclear Security Administration. All rights in the program are
+ reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear
+ Security Administration. The Government is granted for itself and others acting on its behalf a
+ nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare
+ derivative works, distribute copies to the public, perform publicly and display publicly, and
+ to permit others to do so.
+ This program is open source under the BSD-3 License.
+ Redistribution and use in source and binary forms, with or without modification, are permitted
+ provided that the following conditions are met:
+ 1.  Redistributions of source code must retain the above copyright notice, this list of
+ conditions and the following disclaimer.
+ 2.  Redistributions in binary form must reproduce the above copyright notice, this list of
+ conditions and the following disclaimer in the documentation and/or other materials
+ provided with the distribution.
+ 3.  Neither the name of the copyright holder nor the names of its contributors may be used
+ to endorse or promote products derived from this software without specific prior
+ written permission.
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
+ CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+ OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+ OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ **********************************************************************************************/
+
+#ifndef INTERNAL_ENERGY_MINIMIZE_SHAPE_OPT_H
+#define INTERNAL_ENERGY_MINIMIZE_SHAPE_OPT_H
+
+#include "matar.h"
+#include "elements.h"
+#include <string>
+#include <Teuchos_ScalarTraits.hpp>
+#include <Teuchos_RCP.hpp>
+#include <Teuchos_oblackholestream.hpp>
+#include <Teuchos_Tuple.hpp>
+#include <Teuchos_VerboseObject.hpp>
+
+#include <Tpetra_Core.hpp>
+#include <Tpetra_Map.hpp>
+#include <Tpetra_MultiVector.hpp>
+#include "Tpetra_Details_makeColMap.hpp"
+#include "Tpetra_Details_DefaultTypes.hpp"
+
+#include "ROL_Types.hpp"
+#include <ROL_TpetraMultiVector.hpp>
+#include "ROL_Objective.hpp"
+#include "Fierro_Optimization_Objective.hpp"
+#include "ROL_Elementwise_Reduce.hpp"
+#include "FEA_Module_SGH.h"
+#include "Explicit_Solver.h"
+
+class InternalEnergyMinimize_ShapeOpt : public FierroOptimizationObjective
+{
+typedef Tpetra::Map<>::local_ordinal_type LO;
+typedef Tpetra::Map<>::global_ordinal_type GO;
+typedef Tpetra::Map<>::node_type Node;
+typedef Tpetra::Map<LO, GO, Node> Map;
+typedef Tpetra::MultiVector<real_t, LO, GO, Node> MV;
+typedef ROL::Vector<real_t> V;
+typedef ROL::TpetraMultiVector<real_t, LO, GO, Node> ROL_MV;
+
+using traits = Kokkos::ViewTraits<LO*, Kokkos::LayoutLeft, void, void>;
+using array_layout    = typename traits::array_layout;
+using execution_space = typename traits::execution_space;
+using device_type     = typename traits::device_type;
+using memory_traits   = typename traits::memory_traits;
+using global_size_t   = Tpetra::global_size_t;
+
+typedef Kokkos::View<real_t*, Kokkos::LayoutRight, device_type, memory_traits> values_array;
+typedef Kokkos::View<GO*, array_layout, device_type, memory_traits> global_indices_array;
+typedef Kokkos::View<LO*, array_layout, device_type, memory_traits> indices_array;
+
+// typedef Kokkos::DualView<real_t**, Kokkos::LayoutLeft, device_type>::t_dev vec_array;
+typedef MV::dual_view_type::t_dev vec_array;
+typedef MV::dual_view_type::t_host host_vec_array;
+typedef Kokkos::View<const real_t**, array_layout, HostSpace, memory_traits> const_host_vec_array;
+typedef Kokkos::View<const real_t**, array_layout, device_type, memory_traits> const_vec_array;
+typedef MV::dual_view_type dual_vec_array;
+
+private:
+
+    Explicit_Solver* Explicit_Solver_Pointer_;
+    FEA_Module_SGH*  FEM_SGH_;
+    ROL::Ptr<ROL_MV> ROL_Force;
+    ROL::Ptr<ROL_MV> ROL_Velocities;
+    ROL::Ptr<ROL_MV> ROL_Gradients;
+    Teuchos::RCP<MV> previous_gradients;
+    real_t initial_internal_energy;
+    real_t previous_objective_accumulation, objective_sign;
+
+    bool first_init; //prevents ROL from calling init computation twice at start for the AL algorithm
+
+    /////////////////////////////////////////////////////////////////////////////
+    ///
+    /// \fn getVector
+    ///
+    /// \brief Retrieves ROL vector at desired location
+    ///
+    /// \param Pointer to desired ROL vector
+    ///
+    /// \return Returns ROL MV vector
+    ///
+    /////////////////////////////////////////////////////////////////////////////
+    ROL::Ptr<const MV> getVector(const V& x)
+    {
+        return dynamic_cast<const ROL_MV&>(x).getVector();
+    }
+
+    /////////////////////////////////////////////////////////////////////////////
+    ///
+    /// \fn getVector
+    ///
+    /// \brief Retrieves ROL vector at desired location
+    ///
+    /// \param Pointer to desired ROL vector
+    ///
+    /// \return Returns const ROL MV vector
+    ///
+    /////////////////////////////////////////////////////////////////////////////
+    ROL::Ptr<MV> getVector(V& x)
+    {
+        return dynamic_cast<ROL_MV&>(x).getVector();
+    }
+
+public:
+    int    last_comm_step, last_solve_step, current_step;
+    int num_dim;
+    size_t nvalid_modules;
+    size_t nlocal_nodes, num_corners, num_nodes_in_elem, rnum_elem;
+    DViewCArrayKokkos<double> node_mass, node_coords;
+    std::vector<FEA_MODULE_TYPE> valid_fea_modules; // modules that may interface with this objective function
+    FEA_MODULE_TYPE set_module_type;
+    // std::string my_fea_module = "SGH";
+
+    InternalEnergyMinimize_ShapeOpt(Explicit_Solver* Explicit_Solver_Pointer)
+    {
+        Explicit_Solver_Pointer_ = Explicit_Solver_Pointer;
+        first_init = false;
+        valid_fea_modules.push_back(FEA_MODULE_TYPE::SGH);
+        valid_fea_modules.push_back(FEA_MODULE_TYPE::Dynamic_Elasticity);
+        nvalid_modules = valid_fea_modules.size();
+        objective_sign = 1;
+        num_dim = Explicit_Solver_Pointer_->simparam.num_dims;
+        const Simulation_Parameters& simparam = Explicit_Solver_Pointer_->simparam;
+        for (const auto& fea_module : Explicit_Solver_Pointer_->fea_modules) {
+            for (int ivalid = 0; ivalid < nvalid_modules; ivalid++) {
+                if (fea_module->Module_Type == FEA_MODULE_TYPE::SGH) {
+                    FEM_SGH_ = dynamic_cast<FEA_Module_SGH*>(fea_module);
+                    set_module_type = FEA_MODULE_TYPE::SGH;
+                    node_mass = FEM_SGH_->node_mass;
+                    node_coords = FEM_SGH_->node_coords;
+                    nlocal_nodes = FEM_SGH_->nlocal_nodes;
+                    num_corners = FEM_SGH_->num_corners;
+                    num_nodes_in_elem = FEM_SGH_->num_nodes_in_elem;
+                    rnum_elem = FEM_SGH_->rnum_elem;
+                }
+            }
+        }
+        last_comm_step    = last_solve_step = -1;
+        current_step      = 0;
+        time_accumulation = true;
+        
+        previous_gradients = Teuchos::rcp(new MV(Explicit_Solver_Pointer_->map, num_dim));
+        if(Explicit_Solver_Pointer_->simparam.optimization_options.maximize_flag){
+            objective_sign = -1;
+        }
+        // ROL_Force = ROL::makePtr<ROL_MV>(FEM_->Global_Nodal_Forces);
+        if (set_module_type == FEA_MODULE_TYPE::SGH) {
+            ROL_Velocities = ROL::makePtr<ROL_MV>(FEM_SGH_->node_velocities_distributed);
+        }
+    }
+
+  /* --------------------------------------------------------------------------------------
+   Update solver state variables to synchronize with the current design variable vector, z
+  ----------------------------------------------------------------------------------------- */
+
+    void update(const ROL::Vector<real_t>& z, ROL::UpdateType type, int iter = -1)
+    {
+        if (set_module_type == FEA_MODULE_TYPE::SGH) {
+            update_sgh(z, type, iter);
+        }
+    }
+
+  /* --------------------------------------------------------------------------------------
+   Update solver state variables to synchronize with the current design variable vector, z
+  ----------------------------------------------------------------------------------------- */
+
+    void update_sgh(const ROL::Vector<real_t>& z, ROL::UpdateType type, int iter = -1)
+    {
+        // debug
+        std::ostream& out = std::cout;
+        Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
+        bool print_flag = false;
+        ROL::Ptr<const MV>   zp = getVector(z); //tpetra multivector wrapper on design vector
+        const_host_vec_array design_coordinates = zp->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+
+        if (type == ROL::UpdateType::Initial) {
+            if(first_init){
+                // This is the first call to update
+                if (Explicit_Solver_Pointer_->myrank == 0) {
+                    *fos << "called SGH Initial" << std::endl;
+                }
+
+                FEM_SGH_->comm_variables(zp);
+                FEM_SGH_->update_forward_solve(zp);
+                if(Explicit_Solver_Pointer_->myrank == 0){
+                    std::cout << "CURRENT TIME INTEGRAL OF INTERNAL ENERGY " << objective_accumulation << std::endl;
+                }
+                FEM_SGH_->compute_topology_optimization_adjoint_full(zp);
+                previous_objective_accumulation = objective_accumulation;
+                previous_gradients->assign(*(FEM_SGH_->cached_design_gradients_distributed));
+                // initial design density data was already communicated for ghost nodes in init_design()
+                // decide to output current optimization state
+                // FEM_SGH_->Explicit_Solver_Pointer_->write_outputs();
+            }
+            first_init = true;
+        }
+        else if (type == ROL::UpdateType::Accept) {
+            if (Explicit_Solver_Pointer_->myrank == 0) {
+                *fos << "called Accept" << std::endl;
+            }
+            
+            previous_objective_accumulation = objective_accumulation;
+            previous_gradients->assign(*(FEM_SGH_->cached_design_gradients_distributed));
+        }
+        else if (type == ROL::UpdateType::Revert) {
+            // u_ was set to u=S(x) during a trial update
+            // and has been rejected as the new iterate
+            // Revert to cached value
+            // This is a new value of x
+            // communicate density variables for ghosts
+            if (Explicit_Solver_Pointer_->myrank == 0) { *fos << "called Revert" << std::endl; }
+            objective_accumulation = previous_objective_accumulation;
+            FEM_SGH_->cached_design_gradients_distributed->assign(*previous_gradients);
+            if(Explicit_Solver_Pointer_->myrank == 0){
+                std::cout << "CURRENT TIME INTEGRAL OF INTERNAL ENERGY " << objective_accumulation << std::endl;
+            }
+            // FEM_SGH_->comm_variables(zp);
+            // // update deformation variables
+            // FEM_SGH_->update_forward_solve(zp);
+            // FEM_SGH_->compute_topology_optimization_adjoint_full(zp);
+        }
+        else if (type == ROL::UpdateType::Trial) {
+            // This is a new value of x
+            current_step++;
+            if(current_step%FEM_SGH_->simparam->optimization_options.optimization_output_freq==0){
+                print_flag = true;
+            }
+            if (Explicit_Solver_Pointer_->myrank == 0) {
+                *fos << "called Trial" << std::endl;
+            }
+            // communicate density variables for ghosts
+            FEM_SGH_->comm_variables(zp);
+            // update deformation variables
+            FEM_SGH_->update_forward_solve(zp, print_flag);
+            
+            if(Explicit_Solver_Pointer_->myrank == 0){
+                std::cout << "CURRENT TIME INTEGRAL OF INTERNAL ENERGY " << objective_accumulation << std::endl;
+            }
+            FEM_SGH_->compute_topology_optimization_adjoint_full(zp);
+            // decide to output current optimization state
+            // FEM_SGH_->Explicit_Solver_Pointer_->write_outputs();
+        }
+        else{ // ROL::UpdateType::Temp
+            // This is a new value of x used for,
+            // e.g., finite-difference checks
+            if (Explicit_Solver_Pointer_->myrank == 0) {
+                *fos << "called Temp" << std::endl;
+            }
+            FEM_SGH_->comm_variables(zp);
+            FEM_SGH_->update_forward_solve(zp);
+            if(Explicit_Solver_Pointer_->myrank == 0){
+                std::cout << "CURRENT TIME INTEGRAL OF INTERNAL ENERGY " << objective_accumulation << std::endl;
+            }
+            FEM_SGH_->compute_topology_optimization_adjoint_full(zp);
+        }
+    }
+
+  /* --------------------------------------------------------------------------------------
+   Compute time integral contribution for this objective function form
+  ----------------------------------------------------------------------------------------- */
+    void step_accumulation(const real_t& dt, const size_t& cycle, const size_t& rk_level) {
+        
+        const_vec_array current_elem_sie;
+        const_vec_array previous_elem_sie;
+        size_t nlocal_elem_non_overlapping = FEM_SGH_->Explicit_Solver_Pointer_->nlocal_elem_non_overlapping;
+        auto nodes_in_elem = FEM_SGH_->nodes_in_elem;
+        auto elem_mass = FEM_SGH_->elem_mass;
+        bool use_solve_checkpoints    = FEM_SGH_->simparam->optimization_options.use_solve_checkpoints;
+        if(use_solve_checkpoints){
+            current_elem_sie = FEM_SGH_->element_internal_energy_distributed->getLocalView<device_type>(Tpetra::Access::ReadOnly);
+            previous_elem_sie = FEM_SGH_->previous_element_internal_energy_distributed->getLocalView<device_type>(Tpetra::Access::ReadOnly);
+        }
+        else{
+            auto forward_solve_internal_energy_data = FEM_SGH_->forward_solve_internal_energy_data;
+            current_elem_sie = (*forward_solve_internal_energy_data)[cycle + 1]->getLocalView<device_type>(Tpetra::Access::ReadOnly);
+            previous_elem_sie = (*forward_solve_internal_energy_data)[cycle]->getLocalView<device_type>(Tpetra::Access::ReadOnly);
+        }
+
+        double IE_sum = 0.0;
+        double IE_loc_sum = 0.0;
+        // extensive KE
+        if(FEM_SGH_->simparam->optimization_options.optimization_objective_regions.size()){
+            int nobj_volumes = FEM_SGH_->simparam->optimization_options.optimization_objective_regions.size();
+            auto optimization_objective_regions = FEM_SGH_->simparam->optimization_options.optimization_objective_regions;
+            const_vec_array all_initial_node_coords = FEM_SGH_->all_initial_node_coords_distributed->getLocalView<device_type>(Tpetra::Access::ReadOnly);
+            FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+                int node_id;
+                double current_elem_coords[3];
+                bool contained = false;
+                current_elem_coords[0] = 0;
+                current_elem_coords[1] = 0;
+                current_elem_coords[2] = 0;
+                for(int inode=0; inode< num_nodes_in_elem; inode++){
+                    node_id = nodes_in_elem(elem_gid, inode);
+                    current_elem_coords[0] += all_initial_node_coords(node_id, 0)/num_nodes_in_elem;
+                    current_elem_coords[1] += all_initial_node_coords(node_id, 1)/num_nodes_in_elem;
+                    current_elem_coords[2] += all_initial_node_coords(node_id, 2)/num_nodes_in_elem;
+                }
+                for(int ivolume = 0; ivolume < nobj_volumes; ivolume++){
+                    if(optimization_objective_regions(ivolume).contains(current_elem_coords)){
+                        contained = true;
+                    }
+                }
+                if(contained){
+                    IE_loc_sum += elem_mass(elem_gid) *  0.5 * (current_elem_sie(elem_gid,0)+previous_elem_sie(elem_gid,0));
+                }
+
+                
+            }, IE_sum);
+        }
+        else{
+            FOR_REDUCE_SUM_CLASS(elem_gid, 0, nlocal_elem_non_overlapping, IE_loc_sum, {
+                IE_loc_sum += elem_mass(elem_gid) * 0.5 * (current_elem_sie(elem_gid,0)+previous_elem_sie(elem_gid,0));
+            }, IE_sum);
+        }
+
+        Kokkos::fence();
+        objective_accumulation += IE_sum * dt;
+    }
+
+  /* --------------------------------------------------------------------------------------
+   Update objective value with the current design variable vector, z
+  ----------------------------------------------------------------------------------------- */
+
+    real_t value(const ROL::Vector<real_t>& z, real_t& tol)
+    {
+        // std::cout << "Started obj value on task " <<FEM_->myrank  << std::endl;
+        ROL::Ptr<const MV> zp = getVector(z);
+        real_t c = 0.0;
+
+        // debug print
+        // std::ostream &out = std::cout;
+        // Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
+        // if(FEM_->myrank==0)
+        // *fos << "Value function z:" << std::endl;
+        // zp->describe(*fos,Teuchos::VERB_EXTREME);
+        // *fos << std::endl;
+        // std::fflush(stdout);
+
+        const_host_vec_array design_coordinates = zp->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+        // communicate ghosts and solve for nodal degrees of freedom as a function of the current design variables
+        /*
+        if(last_comm_step!=current_step){
+          FEM_->comm_variables(zp);
+          last_comm_step = current_step;
+        }
+
+        if(last_solve_step!=current_step){
+          //std::cout << "UPDATED velocities" << std::endl;
+          FEM_->update_linear_solve(zp);
+          last_solve_step = current_step;
+        }
+        */
+        // debug print of velocities
+        // std::ostream &out = std::cout;
+        // Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
+        // if(FEM_->myrank==0)
+        // *fos << "Displacement data :" << std::endl;
+        // FEM_->node_velocities_distributed->describe(*fos,Teuchos::VERB_EXTREME);
+        // *fos << std::endl;
+        // std::fflush(stdout);
+
+        std::cout.precision(10);
+
+        // std::cout << "Ended obj value on task " <<FEM_->myrank  << std::endl;
+        return objective_sign*objective_accumulation;
+    }
+
+  /* --------------------------------------------------------------------------------------
+   Update gradient vector (g) with the current design variable vector, z
+  ----------------------------------------------------------------------------------------- */
+
+    void gradient(ROL::Vector<real_t>& g, const ROL::Vector<real_t>& z, real_t& tol)
+    {
+        // std::cout << "Started obj gradient on task " <<FEM_->myrank  << std::endl;
+        // get Tpetra multivector pointer from the ROL vector
+        ROL::Ptr<const MV> zp = getVector(z); //pointer to design vector
+        ROL::Ptr<MV> gp = getVector(g); //pointer to gradient vector
+
+        // communicate ghosts and solve for nodal degrees of freedom as a function of the current design variables
+        // FEM_->gradient_print_sync=1;
+        // FEM_->gradient_print_sync=0;
+        // get local view of the data
+
+        if (set_module_type == FEA_MODULE_TYPE::SGH) {
+            FEM_SGH_->compute_shape_optimization_gradient_full(zp, gp);
+        }
+        gp->scale(objective_sign);
+        // debug print of gradient
+        // std::ostream &out = std::cout;
+        // Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
+        // if(FEM_->myrank==0)
+        // *fos << "Gradient data :" << std::endl;
+        // gp->describe(*fos,Teuchos::VERB_EXTREME);
+        // *fos << std::endl;
+        // std::fflush(stdout);
+        // for(int i = 0; i < FEM_->nlocal_nodes; i++){
+        // objective_gradients(i,0) *= -1;
+        // }
+
+        // std::cout << "Objective Gradient called"<< std::endl;
+        // debug print of design variables
+        // std::ostream &out = std::cout;
+        // Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
+        // if(FEM_->myrank==0)
+        // *fos << "Gradient data :" << std::endl;
+        // gp->describe(*fos,Teuchos::VERB_EXTREME);
+
+        // *fos << std::endl;
+        // std::fflush(stdout);
+        // std::cout << "ended obj gradient on task " <<FEM_->myrank  << std::endl;
+    }
+    
+  //contributes to rate of change of adjoint vector due to term with velocity gradient of objective
+    void sie_gradient_adjoint_contribution(vec_array& adjoint_rate_vector, const DViewCArrayKokkos<double>& node_mass,
+                                                const DViewCArrayKokkos<double>& elem_mass, const DViewCArrayKokkos<double>& node_vel,
+                                                const DViewCArrayKokkos<double>& node_coords, const DViewCArrayKokkos<double>& elem_sie,
+                                                const size_t& rk_level){
+        
+
+        FOR_ALL_CLASS(elem_gid, 0, rnum_elem, {
+            adjoint_rate_vector(elem_gid, 0) = elem_mass(elem_gid);
+        }); // end parallel for
+        Kokkos::fence();
+    }
+
+  //contribution to gradient tally from objective w.r.t design coordinate
+    void design_coordinate_gradient_term(vec_array& gradient_vector, const DViewCArrayKokkos<double>& node_mass,
+                               const DViewCArrayKokkos<double>& elem_mass, const DViewCArrayKokkos<double>& node_vel,
+                               const DViewCArrayKokkos<double>& node_coords, const DViewCArrayKokkos<double>& elem_sie,
+                               const size_t& rk_level, const real_t& global_dt = 0){
+        auto optimization_objective_regions = FEM_SGH_->simparam->optimization_options.optimization_objective_regions;
+        auto nodes_in_elem = FEM_SGH_->nodes_in_elem;
+        auto corner_value_storage = FEM_SGH_->corner_value_storage;
+        auto corners_in_node = FEM_SGH_->corners_in_node;
+        auto num_corners_in_node = FEM_SGH_->num_corners_in_node;
+    }
+
+};
+
+#endif // end header guard
diff --git a/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Kinetic_Energy_Minimize.h b/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Kinetic_Energy_Minimize.h
index 932cebb3c..a7e2ee6e8 100644
--- a/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Kinetic_Energy_Minimize.h
+++ b/src/Parallel-Solvers/Parallel-Explicit/Topology_Optimization/Kinetic_Energy_Minimize.h
@@ -399,7 +399,7 @@ typedef MV::dual_view_type dual_vec_array;
             int nobj_volumes = FEM_SGH_->simparam->optimization_options.optimization_objective_regions.size();
             auto optimization_objective_regions = FEM_SGH_->simparam->optimization_options.optimization_objective_regions;
             const_vec_array all_initial_node_coords = FEM_SGH_->all_initial_node_coords_distributed->getLocalView<device_type>(Tpetra::Access::ReadOnly);
-            REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+            FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
                 double ke = 0;
                 double current_node_coords[3];
                 bool contained = false;
@@ -428,7 +428,7 @@ typedef MV::dual_view_type dual_vec_array;
             }, KE_sum);
         }
         else{
-            REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
+            FOR_REDUCE_SUM_CLASS(node_gid, 0, nlocal_nodes, KE_loc_sum, {
                 double ke = 0;
                 for (size_t dim = 0; dim < num_dim; dim++) {
                     // midpoint integration approximation
diff --git a/src/Parallel-Solvers/Simulation_Parameters/Input_Options.h b/src/Parallel-Solvers/Simulation_Parameters/Input_Options.h
index bc16008ad..c00dc9804 100644
--- a/src/Parallel-Solvers/Simulation_Parameters/Input_Options.h
+++ b/src/Parallel-Solvers/Simulation_Parameters/Input_Options.h
@@ -22,6 +22,7 @@ struct Input_Options : Yaml::ValidatedYaml, Yaml::DerivedFields {
   int p_order = 2;
   double unit_scaling = 1.0;
   bool topology_optimization_restart = false;
+  bool shape_optimization_restart = false;
 
   ELEMENT_TYPE element_type = ELEMENT_TYPE::hex8;
   bool zero_index_base = false;
@@ -90,4 +91,4 @@ struct Input_Options : Yaml::ValidatedYaml, Yaml::DerivedFields {
   }
 };
 IMPL_YAML_SERIALIZABLE_FOR(Input_Options, mesh_file_name, mesh_file_format, element_type, zero_index_base, p_order, unit_scaling,
-topology_optimization_restart)
+topology_optimization_restart, shape_optimization_restart)
diff --git a/src/Parallel-Solvers/Simulation_Parameters/Optimization_Options.h b/src/Parallel-Solvers/Simulation_Parameters/Optimization_Options.h
index df209b375..b57a5dd20 100644
--- a/src/Parallel-Solvers/Simulation_Parameters/Optimization_Options.h
+++ b/src/Parallel-Solvers/Simulation_Parameters/Optimization_Options.h
@@ -12,20 +12,23 @@ SERIALIZABLE_ENUM(FUNCTION_TYPE,
 //  VECTOR_INEQUALITY_CONSTRAINT
 )
 
-SERIALIZABLE_ENUM(TO_MODULE_TYPE,
-  Kinetic_Energy_Minimize,
-  Internal_Energy_Minimize,
+SERIALIZABLE_ENUM(OPTIMIZATION_MODULE_TYPE,
+  Kinetic_Energy_Minimize_TopOpt,
+  Internal_Energy_Minimize_TopOpt,
+  Internal_Energy_Minimize_ShapeOpt,
   Multi_Objective,
-  Heat_Capacity_Potential_Minimize,
-  Strain_Energy_Minimize,
-  Mass_Constraint,
-  Moment_of_Inertia_Constraint,
-  Center_of_Mass_Constraint,
-  Heat_Capacity_Potential_Constraint,
+  Heat_Capacity_Potential_Minimize_TopOpt,
+  Strain_Energy_Minimize_TopOpt,
+  Mass_Constraint_TopOpt,
+  Mass_Constraint_ShapeOpt,
+  Moment_of_Inertia_Constraint_TopOpt,
+  Moment_of_Inertia_Constraint_ShapeOpt,
+  Center_of_Mass_Constraint_TopOpt,
+  Heat_Capacity_Potential_Constraint_TopOpt,
   MULTI_OBJECTIVE_TERM,
-  Thermo_Elastic_Strain_Energy_Minimize,
-  Strain_Energy_Constraint,
-  Displacement_Constraint
+  Thermo_Elastic_Strain_Energy_Minimize_TopOpt,
+  Strain_Energy_Constraint_TopOpt,
+  Displacement_Constraint_TopOpt
 )
 
 
@@ -164,6 +167,7 @@ struct Optimization_Options: Yaml::DerivedFields {
   bool use_gradient_tally = false;                            //tallies gradient in tandem with the time sequence solving for the adjoint vectors
   bool optimization_parameters_xml_file = false;
   std::string xml_parameters_file_name = "optimization_parameters.xml";
+  real_t max_coord_move_length = 1;
 
   MULTI_OBJECTIVE_STRUCTURE multi_objective_structure = MULTI_OBJECTIVE_STRUCTURE::linear;
   std::vector<MultiObjectiveModule> multi_objective_modules;
@@ -191,7 +195,7 @@ IMPL_YAML_SERIALIZABLE_FOR(Optimization_Options,
   optimization_process, optimization_objective, 
   constraints, method_of_moving_asymptotes, volume_bound_constraints, objective_regions,
   simp_penalty_power, density_epsilon, thick_condition_boundary,
-  optimization_output_freq, density_filter, minimum_density, maximum_density,
+  optimization_output_freq, density_filter, minimum_density, maximum_density, max_coord_move_length,
   multi_objective_modules, multi_objective_structure, density_filter, retain_outer_shell,
   variable_outer_shell, shell_density, objective_normalization_constant,
   num_solve_checkpoints, use_solve_checkpoints, use_gradient_tally, disable_forward_solve_output,
diff --git a/src/Parallel-Solvers/Simulation_Parameters/Simulation_Parameters.h b/src/Parallel-Solvers/Simulation_Parameters/Simulation_Parameters.h
index 4635dc99c..7130fd015 100644
--- a/src/Parallel-Solvers/Simulation_Parameters/Simulation_Parameters.h
+++ b/src/Parallel-Solvers/Simulation_Parameters/Simulation_Parameters.h
@@ -41,11 +41,11 @@ struct Simulation_Parameters
     bool restart_file = false;
 
     std::set<FEA_MODULE_TYPE> fea_module_must_read;
-    //list of TO functions needed by problem
-    std::vector<TO_MODULE_TYPE> TO_Module_List;
-    std::vector<FUNCTION_TYPE> TO_Function_Type;
-    std::vector<int> TO_Module_My_FEA_Module;
-    std::vector<std::vector<int>> FEA_Module_My_TO_Modules;
+    //list of Optimization functions needed by problem
+    std::vector<OPTIMIZATION_MODULE_TYPE> Optimization_Module_List;
+    std::vector<FUNCTION_TYPE> Optimization_Function_Type;
+    std::vector<int> Optimization_Module_My_FEA_Module;
+    std::vector<std::vector<int>> FEA_Module_My_Optimization_Modules;
     std::vector<std::vector<double>> Function_Arguments;
 
     std::vector<int> Multi_Objective_Modules;
@@ -96,45 +96,64 @@ struct Simulation_Parameters
     }
 
     void derive_objective_module() {
-        switch (optimization_options.optimization_objective) {
-        case OPTIMIZATION_OBJECTIVE::minimize_compliance:
-            add_TO_module(TO_MODULE_TYPE::Strain_Energy_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            optimization_options.normalized_objective = true;
-            break;
-        case OPTIMIZATION_OBJECTIVE::minimize_thermal_resistance:
-            add_TO_module(TO_MODULE_TYPE::Heat_Capacity_Potential_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            optimization_options.normalized_objective = true;
-            break;
-        case OPTIMIZATION_OBJECTIVE::minimize_kinetic_energy:
-            add_TO_module(TO_MODULE_TYPE::Kinetic_Energy_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            break;
-        case OPTIMIZATION_OBJECTIVE::maximize_compliance:
-            add_TO_module(TO_MODULE_TYPE::Strain_Energy_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            optimization_options.normalized_objective = true;
-            optimization_options.maximize_flag = true;
-            break;
-        case OPTIMIZATION_OBJECTIVE::maximize_thermal_resistance:
-            add_TO_module(TO_MODULE_TYPE::Heat_Capacity_Potential_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            optimization_options.normalized_objective = true;
-            optimization_options.maximize_flag = true;
-            break;
-        case OPTIMIZATION_OBJECTIVE::maximize_kinetic_energy:
-            add_TO_module(TO_MODULE_TYPE::Kinetic_Energy_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            optimization_options.maximize_flag = true;
-            break;
-        case OPTIMIZATION_OBJECTIVE::minimize_internal_energy:
-            add_TO_module(TO_MODULE_TYPE::Internal_Energy_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            break;
-        case OPTIMIZATION_OBJECTIVE::maximize_internal_energy:
-            add_TO_module(TO_MODULE_TYPE::Internal_Energy_Minimize, FUNCTION_TYPE::OBJECTIVE, {});
-            optimization_options.maximize_flag = true;
-            break;
-        case OPTIMIZATION_OBJECTIVE::multi_objective:
-            add_TO_module(TO_MODULE_TYPE::Multi_Objective, FUNCTION_TYPE::OBJECTIVE, {});
-            derive_multi_objectives();
-            break;
-        default:
-            throw Yaml::ConfigurationException("Unsupported optimization objective " + to_string(optimization_options.optimization_objective));
+        if(topology_optimization_on){
+            switch (optimization_options.optimization_objective) {
+            case OPTIMIZATION_OBJECTIVE::minimize_compliance:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Strain_Energy_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                optimization_options.normalized_objective = true;
+                break;
+            case OPTIMIZATION_OBJECTIVE::minimize_thermal_resistance:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                optimization_options.normalized_objective = true;
+                break;
+            case OPTIMIZATION_OBJECTIVE::minimize_kinetic_energy:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Kinetic_Energy_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                break;
+            case OPTIMIZATION_OBJECTIVE::maximize_compliance:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Strain_Energy_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                optimization_options.normalized_objective = true;
+                optimization_options.maximize_flag = true;
+                break;
+            case OPTIMIZATION_OBJECTIVE::maximize_thermal_resistance:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                optimization_options.normalized_objective = true;
+                optimization_options.maximize_flag = true;
+                break;
+            case OPTIMIZATION_OBJECTIVE::maximize_kinetic_energy:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Kinetic_Energy_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                optimization_options.maximize_flag = true;
+                break;
+            case OPTIMIZATION_OBJECTIVE::minimize_internal_energy:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                break;
+            case OPTIMIZATION_OBJECTIVE::maximize_internal_energy:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_TopOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                optimization_options.maximize_flag = true;
+                break;
+            case OPTIMIZATION_OBJECTIVE::multi_objective:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Multi_Objective, FUNCTION_TYPE::OBJECTIVE, {});
+                derive_multi_objectives();
+                break;
+            default:
+                throw Yaml::ConfigurationException("Unsupported optimization objective " + to_string(optimization_options.optimization_objective));
+            }
+        }
+        else if(shape_optimization_on){
+            switch (optimization_options.optimization_objective) {
+            case OPTIMIZATION_OBJECTIVE::minimize_internal_energy:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_ShapeOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                break;
+            case OPTIMIZATION_OBJECTIVE::maximize_internal_energy:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_ShapeOpt, FUNCTION_TYPE::OBJECTIVE, {});
+                optimization_options.maximize_flag = true;
+                break;
+            case OPTIMIZATION_OBJECTIVE::multi_objective:
+                add_optimization_module(OPTIMIZATION_MODULE_TYPE::Multi_Objective, FUNCTION_TYPE::OBJECTIVE, {});
+                derive_multi_objectives();
+                break;
+            default:
+                throw Yaml::ConfigurationException("Unsupported optimization objective " + to_string(optimization_options.optimization_objective));
+            }
         }
     }
     void derive_multi_objectives() {
@@ -142,20 +161,20 @@ struct Simulation_Parameters
             return;
         
         for (auto mod : optimization_options.multi_objective_modules) {
-            TO_MODULE_TYPE to_type;
+            OPTIMIZATION_MODULE_TYPE to_type;
             switch (mod.type) {
             case OPTIMIZATION_OBJECTIVE::minimize_compliance:
-                to_type = TO_MODULE_TYPE::Strain_Energy_Minimize;
+                to_type = OPTIMIZATION_MODULE_TYPE::Strain_Energy_Minimize_TopOpt;
                 break;
             case OPTIMIZATION_OBJECTIVE::minimize_thermal_resistance:
-                to_type = TO_MODULE_TYPE::Heat_Capacity_Potential_Minimize;
+                to_type = OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Minimize_TopOpt;
                 break;
             default:
                 throw Yaml::ConfigurationException("Unsupported sub-objective " + to_string(mod.type));
             }
-            Multi_Objective_Modules.push_back(TO_Module_List.size());
+            Multi_Objective_Modules.push_back(Optimization_Module_List.size());
             Multi_Objective_Weights.push_back(mod.weight_coefficient);
-            add_TO_module(to_type, FUNCTION_TYPE::MULTI_OBJECTIVE_TERM, {});
+            add_optimization_module(to_type, FUNCTION_TYPE::MULTI_OBJECTIVE_TERM, {});
         }
     }
     void derive_constraint_modules() {
@@ -168,38 +187,60 @@ struct Simulation_Parameters
             default:
                 throw Yaml::ConfigurationException("Unsupported relation " + to_string(constraint.relation));
             }
-
-            switch (constraint.type) {
-            case CONSTRAINT_TYPE::mass:
-                add_TO_module(
-                    TO_MODULE_TYPE::Mass_Constraint, 
-                    f_type, 
-                    {constraint.value}
-                );
-                break;
-            case CONSTRAINT_TYPE::moment_of_inertia:
-                add_TO_module(
-                    TO_MODULE_TYPE::Moment_of_Inertia_Constraint, 
-                    f_type, 
-                    {constraint.value, (double)component_to_int(constraint.component.value())}
-                );
-                break;
-            case CONSTRAINT_TYPE::center_of_mass:
-                add_TO_module(
-                    TO_MODULE_TYPE::Center_of_Mass_Constraint, 
-                    f_type, 
-                    {constraint.value, (double)component_to_int(constraint.component.value())}
-                );
-                break;
-            case CONSTRAINT_TYPE::displacement:
-                add_TO_module(
-                    TO_MODULE_TYPE::Displacement_Constraint, 
-                    f_type, 
-                    {constraint.value}
-                );
-                break;
-            default:
-                throw Yaml::ConfigurationException("Unsupported constraint type " + to_string(constraint.type));
+            
+            if(topology_optimization_on){
+                switch (constraint.type) {
+                case CONSTRAINT_TYPE::mass:
+                    add_optimization_module(
+                        OPTIMIZATION_MODULE_TYPE::Mass_Constraint_TopOpt, 
+                        f_type, 
+                        {constraint.value}
+                    );
+                    break;
+                case CONSTRAINT_TYPE::moment_of_inertia:
+                    add_optimization_module(
+                        OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt, 
+                        f_type, 
+                        {constraint.value, (double)component_to_int(constraint.component.value())}
+                    );
+                    break;
+                case CONSTRAINT_TYPE::center_of_mass:
+                    add_optimization_module(
+                        OPTIMIZATION_MODULE_TYPE::Center_of_Mass_Constraint_TopOpt, 
+                        f_type, 
+                        {constraint.value, (double)component_to_int(constraint.component.value())}
+                    );
+                    break;
+                case CONSTRAINT_TYPE::displacement:
+                    add_optimization_module(
+                        OPTIMIZATION_MODULE_TYPE::Displacement_Constraint_TopOpt, 
+                        f_type, 
+                        {constraint.value}
+                    );
+                    break;
+                default:
+                    throw Yaml::ConfigurationException("Unsupported constraint type " + to_string(constraint.type));
+                }
+            }
+            else if(shape_optimization_on){
+                switch (constraint.type) {
+                case CONSTRAINT_TYPE::mass:
+                    add_optimization_module(
+                        OPTIMIZATION_MODULE_TYPE::Mass_Constraint_ShapeOpt, 
+                        f_type, 
+                        {constraint.value}
+                    );
+                    break;
+                case CONSTRAINT_TYPE::moment_of_inertia:
+                    add_optimization_module(
+                        OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_ShapeOpt, 
+                        f_type, 
+                        {constraint.value, (double)component_to_int(constraint.component.value())}
+                    );
+                    break;
+                default:
+                    throw Yaml::ConfigurationException("Unsupported constraint type " + to_string(constraint.type));
+                }
             }
         }
     }
@@ -209,22 +250,22 @@ struct Simulation_Parameters
         shape_optimization_on    = optimization_options.optimization_process == OPTIMIZATION_PROCESS::shape_optimization;
     }
 
-    void map_TO_to_FEA() {
+    void map_Opt_to_FEA() {
         // Now we allocate the vectors for each of the currently identified modules.
-        TO_Module_My_FEA_Module.resize(TO_Module_List.size());
-        FEA_Module_My_TO_Modules.resize(fea_module_parameters.size());
+        Optimization_Module_My_FEA_Module.resize(Optimization_Module_List.size());
+        FEA_Module_My_Optimization_Modules.resize(fea_module_parameters.size());
         
         // Finally we can set up the maps.
         // 
         // TODO: This should really use two `std::map<int, std::set<int>>`s
         // instead of this vector stuff.
-        for (size_t to_index = 0; to_index < TO_Module_List.size(); to_index++) {
-            auto to_module = TO_Module_List[to_index];
-            auto fea_index = find_TO_module_dependency(to_module);
+        for (size_t to_index = 0; to_index < Optimization_Module_List.size(); to_index++) {
+            auto to_module = Optimization_Module_List[to_index];
+            auto fea_index = find_Optimization_module_dependency(to_module);
             if (!fea_index.has_value())
                 continue;
-            TO_Module_My_FEA_Module[to_index] = fea_index.value();
-            FEA_Module_My_TO_Modules[fea_index.value()].push_back(to_index);
+            Optimization_Module_My_FEA_Module[to_index] = fea_index.value();
+            FEA_Module_My_Optimization_Modules[fea_index.value()].push_back(to_index);
         }
     }
 
@@ -250,7 +291,7 @@ struct Simulation_Parameters
             derive_objective_module();
             derive_constraint_modules();
         }
-        map_TO_to_FEA();
+        map_Opt_to_FEA();
 
         mtr::from_vector(mat_fill, regions);
         mtr::from_vector(material, materials);
@@ -309,13 +350,13 @@ struct Simulation_Parameters
             throw Yaml::ConfigurationException("Specify exactly one of `input_options` and `mesh_generation_options`.");
         validate_element_type();
         // Check that the FEA module dependencies are satisfied.
-        for (auto to_module : TO_Module_List)
-            find_TO_module_dependency(to_module);
+        for (auto to_module : Optimization_Module_List)
+            find_Optimization_module_dependency(to_module);
     }
 
-    void add_TO_module(TO_MODULE_TYPE type, FUNCTION_TYPE function_type, std::vector<double> arguments) {
-        TO_Module_List.push_back(type);
-        TO_Function_Type.push_back(function_type);
+    void add_optimization_module(OPTIMIZATION_MODULE_TYPE type, FUNCTION_TYPE function_type, std::vector<double> arguments) {
+        Optimization_Module_List.push_back(type);
+        Optimization_Function_Type.push_back(function_type);
         Function_Arguments.push_back(arguments);
     }
 
@@ -326,19 +367,20 @@ struct Simulation_Parameters
      * If dependencies are not satisfied, an exception will be thrown.
      * If there are no dependencies, std::nullopt will be returned.
      */
-    std::optional<size_t> find_TO_module_dependency(TO_MODULE_TYPE type) {
-        const static std::unordered_map<TO_MODULE_TYPE, std::vector<FEA_MODULE_TYPE>> map {
-            {TO_MODULE_TYPE::Heat_Capacity_Potential_Minimize,      {FEA_MODULE_TYPE::Heat_Conduction}},
-            {TO_MODULE_TYPE::Heat_Capacity_Potential_Constraint,    {FEA_MODULE_TYPE::Heat_Conduction}},
-            {TO_MODULE_TYPE::Thermo_Elastic_Strain_Energy_Minimize, {FEA_MODULE_TYPE::Heat_Conduction}},
-            {TO_MODULE_TYPE::Mass_Constraint,                       {FEA_MODULE_TYPE::Inertial       }},
-            {TO_MODULE_TYPE::Center_of_Mass_Constraint,             {FEA_MODULE_TYPE::Inertial       }},
-            {TO_MODULE_TYPE::Moment_of_Inertia_Constraint,          {FEA_MODULE_TYPE::Inertial       }},
-            {TO_MODULE_TYPE::Strain_Energy_Minimize,                {FEA_MODULE_TYPE::Elasticity     }},
-            {TO_MODULE_TYPE::Displacement_Constraint,               {FEA_MODULE_TYPE::Elasticity     }},
-            {TO_MODULE_TYPE::Strain_Energy_Constraint,              {FEA_MODULE_TYPE::Elasticity     }},
-            {TO_MODULE_TYPE::Internal_Energy_Minimize,              {FEA_MODULE_TYPE::SGH            }},
-            {TO_MODULE_TYPE::Kinetic_Energy_Minimize,               {FEA_MODULE_TYPE::SGH, FEA_MODULE_TYPE::Dynamic_Elasticity}},
+    std::optional<size_t> find_Optimization_module_dependency(OPTIMIZATION_MODULE_TYPE type) {
+        const static std::unordered_map<OPTIMIZATION_MODULE_TYPE, std::vector<FEA_MODULE_TYPE>> map {
+            {OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Minimize_TopOpt,      {FEA_MODULE_TYPE::Heat_Conduction}},
+            {OPTIMIZATION_MODULE_TYPE::Heat_Capacity_Potential_Constraint_TopOpt,    {FEA_MODULE_TYPE::Heat_Conduction}},
+            {OPTIMIZATION_MODULE_TYPE::Thermo_Elastic_Strain_Energy_Minimize_TopOpt, {FEA_MODULE_TYPE::Heat_Conduction}},
+            {OPTIMIZATION_MODULE_TYPE::Mass_Constraint_TopOpt,                       {FEA_MODULE_TYPE::Inertial       }},
+            {OPTIMIZATION_MODULE_TYPE::Center_of_Mass_Constraint_TopOpt,             {FEA_MODULE_TYPE::Inertial       }},
+            {OPTIMIZATION_MODULE_TYPE::Moment_of_Inertia_Constraint_TopOpt,          {FEA_MODULE_TYPE::Inertial       }},
+            {OPTIMIZATION_MODULE_TYPE::Strain_Energy_Minimize_TopOpt,                {FEA_MODULE_TYPE::Elasticity     }},
+            {OPTIMIZATION_MODULE_TYPE::Displacement_Constraint_TopOpt,               {FEA_MODULE_TYPE::Elasticity     }},
+            {OPTIMIZATION_MODULE_TYPE::Strain_Energy_Constraint_TopOpt,              {FEA_MODULE_TYPE::Elasticity     }},
+            {OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_TopOpt,              {FEA_MODULE_TYPE::SGH            }},
+            {OPTIMIZATION_MODULE_TYPE::Internal_Energy_Minimize_ShapeOpt,            {FEA_MODULE_TYPE::SGH            }},
+            {OPTIMIZATION_MODULE_TYPE::Kinetic_Energy_Minimize_TopOpt,               {FEA_MODULE_TYPE::SGH, FEA_MODULE_TYPE::Dynamic_Elasticity}},
         };
         if (map.count(type) == 0)
             return {};
diff --git a/src/Parallel-Solvers/Solver.cpp b/src/Parallel-Solvers/Solver.cpp
index f0f94699a..4d8771bd6 100644
--- a/src/Parallel-Solvers/Solver.cpp
+++ b/src/Parallel-Solvers/Solver.cpp
@@ -1031,6 +1031,7 @@ void Solver::read_mesh_vtk(const char* MESH)
     real_t unit_scaling                  = input_options.unit_scaling;
     bool   zero_index_base               = input_options.zero_index_base;
     bool   topology_optimization_restart = input_options.topology_optimization_restart;
+    bool   shape_optimization_restart    = input_options.shape_optimization_restart;
 
     CArrayKokkos<char, array_layout, HostSpace, memory_traits> read_buffer;
 
@@ -1128,7 +1129,7 @@ void Solver::read_mesh_vtk(const char* MESH)
             buffer_iterations++;
         }
 
-        // read coords, also density if restarting
+        // read coords
         read_index_start = 0;
         for (buffer_iteration = 0; buffer_iteration < buffer_iterations; buffer_iteration++)
         {
@@ -1205,7 +1206,6 @@ void Solver::read_mesh_vtk(const char* MESH)
                         dof_value = atof(&read_buffer(scan_loop, 2, 0));
                         node_coords(node_rid, 2) = dof_value * unit_scaling;
                     }
-                    // extract density if restarting
                 }
             }
             read_index_start += BUFFER_LINES;
@@ -1759,6 +1759,161 @@ void Solver::read_mesh_vtk(const char* MESH)
             } // end if
         } // end if(myrank==0)
     }
+
+    //If restarting a topology optimization run, obtain nodal design density data here
+    if(shape_optimization_restart){
+        design_node_coords_distributed = Teuchos::rcp(new MV(map, num_dim));
+        host_vec_array design_node_coords = design_node_coords_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+        if (myrank == 0)
+        {
+            bool found = false;
+            while (found == false&&in->good()) {
+                std::getline(*in, read_line);
+                //std::cout << read_line << std::endl;
+                line_parse.clear();
+                line_parse.str(read_line);
+
+                //stop when the design_density string is reached
+                while (!line_parse.eof()){
+                    line_parse >> substring;
+                    //std::cout << substring << std::endl;
+                    if(!substring.compare("design_coordinates")){
+                        found = true;
+                    }
+                } //while
+
+            } // end while
+
+            if (!found){
+                throw std::runtime_error("ERROR: Failed to find design_coordinates");
+            } // end if
+
+            //skip "LOOKUP_TABLE default" line
+            std::getline(*in, read_line);
+        } // end if(myrank==0)
+        
+        //read in density of each node
+        // allocate read buffer
+        words_per_line = num_dim;
+        read_buffer = CArrayKokkos<char, array_layout, HostSpace, memory_traits>(BUFFER_LINES, words_per_line, MAX_WORD);
+
+        dof_limit = num_nodes;
+        buffer_iterations = dof_limit / BUFFER_LINES;
+        if (dof_limit % BUFFER_LINES != 0)
+        {
+            buffer_iterations++;
+        }
+
+        // read densities
+        read_index_start = 0;
+        for (buffer_iteration = 0; buffer_iteration < buffer_iterations; buffer_iteration++)
+        {
+            // pack buffer on rank 0
+            if (myrank == 0 && buffer_iteration < buffer_iterations - 1)
+            {
+                for (buffer_loop = 0; buffer_loop < BUFFER_LINES; buffer_loop++)
+                {
+                    getline(*in, read_line);
+                    line_parse.clear();
+                    line_parse.str(read_line);
+
+                    for (int iword = 0; iword < words_per_line; iword++)
+                    {
+                        // read portions of the line into the substring variable
+                        line_parse >> substring;
+                        // debug print
+                        // std::cout<<" "<< substring <<std::endl;
+                        // assign the substring variable as a word of the read buffer
+                        strcpy(&read_buffer(buffer_loop, iword, 0), substring.c_str());
+                    }
+                }
+            }
+            else if (myrank == 0)
+            {
+                buffer_loop = 0;
+                while (buffer_iteration * BUFFER_LINES + buffer_loop < num_nodes) {
+                    getline(*in, read_line);
+                    line_parse.clear();
+                    line_parse.str(read_line);
+                    for (int iword = 0; iword < words_per_line; iword++)
+                    {
+                        // read portions of the line into the substring variable
+                        line_parse >> substring;
+                        // debug print
+                        // std::cout<<" "<< substring <<std::endl;
+                        // assign the substring variable as a word of the read buffer
+                        strcpy(&read_buffer(buffer_loop, iword, 0), substring.c_str());
+                    }
+                    buffer_loop++;
+                }
+            }
+
+            // broadcast buffer to all ranks; each rank will determine which nodes in the buffer belong
+            MPI_Bcast(read_buffer.pointer(), BUFFER_LINES * words_per_line * MAX_WORD, MPI_CHAR, 0, world);
+            // broadcast how many nodes were read into this buffer iteration
+            MPI_Bcast(&buffer_loop, 1, MPI_INT, 0, world);
+
+            // debug_print
+            // std::cout << "NODE BUFFER LOOP IS: " << buffer_loop << std::endl;
+            // for(int iprint=0; iprint < buffer_loop; iprint++)
+            // std::cout<<"buffer packing: " << std::string(&read_buffer(iprint,0,0)) << std::endl;
+            // return;
+
+            // determine which data to store in the swage mesh members (the local node data)
+            // loop through read buffer
+            for (scan_loop = 0; scan_loop < buffer_loop; scan_loop++)
+            {
+                // set global node id (ensight specific order)
+                node_gid = read_index_start + scan_loop;
+                // let map decide if this node id belongs locally; if yes store data
+                if (map->isNodeGlobalElement(node_gid))
+                {
+                    // set local node index in this mpi rank
+                    node_rid = map->getLocalElement(node_gid);
+                    // extract nodal position from the read buffer
+                    // for tecplot format this is the three coords in the same line
+                    dof_value = atof(&read_buffer(scan_loop, 0, 0));
+                    design_node_coords(node_rid, 0) = dof_value * unit_scaling;
+                    dof_value = atof(&read_buffer(scan_loop, 1, 0));
+                    design_node_coords(node_rid, 1) = dof_value * unit_scaling;
+                    if(num_dim==3){
+                        dof_value = atof(&read_buffer(scan_loop, 2, 0));
+                        design_node_coords(node_rid, 2) = dof_value * unit_scaling; 
+                    }
+                }
+            }
+            read_index_start += BUFFER_LINES;
+        }
+
+        //Find initial objective value to normalize by
+        if (myrank == 0 && simparam.optimization_options.normalized_objective)
+        {
+            bool found = false;
+            while (found == false&&in->good()) {
+                std::getline(*in, read_line);
+                //std::cout << read_line << std::endl;
+                line_parse.clear();
+                line_parse.str(read_line);
+
+                //stop when the design_density string is reached
+                while (!line_parse.eof()){
+                    line_parse >> substring;
+                    //std::cout << substring << std::endl;
+                    if(!substring.compare("Objective_Normalization_Constant")){
+                        found = true;
+                        line_parse >> substring;
+                        simparam.optimization_options.objective_normalization_constant = stod(substring);
+                    }
+                } //while
+
+            } // end while
+
+            if (!found){
+                throw std::runtime_error("ERROR: Failed to find initial objective value for restart");
+            } // end if
+        } // end if(myrank==0)
+    }
+
     // Close mesh input file
     if (myrank == 0)
     {
diff --git a/src/Parallel-Solvers/Solver.h b/src/Parallel-Solvers/Solver.h
index c3ced4c40..a19e18af0 100644
--- a/src/Parallel-Solvers/Solver.h
+++ b/src/Parallel-Solvers/Solver.h
@@ -225,8 +225,8 @@ class Solver
     Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> map; // map of node indices
     Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> sorted_map; // sorted contiguous map of node indices
     Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> ghost_node_map; // map of node indices with ghosts on each rank
-    Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> all_node_map; // map of node indices with ghosts on each rank
-    Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> nonoverlap_element_node_map; // map of node indices with ghosts on each rank
+    Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> all_node_map; // map of local + ghost node indices on each rank
+    Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> nonoverlap_element_node_map; // map of node indices belonging to unique element map
     Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> element_map; // non overlapping map of elements owned by each rank used in reduction ops
     Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> all_element_map; // overlapping map of elements connected to the local nodes in each rank
     Teuchos::RCP<Tpetra::Map<LO, GO, node_type>> sorted_element_map; // sorted contiguous map of element indices owned by each rank used in parallel IO
@@ -240,12 +240,15 @@ class Solver
     Teuchos::RCP<MV>       all_initial_node_coords_distributed;
     Teuchos::RCP<MV>       all_node_coords_distributed;
     Teuchos::RCP<MV>       design_node_densities_distributed;
+    Teuchos::RCP<MV>       design_node_coords_distributed;
     Teuchos::RCP<MV>       filtered_node_densities_distributed;
     Teuchos::RCP<const MV> test_node_densities_distributed;
     Teuchos::RCP<MV>       all_node_densities_distributed;
+    Teuchos::RCP<MV>       all_design_node_coords_distributed;
     Teuchos::RCP<MV>       all_filtered_node_densities_distributed;
     Teuchos::RCP<MV>       lower_bound_node_densities_distributed, all_lower_bound_node_densities_distributed;
     Teuchos::RCP<MV>       upper_bound_node_densities_distributed;
+    Teuchos::RCP<MV>       lower_bound_node_coordinates_distributed, upper_bound_node_coordinates_distributed;
     Teuchos::RCP<MV>       Global_Element_Densities_Upper_Bound;
     Teuchos::RCP<MV>       Global_Element_Densities_Lower_Bound;
     Teuchos::RCP<MV>       Global_Element_Densities;
diff --git a/src/Voxelizer/src/stl-to-voxelvtk.cpp b/src/Voxelizer/src/stl-to-voxelvtk.cpp
index 5afccb08c..5236fe60d 100644
--- a/src/Voxelizer/src/stl-to-voxelvtk.cpp
+++ b/src/Voxelizer/src/stl-to-voxelvtk.cpp
@@ -400,7 +400,7 @@ void main_function(CArray<bool> &gridOUTPUT, int &gridX, int &gridY, int &gridZ,
     float meshZmin;
     
     // Global maximum x-direction
-    REDUCE_MAX(i,0,n_facets,meshXmax, {
+    FOR_REDUCE_MAX(i,0,n_facets,meshXmax, {
         if (v1X(i) > meshXmax | v2X(i) > meshXmax | v3X(i) > meshXmax) {
             if (v1X(i) > v2X(i) && v1X(i) > v3X(i)) {
                 meshXmax = v1X(i);
@@ -413,7 +413,7 @@ void main_function(CArray<bool> &gridOUTPUT, int &gridX, int &gridY, int &gridZ,
     }, meshXmax);
     
     // Global minimum x-direction
-    REDUCE_MIN(i,0,n_facets,meshXmin, {
+    FOR_REDUCE_MIN(i,0,n_facets,meshXmin, {
         if (v1X(i) < meshXmin | v2X(i) < meshXmin | v3X(i) < meshXmin) {
             if (v1X(i) < v2X(i) && v1X(i) < v3X(i)) {
                 meshXmin = v1X(i);
@@ -426,7 +426,7 @@ void main_function(CArray<bool> &gridOUTPUT, int &gridX, int &gridY, int &gridZ,
     }, meshXmin);
 
     // Global maximum y-direction
-    REDUCE_MAX(i,0,n_facets,meshYmax, {
+    FOR_REDUCE_MAX(i,0,n_facets,meshYmax, {
         if (v1Y(i) > meshYmax | v2Y(i) > meshYmax | v3Y(i) > meshYmax) {
             if (v1Y(i) > v2Y(i) && v1Y(i) > v3Y(i)) {
                 meshYmax = v1Y(i);
@@ -439,7 +439,7 @@ void main_function(CArray<bool> &gridOUTPUT, int &gridX, int &gridY, int &gridZ,
     }, meshYmax);
 
     // Global minimum y-direction
-    REDUCE_MIN(i,0,n_facets,meshYmin, {
+    FOR_REDUCE_MIN(i,0,n_facets,meshYmin, {
         if (v1Y(i) < meshYmin | v2Y(i) < meshYmin | v3Y(i) < meshYmin) {
             if (v1Y(i) < v2Y(i) && v1Y(i) < v3Y(i)) {
                 meshYmin = v1Y(i);
@@ -452,7 +452,7 @@ void main_function(CArray<bool> &gridOUTPUT, int &gridX, int &gridY, int &gridZ,
     }, meshYmin);
 
     // Global maximum z-direction
-    REDUCE_MAX(i,0,n_facets,meshZmax, {
+    FOR_REDUCE_MAX(i,0,n_facets,meshZmax, {
         if (v1Z(i) > meshZmax | v2Z(i) > meshZmax | v3Z(i) > meshZmax) {
             if (v1Z(i) > v2Z(i) && v1Z(i) > v3Z(i)) {
                 meshZmax = v1Z(i);
@@ -465,7 +465,7 @@ void main_function(CArray<bool> &gridOUTPUT, int &gridX, int &gridY, int &gridZ,
     }, meshZmax);
 
     // Global minimum z-direction
-    REDUCE_MIN(i,0,n_facets,meshZmin, {
+    FOR_REDUCE_MIN(i,0,n_facets,meshZmin, {
         if (v1Z(i) < meshZmin | v2Z(i) < meshZmin | v3Z(i) < meshZmin) {
             if (v1Z(i) < v2Z(i) && v1Z(i) < v3Z(i)) {
                 meshZmin = v1Z(i);