Skip to content

Commit

Permalink
Rewrite boundary conditions for multiple geometries
Browse files Browse the repository at this point in the history
  • Loading branch information
@streeve
streeve committed Aug 8, 2024
1 parent 95bc551 commit b50761c
Show file tree
Hide file tree
Showing 3 changed files with 80 additions and 31 deletions.
15 changes: 8 additions & 7 deletions examples/mechanics/crack_branching.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,13 +90,14 @@ void crackBranchingExample( const std::string filename )
double dy = particles->dx[1];
double b0 = sigma0 / dy;

CabanaPD::RegionBoundary plane1( low_corner[0], high_corner[0],
low_corner[1] - dy, low_corner[1] + dy,
low_corner[2], high_corner[2] );
CabanaPD::RegionBoundary 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> planes = { plane1, plane2 };
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 };
auto particles_f = particles->getForce();
auto particles_x = particles->getReferencePosition();
// Create a symmetric force BC in the y-direction.
Expand Down
5 changes: 3 additions & 2 deletions examples/mechanics/kalthoff_winkler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,10 +98,11 @@ void kalthoffWinklerExample( const std::string filename )
// ====================================================
double dx = particles->dx[0];
double x_bc = -0.5 * height;
CabanaPD::RegionBoundary plane(
CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane(
x_bc - dx, x_bc + dx * 1.25, y_prenotch1 - dx * 0.25,
y_prenotch2 + dx * 0.25, -thickness, thickness );
std::vector<CabanaPD::RegionBoundary> planes = { plane };
std::vector<CabanaPD::RegionBoundary<CabanaPD::Rectilinear>> planes = {
plane };
auto bc = createBoundaryCondition( CabanaPD::ForceValueBCTag{}, 0.0,
exec_space{}, *particles, planes );

Expand Down
91 changes: 69 additions & 22 deletions src/CabanaPD_Boundary.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,20 @@
namespace CabanaPD
{

// Define a plane or other rectilinear subset of the system as the boundary.
struct RegionBoundary
struct RectangularPrism
{
};
struct Cylinder
{
};

// Forward declaration.
template <class GeometryType>
struct RegionBoundary;

// Define a subset of the system as the boundary with a rectangular prism.
template <>
struct RegionBoundary<RectangularPrism>
{
double low_x;
double high_x;
Expand All @@ -40,14 +52,53 @@ struct RegionBoundary
, high_y( _high_y )
, low_z( _low_z )
, high_z( _high_z ){};

template <class PositionType>
KOKKOS_INLINE_FUNCTION bool inside( const PositionType& x,
const int pid ) const
{
return ( x( pid, 0 ) >= low_x && x( pid, 0 ) <= high_x &&
x( pid, 1 ) >= low_y && x( pid, 1 ) <= high_y &&
x( pid, 2 ) >= low_z && x( pid, 2 ) <= high_z );
}
};

// Define a subset of the system as the boundary with a cylinder.
template <>
struct RegionBoundary<Cylinder>
{
double radius_in;
double radius_out;
double low_z;
double high_z;
double x_center = 0.0;
double y_center = 0.0;

RegionBoundary( const double _radius_in, const double _radius_out,
const double _low_z, const double _high_z )
: radius_in( _radius_in )
, radius_out( _radius_out )
, low_z( _low_z )
, high_z( _high_z ){};

template <class PositionType>
KOKKOS_INLINE_FUNCTION bool inside( const PositionType& x,
const int pid ) const
{
double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) +
( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center );
return ( rsq >= radius_in * radius_in &&
rsq <= radius_out * radius_out && x( pid, 2 ) >= low_z &&
x( pid, 2 ) <= high_z );
}
};

template <class MemorySpace, class BoundaryType>
struct BoundaryIndexSpace;

// FIXME: fails for some cases if initial guess is not sufficient.
template <class MemorySpace>
struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
template <class MemorySpace, class GeometryType>
struct BoundaryIndexSpace<MemorySpace, RegionBoundary<GeometryType>>
{
using index_view_type = Kokkos::View<std::size_t*, MemorySpace>;
index_view_type _view;
Expand All @@ -60,7 +111,7 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>

template <class ExecSpace, class Particles>
BoundaryIndexSpace( ExecSpace exec_space, Particles particles,
std::vector<RegionBoundary> planes,
std::vector<RegionBoundary<GeometryType>> planes,
const double initial_guess )
{
_timer.start();
Expand All @@ -85,7 +136,8 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
}

template <class ExecSpace, class Particles>
void update( ExecSpace, Particles particles, RegionBoundary plane )
void update( ExecSpace, Particles particles,
RegionBoundary<GeometryType> region )
{
auto count_host =
Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace{}, _count );
Expand All @@ -97,9 +149,7 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
Kokkos::RangePolicy<ExecSpace> policy( 0, particles.n_local );
auto index_functor = KOKKOS_LAMBDA( const std::size_t pid )
{
if ( x( pid, 0 ) >= plane.low_x && x( pid, 0 ) <= plane.high_x &&
x( pid, 1 ) >= plane.low_y && x( pid, 1 ) <= plane.high_y &&
x( pid, 2 ) >= plane.low_z && x( pid, 2 ) <= plane.high_z )
if ( region.inside( x, pid ) )
{
// Resize after count if needed.
auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 );
Expand All @@ -124,9 +174,9 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
auto time() { return _timer.time(); };
};

template <class BoundaryType, class ExecSpace, class Particles>
template <class ExecSpace, class Particles, class BoundaryType>
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;
Expand Down Expand Up @@ -159,9 +209,8 @@ struct BoundaryCondition
{
}

template <class ExecSpace, class Particles>
void update( ExecSpace exec_space, Particles particles,
RegionBoundary plane )
template <class ExecSpace, class Particles, class BoundaryType>
void update( ExecSpace exec_space, Particles particles, BoundaryType plane )
{
_index_space.update( exec_space, particles, plane );
}
Expand Down Expand Up @@ -202,9 +251,8 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
{
}

template <class ExecSpace, class Particles>
void update( ExecSpace exec_space, Particles particles,
RegionBoundary plane )
template <class ExecSpace, class Particles, class BoundaryType>
void update( ExecSpace exec_space, Particles particles, BoundaryType plane )
{
_index_space.update( exec_space, particles, plane );
}
Expand Down Expand Up @@ -247,9 +295,8 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
{
}

template <class ExecSpace, class Particles>
void update( ExecSpace exec_space, Particles particles,
RegionBoundary plane )
template <class ExecSpace, class Particles, class BoundaryType>
void update( ExecSpace exec_space, Particles particles, BoundaryType plane )
{
_index_space.update( exec_space, particles, plane );
}
Expand Down Expand Up @@ -288,7 +335,7 @@ auto createBoundaryCondition( BCTag, const double value, ExecSpace exec_space,
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace<BoundaryType>(
bc_index_type bc_indices = createBoundaryIndexSpace(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
}
Expand All @@ -304,7 +351,7 @@ auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space,
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace<BoundaryType>(
bc_index_type bc_indices = createBoundaryIndexSpace(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, UserFunctor>(
bc_indices, user_functor, force_update );
Expand Down

0 comments on commit b50761c

Please sign in to comment.