CoCalc -- containers

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/solvers/dgsem_t8code/containers_3d.jl
⁵⁵⁹⁰ views
1
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
2
# Since these FMAs can increase the performance of many numerical algorithms,
3
# we need to opt-in explicitly.
4
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
5
@muladd begin
6
#! format: noindent
7

8
# Interpolate tree_node_coordinates to each quadrant at the specified nodes
9
function calc_node_coordinates!(node_coordinates,
10
                                mesh::T8codeMesh{3, RealT},
11
                                nodes::AbstractVector) where {RealT <: Real}
12
    # We use `StrideArray`s here since these buffers are used in performance-critical
13
    # places and the additional information passed to the compiler makes them faster
14
    # than native `Array`s.
15
    tmp1 = StrideArray(undef, real(mesh),
16
                       StaticInt(3), static_length(nodes), static_length(mesh.nodes),
17
                       static_length(mesh.nodes))
18
    matrix1 = StrideArray(undef, real(mesh),
19
                          static_length(nodes), static_length(mesh.nodes))
20
    matrix2 = similar(matrix1)
21
    matrix3 = similar(matrix1)
22
    baryweights_in = barycentric_weights(mesh.nodes)
23

24
    num_local_trees = t8_forest_get_num_local_trees(mesh.forest)
25

26
    current_index = 0
27
    for itree in 0:(num_local_trees - 1)
28
        tree_class = t8_forest_get_tree_class(mesh.forest, itree)
29
        eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class)
30
        num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree)
31
        global_itree = t8_forest_global_tree_id(mesh.forest, itree)
32

33
        for ielement in 0:(num_elements_in_tree - 1)
34
            element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement)
35
            element_level = t8_element_level(eclass_scheme, element)
36

37
            # Note, `t8_hex_len` is encoded as an integer (Morton encoding) in
38
            # relation to `t8_hex_root_len`. This line transforms the
39
            # "integer" length to a float in relation to the unit interval [0,1].
40
            element_length = t8_hex_len(element_level) / t8_hex_root_len
41

42
            element_coords = Vector{RealT}(undef, 3)
43
            t8_element_vertex_reference_coords(eclass_scheme, element, 0,
44
                                               pointer(element_coords))
45

46
            nodes_out_x = (2 *
47
                           (element_length * 0.5f0 * (nodes .+ 1) .+ element_coords[1]) .-
48
                           1)
49
            nodes_out_y = (2 *
50
                           (element_length * 0.5f0 * (nodes .+ 1) .+ element_coords[2]) .-
51
                           1)
52
            nodes_out_z = (2 *
53
                           (element_length * 0.5f0 * (nodes .+ 1) .+ element_coords[3]) .-
54
                           1)
55

56
            polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x,
57
                                             baryweights_in)
58
            polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y,
59
                                             baryweights_in)
60
            polynomial_interpolation_matrix!(matrix3, mesh.nodes, nodes_out_z,
61
                                             baryweights_in)
62

63
            multiply_dimensionwise!(view(node_coordinates, :, :, :, :,
64
                                         current_index += 1),
65
                                    matrix1, matrix2, matrix3,
66
                                    view(mesh.tree_node_coordinates, :, :, :, :,
67
                                         global_itree + 1), tmp1)
68
        end
69
    end
70

71
    return node_coordinates
72
end
73

74
# This routine was copied and adapted from `src/dgsem_p4est/containers_3d.jl`: `orientation_to_indices_p4est`.
75
function init_mortar_neighbor_ids!(mortars::P4estMortarContainer{3}, my_face,
76
                                   other_face, orientation, neighbor_ielements,
77
                                   mortar_id)
78
    # my_face and other_face are the face directions (zero-based)
79
    # of "my side" and "other side" respectively.
80
    # Face corner 0 of the face with the lower face direction connects to a corner of the other face.
81
    # The number of this corner is the orientation code in `p4est`.
82
    lower = my_face <= other_face
83

84
    # x_pos, y_neg, and z_pos are the directions in which the face has right-handed coordinates
85
    # when looked at from the outside.
86
    my_right_handed = my_face in (1, 2, 5)
87
    other_right_handed = other_face in (1, 2, 5)
88

89
    # If both or none are right-handed when looked at from the outside, they will have different
90
    # orientations when looked at from the same side of the interface.
91
    flipped = my_right_handed == other_right_handed
92

93
    # In the following illustrations, the face corner numbering of `p4est` is shown.
94
    # ξ and η are the local coordinates of the respective face.
95
    # We're looking at both faces from the same side of the interface, so that "other side"
96
    # (in the illustrations on the left) has right-handed coordinates.
97
    if !flipped
98
        if orientation == 0
99
            # Corner 0 of other side matches corner 0 of my side
100
            #   2┌──────┐3   2┌──────┐3
101
            #    │      │     │      │
102
            #    │      │     │      │
103
            #   0└──────┘1   0└──────┘1
104
            #     η            η
105
            #     ↑            ↑
106
            #     │            │
107
            #     └───> ξ      └───> ξ
108

109
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1
110
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1
111
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[3] + 1
112
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[4] + 1
113

114
        elseif ((lower && orientation == 2) # Corner 0 of my side matches corner 2 of other side
115
                ||
116
                (!lower && orientation == 1)) # Corner 0 of other side matches corner 1 of my side
117
            #   2┌──────┐3   0┌──────┐2
118
            #    │      │     │      │
119
            #    │      │     │      │
120
            #   0└──────┘1   1└──────┘3
121
            #     η            ┌───> η
122
            #     ↑            │
123
            #     │            ↓
124
            #     └───> ξ      ξ
125

126
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1
127
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[4] + 1
128
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[1] + 1
129
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[3] + 1
130

131
        elseif ((lower && orientation == 1) # Corner 0 of my side matches corner 1 of other side
132
                ||
133
                (!lower && orientation == 2)) # Corner 0 of other side matches corner 2 of my side
134
            #   2┌──────┐3   3┌──────┐1
135
            #    │      │     │      │
136
            #    │      │     │      │
137
            #   0└──────┘1   2└──────┘0
138
            #     η                 ξ
139
            #     ↑                 ↑
140
            #     │                 │
141
            #     └───> ξ     η <───┘
142

143
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[3] + 1
144
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1
145
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[4] + 1
146
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[2] + 1
147

148
        else # orientation == 3
149
            # Corner 0 of my side matches corner 3 of other side and
150
            # corner 0 of other side matches corner 3 of my side.
151
            #   2┌──────┐3   1┌──────┐0
152
            #    │      │     │      │
153
            #    │      │     │      │
154
            #   0└──────┘1   3└──────┘2
155
            #     η           ξ <───┐
156
            #     ↑                 │
157
            #     │                 ↓
158
            #     └───> ξ           η
159

160
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[4] + 1
161
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[3] + 1
162
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[2] + 1
163
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[1] + 1
164
        end
165
    else # flipped
166
        if orientation == 0
167
            # Corner 0 of other side matches corner 0 of my side
168
            #   2┌──────┐3   1┌──────┐3
169
            #    │      │     │      │
170
            #    │      │     │      │
171
            #   0└──────┘1   0└──────┘2
172
            #     η            ξ
173
            #     ↑            ↑
174
            #     │            │
175
            #     └───> ξ      └───> η
176

177
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1
178
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[3] + 1
179
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[2] + 1
180
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[4] + 1
181

182
        elseif orientation == 2
183
            # Corner 0 of my side matches corner 2 of other side and
184
            # corner 0 of other side matches corner 2 of my side.
185
            #   2┌──────┐3   0┌──────┐1
186
            #    │      │     │      │
187
            #    │      │     │      │
188
            #   0└──────┘1   2└──────┘3
189
            #     η            ┌───> ξ
190
            #     ↑            │
191
            #     │            ↓
192
            #     └───> ξ      η
193

194
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[3] + 1
195
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[4] + 1
196
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[1] + 1
197
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[2] + 1
198

199
        elseif orientation == 1
200
            # Corner 0 of my side matches corner 1 of other side and
201
            # corner 0 of other side matches corner 1 of my side.
202
            #   2┌──────┐3   3┌──────┐2
203
            #    │      │     │      │
204
            #    │      │     │      │
205
            #   0└──────┘1   1└──────┘0
206
            #     η                 η
207
            #     ↑                 ↑
208
            #     │                 │
209
            #     └───> ξ     ξ <───┘
210

211
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1
212
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1
213
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[4] + 1
214
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[3] + 1
215

216
        else # orientation == 3
217
            # Corner 0 of my side matches corner 3 of other side and
218
            # corner 0 of other side matches corner 3 of my side.
219
            #   2┌──────┐3   2┌──────┐0
220
            #    │      │     │      │
221
            #    │      │     │      │
222
            #   0└──────┘1   3└──────┘1
223
            #     η           η <───┐
224
            #     ↑                 │
225
            #     │                 ↓
226
            #     └───> ξ           ξ
227

228
            mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[4] + 1
229
            mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1
230
            mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[3] + 1
231
            mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[1] + 1
232
        end
233
    end
234
end
235
end # @muladd
236

237
Product

Resources

Company
