Initial version of multiple pre-notches

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 <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+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<int, 3> 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<double, 3> low_corner = { low_x, low_y, low_z };
+        std::array<double, 3> 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<double, 3> p01 = { low_x, y_prenotch1, low_z };
+        // 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 );
+
+        // 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<double, 3> p01 = { Xc1, Yc1, low_z };
+            Kokkos::Array<double, 3> v1 = { l * cos( theta ), l * sin( theta ),
+                                            0 };
+            Kokkos::Array<double, 3> 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<Kokkos::Array<double, 3>, 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<CabanaPD::PMB, CabanaPD::Fracture>;
+        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<memory_space, typename model_type::base_model>>(
+            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<CabanaPD::RegionBoundary> 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<memory_space>(
+            inputs, particles, force_model, bc, prenotch );
+        cabana_pd->init_force();
+        cabana_pd->run();
+    }
+
+    MPI_Finalize();
+}