1+ #include " utilities.hpp"
2+ #include < stdexcept>
3+
4+ namespace specfem ::assembly::boundaries_impl {
5+
6+ bool is_on_boundary (specfem::mesh_entity::dim3::type type, int iz, int iy,
7+ int ix, int ngllz, int nglly, int ngllx) {
8+ return (type == specfem::mesh_entity::dim3::type::top && iz == ngllz - 1 ) ||
9+ (type == specfem::mesh_entity::dim3::type::bottom && iz == 0 ) ||
10+ (type == specfem::mesh_entity::dim3::type::left && ix == 0 ) ||
11+ (type == specfem::mesh_entity::dim3::type::right && ix == ngllx - 1 ) ||
12+ (type == specfem::mesh_entity::dim3::type::front && iy == nglly - 1 ) ||
13+ (type == specfem::mesh_entity::dim3::type::back && iy == 0 ) ||
14+ // Edges
15+ (type == specfem::mesh_entity::dim3::type::bottom_left && iz == 0 &&
16+ ix == 0 ) ||
17+ (type == specfem::mesh_entity::dim3::type::bottom_right && iz == 0 &&
18+ ix == ngllx - 1 ) ||
19+ (type == specfem::mesh_entity::dim3::type::top_right &&
20+ iz == ngllz - 1 && ix == ngllx - 1 ) ||
21+ (type == specfem::mesh_entity::dim3::type::top_left &&
22+ iz == ngllz - 1 && ix == 0 ) ||
23+ (type == specfem::mesh_entity::dim3::type::front_bottom &&
24+ iy == nglly - 1 && iz == 0 ) ||
25+ (type == specfem::mesh_entity::dim3::type::front_top &&
26+ iy == nglly - 1 && iz == ngllz - 1 ) ||
27+ (type == specfem::mesh_entity::dim3::type::front_left &&
28+ iy == nglly - 1 && ix == 0 ) ||
29+ (type == specfem::mesh_entity::dim3::type::front_right &&
30+ iy == nglly - 1 && ix == ngllx - 1 ) ||
31+ (type == specfem::mesh_entity::dim3::type::back_bottom && iy == 0 &&
32+ iz == 0 ) ||
33+ (type == specfem::mesh_entity::dim3::type::back_top && iy == 0 &&
34+ iz == ngllz - 1 ) ||
35+ (type == specfem::mesh_entity::dim3::type::back_left && iy == 0 &&
36+ ix == 0 ) ||
37+ (type == specfem::mesh_entity::dim3::type::back_right && iy == 0 &&
38+ ix == ngllx - 1 ) ||
39+ // Corners
40+ (type == specfem::mesh_entity::dim3::type::bottom_front_left &&
41+ iz == 0 && iy == nglly - 1 && ix == 0 ) ||
42+ (type == specfem::mesh_entity::dim3::type::bottom_front_right &&
43+ iz == 0 && iy == nglly - 1 && ix == ngllx - 1 ) ||
44+ (type == specfem::mesh_entity::dim3::type::bottom_back_left &&
45+ iz == 0 && iy == 0 && ix == 0 ) ||
46+ (type == specfem::mesh_entity::dim3::type::bottom_back_right &&
47+ iz == 0 && iy == 0 && ix == ngllx - 1 ) ||
48+ (type == specfem::mesh_entity::dim3::type::top_front_left &&
49+ iz == ngllz - 1 && iy == nglly - 1 && ix == 0 ) ||
50+ (type == specfem::mesh_entity::dim3::type::top_front_right &&
51+ iz == ngllz - 1 && iy == nglly - 1 && ix == ngllx - 1 ) ||
52+ (type == specfem::mesh_entity::dim3::type::top_back_left &&
53+ iz == ngllz - 1 && iy == 0 && ix == 0 ) ||
54+ (type == specfem::mesh_entity::dim3::type::top_back_right &&
55+ iz == ngllz - 1 && iy == 0 && ix == ngllx - 1 );
56+ }
57+
58+ std::tuple<std::array<type_real, 3 >, type_real> get_boundary_face_and_weight (
59+ specfem::mesh_entity::dim3::type type,
60+ const std::array<type_real, 3 > &weights,
61+ const specfem::point::jacobian_matrix<specfem::dimension::type::dim3, true ,
62+ false > &point_jacobian_matrix) {
63+
64+ if (type == specfem::mesh_entity::dim3::type::bottom_front_left ||
65+ type == specfem::mesh_entity::dim3::type::bottom_back_left ||
66+ type == specfem::mesh_entity::dim3::type::top_front_left ||
67+ type == specfem::mesh_entity::dim3::type::top_back_left ||
68+ type == specfem::mesh_entity::dim3::type::front_left ||
69+ type == specfem::mesh_entity::dim3::type::back_left ||
70+ type == specfem::mesh_entity::dim3::type::bottom_left ||
71+ type == specfem::mesh_entity::dim3::type::top_left ||
72+ type == specfem::mesh_entity::dim3::type::left) {
73+ const auto normal = point_jacobian_matrix.compute_normal (
74+ specfem::mesh_entity::dim3::type::left);
75+ const std::array<type_real, 3 > face_normal = { normal (0 ), normal (1 ),
76+ normal (2 ) };
77+ return std::make_tuple (face_normal, weights[1 ] * weights[2 ]);
78+ }
79+
80+ if (type == specfem::mesh_entity::dim3::type::bottom_front_right ||
81+ type == specfem::mesh_entity::dim3::type::bottom_back_right ||
82+ type == specfem::mesh_entity::dim3::type::top_front_right ||
83+ type == specfem::mesh_entity::dim3::type::top_back_right ||
84+ type == specfem::mesh_entity::dim3::type::front_right ||
85+ type == specfem::mesh_entity::dim3::type::back_right ||
86+ type == specfem::mesh_entity::dim3::type::bottom_right ||
87+ type == specfem::mesh_entity::dim3::type::top_right ||
88+ type == specfem::mesh_entity::dim3::type::right) {
89+ const auto normal = point_jacobian_matrix.compute_normal (
90+ specfem::mesh_entity::dim3::type::right);
91+ const std::array<type_real, 3 > face_normal = { normal (0 ), normal (1 ),
92+ normal (2 ) };
93+ return std::make_tuple (face_normal, weights[1 ] * weights[2 ]);
94+ }
95+
96+ if (type == specfem::mesh_entity::dim3::type::top_front_left ||
97+ type == specfem::mesh_entity::dim3::type::top_front_right ||
98+ type == specfem::mesh_entity::dim3::type::top_back_left ||
99+ type == specfem::mesh_entity::dim3::type::top_back_right ||
100+ type == specfem::mesh_entity::dim3::type::front_top ||
101+ type == specfem::mesh_entity::dim3::type::back_top ||
102+ type == specfem::mesh_entity::dim3::type::top_left ||
103+ type == specfem::mesh_entity::dim3::type::top_right ||
104+ type == specfem::mesh_entity::dim3::type::top) {
105+ const auto normal = point_jacobian_matrix.compute_normal (
106+ specfem::mesh_entity::dim3::type::top);
107+ const std::array<type_real, 3 > face_normal = { normal (0 ), normal (1 ),
108+ normal (2 ) };
109+ return std::make_tuple (face_normal, weights[0 ] * weights[2 ]);
110+ }
111+
112+ if (type == specfem::mesh_entity::dim3::type::bottom_front_left ||
113+ type == specfem::mesh_entity::dim3::type::bottom_front_right ||
114+ type == specfem::mesh_entity::dim3::type::bottom_back_left ||
115+ type == specfem::mesh_entity::dim3::type::bottom_back_right ||
116+ type == specfem::mesh_entity::dim3::type::front_bottom ||
117+ type == specfem::mesh_entity::dim3::type::back_bottom ||
118+ type == specfem::mesh_entity::dim3::type::bottom_left ||
119+ type == specfem::mesh_entity::dim3::type::bottom_right ||
120+ type == specfem::mesh_entity::dim3::type::bottom) {
121+ const auto normal = point_jacobian_matrix.compute_normal (
122+ specfem::mesh_entity::dim3::type::bottom);
123+ const std::array<type_real, 3 > face_normal = { normal (0 ), normal (1 ),
124+ normal (2 ) };
125+ return std::make_tuple (face_normal, weights[0 ] * weights[2 ]);
126+ }
127+
128+ if (type == specfem::mesh_entity::dim3::type::bottom_front_left ||
129+ type == specfem::mesh_entity::dim3::type::bottom_front_right ||
130+ type == specfem::mesh_entity::dim3::type::top_front_left ||
131+ type == specfem::mesh_entity::dim3::type::top_front_right ||
132+ type == specfem::mesh_entity::dim3::type::front_left ||
133+ type == specfem::mesh_entity::dim3::type::front_right ||
134+ type == specfem::mesh_entity::dim3::type::front_bottom ||
135+ type == specfem::mesh_entity::dim3::type::front_top ||
136+ type == specfem::mesh_entity::dim3::type::front) {
137+ const auto normal = point_jacobian_matrix.compute_normal (
138+ specfem::mesh_entity::dim3::type::front);
139+ const std::array<type_real, 3 > face_normal = { normal (0 ), normal (1 ),
140+ normal (2 ) };
141+ return std::make_tuple (face_normal, weights[0 ] * weights[1 ]);
142+ }
143+
144+ if (type == specfem::mesh_entity::dim3::type::bottom_back_left ||
145+ type == specfem::mesh_entity::dim3::type::bottom_back_right ||
146+ type == specfem::mesh_entity::dim3::type::top_back_left ||
147+ type == specfem::mesh_entity::dim3::type::top_back_right ||
148+ type == specfem::mesh_entity::dim3::type::back_left ||
149+ type == specfem::mesh_entity::dim3::type::back_right ||
150+ type == specfem::mesh_entity::dim3::type::back_bottom ||
151+ type == specfem::mesh_entity::dim3::type::back_top ||
152+ type == specfem::mesh_entity::dim3::type::back) {
153+ const auto normal = point_jacobian_matrix.compute_normal (
154+ specfem::mesh_entity::dim3::type::back);
155+ const std::array<type_real, 3 > face_normal = { normal (0 ), normal (1 ),
156+ normal (2 ) };
157+ return std::make_tuple (face_normal, weights[0 ] * weights[1 ]);
158+ }
159+
160+ throw std::invalid_argument (" Error: Unknown boundary type" );
161+ }
162+
163+ } // namespace specfem::assembly::boundaries_impl
0 commit comments