diff --git a/data/simple_cube.g b/data/simple_cube.g
new file mode 100644
index 000000000..07ca3fec8
Binary files /dev/null and b/data/simple_cube.g differ
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 1f8dda0c6..252623489 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -3,4 +3,5 @@ add_subdirectory(team7)
 add_subdirectory(complex_team7)
 add_subdirectory(closed_coil)
 add_subdirectory(open_coil)
-add_subdirectory(magnetostatic)
\ No newline at end of file
+add_subdirectory(magnetostatic)
+add_subdirectory(thermal_expansion)
diff --git a/examples/thermal_expansion/CMakeLists.txt b/examples/thermal_expansion/CMakeLists.txt
new file mode 100644
index 000000000..1c612f572
--- /dev/null
+++ b/examples/thermal_expansion/CMakeLists.txt
@@ -0,0 +1,23 @@
+file(GLOB_RECURSE example_src "*.cpp")
+
+file(GLOB_RECURSE test_src_files "${PROJECT_SOURCE_DIR}/src/*.h" "${PROJECT_SOURCE_DIR}/src/*.hpp" "${PROJECT_SOURCE_DIR}/src/*.cpp")
+
+set (${PROJECT_NAME}_INCLUDE_DIRS "")
+foreach (_srcFile ${test_src_files})
+    get_filename_component(_dir ${_srcFile} PATH)
+    list (APPEND ${PROJECT_NAME}_INCLUDE_DIRS ${_dir})
+endforeach()
+list (REMOVE_DUPLICATES ${PROJECT_NAME}_INCLUDE_DIRS)
+
+
+add_executable(thermalexpansion ${example_src})
+add_compile_options(thermalexpansion ${BUILD_TYPE_COMPILER_FLAGS})
+target_include_directories(thermalexpansion PUBLIC ${MFEM_COMMON_INCLUDES} ${MFEM_INCLUDE_DIRS})
+target_include_directories(thermalexpansion PUBLIC ${${PROJECT_NAME}_INCLUDE_DIRS})
+target_include_directories(thermalexpansion PUBLIC ${PROJECT_SOURCE_DIR}/data/)
+
+set_property(TARGET thermalexpansion PROPERTY CXX_STANDARD 17)
+
+target_link_libraries(thermalexpansion ${GTEST_LIBRARY} ${GTEST_MAIN_LIBRARY} pthread)
+target_link_libraries(thermalexpansion ${CMAKE_LIBRARY_OUTPUT_DIRECTORY}/lib${PROJECT_NAME}.so)
+target_link_libraries(thermalexpansion ${MFEM_LIBRARIES} ${MFEM_COMMON_LIBRARY} -lrt)
diff --git a/examples/thermal_expansion/Main.cpp b/examples/thermal_expansion/Main.cpp
new file mode 100644
index 000000000..cc96e3feb
--- /dev/null
+++ b/examples/thermal_expansion/Main.cpp
@@ -0,0 +1,120 @@
+#include "hephaestus.hpp"
+
+const char *DATA_DIR = "../../data/";
+
+hephaestus::Coefficients defineCoefficients(){
+    
+  hephaestus::Coefficients coefficients;
+
+  coefficients.scalars.Register("thermal_expansion_coef",
+                                new mfem::ConstantCoefficient(1e-04),
+                                true);
+
+  coefficients.scalars.Register("thermal_conductivity",
+                                new mfem::ConstantCoefficient(300),
+                                true);
+
+  coefficients.scalars.Register("lame_param",
+                                new mfem::ConstantCoefficient(10),
+                                true);
+
+  coefficients.scalars.Register("shear_modulus",
+                                new mfem::ConstantCoefficient(10),
+                                true);                        
+
+  coefficients.scalars.Register("stress_free_temp",
+                                new mfem::ConstantCoefficient(250),
+                                true);
+
+  coefficients.scalars.Register("zero",
+                                new mfem::ConstantCoefficient(0.0),
+                                true);
+                                                                            
+
+  return coefficients;
+}
+
+hephaestus::Outputs defineOutputs() {
+  // std::map<std::string, mfem::DataCollection *> data_collections;
+  // data_collections["ParaViewDataCollection"] =
+  //     new mfem::ParaViewDataCollection("ThermalExpansionExample");
+  hephaestus::Outputs outputs;
+  outputs.Register("Paraview", new mfem::ParaViewDataCollection("ThermalExpansionExample"), true);
+  return outputs;
+}
+
+
+hephaestus::BCMap defineBoundaries() {
+
+  hephaestus::BCMap boundaries;
+  mfem::ConstantCoefficient *cold = new mfem::ConstantCoefficient(300.00);
+  mfem::ConstantCoefficient *hot = new mfem::ConstantCoefficient(500.00);
+  // mfem::VectorConstantCoefficient *zero = new mfem::VectorConstantCoefficient(mfem::Vector({0, 0, 0}));
+  mfem::Array<int> therm_bound_one_arr(1);
+  mfem::Array<int> therm_bound_two_arr(1);
+  mfem::Array<int> fixed_disp_arr(1);
+
+  therm_bound_one_arr = 1;
+  therm_bound_two_arr = 2;
+  fixed_disp_arr = 1;
+
+  boundaries.Register("thermal_boundary_one", new hephaestus::DirichletBC("temperature", therm_bound_one_arr, cold), true);
+  boundaries.Register("thermal_boundary_two", new hephaestus::DirichletBC("temperature", therm_bound_two_arr, hot), true);
+  boundaries.Register("fixed_displacement", new hephaestus::VectorDirichletBC("displacement", fixed_disp_arr,  new mfem::VectorConstantCoefficient(mfem::Vector({0, 0, 0})), nullptr, hephaestus::VectorDirichletBC::APPLY_TYPE::STANDARD), true);
+
+  return boundaries;
+}
+
+
+int main(int argc, char *argv[]) {
+  mfem::OptionsParser args(argc, argv);
+  args.AddOption(&DATA_DIR, "-dataDir", "--data_directory",
+                 "Directory storing input data for tests.");
+  args.Parse();
+  MPI_Init(&argc, &argv);
+
+  // Create Formulation
+  hephaestus::SteadyStateProblemBuilder *problem_builder = new hephaestus::ThermalExpansionFormulation();
+  // Set Mesh
+  mfem::Mesh mesh((std::string(DATA_DIR) + std::string("./simple_cube.g")).c_str(), 1,
+                  1);
+  std::shared_ptr<mfem::ParMesh> pmesh =
+      std::make_shared<mfem::ParMesh>(mfem::ParMesh(MPI_COMM_WORLD, mesh));
+  problem_builder->SetMesh(pmesh);
+  problem_builder->AddFESpace(std::string("H1"), std::string("H1_3D_P1"));
+  problem_builder->AddFESpace(std::string("H1_3"), std::string("H1_3D_P1"), 3, mfem::Ordering::byVDIM);
+  
+  problem_builder->AddGridFunction(std::string("temperature"),
+                                   std::string("H1"));
+  problem_builder->AddGridFunction(std::string("displacement"),
+                                   std::string("H1_3"));
+
+  hephaestus::Coefficients coefficients = defineCoefficients();
+  problem_builder->SetCoefficients(coefficients);
+
+  hephaestus::Outputs outputs = defineOutputs();
+  problem_builder->SetOutputs(outputs);
+
+  hephaestus::BCMap boundaries = defineBoundaries();
+  problem_builder->SetBoundaryConditions(boundaries);
+  
+  hephaestus::InputParameters solver_options;
+  solver_options.SetParam("Tolerance", float(1.0e-5));
+  solver_options.SetParam("MaxIter", (unsigned int)1000);
+  solver_options.SetParam("PrintLevel", 0);
+  problem_builder->SetSolverOptions(solver_options);
+
+  hephaestus::ProblemBuildSequencer sequencer(problem_builder);
+  sequencer.ConstructEquationSystemProblem();
+  std::unique_ptr<hephaestus::SteadyStateProblem> problem =
+      problem_builder->ReturnProblem();
+  hephaestus::InputParameters exec_params;
+  exec_params.SetParam("Problem", problem.get());
+  hephaestus::SteadyExecutioner *executioner =
+      new hephaestus::SteadyExecutioner(exec_params);
+
+  mfem::out << "Created executioner";
+  executioner->Execute();
+
+  MPI_Finalize();
+}
diff --git a/src/CMakeCache.txt b/src/CMakeCache.txt
new file mode 100644
index 000000000..97d0e3467
--- /dev/null
+++ b/src/CMakeCache.txt
@@ -0,0 +1,371 @@
+# This is the CMakeCache file.
+# For build in directory: /home/bill/Projects/hephaestus/src
+# It was generated by CMake: /usr/bin/cmake
+# You can edit this file to change values found and used by cmake.
+# If you do not want to change any of the values, simply exit the editor.
+# If you do want to change a value, simply edit, save, and exit the editor.
+# The syntax for the file is as follows:
+# KEY:TYPE=VALUE
+# KEY is the name of a variable in the cache.
+# TYPE is a hint to GUIs for the type of VALUE, DO NOT EDIT TYPE!.
+# VALUE is the current value for the KEY.
+
+########################
+# EXTERNAL cache entries
+########################
+
+//The directory containing a CMake configuration file for Boost.
+Boost_DIR:PATH=/usr/lib/x86_64-linux-gnu/cmake/Boost-1.74.0
+
+//Path to a file.
+Boost_INCLUDE_DIR:PATH=/usr/include
+
+//Path to a program.
+CMAKE_ADDR2LINE:FILEPATH=/usr/bin/addr2line
+
+//Path to a program.
+CMAKE_AR:FILEPATH=/usr/bin/ar
+
+//Choose the type of build, options are: None Debug Release RelWithDebInfo
+// MinSizeRel ...
+CMAKE_BUILD_TYPE:STRING=
+
+//Enable/Disable color output during build.
+CMAKE_COLOR_MAKEFILE:BOOL=ON
+
+//CXX compiler
+CMAKE_CXX_COMPILER:FILEPATH=/usr/bin/c++
+
+//A wrapper around 'ar' adding the appropriate '--plugin' option
+// for the GCC compiler
+CMAKE_CXX_COMPILER_AR:FILEPATH=/usr/bin/gcc-ar-12
+
+//A wrapper around 'ranlib' adding the appropriate '--plugin' option
+// for the GCC compiler
+CMAKE_CXX_COMPILER_RANLIB:FILEPATH=/usr/bin/gcc-ranlib-12
+
+//Flags used by the CXX compiler during all build types.
+CMAKE_CXX_FLAGS:STRING=
+
+//Flags used by the CXX compiler during DEBUG builds.
+CMAKE_CXX_FLAGS_DEBUG:STRING=-g
+
+//Flags used by the CXX compiler during MINSIZEREL builds.
+CMAKE_CXX_FLAGS_MINSIZEREL:STRING=-Os -DNDEBUG
+
+//Flags used by the CXX compiler during RELEASE builds.
+CMAKE_CXX_FLAGS_RELEASE:STRING=-O3 -DNDEBUG
+
+//Flags used by the CXX compiler during RELWITHDEBINFO builds.
+CMAKE_CXX_FLAGS_RELWITHDEBINFO:STRING=-O2 -g -DNDEBUG
+
+//Path to a program.
+CMAKE_DLLTOOL:FILEPATH=CMAKE_DLLTOOL-NOTFOUND
+
+//Flags used by the linker during all build types.
+CMAKE_EXE_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during DEBUG builds.
+CMAKE_EXE_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during MINSIZEREL builds.
+CMAKE_EXE_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during RELEASE builds.
+CMAKE_EXE_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during RELWITHDEBINFO builds.
+CMAKE_EXE_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//Enable/Disable output of compile commands during generation.
+CMAKE_EXPORT_COMPILE_COMMANDS:BOOL=
+
+//Value Computed by CMake.
+CMAKE_FIND_PACKAGE_REDIRECTS_DIR:STATIC=/home/bill/Projects/hephaestus/src/CMakeFiles/pkgRedirects
+
+//Install path prefix, prepended onto install directories.
+CMAKE_INSTALL_PREFIX:PATH=/usr/local
+
+//Path to a program.
+CMAKE_LINKER:FILEPATH=/usr/bin/ld
+
+//Path to a program.
+CMAKE_MAKE_PROGRAM:FILEPATH=/usr/bin/gmake
+
+//Flags used by the linker during the creation of modules during
+// all build types.
+CMAKE_MODULE_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during the creation of modules during
+// DEBUG builds.
+CMAKE_MODULE_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during the creation of modules during
+// MINSIZEREL builds.
+CMAKE_MODULE_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during the creation of modules during
+// RELEASE builds.
+CMAKE_MODULE_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during the creation of modules during
+// RELWITHDEBINFO builds.
+CMAKE_MODULE_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//Path to a program.
+CMAKE_NM:FILEPATH=/usr/bin/nm
+
+//Path to a program.
+CMAKE_OBJCOPY:FILEPATH=/usr/bin/objcopy
+
+//Path to a program.
+CMAKE_OBJDUMP:FILEPATH=/usr/bin/objdump
+
+//Value Computed by CMake
+CMAKE_PROJECT_DESCRIPTION:STATIC=
+
+//Value Computed by CMake
+CMAKE_PROJECT_HOMEPAGE_URL:STATIC=
+
+//Value Computed by CMake
+CMAKE_PROJECT_NAME:STATIC=hephaestus
+
+//Value Computed by CMake
+CMAKE_PROJECT_VERSION:STATIC=0.0.0
+
+//Value Computed by CMake
+CMAKE_PROJECT_VERSION_MAJOR:STATIC=0
+
+//Value Computed by CMake
+CMAKE_PROJECT_VERSION_MINOR:STATIC=0
+
+//Value Computed by CMake
+CMAKE_PROJECT_VERSION_PATCH:STATIC=0
+
+//Value Computed by CMake
+CMAKE_PROJECT_VERSION_TWEAK:STATIC=
+
+//Path to a program.
+CMAKE_RANLIB:FILEPATH=/usr/bin/ranlib
+
+//Path to a program.
+CMAKE_READELF:FILEPATH=/usr/bin/readelf
+
+//Flags used by the linker during the creation of shared libraries
+// during all build types.
+CMAKE_SHARED_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during the creation of shared libraries
+// during DEBUG builds.
+CMAKE_SHARED_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during the creation of shared libraries
+// during MINSIZEREL builds.
+CMAKE_SHARED_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during the creation of shared libraries
+// during RELEASE builds.
+CMAKE_SHARED_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during the creation of shared libraries
+// during RELWITHDEBINFO builds.
+CMAKE_SHARED_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//If set, runtime paths are not added when installing shared libraries,
+// but are added when building.
+CMAKE_SKIP_INSTALL_RPATH:BOOL=NO
+
+//If set, runtime paths are not added when using shared libraries.
+CMAKE_SKIP_RPATH:BOOL=NO
+
+//Flags used by the linker during the creation of static libraries
+// during all build types.
+CMAKE_STATIC_LINKER_FLAGS:STRING=
+
+//Flags used by the linker during the creation of static libraries
+// during DEBUG builds.
+CMAKE_STATIC_LINKER_FLAGS_DEBUG:STRING=
+
+//Flags used by the linker during the creation of static libraries
+// during MINSIZEREL builds.
+CMAKE_STATIC_LINKER_FLAGS_MINSIZEREL:STRING=
+
+//Flags used by the linker during the creation of static libraries
+// during RELEASE builds.
+CMAKE_STATIC_LINKER_FLAGS_RELEASE:STRING=
+
+//Flags used by the linker during the creation of static libraries
+// during RELWITHDEBINFO builds.
+CMAKE_STATIC_LINKER_FLAGS_RELWITHDEBINFO:STRING=
+
+//Path to a program.
+CMAKE_STRIP:FILEPATH=/usr/bin/strip
+
+//If this value is on, makefiles will be generated without the
+// .SILENT directive, and all commands will be echoed to the console
+// during the make.  This is useful for debugging only. With Visual
+// Studio IDE projects all commands are done without /nologo.
+CMAKE_VERBOSE_MAKEFILE:BOOL=FALSE
+
+//Path to the MFEM common miniapp headers.
+MFEM_COMMON_INCLUDES:PATH=/home/bill/Projects/hephaestus/../mfem/miniapps/common
+
+//Path to the MFEM build or install prefix.
+MFEM_DIR:PATH=/home/bill/Projects/hephaestus/../mfem_build
+
+//The directory containing a CMake configuration file for boost_headers.
+boost_headers_DIR:PATH=/usr/lib/x86_64-linux-gnu/cmake/boost_headers-1.74.0
+
+//Value Computed by CMake
+hephaestus_BINARY_DIR:STATIC=/home/bill/Projects/hephaestus/src
+
+//Value Computed by CMake
+hephaestus_IS_TOP_LEVEL:STATIC=ON
+
+//Value Computed by CMake
+hephaestus_SOURCE_DIR:STATIC=/home/bill/Projects/hephaestus
+
+//The directory containing a CMake configuration file for mfem.
+mfem_DIR:PATH=mfem_DIR-NOTFOUND
+
+
+########################
+# INTERNAL cache entries
+########################
+
+//ADVANCED property for variable: Boost_DIR
+Boost_DIR-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_ADDR2LINE
+CMAKE_ADDR2LINE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_AR
+CMAKE_AR-ADVANCED:INTERNAL=1
+//This is the directory where this CMakeCache.txt was created
+CMAKE_CACHEFILE_DIR:INTERNAL=/home/bill/Projects/hephaestus/src
+//Major version of cmake used to create the current loaded cache
+CMAKE_CACHE_MAJOR_VERSION:INTERNAL=3
+//Minor version of cmake used to create the current loaded cache
+CMAKE_CACHE_MINOR_VERSION:INTERNAL=25
+//Patch version of cmake used to create the current loaded cache
+CMAKE_CACHE_PATCH_VERSION:INTERNAL=1
+//ADVANCED property for variable: CMAKE_COLOR_MAKEFILE
+CMAKE_COLOR_MAKEFILE-ADVANCED:INTERNAL=1
+//Path to CMake executable.
+CMAKE_COMMAND:INTERNAL=/usr/bin/cmake
+//Path to cpack program executable.
+CMAKE_CPACK_COMMAND:INTERNAL=/usr/bin/cpack
+//Path to ctest program executable.
+CMAKE_CTEST_COMMAND:INTERNAL=/usr/bin/ctest
+//ADVANCED property for variable: CMAKE_CXX_COMPILER
+CMAKE_CXX_COMPILER-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_COMPILER_AR
+CMAKE_CXX_COMPILER_AR-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_COMPILER_RANLIB
+CMAKE_CXX_COMPILER_RANLIB-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS
+CMAKE_CXX_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_DEBUG
+CMAKE_CXX_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_MINSIZEREL
+CMAKE_CXX_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_RELEASE
+CMAKE_CXX_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_CXX_FLAGS_RELWITHDEBINFO
+CMAKE_CXX_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_DLLTOOL
+CMAKE_DLLTOOL-ADVANCED:INTERNAL=1
+//Executable file format
+CMAKE_EXECUTABLE_FORMAT:INTERNAL=ELF
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS
+CMAKE_EXE_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_DEBUG
+CMAKE_EXE_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_MINSIZEREL
+CMAKE_EXE_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_RELEASE
+CMAKE_EXE_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXE_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_EXE_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_EXPORT_COMPILE_COMMANDS
+CMAKE_EXPORT_COMPILE_COMMANDS-ADVANCED:INTERNAL=1
+//Name of external makefile project generator.
+CMAKE_EXTRA_GENERATOR:INTERNAL=
+//Name of generator.
+CMAKE_GENERATOR:INTERNAL=Unix Makefiles
+//Generator instance identifier.
+CMAKE_GENERATOR_INSTANCE:INTERNAL=
+//Name of generator platform.
+CMAKE_GENERATOR_PLATFORM:INTERNAL=
+//Name of generator toolset.
+CMAKE_GENERATOR_TOOLSET:INTERNAL=
+//Source directory with the top level CMakeLists.txt file for this
+// project
+CMAKE_HOME_DIRECTORY:INTERNAL=/home/bill/Projects/hephaestus
+//Install .so files without execute permission.
+CMAKE_INSTALL_SO_NO_EXE:INTERNAL=1
+//ADVANCED property for variable: CMAKE_LINKER
+CMAKE_LINKER-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MAKE_PROGRAM
+CMAKE_MAKE_PROGRAM-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS
+CMAKE_MODULE_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_DEBUG
+CMAKE_MODULE_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_MINSIZEREL
+CMAKE_MODULE_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_RELEASE
+CMAKE_MODULE_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_MODULE_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_MODULE_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_NM
+CMAKE_NM-ADVANCED:INTERNAL=1
+//number of local generators
+CMAKE_NUMBER_OF_MAKEFILES:INTERNAL=1
+//ADVANCED property for variable: CMAKE_OBJCOPY
+CMAKE_OBJCOPY-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_OBJDUMP
+CMAKE_OBJDUMP-ADVANCED:INTERNAL=1
+//Platform information initialized
+CMAKE_PLATFORM_INFO_INITIALIZED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_RANLIB
+CMAKE_RANLIB-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_READELF
+CMAKE_READELF-ADVANCED:INTERNAL=1
+//Path to CMake installation.
+CMAKE_ROOT:INTERNAL=/usr/share/cmake-3.25
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS
+CMAKE_SHARED_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_DEBUG
+CMAKE_SHARED_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_MINSIZEREL
+CMAKE_SHARED_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_RELEASE
+CMAKE_SHARED_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SHARED_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_SHARED_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SKIP_INSTALL_RPATH
+CMAKE_SKIP_INSTALL_RPATH-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_SKIP_RPATH
+CMAKE_SKIP_RPATH-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS
+CMAKE_STATIC_LINKER_FLAGS-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_DEBUG
+CMAKE_STATIC_LINKER_FLAGS_DEBUG-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_MINSIZEREL
+CMAKE_STATIC_LINKER_FLAGS_MINSIZEREL-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_RELEASE
+CMAKE_STATIC_LINKER_FLAGS_RELEASE-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STATIC_LINKER_FLAGS_RELWITHDEBINFO
+CMAKE_STATIC_LINKER_FLAGS_RELWITHDEBINFO-ADVANCED:INTERNAL=1
+//ADVANCED property for variable: CMAKE_STRIP
+CMAKE_STRIP-ADVANCED:INTERNAL=1
+//uname command
+CMAKE_UNAME:INTERNAL=/usr/bin/uname
+//ADVANCED property for variable: CMAKE_VERBOSE_MAKEFILE
+CMAKE_VERBOSE_MAKEFILE-ADVANCED:INTERNAL=1
+//Details about finding Boost
+FIND_PACKAGE_MESSAGE_DETAILS_Boost:INTERNAL=[/usr/lib/x86_64-linux-gnu/cmake/Boost-1.74.0/BoostConfig.cmake][c ][v1.74.0()]
+//linker supports push/pop state
+_CMAKE_LINKER_PUSHPOP_STATE_SUPPORTED:INTERNAL=TRUE
+//ADVANCED property for variable: boost_headers_DIR
+boost_headers_DIR-ADVANCED:INTERNAL=1
+
diff --git a/src/boundary_conditions/Essential/dirichlet_bc.cpp b/src/boundary_conditions/Essential/dirichlet_bc.cpp
new file mode 100644
index 000000000..207b971e1
--- /dev/null
+++ b/src/boundary_conditions/Essential/dirichlet_bc.cpp
@@ -0,0 +1,22 @@
+#include "dirichlet_bc.hpp"
+
+namespace hephaestus {
+
+DirichletBC::DirichletBC(const std::string &name_,
+                                         mfem::Array<int> bdr_attributes_)
+    : EssentialBC(name_, bdr_attributes_) {}
+
+DirichletBC::DirichletBC(const std::string &name_,
+                                         mfem::Array<int> bdr_attributes_,
+                                         mfem::ConstantCoefficient *coeff_,
+                                         mfem::ConstantCoefficient *coeff_im_)
+    : EssentialBC(name_, bdr_attributes_), coeff(coeff_), coeff_im(coeff_im_) {}
+
+void DirichletBC::applyBC(mfem::GridFunction &gridfunc,
+                                  mfem::Mesh *mesh_) {
+  mfem::Array<int> ess_bdrs(mesh_->bdr_attributes.Max());
+  ess_bdrs = this->getMarkers(*mesh_);
+  gridfunc.ProjectBdrCoefficient(*(this->coeff), ess_bdrs);
+}
+
+} // namespace hephaestus
diff --git a/src/boundary_conditions/Essential/dirichlet_bc.hpp b/src/boundary_conditions/Essential/dirichlet_bc.hpp
new file mode 100644
index 000000000..dc48f77fa
--- /dev/null
+++ b/src/boundary_conditions/Essential/dirichlet_bc.hpp
@@ -0,0 +1,22 @@
+#pragma once
+#include "essential_bc_base.hpp"
+
+namespace hephaestus{
+
+class DirichletBC : public EssentialBC{
+public:
+  DirichletBC(const std::string &name_,
+                      mfem::Array<int> bdr_attributes_);
+                      
+  DirichletBC(const std::string &name_,
+                      mfem::Array<int> bdr_attributes_,
+                      mfem::ConstantCoefficient *coeff_,
+                      mfem::ConstantCoefficient *coeff_im_ = nullptr);
+
+  virtual void applyBC(mfem::GridFunction &gridfunc,
+                       mfem::Mesh *mesh_) override;
+
+  mfem::ConstantCoefficient *coeff;
+  mfem::ConstantCoefficient *coeff_im;
+};
+} // namespace hephaestus
\ No newline at end of file
diff --git a/src/boundary_conditions/Essential/essential_bcs.hpp b/src/boundary_conditions/Essential/essential_bcs.hpp
index 7b2f406d0..245899c32 100644
--- a/src/boundary_conditions/Essential/essential_bcs.hpp
+++ b/src/boundary_conditions/Essential/essential_bcs.hpp
@@ -1,3 +1,4 @@
 #pragma once
 #include "function_dirichlet_bc.hpp"
 #include "vector_dirichlet_bc.hpp"
+#include "dirichlet_bc.hpp"
diff --git a/src/boundary_conditions/Essential/vector_dirichlet_bc.hpp b/src/boundary_conditions/Essential/vector_dirichlet_bc.hpp
index bcd3ecd0d..196dca574 100644
--- a/src/boundary_conditions/Essential/vector_dirichlet_bc.hpp
+++ b/src/boundary_conditions/Essential/vector_dirichlet_bc.hpp
@@ -5,7 +5,7 @@ namespace hephaestus {
 
 class VectorDirichletBC : public EssentialBC {
 
-protected:
+public:
   enum APPLY_TYPE
   {
     STANDARD,
diff --git a/src/factory/factory.hpp b/src/factory/factory.hpp
index 5cab10a0e..ad83e532f 100644
--- a/src/factory/factory.hpp
+++ b/src/factory/factory.hpp
@@ -7,6 +7,7 @@
 #include "e_formulation.hpp"
 #include "eb_dual_formulation.hpp"
 #include "h_formulation.hpp"
+#include "thermal_expansion_formulation.hpp"
 #include "inputs.hpp"
 #include "magnetostatic_formulation.hpp"
 
diff --git a/src/formulations/ThermalExpansion/thermal_expansion_formulation.cpp b/src/formulations/ThermalExpansion/thermal_expansion_formulation.cpp
new file mode 100644
index 000000000..a505ad05b
--- /dev/null
+++ b/src/formulations/ThermalExpansion/thermal_expansion_formulation.cpp
@@ -0,0 +1,253 @@
+#include "thermal_expansion_formulation.hpp"
+
+namespace hephaestus {
+
+ThermalExpansionFormulation::ThermalExpansionFormulation() : SteadyStateFormulation() {
+
+  temp_var_name = std::string("temperature");
+  displacement_var_name = std::string("displacement");
+  lame_param_coef_name = std::string("lame_param");
+  shear_modulus_coef_name = std::string("shear_modulus");
+  thermal_expansion_coef_name = std::string("thermal_expansion_coef");
+  thermal_conductivity_coef_name = std::string("thermal_conductivity");
+  stress_free_temp_coef_name = std::string("stress_free_temp");
+  thermal_expansion_bilin_coef_name = std::string("thermal_expansion_bilin_coef");
+  thermal_expansion_lin_coef_name = std::string("thermal_expansion_lin_coef");
+  zero_coef_name = std::string("zero");
+}
+
+void ThermalExpansionFormulation::RegisterCoefficients() {
+  std::cout << "REGISTER" << std::endl;
+  hephaestus::Coefficients &coefficients = this->GetProblem()->coefficients;
+
+  if (!coefficients.scalars.Has(lame_param_coef_name)) {
+    MFEM_ABORT(lame_param_coef_name + " coefficient not found.");
+  }
+
+  if (!coefficients.scalars.Has(shear_modulus_coef_name)) {
+    MFEM_ABORT(shear_modulus_coef_name + " coefficient not found.");
+  }
+
+  if (!coefficients.scalars.Has(stress_free_temp_coef_name)) {
+    MFEM_ABORT(stress_free_temp_coef_name + " coefficient not found.");
+  }
+
+  if (!coefficients.scalars.Has(thermal_expansion_coef_name)) {
+    MFEM_ABORT(thermal_expansion_coef_name + " coefficient not found.");
+  }
+
+  if (!coefficients.scalars.Has(thermal_conductivity_coef_name)) {
+    MFEM_ABORT(thermal_conductivity_coef_name + " coefficient not found.");
+  }
+
+  if (!coefficients.scalars.Has(zero_coef_name)) {
+    MFEM_ABORT(zero_coef_name + " coefficient not found.");
+  }
+
+  if (!coefficients.scalars.Has("materialTerm")) {
+    mfem::SumCoefficient* materialTerm = new mfem::SumCoefficient(*(coefficients.scalars.Get(lame_param_coef_name)), *(coefficients.scalars.Get(shear_modulus_coef_name)), 3, 2);
+    coefficients.scalars.Register("materialTerm", materialTerm, true);
+  }
+
+  if (!coefficients.scalars.Has("thexpStressFreeTemp")) {
+    mfem::SumCoefficient* thexpStressFreeTemp = new mfem::SumCoefficient(*(coefficients.scalars.Get(thermal_expansion_coef_name)), *(coefficients.scalars.Get(stress_free_temp_coef_name)));
+    coefficients.scalars.Register("thexpStressFreeTemp", thexpStressFreeTemp, true);
+  }
+
+  if (!coefficients.scalars.Has("bilinearFormCoefPositive")) {
+    mfem::ProductCoefficient *bilinearFormCoefPositive = new mfem::ProductCoefficient(*(coefficients.scalars.Get(thermal_expansion_coef_name)),
+                                                    *(coefficients.scalars.Get("materialTerm")));
+    coefficients.scalars.Register("bilinearFormCoefPositive", bilinearFormCoefPositive, true);
+  }
+                                            
+  if (!coefficients.scalars.Has(thermal_expansion_bilin_coef_name)) {
+    coefficients.scalars.Register(thermal_expansion_bilin_coef_name,
+                                    new mfem::ProductCoefficient(-1.0, *(coefficients.scalars.Get("bilinearFormCoefPositive"))), true);
+  }
+
+  if (!coefficients.scalars.Has(thermal_expansion_lin_coef_name)) {
+    coefficients.scalars.Register(thermal_expansion_lin_coef_name,
+                                    new mfem::ProductCoefficient(*(coefficients.scalars.Get("thexpStressFreeTemp")),
+                                    *(coefficients.scalars.Get("materialTerm"))), true);
+  }
+}
+
+void ThermalExpansionFormulation::ConstructEquationSystem() {
+  hephaestus::InputParameters weak_form_params;
+  weak_form_params.SetParam("TempVarName", temp_var_name);
+  weak_form_params.SetParam("DisplacementVarName", displacement_var_name);
+  weak_form_params.SetParam("StressFreeTempCoefName", stress_free_temp_coef_name);
+  weak_form_params.SetParam("LameParamCoefName", lame_param_coef_name);
+  weak_form_params.SetParam("ShearModulusCoefName", shear_modulus_coef_name);
+  weak_form_params.SetParam("ThermalExpansionCoefName", thermal_expansion_coef_name);
+  weak_form_params.SetParam("ThermalConductivityCoefName", thermal_conductivity_coef_name);
+  weak_form_params.SetParam("ThermalExpansionBilinCoefName", thermal_expansion_bilin_coef_name);
+  weak_form_params.SetParam("ThermalExpansionLinCoefName", thermal_expansion_lin_coef_name);
+  weak_form_params.SetParam("ZeroCoefName", zero_coef_name);
+  this->GetProblem()->eq_sys =
+      std::make_unique<hephaestus::ThermalExpansionEquationSystem>(weak_form_params);
+}
+
+
+void ThermalExpansionFormulation::ConstructOperator() {
+  hephaestus::InputParameters &solver_options =
+
+      this->GetProblem()->solver_options;
+  solver_options.SetParam("TempVarName", temp_var_name);
+  solver_options.SetParam("DisplacementVarName", displacement_var_name);
+  solver_options.SetParam("LameCoefName", lame_param_coef_name);
+  solver_options.SetParam("ShearModulusCoefName", shear_modulus_coef_name);
+  solver_options.SetParam("ThermalExpansionCoefName", thermal_expansion_coef_name);
+  solver_options.SetParam("ThermalConductivityCoefName", thermal_conductivity_coef_name);
+  solver_options.SetParam("StressFreeTempCoefName", stress_free_temp_coef_name);
+  solver_options.SetParam("ZeroCoefName", zero_coef_name);
+
+  this->problem->eq_sys_operator = std::make_unique<hephaestus::ThermalExpansionOperator>(
+      *(this->problem->pmesh), this->problem->fespaces,
+      this->problem->gridfunctions, this->problem->bc_map,
+      this->problem->coefficients, this->problem->sources,
+      this->problem->solver_options);
+  this->problem->eq_sys_operator->SetEquationSystem(
+      this->problem->eq_sys.get());
+  this->problem->eq_sys_operator->SetGridFunctions();
+};
+
+void ThermalExpansionFormulation::RegisterGridFunctions() {
+  int &myid = this->GetProblem()->myid_;
+  hephaestus::GridFunctions &gridfunctions = this->GetProblem()->gridfunctions;
+  hephaestus::FESpaces &fespaces = this->GetProblem()->fespaces;
+
+  // Register default ParGridFunctions of state gridfunctions if not provided
+  if (!gridfunctions.Has(temp_var_name)) {
+    if (myid == 0) {
+      MFEM_WARNING(temp_var_name << " not found in gridfunctions: building "
+                                      "gridfunction from defaults");
+    }
+    AddFESpace(std::string("_TempFESpace"), std::string("H1_"));
+    AddGridFunction(temp_var_name, std::string("_TempFESpace"));
+  };
+
+  if (!gridfunctions.Has(displacement_var_name)) {
+      if (myid == 0) {
+        MFEM_WARNING(displacement_var_name << " not found in gridfunctions: building "
+                                        "gridfunction from defaults");
+      }
+      AddFESpace(std::string("_DisplacementFESpace"), std::string("H1_"), 3, mfem::Ordering::byVDIM);
+      AddGridFunction(displacement_var_name, std::string("_DisplacementFESpace"));
+  };
+    // Register time derivatives
+  SteadyStateProblemBuilder::RegisterGridFunctions();
+};
+
+
+
+/**
+ * OPERATOR
+ * 
+ * 
+ * 
+ * 
+ * 
+ * 
+*/
+ThermalExpansionOperator::ThermalExpansionOperator(
+    mfem::ParMesh &pmesh, hephaestus::FESpaces &fespaces,
+    hephaestus::GridFunctions &gridfunctions, hephaestus::BCMap &bc_map,
+    hephaestus::Coefficients &coefficients, hephaestus::Sources &sources,
+    hephaestus::InputParameters &solver_options)    
+    : EquationSystemOperator(pmesh, fespaces, gridfunctions,
+                                            bc_map, coefficients, sources,
+                                            solver_options)
+      {}
+
+
+
+// This method no longer needs to be overridden now that eq is being used
+// void ThermalExpansionOperator::SetGridFunctions() {
+//   state_var_names.push_back(temp_var_name);
+//   state_var_names.push_back(displacement_var_name);
+  
+//   EquationSystemOperator::SetGridFunctions();
+
+//   t_ = new mfem::ParGridFunction(local_test_vars.at(0)->ParFESpace());
+//   u_ = new mfem::ParGridFunction(local_test_vars.at(1)->ParFESpace());
+// };
+
+void ThermalExpansionOperator::Init(mfem::Vector &X) {
+  EquationSystemOperator::Init(X);
+}
+
+void ThermalExpansionOperator::Solve(mfem::Vector &X) {  
+
+  _equation_system->FormLinearSystem(blockA, trueX, trueRhs);
+
+  preconditioner = new mfem::HypreBoomerAMG(*blockA.As<mfem::HypreParMatrix>());
+  preconditioner->SetSystemsOptions(this->pmesh_->Dimension());
+  solver = new mfem::HyprePCG(*blockA.As<mfem::HypreParMatrix>());
+  
+  solver->SetPreconditioner(*preconditioner);
+  solver->SetTol(1e-11);
+  solver->Mult(trueRhs, trueX);
+  _equation_system->RecoverFEMSolution(trueX, _gridfunctions);
+}
+
+/**
+ * EQUATION SYSTEM
+ * 
+ * 
+ * 
+ * 
+ * 
+ * 
+*/
+ThermalExpansionEquationSystem::ThermalExpansionEquationSystem(
+  const hephaestus::InputParameters &params)
+  : EquationSystem(params),
+    temp_var_name(params.GetParam<std::string>("TempVarName")), 
+    displacement_var_name(params.GetParam<std::string>("DisplacementVarName")), 
+    stress_free_temp_coef_name(params.GetParam<std::string>("StressFreeTempCoefName")),
+    thermal_conductivity_coef_name(params.GetParam<std::string>("ThermalConductivityCoefName")),
+    lame_param_coef_name(params.GetParam<std::string>("LameParamCoefName")),
+    shear_modulus_coef_name(params.GetParam<std::string>("ShearModulusCoefName")),
+    thermal_expansion_bilin_coef_name(params.GetParam<std::string>("ThermalExpansionBilinCoefName")),
+    thermal_expansion_lin_coef_name(params.GetParam<std::string>("ThermalExpansionLinCoefName"))
+     {}
+
+void ThermalExpansionEquationSystem::Init(hephaestus::GridFunctions &gridfunctions,
+                                  const hephaestus::FESpaces &fespaces,
+                                  hephaestus::BCMap &bc_map,
+                                  hephaestus::Coefficients &coefficients) {
+
+  EquationSystem::Init(gridfunctions, fespaces, bc_map, coefficients);
+}
+
+void ThermalExpansionEquationSystem::addKernels() {
+  // Add missing variable names
+  addVariableNameIfMissing(temp_var_name);
+  addVariableNameIfMissing(displacement_var_name);
+
+
+  hephaestus::InputParameters diffusionIntegratorParams;
+  diffusionIntegratorParams.SetParam("CoefficientName", thermal_conductivity_coef_name);
+  addKernel(temp_var_name, new hephaestus::DiffusionKernel(diffusionIntegratorParams));
+
+  hephaestus::InputParameters tempLFParams;
+  tempLFParams.SetParam("CoefficientName", zero_coef_name);
+  addKernel(temp_var_name, new hephaestus::DomainLFKernel(tempLFParams));
+
+  hephaestus::InputParameters elasticityIntegratorParams;
+  elasticityIntegratorParams.SetParam("LameParameterCoefName", lame_param_coef_name);
+  elasticityIntegratorParams.SetParam("ShearModulusCoefName", shear_modulus_coef_name);
+  addKernel(displacement_var_name, new hephaestus::LinearElasticityKernel(elasticityIntegratorParams));
+
+  hephaestus::InputParameters domainDivergenceLFParams;
+  domainDivergenceLFParams.SetParam("CoefficientName", thermal_expansion_lin_coef_name);
+  addKernel(displacement_var_name, new hephaestus::DomainDivergenceLFKernel(domainDivergenceLFParams));
+
+  hephaestus::InputParameters mixedWeakDivergenceParams;
+  mixedWeakDivergenceParams.SetParam("CoefficientName", thermal_expansion_bilin_coef_name);
+  addKernel(temp_var_name, displacement_var_name,
+            new hephaestus::MixedWeakDivergenceKernel(mixedWeakDivergenceParams));
+}
+
+} //hephaestus
diff --git a/src/formulations/ThermalExpansion/thermal_expansion_formulation.hpp b/src/formulations/ThermalExpansion/thermal_expansion_formulation.hpp
new file mode 100644
index 000000000..e0eb8d20f
--- /dev/null
+++ b/src/formulations/ThermalExpansion/thermal_expansion_formulation.hpp
@@ -0,0 +1,68 @@
+#pragma once
+#include "../common/pfem_extras.hpp"
+#include "formulation.hpp"
+#include "kernels.hpp"
+#include "inputs.hpp"
+#include "sources.hpp"
+#include "integrators.hpp"
+
+namespace hephaestus {
+
+class ThermalExpansionFormulation : public SteadyStateFormulation {
+public:
+  ThermalExpansionFormulation();
+
+  virtual void ConstructEquationSystem() override;
+
+  virtual void ConstructOperator() override;
+
+  virtual void RegisterGridFunctions() override;
+
+  virtual void RegisterCoefficients() override;
+
+protected:
+  std::string temp_var_name, displacement_var_name, 
+                stress_free_temp_coef_name, lame_param_coef_name, shear_modulus_coef_name, thermal_expansion_coef_name, thermal_conductivity_coef_name,
+                  thermal_expansion_bilin_coef_name, thermal_expansion_lin_coef_name, zero_coef_name;
+};
+
+
+// Equation System
+class ThermalExpansionEquationSystem : public EquationSystem {
+public:
+  ThermalExpansionEquationSystem(const hephaestus::InputParameters &params);
+
+  virtual void Init(hephaestus::GridFunctions &gridfunctions,
+                    const hephaestus::FESpaces &fespaces,
+                    hephaestus::BCMap &bc_map,
+                    hephaestus::Coefficients &coefficients) override;
+  virtual void addKernels() override;
+
+  std::string temp_var_name, displacement_var_name, 
+                stress_free_temp_coef_name, lame_param_coef_name, shear_modulus_coef_name,
+                thermal_conductivity_coef_name, thermal_expansion_bilin_coef_name, thermal_expansion_lin_coef_name, zero_coef_name;
+};
+
+class ThermalExpansionOperator : public EquationSystemOperator {
+public:
+  ThermalExpansionOperator(mfem::ParMesh &pmesh, hephaestus::FESpaces &fespaces,
+                hephaestus::GridFunctions &gridfunctions,
+                hephaestus::BCMap &bc_map,
+                hephaestus::Coefficients &coefficients,
+                hephaestus::Sources &sources,
+                hephaestus::InputParameters &solver_options);
+
+  ~ThermalExpansionOperator(){};
+
+  // virtual void SetGridFunctions() override;
+  virtual void Init(mfem::Vector &X) override;
+  virtual void Solve(mfem::Vector &X) override;
+
+  // Method for manipulating existing coefficients to get the terms needed for the thermal expansion solve
+  void MakeCoefficients();
+
+  mfem::HypreBoomerAMG *preconditioner;
+  mfem::HyprePCG *solver;
+}; //
+
+} // namespace hephaestus
diff --git a/src/formulations/equation_system_operator.cpp b/src/formulations/equation_system_operator.cpp
index 9f5e2ed4f..8a4fb0264 100644
--- a/src/formulations/equation_system_operator.cpp
+++ b/src/formulations/equation_system_operator.cpp
@@ -3,6 +3,7 @@
 namespace hephaestus {
 
 void EquationSystemOperator::SetGridFunctions() {
+  state_var_names = _equation_system->var_names;
   local_test_vars = populateVectorFromNamedFieldsMap<mfem::ParGridFunction>(
       _gridfunctions, state_var_names);
 
@@ -36,6 +37,12 @@ void EquationSystemOperator::SetGridFunctions() {
   }
 };
 
+void EquationSystemOperator::SetEquationSystem(
+  hephaestus::EquationSystem *equation_system) {
+    _equation_system = equation_system;
+  
+}
+
 void EquationSystemOperator::Init(mfem::Vector &X) {
   // Define material property coefficients
   for (unsigned int ind = 0; ind < local_test_vars.size(); ++ind) {
@@ -43,6 +50,7 @@ void EquationSystemOperator::Init(mfem::Vector &X) {
                                      const_cast<mfem::Vector &>(X),
                                      true_offsets[ind]);
   }
+  _equation_system->buildEquationSystem(_bc_map, _sources);
 }
 
 void EquationSystemOperator::Solve(mfem::Vector &X) {}
diff --git a/src/formulations/equation_system_operator.hpp b/src/formulations/equation_system_operator.hpp
index ef64a2075..553b46b23 100644
--- a/src/formulations/equation_system_operator.hpp
+++ b/src/formulations/equation_system_operator.hpp
@@ -20,8 +20,9 @@ class EquationSystemOperator : public mfem::Operator {
         _coefficients(coefficients), _solver_options(solver_options){};
 
   ~EquationSystemOperator(){};
-
+  
   virtual void SetGridFunctions();
+  virtual void SetEquationSystem(hephaestus::EquationSystem *equation_system);
   virtual void Init(mfem::Vector &X);
   virtual void Solve(mfem::Vector &X);
   void Mult(const mfem::Vector &x, mfem::Vector &y) const override{};
@@ -39,6 +40,8 @@ class EquationSystemOperator : public mfem::Operator {
 
   std::vector<mfem::ParGridFunction *> local_test_vars;
 
+  hephaestus::EquationSystem *_equation_system;
+
   int myid_;
   int num_procs_;
   mfem::ParMesh *pmesh_;
diff --git a/src/formulations/formulation.hpp b/src/formulations/formulation.hpp
index 04eed44ed..4f667d10b 100644
--- a/src/formulations/formulation.hpp
+++ b/src/formulations/formulation.hpp
@@ -2,3 +2,4 @@
 #include "frequency_domain_em_formulation.hpp"
 #include "steady_state_em_formulation.hpp"
 #include "time_domain_em_formulation.hpp"
+#include "steady_state_formulation.hpp"
diff --git a/src/formulations/steady_state_formulation.cpp b/src/formulations/steady_state_formulation.cpp
new file mode 100644
index 000000000..f6c6b7a6e
--- /dev/null
+++ b/src/formulations/steady_state_formulation.cpp
@@ -0,0 +1,7 @@
+#include "steady_state_formulation.hpp"
+
+namespace hephaestus {
+
+SteadyStateFormulation::SteadyStateFormulation(){};
+
+} // namespace hephaestus
diff --git a/src/formulations/steady_state_formulation.hpp b/src/formulations/steady_state_formulation.hpp
new file mode 100644
index 000000000..d45319720
--- /dev/null
+++ b/src/formulations/steady_state_formulation.hpp
@@ -0,0 +1,13 @@
+#pragma once
+#include "steady_state_problem_builder.hpp"
+
+namespace hephaestus {
+
+// 
+class SteadyStateFormulation : public hephaestus::SteadyStateProblemBuilder {
+    
+    
+public:
+    SteadyStateFormulation();
+}; 
+} // namespace hephaestus 
\ No newline at end of file
diff --git a/src/integrators/CMakeLists.txt b/src/integrators/CMakeLists.txt
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/integrators/DomainLFH1DivIntegrator.cpp b/src/integrators/DomainLFH1DivIntegrator.cpp
new file mode 100644
index 000000000..aedcab492
--- /dev/null
+++ b/src/integrators/DomainLFH1DivIntegrator.cpp
@@ -0,0 +1,55 @@
+#include "DomainLFH1DivIntegrator.hpp"
+
+
+namespace mfem {
+void DomainLFH1DivIntegrator ::AssembleRHSElementVect(
+   const mfem::FiniteElement &el, mfem::ElementTransformation &Tr, mfem::Vector &elvect)
+{
+    const int dim = el.GetDim();
+    const int dof = el.GetDof();
+    double w, coeff;
+    const int sdim = Tr.GetSpaceDim();
+    
+    dshape.SetSize(dof, dim);
+    gshape.SetSize(dof, dim);
+    divShape.SetSize(dof * dim);
+    
+    dshape = 0.0;
+    gshape = 0.0;
+    divShape = 0.0;
+
+    elvect.SetSize(dof * dim);
+    elvect = 0.0;
+
+    const mfem::IntegrationRule *ir = IntRule;
+    if (ir == NULL)
+    {
+        int intorder = 2 * el.GetOrder();
+        ir = &mfem::IntRules.Get(el.GetGeomType(), intorder);
+    }
+
+    for (int q = 0; q < ir->GetNPoints(); q++)
+    {
+        const mfem::IntegrationPoint &ip = ir->IntPoint(q);
+
+        Tr.SetIntPoint(&ip);
+        el.CalcDShape(ip, dshape);
+        
+        w = ip.weight;
+        
+        coeff = Q->Eval(Tr, ip);  
+        
+        Mult(dshape, Tr.AdjugateJacobian(), gshape);
+       
+        gshape.GradToDiv(divShape);
+        add(elvect, w * coeff, divShape, elvect);
+    }
+}
+} //namespace mfem
+
+// void DomainLFH1DivIntegrator ::AssembleDeltaElementVect(const mfem::FiniteElement &fe,
+//                                         mfem::ElementTransformation &Trans,
+//                                         mfem::Vector &elvect)
+// {
+//    MFEM_ABORT("Not implemented!");
+// }
\ No newline at end of file
diff --git a/src/integrators/DomainLFH1DivIntegrator.hpp b/src/integrators/DomainLFH1DivIntegrator.hpp
new file mode 100644
index 000000000..5933456d5
--- /dev/null
+++ b/src/integrators/DomainLFH1DivIntegrator.hpp
@@ -0,0 +1,43 @@
+#pragma once
+#include "mfem.hpp"
+
+// MFEM integrator to calculate the Rhs term, when your test function is a vector value
+//  represented by multiple copies of a scalar field i.e. H1^3
+// (f, (∇ . v))
+
+
+
+namespace mfem{
+class DomainLFH1DivIntegrator : public mfem::LinearFormIntegrator
+{
+protected:
+    mfem::Coefficient *Q;
+
+private:
+   mfem::Vector shape, divShape;
+   mfem::DenseMatrix dshape, gshape;
+
+public:
+    /// Constructs the domain integrator (Q, grad v)
+    DomainLFH1DivIntegrator (mfem::Coefficient &q)
+        : Q(&q)
+    {}
+
+    ~DomainLFH1DivIntegrator () {};
+
+    // virtual bool SupportsDevice() override { return true; }
+    
+    /** Given a particular Finite Element and a transformation (Tr)
+         computes the element right hand side element vector, elvect. */
+    virtual void AssembleRHSElementVect(const mfem::FiniteElement &el,
+                                        mfem::ElementTransformation &Tr,
+                                        mfem::Vector &elvect) override;
+
+    // virtual void AssembleDeltaElementVect(const mfem::FiniteElement &fe,
+    //                                      mfem::ElementTransformation &Trans,
+    //                                      mfem::Vector &elvect) override;
+
+    // using mfem::LinearFormIntegrator::AssembleRHSElementVect;
+
+};
+} //namespace mfem
\ No newline at end of file
diff --git a/src/integrators/MixedWeakDivergenceIntegrator.cpp b/src/integrators/MixedWeakDivergenceIntegrator.cpp
new file mode 100644
index 000000000..33fe7f868
--- /dev/null
+++ b/src/integrators/MixedWeakDivergenceIntegrator.cpp
@@ -0,0 +1,69 @@
+#include "MixedWeakDivergenceIntegrator.hpp"
+
+
+namespace mfem {
+    
+void MixedWeakDivergenceIntegrator::AssembleElementMatrix2(const mfem::FiniteElement &trial_fe,
+                                                            const mfem::FiniteElement &test_fe,
+                                                            mfem::ElementTransformation &Trans,
+                                                            mfem::DenseMatrix &elmat)
+{
+    //  Get number of dofs/ basis functions
+    dim = trial_fe.GetDim();
+    int trial_dof = trial_fe.GetDof();
+    int test_dof = test_fe.GetDof();
+    double w, coeff;
+    Jadj.SetSize(dim);
+
+    // vector to store values of basis functions
+    trial_shape.SetSize(trial_dof);
+    test_div_shape.SetSize(dim * test_dof);
+
+    // Derivitive placeholder for getting div
+    test_d_shape.SetSize(test_dof, dim);
+    test_g_shape.SetSize(test_dof, dim);
+    
+
+    elmat.SetSize(dim*test_dof, trial_dof);
+    elmat = 0.0;
+    pelmat.SetSize(dim*test_dof, trial_dof);
+    
+    mfem::Geometry::Type geom = Trans.GetGeometryType();
+    
+    const mfem::IntegrationRule *ir = &GetIntRule(trial_fe, test_fe, Trans);
+   
+    for(int i = 0; i < ir->GetNPoints(); i++)
+    {
+        const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+        Trans.SetIntPoint(&ip);
+
+        trial_fe.CalcShape(ip, trial_shape);
+        test_fe.CalcDShape(ip, test_d_shape);
+
+        w = ip.weight;
+
+        if (a)
+        {
+            coeff = a -> Eval(Trans, ip);
+        }
+        Mult(test_d_shape, Trans.AdjugateJacobian(), test_g_shape);
+        trial_shape *= w;
+        trial_shape *= coeff;
+
+        test_g_shape.GradToDiv(test_div_shape);
+        AddMultVWt(test_div_shape, trial_shape, elmat);
+    }
+}
+
+const IntegrationRule& mfem::MixedWeakDivergenceIntegrator::GetIntRule(const mfem::FiniteElement &trial_fe,
+                                   const mfem::FiniteElement &test_fe,
+                                   mfem::ElementTransformation &Trans)
+{
+    int order = trial_fe.GetOrder()
+                + Trans.OrderGrad(&test_fe)
+                + Trans.OrderJ();
+    
+    return mfem::IntRules.Get(trial_fe.GetGeomType(), order);
+} 
+
+} //namespace mfem
\ No newline at end of file
diff --git a/src/integrators/MixedWeakDivergenceIntegrator.hpp b/src/integrators/MixedWeakDivergenceIntegrator.hpp
new file mode 100644
index 000000000..f2023ad0d
--- /dev/null
+++ b/src/integrators/MixedWeakDivergenceIntegrator.hpp
@@ -0,0 +1,48 @@
+#pragma once
+#include "mfem.hpp"
+
+// MFEM mixed integrator to calculate the bilinear form,
+// (u, (∇ . v))
+
+//  For when your test function is a vector value represented by multiple copies of a scalar field i.e. H1^3,
+//  and your trial function is a scalar in H1
+
+
+namespace mfem {
+    
+class MixedWeakDivergenceIntegrator : public mfem::BilinearFormIntegrator
+{
+protected:
+    mfem::Coefficient *a;
+    
+private:
+    mfem::Vector trial_shape;
+    mfem::Vector test_div_shape;
+    mfem::DenseMatrix test_d_shape;
+    mfem::DenseMatrix pelmat;
+    mfem::DenseMatrix test_g_shape;
+    mfem::DenseMatrix Jadj;
+
+    const mfem::DofToQuad *trial_maps, *test_maps; ///< Not owned
+    const mfem::GeometricFactors *geom;            ///< Not owned
+    int dim;
+
+public:
+
+    MixedWeakDivergenceIntegrator(mfem::Coefficient &a) :
+        a(&a), trial_maps(NULL), test_maps(NULL), geom(NULL) 
+        {}        
+
+    ~MixedWeakDivergenceIntegrator() {;}
+
+    void AssembleElementMatrix2(const mfem::FiniteElement &trial_fe,
+                                const mfem::FiniteElement &test_fe,
+                                mfem::ElementTransformation &Trans,
+                                mfem::DenseMatrix &elmat) override;
+
+    const mfem::IntegrationRule& GetIntRule(const mfem::FiniteElement &trial_fe,
+                                    const mfem::FiniteElement &test_fe,
+                                    mfem::ElementTransformation &Trans);
+};
+
+}  // namespace mfem
\ No newline at end of file
diff --git a/src/integrators/integrators.hpp b/src/integrators/integrators.hpp
new file mode 100644
index 000000000..704ab7e48
--- /dev/null
+++ b/src/integrators/integrators.hpp
@@ -0,0 +1,3 @@
+#pragma once
+#include "MixedWeakDivergenceIntegrator.hpp"
+#include "DomainLFH1DivIntegrator.hpp"
\ No newline at end of file
diff --git a/src/kernels/domain_divergence_lf_kernel.cpp b/src/kernels/domain_divergence_lf_kernel.cpp
new file mode 100644
index 000000000..11c9b041b
--- /dev/null
+++ b/src/kernels/domain_divergence_lf_kernel.cpp
@@ -0,0 +1,21 @@
+#include "domain_divergence_lf_kernel.hpp"
+
+namespace hephaestus {
+
+DomainDivergenceLFKernel::DomainDivergenceLFKernel(const hephaestus::InputParameters &params)
+    : Kernel(params),
+      coef_name(params.GetParam<std::string>("CoefficientName")) {}
+
+void DomainDivergenceLFKernel::Init(hephaestus::GridFunctions &gridfunctions,
+                          const hephaestus::FESpaces &fespaces,
+                          hephaestus::BCMap &bc_map,
+                          hephaestus::Coefficients &coefficients) {
+
+  coef = coefficients.scalars.Get(coef_name);
+}
+
+void DomainDivergenceLFKernel::Apply(mfem::ParLinearForm *lf) {
+  lf->AddDomainIntegrator(new mfem::DomainLFH1DivIntegrator(*coef));
+}
+
+} // namespace hephaestus
diff --git a/src/kernels/domain_divergence_lf_kernel.hpp b/src/kernels/domain_divergence_lf_kernel.hpp
new file mode 100644
index 000000000..cdb49ed7a
--- /dev/null
+++ b/src/kernels/domain_divergence_lf_kernel.hpp
@@ -0,0 +1,25 @@
+#pragma once
+#include "kernel_base.hpp"
+#include "DomainLFH1DivIntegrator.hpp"
+
+namespace hephaestus {
+
+/*
+(∇∙v, f)
+*/
+class DomainDivergenceLFKernel : public Kernel<mfem::ParLinearForm> {
+public:
+  DomainDivergenceLFKernel(const hephaestus::InputParameters &params);
+  ~DomainDivergenceLFKernel();
+  virtual void Init(hephaestus::GridFunctions &gridfunctions,
+                    const hephaestus::FESpaces &fespaces,
+                    hephaestus::BCMap &bc_map,
+                    hephaestus::Coefficients &coefficients) override;
+  virtual void Apply(mfem::ParLinearForm *lf) override;
+
+  std::string coef_name;
+  mfem::Coefficient *coef;
+  
+};
+
+}; // namespace hephaestus
diff --git a/src/kernels/domain_lf_kernel.cpp b/src/kernels/domain_lf_kernel.cpp
new file mode 100644
index 000000000..c63218d9f
--- /dev/null
+++ b/src/kernels/domain_lf_kernel.cpp
@@ -0,0 +1,21 @@
+#include "domain_lf_kernel.hpp"
+
+namespace hephaestus {
+
+DomainLFKernel::DomainLFKernel(const hephaestus::InputParameters &params)
+    : Kernel(params),
+      coef_name(params.GetParam<std::string>("CoefficientName")) {}
+
+void DomainLFKernel::Init(hephaestus::GridFunctions &gridfunctions,
+                           const hephaestus::FESpaces &fespaces,
+                           hephaestus::BCMap &bc_map,
+                           hephaestus::Coefficients &coefficients) {
+
+  coef = coefficients.scalars.Get(coef_name);
+}
+
+void DomainLFKernel::Apply(mfem::ParLinearForm *lf) {
+  lf->AddDomainIntegrator(new mfem::DomainLFIntegrator(*coef));
+};
+
+} // namespace hephaestus
diff --git a/src/kernels/domain_lf_kernel.hpp b/src/kernels/domain_lf_kernel.hpp
new file mode 100644
index 000000000..2d498b6bd
--- /dev/null
+++ b/src/kernels/domain_lf_kernel.hpp
@@ -0,0 +1,21 @@
+#pragma once
+#include "kernel_base.hpp"
+
+namespace hephaestus {
+
+/*
+(σ ∇ V, ∇ V')
+*/
+class DomainLFKernel : public Kernel<mfem::ParLinearForm> {
+public:
+  DomainLFKernel(const hephaestus::InputParameters &params);
+  virtual void Init(hephaestus::GridFunctions &gridfunctions,
+                    const hephaestus::FESpaces &fespaces,
+                    hephaestus::BCMap &bc_map,
+                    hephaestus::Coefficients &coefficients) override;
+  virtual void Apply(mfem::ParLinearForm *lf) override;
+  std::string coef_name;
+  mfem::Coefficient *coef;
+};
+
+}; // namespace hephaestus
diff --git a/src/kernels/kernels.hpp b/src/kernels/kernels.hpp
index 61d683dce..1ba8565bb 100644
--- a/src/kernels/kernels.hpp
+++ b/src/kernels/kernels.hpp
@@ -6,3 +6,9 @@
 #include "vector_fe_weak_divergence_kernel.hpp"
 #include "weak_curl_curl_kernel.hpp"
 #include "weak_curl_kernel.hpp"
+#include "mixed_weak_divergence_kernel.hpp"
+#include "domain_divergence_lf_kernel.hpp"
+#include "domain_lf_kernel.hpp"
+#include "linear_elasticity_kernel.hpp"
+
+
diff --git a/src/kernels/linear_elasticity_kernel.cpp b/src/kernels/linear_elasticity_kernel.cpp
new file mode 100644
index 000000000..eab170134
--- /dev/null
+++ b/src/kernels/linear_elasticity_kernel.cpp
@@ -0,0 +1,25 @@
+#include "linear_elasticity_kernel.hpp"
+
+namespace hephaestus {
+
+LinearElasticityKernel::LinearElasticityKernel(
+    const hephaestus::InputParameters &params)
+    : Kernel(params),
+      lame_parameter_name(params.GetParam<std::string>("LameParameterCoefName")),
+      shear_modulus_name(params.GetParam<std::string>("ShearModulusCoefName"))
+    {}
+
+void LinearElasticityKernel::Init(hephaestus::GridFunctions &gridfunctions,
+                                     const hephaestus::FESpaces &fespaces,
+                                     hephaestus::BCMap &bc_map,
+                                     hephaestus::Coefficients &coefficients) {
+
+  lame_parameter = coefficients.scalars.Get(lame_parameter_name);
+  shear_modulus = coefficients.scalars.Get(shear_modulus_name);
+}
+
+void LinearElasticityKernel::Apply(mfem::ParBilinearForm *blf) {
+  blf->AddDomainIntegrator(new mfem::ElasticityIntegrator(*lame_parameter, *shear_modulus));
+}
+
+} // namespace hephaestus
diff --git a/src/kernels/linear_elasticity_kernel.hpp b/src/kernels/linear_elasticity_kernel.hpp
new file mode 100644
index 000000000..ff1aaa519
--- /dev/null
+++ b/src/kernels/linear_elasticity_kernel.hpp
@@ -0,0 +1,24 @@
+#pragma once
+#include "kernel_base.hpp"
+
+namespace hephaestus {
+
+/*
+Implements the provided mfem linear elasticity integrator into a hephaestus kernel.
+a(u,v) = (λ div(u), div(v)) + (2 μ e(u), e(v)),
+where e(v) = (1/2) (grad(v) + grad(v)^T).
+*/
+class LinearElasticityKernel : public Kernel<mfem::ParBilinearForm> {
+public:
+  LinearElasticityKernel(const hephaestus::InputParameters &params);
+  virtual void Init(hephaestus::GridFunctions &gridfunctions,
+                    const hephaestus::FESpaces &fespaces,
+                    hephaestus::BCMap &bc_map,
+                    hephaestus::Coefficients &coefficients) override;
+  virtual void Apply(mfem::ParBilinearForm *blf) override;
+  std::string lame_parameter_name, shear_modulus_name;
+  mfem::Coefficient *lame_parameter, *shear_modulus;
+
+};
+
+}; // namespace hephaestus
diff --git a/src/kernels/mixed_weak_divergence_kernel.cpp b/src/kernels/mixed_weak_divergence_kernel.cpp
new file mode 100644
index 000000000..759cd87bd
--- /dev/null
+++ b/src/kernels/mixed_weak_divergence_kernel.cpp
@@ -0,0 +1,22 @@
+#include "mixed_weak_divergence_kernel.hpp"
+
+namespace hephaestus {
+
+MixedWeakDivergenceKernel::MixedWeakDivergenceKernel(
+    const hephaestus::InputParameters &params)
+    : Kernel(params),
+      coef_name(params.GetParam<std::string>("CoefficientName")) {}
+
+void MixedWeakDivergenceKernel::Init(hephaestus::GridFunctions &gridfunctions,
+                                     const hephaestus::FESpaces &fespaces,
+                                     hephaestus::BCMap &bc_map,
+                                     hephaestus::Coefficients &coefficients) {
+
+  coef = coefficients.scalars.Get(coef_name);
+}
+
+void MixedWeakDivergenceKernel::Apply(mfem::ParMixedBilinearForm *mblf) {
+  mblf->AddDomainIntegrator(new mfem::MixedWeakDivergenceIntegrator(*coef));
+}
+
+} // namespace hephaestus
diff --git a/src/kernels/mixed_weak_divergence_kernel.hpp b/src/kernels/mixed_weak_divergence_kernel.hpp
new file mode 100644
index 000000000..61bf6d5f2
--- /dev/null
+++ b/src/kernels/mixed_weak_divergence_kernel.hpp
@@ -0,0 +1,22 @@
+#pragma once
+#include "kernel_base.hpp"
+#include "MixedWeakDivergenceIntegrator.hpp"
+
+namespace hephaestus {
+
+/*
+(σ ∇ V, u')
+*/
+class MixedWeakDivergenceKernel : public Kernel<mfem::ParMixedBilinearForm> {
+public:
+  MixedWeakDivergenceKernel(const hephaestus::InputParameters &params);
+  virtual void Init(hephaestus::GridFunctions &gridfunctions,
+                    const hephaestus::FESpaces &fespaces,
+                    hephaestus::BCMap &bc_map,
+                    hephaestus::Coefficients &coefficients) override;
+  virtual void Apply(mfem::ParMixedBilinearForm *mblf) override;
+  std::string coef_name;
+  mfem::Coefficient *coef;
+};
+
+}; // namespace hephaestus
diff --git a/src/problem_builders/steady_state_problem_builder.cpp b/src/problem_builders/steady_state_problem_builder.cpp
index c74d91558..98382e06b 100644
--- a/src/problem_builders/steady_state_problem_builder.cpp
+++ b/src/problem_builders/steady_state_problem_builder.cpp
@@ -2,13 +2,6 @@
 
 namespace hephaestus {
 
-void SteadyStateProblemBuilder::InitializeKernels() {
-  this->problem->preprocessors.Init(this->problem->gridfunctions,
-                                    this->problem->coefficients);
-  this->problem->sources.Init(this->problem->gridfunctions,
-                              this->problem->fespaces, this->problem->bc_map,
-                              this->problem->coefficients);
-}
 void SteadyStateProblemBuilder::ConstructOperator() {
   this->problem->eq_sys_operator =
       std::make_unique<hephaestus::EquationSystemOperator>(
@@ -19,10 +12,34 @@ void SteadyStateProblemBuilder::ConstructOperator() {
   this->problem->eq_sys_operator->SetGridFunctions();
 }
 
+void SteadyStateProblemBuilder::InitializeKernels() {
+  this->problem->eq_sys->Init(
+      this->problem->gridfunctions, this->problem->fespaces,
+      this->problem->bc_map, this->problem->coefficients);
+  this->problem->preprocessors.Init(this->problem->gridfunctions,
+                                    this->problem->coefficients);
+  this->problem->sources.Init(this->problem->gridfunctions,
+                              this->problem->fespaces, this->problem->bc_map,
+                              this->problem->coefficients);
+}
+
 void SteadyStateProblemBuilder::ConstructState() {
   this->problem->F = new mfem::BlockVector(
       this->problem->eq_sys_operator->true_offsets); // Vector of dofs
   this->problem->eq_sys_operator->Init(
       *(this->problem->F)); // Set up initial conditions
 }
+
+void SteadyStateProblemBuilder::ConstructEquationSystem() {
+  hephaestus::InputParameters params;
+  this->problem->eq_sys = 
+      std::make_unique<hephaestus::EquationSystem>(params);
+}
+
+void SteadyStateProblemBuilder::RegisterGridFunctions() {
+  std::vector<std::string> gridfunction_names;
+  for (auto const &[name, gf] : this->problem->gridfunctions) {
+    gridfunction_names.push_back(name);
+  }
+}
 } // namespace hephaestus
diff --git a/src/problem_builders/steady_state_problem_builder.hpp b/src/problem_builders/steady_state_problem_builder.hpp
index d2dccefde..a0d559a6c 100644
--- a/src/problem_builders/steady_state_problem_builder.hpp
+++ b/src/problem_builders/steady_state_problem_builder.hpp
@@ -38,15 +38,15 @@ class SteadyStateProblemBuilder : public hephaestus::ProblemBuilder {
 
   virtual void RegisterFESpaces() override{};
 
-  virtual void RegisterGridFunctions() override{};
+  virtual void RegisterGridFunctions() override;
 
   virtual void RegisterAuxSolvers() override{};
 
   virtual void RegisterCoefficients() override{};
 
-  virtual void InitializeKernels() override;
+  virtual void ConstructEquationSystem() override;
 
-  virtual void ConstructEquationSystem() override{};
+  virtual void InitializeKernels() override;
 
   virtual void ConstructOperator() override;