diff --git a/lib/Elements b/lib/Elements
index cf521e627..0de655cd2 160000
--- a/lib/Elements
+++ b/lib/Elements
@@ -1 +1 @@
-Subproject commit cf521e627c76926425f346a89f41d61087631468
+Subproject commit 0de655cd260e60e4f037a78561dddfb9c8757fc7
diff --git a/single-node-refactor/CMakeLists.txt b/single-node-refactor/CMakeLists.txt
index 78d9f53e7..817595fe6 100755
--- a/single-node-refactor/CMakeLists.txt
+++ b/single-node-refactor/CMakeLists.txt
@@ -33,9 +33,20 @@ if(NOT CMAKE_BUILD_TYPE)
  # set(CMAKE_BUILD_TYPE Release)
 endif(NOT CMAKE_BUILD_TYPE)
 
-include_directories(../lib/Elements/matar)
-include_directories(../lib/Elements/matar/src/include)
-add_subdirectory(../lib/Elements/matar cbin)
+###DANinclude_directories(../lib/Elements/matar)
+###DANinclude_directories(../lib/Elements/matar/src/include)
+
+###DANadd_subdirectory(../lib/Elements/matar cbin)
+if (MPI)
+    find_package(MPI REQUIRED)
+    add_definitions(-DHAVE_MPI=1)
+    list(APPEND LINKING_LIBRARIES MPI::MPI_CXX)
+    if (GPU_AWARE_MPI)
+        add_definitions(-DHAVE_GPU_AWARE_MPI)
+#    elseif (GPU_SHARED_MEM)
+#        add_definitions(-DHAVE_GPU_SHARED_MEM)
+    endif()
+endif()
 
 # include_directories(Mesh-Builder)
 # add_subdirectory(Mesh-Builder)
diff --git a/single-node-refactor/scripts/cmake_build.sh b/single-node-refactor/scripts/cmake_build.sh
index d0c1e8442..9d4b46ef3 100755
--- a/single-node-refactor/scripts/cmake_build.sh
+++ b/single-node-refactor/scripts/cmake_build.sh
@@ -19,6 +19,9 @@ if [ "$solver" = "SGH" ]; then
 else
     echo "Error: Solver not supported."
 fi
+cmake_options+=(
+    -D MPI=ON
+)
 
 # Print CMake options for reference
 echo "CMake Options: ${cmake_options[@]}"
diff --git a/single-node-refactor/scripts/matar-install.sh b/single-node-refactor/scripts/matar-install.sh
index 6ec120738..8b5c6ba78 100755
--- a/single-node-refactor/scripts/matar-install.sh
+++ b/single-node-refactor/scripts/matar-install.sh
@@ -19,6 +19,9 @@ else
         -D Matar_ENABLE_KOKKOS=ON
     )
 fi
+cmake_options+=(
+    -D Matar_ENABLE_MPI=ON
+)
 
 # Print CMake options for reference
 echo "CMake Options: ${cmake_options[@]}"
diff --git a/single-node-refactor/src/CMakeLists.txt b/single-node-refactor/src/CMakeLists.txt
index 859a059e9..06bdf2ff2 100755
--- a/single-node-refactor/src/CMakeLists.txt
+++ b/single-node-refactor/src/CMakeLists.txt
@@ -37,6 +37,10 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${VECTOR_CXX_FLAGS}")
 
 find_package(Kokkos REQUIRED)
 find_package(Matar REQUIRED)
+if (MPI)
+    find_package(MPI REQUIRED)
+    add_definitions(-DHAVE_MPI=1)
+endif()
 
 
 # if(BUILD_KOKKOS_SGH)
@@ -77,3 +81,6 @@ include_directories(common)
 
 add_executable(Fierro main.cpp solver.cpp )
 target_link_libraries(Fierro PRIVATE matar parse_yaml sgh_solver Kokkos::kokkos)
+if (MPI)
+    target_link_libraries(Fierro PRIVATE MPI::MPI_CXX)
+endif()
diff --git a/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h b/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h
index 4c93c4ff7..7cc7fb228 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h
+++ b/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h
@@ -136,13 +136,13 @@ class SGH : public Solver
     void boundary_velocity(
         const mesh_t& mesh,
         const BoundaryCondition_t& Boundary,
-        DCArrayKokkos<double>& node_vel,
+        MPIArrayKokkos<double>& node_vel,
         const double time_value) const;
 
     void boundary_contact(
         const mesh_t& mesh,
         const BoundaryCondition_t& Boundary,
-        DCArrayKokkos<double>& node_vel,
+        MPIArrayKokkos<double>& node_vel,
         const double time_value) const;
 
     // **** Functions defined in energy_sgh.cpp **** //
@@ -150,7 +150,7 @@ class SGH : public Solver
         double rk_alpha,
         double dt,
         const mesh_t& mesh,
-        const DCArrayKokkos<double>& node_vel,
+        const MPIArrayKokkos<double>& node_vel,
         const DCArrayKokkos<double>& node_coords,
         DCArrayKokkos<double>& MaterialPoints_sie,
         const DCArrayKokkos<double>& MaterialPoints_mass,
@@ -161,7 +161,7 @@ class SGH : public Solver
         const Material_t& Materials,
         const mesh_t& mesh,
         const DCArrayKokkos<double>& node_coords,
-        const DCArrayKokkos<double>& node_vel,
+        const MPIArrayKokkos<double>& node_vel,
         const DCArrayKokkos<double>& MaterialPoints_den,
         const DCArrayKokkos<double>& MaterialPoints_sie,
         const DCArrayKokkos<double>& MaterialPoints_pres,
@@ -182,7 +182,7 @@ class SGH : public Solver
         const Material_t& Materials,
         const mesh_t& mesh,
         const DCArrayKokkos<double>& node_coords,
-        const DCArrayKokkos<double>& node_vel,
+        const MPIArrayKokkos<double>& node_vel,
         const DCArrayKokkos<double>& MaterialPoints_den,
         const DCArrayKokkos<double>& MaterialPoints_sie,
         const DCArrayKokkos<double>& MaterialPoints_pres,
@@ -205,14 +205,14 @@ class SGH : public Solver
         const size_t num_dims,
         const size_t num_nodes,
         DCArrayKokkos<double>& node_coords,
-        const DCArrayKokkos<double>& node_vel) const;
+        const MPIArrayKokkos<double>& node_vel) const;
 
     // **** Functions defined in momentum.cpp **** //
     void update_velocity(
         double rk_alpha,
         double dt,
         const mesh_t& mesh,
-        DCArrayKokkos<double>& node_vel,
+        MPIArrayKokkos<double>& node_vel,
         const DCArrayKokkos<double>& node_mass,
         const DCArrayKokkos<double>& corner_force) const;
 
@@ -220,7 +220,7 @@ class SGH : public Solver
     void get_velgrad(
         ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const ViewCArrayKokkos<double>& b_matrix,
         const double GaussPoints_vol,
         const size_t elem_gid) const;
@@ -229,7 +229,7 @@ class SGH : public Solver
     void get_velgrad2D(
         ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const ViewCArrayKokkos<double>& b_matrix,
         const double GaussPoints_vol,
         const double elem_area,
@@ -239,14 +239,14 @@ class SGH : public Solver
         DCArrayKokkos<double>& GaussPoints_div,
         const mesh_t mesh,
         const DCArrayKokkos<double>& node_coords,
-        const DCArrayKokkos<double>& node_vel,
+        const MPIArrayKokkos<double>& node_vel,
         const DCArrayKokkos<double>& GaussPoints_vol) const;
 
     void get_divergence2D(
         DCArrayKokkos<double>& GaussPoints_div,
         const mesh_t mesh,
         const DCArrayKokkos<double>& node_coords,
-        const DCArrayKokkos<double>& node_vel,
+        const MPIArrayKokkos<double>& node_vel,
         const DCArrayKokkos<double>& GaussPoints_vol) const;
 
     KOKKOS_FUNCTION
@@ -260,7 +260,7 @@ class SGH : public Solver
         const Material_t& Materials,
         const mesh_t& mesh,
         const DCArrayKokkos<double>& node_coords,
-        const DCArrayKokkos<double>& node_vel,
+        const MPIArrayKokkos<double>& node_vel,
         DCArrayKokkos<double>& MaterialPoints_den,
         DCArrayKokkos<double>& MaterialPoints_pres,
         DCArrayKokkos<double>& MaterialPoints_stress,
@@ -278,7 +278,7 @@ class SGH : public Solver
         const Material_t& Materials,
         const mesh_t& mesh,
         const DCArrayKokkos<double>& node_coords,
-        const DCArrayKokkos<double>& node_vel,
+        const MPIArrayKokkos<double>& node_vel,
         DCArrayKokkos<double>& MaterialPoints_den,
         DCArrayKokkos<double>& MaterialPoints_pres,
         DCArrayKokkos<double>& MaterialPoints_stress,
@@ -295,7 +295,7 @@ class SGH : public Solver
     // NOTE: Consider pulling up
     void rk_init(
         DCArrayKokkos<double>& node_coords,
-        DCArrayKokkos<double>& node_vel,
+        MPIArrayKokkos<double>& node_vel,
         DCArrayKokkos<double>& MaterialPoints_sie,
         DCArrayKokkos<double>& MaterialPoints_stress,
         const size_t num_dims,
@@ -305,7 +305,7 @@ class SGH : public Solver
     void get_timestep(
         mesh_t& mesh,
         DCArrayKokkos<double>& node_coords,
-        DCArrayKokkos<double>& node_vel,
+        MPIArrayKokkos<double>& node_vel,
         DCArrayKokkos<double>& MaterialPoints_sspd,
         DCArrayKokkos<double>& GaussPoints_vol,
         double time_value,
@@ -320,7 +320,7 @@ class SGH : public Solver
     void get_timestep2D(
         mesh_t& mesh,
         DCArrayKokkos<double>& node_coords,
-        DCArrayKokkos<double>& node_vel,
+        MPIArrayKokkos<double>& node_vel,
         DCArrayKokkos<double>& MaterialPoints_sspd,
         DCArrayKokkos<double>& GaussPoints_vol,
         double time_value,
@@ -358,7 +358,7 @@ class SGH : public Solver
         const ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
         const DCArrayKokkos<double>&    node_coords,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const double vol,
         const double dt,
         const double rk_alpha);
diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp
index 4253b4a1f..b23a45a40 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp
@@ -48,7 +48,7 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 /////////////////////////////////////////////////////////////////////////////
 void SGH::boundary_velocity(const mesh_t&      mesh,
                             const BoundaryCondition_t& BoundaryConditions,
-                            DCArrayKokkos<double>& node_vel,
+                            MPIArrayKokkos<double>& node_vel,
                             const double time_value) const
 {
 
@@ -94,7 +94,7 @@ void SGH::boundary_velocity(const mesh_t&      mesh,
 /////////////////////////////////////////////////////////////////////////////
 void SGH::boundary_contact(const mesh_t& mesh,
                            const BoundaryCondition_t& BoundaryConditions,
-                           DCArrayKokkos<double>& node_vel,
+                           MPIArrayKokkos<double>& node_vel,
                            const double time_value) const
 {
     return;
diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp
index bb28f30e6..62d373e60 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp
@@ -53,7 +53,7 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 void SGH::update_energy(double rk_alpha,
     double dt,
     const mesh_t& mesh,
-    const DCArrayKokkos<double>& node_vel,
+    const MPIArrayKokkos<double>& node_vel,
     const DCArrayKokkos<double>& node_coords,
     DCArrayKokkos<double>& MaterialPoints_sie,
     const DCArrayKokkos<double>& MaterialPoints_mass,
diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp
index 6faed136e..c03fe7d2c 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp
@@ -62,7 +62,7 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 void SGH::get_force(const Material_t& Materials,
     const mesh_t& mesh,
     const DCArrayKokkos<double>& node_coords,
-    const DCArrayKokkos<double>& node_vel,
+    const MPIArrayKokkos<double>& node_vel,
     const DCArrayKokkos<double>& MaterialPoints_den,
     const DCArrayKokkos<double>& MaterialPoints_sie,
     const DCArrayKokkos<double>& MaterialPoints_pres,
@@ -436,7 +436,7 @@ void SGH::get_force(const Material_t& Materials,
 void SGH::get_force_2D(const Material_t& Materials,
     const mesh_t& mesh,
     const DCArrayKokkos<double>& node_coords,
-    const DCArrayKokkos<double>& node_vel,
+    const MPIArrayKokkos<double>& node_vel,
     const DCArrayKokkos<double>& MaterialPoints_den,
     const DCArrayKokkos<double>& MaterialPoints_sie,
     const DCArrayKokkos<double>& MaterialPoints_pres,
diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp
index 86f2b6800..265a79a73 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp
@@ -50,7 +50,7 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 void SGH::update_velocity(double rk_alpha,
     double dt,
     const mesh_t& mesh,
-    DCArrayKokkos<double>& node_vel,
+    MPIArrayKokkos<double>& node_vel,
     const DCArrayKokkos<double>& node_mass,
     const DCArrayKokkos<double>& corner_force) const
 {
@@ -101,7 +101,7 @@ void SGH::update_velocity(double rk_alpha,
 KOKKOS_FUNCTION
 void SGH::get_velgrad(ViewCArrayKokkos<double>& vel_grad,
     const ViewCArrayKokkos<size_t>& elem_node_gids,
-    const DCArrayKokkos<double>&    node_vel,
+    const MPIArrayKokkos<double>&    node_vel,
     const ViewCArrayKokkos<double>& b_matrix,
     const double elem_vol,
     const size_t elem_gid) const
@@ -196,7 +196,7 @@ void SGH::get_velgrad(ViewCArrayKokkos<double>& vel_grad,
 KOKKOS_FUNCTION
 void SGH::get_velgrad2D(ViewCArrayKokkos<double>& vel_grad,
     const ViewCArrayKokkos<size_t>& elem_node_gids,
-    const DCArrayKokkos<double>&    node_vel,
+    const MPIArrayKokkos<double>&    node_vel,
     const ViewCArrayKokkos<double>& b_matrix,
     const double elem_vol,
     const double elem_area,
@@ -266,7 +266,7 @@ void SGH::get_velgrad2D(ViewCArrayKokkos<double>& vel_grad,
 void SGH::get_divergence(DCArrayKokkos<double>& elem_div,
     const mesh_t mesh,
     const DCArrayKokkos<double>& node_coords,
-    const DCArrayKokkos<double>& node_vel,
+    const MPIArrayKokkos<double>& node_vel,
     const DCArrayKokkos<double>& elem_vol) const
 {
     // --- calculate the forces acting on the nodes from the element ---
@@ -342,7 +342,7 @@ void SGH::get_divergence(DCArrayKokkos<double>& elem_div,
 void SGH::get_divergence2D(DCArrayKokkos<double>& elem_div,
     const mesh_t mesh,
     const DCArrayKokkos<double>& node_coords,
-    const DCArrayKokkos<double>& node_vel,
+    const MPIArrayKokkos<double>& node_vel,
     const DCArrayKokkos<double>& elem_vol) const
 {
     // --- calculate the forces acting on the nodes from the element ---
diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/position.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/position.cpp
index 2dd0d59f9..1013777a6 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/src/position.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver/src/position.cpp
@@ -52,7 +52,7 @@ void SGH::update_position(double rk_alpha,
     const size_t num_dims,
     const size_t num_nodes,
     DCArrayKokkos<double>& node_coords,
-    const DCArrayKokkos<double>& node_vel) const
+    const MPIArrayKokkos<double>& node_vel) const
 {
     // loop over all the nodes in the mesh
     FOR_ALL(node_gid, 0, num_nodes, {
@@ -61,4 +61,4 @@ void SGH::update_position(double rk_alpha,
             node_coords(1, node_gid, dim) = node_coords(0, node_gid, dim) + rk_alpha * dt * half_vel;
         }
     }); // end parallel for over nodes
-} // end subroutine
\ No newline at end of file
+} // end subroutine
diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp
index 1a748ab08..b27baa8f3 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp
@@ -60,7 +60,7 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 void SGH::update_state(const Material_t& Materials,
     const mesh_t& mesh,
     const DCArrayKokkos<double>& node_coords,
-    const DCArrayKokkos<double>& node_vel,
+    const MPIArrayKokkos<double>& node_vel,
     DCArrayKokkos<double>& MaterialPoints_den,
     DCArrayKokkos<double>& MaterialPoints_pres,
     DCArrayKokkos<double>& MaterialPoints_stress,
@@ -220,7 +220,7 @@ void SGH::update_state(const Material_t& Materials,
 void SGH::update_state2D(const Material_t& Materials,
     const mesh_t& mesh,
     const DCArrayKokkos<double>& node_coords,
-    const DCArrayKokkos<double>& node_vel,
+    const MPIArrayKokkos<double>& node_vel,
     DCArrayKokkos<double>& MaterialPoints_den,
     DCArrayKokkos<double>& MaterialPoints_pres,
     DCArrayKokkos<double>& MaterialPoints_stress,
diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp
index 1cfe57695..df3ce286e 100644
--- a/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp
+++ b/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp
@@ -50,7 +50,7 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ///
 /////////////////////////////////////////////////////////////////////////////
 void SGH::rk_init(DCArrayKokkos<double>& node_coords,
-    DCArrayKokkos<double>& node_vel,
+    MPIArrayKokkos<double>& node_vel,
     DCArrayKokkos<double>& MaterialPoints_sie,
     DCArrayKokkos<double>& MaterialPoints_stress,
     const size_t num_dims,
@@ -101,7 +101,7 @@ void SGH::rk_init(DCArrayKokkos<double>& node_coords,
 /////////////////////////////////////////////////////////////////////////////
 void SGH::get_timestep(mesh_t& mesh,
     DCArrayKokkos<double>&     node_coords,
-    DCArrayKokkos<double>&     node_vel,
+    MPIArrayKokkos<double>&     node_vel,
     DCArrayKokkos<double>&     MaterialPoints_sspd,
     DCArrayKokkos<double>&     GaussPoints_vol,
     double time_value,
@@ -211,7 +211,7 @@ void SGH::get_timestep(mesh_t& mesh,
 /////////////////////////////////////////////////////////////////////////////
 void SGH::get_timestep2D(mesh_t& mesh,
     DCArrayKokkos<double>& node_coords,
-    DCArrayKokkos<double>& node_vel,
+    MPIArrayKokkos<double>& node_vel,
     DCArrayKokkos<double>& MaterialPoints_sspd,
     DCArrayKokkos<double>& GaussPoints_vol,
     double time_value,
diff --git a/single-node-refactor/src/boundary_conditions/velocity/constant_velocity_bc.h b/single-node-refactor/src/boundary_conditions/velocity/constant_velocity_bc.h
index e729022b8..cd09dd7fb 100644
--- a/single-node-refactor/src/boundary_conditions/velocity/constant_velocity_bc.h
+++ b/single-node-refactor/src/boundary_conditions/velocity/constant_velocity_bc.h
@@ -65,7 +65,7 @@ namespace ConstantVelocityBC {
                          const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                          const DCArrayKokkos<double>& bc_global_vars,
                          const DCArrayKokkos<double>& bc_state_vars,
-                         const DCArrayKokkos<double>& node_vel,
+                         const MPIArrayKokkos<double>& node_vel,
                          const double time_value,
                          const size_t bdy_node_gid,
                          const size_t bdy_set)
@@ -82,4 +82,4 @@ namespace ConstantVelocityBC {
 
 } // end namespace
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/boundary_conditions/velocity/no_velocity_bc.h b/single-node-refactor/src/boundary_conditions/velocity/no_velocity_bc.h
index 63b92ef6e..59da6c25a 100644
--- a/single-node-refactor/src/boundary_conditions/velocity/no_velocity_bc.h
+++ b/single-node-refactor/src/boundary_conditions/velocity/no_velocity_bc.h
@@ -63,7 +63,7 @@ namespace NoVelocityBC {
                          const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                          const DCArrayKokkos<double>& bc_global_vars,
                          const DCArrayKokkos<double>& bc_state_vars,
-                         const DCArrayKokkos<double>& node_vel,
+                         const MPIArrayKokkos<double>& node_vel,
                          const double time_value,
                          const size_t bdy_node_gid,
                          const size_t bdy_set)
@@ -75,4 +75,4 @@ namespace NoVelocityBC {
 
 } // end namespace
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/boundary_conditions/velocity/piston_velocity_bc.h b/single-node-refactor/src/boundary_conditions/velocity/piston_velocity_bc.h
index 35d374e5e..2dcc958c7 100644
--- a/single-node-refactor/src/boundary_conditions/velocity/piston_velocity_bc.h
+++ b/single-node-refactor/src/boundary_conditions/velocity/piston_velocity_bc.h
@@ -67,7 +67,7 @@ namespace PistonVelocityBC {
                          const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                          const DCArrayKokkos<double>& bc_global_vars,
                          const DCArrayKokkos<double>& bc_state_vars,
-                         const DCArrayKokkos<double>& node_vel,
+                         const MPIArrayKokkos<double>& node_vel,
                          const double time_value,
                          const size_t bdy_node_gid,
                          const size_t bdy_set)
@@ -87,4 +87,4 @@ namespace PistonVelocityBC {
 
 } // end namespace
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/boundary_conditions/velocity/reflected_velocity_bc.h b/single-node-refactor/src/boundary_conditions/velocity/reflected_velocity_bc.h
index 3535f81b6..912461f81 100644
--- a/single-node-refactor/src/boundary_conditions/velocity/reflected_velocity_bc.h
+++ b/single-node-refactor/src/boundary_conditions/velocity/reflected_velocity_bc.h
@@ -64,7 +64,7 @@ namespace ReflectedVelocityBC {
                          const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                          const DCArrayKokkos<double>& bc_global_vars,
                          const DCArrayKokkos<double>& bc_state_vars,
-                         const DCArrayKokkos<double>& node_vel,
+                         const MPIArrayKokkos<double>& node_vel,
                          const double time_value,
                          const size_t bdy_node_gid,
                          const size_t bdy_set)
@@ -84,4 +84,4 @@ namespace ReflectedVelocityBC {
 } // end namespace
 
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/boundary_conditions/velocity/time_varying_velocity_bc.h b/single-node-refactor/src/boundary_conditions/velocity/time_varying_velocity_bc.h
index 279e69db0..ca12152f5 100644
--- a/single-node-refactor/src/boundary_conditions/velocity/time_varying_velocity_bc.h
+++ b/single-node-refactor/src/boundary_conditions/velocity/time_varying_velocity_bc.h
@@ -75,7 +75,7 @@ namespace TimeVaryingVelocityBC {
                          const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                          const DCArrayKokkos<double>& bc_global_vars,
                          const DCArrayKokkos<double>& bc_state_vars,
-                         const DCArrayKokkos<double>& node_vel,
+                         const MPIArrayKokkos<double>& node_vel,
                          const double time_value,
                          const size_t bdy_node_gid,
                          const size_t bdy_set)
@@ -108,4 +108,4 @@ namespace TimeVaryingVelocityBC {
 
 } // end namespace
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/boundary_conditions/velocity/user_defined_velocity_bc.h b/single-node-refactor/src/boundary_conditions/velocity/user_defined_velocity_bc.h
index 3f8ecdfb3..46aeb80e5 100644
--- a/single-node-refactor/src/boundary_conditions/velocity/user_defined_velocity_bc.h
+++ b/single-node-refactor/src/boundary_conditions/velocity/user_defined_velocity_bc.h
@@ -65,7 +65,7 @@ namespace UserDefinedVelocityBC {
                          const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                          const DCArrayKokkos<double>& bc_global_vars,
                          const DCArrayKokkos<double>& bc_state_vars,
-                         const DCArrayKokkos<double>& node_vel,
+                         const MPIArrayKokkos<double>& node_vel,
                          const double time_value,
                          const size_t bdy_node_gid,
                          const size_t bdy_set)
@@ -78,4 +78,4 @@ namespace UserDefinedVelocityBC {
 
 } // end namespace
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/boundary_conditions/velocity/zero_velocity_bc.h b/single-node-refactor/src/boundary_conditions/velocity/zero_velocity_bc.h
index 10a893aeb..b07808f9c 100644
--- a/single-node-refactor/src/boundary_conditions/velocity/zero_velocity_bc.h
+++ b/single-node-refactor/src/boundary_conditions/velocity/zero_velocity_bc.h
@@ -66,7 +66,7 @@ namespace ZeroVelocityBC {
                          const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                          const DCArrayKokkos<double>& bc_global_vars,
                          const DCArrayKokkos<double>& bc_state_vars,
-                         const DCArrayKokkos<double>& node_vel,
+                         const MPIArrayKokkos<double>& node_vel,
                          const double time_value,
                          const size_t bdy_node_gid,
                          const size_t bdy_set)
@@ -82,4 +82,4 @@ namespace ZeroVelocityBC {
 
 } // end namespace
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/common/boundary_conditions.h b/single-node-refactor/src/common/boundary_conditions.h
index d26272fea..8919a5677 100644
--- a/single-node-refactor/src/common/boundary_conditions.h
+++ b/single-node-refactor/src/common/boundary_conditions.h
@@ -189,7 +189,7 @@ struct BoundaryConditionFunctions_t
                           const DCArrayKokkos <BoundaryConditionEnums_t>& BoundaryConditionEnums,
                           const DCArrayKokkos<double>& bc_global_vars,
                           const DCArrayKokkos<double>& bc_state_vars,
-                          const DCArrayKokkos<double>& node_vel,
+                          const MPIArrayKokkos<double>& node_vel,
                           const double time_value,
                           const size_t bdy_node_gid,
                           const size_t bdy_set) = NULL;
@@ -259,4 +259,4 @@ static std::vector<std::string> bc_required_inps
 {
 };
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/common/io_utils.h b/single-node-refactor/src/common/io_utils.h
index c39d4769a..ad01f3e58 100644
--- a/single-node-refactor/src/common/io_utils.h
+++ b/single-node-refactor/src/common/io_utils.h
@@ -354,6 +354,7 @@ class MeshBuilder
         const int num_elems = num_elems_i * num_elems_j;
 
         std::vector<double> origin(num_dim);
+
         //SimulationParamaters.mesh_input.origin.update_host();
         for (int i = 0; i < num_dim; i++) { origin[i] = SimulationParamaters.mesh_input.origin[i]; }
 
@@ -378,13 +379,35 @@ class MeshBuilder
         mesh.initialize_nodes(num_nodes);
         node.initialize(rk_num_bins, num_nodes, num_dim);
 
+        // MPI data
+        int world_size, rank;
+        MPI_Comm_size(MPI_COMM_WORLD, &world_size);
+        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+        int grid_world_size, jrem, irem, col, row, jadd, iadd, jsub, isub, jchunk, ichunk, jstart, jend, istart, iend;
+        grid_world_size = sqrt(world_size);
+        jrem = num_points_j % grid_world_size; // how many extra points we have after a division of ranks
+        irem = num_points_i % grid_world_size;
+        col = rank % grid_world_size;
+        row = rank / grid_world_size;
+        jadd = col / (grid_world_size - jrem); // 1 if you have an extra piece, 0 otherwise
+        iadd = row / (grid_world_size - irem); // 1 if you have an extra piece, 0 otherwise
+        jsub = jadd * (grid_world_size - jrem); // subtraction from your start based on how many other ranks have an extra piece
+        isub = iadd * (grid_world_size - irem);
+        jchunk = num_points_j / grid_world_size + jadd;
+        ichunk = num_points_i / grid_world_size + iadd;
+        jstart = col * jchunk - jsub;
+        jend = jstart +  jchunk;
+        istart = row * ichunk - isub;
+        iend = istart +  ichunk;
+
         // --- Build nodes ---
 
         // populate the point data structures
-        for (int j = 0; j < num_points_j; j++) {
-            for (int i = 0; i < num_points_i; i++) {
+        for (int j = jstart, j_loc = 0; j < jend; j++, j_loc++) {
+            for (int i = istart, i_loc = 0; i < iend; i++, i_loc++) {
                 // global id for the point
-                int node_gid = get_id(i, j, 0, num_points_i, num_points_j);
+                int node_gid = get_id(i_loc, j_loc, 0, ichunk, jchunk);
 
                 // store the point coordinates
                 node.coords.host(0, node_gid, 0) = origin[0] + (double)i * dx;
@@ -393,7 +416,7 @@ class MeshBuilder
         } // end for j
 
         for (int rk_level = 1; rk_level < rk_num_bins; rk_level++) {
-            for (int node_gid = 0; node_gid < num_nodes; node_gid++) {
+            for (int node_gid = 0; node_gid < jchunk * ichunk; node_gid++) {
                 node.coords.host(rk_level, node_gid, 0) = node.coords.host(0, node_gid, 0);
                 node.coords.host(rk_level, node_gid, 1) = node.coords.host(0, node_gid, 1);
             }
@@ -401,15 +424,31 @@ class MeshBuilder
         node.coords.update_device();
 
         // intialize elem variables
+        
+        jrem = (int) num_elems_j % grid_world_size; // how many extra points we have after a division of ranks
+        irem = (int) num_elems_i % grid_world_size;
+        col = rank % grid_world_size;
+        row = rank / grid_world_size;
+        jadd = col / (grid_world_size - jrem); // 1 if you have an extra piece, 0 otherwise
+        iadd = row / (grid_world_size - irem); // 1 if you have an extra piece, 0 otherwise
+        jsub = jadd * (grid_world_size - jrem); // subtraction from your start based on how many other ranks have an extra piece
+        isub = iadd * (grid_world_size - irem);
+        jchunk = (int) num_elems_j / grid_world_size + jadd;
+        ichunk = (int) num_elems_i / grid_world_size + iadd;
+        jstart = col * jchunk - jsub;
+        jend = jstart +  jchunk;
+        istart = row * ichunk - isub;
+        iend = istart +  ichunk;
+
         mesh.initialize_elems(num_elems, num_dim);
         MaterialPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D
         GaussPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D
 
         // populate the elem center data structures
-        for (int j = 0; j < num_elems_j; j++) {
-            for (int i = 0; i < num_elems_i; i++) {
+        for (int j = jstart, j_loc = 0; j < jend; j++, j_loc++) {
+            for (int i = istart, i_loc = 0; i < i<iend; i++, i_loc++) {
                 // global id for the elem
-                int elem_gid = get_id(i, j, 0, num_elems_i, num_elems_j);
+                int elem_gid = get_id(i_loc, j_loc, 0, ichunk, jchunk);
 
                 // store the point IDs for this elem where the range is
                 // (i:i+1, j:j+1 for a linear quad
@@ -1288,4 +1327,4 @@ class MeshWriter
     }
 };
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/common/material.h b/single-node-refactor/src/common/material.h
index 2b82de1c9..e1933461f 100644
--- a/single-node-refactor/src/common/material.h
+++ b/single-node-refactor/src/common/material.h
@@ -222,7 +222,7 @@ struct MaterialFunctions_t
                         const ViewCArrayKokkos<double>& vel_grad,
                         const ViewCArrayKokkos<size_t>& MaterialPoints_node_gids,
                         const DCArrayKokkos<double>&    node_coords,
-                        const DCArrayKokkos<double>&    node_vel,
+                        const MPIArrayKokkos<double>&    node_vel,
                         const double vol,
                         const double dt,
                         const double rk_alpha) = NULL;
@@ -359,4 +359,4 @@ static std::vector<std::string> material_thermal_statics_required_inps
 };
 
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/common/mesh.h b/single-node-refactor/src/common/mesh.h
index a491619ed..cf4109473 100644
--- a/single-node-refactor/src/common/mesh.h
+++ b/single-node-refactor/src/common/mesh.h
@@ -1537,7 +1537,7 @@ void decompose_vel_grad(ViewCArrayKokkos<double>& D_tensor,
                         const ViewCArrayKokkos<size_t>& elem_node_gids,
                         const size_t elem_gid,
                         const DCArrayKokkos<double>& node_coords,
-                        const DCArrayKokkos<double>& node_vel,
+                        const MPIArrayKokkos<double>& node_vel,
                         const double vol);
 
-#endif
\ No newline at end of file
+#endif
diff --git a/single-node-refactor/src/common/mesh_old.h b/single-node-refactor/src/common/mesh_old.h
index e79041ea9..6ca4f5eae 100644
--- a/single-node-refactor/src/common/mesh_old.h
+++ b/single-node-refactor/src/common/mesh_old.h
@@ -919,7 +919,7 @@ struct mesh_t
 
 void ensight(const mesh_t& mesh,
              const DCArrayKokkos<double>& node_coords,
-             const DCArrayKokkos<double>& node_vel,
+             const MPIArrayKokkos<double>& node_vel,
              const DCArrayKokkos<double>& node_mass,
              const DCArrayKokkos<double>& elem_den,
              const DCArrayKokkos<double>& elem_pres,
@@ -935,7 +935,7 @@ void ensight(const mesh_t& mesh,
 
 void state_file(const mesh_t& mesh,
                 const DCArrayKokkos<double>& node_coords,
-                const DCArrayKokkos<double>& node_vel,
+                const MPIArrayKokkos<double>& node_vel,
                 const DCArrayKokkos<double>& node_mass,
                 const DCArrayKokkos<double>& elem_den,
                 const DCArrayKokkos<double>& elem_pres,
@@ -959,7 +959,7 @@ void decompose_vel_grad(ViewCArrayKokkos<double>& D_tensor,
                         const ViewCArrayKokkos<size_t>& elem_node_gids,
                         const size_t elem_gid,
                         const DCArrayKokkos<double>& node_coords,
-                        const DCArrayKokkos<double>& node_vel,
+                        const MPIArrayKokkos<double>& node_vel,
                         const double vol);
 
 #endif
diff --git a/single-node-refactor/src/common/state.h b/single-node-refactor/src/common/state.h
index 39057b8ff..9ac2f1661 100644
--- a/single-node-refactor/src/common/state.h
+++ b/single-node-refactor/src/common/state.h
@@ -48,14 +48,14 @@ using namespace mtr;
 struct node_t
 {
     DCArrayKokkos<double> coords; ///< Nodal coordinates
-    DCArrayKokkos<double> vel;  ///< Nodal velocity
+    MPIArrayKokkos<double> vel;  ///< Nodal velocity
     DCArrayKokkos<double> mass; ///< Nodal mass
 
     // initialization method (num_rk_storage_bins, num_nodes, num_dims)
     void initialize(size_t num_rk, size_t num_nodes, size_t num_dims)
     {
         this->coords = DCArrayKokkos<double>(num_rk, num_nodes, num_dims, "node_coordinates");
-        this->vel    = DCArrayKokkos<double>(num_rk, num_nodes, num_dims, "node_velocity");
+        this->vel    = MPIArrayKokkos<double>(num_rk, num_nodes, num_dims, "node_velocity");
         this->mass   = DCArrayKokkos<double>(num_nodes, "node_mass");
     }; // end method
 }; // end node_t
diff --git a/single-node-refactor/src/main.cpp b/single-node-refactor/src/main.cpp
index 01ce598e2..d68748eef 100644
--- a/single-node-refactor/src/main.cpp
+++ b/single-node-refactor/src/main.cpp
@@ -60,8 +60,10 @@ int main(int argc, char* argv[])
         return 0;
     } // end if
 
-    Kokkos::initialize();
-    {
+    //Kokkos::initialize(argc, argv);
+    MATAR_MPI_INIT
+    MATAR_KOKKOS_INIT
+    { // kokkos scope
 
         // Create driver
         Driver* driver = new Driver(argv[1]);
@@ -73,9 +75,10 @@ int main(int argc, char* argv[])
 
         // Delete driver
         delete driver;
-    }
-
-    Kokkos::finalize();
+    } // end kokkos scope
+    MATAR_KOKKOS_FINALIZE
+    MATAR_MPI_FINALIZE
+    //Kokkos::finalize();
 
     std::cout << "**** End of main **** " << std::endl;
     return 0;
diff --git a/single-node-refactor/src/material_models/strength/no_strength.h b/single-node-refactor/src/material_models/strength/no_strength.h
index 957c9c37f..488995d65 100644
--- a/single-node-refactor/src/material_models/strength/no_strength.h
+++ b/single-node-refactor/src/material_models/strength/no_strength.h
@@ -76,7 +76,7 @@ namespace UserDefinedStrengthModel {
         const ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
         const DCArrayKokkos<double>&    node_coords,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const double vol,
         const double dt,
         const double rk_alpha)
@@ -111,7 +111,7 @@ namespace NoStrengthModel {
         const ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
         const DCArrayKokkos<double>&    node_coords,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const double vol,
         const double dt,
         const double rk_alpha)
@@ -155,7 +155,7 @@ namespace NotionalStrengthModel {
         const ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
         const DCArrayKokkos<double>&    node_coords,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const double vol,
         const double dt,
         const double rk_alpha)
@@ -168,4 +168,4 @@ namespace NotionalStrengthModel {
 
 
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard
diff --git a/single-node-refactor/src/material_models/strength/user_defined_strength.h b/single-node-refactor/src/material_models/strength/user_defined_strength.h
index 8922cb007..6a18b0d74 100644
--- a/single-node-refactor/src/material_models/strength/user_defined_strength.h
+++ b/single-node-refactor/src/material_models/strength/user_defined_strength.h
@@ -76,7 +76,7 @@ namespace UserDefinedStrengthModel {
         const ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
         const DCArrayKokkos<double>&    node_coords,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const double vol,
         const double dt,
         const double rk_alpha)
@@ -129,7 +129,7 @@ namespace NotionalStrengthModel {
         const ViewCArrayKokkos<double>& vel_grad,
         const ViewCArrayKokkos<size_t>& elem_node_gids,
         const DCArrayKokkos<double>&    node_coords,
-        const DCArrayKokkos<double>&    node_vel,
+        const MPIArrayKokkos<double>&    node_vel,
         const double vol,
         const double dt,
         const double rk_alpha)
@@ -142,4 +142,4 @@ namespace NotionalStrengthModel {
 
 
 
-#endif // end Header Guard
\ No newline at end of file
+#endif // end Header Guard