Enable pre-notches with differing orientations

README.md

@@ -151,6 +151,12 @@ The third example is crack branching in a pre-notched soda-lime glass plate due
 ./CabanaPD/build/install/bin/CrackBranching CabanaPD/examples/mechanics/inputs/crack_branching.json
 ```
 
+A similar case but with multiple random pre-notches can be run with:
+
+```
+./CabanaPD/build/install/bin/RandomCracks CabanaPD/examples/mechanics/inputs/random_cracks.json
+```
+
 The fourth example is a fragmenting cylinder due to internal pressure [4]. The example can be run with:
 
 ```

examples/mechanics/CMakeLists.txt

@@ -10,4 +10,7 @@ target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD)
 add_executable(FragmentingCylinder fragmenting_cylinder.cpp)
 target_link_libraries(FragmentingCylinder LINK_PUBLIC CabanaPD)
 
-install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DESTINATION ${CMAKE_INSTALL_BINDIR})
+add_executable(RandomCracks random_cracks.cpp)
+target_link_libraries(RandomCracks LINK_PUBLIC CabanaPD)
+
+install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder RandomCracks DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/mechanics/inputs/random_cracks.json

@@ -0,0 +1,16 @@
+{
+    "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"},
+    "minimum_prenotch_length" : {"value": 0.001, "unit": "m"},
+    "maximum_prenotch_length" : {"value": 0.0025, "unit": "m"},
+    "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/random_cracks.cpp

@@ -0,0 +1,205 @@
+/****************************************************************************
+ * 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 propagation from multiple pre-notches.
+void randomCracksExample( const std::string filename )
+{
+    // ====================================================
+    //             Use default 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 delta = inputs["horizon"];
+    delta += 1e-10;
+
+    // ====================================================
+    //                  Discretization
+    // ====================================================
+    std::array<double, 3> low_corner = inputs["low_corner"];
+    std::array<double, 3> high_corner = inputs["high_corner"];
+
+    // ====================================================
+    //                    Pre-notches
+    // ====================================================
+    // Number of pre-notches
+    constexpr int Npn = 200;
+
+    double minl = inputs["minimum_prenotch_length"];
+    double maxl = inputs["maximum_prenotch_length"];
+    double thickness = high_corner[2] - low_corner[2];
+
+    // Initialize pre-notch arrays
+    Kokkos::Array<Kokkos::Array<double, 3>, Npn> notch_positions;
+    Kokkos::Array<Kokkos::Array<double, 3>, Npn> notch_v1;
+    Kokkos::Array<Kokkos::Array<double, 3>, Npn> notch_v2;
+
+    // Changing this seed will re-randomize the cracks.
+    std::size_t seed = 44758454;
+    // Random number generator
+    std::mt19937 gen( seed );
+    std::uniform_real_distribution<> dis( 0.0, 1.0 );
+
+    // Loop over pre-notches
+    for ( int n = 0; n < Npn; n++ )
+    {
+        // Random numbers for pre-notch position
+        double random_number_x = dis( gen );
+        double random_number_y = dis( gen );
+
+        // Random coordinates of one endpoint of the pre-notch
+        // Note: the addition and subtraction of "maxl" ensures the prenotch
+        // does not extend outside the domain
+        double Xc1 = ( low_corner[0] + maxl ) +
+                     ( ( high_corner[0] - maxl ) - ( low_corner[0] + maxl ) ) *
+                         random_number_x;
+        double Yc1 = ( low_corner[1] + maxl ) +
+                     ( ( high_corner[1] - maxl ) - ( low_corner[1] + maxl ) ) *
+                         random_number_y;
+        Kokkos::Array<double, 3> p0 = { Xc1, Yc1, low_corner[2] };
+
+        // Assign pre-notch position
+        notch_positions[n] = p0;
+
+        // Random pre-notch length on XY-plane
+        double random_number_l = dis( gen );
+        double l = minl + ( maxl - minl ) * random_number_l;
+
+        // Random pre-notch orientation on XY-plane
+        double random_number_theta = dis( gen );
+        double theta = CabanaPD::pi * random_number_theta;
+
+        // Pre-notch v1 vector on XY-plane
+        Kokkos::Array<double, 3> v1 = { l * cos( theta ), l * sin( theta ), 0 };
+        notch_v1[n] = v1;
+
+        // Random number for y-component of v2 vector: the angle of v2 in the
+        // YZ-plane is between -45 and 45 deg.
+        double random_number_v2_y = dis( gen );
+
+        // Pre-notch v2 vector on YZ-plane
+        Kokkos::Array<double, 3> v2 = {
+            0, ( -1.0 + 2.0 * random_number_v2_y ) * thickness, thickness };
+        notch_v2[n] = v2;
+    }
+
+    CabanaPD::Prenotch<Npn> prenotch( notch_v1, notch_v2, notch_positions );
+
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    using model_type = CabanaPD::ForceModel<CabanaPD::PMB>;
+    model_type force_model( delta, K, G0 );
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    // Note that individual inputs can be passed instead (see other examples).
+    auto particles = CabanaPD::createParticles<memory_space, model_type>(
+        exec_space{}, inputs );
+
+    // ====================================================
+    //                Boundary conditions planes
+    // ====================================================
+    double dy = particles->dx[1];
+    CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane1(
+        low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy,
+        low_corner[2], high_corner[2] );
+    CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane2(
+        low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy,
+        low_corner[2], high_corner[2] );
+    std::vector<CabanaPD::RegionBoundary<CabanaPD::RectangularPrism>> planes = {
+        plane1, plane2 };
+
+    // ====================================================
+    //            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_y + delta + 1e-10 ||
+             x( pid, 1 ) >= plane2.high_y - delta - 1e-10 )
+            nofail( pid ) = 1;
+    };
+    particles->updateParticles( exec_space{}, init_functor );
+
+    // ====================================================
+    //                   Create solver
+    // ====================================================
+    auto cabana_pd =
+        CabanaPD::createSolver<memory_space>( inputs, particles, force_model );
+
+    // ====================================================
+    //                Boundary conditions
+    // ====================================================
+    // Create BC last to ensure ghost particles are included.
+    double sigma0 = inputs["traction"];
+    double b0 = sigma0 / dy;
+    f = particles->sliceForce();
+    x = 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;
+    };
+    auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes,
+                                       true );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+    cabana_pd->init( bc, prenotch );
+    cabana_pd->run( bc );
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    randomCracksExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

src/CabanaPD_Prenotch.hpp

@@ -171,24 +171,45 @@ int bondPrenotchIntersection( const Kokkos::Array<double, 3> v1,
     return keep_bond;
 }
 
-template <std::size_t NumNotch>
+template <std::size_t NumNotch, std::size_t NumVector = 1>
 struct Prenotch
 {
     static constexpr std::size_t num_notch = NumNotch;
-    Kokkos::Array<double, 3> _v1;
-    Kokkos::Array<double, 3> _v2;
-    Kokkos::Array<Kokkos::Array<double, 3>, num_notch> _p0_list;
+    static constexpr std::size_t num_vector = NumNotch;
+    Kokkos::Array<Kokkos::Array<double, 3>, num_vector> _v1;
+    Kokkos::Array<Kokkos::Array<double, 3>, num_vector> _v2;
+    Kokkos::Array<Kokkos::Array<double, 3>, num_notch> _p0;
+    bool fixed_orientation;
 
     Timer _timer;
 
+    // Default constructor
     Prenotch() {}
 
+    // Constructor if all pre-notches are oriented the same way (e.g.
+    // Kalthoff-Winkler).
     Prenotch( Kokkos::Array<double, 3> v1, Kokkos::Array<double, 3> v2,
-              Kokkos::Array<Kokkos::Array<double, 3>, num_notch> p0_list )
+              Kokkos::Array<Kokkos::Array<double, 3>, num_notch> p0 )
+        : _v1( { v1 } )
+        , _v2( { v2 } )
+        , _p0( p0 )
+    {
+        fixed_orientation = true;
+    }
+
+    // Constructor for general case of any orientation for any number of
+    // pre-notches.
+    Prenotch( Kokkos::Array<Kokkos::Array<double, 3>, num_vector> v1,
+              Kokkos::Array<Kokkos::Array<double, 3>, num_vector> v2,
+              Kokkos::Array<Kokkos::Array<double, 3>, num_notch> p0 )
         : _v1( v1 )
         , _v2( v2 )
-        , _p0_list( p0_list )
+        , _p0( p0 )
     {
+        static_assert(
+            num_vector == num_notch,
+            "Number of orientation vectors must match number of pre-notches." );
+        fixed_orientation = false;
     }
 
     template <class ExecSpace, class NeighborView, class Particles,
@@ -202,11 +223,14 @@ struct Prenotch
         // TODO: decide whether to disallow prenotches in frozen particles.
         Kokkos::RangePolicy<ExecSpace> policy( 0, particles.localOffset() );
 
-        auto v1 = _v1;
-        auto v2 = _v2;
-        for ( std::size_t p = 0; p < _p0_list.size(); p++ )
+        for ( std::size_t p = 0; p < _p0.size(); p++ )
         {
-            auto p0 = _p0_list[p];
+            // These will always be different positions.
+            auto p0 = _p0[p];
+            // These may all have the same orientation or all different
+            // orientations.
+            auto v1 = getV1( p );
+            auto v2 = getV2( p );
             auto notch_functor = KOKKOS_LAMBDA( const int i )
             {
                 std::size_t num_neighbors =
@@ -235,6 +259,22 @@ struct Prenotch
         _timer.stop();
     }
     auto time() { return _timer.time(); };
+
+    auto getV1( const int p )
+    {
+        if ( fixed_orientation )
+            return _v1[0];
+        else
+            return _v1[p];
+    }
+
+    auto getV2( const int p )
+    {
+        if ( fixed_orientation )
+            return _v2[0];
+        else
+            return _v2[p];
+    }
 };
 
 } // namespace CabanaPD