Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 266:9ad447176ba1 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Minor work on implementation of SATs for laplace
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 05 Dec 2019 09:27:31 +0100
parents f1e90a92ad74
children f67ce2eb6019
comparison
equal deleted inserted replaced
265:4308b500d6e7 266:9ad447176ba1
179 179
180 180
181 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end 181 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
182 182
183 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} 183 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
184 e = boundary_value(L.op, Bid()) 184 e = boundary_value(L, Bid())
185 d = normal_derivative(L.op, Bid()) 185 d = normal_derivative(L, Bid())
186 Hᵧ = boundary_quadrature(L.op, Bid()) 186 Hᵧ = boundary_quadrature(L, Bid())
187 187 H⁻¹ = inverse_quadrature(L)
188 return -L.Hi*e*Hᵧ*(d'*v - g) 188 return (-H⁻¹*e*Hᵧ*(d'*v - g))[I]
189 # Need to handle d'*v - g so that it is an AbstractArray that TensorMappings can act on
190 end 189 end
191 190
192 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition 191 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
193 tau::Float64 192 tau::Float64
194 end 193 end
195 194
196 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} 195 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
197 e = boundary_value(L.op, Bid()) 196 e = boundary_value(L, Bid())
198 d = normal_derivative(L.op, Bid()) 197 d = normal_derivative(L, Bid())
199 Hᵧ = boundary_quadrature(L.op, Bid()) 198 Hᵧ = boundary_quadrature(L, Bid())
200 199 H⁻¹ = inverse_quadrature(L)
201 return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g) 200 return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I]
202 # Need to handle scalar multiplication and addition of TensorMapping 201 # Need to handle scalar multiplication and addition of TensorMapping
203 end 202 end
204 203
205 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D 204 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
206 # return apply(s.L, v, i) + 205 # return apply(s.L, v, i) +