diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 17cc5f92..4faf395d 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -118,7 +118,7 @@ jobs:
       matrix:
         cxx: ['hipcc']
         cmake_build_type: ['Release']
-        kokkos_ver: ['3.6.01']
+        kokkos_ver: ['4.0.01']
     runs-on: ubuntu-20.04
     container: ghcr.io/ecp-copa/ci-containers/rocm:latest
     steps:
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index df1b8e15..79a697a6 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -7,4 +7,7 @@ target_link_libraries(KalthoffWinkler LINK_PUBLIC CabanaPD)
 add_executable(CrackBranching crack_branching.cpp)
 target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD)
 
+add_executable(ReadFile read_file.cpp)
+target_link_libraries(ReadFile LINK_PUBLIC CabanaPD)
+
 install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR})
diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp
index e4c8ddd0..cc045cc3 100644
--- a/examples/crack_branching.cpp
+++ b/examples/crack_branching.cpp
@@ -106,9 +106,11 @@ int main( int argc, char* argv[] )
         CabanaPD::RegionBoundary plane2( low_x, high_x, high_y - dy,
                                          high_y + dy, low_z, high_z );
         std::vector<CabanaPD::RegionBoundary> planes = { plane1, plane2 };
-        auto bc =
-            createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{},
-                                     exec_space{}, *particles, planes, b0 );
+        int bc_dim = 1;
+        double center = 0.0;
+        auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{},
+                                           exec_space{}, *particles, planes, b0,
+                                           bc_dim, center );
 
         auto init_functor = KOKKOS_LAMBDA( const int pid )
         {
diff --git a/examples/read_file.cpp b/examples/read_file.cpp
new file mode 100644
index 00000000..e04f95b6
--- /dev/null
+++ b/examples/read_file.cpp
@@ -0,0 +1,201 @@
+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+template <class ParticleType>
+void read_particles( const std::string filename, ParticleType& particles )
+{
+    // Read particles from csv file.
+    std::vector<double> csv_x;
+    std::vector<double> csv_y;
+    std::vector<double> csv_v;
+
+    std::vector<std::string> row;
+    std::string line, word;
+
+    std::fstream file( filename, std::ios::in );
+    if ( file.is_open() )
+    {
+        std::getline( file, line );
+        while ( std::getline( file, line ) )
+        {
+            row.clear();
+
+            std::stringstream str( line );
+
+            while ( std::getline( str, word, ',' ) )
+            {
+                row.push_back( word );
+            }
+            csv_x.push_back( std::stod( row[0] ) );
+            csv_y.push_back( std::stod( row[1] ) );
+            csv_v.push_back( std::stod( row[2] ) );
+        }
+    }
+    else
+        throw( "Could not open file." + filename );
+
+    // Create unmanaged Views in order to copy to device.
+    Kokkos::View<double*, Kokkos::HostSpace> x_host( csv_x.data(),
+                                                     csv_x.size() );
+    Kokkos::View<double*, Kokkos::HostSpace> y_host( csv_y.data(),
+                                                     csv_y.size() );
+    Kokkos::View<double*, Kokkos::HostSpace> vol_host( csv_v.data(),
+                                                       csv_v.size() );
+
+    // Copy to the device.
+    using memory_space = typename ParticleType::memory_space;
+    auto x = Kokkos::create_mirror_view_and_copy( memory_space(), x_host );
+    auto y = Kokkos::create_mirror_view_and_copy( memory_space(), y_host );
+    auto vol = Kokkos::create_mirror_view_and_copy( memory_space(), vol_host );
+    // Resize internal variables (no ghosts initially).
+    particles.resize( x.size(), 0 );
+
+    auto px = particles.sliceReferencePosition();
+    auto pvol = particles.sliceVolume();
+    auto v = particles.sliceVelocity();
+    auto f = particles.sliceForce();
+    auto type = particles.sliceType();
+    auto rho = particles.sliceDensity();
+    auto u = particles.sliceDisplacement();
+    auto nofail = particles.sliceNoFail();
+    auto damage = particles.sliceDamage();
+
+    // Using atomic memory for the count.
+    using CountView =
+        typename Kokkos::View<int, Kokkos::LayoutRight, memory_space,
+                              Kokkos::MemoryTraits<Kokkos::Atomic>>;
+    CountView count( "count" );
+
+    using exec_space = typename memory_space::execution_space;
+    Kokkos::parallel_for(
+        "copy_to_particles", Kokkos::RangePolicy<exec_space>( 0, x.size() ),
+        KOKKOS_LAMBDA( const int pid ) {
+            // Guard against duplicate points. This threaded loop should be
+            // faster than checking during the host read.
+            for ( int p = 0; p < pid; p++ )
+            {
+                auto xdiff = x( p ) - x( pid );
+                auto ydiff = y( p ) - y( pid );
+                if ( ( Kokkos::abs( xdiff ) < 1e-14 ) &&
+                     ( Kokkos::abs( ydiff ) < 1e-14 ) )
+                    return;
+            }
+            const std::size_t c = count()++;
+
+            // Set the particle position and volume.
+            pvol( c ) = vol( pid );
+            px( c, 0 ) = x( pid );
+            px( c, 1 ) = y( pid );
+
+            // Initialize everything else to zero.
+            // Currently psuedo-2d.
+            px( c, 2 ) = 0.0;
+            for ( int d = 0; d < 3; d++ )
+            {
+                u( c, d ) = 0.0;
+                v( c, d ) = 0.0;
+                f( c, d ) = 0.0;
+            }
+            type( c ) = 0;
+            nofail( c ) = 0;
+            rho( c ) = 1.0;
+            damage( c ) = 0;
+        } );
+
+    int final_count = 0;
+    Kokkos::deep_copy( final_count, count );
+    if ( final_count < static_cast<int>( x.size() ) )
+        particles.resize( final_count, 0 );
+}
+
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    {
+        Kokkos::ScopeGuard scope_guard( argc, argv );
+
+        // FIXME: change backend at compile time for now.
+        using exec_space = Kokkos::DefaultExecutionSpace;
+        using memory_space = typename exec_space::memory_space;
+
+        // Time
+        double t_final = 1e-1;
+        double t_ramp = 1e-3;
+        double dt = 1e-7;
+        double output_frequency = 100;
+
+        // Material constants
+        double E = 72e+9;                      // [Pa]
+        double nu = 0.25;                      // unitless
+        double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa]
+        double G0 = 3.8;                       // [J/m^2]
+
+        // PD horizon
+        int m = 3;
+        double delta = 4.0;
+        double dx = delta / m;
+
+        // Choose force model type.
+        using model_type =
+            CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
+        model_type force_model( delta, K, G0 );
+        CabanaPD::Inputs<3> inputs( t_final, dt, output_frequency, true );
+        inputs.read_args( argc, argv );
+
+        // Default construct to then read particles.
+        auto particles = std::make_shared<CabanaPD::Particles<
+            memory_space, typename model_type::base_model, 3>>();
+
+        // Read particles from file.
+        read_particles( inputs.input_file, *particles );
+        // Update after reading. Currently requires fixed cell spacing.
+        particles->updateAfterRead( exec_space(), dx );
+
+        // Do this separately so that the mesh is already set up.
+        auto x = particles->sliceReferencePosition();
+        auto nofail = particles->sliceNoFail();
+        auto f = particles->sliceForce();
+        double max_bc = particles->local_mesh_hi[1] - delta;
+        double min_bc = particles->local_mesh_lo[1] + delta;
+        auto init_functor = KOKKOS_LAMBDA( const int pid )
+        {
+            if ( x( pid, 1 ) <= min_bc + 1e-10 ||
+                 x( pid, 1 ) >= max_bc - 1e-10 )
+                nofail( pid ) = 1;
+        };
+        particles->updateParticles( exec_space{}, init_functor );
+
+        CabanaPD::Prenotch<0> prenotch;
+
+        CabanaPD::RegionBoundary planeUpper(
+            particles->local_mesh_lo[0], particles->local_mesh_hi[0], max_bc,
+            particles->local_mesh_hi[1], particles->local_mesh_lo[2],
+            particles->local_mesh_hi[2] );
+        CabanaPD::RegionBoundary planeLower(
+            particles->local_mesh_lo[0], particles->local_mesh_hi[0],
+            particles->local_mesh_lo[1], min_bc, particles->local_mesh_lo[2],
+            particles->local_mesh_hi[2] );
+
+        std::vector<CabanaPD::RegionBoundary> planes = { planeUpper,
+                                                         planeLower };
+        int bc_dim = 1;
+        double center = particles->local_mesh_ext[bc_dim] / 2.0 +
+                        particles->local_mesh_lo[bc_dim];
+        auto bc = createBoundaryCondition(
+            CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes,
+            2e4, 0.0, bc_dim, center, 0.0, t_ramp );
+
+        auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
+            inputs, particles, force_model, bc, prenotch );
+        cabana_pd->init_force();
+        cabana_pd->run();
+    }
+
+    MPI_Finalize();
+}
diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp
index 2e45522f..32912c42 100644
--- a/src/CabanaPD_Boundary.hpp
+++ b/src/CabanaPD_Boundary.hpp
@@ -42,7 +42,12 @@ struct RegionBoundary
         , low_y( _low_y )
         , high_y( _high_y )
         , low_z( _low_z )
-        , high_z( _high_z ){};
+        , high_z( _high_z )
+    {
+        assert( low_x < high_x );
+        assert( low_y < high_y );
+        assert( low_z < high_z );
+    };
 };
 
 template <class MemorySpace, class BoundaryType>
@@ -118,7 +123,7 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
 
 template <class BoundaryType, class ExecSpace, class Particles>
 auto createBoundaryIndexSpace( ExecSpace exec_space, Particles particles,
-                               std::vector<RegionBoundary> planes,
+                               std::vector<BoundaryType> planes,
                                const double initial_guess )
 {
     using memory_space = typename Particles::memory_space;
@@ -133,22 +138,32 @@ struct ForceUpdateBCTag
 {
 };
 
-struct ForceCrackBranchBCTag
+struct ForceSymmetric1dBCTag
 {
 };
 
 template <class BCIndexSpace, class BCTag>
 struct BoundaryCondition;
 
+// Reset BC (all dimensions).
 template <class BCIndexSpace>
 struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
 {
     double _value;
     BCIndexSpace _index_space;
-
-    BoundaryCondition( const double value, BCIndexSpace bc_index_space )
+    double _time_start;
+    double _time_ramp;
+    double _time_end;
+
+    BoundaryCondition( BCIndexSpace bc_index_space, const double value,
+                       const double start = 0.0,
+                       const double ramp = std::numeric_limits<double>::max(),
+                       const double end = std::numeric_limits<double>::max() )
         : _value( value )
         , _index_space( bc_index_space )
+        , _time_start( start )
+        , _time_ramp( ramp )
+        , _time_end( end )
     {
     }
 
@@ -159,33 +174,56 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
         _index_space.update( exec_space, particles, plane );
     }
 
-    template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType& particles )
+    template <class ExecSpace, class FieldType, class PositionType>
+    void apply( ExecSpace, FieldType& f, PositionType&, const double t )
     {
-        auto f = particles.sliceForce();
-        auto x = particles.sliceReferencePosition();
         auto index_space = _index_space._view;
         Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
         auto value = _value;
+        auto start = _time_start;
+        auto end = _time_end;
+        auto ramp = _time_ramp;
         Kokkos::parallel_for(
             "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) {
-                auto pid = index_space( b );
-                for ( int d = 0; d < 3; d++ )
-                    f( pid, d ) = value;
+                double current = value;
+                // Initial linear ramp.
+                if ( t < ramp )
+                {
+                    current = value * ( t - start ) / ( ramp - start );
+                }
+                if ( t > start && t < end )
+                {
+                    auto pid = index_space( b );
+                    for ( int d = 0; d < 3; d++ )
+                        f( pid, d ) = current;
+                }
             } );
     }
 };
 
+// Increment BC (all dimensions).
 template <class BCIndexSpace>
 struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
 {
     double _value;
     BCIndexSpace _index_space;
-
-    BoundaryCondition( const double value, BCIndexSpace bc_index_space )
+    double _time_start;
+    double _time_ramp;
+    double _time_end;
+
+    BoundaryCondition( BCIndexSpace bc_index_space, const double value,
+                       const double start = 0.0,
+                       const double ramp = std::numeric_limits<double>::max(),
+                       const double end = std::numeric_limits<double>::max() )
         : _value( value )
         , _index_space( bc_index_space )
+        , _time_start( start )
+        , _time_ramp( ramp )
+        , _time_end( end )
     {
+        assert( _time_ramp >= _time_start );
+        assert( _time_end >= _time_ramp );
+        assert( _time_end >= _time_start );
     }
 
     template <class ExecSpace, class Particles>
@@ -195,32 +233,88 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
         _index_space.update( exec_space, particles, plane );
     }
 
-    template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType& particles )
+    template <class ExecSpace, class FieldType, class PositionType>
+    void apply( ExecSpace, FieldType& f, PositionType&, const double t )
     {
-        auto f = particles.sliceForce();
         auto index_space = _index_space._view;
         Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
         auto value = _value;
+        auto start = _time_start;
+        auto end = _time_end;
+        auto ramp = _time_ramp;
         Kokkos::parallel_for(
             "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) {
-                auto pid = index_space( b );
-                for ( int d = 0; d < 3; d++ )
-                    f( pid, d ) += value;
+                double current = value;
+                // Initial linear ramp.
+                if ( t < ramp )
+                {
+                    current = value * ( t - start ) / ( ramp - start );
+                }
+                if ( t > start && t < end )
+                {
+                    auto pid = index_space( b );
+                    for ( int d = 0; d < 3; d++ )
+                        f( pid, d ) += current;
+                }
             } );
     }
 };
 
+// Symmetric 1d BC applied with opposite sign based on position from the
+// midplane in the specified direction.
 template <class BCIndexSpace>
-struct BoundaryCondition<BCIndexSpace, ForceCrackBranchBCTag>
+struct BoundaryCondition<BCIndexSpace, ForceSymmetric1dBCTag>
 {
-    double _value;
     BCIndexSpace _index_space;
+    double _value_top;
+    double _value_bottom;
+    double _time_start;
+    double _time_ramp;
+    double _time_end;
+
+    int _dim;
+    double _center;
+
+    BoundaryCondition( BCIndexSpace bc_index_space, const double value,
+                       const int dim, const double center,
+                       const double start = 0.0,
+                       const double ramp = std::numeric_limits<double>::max(),
+                       const double end = std::numeric_limits<double>::max() )
+        : _index_space( bc_index_space )
+        , _value_top( value )
+        , _value_bottom( value )
+        , _time_start( start )
+        , _time_ramp( ramp )
+        , _time_end( end )
+        , _dim( dim )
+        , _center( center )
+    {
+        assert( _time_ramp >= _time_start );
+        assert( _time_end >= _time_ramp );
+        assert( _time_end >= _time_start );
+        assert( _dim >= 0 );
+        assert( _dim < 3 );
+    }
 
-    BoundaryCondition( const double value, BCIndexSpace bc_index_space )
-        : _value( value )
-        , _index_space( bc_index_space )
+    BoundaryCondition( BCIndexSpace bc_index_space, const double value_top,
+                       const double value_bottom, const int dim,
+                       const double center, const double start = 0.0,
+                       const double ramp = std::numeric_limits<double>::max(),
+                       const double end = std::numeric_limits<double>::max() )
+        : _index_space( bc_index_space )
+        , _value_top( value_top )
+        , _value_bottom( value_bottom )
+        , _time_start( start )
+        , _time_ramp( ramp )
+        , _time_end( end )
+        , _dim( dim )
+        , _center( center )
     {
+        assert( _time_ramp >= _time_start );
+        assert( _time_end >= _time_ramp );
+        assert( _time_end >= _time_start );
+        assert( _dim >= 0 );
+        assert( _dim < 3 );
     }
 
     template <class ExecSpace, class Particles>
@@ -230,38 +324,96 @@ struct BoundaryCondition<BCIndexSpace, ForceCrackBranchBCTag>
         _index_space.update( exec_space, particles, plane );
     }
 
-    template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType& particles )
+    template <class ExecSpace, class FieldType, class PositionType>
+    void apply( ExecSpace, FieldType& f, PositionType& x, const double t )
     {
-        auto f = particles.sliceForce();
-        auto x = particles.sliceReferencePosition();
         auto index_space = _index_space._view;
         Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
-        auto value = _value;
+        auto dim = _dim;
+        auto center = _center;
+        auto value_top = _value_top;
+        auto value_bottom = _value_bottom;
+        auto start = _time_start;
+        auto end = _time_end;
+        auto ramp = _time_ramp;
         Kokkos::parallel_for(
             "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) {
-                auto pid = index_space( b );
-                // This is specifically for the crack branching.
-                auto sign = std::abs( x( pid, 1 ) ) / x( pid, 1 );
-                f( pid, 1 ) += value * sign;
+                auto current_top = value_top;
+                auto current_bottom = value_bottom;
+                // Initial linear ramp.
+                if ( t < ramp )
+                {
+                    auto t_factor = ( t - start ) / ( ramp - start );
+                    current_top = value_top * t_factor;
+                    current_bottom = value_bottom * t_factor;
+                }
+                if ( t > start && t < end )
+                {
+                    auto pid = index_space( b );
+                    if ( x( pid, dim ) > center )
+                        f( pid, dim ) += current_top;
+                    else
+                        f( pid, dim ) += current_bottom;
+                }
             } );
     }
 };
 
 // FIXME: relatively large initial guess for allocation.
-template <class BoundaryType, class BCTag, class ExecSpace, class Particles>
-auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles,
-                              std::vector<BoundaryType> planes,
-                              const double value,
-                              const double initial_guess = 0.5 )
+template <class BoundaryType, class BCTag, class ExecSpace, class Particles,
+          class... Args>
+auto createBoundaryCondition(
+    BCTag, ExecSpace exec_space, Particles particles,
+    std::vector<BoundaryType> planes, const double value,
+    const double start = 0.0,
+    const double ramp = std::numeric_limits<double>::max(),
+    const double end = std::numeric_limits<double>::max(),
+    const double initial_guess = 0.5 )
 {
     using memory_space = typename Particles::memory_space;
     using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
     bc_index_type bc_indices = createBoundaryIndexSpace<BoundaryType>(
         exec_space, particles, planes, initial_guess );
-    return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
+    return BoundaryCondition<bc_index_type, BCTag>( bc_indices, value, start,
+                                                    ramp, end );
 }
 
+// FIXME: relatively large initial guess for allocation.
+template <class BoundaryType, class ExecSpace, class Particles>
+auto createBoundaryCondition(
+    ForceSymmetric1dBCTag, ExecSpace exec_space, Particles particles,
+    std::vector<BoundaryType> planes, const double value, const int dim,
+    const double center, const double start = 0.0,
+    const double ramp = std::numeric_limits<double>::max(),
+    const double end = std::numeric_limits<double>::max(),
+    const double initial_guess = 0.5 )
+{
+    using memory_space = typename Particles::memory_space;
+    using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
+    bc_index_type bc_indices = createBoundaryIndexSpace(
+        exec_space, particles, planes, initial_guess );
+    return BoundaryCondition<bc_index_type, ForceSymmetric1dBCTag>(
+        bc_indices, value, dim, center, start, ramp, end );
+}
+
+// FIXME: relatively large initial guess for allocation.
+template <class BoundaryType, class ExecSpace, class Particles>
+auto createBoundaryCondition(
+    ForceSymmetric1dBCTag, ExecSpace exec_space, Particles particles,
+    std::vector<BoundaryType> planes, const double value_top,
+    const double value_bottom, const int dim, const double center,
+    const double start = 0.0,
+    const double ramp = std::numeric_limits<double>::max(),
+    const double end = std::numeric_limits<double>::max(),
+    const double initial_guess = 0.5 )
+{
+    using memory_space = typename Particles::memory_space;
+    using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
+    bc_index_type bc_indices = createBoundaryIndexSpace(
+        exec_space, particles, planes, initial_guess );
+    return BoundaryCondition<bc_index_type, ForceSymmetric1dBCTag>(
+        bc_indices, value_top, value_bottom, dim, center, start, ramp, end );
+}
 } // namespace CabanaPD
 
 #endif
diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp
index 51b28f59..544feabe 100644
--- a/src/CabanaPD_Input.hpp
+++ b/src/CabanaPD_Input.hpp
@@ -24,6 +24,7 @@ class Inputs
   public:
     std::string output_file = "cabanaPD.out";
     std::string error_file = "cabanaPD.err";
+    std::string input_file = "particles.csv";
     std::string device_type = "SERIAL";
 
     std::array<int, Dim> num_cells;
@@ -51,6 +52,15 @@ class Inputs
     {
         num_steps = final_time / timestep;
     }
+    Inputs( const double t_f, const double dt, const int of,
+            const bool output_ref )
+        : final_time( t_f )
+        , timestep( dt )
+        , output_frequency( of )
+        , output_reference( output_ref )
+    {
+        num_steps = final_time / timestep;
+    }
     ~Inputs(){};
 
     void read_args( int argc, char* argv[] )
@@ -87,6 +97,12 @@ class Inputs
                 output_file = argv[i + 1];
                 ++i;
             }
+            else if ( ( strcmp( argv[i], "-i" ) == 0 ) ||
+                      ( strcmp( argv[i], "--input-file" ) == 0 ) )
+            {
+                input_file = argv[i + 1];
+                ++i;
+            }
             else if ( ( strcmp( argv[i], "-e" ) == 0 ) ||
                       ( strcmp( argv[i], "--error-file" ) == 0 ) )
             {
diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp
index 5ec97f1f..55e33a29 100644
--- a/src/CabanaPD_Particles.hpp
+++ b/src/CabanaPD_Particles.hpp
@@ -281,6 +281,136 @@ class Particles<MemorySpace, PMB, Dimension>
         resize( n_local, 0 );
         size = _aosoa_x.size();
 
+        update_global();
+    }
+
+    // This is necessary after reading in particles from file for a consistent
+    // state.
+    template <class ExecSpace, std::size_t NSD = dim>
+    std::enable_if_t<3 == NSD, void>
+    updateAfterRead( const ExecSpace& exec_space, const double dx )
+    {
+        // Because this may be a non-uniform mesh, build a background mesh with
+        // dx=cutoff and force halo_width=1
+        halo_width = 1;
+        auto x = sliceReferencePosition();
+        n_local = x.size();
+        n_ghost = 0;
+        size = n_local;
+        update_global();
+
+        double max_x;
+        Kokkos::Max<double> max_x_reducer( max_x );
+        double min_x;
+        Kokkos::Min<double> min_x_reducer( min_x );
+        double max_y;
+        Kokkos::Max<double> max_y_reducer( max_y );
+        double min_y;
+        Kokkos::Min<double> min_y_reducer( min_y );
+        double max_z;
+        Kokkos::Max<double> max_z_reducer( max_z );
+        double min_z;
+        Kokkos::Min<double> min_z_reducer( min_z );
+        Kokkos::parallel_reduce(
+            "CabanaPD::Particles::min_max_positions",
+            Kokkos::RangePolicy<ExecSpace>( exec_space, 0, n_local ),
+            KOKKOS_LAMBDA( const int i, double& min_x, double& min_y,
+                           double& min_z, double& max_x, double& max_y,
+                           double& max_z ) {
+                if ( x( i, 0 ) > max_x )
+                    max_x = x( i, 0 );
+                else if ( x( i, 0 ) < min_x )
+                    min_x = x( i, 0 );
+                if ( x( i, 1 ) > max_y )
+                    max_y = x( i, 1 );
+                else if ( x( i, 1 ) < min_y )
+                    min_y = x( i, 1 );
+                if ( x( i, 2 ) > max_z )
+                    max_z = x( i, 2 );
+                else if ( x( i, 2 ) < min_z )
+                    min_z = x( i, 2 );
+            },
+            min_x_reducer, min_y_reducer, min_z_reducer, max_x_reducer,
+            max_y_reducer, max_z_reducer );
+
+        std::array<double, 3> min_corner = { min_x, min_y, min_z };
+        std::array<double, 3> max_corner = { max_x, max_y, max_z };
+
+        std::array<int, 3> num_cells;
+        for ( int d = 0; d < 3; d++ )
+        {
+            num_cells[d] =
+                static_cast<int>( ( max_corner[d] - min_corner[d] ) / dx );
+            if ( num_cells[d] == 0 )
+            {
+                num_cells[d]++;
+                // Issues with pseudo-2d with only one bin in a given direction.
+                max_corner[d] += dx * 3;
+                min_corner[d] -= dx * 3;
+            }
+        }
+        createDomain( min_corner, max_corner, num_cells );
+    }
+
+    template <class ExecSpace, std::size_t NSD = dim>
+    std::enable_if_t<2 == NSD, void>
+    updateAfterRead( const ExecSpace& exec_space, double dx )
+    {
+        // Because this may be a non-uniform mesh, build a background mesh with
+        // dx=cutoff and force halo_width=1
+        halo_width = 1;
+        auto x = sliceReferencePosition();
+        n_local = x.size();
+        n_ghost = 0;
+        size = n_local;
+        update_global();
+
+        double max_x;
+        Kokkos::Max<double> max_x_reducer( max_x );
+        double min_x;
+        Kokkos::Min<double> min_x_reducer( min_x );
+        double max_y;
+        Kokkos::Max<double> max_y_reducer( max_y );
+        double min_y;
+        Kokkos::Min<double> min_y_reducer( min_y );
+        Kokkos::parallel_reduce(
+            "CabanaPD::Particles::min_max_positions",
+            Kokkos::RangePolicy<ExecSpace>( exec_space, 0, n_local ),
+            KOKKOS_LAMBDA( const int i, double& min_x, double& min_y,
+                           double& max_x, double& max_y ) {
+                if ( x( i, 0 ) > max_x )
+                    max_x = x( i, 0 );
+                else if ( x( i, 0 ) < min_x )
+                    min_x = x( i, 0 );
+                if ( x( i, 1 ) > max_y )
+                    max_y = x( i, 1 );
+                else if ( x( i, 1 ) < min_y )
+                    min_y = x( i, 1 );
+            },
+            min_x_reducer, min_y_reducer, max_x_reducer, max_y_reducer );
+
+        std::array<double, 2> min_corner = { min_x, min_y };
+        std::array<double, 2> max_corner = { max_x, max_y };
+
+        std::array<int, 2> num_cells;
+        for ( int d = 0; d < 2; d++ )
+        {
+            num_cells[d] =
+                static_cast<int>( ( max_corner[d] - min_corner[d] ) / dx );
+            if ( num_cells[d] == 0 )
+            {
+                num_cells[d]++;
+                // Potentially issues with pseudo-1d with only one bin in a
+                // given direction.
+                max_corner[d] += dx * 3;
+                min_corner[d] -= dx * 3;
+            }
+            createDomain( min_corner, max_corner, num_cells );
+        }
+    }
+
+    void update_global()
+    {
         // Not using Allreduce because global count is only used for printing.
         MPI_Reduce( &n_local, &n_global, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0,
                     MPI_COMM_WORLD );
@@ -519,6 +649,16 @@ class Particles<MemorySpace, LPS, Dimension>
         _aosoa_theta.resize( 0 );
     }
 
+    template <class ExecSpace>
+    void updateAfterRead( const ExecSpace& exec_space, const int hw, double dx )
+    {
+        base_type::updateAfterRead( exec_space, hw, dx );
+
+        // Only need to resize LPS variables
+        _aosoa_theta.resize( n_local + n_ghost );
+        _aosoa_m.resize( n_local + n_ghost );
+    }
+
     auto sliceDilatation()
     {
         return Cabana::slice<0>( _aosoa_theta, "dilatation" );
diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp
index 07407f86..800e78d3 100644
--- a/src/CabanaPD_Solver.hpp
+++ b/src/CabanaPD_Solver.hpp
@@ -357,9 +357,11 @@ class SolverFracture
     SolverFracture( input_type _inputs,
                     std::shared_ptr<particle_type> _particles,
                     force_model_type force_model, bc_type bc,
-                    prenotch_type prenotch )
+                    prenotch_type prenotch,
+                    const bool force_boundary_condition )
         : base_type( _inputs, _particles, force_model )
         , boundary_condition( bc )
+        , force_bc( force_boundary_condition )
     {
         init_timer.reset();
 
@@ -391,8 +393,18 @@ class SolverFracture
         computeEnergy( *force, *particles, *neighbors, mu, neigh_iter_tag() );
 
         // Add boundary condition.
-        boundary_condition.apply( exec_space(), *particles );
-
+        auto x = particles->sliceReferencePosition();
+        // FIXME: this needs to be generalized to any field.
+        if ( force_bc )
+        {
+            auto f = particles->sliceForce();
+            boundary_condition.apply( exec_space(), f, x, 0.0 );
+        }
+        else
+        {
+            auto u = particles->sliceDisplacement();
+            boundary_condition.apply( exec_space(), u, x, 0.0 );
+        }
         particles->output( 0, 0.0, output_reference );
         init_time += init_timer.seconds();
     }
@@ -436,7 +448,19 @@ class SolverFracture
             force_time += force_timer.seconds();
 
             // Add boundary condition.
-            boundary_condition.apply( exec_space{}, *particles );
+            auto x = particles->sliceReferencePosition();
+            // FIXME: this needs to be generalized to any field.
+            auto time = step * inputs->timestep;
+            if ( force_bc )
+            {
+                auto f = particles->sliceForce();
+                boundary_condition.apply( exec_space(), f, x, time );
+            }
+            else
+            {
+                auto u = particles->sliceDisplacement();
+                boundary_condition.apply( exec_space(), u, x, time );
+            }
 
             // Integrate - velocity Verlet second half.
             integrate_timer.reset();
@@ -473,6 +497,7 @@ class SolverFracture
     using base_type::neighbors;
     using base_type::particles;
     bc_type boundary_condition;
+    bool force_bc;
 
     using NeighborView = typename Kokkos::View<int**, memory_space>;
     NeighborView mu;
@@ -510,12 +535,13 @@ template <class MemorySpace, class InputsType, class ParticleType,
           class ForceModel, class BCType, class PrenotchType>
 auto createSolverFracture( InputsType inputs,
                            std::shared_ptr<ParticleType> particles,
-                           ForceModel model, BCType bc, PrenotchType prenotch )
+                           ForceModel model, BCType bc, PrenotchType prenotch,
+                           const bool force_bc = true )
 {
     return std::make_shared<
         SolverFracture<MemorySpace, InputsType, ParticleType, ForceModel,
                        BCType, PrenotchType>>( inputs, particles, model, bc,
-                                               prenotch );
+                                               prenotch, force_bc );
 }
 
 /*