Initial Condition

InitialCondition(; coordinates, density, velocity=zeros(size(coordinates, 1)),
                 mass=nothing, pressure=0.0, particle_spacing=-1.0)

Struct to hold the initial configuration of the particles.

The following setups return InitialConditions for commonly used setups:

• RectangularShape
• SphereShape
• RectangularTank
• ComplexShape
• extrude_geometry

InitialConditions support the set operations union, setdiff and intersect in order to build more complex geometries.

Arguments

• coordinates: An array where the $i$-th column holds the coordinates of particle $i$.
• density: Either a vector holding the density of each particle, or a function mapping each particle's coordinates to its density, or a scalar for a constant density over all particles.

Keywords

• velocity: Either an array where the $i$-th column holds the velocity of particle $i$, or a function mapping each particle's coordinates to its velocity, or, for a constant fluid velocity, a vector holding this velocity. Velocity is constant zero by default.
• mass: Either nothing (default) to automatically compute particle mass from particle density and spacing, or a vector holding the mass of each particle, or a function mapping each particle's coordinates to its mass, or a scalar for a constant mass over all particles.
• pressure: Either a vector holding the pressure of each particle, or a function mapping each particle's coordinates to its pressure, or a scalar for a constant pressure over all particles. This is optional and only needed when using the EntropicallyDampedSPHSystem.
• particle_spacing: The spacing between the particles. This is a scalar, as the spacing is assumed to be uniform. This is only needed when using set operations on the InitialCondition or for automatic mass calculation.

Examples

# Rectangle filled with particles
initial_condition = RectangularShape(0.1, (3, 4), (-1.0, 1.0), density=1.0)

# Two spheres in one initial condition
initial_condition = union(SphereShape(0.15, 0.5, (-1.0, 1.0), 1.0),
                          SphereShape(0.15, 0.2, (0.0, 1.0), 1.0))

# Rectangle with a spherical hole
shape1 = RectangularShape(0.1, (16, 13), (-0.8, 0.0), density=1.0)
shape2 = SphereShape(0.1, 0.35, (0.0, 0.6), 1.0, sphere_type=RoundSphere())
initial_condition = setdiff(shape1, shape2)

# Intersect of a rectangle with a sphere. Note that this keeps the particles of the
# rectangle that are in the intersect, while `intersect(shape2, shape1)` would consist of
# the particles of the sphere that are in the intersect.
shape1 = RectangularShape(0.1, (16, 13), (-0.8, 0.0), density=1.0)
shape2 = SphereShape(0.1, 0.35, (0.0, 0.6), 1.0, sphere_type=RoundSphere())
initial_condition = intersect(shape1, shape2)

# Build `InitialCondition` manually
coordinates = [0.0 1.0 1.0
               0.0 0.0 1.0]
velocity = zero(coordinates)
mass = ones(3)
density = 1000 * ones(3)
initial_condition = InitialCondition(; coordinates, velocity, mass, density)

# With functions
initial_condition = InitialCondition(; coordinates, velocity=x -> 2x, mass=1.0, density=1000.0)

ComplexShape(geometry::Union{TriangleMesh, Polygon}; particle_spacing, density,
             pressure=0.0, mass=nothing, velocity=zeros(ndims(geometry)),
             point_in_geometry_algorithm=WindingNumberJacobson(; geometry,
                                                               hierarchical_winding=false,
                                                               winding_number_factor=sqrt(eps())),
             grid_offset::Real=0.0, max_nparticles=10^7,
             pad_initial_particle_grid=2particle_spacing)

Sample a complex geometry with particles. Returns an InitialCondition. Note that an initial particle grid is generated inside the bounding box of the geometry. A point_in_geometry_algorithm checks if particles are inside the geometry or not. For more information about the method see WindingNumberJacobson or WindingNumberHorman.

Arguments

• geometry: Geometry returned by load_geometry.

Keywords

• particle_spacing: Spacing between the particles.
• density: Either a function mapping each particle's coordinates to its density, or a scalar for a constant density over all particles.
• velocity: Either a function mapping each particle's coordinates to its velocity, or, for a constant fluid velocity, a vector holding this velocity. Velocity is constant zero by default.
• mass: Either nothing (default) to automatically compute particle mass from particle density and spacing, or a function mapping each particle's coordinates to its mass, or a scalar for a constant mass over all particles.
• pressure: Scalar to set the pressure of all particles to this value. This is only used by the EntropicallyDampedSPHSystem and will be overwritten when using an initial pressure function in the system.
• point_in_geometry_algorithm: Algorithm for sampling the complex geometry with particles. It basically checks whether a particle is inside an object or not. For more information see WindingNumberJacobson or WindingNumberHorman
• grid_offset: Offset of the initial particle grid of the bounding box of the geometry.
• max_nparticles: Maximum number of particles in the initial particle grid. This is only used to avoid accidentally choosing a particle_spacing that is too small for the scale of the geometry.
• pad_initial_particle_grid: Padding of the initial particle grid.

extrude_geometry(geometry; particle_spacing, direction, n_extrude::Integer,
                 velocity=zeros(length(direction)),
                 mass=nothing, density=nothing, pressure=0.0)

Extrude either a line, a plane or a shape along a specific direction. Returns an InitialCondition.

Arguments

• geometry: Either particle coordinates or an InitialCondition defining a 2D shape to extrude to a 3D volume, or two 2D points $(A, B)$ defining the interval $[A, B]$ to extrude to a plane in 2D, or three 3D points $(A, B, C)$ defining the parallelogram spanned by the vectors $\widehat{AB}$ and $\widehat {AC}$ to extrude to a parallelepiped.

Keywords

• particle_spacing: Spacing between the particles. Can be omitted when geometry is an InitialCondition (unless geometry.particle_spacing == -1).
• direction: A vector that specifies the direction in which to extrude.
• n_extrude: Number of layers of particles created in the direction of extrusion.
• velocity: Either a function mapping each particle's coordinates to its velocity, or, for a constant fluid velocity, a vector holding this velocity. Velocity is constant zero by default.
• mass: Either nothing (default) to automatically compute particle mass from particle density and spacing, or a function mapping each particle's coordinates to its mass, or a scalar for a constant mass over all particles.
• density: Either a function mapping each particle's coordinates to its density, or a scalar for a constant density over all particles.
• pressure: Scalar to set the pressure of all particles to this value. This is only used by the EntropicallyDampedSPHSystem and will be overwritten when using an initial pressure function in the system.
• tlsph: With the TotalLagrangianSPHSystem, particles need to be placed on the boundary of the shape and not one particle radius away, as for fluids. When tlsph=true, particles will be placed on the boundary of the shape.

Examples

# Extrude a line in 2D to a plane in 2D
p1 = [0.0, 0.0]
p2 = [1.0, 1.0]

direction = [-1.0, 1.0]

shape = extrude_geometry((p1, p2); direction, particle_spacing=0.1, n_extrude=4, density=1000.0)

# Extrude a parallelogram in 3D to a parallelepiped in 3D
p1 = [0.0, 0.0, 0.0]
p2 = [0.5, 1.0, 0.0]
p3 = [1.0, 0.2, 0.0]

direction = [0.0, 0.0, 1.0]

shape = extrude_geometry((p1, p2, p3); direction, particle_spacing=0.1, n_extrude=4, density=1000.0)

# Extrude a 2D shape (here: a disc) to a 3D shape (here: a cylinder)
shape = SphereShape(0.1, 0.5, (0.2, 0.4), 1000.0, n_layers=3,
                    sphere_type=RoundSphere(end_angle=pi))

direction = [0.0, 0.0, 1.0]

shape = extrude_geometry(shape; direction, particle_spacing=0.1, n_extrude=4, density=1000.0)

RectangularShape(particle_spacing, n_particles_per_dimension, min_coordinates;
                 velocity=zeros(length(n_particles_per_dimension)),
                 mass=nothing, density=nothing, pressure=0.0,
                 acceleration=nothing, state_equation=nothing,
                 tlsph=false, loop_order=nothing)

Rectangular shape filled with particles. Returns an InitialCondition.

Arguments

• particle_spacing: Spacing between the particles.
• n_particles_per_dimension: Tuple containing the number of particles in x, y and z (only 3D) direction, respectively.
• min_coordinates: Coordinates of the corner in negative coordinate directions.

Keywords

• velocity: Either a function mapping each particle's coordinates to its velocity, or, for a constant fluid velocity, a vector holding this velocity. Velocity is constant zero by default.
• mass: Either nothing (default) to automatically compute particle mass from particle density and spacing, or a function mapping each particle's coordinates to its mass, or a scalar for a constant mass over all particles.
• density: Either a function mapping each particle's coordinates to its density, or a scalar for a constant density over all particles. Obligatory when not using a state equation. Cannot be used together with state_equation.
• pressure: Scalar to set the pressure of all particles to this value. This is only used by the EntropicallyDampedSPHSystem and will be overwritten when using an initial pressure function in the system. Cannot be used together with hydrostatic pressure gradient.
• acceleration: In order to initialize particles with a hydrostatic pressure gradient, an acceleration vector can be passed. Note that only accelerations in one coordinate direction and no diagonal accelerations are supported. This will only change the pressure of the particles. When using the WeaklyCompressibleSPHSystem, pass a state_equation as well to initialize the particles with the corresponding density and mass. When using the EntropicallyDampedSPHSystem, the pressure will be overwritten when using an initial pressure function in the system. This cannot be used together with the pressure keyword argument.
• state_equation: When calculating a hydrostatic pressure gradient by setting acceleration, the state_equation will be used to set the corresponding density. Cannot be used together with density.
• tlsph: With the TotalLagrangianSPHSystem, particles need to be placed on the boundary of the shape and not one particle radius away, as for fluids. When tlsph=true, particles will be placed on the boundary of the shape.

Examples

# 2D
rectangular = RectangularShape(particle_spacing, (5, 4), (1.0, 2.0), density=1000.0)

# 2D with hydrostatic pressure gradient.
# `state_equation` has to be the same as for the WCSPH system.
state_equation = StateEquationCole(sound_speed=20.0, exponent=7, reference_density=1000.0)
rectangular = RectangularShape(particle_spacing, (5, 4), (1.0, 2.0),
                               acceleration=(0.0, -9.81), state_equation=state_equation)

# 3D
rectangular = RectangularShape(particle_spacing, (5, 4, 7), (1.0, 2.0, 3.0), density=1000.0)

RectangularTank(particle_spacing, fluid_size, tank_size, fluid_density;
                velocity=zeros(length(fluid_size)), fluid_mass=nothing,
                pressure=0.0,
                acceleration=nothing, state_equation=nothing,
                boundary_density=fluid_density,
                n_layers=1, spacing_ratio=1.0,
                min_coordinates=zeros(length(fluid_size)),
                faces=Tuple(trues(2 * length(fluid_size))))

Rectangular tank filled with a fluid to set up dam-break-style simulations.

Arguments

• particle_spacing: Spacing between the fluid particles.
• fluid_size: The dimensions of the fluid as (x, y) (or (x, y, z) in 3D).
• tank_size: The dimensions of the tank as (x, y) (or (x, y, z) in 3D).
• fluid_density: The rest density of the fluid. Will only be used as default for boundary_density when using a state equation.

Keywords

• velocity: Either a function mapping each particle's coordinates to its velocity, or, for a constant fluid velocity, a vector holding this velocity. Velocity is constant zero by default.
• fluid_mass: Either nothing (default) to automatically compute particle mass from particle density and spacing, or a function mapping each particle's coordinates to its mass, or a scalar for a constant mass over all particles.
• pressure: Scalar to set the pressure of all particles to this value. This is only used by the EntropicallyDampedSPHSystem and will be overwritten when using an initial pressure function in the system. Cannot be used together with hydrostatic pressure gradient.
• acceleration: In order to initialize particles with a hydrostatic pressure gradient, an acceleration vector can be passed. Note that only accelerations in one coordinate direction and no diagonal accelerations are supported. This will only change the pressure of the particles. When using the WeaklyCompressibleSPHSystem, pass a state_equation as well to initialize the particles with the corresponding density and mass. When using the EntropicallyDampedSPHSystem, the pressure will be overwritten when using an initial pressure function in the system. This cannot be used together with the pressure keyword argument.
• state_equation: When calculating a hydrostatic pressure gradient by setting acceleration, the state_equation will be used to set the corresponding density. Cannot be used together with density.
• boundary_density: Density of each boundary particle (by default set to the fluid density)
• n_layers: Number of boundary layers.
• spacing_ratio: Ratio of particle_spacing to boundary particle spacing. A value of 2 means that the boundary particle spacing will be half the fluid particle spacing.
• min_coordinates: Coordinates of the corner in negative coordinate directions.
• faces: By default all faces are generated. Set faces by passing a bit-array of length 4 (2D) or 6 (3D) to generate the faces in the normal direction: -x,+x,-y,+y,-z,+z.

Fields

• fluid::InitialCondition: InitialCondition for the fluid.
• boundary::InitialCondition: InitialCondition for the boundary.
• fluid_size::Tuple: Tuple containing the size of the fluid in each dimension after rounding.
• tank_size::Tuple: Tuple containing the size of the tank in each dimension after rounding.

Examples

# 2D
setup = RectangularTank(particle_spacing, (water_width, water_height),
                        (container_width, container_height), fluid_density,
                        n_layers=2, spacing_ratio=3)

# 2D with hydrostatic pressure gradient.
# `state_equation` has to be the same as for the WCSPH system.
state_equation = StateEquationCole(sound_speed=10.0, exponent=1, reference_density=1000.0)
setup = RectangularTank(particle_spacing, (water_width, water_height),
                        (container_width, container_height), fluid_density,
                        acceleration=(0.0, -9.81), state_equation=state_equation)

# 3D
setup = RectangularTank(particle_spacing, (water_width, water_height, water_depth),
                        (container_width, container_height, container_depth), fluid_density,
                        n_layers=2)

See also: reset_wall!.

reset_wall!(rectangular_tank::RectangularTank, reset_faces, positions)

The selected walls of the tank will be placed at the new positions.

Arguments

• reset_faces: Boolean tuple of 4 (in 2D) or 6 (in 3D) dimensions, similar to faces in RectangularTank.
• positions: Tuple of new positions

RoundSphere(; start_angle=0.0, end_angle=2π)

Construct a sphere (or sphere segment) by nesting perfectly round concentric spheres. The resulting ball will be perfectly round, but will not have a regular inner structure.

Keywords

• start_angle: The starting angle of the sphere segment in radians. It determines the beginning point of the segment. The default is set to 0.0 representing the positive x-axis.
• end_angle: The ending angle of the sphere segment in radians. It defines the termination point of the segment. The default is set to 2pi, completing a full sphere.

VoxelSphere()

Construct a sphere of voxels (where particles are placed in the voxel center) with a regular inner structure but corners on the surface. Essentially, a grid of particles is generated and all particles outside the sphere are removed. The resulting sphere will have a perfect inner structure, but is not perfectly round, as it will have corners (like a sphere in Minecraft).

SphereShape(particle_spacing, radius, center_position, density;
            sphere_type=VoxelSphere(), n_layers=-1, layer_outwards=false,
            cutout_min=(0.0, 0.0), cutout_max=(0.0, 0.0), tlsph=false,
            velocity=zeros(length(center_position)), mass=nothing, pressure=0.0)

Generate a sphere that is either completely filled (by default) or hollow (by passing n_layers).

With the sphere type VoxelSphere, a sphere of voxels (where particles are placed in the voxel center) with a regular inner structure but corners on the surface is created. Essentially, a grid of particles is generated and all particles outside the sphere are removed. With the sphere type RoundSphere, a perfectly round sphere with an imperfect inner structure is created.

A cuboid can be cut out of the sphere by specifying the two corners in negative and positive coordinate directions as cutout_min and cutout_max.

Arguments

• particle_spacing: Spacing between the particles.
• radius: Radius of the sphere.
• center_position: The coordinates of the center of the sphere.
• density: Either a function mapping each particle's coordinates to its density, or a scalar for a constant density over all particles.

Keywords

• sphere_type: Either VoxelSphere or RoundSphere (see explanation above).
• n_layers: Set to an integer greater than zero to generate a hollow sphere, where the shell consists of n_layers layers.
• layer_outwards: When set to false (by default), radius is the outer radius of the sphere. When set to true, radius is the inner radius of the sphere. This is only used when n_layers > 0.
• cutout_min: Corner in negative coordinate directions of a cuboid that is to be cut out of the sphere.
• cutout_max: Corner in positive coordinate directions of a cuboid that is to be cut out of the sphere.
• tlsph: With the TotalLagrangianSPHSystem, particles need to be placed on the boundary of the shape and not one particle radius away, as for fluids. When tlsph=true, particles will be placed on the boundary of the shape.
• velocity: Either a function mapping each particle's coordinates to its velocity, or, for a constant fluid velocity, a vector holding this velocity. Velocity is constant zero by default.
• mass: Either nothing (default) to automatically compute particle mass from particle density and spacing, or a function mapping each particle's coordinates to its mass, or a scalar for a constant mass over all particles.
• pressure: Either a function mapping each particle's coordinates to its pressure, or a scalar for a constant pressure over all particles. This is optional and only needed when using the EntropicallyDampedSPHSystem.

Examples

# Filled circle with radius 0.5, center in (0.2, 0.4) and a particle spacing of 0.1
SphereShape(0.1, 0.5, (0.2, 0.4), 1000.0)

# Same as before, but perfectly round
SphereShape(0.1, 0.5, (0.2, 0.4), 1000.0, sphere_type=RoundSphere())

# Hollow circle with ~3 layers, outer radius 0.5, center in (0.2, 0.4) and a particle
# spacing of 0.1.
SphereShape(0.1, 0.5, (0.2, 0.4), 1000.0, n_layers=3)

# Same as before, but perfectly round
SphereShape(0.1, 0.5, (0.2, 0.4), 1000.0, n_layers=3, sphere_type=RoundSphere())

# Hollow circle with 3 layers, inner radius 0.5, center in (0.2, 0.4) and a particle spacing
# of 0.1.
SphereShape(0.1, 0.5, (0.2, 0.4), 1000.0, n_layers=3, layer_outwards=true)

# Filled circle with radius 0.1, center in (0.0, 0.0), particle spacing 0.1, but the
# rectangle [0, 1] x [-0.2, 0.2] is cut out.
SphereShape(0.1, 1.0, (0.0, 0.0), 1000.0, cutout_min=(0.0, -0.2), cutout_max=(1.0, 0.2))

# Filled 3D sphere with radius 0.5, center in (0.2, 0.4, 0.3) and a particle spacing of 0.1
SphereShape(0.1, 0.5, (0.2, 0.4, 0.3), 1000.0)

# Same as before, but perfectly round
SphereShape(0.1, 0.5, (0.2, 0.4, 0.3), 1000.0, sphere_type=RoundSphere())