Mercurial > repos > public > sbplib_julia

comparison src/Grids/mapped_grid.jl @ 1566:b9c7bab94241 feature/grids/manifolds

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/curvilinear
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 25 Apr 2024 13:22:18 +0200
parents 5d32ecb98db8
children 64baaf29ae4e 063a2bfb03da
comparison
equal deleted inserted replaced
1565:c3425b4302b8 1566:b9c7bab94241
1 struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D}
2 logicalgrid::GT
3 physicalcoordinates::CT
4 jacobian::JT
5 end
6
7 jacobian(g::MappedGrid) = g.jacobian
8 logicalgrid(g::MappedGrid) = g.logicalgrid
9
10
11 # Indexing interface
12 Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...]
13 Base.eachindex(g::MappedGrid) = eachindex(g.logicalgrid)
14
15 Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d)
16 Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d)
17
18 # Iteration interface
19
20 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates)
21 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state)
22
23 Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}()
24 Base.length(g::MappedGrid) = length(g.logicalgrid)
25 Base.size(g::MappedGrid) = size(g.logicalgrid)
26 Base.size(g::MappedGrid, d) = size(g.logicalgrid, d)
27
28 boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logicalgrid)
29 boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logicalgrid, id)
30
31 function boundary_grid(g::MappedGrid, id::TensorGridBoundary)
32 b_indices = boundary_indices(g.logicalgrid, id)
33
34 # Calculate indices of needed jacobian components
35 D = ndims(g)
36 all_indices = SVector{D}(1:D)
37 free_variable_indices = deleteat(all_indices, grid_id(id))
38 jacobian_components = (:, free_variable_indices)
39
40 # Create grid function for boundary grid jacobian
41 boundary_jacobian = componentview((@view g.jacobian[b_indices...]) , jacobian_components...)
42 boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...]
43
44 return MappedGrid(
45 boundary_grid(g.logicalgrid, id),
46 boundary_physicalcoordinates,
47 boundary_jacobian,
48 )
49 end
50
51 # TBD: refine and coarsen could be implemented once we have a simple manifold implementation.
52 # Before we do, we should consider the overhead of including such a field in the mapped grid struct.
53
54 function mapped_grid(x, J, size...)
55 D = length(size)
56 lg = equidistant_grid(size, ntuple(i->0., D), ntuple(i->1., D))
57 return MappedGrid(
58 lg,
59 map(x,lg),
60 map(J,lg),
61 )
62 end
63
64 function jacobian_determinant(g::MappedGrid)
65 return map(jacobian(g)) do ∂x∂ξ
66 det(∂x∂ξ)
67 end
68 end
69
70 function geometric_tensor(g::MappedGrid)
71 return map(jacobian(g)) do ∂x∂ξ
72 ∂x∂ξ'*∂x∂ξ
73 end
74 end
75
76 function geometric_tensor_inverse(g::MappedGrid)
77 return map(jacobian(g)) do ∂x∂ξ
78 inv(∂x∂ξ'*∂x∂ξ)
79 end
80 end
81