Simple boundary conditions

README.md

@@ -10,3 +10,108 @@ Julia package for solving PDEs via meshless methods.
 > [!CAUTION]
 > This is under heavy development and does not have robust testing for all methods
 yet. Use at your own risk.
+
+## Contributing
+
+We welcome contributions! The package is designed to be easily extensible, especially for adding new physics domains and boundary conditions.
+
+### Adding New Physics Domains
+
+The boundary condition system uses a trait-based design that separates mathematical BC types (Dirichlet, Neumann, Robin) from physics domains (Energy, Fluids, Mechanics, etc.). This makes it easy to add new physics without code duplication.
+
+#### Boundary Condition System Structure
+
+```
+src/boundary_conditions/
+├── core/                          # Core infrastructure (reusable)
+│   ├── physics_traits.jl          # Physics domain trait system
+│   ├── bc_hierarchy.jl            # Mathematical BC type hierarchy
+│   └── generic_types.jl           # Generic BC implementations
+├── numerical/                     # Numerical methods
+│   └── derivatives.jl             # Derivative computation
+├── energy.jl                      # Energy/thermal BCs
+├── fluids.jl                      # Fluid dynamics BCs
+├── walls.jl                       # Wall BCs
+└── boundary_conditions.jl         # Main orchestrator
+```
+
+#### Example: Adding Structural Mechanics
+
+To add a new physics domain (e.g., Structural Mechanics), follow these steps:
+
+**1. Define the physics domain in `core/physics_traits.jl`:**
+
+```julia
+"""
+    MechanicsPhysics <: PhysicsDomain
+
+Physics domain for structural mechanics models and boundary conditions.
+"""
+struct MechanicsPhysics <: PhysicsDomain end
+
+# Define compatibility rules
+is_compatible(::MechanicsPhysics, ::MechanicsPhysics) = true
+```
+
+**2. Create a new file `src/boundary_conditions/mechanics.jl`:**
+
+```julia
+"""
+    Displacement{T} <: Dirichlet
+
+Prescribed displacement boundary condition.
+"""
+const Displacement{T} = FixedValue{MechanicsPhysics, T}
+
+# Constructor for convenience
+Displacement(value::T) where {T} = FixedValue{MechanicsPhysics, T}(value)
+
+Base.show(io::IO, bc::Displacement) = print(io, "Displacement: $(bc.value)")
+
+"""
+    Traction{Q} <: Neumann
+
+Prescribed traction (force per unit area) boundary condition.
+"""
+const Traction{Q} = Flux{MechanicsPhysics, Q}
+
+Traction(flux::Q) where {Q} = Flux{MechanicsPhysics, Q}(flux)
+
+Base.show(io::IO, bc::Traction) = print(io, "Traction: $(bc.flux)")
+
+"""
+    ZeroStress <: Neumann
+
+Stress-free boundary: ∂u/∂n = 0
+"""
+const ZeroStress = ZeroGradient{MechanicsPhysics}
+
+Base.show(io::IO, ::ZeroStress) = print(io, "ZeroStress")
+```
+
+**3. Add the physics domain trait to your model in `src/models/mechanics.jl`:**
+
+```julia
+struct ElasticSolid{E, NU} <: AbstractModel
+    E::E    # Young's modulus
+    ν::NU   # Poisson's ratio
+end
+
+# Physics domain trait
+physics_domain(::Type{<:ElasticSolid}) = MechanicsPhysics()
+```
+
+**4. Update exports in `src/MeshlessMultiphysics.jl`:**
+
+```julia
+# Add to physics domain exports
+export MechanicsPhysics
+
+# Add mechanics BCs
+include("boundary_conditions/mechanics.jl")
+export Displacement, Traction, ZeroStress
+```
+
+**That's it!** The generic implementations (`FixedValue`, `Flux`, `ZeroGradient`) provide all the mathematical machinery. You just define type aliases with your physics domain and implement the display methods.
+
+For more examples, see the existing implementations in `src/boundary_conditions/energy.jl` and `src/boundary_conditions/fluids.jl`.

examples/heat_eq_2d.jl

@@ -0,0 +1,85 @@
+using MeshlessMultiphysics
+import MeshlessMultiphysics as MM
+using RadialBasisFunctions
+using WhatsThePoint
+import WhatsThePoint as WTP
+using StaticArrays
+using LinearAlgebra
+using DifferentialEquations
+using LinearSolve
+using CUDA
+using CUDA.CUSPARSE
+using GLMakie
+using UnicodePlots
+using Unitful: m, °, ustrip
+using JLD2
+
+include("visualize_results.jl")
+
+##
+# create boundary points
+
+L = (1m, 1m)
+
+dx = 1 / 129 * m # boundary point spacing
+S = ConstantSpacing(dx)
+rx = dx:dx:(L[1] - dx)
+ry = dx:dx:(L[2] - dx)
+
+p_bot = map(i -> WTP.Point(i, 0m), rx)
+p_right = map(i -> WTP.Point(L[1], i), ry)
+p_top = map(i -> WTP.Point(i, L[2]), reverse(rx))
+p_left = map(i -> WTP.Point(0m, i), reverse(ry))
+
+n_bot = map(i -> WTP.Vec(0.0, -1.0), rx)
+n_right = map(i -> WTP.Vec(1.0, 0.0), ry)
+n_top = map(i -> WTP.Vec(0.0, 1.0), rx)
+n_left = map(i -> WTP.Vec(-1.0, 0.0), ry)
+
+p = vcat(p_bot, p_right, p_top, p_left) # points
+n = vcat(n_bot, n_right, n_top, n_left) # normals
+a = fill(dx, length(p)) # areas
+
+part = PointBoundary(p, n, a)
+split_surface!(part, 75°)
+combine_surfaces!(part, :surface3, :surface4)
+
+figsize = (1500, 1500)
+markersize = 0.0025
+visualize(part; markersize = 1.5markersize, size = figsize)
+
+##
+
+Δ = dx
+cloud = WhatsThePoint.discretize(part, ConstantSpacing(Δ), alg = VanDerSandeFornberg())
+
+conv = repel!(cloud, ConstantSpacing(Δ); α = Δ / 20, max_iters = 1000)
+display(lineplot(conv))
+
+visualize(cloud; markersize = markersize, size = figsize)
+#save(joinpath(@__DIR__, "rectangle-0.04.jld2"), Dict("cloud"=>cloud))
+
+##
+
+# physics models and boundary conditions
+h = 250 / 1e6 # W / (mm^2 K)
+T∞ = 25 + 273.15 # K
+k = 40 / 1e3 # W / (mm K)
+ρ = 7833 / 1e9  # kg/mm^3
+cₚ = 0.465 * 1e3 # J / (kg K)]
+α = k / (cₚ * ρ) # mm^2 / s
+
+bcs = Dict(
+    :surface1 => MM.Temperature(10), :surface2 => MM.Temperature(0), :surface3 => MM.Temperature(5))
+
+domain = MM.Domain(cloud, bcs, SolidEnergy(k = k, ρ = ρ, cₚ = cₚ))
+
+u0 = zeros(length(domain.cloud))
+
+# solve using LinearSolve.jl
+prob = MM.LinearProblem(domain)
+@time sol = solve(prob)
+T = sol.u
+
+#exportvtk("heat-equation-2d", pointify(cloud), [T], ["T"])
+viz_2d(domain, T; markersize = markersize, size = figsize, levels = 32)

examples/test_shadow_bc.jl

@@ -0,0 +1,183 @@
+using MeshlessMultiphysics
+import MeshlessMultiphysics as MM
+using RadialBasisFunctions
+using WhatsThePoint
+import WhatsThePoint as WTP
+using StaticArrays
+using LinearAlgebra
+using LinearSolve
+using Unitful: m, °, ustrip
+
+println("="^60)
+println("Testing Shadow Point Boundary Conditions")
+println("="^60)
+
+##
+# Create a simple square domain
+println("\n1. Creating square domain...")
+
+L = (1m, 1m)
+dx = 0.1m  # boundary point spacing
+S = ConstantSpacing(dx)
+
+# Create boundary points for a square
+rx = dx:dx:(L[1] - dx)
+ry = dx:dx:(L[2] - dx)
+
+p_bot = map(i -> WTP.Point(i, 0m), rx)
+p_right = map(i -> WTP.Point(L[1], i), ry)
+p_top = map(i -> WTP.Point(i, L[2]), reverse(rx))
+p_left = map(i -> WTP.Point(0m, i), reverse(ry))
+
+n_bot = map(i -> WTP.Vec(0.0, -1.0), rx)
+n_right = map(i -> WTP.Vec(1.0, 0.0), ry)
+n_top = map(i -> WTP.Vec(0.0, 1.0), rx)
+n_left = map(i -> WTP.Vec(-1.0, 0.0), ry)
+
+p = vcat(p_bot, p_right, p_top, p_left)
+n = vcat(n_bot, n_right, n_top, n_left)
+a = fill(dx, length(p))
+
+part = PointBoundary(p, n, a)
+split_surface!(part, 75°)
+
+println("   Created boundary with $(length(p)) points")
+println("   Surfaces: ", keys(part.surfaces))
+
+##
+# Create cloud
+println("\n2. Creating point cloud...")
+
+Δ = dx
+cloud = WhatsThePoint.discretize(part, ConstantSpacing(Δ), alg = VanDerSandeFornberg())
+
+println("   Total points in cloud: ", size(cloud))
+println("   Volume points: ", length(cloud.volume.points))
+for (key, val) in cloud.boundary.surfaces
+    println("   Surface $key: ", length(val), " points")
+end
+
+##
+# Set up boundary conditions with shadow points
+println("\n3. Setting up boundary conditions...")
+
+# Material properties
+k = 1.0      # thermal conductivity [W/(m·K)]
+ρ = 1.0      # density [kg/m³]
+cₚ = 1.0     # specific heat [J/(kg·K)]
+
+# Test Case 1: Neumann BC with shadow points
+println("\n   Test Case: Neumann (Adiabatic) with shadow points")
+bcs = Dict(
+    :surface1 => MM.Temperature(100.0),  # Dirichlet: T = 100 (bottom)
+    :surface2 => MM.Temperature(0.0),    # Dirichlet: T = 0 (right)
+    :surface3 => MM.Temperature(50.0),   # Dirichlet: T = 50 (top)
+    :surface4 => MM.Adiabatic()          # Neumann: ∂T/∂n = 0 (left, adiabatic)
+)
+
+println("   Boundary conditions:")
+for (surf, bc) in bcs
+    println("      $surf => $bc")
+end
+
+domain = MM.Domain(cloud, bcs, MM.SolidEnergy(k = k, ρ = ρ, cₚ = cₚ))
+
+##
+# Solve with standard derivative (no shadow points)
+println("\n4. Solving with standard derivative (no shadow points)...")
+
+prob = MM.LinearProblem(domain)
+println("   System size: ", size(prob.A))
+println("   Solving...")
+sol = solve(prob)
+T_standard = sol.u
+println("   ✓ Solution obtained")
+println("   Temperature range: [$(minimum(T_standard)), $(maximum(T_standard))]")
+
+##
+# Now test with shadow points
+println("\n5. Testing with FIRST ORDER shadow points...")
+
+# Same BCs as before, but now we pass the scheme to LinearProblem
+bcs_shadow1 = Dict(
+    :surface1 => MM.Temperature(100.0),
+    :surface2 => MM.Temperature(0.0),
+    :surface3 => MM.Temperature(50.0),
+    :surface4 => MM.Adiabatic()  # Adiabatic BC - scheme is specified in LinearProblem
+)
+
+println("   Boundary conditions:")
+for (surf, bc) in bcs_shadow1
+    println("      $surf => $bc")
+end
+println("   Scheme: ShadowPoints(dx, 1)")
+
+domain_shadow1 = MM.Domain(cloud, bcs_shadow1, MM.SolidEnergy(k = k, ρ = ρ, cₚ = cₚ))
+prob_shadow1 = MM.LinearProblem(
+    domain_shadow1; scheme = WTP.ShadowPoints(WTP.ConstantSpacing(dx), 1))
+println("   System size: ", size(prob_shadow1.A))
+println("   Solving...")
+sol_shadow1 = solve(prob_shadow1)
+T_shadow1 = sol_shadow1.u
+println("   ✓ Solution obtained")
+println("   Temperature range: [$(minimum(T_shadow1)), $(maximum(T_shadow1))]")
+
+##
+# Test with second order shadow points
+println("\n6. Testing with SECOND ORDER shadow points...")
+
+bcs_shadow2 = Dict(
+    :surface1 => MM.Temperature(100.0),
+    :surface2 => MM.Temperature(0.0),
+    :surface3 => MM.Temperature(50.0),
+    :surface4 => MM.Adiabatic()  # Adiabatic BC - scheme is specified in LinearProblem
+)
+
+println("   Boundary conditions:")
+for (surf, bc) in bcs_shadow2
+    println("      $surf => $bc")
+end
+println("   Scheme: ShadowPoints(dx, 2)")
+
+domain_shadow2 = MM.Domain(cloud, bcs_shadow2, MM.SolidEnergy(k = k, ρ = ρ, cₚ = cₚ))
+prob_shadow2 = MM.LinearProblem(
+    domain_shadow2; scheme = WTP.ShadowPoints(WTP.ConstantSpacing(dx), 2))
+println("   System size: ", size(prob_shadow2.A))
+println("   Solving...")
+sol_shadow2 = solve(prob_shadow2)
+T_shadow2 = sol_shadow2.u
+println("   ✓ Solution obtained")
+println("   Temperature range: [$(minimum(T_shadow2)), $(maximum(T_shadow2))]")
+
+##
+# Test Robin BC (Convection)
+println("\n7. Testing Robin BC (Convection) with shadow points...")
+
+h = 10.0  # heat transfer coefficient [W/(m²·K)]
+T_inf = 25.0  # ambient temperature [K]
+
+bcs_robin = Dict(
+    :surface1 => MM.Temperature(100.0),
+    :surface2 => MM.Temperature(0.0),
+    :surface3 => MM.Temperature(50.0),
+    :surface4 => MM.Convection(h, k, T_inf)  # Robin BC: h*T + k*∂T/∂n = h*T_inf
+)
+
+println("   Boundary conditions:")
+for (surf, bc) in bcs_robin
+    println("      $surf => $bc")
+end
+
+domain_robin = MM.Domain(cloud, bcs_robin, MM.SolidEnergy(k = k, ρ = ρ, cₚ = cₚ))
+prob_robin = MM.LinearProblem(
+    domain_robin; scheme = WTP.ShadowPoints(WTP.ConstantSpacing(dx), 1))
+println("   System size: ", size(prob_robin.A))
+println("   Solving...")
+sol_robin = solve(prob_robin)
+T_robin = sol_robin.u
+println("   ✓ Solution obtained")
+println("   Temperature range: [$(minimum(T_robin)), $(maximum(T_robin))]")
+
+println("\n" * "="^60)
+println("Test completed!")
+println("="^60)