Simple boundary conditions

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 => Temperature(10), :surface2 => Temperature(0), :surface3 => 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/visualize_results.jl

@@ -0,0 +1,119 @@
+function viz_2d(
+        domain,
+        labels;
+        size = (1000, 1000),
+        colorrange = WhatsThePoint._get_colorrange(labels),
+        colormap = :Spectral,
+        levels = 32,
+        kwargs...
+)
+    fig = Figure(; size = size)
+    ax = Axis(fig[1, 1]; aspect = DataAspect())
+
+    cmap = Makie.cgrad(colormap, levels; categorical = true)
+
+    #boundary
+    for b in domain.boundaries
+        ids = b.second[1]
+        points = pointify(domain.cloud[b.first])
+        c = coords.(points)
+        x = map(c -> ustrip(c.x), c)
+        y = map(c -> ustrip(c.y), c)
+        meshscatter!(
+            ax,
+            ustrip.(x),
+            ustrip.(y);
+            color = labels[ids],
+            shading = Makie.NoShading,
+            colorrange = colorrange,
+            colormap = cmap,
+            kwargs...
+        )
+    end
+
+    # volume
+    start = maximum(domain.boundaries) do b
+        b[2][1][end]
+    end + 1
+    ids = start:(start + length(domain.cloud.volume) - 1)
+    c = coords.(domain.cloud.volume.points)
+    x = map(c -> ustrip(c.x), c)
+    y = map(c -> ustrip(c.y), c)
+    meshscatter!(
+        ax,
+        ustrip.(x),
+        ustrip.(y);
+        color = labels[ids],
+        shading = Makie.NoShading,
+        colorrange = colorrange,
+        colormap = cmap,
+        kwargs...
+    )
+    Colorbar(fig[1, 2]; colorrange = colorrange, colormap = cmap)
+    save("visualization.png", fig)
+    return nothing
+end
+
+function viz_3d(
+        domain,
+        labels,
+        f_filter = x -> true;
+        size = (1000, 1000),
+        colorrange = WhatsThePoint._get_colorrange(labels),
+        azimuth = 1.275π,
+        elevation = π / 8,
+        colormap = :Spectral,
+        levels = 32,
+        kwargs...
+)
+    fig = Figure(; size = size)
+    ax = Axis3(fig[1, 1]; azimuth = azimuth, elevation = elevation)
+    ax.aspect = :data
+
+    cmap = Makie.cgrad(colormap, levels; categorical = true)
+
+    for b in domain.boundaries
+        ids = b.second[1]
+        points = pointify(domain.cloud[b.first])
+        c = coords.(points)
+        filtered_ids = findall(f_filter, c)
+        c = c[filtered_ids]
+        x = map(c -> ustrip(c.x), c)
+        y = map(c -> ustrip(c.y), c)
+        z = map(c -> ustrip(c.z), c)
+        meshscatter!(
+            ax,
+            ustrip.(x),
+            ustrip.(y),
+            ustrip.(z);
+            color = labels[ids][filtered_ids],
+            colorrange = colorrange,
+            colormap = cmap,
+            kwargs...
+        )
+    end
+
+    # volume
+    start = maximum(domain.boundaries) do b
+        b[2][1][end]
+    end + 1
+    ids = start:(start + length(domain.cloud.volume) - 1)
+    c = coords.(domain.cloud.volume.points)
+    filtered_ids = findall(f_filter, c)
+    c = c[filtered_ids]
+    x = map(c -> ustrip(c.x), c)
+    y = map(c -> ustrip(c.y), c)
+    z = map(c -> ustrip(c.z), c)
+    meshscatter!(
+        ax,
+        ustrip.(x),
+        ustrip.(y),
+        ustrip.(z);
+        color = labels[ids][filtered_ids],
+        colorrange = colorrange,
+        colormap = cmap,
+        kwargs...
+    )
+    Colorbar(fig[1, 2]; colorrange = colorrange, colormap = cmap)
+    return fig
+end

src/MeshlessMultiphysics.jl

@@ -22,25 +22,25 @@ using WriteVTK
 include("utils.jl")
 
 #################### Abstract Types ####################
-abstract type AbstractBoundaryCondition end
 abstract type AbstractModel end
 
+#################### Boundary Conditions ####################
+include("boundary_conditions/boundary_derivatives.jl")
+include("boundary_conditions/boundary_conditions.jl")
+
+export AbstractBoundaryCondition
+
 #################### Domains ####################
 include("domain.jl")
 export Domain
 
-#################### Boundary Conditions ####################
-export AbstractBoundaryCondition
-
 include("boundary_conditions/walls.jl")
 export Wall
 
 include("boundary_conditions/fluids.jl")
-export FluidBoundaryCondition
 export VelocityInlet, PressureOutlet
 
 include("boundary_conditions/energy.jl")
-export EnergyBoundaryCondition
 export Adiabatic, Temperature, HeatFlux, Convection
 
 export make_bc, make_bc!

src/boundary_conditions/boundary_conditions.jl

@@ -0,0 +1,92 @@
+abstract type AbstractBoundaryCondition end
+abstract type DerivativeBoundaryCondition <: AbstractBoundaryCondition end
+
+abstract type Dirichlet <: AbstractBoundaryCondition end
+abstract type Neumann <: DerivativeBoundaryCondition end
+abstract type Robin <: DerivativeBoundaryCondition end
+
+bc_family(::Type{<:Dirichlet}) = Dirichlet
+bc_family(::Type{<:DerivativeBoundaryCondition}) = DerivativeBoundaryCondition
+
+bc_type(::Type{<:Dirichlet}) = Dirichlet
+bc_type(::Type{<:Neumann}) = Neumann
+bc_type(::Type{<:Robin}) = Robin
+
+function make_bc!(A, b, boundary::T, surf, domain, ids; kwargs...) where {T}
+    return make_bc!(bc_family(T), A, b, boundary, surf, domain, ids; kwargs...)
+end
+
+function make_bc!(::Type{Dirichlet}, A, b, boundary, surf, domain, ids; kwargs...)
+    return write_bc_dirichlet!(A, b, ids, boundary, surf)
+end
+
+function make_bc!(::Type{DerivativeBoundaryCondition}, A, b, boundary, surf, domain, ids;
+        scheme = nothing, kwargs...)
+    # scheme comes from kwargs, passed down from LinearProblem
+    weights = compute_derivative_weights(surf, domain, scheme; kwargs...)
+    return write_bc_derivative!(
+        bc_type(typeof(boundary)), A, b, ids, weights, boundary, surf; kwargs...)
+end
+
+function write_bc_derivative!(
+        ::Type{Neumann}, A, b, ids, weights, boundary, surf; kwargs...)
+    return write_bc_neumann!(A, b, ids, weights, boundary, surf; kwargs...)
+end
+
+function write_bc_derivative!(::Type{Robin}, A, b, ids, weights, boundary, surf; kwargs...)
+    return write_bc_robin!(A, b, ids, weights, boundary, surf; kwargs...)
+end
+
+function get_bc_value_at_index(value::Number, surf, ids, i)
+    # Scalar value - return as-is
+    return value
+end
+
+function get_bc_value_at_index(value::AbstractVector, surf, ids, i)
+    # Vector of values - index into it
+    local_idx = i - first(ids) + 1
+    return value[local_idx]
+end
+
+function get_bc_value_at_index(value::Function, surf, ids, i)
+    # Function - evaluate at surface point
+    local_idx = i - first(ids) + 1
+    point = coordinates(surf)[local_idx]
+    return value(point)
+end
+
+function write_bc_dirichlet!(A::AbstractMatrix{TA}, b::AbstractVector{TB},
+        ids, boundary, surf) where {TA, TB}
+    bc_value = boundary()  # Can be scalar, vector, or function
+
+    for i in ids
+        A[i, :] .= zero(TA)
+        A[i, i] = one(TA)
+        b[i] = convert(TB, get_bc_value_at_index(bc_value, surf, ids, i))
+    end
+end
+
+function write_bc_neumann!(A::AbstractMatrix{TA}, b::AbstractVector{TB},
+        ids, weights, boundary, surf; kwargs...) where {TA, TB}
+    bc_value = boundary()  # Can be scalar, vector, or function
+
+    for i in ids
+        local_idx = i - first(ids) + 1
+        A[i, :] .= weights[local_idx, :]
+        b[i] = convert(TB, get_bc_value_at_index(bc_value, surf, ids, i))
+    end
+end
+
+function write_bc_robin!(A::AbstractMatrix{TA}, b::AbstractVector{TB},
+        ids, weights, boundary, surf; kwargs...) where {TA, TB}
+    α_val = convert(TA, α(boundary))
+    β_val = convert(TA, β(boundary))
+    bc_value = boundary()  # Can be scalar, vector, or function
+
+    for i in ids
+        local_idx = i - first(ids) + 1
+        A[i, :] .= convert(TA, β_val) .* weights[local_idx, :]
+        A[i, i] += convert(TA, α_val)
+        b[i] = convert(TB, get_bc_value_at_index(bc_value, surf, ids, i))
+    end
+end