changeset 1660:6d196fb85133 feature/sbp_operators/laplace_curvilinear

Merge feature/grids/curvilinear
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 28 Jun 2024 17:04:05 +0200
parents cc9d18a5ff2d (current diff) 3bbcd496e021 (diff)
children bdb4becac704
files src/Grids/Grids.jl src/Grids/mapped_grid.jl test/Grids/mapped_grid_test.jl
diffstat 3 files changed, 36 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/src/Grids/Grids.jl	Fri Jun 28 17:00:57 2024 +0200
+++ b/src/Grids/Grids.jl	Fri Jun 28 17:04:05 2024 +0200
@@ -44,6 +44,7 @@
 export eval_on
 export componentview
 export ArrayComponentView
+export normal
 
 export BoundaryIdentifier
 export TensorGridBoundary
--- a/src/Grids/mapped_grid.jl	Fri Jun 28 17:00:57 2024 +0200
+++ b/src/Grids/mapped_grid.jl	Fri Jun 28 17:04:05 2024 +0200
@@ -81,3 +81,27 @@
     end
 end
 
+"""
+    normal(g::MappedGrid, boundary)
+
+The outward pointing normal as a grid function on the boundary
+"""
+function normal(g::MappedGrid, boundary)
+    b_indices = boundary_indices(g, boundary)
+    σ =_boundary_sign(component_type(g), boundary)
+    return map(jacobian(g)[b_indices...]) do ∂x∂ξ
+        ∂ξ∂x = inv(∂x∂ξ)
+        k = grid_id(boundary)
+        σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
+    end
+end
+
+function _boundary_sign(T, boundary)
+    if boundary_id(boundary) == Upper()
+        return one(T)
+    elseif boundary_id(boundary) == Lower()
+        return -one(T)
+    else
+        throw(ArgumentError("The boundary identifier must be either `Lower()` or `Upper()`"))
+    end
+end
--- a/test/Grids/mapped_grid_test.jl	Fri Jun 28 17:00:57 2024 +0200
+++ b/test/Grids/mapped_grid_test.jl	Fri Jun 28 17:04:05 2024 +0200
@@ -183,4 +183,15 @@
     lg = equidistant_grid((0,0), (1,1), 10, 11)
     @test logicalgrid(mg) == lg
     @test collect(mg) == map(x̄, lg)
+
+
+    @testset "normal" begin
+        @test normal(mg, CartesianBoundary{1,Lower}()) == fill(@SVector[-1,0], 11)
+        @test normal(mg, CartesianBoundary{1,Upper}()) == fill(@SVector[1,0], 11)
+        @test normal(mg, CartesianBoundary{2,Lower}()) == fill(@SVector[0,-1], 10)
+        @test normal(mg, CartesianBoundary{2,Upper}()) ≈ map(boundary_grid(mg,CartesianBoundary{2,Upper}())|>logicalgrid) do ξ̄
+            α = 1-2ξ̄[1]
+            @SVector[α,1]/√(α^2 + 1)
+        end
+    end
 end