Macchiato.jl/src/boundary_conditions/boundary_conditions.jl at 33b1564374423686ac4e69bc3f26e5432e7dd913 · JuliaMeshless/Macchiato.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
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