Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/boundaryops/normal_derivative.jl @ 1656:89456aa6fa80 feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Flesh out normal_derivative implementation some more.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 27 Jun 2024 09:06:01 +0200
parents f4dc17cfafce
children 9c84e97de895
comparison
equal deleted inserted replaced
1655:51f23a0b07fa 1656:89456aa6fa80
28 scaled_stencil = scale(closure_stencil,h_inv) 28 scaled_stencil = scale(closure_stencil,h_inv)
29 return BoundaryOperator(g, scaled_stencil, boundary) 29 return BoundaryOperator(g, scaled_stencil, boundary)
30 end 30 end
31 31
32 function normal_derivative(g::MappedGrid, stencil_set::StencilSet, boundary) 32 function normal_derivative(g::MappedGrid, stencil_set::StencilSet, boundary)
33 g⁻¹ = geometric_tensor_inverse(g) # Extract boundary part 33 b_indices = boundary_indices(logicalgrid(g), boundary)
34 k = NaN # Dimension of boundary 34
35 mapreduce(1:ndims(g)) do i 35 k = grid_id(boundary)
36 gᵏⁱ = componentview(g⁻¹,k,i) 36
37 gᵏᵏ = componentview(g⁻¹,k,k) 37
38 # ∂ξᵢ = ... 38 # Compute the weights for the logival derivatives
39 DiagonalTensor(gᵏⁱ./sqrt.(gᵏᵏ)) * ∂ξᵢ # Should the metric expression be mapped lazily? 39 g⁻¹ = geometric_tensor_inverse(g)
40 α = map(CartesianIndices(g⁻¹)[b_indices...]) do I # TODO: Fix iterator here
41 gᵏⁱ = g⁻¹[I][k,:]
42 gᵏᵏ = g⁻¹[I][k,k]
43
44 gᵏⁱ./sqrt(gᵏᵏ)
45 end
46
47
48 mapreduce(+,1:ndims(g)) do i
49 if i == k
50 ∂ξᵢ = normal_derivative(logicalgrid(g), stencil_set, boundary)
51 else
52 ∂ξᵢ = first_derivative(logicalgrid(g), stencil_set, i)
53 end
54
55 αᵢ = componentview(α,i)
56 DiagonalTensor(αᵢ) ∘ ∂ξᵢ
40 end 57 end
41 end 58 end