changeset 1587:aef3827ef522 feature/sbp_operators/laplace_curvilinear

Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 26 Apr 2024 22:49:36 +0200
parents d359d0d7fb12 (diff) d4a6f9effcdd (current diff)
children f6774e98d223
files src/Grids/Grids.jl
diffstat 2 files changed, 28 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Fri Apr 26 22:49:01 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -52,3 +52,25 @@
     return Δ
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
+
+
+function laplace(grid::MappedGrid, stencil_set)
+    J = jacobian_determinant(grid)
+    J⁻¹ = DiagonalTensor(map(inv, J))
+
+    Jg = map(*, J, ggeometric_tensor_inverse(grid))
+    lg = logicalgrid(grid)
+
+    return mapreduce(+, CartesianIndices(first(Jg))) do I
+        i,j = I[1], I[2]
+        Jgⁱʲ = componentview(Jg, I[1], I[2])
+
+        if i == j
+            J⁻¹∘second_derivative_variable(lg, Jgⁱʲ, stencil_set, i)
+        else
+            Dᵢ = first_derivative(lg, stencil_set, i)
+            Dⱼ = first_derivative(lg, stencil_set, j)
+            J⁻¹∘Dᵢ∘DiagonalTensor(Jgⁱʲ)∘Dⱼ
+        end
+    end
+end
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Fri Apr 26 22:49:01 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -72,12 +72,12 @@
     g_1D = equidistant_grid(0.0, 1., 101)
     g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
 
-    @testset "1D" begin
+    @testset "EquidistantGrid" begin
         Δ = laplace(g_1D, stencil_set)
         @test Δ == second_derivative(g_1D, stencil_set)
         @test Δ isa LazyTensor{Float64,1,1}
     end
-    @testset "3D" begin
+    @testset "TensorGrid" begin
         Δ = laplace(g_3D, stencil_set)
         @test Δ isa LazyTensor{Float64,3,3}
         Dxx = second_derivative(g_3D, stencil_set, 1)
@@ -86,5 +86,9 @@
         @test Δ == Dxx + Dyy + Dzz
         @test Δ isa LazyTensor{Float64,3,3}
     end
+
+    @testset "MappedGrid" begin
+        @test_broken false
+    end
 end