API Reference

Domain{M, C}

Central container that ties together a point cloud, boundary conditions, and physics models.

Fields

• cloud::PointCloud{M, C}: The discretized geometry (boundary + interior points)
• boundaries::Dict{Symbol, Tuple{UnitRange, AbstractBoundaryCondition}}: Mapping from surface names to (index_range, bc) pairs, where index_range gives the global indices of that surface's points in the assembled system
• models::AbstractVector{<:AbstractModel}: Physics models (e.g., SolidEnergy, LinearElasticity)
• name::Symbol: Domain identifier (defaults to :domain1)

Constructors

Domain(cloud, boundaries, model)    # cloud + BCs + model(s)
Domain(cloud, model)                # cloud + model(s), no BCs

At construction, the Domain validates that:

1. Every boundary condition key matches a surface name in the point cloud
2. Every surface in the point cloud has a corresponding boundary condition entry

add!(domain::Domain, model::AbstractModel)

Append a physics model to the domain's model list.

add!(domain::Domain, boundary::AbstractBoundaryCondition, name::Symbol)

Attach a boundary condition to the named surface on the domain.

delete!(domain::Domain, model::AbstractModel)

Remove a physics model from the domain.

SolidEnergy(; k, ρ, cₚ, source=nothing)

Solid-body energy (heat) transport model.

Solves the heat equation in a solid medium:

• Steady-state: k ∇²T = -f (Poisson equation)
• Transient: ρ cₚ ∂T/∂t = k ∇²T + f

Fields

• k: Thermal conductivity
• ρ: Density
• cₚ: Specific heat capacity
• source: Optional volumetric source term f(x, t) -> value (default: nothing)

Example

model = SolidEnergy(k=50.0, ρ=7800.0, cₚ=500.0)
model = SolidEnergy(k=1.0, ρ=1.0, cₚ=1.0, source=(x, t) -> -4.0)

LinearElasticity{E, Nu, Rho, F} <: Solid

Linear isotropic elasticity model for solid mechanics (Navier-Cauchy equations).

Supports 2D plane stress formulation. The governing equations in displacement form:

(λ*+2μ) ∂²u/∂x² + μ ∂²u/∂y² + (λ*+μ) ∂²v/∂x∂y + fₓ = 0
(λ*+μ) ∂²u/∂x∂y + μ ∂²v/∂x² + (λ*+2μ) ∂²v/∂y² + fᵧ = 0

Fields

• E: Young's modulus
• ν: Poisson's ratio
• ρ: Density (optional, for body forces / dynamics)
• body_force: Body force function f(x) -> (fx, fy) (optional)

Example

model = LinearElasticity(E=200e3, ν=0.3)
model = LinearElasticity(E=200e3, ν=0.3, body_force=x -> (0.0, -9.81))

lame_parameters(model::LinearElasticity)

Compute Lamé parameters for plane stress:

• μ = E / (2(1+ν))
• λ = Eν / ((1+ν)(1-2ν))
• λ* = 2μλ / (λ+2μ) (plane stress modification)

Returns (μ, λstar).

IncompressibleNavierStokes(; μ, ρ)
IncompressibleNavierStokes(μ::Real, ρ)

Incompressible Navier-Stokes fluid model.

When μ is a plain number, it is automatically wrapped in NewtonianViscosity. Pass an AbstractViscosity subtype (e.g., CarreauYasudaViscosity) for non-Newtonian behavior.

Fields

• μ::AbstractViscosity: Dynamic viscosity model
• ρ: Fluid density

Example

model = IncompressibleNavierStokes(μ=0.001, ρ=1000.0)
model = IncompressibleNavierStokes(μ=CarreauYasudaViscosity(μ_inf=0.0035, μ_0=0.056), ρ=1060.0)

AbstractViscosity

Abstract supertype for viscosity models used with IncompressibleNavierStokes. Concrete subtypes must be callable as μ(γ̇) returning the dynamic viscosity at shear rate γ̇.

NewtonianViscosity(μ)

Constant (Newtonian) viscosity model. Returns the same viscosity μ regardless of shear rate.

Example

visc = NewtonianViscosity(0.001)  # water-like viscosity
visc(100.0)  # => 0.001

CarreauYasudaViscosity(; μ_inf, μ_0, n=0.333, λ=0.31, a=2)

Generalized Newtonian viscosity model for shear-thinning (or shear-thickening) fluids.

The Carreau-Yasuda model:

μ(γ̇) = μ_inf + (μ_0 - μ_inf) * (1 + (λ γ̇)^a)^((n - 1) / a)

Fields

• μ_inf: Infinite-shear-rate viscosity
• μ_0: Zero-shear-rate viscosity
• n: Power-law index (< 1 for shear-thinning)
• λ: Relaxation time
• a: Yasuda parameter (a = 2 recovers the Carreau model)

Example

visc = CarreauYasudaViscosity(μ_inf=0.0035, μ_0=0.056, n=0.333, λ=0.31)

AbstractSimulationMode

Abstract supertype for simulation modes. See Steady and Transient.

Steady()

Mode for steady-state simulations solved via LinearSolve.

Transient(; Δt, stop_time, solver=Tsit5())

Mode for transient (time-dependent) simulations solved via OrdinaryDiffEq.solve.

Arguments

• Δt: Time step size
• stop_time: End time for the simulation
• solver: ODE solver (default: Tsit5())

Examples

Transient(Δt=0.001, stop_time=1.0)
Transient(Δt=0.001, stop_time=1.0, solver=RK4())

AbstractModel

Abstract supertype for all PDE models. Subtype this to define a custom PDE.

See Custom PDEs for a complete walkthrough.

_num_vars(model::AbstractModel, dim) -> Int

Return the number of solution variables per point for model in dim dimensions.

Examples: 1 for scalar PDEs, dim for vector PDEs, dim + 1 for velocity + pressure.

make_system(model::AbstractModel, domain; kwargs...) -> (A, b)

Assemble the system matrix A and right-hand side b for steady-state solving.

Required for steady-state simulations. Macchiato applies boundary conditions and solves Ax = b.

make_f(model::AbstractModel, domain; kwargs...) -> f

Return an in-place ODE function f(du, u, p, t) for transient integration.

Required for transient simulations. Macchiato passes the returned function to OrdinaryDiffEq.jl.

AbstractBoundaryCondition

Base abstract type for all boundary conditions. All BCs must subtype one of: Dirichlet, Neumann, or Robin.

Dirichlet <: AbstractBoundaryCondition

Essential boundary conditions that prescribe values at the boundary. Examples: Temperature, VelocityInlet, Displacement.

Neumann <: DerivativeBoundaryCondition

Natural boundary conditions that prescribe normal derivatives: ∂u/∂n = g. Examples: HeatFlux, Adiabatic, VelocityOutlet.

Robin <: DerivativeBoundaryCondition

Mixed boundary conditions combining value and derivative: β·∂u/∂n + α·u = g. Example: Convection.

DerivativeBoundaryCondition <: AbstractBoundaryCondition

Abstract type for BCs involving derivatives (Neumann and Robin).

PrescribedValue{F<:Function} <: Dirichlet

Generic Dirichlet BC that prescribes a value via a function.

The function has signature f(x, t) -> value where:

• x: spatial coordinate of the boundary point
• t: time
• Returns the prescribed value at that location and time

Fields

• f: Function with signature (x, t) -> value
• name: Display name (e.g., :Temperature, :PrescribedValue)

For built-in physics, use named constructors (e.g., Temperature, VelocityInlet). For custom PDEs, use the unparameterized constructors directly: PrescribedValue(0.0).

PrescribedFlux{F<:Function} <: Neumann

Generic Neumann BC that prescribes a flux (normal derivative) via a function.

The flux condition is: ∂u/∂n = f(x, t)

The function has signature f(x, t) -> flux_value where:

• x: spatial coordinate of the boundary point
• t: time
• Returns the prescribed flux value at that location and time

Fields

• f: Function with signature (x, t) -> flux
• name: Display name (e.g., :HeatFlux, :PrescribedFlux)

For built-in physics, use named constructors (e.g., HeatFlux, Traction). For custom PDEs, use the unparameterized constructors directly: PrescribedFlux(1.0).

ZeroFlux <: Neumann

Generic Neumann BC with zero flux: ∂u/∂n = 0.

Represents symmetry, insulation, or fully-developed flow depending on context.

Fields

• name: Display name (e.g., :Adiabatic, :VelocityOutlet, :ZeroFlux)

For built-in physics, use named constructors (e.g., Adiabatic, VelocityOutlet). For custom PDEs, use the unparameterized constructor directly: ZeroFlux().

Temperature(value)

Prescribed temperature BC. Value can be a Number or Function (x, t) -> value.

HeatFlux(flux)

Prescribed heat flux BC: ∂T/∂n = q. Flux can be a Number or Function (x, t) -> flux.

Adiabatic()

Thermally insulated boundary: ∂T/∂n = 0.

Convection(h, k, T∞)

Convective heat transfer: h·T + k·∂T/∂n = h·T∞(x,t). T∞ can be a Number or Function (x, t) -> ambient_temp.

Displacement{F<:Function} <: Dirichlet

Prescribed displacement BC for solid mechanics. The function returns a tuple of displacement components: f(x, t) -> (ux, uy) for 2D.

Constructors

Displacement((x, t) -> (0.0, 0.0))          # Function returning tuple
Displacement(0.0, 0.0)                       # Constant displacement
Displacement(ux::Function, uy::Function)     # Per-component functions

Traction{F<:Function} <: Neumann

Prescribed traction BC for solid mechanics. The function returns a tuple of traction components: f(x, t) -> (tx, ty) for 2D.

In terms of stress: t = σ·n where n is the outward normal.

Constructors

Traction((x, t) -> (0.0, -1000.0))          # Function returning tuple
Traction(tx::Number, ty::Number)             # Constant traction
Traction(tx::Function, ty::Function)         # Per-component functions

TractionFree()

Zero-traction (free surface) BC: σ·n = 0. Convenience for Traction(0.0, 0.0).

VelocityInlet(velocity)

Prescribed velocity at inlet. Value can be a Number or Function (x, t) -> velocity.

PressureOutlet(pressure)

Prescribed pressure at outlet. Value can be a Number or Function (x, t) -> pressure.

VelocityOutlet()

Zero-gradient velocity outlet: ∂v/∂n = 0. Used for fully developed outflow.

Wall(velocity)
Wall()

No-slip wall or moving wall BC. Value can be a Number or Function (x, t) -> velocity. No arguments creates stationary wall (v=0).

Simulation{M, C, Mode}

High-level simulation container that manages solving and solution storage.

Fields

• domain::Domain{M, C}: The computational domain with models and boundary conditions
• mode::Mode: Simulation mode (Steady or Transient)
• u0: Initial condition vector (transient only)
• time: Current simulation time
• running: Whether simulation is currently running
• _solution: Solution vector

Constructors

Steady-state simulation

sim = Simulation(domain)
sim = Simulation(domain, Steady())

Transient simulation

sim = Simulation(domain, Transient(Δt=0.001, stop_time=1.0))

run!(sim::Simulation; kwargs...)

Execute the simulation. Dispatches to transient or steady-state path based on mode. Extra kwargs are forwarded to OrdinaryDiffEq.solve (transient) or LinearSolve.LinearProblem (steady).

Returns the simulation object.

set!(sim::Simulation; kwargs...)

Set initial conditions for simulation fields.

Arguments

• sim: Simulation to set initial conditions for
• kwargs: Field name/value pairs

Supported value types

• Number: Uniform value for entire field
• Function: Called as f(x) where x is coordinate vector [x, y] or [x, y, z]
• Vector: Direct assignment (must match field length)

Examples

set!(sim, T=300.0)                           # Uniform temperature
set!(sim, T=x -> 300 + 10*x[1])              # Temperature function of position
set!(sim, u=0.0, v=0.0, p=0.0)               # Multiple fields

temperature(sim) -> Vector{Float64}

Extract temperature field from simulation.

velocity(sim) -> Tuple{Vector{Float64}, ...}

Extract velocity components from simulation. Returns (u, v) for 2D or (u, v, w) for 3D.

pressure(sim) -> Vector{Float64}

Extract pressure field from simulation.

displacement(sim) -> Tuple{Vector{Float64}, ...}

Extract displacement components from simulation. Returns (ux, uy) for 2D or (ux, uy, uz) for 3D.

LinearSolve.LinearProblem(domain::Domain; scheme=nothing, kwargs...)

Construct a LinearProblem for steady-state simulation from a Domain.

Assembles the system matrix A and right-hand side b from the physics model, then applies boundary conditions by modifying the appropriate rows of A and b. The resulting system Ax = b is solved with LinearSolve.solve.

exportvtk(filename, points, data, names)

Export point-based simulation results to a VTK file.

Arguments

• filename::String: Output file path (without .vtu extension)
• points::AbstractVector: Point cloud or coordinate vector
• data::AbstractVector{<:AbstractVector}: Field data arrays to export
• names::AbstractVector: Corresponding field names (e.g., ["T", "u"])

Example

exportvtk("results/temperature", points(cloud), [T_values], ["T"])

savevtk!(vtkfile)

Write the VTK file to disk.

upwind(data, eval_points, dim[, basis]; Δ=nothing, k=autoselect_k(data, basis))
upwind(data, dim[, basis]; Δ=nothing, k=autoselect_k(data, basis))

Build an upwind finite-difference-style operator using RBF interpolation.

Computes backward, forward, and centered partial derivatives with respect to dimension dim, then returns a function (ϕ, v, θ) that blends them based on flow direction v and upwind parameter θ ∈ [0, 1] (1 = full upwind, 0 = centered).

The single-argument form upwind(data, dim) evaluates at the data points themselves.

Arguments

• data: Stencil points for RBF approximation
• eval_points: Points where the derivative is evaluated
• dim: Spatial dimension (1 = x, 2 = y, …)
• basis: Radial basis function (default: PHS(3; poly_deg=2))
• Δ: Virtual node offset distance (auto-detected if nothing)
• k: Number of nearest neighbors for stencil