GitHub

Getting started

Info

Step-by-step instructions are provided after the code.

using Tesserae
import Plots

function main()

    # Material constants
    E  = 1000.0                 # Young's modulus
    ν  = 0.3                    # Poisson's ratio
    λ  = (E*ν) / ((1+ν)*(1-2ν)) # Lame's first parameter
    μ  = E / 2(1 + ν)           # Shear modulus
    ρ⁰ = 1000.0                 # Initial density

    # Grid and particle properties
    GridProp = @NamedTuple begin
        x  :: Vec{2, Float64} # Position
        m  :: Float64         # Mass
        mv :: Vec{2, Float64} # Momentum
        f  :: Vec{2, Float64} # Force
        v  :: Vec{2, Float64} # Velocity
        vⁿ :: Vec{2, Float64} # Velocity at t = tⁿ
    end
    ParticleProp = @NamedTuple begin
        x  :: Vec{2, Float64}                           # Position
        m  :: Float64                                   # Mass
        V⁰ :: Float64                                   # Initial volume
        v  :: Vec{2, Float64}                           # Velocity
        ∇v :: SecondOrderTensor{2, Float64, 4}          # Velocity gradient
        F  :: SecondOrderTensor{2, Float64, 4}          # Deformation gradient
        σ  :: SymmetricSecondOrderTensor{2, Float64, 3} # Cauchy stress
    end

    # Mesh
    mesh = CartesianMesh(0.05, (0,1), (0,1))

    # Background grid
    grid = generate_grid(GridProp, mesh)

    # Particles
    particles = let
        pts = generate_particles(ParticleProp, mesh; alg=GridSampling())
        pts.V⁰ .= volume(mesh) / length(pts) # Set initial volume

        # Left and right disks
        r = 0.2 # Radius
        lhs = findall(x -> norm(@. x-r    ) < r, pts.x)
        rhs = findall(x -> norm(@. x-(1-r)) < r, pts.x)

        # Set initial velocities
        pts.v[lhs] .= Ref(Vec( 0.1, 0.1))
        pts.v[rhs] .= Ref(Vec(-0.1,-0.1))

        pts[[lhs; rhs]]
    end
    @. particles.m = ρ⁰ * particles.V⁰
    @. particles.F = one(particles.F)

    # Interpolation
    mpvalues = map(p -> MPValue(BSpline(Linear()), mesh), eachindex(particles))

    # Create animation by `Plots.@gif`
    Δt = 0.001
    Plots.@gif for t in range(0, 4, step=Δt)

        # Update basis function values
        for p in eachindex(particles)
            update!(mpvalues[p], particles.x[p], mesh)
        end

        @P2G grid=>i particles=>p mpvalues=>ip begin
            m[i]  = @∑ w[ip] * m[p]
            mv[i] = @∑ w[ip] * m[p] * v[p]
            f[i]  = @∑ -V⁰[p] * det(F[p]) * σ[p] * ∇w[ip]
        end

        @. grid.vⁿ = grid.mv / grid.m
        @. grid.v  = grid.vⁿ + (grid.f / grid.m) * Δt

        @G2P grid=>i particles=>p mpvalues=>ip begin
            v[p] += @∑ w[ip] * (v[i] - vⁿ[i])
            ∇v[p] = @∑ v[i] ⊗ ∇w[ip]
            x[p] += @∑ w[ip] * v[i] * Δt
        end

        for p in eachindex(particles)
            Δϵₚ = symmetric(particles.∇v[p]) * Δt
            particles.F[p]  = (I + particles.∇v[p]*Δt) * particles.F[p]
            particles.σ[p] += λ*tr(Δϵₚ)*I + 2μ*Δϵₚ # Linear elastic material
        end

        # Plot results
        Plots.scatter(
            reinterpret(Tuple{Float64,Float64}, particles.x),
            lims = (0,1),
            ticks = 0:0.2:1,
            minorgrid = true,
            minorticks = 4,
            aspect_ratio = :equal,
            legend = false,
        )
    end every 100

end

main()
Example block output

Grid and particle generation

1. Grid and particle properties

Before generating grid and particles, we must define their properties. This can be done by simply defining a NamedTuple for the grid and particles, respectively, as follows:

GridProp = @NamedTuple begin
    x  :: Vec{2, Float64} # Position
    m  :: Float64         # Mass
    mv :: Vec{2, Float64} # Momentum
    f  :: Vec{2, Float64} # Force
    v  :: Vec{2, Float64} # Velocity
    vⁿ :: Vec{2, Float64} # Velocity at t = tⁿ
end
ParticleProp = @NamedTuple begin
    x  :: Vec{2, Float64}                           # Position
    m  :: Float64                                   # Mass
    V⁰ :: Float64                                   # Initial volume
    v  :: Vec{2, Float64}                           # Velocity
    ∇v :: SecondOrderTensor{2, Float64, 4}          # Velocity gradient
    F  :: SecondOrderTensor{2, Float64, 4}          # Deformation gradient
    σ  :: SymmetricSecondOrderTensor{2, Float64, 3} # Cauchy stress
end

These properties are fully customizable, allowing users to define any variables. However, two conditions must be met: (1) the property type must be of isbitstype, and (2) the first variable must represent the position of the grid nodes and particles. These position values are automatically set during the grid and particle generation process.

julia> isbitstype(GridProp)true
julia> isbitstype(ParticleProp)true
Info

The struct can also be used instead of NamedTuple as follows:

struct GridProp
    x  :: Vec{2, Float64} # Position
    m  :: Float64         # Mass
    mv :: Vec{2, Float64} # Momentum
    f  :: Vec{2, Float64} # Force
    v  :: Vec{2, Float64} # Velocity
    vⁿ :: Vec{2, Float64} # Velocity at t = tⁿ
end

2. Mesh and grid generation

In Tesserae, a mesh and a grid differ in that a mesh only contains information about the positions of the nodes, whereas a grid includes the mesh (i.e., nodal positions) and additional user-defined variables (as defined in GridProp above).

To create Cartesian mesh, use CartesianMesh(spacing, (xmin, xmax), (ymin, ymax)...) as

mesh = CartesianMesh(0.05, (0,1), (0,1))
21×21 CartesianMesh{2, Float64, Vector{Float64}}:
 [0.0, 0.0]   [0.0, 0.05]   [0.0, 0.1]   …  [0.0, 0.95]   [0.0, 1.0]
 [0.05, 0.0]  [0.05, 0.05]  [0.05, 0.1]     [0.05, 0.95]  [0.05, 1.0]
 [0.1, 0.0]   [0.1, 0.05]   [0.1, 0.1]      [0.1, 0.95]   [0.1, 1.0]
 [0.15, 0.0]  [0.15, 0.05]  [0.15, 0.1]     [0.15, 0.95]  [0.15, 1.0]
 [0.2, 0.0]   [0.2, 0.05]   [0.2, 0.1]      [0.2, 0.95]   [0.2, 1.0]
 [0.25, 0.0]  [0.25, 0.05]  [0.25, 0.1]  …  [0.25, 0.95]  [0.25, 1.0]
 [0.3, 0.0]   [0.3, 0.05]   [0.3, 0.1]      [0.3, 0.95]   [0.3, 1.0]
 [0.35, 0.0]  [0.35, 0.05]  [0.35, 0.1]     [0.35, 0.95]  [0.35, 1.0]
 [0.4, 0.0]   [0.4, 0.05]   [0.4, 0.1]      [0.4, 0.95]   [0.4, 1.0]
 [0.45, 0.0]  [0.45, 0.05]  [0.45, 0.1]     [0.45, 0.95]  [0.45, 1.0]
 ⋮                                       ⋱                ⋮
 [0.6, 0.0]   [0.6, 0.05]   [0.6, 0.1]      [0.6, 0.95]   [0.6, 1.0]
 [0.65, 0.0]  [0.65, 0.05]  [0.65, 0.1]     [0.65, 0.95]  [0.65, 1.0]
 [0.7, 0.0]   [0.7, 0.05]   [0.7, 0.1]      [0.7, 0.95]   [0.7, 1.0]
 [0.75, 0.0]  [0.75, 0.05]  [0.75, 0.1]  …  [0.75, 0.95]  [0.75, 1.0]
 [0.8, 0.0]   [0.8, 0.05]   [0.8, 0.1]      [0.8, 0.95]   [0.8, 1.0]
 [0.85, 0.0]  [0.85, 0.05]  [0.85, 0.1]     [0.85, 0.95]  [0.85, 1.0]
 [0.9, 0.0]   [0.9, 0.05]   [0.9, 0.1]      [0.9, 0.95]   [0.9, 1.0]
 [0.95, 0.0]  [0.95, 0.05]  [0.95, 0.1]     [0.95, 0.95]  [0.95, 1.0]
 [1.0, 0.0]   [1.0, 0.05]   [1.0, 0.1]   …  [1.0, 0.95]   [1.0, 1.0]

Using this mesh, the grid can be generated as

grid = generate_grid(GridProp, mesh)
21×21 StructArray(::CartesianMesh{2, Float64, Vector{Float64}}, ::Matrix{Float64}, ::Matrix{Vec{2, Float64}}, ::Matrix{Vec{2, Float64}}, ::Matrix{Vec{2, Float64}}, ::Matrix{Vec{2, Float64}}) with eltype @NamedTuple{x::Vec{2, Float64}, m::Float64, mv::Vec{2, Float64}, f::Vec{2, Float64}, v::Vec{2, Float64}, vⁿ::Vec{2, Float64}}:
 (x = [0.0, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])   …  (x = [0.0, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.05, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])     (x = [0.05, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.1, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.1, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.15, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])     (x = [0.15, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.2, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.2, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.25, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])  …  (x = [0.25, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.3, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.3, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.35, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])     (x = [0.35, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.4, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.4, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.45, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])     (x = [0.45, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 ⋮                                                                                             ⋱  ⋮
 (x = [0.6, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.6, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.65, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])     (x = [0.65, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.7, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.7, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.75, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])  …  (x = [0.75, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.8, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.8, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.85, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])     (x = [0.85, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.9, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])      (x = [0.9, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [0.95, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])     (x = [0.95, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])
 (x = [1.0, 0.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])   …  (x = [1.0, 1.0], m = 0.0, mv = [0.0, 0.0], f = [0.0, 0.0], v = [0.0, 0.0], vⁿ = [0.0, 0.0])

This grid is a StructArray with an element type of GridProp. Thus, each variable defined in GridProp can be accessed using the . notation as follows:

julia> grid.v21×21 Matrix{Vec{2, Float64}}:
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]  …  [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]  …  [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 ⋮                                   ⋱                          ⋮
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]  …  [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]     [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
 [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]  …  [0.0, 0.0]  [0.0, 0.0]  [0.0, 0.0]
Info

Note that grid.x simply returns the mesh.

julia> grid.x === meshtrue

3. Particle generation

The particles can be generated using generate_particles function:

pts = generate_particles(ParticleProp, mesh; alg=GridSampling())
1600-element StructArray(::Vector{Vec{2, Float64}}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Vec{2, Float64}}, ::Vector{Tensor{Tuple{2, 2}, Float64, 2, 4}}, ::Vector{Tensor{Tuple{2, 2}, Float64, 2, 4}}, ::Vector{SymmetricSecondOrderTensor{2, Float64, 3}}) with eltype @NamedTuple{x::Vec{2, Float64}, m::Float64, V⁰::Float64, v::Vec{2, Float64}, ∇v::Tensor{Tuple{2, 2}, Float64, 2, 4}, F::Tensor{Tuple{2, 2}, Float64, 2, 4}, σ::SymmetricSecondOrderTensor{2, Float64, 3}}:
 (x = [0.0125, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.0375, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.0625, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.0875, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.1125, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.1375, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.1625, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.1875, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.2125, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.2375, 0.0125], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 ⋮
 (x = [0.7875, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.8125, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.8375, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.8625, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.8875, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.9125, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.9375, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.9625, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])
 (x = [0.9875, 0.9875], m = 0.0, V⁰ = 0.0, v = [0.0, 0.0], ∇v = [0.0 0.0; 0.0 0.0], F = [0.0 0.0; 0.0 0.0], σ = [0.0 0.0; 0.0 0.0])

This pts is also a StructArray, similar to grid.

Info

In the generate_particles function, particles are generated across the entire mesh domain. Consequently, any unnecessary particles need to be removed.

Basis function values

In Tesserae, the basis function values are stored in MPValue. For example, MPValue with the linear basis function can be constructed as

julia> mp = MPValue(BSpline(Linear()), mesh)MPValue:
  Interpolation: BSpline(Linear())
  Property names: w::Matrix{Float64}, ∇w::Matrix{Vec{2, Float64}}
  Neighboring nodes: CartesianIndices((1:0, 1:0))

This mp can be updated by passing the particle position to the update! function:

julia> update!(mp, particles.x[1], mesh)MPValue:
  Interpolation: BSpline(Linear())
  Property names: w::Matrix{Float64}, ∇w::Matrix{Vec{2, Float64}}
  Neighboring nodes: CartesianIndices((3:4, 1:2))
Info

After updating mp, you can check the partition of unity $\sum_i w_{ip} = 1$:

julia> sum(mp.w)1.0

and the linear field reproduction $\sum_i w_{ip} \bm{x}_i = \bm{x}_p$:

julia> nodeindices = neighboringnodes(mp)CartesianIndices((3:4, 1:2))
julia> sum(eachindex(nodeindices)) do ip
           i = nodeindices[ip]
           mp.w[ip] * mesh[i]
       end2-element Vec{2, Float64}:
 0.1375
 0.0125

For the sake of performance, it's best to prepare the same number of MPValues as there are particles. This means that each particle has its own storage for the basis function values.

mpvalues = map(p -> MPValue(BSpline(Linear()), mesh), eachindex(particles))
Info

It is also possible to construct MPValues with Structure-Of-Arrays (SOA) layout using generate_mpvalues.

julia> mpvalues = generate_mpvalues(BSpline(Linear()), mesh, length(particles))416-element MPValueVector:
  Interpolation: BSpline(Linear())
  Property names: w, ∇w

This SoA layout for MPValues is generally preferred for performance, although it cannot be resized.

Transfer between grid and particles

Particle-to-grid transfer

For the particle-to-grid transfer, the @P2G macro is useful:

@P2G grid=>i particles=>p mpvalues=>ip begin
    m[i]  = @∑ w[ip] * m[p]
    mv[i] = @∑ w[ip] * m[p] * v[p]
    f[i]  = @∑ -V⁰[p] * det(F[p]) * σ[p] * ∇w[ip]
end

This macro expands to roughly the following code:

@. grid.m  = zero(grid.m)
@. grid.mv = zero(grid.mv)
@. grid.f  = zero(grid.f)
for p in eachindex(particles)
    mp = mpvalues[p]
    nodeindices = neighboringnodes(mp)
    for ip in eachindex(nodeindices)
        i = nodeindices[ip]
        grid.m[i]  += mp.w[ip] * particles.m[p]
        grid.mv[i] += mp.w[ip] * particles.m[p] * particles.v[p]
        grid.f[i]  += -particles.V⁰[p] * det(particles.F[p]) * particles.σ[p] * mp.∇w[ip]
    end
end

Grid-to-particle transfer

Similar to the particle-to-grid transfer, the @G2P macro is provided for grid-to-particle transfer:

@G2P grid=>i particles=>p mpvalues=>ip begin
    v[p] += @∑ w[ip] * (v[i] - vⁿ[i])
    ∇v[p] = @∑ v[i] ⊗ ∇w[ip]
    x[p] += @∑ w[ip] * v[i] * Δt
end

This macro expands to roughly the following code:

for p in eachindex(particles)
    mp = mpvalues[p]
    nodeindices = neighboringnodes(mp)
    Δvₚ = zero(eltype(particles.v))
    ∇vₚ = zero(eltype(particles.∇v))
    Δxₚ = zero(eltype(particles.x))
    for ip in eachindex(nodeindices)
        i = nodeindices[ip]
        Δvₚ += mp.w[ip] * (grid.v[i] - grid.vⁿ[i])
        ∇vₚ += grid.v[i] ⊗ mp.∇w[ip]
        Δxₚ += mp.w[ip] * grid.v[i] * Δt
    end
    particles.v[p] += Δvₚ
    particles.∇v[p] = ∇vₚ
    particles.x[p] += Δxₚ
end