Mercurial > repos > public > sbplib_julia
comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1562:efa994405c38 feature/sbp_operators/laplace_curvilinear
Add attempt at laplace for mapped grids
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 25 Apr 2024 09:34:34 +0200 |
| parents | 08f06bfacd5c |
| children | d359d0d7fb12 |
comparison
equal
deleted
inserted
replaced
| 1561:6fdc81860b0c | 1562:efa994405c38 |
|---|---|
| 50 Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d) | 50 Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d) |
| 51 end | 51 end |
| 52 return Δ | 52 return Δ |
| 53 end | 53 end |
| 54 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set) | 54 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set) |
| 55 | |
| 56 | |
| 57 function laplace(g::MappedGrid, stencil_set) | |
| 58 J = jacobian_determinant(g) | |
| 59 J⁻¹ = map(inv, J) | |
| 60 | |
| 61 Jḡ = map(*, J, ggeometric_tensor_inverse(g)) | |
| 62 lg = logicalgrid(g) | |
| 63 | |
| 64 return mapreduce(+, CartesianIndices(first(ḡ))) do I | |
| 65 i,j = I[1], I[2] | |
| 66 Jgⁱʲ = componentview(Jḡ, I[1], I[2]) | |
| 67 | |
| 68 if i == j | |
| 69 J⁻¹∘second_derivative_variable(lg, Jgⁱʲ, stencil_set, i) | |
| 70 else | |
| 71 Dᵢ = first_derivative(lg, stencil_set, i) | |
| 72 Dⱼ = first_derivative(lg, stencil_set, j) | |
| 73 J⁻¹∘Dᵢ∘DiagonalTensor(Jgⁱʲ)∘Dⱼ | |
| 74 end | |
| 75 end | |
| 76 end |