deal.II version 9.8.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
portable_hanging_nodes_internal.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2021 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#ifndef dealii_portable_hanging_nodes_internal_h
14#define dealii_portable_hanging_nodes_internal_h
15
16#include <deal.II/base/config.h>
17
18#include <deal.II/matrix_free/hanging_nodes_internal.h>
19
20#include <Kokkos_Core.hpp>
21
22
23DEAL_II_NAMESPACE_OPEN
24namespace Portable
25{
26 namespace internal
27 {
28 //------------------------------------------------------------------------//
29 // Functions for resolving the hanging node constraints on the GPU //
30 //------------------------------------------------------------------------//
31 template <unsigned int size>
32 DEAL_II_HOST_DEVICE inline unsigned int
33 index2(unsigned int i, unsigned int j)
34 {
35 return i + size * j;
36 }
37
38
39
40 template <unsigned int size>
41 DEAL_II_HOST_DEVICE inline unsigned int
42 index3(unsigned int i, unsigned int j, unsigned int k)
43 {
44 return i + size * j + size * size * k;
45 }
46
47
48
49 template <unsigned int fe_degree, unsigned int direction>
50 DEAL_II_HOST_DEVICE inline bool
51 is_constrained_dof_2d(
52 const ::internal::MatrixFreeFunctions::ConstraintKinds
53 &constraint_mask,
54 const unsigned int x_idx,
55 const unsigned int y_idx)
56 {
57 return ((direction == 0) &&
58 (((constraint_mask & ::internal::MatrixFreeFunctions::
59 ConstraintKinds::subcell_y) !=
60 ::internal::MatrixFreeFunctions::ConstraintKinds::
61 unconstrained) ?
62 (y_idx == 0) :
63 (y_idx == fe_degree))) ||
64 ((direction == 1) &&
65 (((constraint_mask & ::internal::MatrixFreeFunctions::
66 ConstraintKinds::subcell_x) !=
67 ::internal::MatrixFreeFunctions::ConstraintKinds::
68 unconstrained) ?
69 (x_idx == 0) :
70 (x_idx == fe_degree)));
71 }
72
73 template <unsigned int fe_degree, unsigned int direction>
74 DEAL_II_HOST_DEVICE inline bool
75 is_constrained_dof_3d(
76 const ::internal::MatrixFreeFunctions::ConstraintKinds
77 &constraint_mask,
78 const unsigned int x_idx,
79 const unsigned int y_idx,
80 const unsigned int z_idx,
81 const ::internal::MatrixFreeFunctions::ConstraintKinds face1_type,
82 const ::internal::MatrixFreeFunctions::ConstraintKinds face2_type,
83 const ::internal::MatrixFreeFunctions::ConstraintKinds face1,
84 const ::internal::MatrixFreeFunctions::ConstraintKinds face2,
85 const ::internal::MatrixFreeFunctions::ConstraintKinds edge)
86 {
87 const unsigned int face1_idx = (direction == 0) ? y_idx :
88 (direction == 1) ? z_idx :
89 x_idx;
90 const unsigned int face2_idx = (direction == 0) ? z_idx :
91 (direction == 1) ? x_idx :
92 y_idx;
93
94 const bool on_face1 = ((constraint_mask & face1_type) !=
95 ::internal::MatrixFreeFunctions::
96 ConstraintKinds::unconstrained) ?
97 (face1_idx == 0) :
98 (face1_idx == fe_degree);
99 const bool on_face2 = ((constraint_mask & face2_type) !=
100 ::internal::MatrixFreeFunctions::
101 ConstraintKinds::unconstrained) ?
102 (face2_idx == 0) :
103 (face2_idx == fe_degree);
104 return (
105 (((constraint_mask & face1) != ::internal::MatrixFreeFunctions::
106 ConstraintKinds::unconstrained) &&
107 on_face1) ||
108 (((constraint_mask & face2) != ::internal::MatrixFreeFunctions::
109 ConstraintKinds::unconstrained) &&
110 on_face2) ||
111 (((constraint_mask & edge) != ::internal::MatrixFreeFunctions::
112 ConstraintKinds::unconstrained) &&
113 on_face1 && on_face2));
114 }
115
116
117
118 template <unsigned int fe_degree,
119 unsigned int direction,
120 bool transpose,
121 typename Number,
122 typename ViewType>
123 DEAL_II_HOST_DEVICE inline void
124 interpolate_boundary_2d(
125 const Kokkos::TeamPolicy<
126 MemorySpace::Default::kokkos_space::execution_space>::member_type
127 &team_member,
128 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
129 constraint_weights,
130 const ::internal::MatrixFreeFunctions::ConstraintKinds
131 &constraint_mask,
132 ViewType values)
133 {
134 constexpr unsigned int n_q_points_1d = fe_degree + 1;
135 constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 2);
136
137 // Flag is true if dof is constrained for the given direction and the
138 // given face.
139 const bool constrained_face =
140 (constraint_mask &
141 (((direction == 0) ?
142 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
143 ::internal::MatrixFreeFunctions::ConstraintKinds::
144 unconstrained) |
145 ((direction == 1) ?
146 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
147 ::internal::MatrixFreeFunctions::ConstraintKinds::
148 unconstrained))) !=
149 ::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained;
150
151 Number tmp[n_q_points];
152 Kokkos::parallel_for(
153 Kokkos::TeamThreadRange(team_member, n_q_points),
154 [&](const int &q_point) {
155 const unsigned int x_idx = q_point % n_q_points_1d;
156 const unsigned int y_idx = q_point / n_q_points_1d;
157
158 const auto this_type =
159 (direction == 0) ?
160 ::internal::MatrixFreeFunctions::ConstraintKinds::
161 subcell_x :
162 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y;
163
164 const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
165 tmp[q_point] = 0;
166
167 // Flag is true if for the given direction, the dof is constrained
168 // with the right type and is on the correct side (left (= 0) or right
169 // (= fe_degree))
170 const bool constrained_dof =
171 is_constrained_dof_2d<fe_degree, direction>(constraint_mask,
172 x_idx,
173 y_idx);
174
175 if (constrained_face && constrained_dof)
176 {
177 const bool type = (constraint_mask & this_type) !=
178 ::internal::MatrixFreeFunctions::
179 ConstraintKinds::unconstrained;
180
181 if (type)
182 {
183 for (unsigned int i = 0; i <= fe_degree; ++i)
184 {
185 const unsigned int real_idx =
186 (direction == 0) ? index2<n_q_points_1d>(i, y_idx) :
187 index2<n_q_points_1d>(x_idx, i);
188
189 const Number w =
190 transpose ?
191 constraint_weights[i * n_q_points_1d + interp_idx] :
192 constraint_weights[interp_idx * n_q_points_1d + i];
193 tmp[q_point] += w * values[real_idx];
194 }
195 }
196 else
197 {
198 for (unsigned int i = 0; i <= fe_degree; ++i)
199 {
200 const unsigned int real_idx =
201 (direction == 0) ? index2<n_q_points_1d>(i, y_idx) :
202 index2<n_q_points_1d>(x_idx, i);
203
204 const Number w =
205 transpose ?
206 constraint_weights[(fe_degree - i) * n_q_points_1d +
207 fe_degree - interp_idx] :
208 constraint_weights[(fe_degree - interp_idx) *
209 n_q_points_1d +
210 fe_degree - i];
211 tmp[q_point] += w * values[real_idx];
212 }
213 }
214 }
215 });
216
217 // The synchronization is done for all the threads in one team with
218 // each team being assigned to one element.
219 team_member.team_barrier();
220 Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, n_q_points),
221 [&](const int &q_point) {
222 const unsigned int x_idx = q_point % n_q_points_1d;
223 const unsigned int y_idx = q_point / n_q_points_1d;
224 const bool constrained_dof =
225 is_constrained_dof_2d<fe_degree, direction>(
226 constraint_mask, x_idx, y_idx);
227 if (constrained_face && constrained_dof)
228 values[index2<fe_degree + 1>(x_idx, y_idx)] =
229 tmp[q_point];
230 });
231
232 team_member.team_barrier();
233 }
234
235
236
237 template <unsigned int fe_degree,
238 unsigned int direction,
239 bool transpose,
240 typename Number,
241 typename ViewType>
242 DEAL_II_HOST_DEVICE inline void
243 interpolate_boundary_3d(
244 const Kokkos::TeamPolicy<
245 MemorySpace::Default::kokkos_space::execution_space>::member_type
246 &team_member,
247 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
248 constraint_weights,
249 const ::internal::MatrixFreeFunctions::ConstraintKinds
250 constraint_mask,
251 ViewType values)
252 {
253 constexpr unsigned int n_q_points_1d = fe_degree + 1;
254 constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 3);
255
256 const auto this_type =
257 (direction == 0) ?
258 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x :
259 (direction == 1) ?
260 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
261 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z;
262 const auto face1_type =
263 (direction == 0) ?
264 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
265 (direction == 1) ?
266 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
267 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x;
268 const auto face2_type =
269 (direction == 0) ?
270 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
271 (direction == 1) ?
272 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x :
273 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y;
274
275 // If computing in x-direction, need to match against face_y or
276 // face_z
277 const auto face1 =
278 (direction == 0) ?
279 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
280 (direction == 1) ?
281 ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
282 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x;
283 const auto face2 =
284 (direction == 0) ?
285 ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
286 (direction == 1) ?
287 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
288 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y;
289 const auto edge =
290 (direction == 0) ?
291 ::internal::MatrixFreeFunctions::ConstraintKinds::edge_x :
292 (direction == 1) ?
293 ::internal::MatrixFreeFunctions::ConstraintKinds::edge_y :
294 ::internal::MatrixFreeFunctions::ConstraintKinds::edge_z;
295 const auto constrained_face = constraint_mask & (face1 | face2 | edge);
296
297 Number tmp[n_q_points];
298 Kokkos::parallel_for(
299 Kokkos::TeamThreadRange(team_member, n_q_points),
300 [&](const int &q_point) {
301 const unsigned int x_idx = q_point % n_q_points_1d;
302 const unsigned int y_idx = (q_point / n_q_points_1d) % n_q_points_1d;
303 const unsigned int z_idx = q_point / (n_q_points_1d * n_q_points_1d);
304
305 const unsigned int interp_idx = (direction == 0) ? x_idx :
306 (direction == 1) ? y_idx :
307 z_idx;
308 const bool constrained_dof =
309 is_constrained_dof_3d<fe_degree, direction>(constraint_mask,
310 x_idx,
311 y_idx,
312 z_idx,
313 face1_type,
314 face2_type,
315 face1,
316 face2,
317 edge);
318 tmp[q_point] = 0;
319 if ((constrained_face != ::internal::MatrixFreeFunctions::
320 ConstraintKinds::unconstrained) &&
321 constrained_dof)
322 {
323 const bool type = (constraint_mask & this_type) !=
324 ::internal::MatrixFreeFunctions::
325 ConstraintKinds::unconstrained;
326 if (type)
327 {
328 for (unsigned int i = 0; i <= fe_degree; ++i)
329 {
330 const unsigned int real_idx =
331 (direction == 0) ?
332 index3<fe_degree + 1>(i, y_idx, z_idx) :
333 (direction == 1) ?
334 index3<fe_degree + 1>(x_idx, i, z_idx) :
335 index3<fe_degree + 1>(x_idx, y_idx, i);
336
337 const Number w =
338 transpose ?
339 constraint_weights[i * n_q_points_1d + interp_idx] :
340 constraint_weights[interp_idx * n_q_points_1d + i];
341 tmp[q_point] += w * values[real_idx];
342 }
343 }
344 else
345 {
346 for (unsigned int i = 0; i <= fe_degree; ++i)
347 {
348 const unsigned int real_idx =
349 (direction == 0) ?
350 index3<n_q_points_1d>(i, y_idx, z_idx) :
351 (direction == 1) ?
352 index3<n_q_points_1d>(x_idx, i, z_idx) :
353 index3<n_q_points_1d>(x_idx, y_idx, i);
354
355 const Number w =
356 transpose ?
357 constraint_weights[(fe_degree - i) * n_q_points_1d +
358 fe_degree - interp_idx] :
359 constraint_weights[(fe_degree - interp_idx) *
360 n_q_points_1d +
361 fe_degree - i];
362 tmp[q_point] += w * values[real_idx];
363 }
364 }
365 }
366 });
367
368 // The synchronization is done for all the threads in one team with
369 // each team being assigned to one element.
370 team_member.team_barrier();
371
372 Kokkos::parallel_for(
373 Kokkos::TeamThreadRange(team_member, n_q_points),
374 [&](const int &q_point) {
375 const unsigned int x_idx = q_point % n_q_points_1d;
376 const unsigned int y_idx = (q_point / n_q_points_1d) % n_q_points_1d;
377 const unsigned int z_idx = q_point / (n_q_points_1d * n_q_points_1d);
378 const bool constrained_dof =
379 is_constrained_dof_3d<fe_degree, direction>(constraint_mask,
380 x_idx,
381 y_idx,
382 z_idx,
383 face1_type,
384 face2_type,
385 face1,
386 face2,
387 edge);
388 if ((constrained_face != ::internal::MatrixFreeFunctions::
389 ConstraintKinds::unconstrained) &&
390 constrained_dof)
391 values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = tmp[q_point];
392 });
393
394 team_member.team_barrier();
395 }
396
397
398
406 template <int dim,
407 int fe_degree,
408 bool transpose,
409 typename Number,
410 typename ViewType>
411 DEAL_II_HOST_DEVICE void
412 resolve_hanging_nodes(
413 const Kokkos::TeamPolicy<
414 MemorySpace::Default::kokkos_space::execution_space>::member_type
415 &team_member,
416 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
417 constraint_weights,
418 const ::internal::MatrixFreeFunctions::ConstraintKinds
419 constraint_mask,
420 ViewType values)
421 {
422 if constexpr (dim == 2)
423 {
424 interpolate_boundary_2d<fe_degree, 0, transpose>(team_member,
425 constraint_weights,
426 constraint_mask,
427 values);
428
429 interpolate_boundary_2d<fe_degree, 1, transpose>(team_member,
430 constraint_weights,
431 constraint_mask,
432 values);
433 }
434 else if constexpr (dim == 3)
435 {
436 // Interpolate y and z faces (x-direction)
437 interpolate_boundary_3d<fe_degree, 0, transpose>(team_member,
438 constraint_weights,
439 constraint_mask,
440 values);
441 // Interpolate x and z faces (y-direction)
442 interpolate_boundary_3d<fe_degree, 1, transpose>(team_member,
443 constraint_weights,
444 constraint_mask,
445 values);
446 // Interpolate x and y faces (z-direction)
447 interpolate_boundary_3d<fe_degree, 2, transpose>(team_member,
448 constraint_weights,
449 constraint_mask,
450 values);
451 }
452 }
453 } // namespace internal
454} // namespace Portable
455
456DEAL_II_NAMESPACE_CLOSE
457#endif
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_HOST_DEVICE
#define DEAL_II_HOST_DEVICE
Definition config.h:170
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:629
hanging_nodes_internal.h
Portable::internal::is_constrained_dof_2d
bool is_constrained_dof_2d(const ::internal::MatrixFreeFunctions::ConstraintKinds &constraint_mask, const unsigned int x_idx, const unsigned int y_idx)
Definition portable_hanging_nodes_internal.h:51
Portable::internal::interpolate_boundary_3d
void interpolate_boundary_3d(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights, const ::internal::MatrixFreeFunctions::ConstraintKinds constraint_mask, ViewType values)
Definition portable_hanging_nodes_internal.h:243
Portable::internal::is_constrained_dof_3d
bool is_constrained_dof_3d(const ::internal::MatrixFreeFunctions::ConstraintKinds &constraint_mask, const unsigned int x_idx, const unsigned int y_idx, const unsigned int z_idx, const ::internal::MatrixFreeFunctions::ConstraintKinds face1_type, const ::internal::MatrixFreeFunctions::ConstraintKinds face2_type, const ::internal::MatrixFreeFunctions::ConstraintKinds face1, const ::internal::MatrixFreeFunctions::ConstraintKinds face2, const ::internal::MatrixFreeFunctions::ConstraintKinds edge)
Definition portable_hanging_nodes_internal.h:75
Portable::internal::index3
unsigned int index3(unsigned int i, unsigned int j, unsigned int k)
Definition portable_hanging_nodes_internal.h:42
Portable::internal::resolve_hanging_nodes
void resolve_hanging_nodes(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights, const ::internal::MatrixFreeFunctions::ConstraintKinds constraint_mask, ViewType values)
Definition portable_hanging_nodes_internal.h:412
Portable::internal::index2
unsigned int index2(unsigned int i, unsigned int j)
Definition portable_hanging_nodes_internal.h:33
Portable::internal::interpolate_boundary_2d
void interpolate_boundary_2d(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights, const ::internal::MatrixFreeFunctions::ConstraintKinds &constraint_mask, ViewType values)
Definition portable_hanging_nodes_internal.h:124
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
n_q_points
*  const unsigned int n_q_points
Definition nonlinear-heat_transfer_with_AD_NOX.h:1386
Generated by doxygen 1.16.1