CoCalc -- transfinite_mappings

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/meshes/transfinite_mappings_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
# Illustration of the corner (circled), edge (braces), and face index numbering convention
9
# used in these functions.
10
#
11
#                    ⑧────────────────────────{7}────────────────────────⑦
12
#                   ╱│                                                   ╱│
13
#                  ╱ │                                                  ╱ │
14
#                 ╱  │                                                 ╱  │
15
#                ╱   │                5 (+z)                          ╱   │
16
#               ╱    │                                               ╱    │
17
#              ╱     │                                              ╱     │
18
#            {12}    │                                            {11}    │
19
#            ╱       │                                            ╱       │
20
#           ╱        │                                           ╱        │
21
#          ╱         │                    2 (+y)                ╱         │
22
#         ╱          │                                         ╱          │
23
#        ╱          {8}                                       ╱          {6}
24
#       ╱            │                                       ╱            │
25
#      ⑤─────────────────────────{3}───────────────────────⑥    4 (+x)   │
26
#      │             │                                      │             │
27
#      │             │                                      │             │
28
#      │             │                                      │             │
29
#      │    6 (-x)   │                                      │             │
30
#      │             │                                      │             │
31
#      │             │                                      │             │
32
#      │             │                                      │             │
33
#      │             │                                      │             │
34
#      │             ④────────────────────────{5}──────────│─────────────③
35
#      │            ╱                                       │            ╱
36
#      │           ╱              1 (-y)                    │           ╱
37
#     {4}         ╱                                        {2}         ╱
38
#      │         ╱                                          │         ╱
39
#      │        ╱                                           │        ╱
40
#      │      {9}                                           │      {10}
41
#      │      ╱                                             │      ╱
42
#      │     ╱                                              │     ╱  Global coordinates:
43
#      │    ╱                                               │    ╱       z
44
#      │   ╱                             3 (-z)             │   ╱        ↑   y
45
#      │  ╱                                                 │  ╱         │  ╱
46
#      │ ╱                                                  │ ╱          │ ╱
47
#      │╱                                                   │╱           └─────> x
48
#      ①───────────────────────{1}─────────────────────────②
49

50
# Transfinite mapping formula from a point (xi, eta, zeta) in reference space [-1,1]^3 to a
51
# physical coordinate (x, y, z) for a hexahedral element with straight sides
52
function straight_side_hex_map(xi, eta, zeta, corner_points)
53
    coordinate = zeros(eltype(xi), 3)
54
    for j in 1:3
55
        coordinate[j] += (0.125f0 *
56
                          (corner_points[j, 1] * (1 - xi) * (1 - eta) * (1 - zeta)
57
                           + corner_points[j, 2] * (1 + xi) * (1 - eta) * (1 - zeta)
58
                           + corner_points[j, 3] * (1 + xi) * (1 + eta) * (1 - zeta)
59
                           + corner_points[j, 4] * (1 - xi) * (1 + eta) * (1 - zeta)
60
                           + corner_points[j, 5] * (1 - xi) * (1 - eta) * (1 + zeta)
61
                           + corner_points[j, 6] * (1 + xi) * (1 - eta) * (1 + zeta)
62
                           + corner_points[j, 7] * (1 + xi) * (1 + eta) * (1 + zeta)
63
                           + corner_points[j, 8] * (1 - xi) * (1 + eta) * (1 + zeta)))
64
    end
65

66
    return coordinate
67
end
68

69
# Construct the (x, y, z) node coordinates in the volume of a straight sided hexahedral element
70
function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, element,
71
                                nodes, corners)
72
    for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes)
73
        node_coordinates[:, i, j, k, element] .= straight_side_hex_map(nodes[i],
74
                                                                       nodes[j],
75
                                                                       nodes[k],
76
                                                                       corners)
77
    end
78

79
    return node_coordinates
80
end
81

82
# Transfinite mapping formula from a point (xi, eta, zeta) in reference space [-1,1]^3 to a point
83
# (x,y,z) in physical coordinate space for a hexahedral element with general curved sides
84
# See Section 4.3
85
# - Andrew R. Winters (2014)
86
#   Discontinuous Galerkin spectral element approximations for the reflection and
87
#   transmission of waves from moving material interfaces
88
#   [PhD thesis, Florida State University](https://diginole.lib.fsu.edu/islandora/object/fsu%3A185342)
89
function transfinite_hex_map(xi, eta, zeta, face_curves::AbstractVector{<:CurvedFace})
90
    coordinate = zeros(eltype(xi), 3)
91
    face_values = zeros(eltype(xi), (3, 6))
92
    edge_values = zeros(eltype(xi), (3, 12))
93
    corners = zeros(eltype(xi), (3, 8))
94

95
    # Compute values along the face edges
96
    edge_values[:, 1] .= evaluate_at(SVector(xi, -1), face_curves[1])
97
    edge_values[:, 2] .= evaluate_at(SVector(1, zeta), face_curves[1])
98
    edge_values[:, 3] .= evaluate_at(SVector(xi, 1), face_curves[1])
99
    edge_values[:, 4] .= evaluate_at(SVector(-1, zeta), face_curves[1])
100

101
    edge_values[:, 5] .= evaluate_at(SVector(xi, -1), face_curves[2])
102
    edge_values[:, 6] .= evaluate_at(SVector(1, zeta), face_curves[2])
103
    edge_values[:, 7] .= evaluate_at(SVector(xi, 1), face_curves[2])
104
    edge_values[:, 8] .= evaluate_at(SVector(-1, zeta), face_curves[2])
105

106
    edge_values[:, 9] .= evaluate_at(SVector(eta, -1), face_curves[6])
107
    edge_values[:, 10] .= evaluate_at(SVector(eta, -1), face_curves[4])
108
    edge_values[:, 11] .= evaluate_at(SVector(eta, 1), face_curves[4])
109
    edge_values[:, 12] .= evaluate_at(SVector(eta, 1), face_curves[6])
110

111
    # Compute values on the face
112
    face_values[:, 1] .= evaluate_at(SVector(xi, zeta), face_curves[1])
113
    face_values[:, 2] .= evaluate_at(SVector(xi, zeta), face_curves[2])
114
    face_values[:, 3] .= evaluate_at(SVector(xi, eta), face_curves[3])
115
    face_values[:, 4] .= evaluate_at(SVector(eta, zeta), face_curves[4])
116
    face_values[:, 5] .= evaluate_at(SVector(xi, eta), face_curves[5])
117
    face_values[:, 6] .= evaluate_at(SVector(eta, zeta), face_curves[6])
118

119
    # Pull the eight corner values and compute the straight sided hex mapping
120
    corners[:, 1] .= face_curves[1].coordinates[:, 1, 1]
121
    corners[:, 2] .= face_curves[1].coordinates[:, end, 1]
122
    corners[:, 3] .= face_curves[2].coordinates[:, end, 1]
123
    corners[:, 4] .= face_curves[2].coordinates[:, 1, 1]
124
    corners[:, 5] .= face_curves[1].coordinates[:, 1, end]
125
    corners[:, 6] .= face_curves[1].coordinates[:, end, end]
126
    corners[:, 7] .= face_curves[2].coordinates[:, end, end]
127
    corners[:, 8] .= face_curves[2].coordinates[:, 1, end]
128

129
    coordinate_straight = straight_side_hex_map(xi, eta, zeta, corners)
130

131
    # Compute the transfinite mapping
132
    for j in 1:3
133
        # Linear interpolation between opposite faces
134
        coordinate[j] = (0.5f0 *
135
                         (face_values[j, 6] * (1 - xi) + face_values[j, 4] * (1 + xi)
136
                          + face_values[j, 1] * (1 - eta) +
137
                          face_values[j, 2] * (1 + eta)
138
                          + face_values[j, 3] * (1 - zeta) +
139
                          face_values[j, 5] * (1 + zeta)))
140

141
        # Edge corrections to ensure faces match
142
        coordinate[j] -= (0.25f0 * (edge_values[j, 1] * (1 - eta) * (1 - zeta)
143
                           + edge_values[j, 2] * (1 + xi) * (1 - eta)
144
                           + edge_values[j, 3] * (1 - eta) * (1 + zeta)
145
                           + edge_values[j, 4] * (1 - xi) * (1 - eta)
146
                           + edge_values[j, 5] * (1 + eta) * (1 - zeta)
147
                           + edge_values[j, 6] * (1 + xi) * (1 + eta)
148
                           + edge_values[j, 7] * (1 + eta) * (1 + zeta)
149
                           + edge_values[j, 8] * (1 - xi) * (1 + eta)
150
                           + edge_values[j, 9] * (1 - xi) * (1 - zeta)
151
                           + edge_values[j, 10] * (1 + xi) * (1 - zeta)
152
                           + edge_values[j, 11] * (1 + xi) * (1 + zeta)
153
                           + edge_values[j, 12] * (1 - xi) * (1 + zeta)))
154

155
        # Subtracted interior twice, so add back the straight-sided hexahedral mapping
156
        coordinate[j] += coordinate_straight[j]
157
    end
158

159
    return coordinate
160
end
161

162
# Construct the (x, y, z) node coordinates in the volume of a curved sided hexahedral element
163
function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, element,
164
                                nodes,
165
                                face_curves::AbstractVector{<:CurvedFace})
166
    for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes)
167
        node_coordinates[:, i, j, k, element] .= transfinite_hex_map(nodes[i], nodes[j],
168
                                                                     nodes[k],
169
                                                                     face_curves)
170
    end
171

172
    return node_coordinates
173
end
174
end # @muladd
175

176
Product

Resources

Company
