GitHub

Examples

2D Heat Conduction

Steady-state heat conduction on a 1m × 1m square with fixed temperatures on each edge.

Steady-State

using WhatsThePoint
using WhatsThePoint: coords
import WhatsThePoint as WTP
using Macchiato
using Unitful: m, °, ustrip
using CairoMakie

# Geometry
dx = 1/33 * m
part = PointBoundary(rectangle(1m, 1m)...)
split_surface!(part, 75°)
cloud = discretize(part, ConstantSpacing(dx), alg=VanDerSandeFornberg())

# Boundary conditions & model
bcs = Dict(
    :surface1 => Temperature(0.0),    # bottom
    :surface2 => Temperature(0.0),    # right
    :surface3 => Temperature(100.0),  # top
    :surface4 => Temperature(0.0),    # left
)
domain = Domain(cloud, bcs, SolidEnergy(k=1.0, ρ=1.0, cₚ=1.0))

# Solve
sim = Simulation(domain)
run!(sim)
T = temperature(sim)

# Visualize
pts = points(cloud)
x = [ustrip(coords(pt).x) for pt in pts]
y = [ustrip(coords(pt).y) for pt in pts]

fig = Figure(; size=(800, 700))
ax = Axis(fig[1, 1]; title="Temperature", xlabel="x [m]", ylabel="y [m]", aspect=DataAspect())
sc = scatter!(ax, x, y; color=T, colormap=:inferno, markersize=8)
Colorbar(fig[1, 2], sc; label="T")
fig
Example block output

Transient

The same geometry and BCs can be run as a transient simulation by providing a time step and stop time:

sim = Simulation(domain, Transient(Δt=0.001, stop_time=0.01))
set!(sim, T=0.0)
run!(sim)

T_final = temperature(sim)
1560-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.8947562017443811
 0.8711682133179561
 0.9074225760886687
 0.9190665696551967
 0.917163078579284
 0.9419273137658155
 0.8848407440212035
 0.9329468938632319
 0.9737792580648601

2D Cantilever Beam (Linear Elasticity)

A cantilever beam under end shear, validated against the Timoshenko analytical solution.

Problem Setup

Geometry: L × 2D beam, x ∈ [0, L], y ∈ [-D, D]. The left end is clamped (prescribed displacement from the exact solution), the right end has a parabolic shear traction, and the top/bottom surfaces are traction-free.

Timoshenko beam solution (plane stress):

u(x,y) = -P/(6EI) [y((6L-3x)x + (2+ν)(y²-D²))]
v(x,y) =  P/(6EI) [3νy²(L-x) + (4+5ν)D²x + (3L-x)x²]

where I = 2D³/3 is the second moment of area.

Full Example

using WhatsThePoint
using WhatsThePoint: coords
import WhatsThePoint as WTP
using Macchiato
using RadialBasisFunctions: PHS
using Unitful: m, °, ustrip
using LinearAlgebra
using Statistics: mean
using CairoMakie

# Problem parameters
L = 8.0    # Beam length
D = 1.0    # Half-height
P = 1000.0 # Applied load
E_val = 1e7
ν_val = 0.3
I = 2D^3 / 3

# Analytical solution
u_exact(x, y) = -P / (6E_val * I) * y * ((6L - 3x) * x + (2 + ν_val) * (y^2 - D^2))
v_exact(x, y) = P / (6E_val * I) * (3ν_val * y^2 * (L - x) + (4 + 5ν_val) * D^2 * x + (3L - x) * x^2)

# Geometry
dx = 0.1 * m
rx = dx:dx:((L * m) - dx)
ry = dx:dx:((2D * m) - dx)

p_bot = [WTP.Point(i, -D * m) for i in rx]
n_bot = [WTP.Vec(0.0, -1.0) for _ in rx]
p_right = [WTP.Point(L * m, -D * m + i) for i in ry]
n_right = [WTP.Vec(1.0, 0.0) for _ in ry]
p_top = [WTP.Point(i, D * m) for i in reverse(rx)]
n_top = [WTP.Vec(0.0, 1.0) for _ in rx]
p_left = [WTP.Point(0.0m, -D * m + i) for i in reverse(ry)]
n_left = [WTP.Vec(-1.0, 0.0) for _ in ry]

pts = vcat(p_bot, p_right, p_top, p_left)
nrms = vcat(n_bot, n_right, n_top, n_left)
areas = fill(dx, length(pts))

part = PointBoundary(pts, nrms, areas)
split_surface!(part, 75°)

cloud = WTP.discretize(part, ConstantSpacing(dx), alg=VanDerSandeFornberg())

# Boundary conditions
bc_left(x, t) = (u_exact(x[1], x[2]), v_exact(x[1], x[2]))
bc_right(x, t) = (0.0, P * (D^2 - x[2]^2) / (2I))

bcs = Dict(
    :surface1 => TractionFree(),          # bottom: free surface
    :surface2 => Traction(bc_right),       # right: parabolic shear
    :surface3 => TractionFree(),          # top: free surface
    :surface4 => Displacement(bc_left),   # left: clamped
)

# Model and solve
model = LinearElasticity(E=E_val, ν=ν_val)
domain = Domain(cloud, bcs, model)

sim = Simulation(domain)
run!(sim; basis=PHS(3; poly_deg=3))

# Extract results
ux_sim, uy_sim = displacement(sim)
N = length(cloud)

# Compare with analytical solution
pts = points(cloud)
ux_ana = [u_exact(ustrip(coords(pt).x), ustrip(coords(pt).y)) for pt in pts]
uy_ana = [v_exact(ustrip(coords(pt).x), ustrip(coords(pt).y)) for pt in pts]

println("Mean absolute error uₓ: ", mean(abs.(ux_sim .- ux_ana)))
println("Mean absolute error uᵧ: ", mean(abs.(uy_sim .- uy_ana)))

# Visualize displacement magnitude
x = [ustrip(coords(pt).x) for pt in pts]
y = [ustrip(coords(pt).y) for pt in pts]
displacement_mag = sqrt.(ux_sim .^ 2 .+ uy_sim .^ 2)

fig = Figure(; size=(1000, 500))
ax = Axis(fig[1, 1]; title="Displacement Magnitude ‖u‖", xlabel="x", ylabel="y", aspect=DataAspect())
sc = scatter!(ax, x, y; color=displacement_mag, colormap=:viridis, markersize=6)
Colorbar(fig[1, 2], sc; label="‖u‖")
fig
Example block output