Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1136:470a70a6c1e6 feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Reformulate boundary routines (sat tensors) for Laplace and fix export
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 11 Oct 2022 18:16:42 +0200
parents a8c8517a310f
children bdcdbd4ea9cd
comparison
equal deleted inserted replaced
1135:05b1d6fd6401 1136:470a70a6c1e6
52 Δ += second_derivative(grid, inner_stencil, closure_stencils, d) 52 Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
53 end 53 end
54 return Δ 54 return Δ
55 end 55 end
56 56
57 sat(Δ::Laplace, g::EquidistantGrid, bc::BoundaryCondition{Neumann}) = neumann_sat(Δ, g, bc.id) 57 """
58 sat_tensors(Δ::Laplace, g::EquidistantGrid, bc::NeumannCondition)
58 59
59 function neumann_sat(Δ::Laplace, g::EquidistantGrid{1}, id::BoundaryIdentifier) 60 Returns anonymous functions for construction the `LazyTensorApplication`s
61 recuired in order to impose a Neumann boundary condition.
62
63 See also: [`sat`,`NeumannCondition`](@ref).
64 """
65 function BoundaryConditions.sat_tensors(Δ::Laplace, g::EquidistantGrid, bc::NeumannCondition)
66 id = bc.id
60 set = Δ.stencil_set 67 set = Δ.stencil_set
61 H⁻¹ = inverse_inner_product(g,set) 68 H⁻¹ = inverse_inner_product(g,set)
69 Hᵧ = inner_product(boundary_grid(g, id), set)
62 e = boundary_restriction(g, set, id) 70 e = boundary_restriction(g, set, id)
63 d = normal_derivative(g, set, id) 71 d = normal_derivative(g, set, id)
64 return f(u,data) = H⁻¹∘e'∘(d*u .- data) 72
65 #closure = H⁻¹∘e'∘d 73 closure(u) = H⁻¹*e'*Hᵧ*d*u
66 #penalty = -H⁻¹∘e' 74 penalty(g) = -H⁻¹*e'*Hᵧ*g
67 #return closure, penalty 75 return closure, penalty
68 end 76 end
69
70
71
72 function neumann_sat(Δ::Laplace, g::EquidistantGrid, id::BoundaryIdentifier)
73 set = Δ.stencil_set
74 H⁻¹ = inverse_inner_product(g, set)
75 orth_dims = Tuple(filter(i -> i != dimension(g), 1:dimension(g)))
76 Hᵧ = inner_product(restrict(g, orth_dims...), set)
77 e = boundary_restriction(g, set, id)
78 d = normal_derivative(g, set, id)
79 return f(u,data) = H⁻¹∘e'∘Hᵧ∘(d*u .- data)
80 #closure = H⁻¹∘e'∘Hᵧ∘d
81 #penalty = -H⁻¹∘e'∘Hᵧ
82 #return closure, penalty
83 end