Mercurial > repos > public > sbplib_julia
comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 618:c64793f77509 feature/volume_and_boundary_operators
Move Laplace and SecondDerivative into the volumeops directory
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 07 Dec 2020 12:16:09 +0100 |
| parents | src/SbpOperators/laplace/laplace.jl@d9324671b412 |
| children | a85db383484f |
comparison
equal
deleted
inserted
replaced
| 617:f59e1732eacc | 618:c64793f77509 |
|---|---|
| 1 # """ | |
| 2 # Laplace{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} | |
| 3 # | |
| 4 # Implements the Laplace operator `L` in Dim dimensions as a tensor operator | |
| 5 # The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative | |
| 6 # tensor operators. | |
| 7 # """ | |
| 8 function Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim | |
| 9 Δ = SecondDerivative(grid, inner_stencil, closure_stencils, 1) | |
| 10 for d = 2:Dim | |
| 11 Δ += SecondDerivative(grid, inner_stencil, closure_stencils, d) | |
| 12 end | |
| 13 return Δ | |
| 14 end | |
| 15 export Laplace | |
| 16 | |
| 17 # quadrature(L::Laplace) = Quadrature(L.op, L.grid) | |
| 18 # inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid) | |
| 19 # boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId) | |
| 20 # normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId) | |
| 21 # boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId) | |
| 22 # export NormalDerivative | |
| 23 # """ | |
| 24 # NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} | |
| 25 # | |
| 26 # Implements the boundary operator `d` as a TensorMapping | |
| 27 # """ | |
| 28 # struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} | |
| 29 # op::D2{T,N,M,K} | |
| 30 # grid::EquidistantGrid{2} | |
| 31 # bId::CartesianBoundary | |
| 32 # end | |
| 33 # | |
| 34 # # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue? | |
| 35 # # Can we give special treatment to TensorMappings that go to a higher dim? | |
| 36 # function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer}) | |
| 37 # if dim(e.bId) == 1 | |
| 38 # return (UnknownDim, domain_size[1]) | |
| 39 # elseif dim(e.bId) == 2 | |
| 40 # return (domain_size[1], UnknownDim) | |
| 41 # end | |
| 42 # end | |
| 43 # LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],) | |
| 44 # | |
| 45 # # TODO: Not type stable D:< | |
| 46 # # TODO: Make this independent of dimension | |
| 47 # function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T | |
| 48 # i = I[dim(d.bId)] | |
| 49 # j = I[3-dim(d.bId)] | |
| 50 # N_i = size(d.grid)[dim(d.bId)] | |
| 51 # h_inv = inverse_spacing(d.grid)[dim(d.bId)] | |
| 52 # return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId)) | |
| 53 # end | |
| 54 # | |
| 55 # function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T | |
| 56 # u = selectdim(v,3-dim(d.bId),Int(I[1])) | |
| 57 # return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId)) | |
| 58 # end | |
| 59 # | |
| 60 # """ | |
| 61 # BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} | |
| 62 # | |
| 63 # Implements the boundary operator `q` as a TensorOperator | |
| 64 # """ | |
| 65 # export BoundaryQuadrature | |
| 66 # struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} | |
| 67 # op::D2{T,N,M,K} | |
| 68 # grid::EquidistantGrid{2} | |
| 69 # bId::CartesianBoundary | |
| 70 # end | |
| 71 # | |
| 72 # | |
| 73 # # TODO: Make this independent of dimension | |
| 74 # function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T | |
| 75 # h = spacing(q.grid)[3-dim(q.bId)] | |
| 76 # N = size(v) | |
| 77 # return apply_quadrature(q.op, h, v[I[1]], I[1], N[1]) | |
| 78 # end | |
| 79 # | |
| 80 # LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I) | |
| 81 # | |
| 82 # | |
| 83 # | |
| 84 # | |
| 85 # struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end | |
| 86 # | |
| 87 # function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} | |
| 88 # e = boundary_value(L, Bid()) | |
| 89 # d = normal_derivative(L, Bid()) | |
| 90 # Hᵧ = boundary_quadrature(L, Bid()) | |
| 91 # H⁻¹ = inverse_quadrature(L) | |
| 92 # return (-H⁻¹*e*Hᵧ*(d'*v - g))[I] | |
| 93 # end | |
| 94 # | |
| 95 # struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition | |
| 96 # tau::Float64 | |
| 97 # end | |
| 98 # | |
| 99 # function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} | |
| 100 # e = boundary_value(L, Bid()) | |
| 101 # d = normal_derivative(L, Bid()) | |
| 102 # Hᵧ = boundary_quadrature(L, Bid()) | |
| 103 # H⁻¹ = inverse_quadrature(L) | |
| 104 # return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I] | |
| 105 # # Need to handle scalar multiplication and addition of TensorMapping | |
| 106 # end | |
| 107 | |
| 108 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D | |
| 109 # return apply(s.L, v, i) + | |
| 110 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + | |
| 111 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + | |
| 112 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + | |
| 113 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) | |
| 114 # end |