diff --git a/examples/multiple_cracks.cpp b/examples/multiple_cracks.cpp new file mode 100644 index 00000000..d3851b27 --- /dev/null +++ b/examples/multiple_cracks.cpp @@ -0,0 +1,189 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 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 +#include + +#include "mpi.h" + +#include + +#include + +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; + + // Plate dimensions (m) + double height = 0.1; + double width = 0.04; + double thickness = 0.002; + + // Domain + std::array num_cell = { 400, 160, 8 }; // 400 x 160 x 8 + double low_x = -0.5 * height; + double low_y = -0.5 * width; + double low_z = -0.5 * thickness; + double high_x = 0.5 * height; + double high_y = 0.5 * width; + double high_z = 0.5 * thickness; + std::array low_corner = { low_x, low_y, low_z }; + std::array high_corner = { high_x, high_y, high_z }; + + // Time + double t_final = 43e-6; + double dt = 5e-8; + double output_frequency = 5; + bool output_reference = true; + + // Material constants + double E = 72e+9; // [Pa] + double nu = 0.25; // unitless + double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa] + double rho0 = 2440; // [kg/m^3] + double G0 = 3.8; // [J/m^2] + + // PD horizon + double delta = 0.001 + 1e-10; + + // FIXME: set halo width based on delta + int m = std::floor( + delta / ( ( high_corner[0] - low_corner[0] ) / num_cell[0] ) ); + int halo_width = m + 1; + + // Prenotch + // double L_prenotch = height / 2.0; + // double y_prenotch1 = 0.0; + // Kokkos::Array p01 = { low_x, y_prenotch1, low_z }; + // Kokkos::Array v1 = { L_prenotch, 0, 0 }; + // Kokkos::Array v2 = { 0, 0, thickness }; + // Kokkos::Array, 1> notch_positions = { p01 }; + // CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions ); + + // Multiple random pre-notches + + // Number of pre-notches + int Npn = 15; + + // Minimum possible pre-notch length + minl = 0.1; + + // Maximim possible pre-notch length + maxl = 0.25; + + // Initialize pre-notch endpoint + double Xc1 = 0.0; + double Yc1 = 0.0; + + // Define Pi constant + double Pi = 3.141592653589793; + + // Initialize pre-notch length and orientation angle + double l = 0.0; + double theta = 0.0; + + // Initialize random number generator + srand( (unsigned)time( &t ) ); + + for ( int n = 0; n < Npn; n++ ) + { + // Coordinates of one endpoint of the pre-notch (random) + // Note: the addition and subtraction of "maxl" ensures + // the prenotch does not extend outside the domain + Xc1 = ( low_x + maxl ) + + ( ( high_x - maxl ) - ( low_x + maxl ) ) * rand(); + Yc1 = ( low_y + maxl ) + + ( ( high_y - maxl ) - ( low_y + maxl ) ) * rand(); + + // Random pre-notch length + l = minl + ( maxl - minl ) * rand(); + + // Random orientation + theta = ( 2 * Pi ) * rand(); + + Kokkos::Array p01 = { Xc1, Yc1, low_z }; + Kokkos::Array v1 = { l * cos( theta ), l * sin( theta ), + 0 }; + Kokkos::Array v2 = { 0, 0, thickness }; + + CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions ); + + // This line needs to be extended to accummulate notches within the + // loop + Kokkos::Array, 2> notch_positions = { + p01, p02 }; + } + + // This line needs to be extended to have a list of v1, so that each + // pre-notch has a different orientation + CabanaPD::Prenotch<2> prenotch( v1, v2, notch_positions ); + + // Choose force model type. + using model_type = + CabanaPD::ForceModel; + model_type force_model( delta, K, G0 ); + CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final, + dt, output_frequency, output_reference ); + inputs.read_args( argc, argv ); + + // Create particles from mesh. + // Does not set displacements, velocities, etc. + auto particles = std::make_shared< + CabanaPD::Particles>( + exec_space(), inputs.low_corner, inputs.high_corner, + inputs.num_cells, halo_width ); + + // Define particle initialization. + auto x = particles->sliceReferencePosition(); + auto v = particles->sliceVelocity(); + auto f = particles->sliceForce(); + auto rho = particles->sliceDensity(); + auto nofail = particles->sliceNoFail(); + + // Relying on uniform grid here. + double dy = particles->dx[1]; + double b0 = 2e6 / dy; // Pa/m + + CabanaPD::RegionBoundary plane1( low_x, high_x, low_y - dy, low_y + dy, + low_z, high_z ); + CabanaPD::RegionBoundary plane2( low_x, high_x, high_y - dy, + high_y + dy, low_z, high_z ); + std::vector planes = { plane1, plane2 }; + auto bc = + createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, + exec_space{}, *particles, planes, b0 ); + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + rho( pid ) = rho0; + // Set the 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 ); + + // FIXME: use createSolver to switch backend at runtime. + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model, bc, prenotch ); + cabana_pd->init_force(); + cabana_pd->run(); + } + + MPI_Finalize(); +}