Mercurial > repos > public > sbplib_julia
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