Draft - Initial indentation example

examples/mechanics/CMakeLists.txt

@@ -22,4 +22,13 @@ target_link_libraries(PlateWithHole LINK_PUBLIC CabanaPD)
 add_executable(CrackInclusion crack_inclusion.cpp)
 target_link_libraries(CrackInclusion LINK_PUBLIC CabanaPD)
 
-install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder RandomCracks DogboneTensileTest PlateWithHole CrackInclusion DESTINATION ${CMAKE_INSTALL_BINDIR})
+add_executable(FlatIndenter flat_indenter.cpp)
+target_link_libraries(FlatIndenter LINK_PUBLIC CabanaPD)
+
+add_executable(SphericalIndenter spherical_indenter.cpp)
+target_link_libraries(SphericalIndenter LINK_PUBLIC CabanaPD)
+
+add_executable(ConicalIndenter conical_indenter.cpp)
+target_link_libraries(ConicalIndenter LINK_PUBLIC CabanaPD)
+
+install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder RandomCracks DogboneTensileTest PlateWithHole CrackInclusion FlatIndenter SphericalIndenter ConicalIndenter DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/mechanics/flat_indenter.cpp

@@ -0,0 +1,164 @@
+/****************************************************************************
+ * Copyright (c) 2022 by Oak Ridge National Laboratory                      *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of CabanaPD. CabanaPD is distributed under a           *
+ * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
+ * the top-level directory.                                                 *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+
+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+// Simulate crack branching from an pre-crack.
+void flatIndenterExample( const std::string filename )
+{
+    // ====================================================
+    //               Choose Kokkos spaces
+    // ====================================================
+    using exec_space = Kokkos::DefaultExecutionSpace;
+    using memory_space = typename exec_space::memory_space;
+
+    // ====================================================
+    //                   Read inputs
+    // ====================================================
+    CabanaPD::Inputs inputs( filename );
+
+    // ====================================================
+    //                Material parameters
+    // ====================================================
+    double rho0 = inputs["density"];
+    double E = inputs["elastic_modulus"];
+    double nu = 0.25; // Use bond-based model
+    double K = E / ( 3 * ( 1 - 2 * nu ) );
+    double G0 = inputs["fracture_energy"];
+    double horizon = inputs["horizon"];
+    horizon += 1e-10;
+
+    // ====================================================
+    //                  Discretization
+    // ====================================================
+    std::array<double, 3> low_corner = inputs["low_corner"];
+    std::array<double, 3> high_corner = inputs["high_corner"];
+
+    // ====================================================
+    //                    Pre-notch
+    // ====================================================
+/*	
+    double height = inputs["system_size"][0];
+    double thickness = inputs["system_size"][2];
+    double L_prenotch = height / 2.0;
+    double y_prenotch = 0.0;
+    Kokkos::Array<double, 3> p01 = { low_corner[0], y_prenotch, low_corner[2] };
+    Kokkos::Array<double, 3> v1 = { L_prenotch, 0, 0 };
+    Kokkos::Array<double, 3> v2 = { 0, 0, thickness };
+    Kokkos::Array<Kokkos::Array<double, 3>, 1> notch_positions = { p01 };
+    CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions );
+*/
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    using model_type = CabanaPD::PMB;
+    CabanaPD::ForceModel force_model( model_type{}, horizon, K, G0 );
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    CabanaPD::Particles particles( memory_space{}, model_type{} );
+
+    // Note that individual inputs can be passed instead (see other examples).
+    particles.domain( inputs );
+    particles.create( exec_space{} );
+
+    // ====================================================
+    //                Boundary conditions planes
+    // ====================================================
+
+    double dz = particles.dx[2];
+    CabanaPD::Region<CabanaPD::RectangularPrism> square_pressure(
+        0.5 * low_corner[0], 0.5 * high_corner[0], 0.5 * low_corner[1],
+        0.5 * high_corner[1], high_corner[2] - dz, high_corner[2] + dz );
+
+    // ====================================================
+    //            Custom particle initialization
+    // ====================================================
+    auto rho = particles.sliceDensity();
+    auto x = particles.sliceReferencePosition();
+    auto v = particles.sliceVelocity();
+    auto f = particles.sliceForce();
+    auto nofail = particles.sliceNoFail();
+
+
+    auto init_functor = KOKKOS_LAMBDA( const int pid )
+    {
+
+        // Density
+        rho( pid ) = rho0;
+/*		
+        // No-fail zone
+        if ( x( pid, 1 ) <= plane1.low[1] + horizon + 1e-10 ||
+             x( pid, 1 ) >= plane2.high[1] - horizon - 1e-10 )
+            nofail( pid ) = 1;
+*/
+    };
+    particles.update( exec_space{}, init_functor );
+
+    // ====================================================
+    //                   Create solver
+    // ====================================================
+    CabanaPD::Solver solver( inputs, particles, force_model );
+
+    // ====================================================
+    //                Boundary conditions
+    // ====================================================
+    // Create BC last to ensure ghost particles are included.
+    double sigma0 = inputs["traction"];
+    double b0 = -sigma0 / dz;
+	auto indent_force = 0;
+    f = solver.particles.sliceForce();
+    x = solver.particles.sliceReferencePosition();
+    // Create a symmetric force BC in the z-direction.
+
+    auto bc_op = KOKKOS_LAMBDA( const int pid, const double )
+    {		
+	  double xsq = x(pid,0) * x(pid,0);
+	  double ysq = x(pid,1) * x(pid,1);
+	  double rsq = xsq + ysq;
+
+	        if ( xsq + ysq < 9e-6 )
+	        {
+	          f( pid, 2 ) += b0;
+	        }  	
+
+
+    };
+    auto bc = createBoundaryCondition( bc_op, exec_space{}, solver.particles,
+                                       true, square_pressure );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+    //solver.init( bc, prenotch );
+	solver.init( bc );
+    solver.run( bc );
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    flatIndenterExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

examples/mechanics/indentation_test.cpp

@@ -0,0 +1,164 @@
+/****************************************************************************
+ * Copyright (c) 2022 by Oak Ridge National Laboratory                      *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of CabanaPD. CabanaPD is distributed under a           *
+ * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
+ * the top-level directory.                                                 *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+
+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+// Simulate crack branching from an pre-crack.
+void crackBranchingExample( const std::string filename )
+{
+    // ====================================================
+    //               Choose Kokkos spaces
+    // ====================================================
+    using exec_space = Kokkos::DefaultExecutionSpace;
+    using memory_space = typename exec_space::memory_space;
+
+    // ====================================================
+    //                   Read inputs
+    // ====================================================
+    CabanaPD::Inputs inputs( filename );
+
+    // ====================================================
+    //                Material parameters
+    // ====================================================
+    double rho0 = inputs["density"];
+    double E = inputs["elastic_modulus"];
+    double nu = 0.25; // Use bond-based model
+    double K = E / ( 3 * ( 1 - 2 * nu ) );
+    double G0 = inputs["fracture_energy"];
+    double horizon = inputs["horizon"];
+    horizon += 1e-10;
+
+    // ====================================================
+    //                  Discretization
+    // ====================================================
+    std::array<double, 3> low_corner = inputs["low_corner"];
+    std::array<double, 3> high_corner = inputs["high_corner"];
+
+    // ====================================================
+    //                    Pre-notch
+    // ====================================================
+    /*
+    double height = inputs["system_size"][0];
+    double thickness = inputs["system_size"][2];
+    double L_prenotch = height / 2.0;
+    double y_prenotch = 0.0;
+    Kokkos::Array<double, 3> p01 = { low_corner[0], y_prenotch, low_corner[2] };
+    Kokkos::Array<double, 3> v1 = { L_prenotch, 0, 0 };
+    Kokkos::Array<double, 3> v2 = { 0, 0, thickness };
+    Kokkos::Array<Kokkos::Array<double, 3>, 1> notch_positions = { p01 };
+    CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions );
+    */
+
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    using model_type = CabanaPD::PMB;
+    CabanaPD::ForceModel force_model( model_type{}, horizon, K, G0 );
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    CabanaPD::Particles particles( memory_space{}, model_type{} );
+
+    // Note that individual inputs can be passed instead (see other examples).
+    particles.domain( inputs );
+    particles.create( exec_space{} );
+
+    // ====================================================
+    //                Boundary conditions planes
+    // ====================================================
+    // double dy = particles.dx[1];
+    // CabanaPD::Region<CabanaPD::RectangularPrism> plane1(
+    //    low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy,
+    //    low_corner[2], high_corner[2] );
+    // CabanaPD::Region<CabanaPD::RectangularPrism> plane2(
+    //    low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] +
+    //    dy, low_corner[2], high_corner[2] );
+
+    double dz = particles.dx[2];
+    CabanaPD::Region<CabanaPD::RectangularPrism> square_pressure(
+        0.5 * low_corner[0], 0.5 * high_corner[0], 0.5 * low_corner[1],
+        0.5 * high_corner[1], high_corner[2] - dz, high_corner[2] + dz );
+
+    // ====================================================
+    //            Custom particle initialization
+    // ====================================================
+    auto rho = particles.sliceDensity();
+    auto x = particles.sliceReferencePosition();
+    auto v = particles.sliceVelocity();
+    auto f = particles.sliceForce();
+    auto nofail = particles.sliceNoFail();
+
+    auto init_functor = KOKKOS_LAMBDA( const int pid )
+    {
+        // Density
+        rho( pid ) = rho0;
+        // No-fail zone
+        // if ( x( pid, 1 ) <= plane1.low[1] + horizon + 1e-10 ||
+        //     x( pid, 1 ) >= plane2.high[1] - horizon - 1e-10 )
+        //    nofail( pid ) = 1;
+    };
+    particles.update( exec_space{}, init_functor );
+
+    // ====================================================
+    //                   Create solver
+    // ====================================================
+    CabanaPD::Solver solver( inputs, particles, force_model );
+
+    // ====================================================
+    //                Boundary conditions
+    // ====================================================
+    // Create BC last to ensure ghost particles are included.
+    double sigma0 = inputs["traction"];
+    // double b0 = sigma0 / dy;
+    double b0 = -sigma0 / dz;
+    f = solver.particles.sliceForce();
+    x = solver.particles.sliceReferencePosition();
+    // Create a symmetric force BC in the y-direction.
+    auto bc_op = KOKKOS_LAMBDA( const int pid, const double )
+    {
+        // auto ypos = x( pid, 1 );
+        // auto sign = std::abs( ypos ) / ypos;
+        // f( pid, 1 ) += b0 * sign;
+        f( pid, 2 ) += b0;
+    };
+    // auto bc = createBoundaryCondition( bc_op, exec_space{}, solver.particles,
+    //                                    true, plane1, plane2 );
+
+    auto bc = createBoundaryCondition( bc_op, exec_space{}, solver.particles,
+                                       true, square_pressure );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+    // solver.init( bc, prenotch );
+    solver.init( bc );
+    solver.run( bc );
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    crackBranchingExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

examples/mechanics/inputs/flat_indenter.json

@@ -0,0 +1,14 @@
+{
+    "num_cells"              : {"value": [400, 160, 8]},
+    "system_size"            : {"value": [0.1, 0.04, 0.002], "unit": "m"},
+    "density"                : {"value": 2440,  "unit": "kg/m^3"},
+    "elastic_modulus"        : {"value": 72e+9, "unit": "Pa"},
+    "fracture_energy"        : {"value": 3.8,   "unit": "J/m^2"},
+    "horizon"                : {"value": 0.001, "unit": "m"},
+    "traction"               : {"value": 2e6,   "unit": "Pa"},
+    "final_time"             : {"value": 43e-6, "unit": "s"},
+    "timestep"               : {"value": 4.5e-8,  "unit": "s"},
+    "timestep_safety_factor" : {"value": 0.85},
+    "output_frequency"       : {"value": 5},
+    "output_reference"       : {"value": true}
+}

examples/mechanics/inputs/indentation_test.json

@@ -0,0 +1,19 @@
+{
+    "num_cells_OLD"              : {"value": [400, 400, 80]},
+    "num_cells"              : {"value": [200, 200, 40]},
+    "system_size_OLD"            : {"value": [0.1, 0.04, 0.002], "unit": "m"},
+    "system_size"            : {"value": [0.1, 0.1, 0.02], "unit": "m"},
+    "density"                : {"value": 2440,  "unit": "kg/m^3"},
+    "elastic_modulus"        : {"value": 72e+9, "unit": "Pa"},
+    "fracture_energy"        : {"value": 3.8,   "unit": "J/m^2"},
+    "horizon_ref"                : {"value": 0.001, "unit": "m"},  
+    "horizon"                : {"value": 0.0015, "unit": "m"},  
+    "traction"               : {"value": 2e6,   "unit": "Pa"},
+    "final_time"             : {"value": 43e-6, "unit": "s"},
+    "timestep_OLD"               : {"value": 4.5e-8,  "unit": "s"},
+    "timestep"               : {"value": 6.75e-8,  "unit": "s"},
+    "timestep_safety_factor" : {"value": 0.85},
+    "output_frequency"       : {"value": 5},
+    "output_reference_OLD"       : {"value": true},
+    "output_reference"       : {"value": false}
+}