deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+00:00
\(\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
166 // Flag is true if for the given direction, the dof is constrained
167 // with the right type and is on the correct side (left (= 0) or right
168 // (= fe_degree))
169 const bool constrained_dof =
170 is_constrained_dof_2d<fe_degree, direction>(constraint_mask,
171 x_idx,
172 y_idx);
173
174 if (constrained_face && constrained_dof)
175 {
176 Number sum = 0.;
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 sum += 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 sum += w * values[real_idx];
212 }
213 }
214 tmp[q_point] = sum;
215 }
216 });
217
218 // The synchronization is done for all the threads in one team with
219 // each team being assigned to one element.
220 team_member.team_barrier();
221 Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, n_q_points),
222 [&](const int &q_point) {
223 const unsigned int x_idx = q_point % n_q_points_1d;
224 const unsigned int y_idx = q_point / n_q_points_1d;
225 const bool constrained_dof =
226 is_constrained_dof_2d<fe_degree, direction>(
227 constraint_mask, x_idx, y_idx);
228 if (constrained_face && constrained_dof)
229 values[index2<fe_degree + 1>(x_idx, y_idx)] =
230 tmp[q_point];
231 });
232
233 team_member.team_barrier();
234 }
235
236
237
238 template <unsigned int fe_degree,
239 unsigned int direction,
240 bool transpose,
241 typename Number,
242 typename ViewType>
243 DEAL_II_HOST_DEVICE inline void
244 interpolate_boundary_3d(
245 const Kokkos::TeamPolicy<
246 MemorySpace::Default::kokkos_space::execution_space>::member_type
247 &team_member,
248 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
249 constraint_weights,
250 const ::internal::MatrixFreeFunctions::ConstraintKinds
251 constraint_mask,
252 ViewType values)
253 {
254 constexpr unsigned int n_q_points_1d = fe_degree + 1;
255 constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 3);
256
257 const auto this_type =
258 (direction == 0) ?
259 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x :
260 (direction == 1) ?
261 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
262 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z;
263 const auto face1_type =
264 (direction == 0) ?
265 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
266 (direction == 1) ?
267 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
268 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x;
269 const auto face2_type =
270 (direction == 0) ?
271 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
272 (direction == 1) ?
273 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x :
274 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y;
275
276 // If computing in x-direction, need to match against face_y or
277 // face_z
278 const auto face1 =
279 (direction == 0) ?
280 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
281 (direction == 1) ?
282 ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
283 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x;
284 const auto face2 =
285 (direction == 0) ?
286 ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
287 (direction == 1) ?
288 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
289 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y;
290 const auto edge =
291 (direction == 0) ?
292 ::internal::MatrixFreeFunctions::ConstraintKinds::edge_x :
293 (direction == 1) ?
294 ::internal::MatrixFreeFunctions::ConstraintKinds::edge_y :
295 ::internal::MatrixFreeFunctions::ConstraintKinds::edge_z;
296 const auto constrained_face = constraint_mask & (face1 | face2 | edge);
297
298 Number tmp[n_q_points] = {};
299 Kokkos::parallel_for(
300 Kokkos::TeamThreadRange(team_member, n_q_points),
301 [&](const int &q_point) {
302 const unsigned int x_idx = q_point % n_q_points_1d;
303 const unsigned int y_idx = (q_point / n_q_points_1d) % n_q_points_1d;
304 const unsigned int z_idx = q_point / (n_q_points_1d * n_q_points_1d);
305
306 const unsigned int interp_idx = (direction == 0) ? x_idx :
307 (direction == 1) ? y_idx :
308 z_idx;
309 const bool constrained_dof =
310 is_constrained_dof_3d<fe_degree, direction>(constraint_mask,
311 x_idx,
312 y_idx,
313 z_idx,
314 face1_type,
315 face2_type,
316 face1,
317 face2,
318 edge);
319 if ((constrained_face != ::internal::MatrixFreeFunctions::
320 ConstraintKinds::unconstrained) &&
321 constrained_dof)
322 {
323 Number sum = 0.;
324 const bool type = (constraint_mask & this_type) !=
325 ::internal::MatrixFreeFunctions::
326 ConstraintKinds::unconstrained;
327 if (type)
328 {
329 for (unsigned int i = 0; i <= fe_degree; ++i)
330 {
331 const unsigned int real_idx =
332 (direction == 0) ?
333 index3<fe_degree + 1>(i, y_idx, z_idx) :
334 (direction == 1) ?
335 index3<fe_degree + 1>(x_idx, i, z_idx) :
336 index3<fe_degree + 1>(x_idx, y_idx, i);
337
338 const Number w =
339 transpose ?
340 constraint_weights[i * n_q_points_1d + interp_idx] :
341 constraint_weights[interp_idx * n_q_points_1d + i];
342 sum += w * values[real_idx];
343 }
344 }
345 else
346 {
347 for (unsigned int i = 0; i <= fe_degree; ++i)
348 {
349 const unsigned int real_idx =
350 (direction == 0) ?
351 index3<n_q_points_1d>(i, y_idx, z_idx) :
352 (direction == 1) ?
353 index3<n_q_points_1d>(x_idx, i, z_idx) :
354 index3<n_q_points_1d>(x_idx, y_idx, i);
355
356 const Number w =
357 transpose ?
358 constraint_weights[(fe_degree - i) * n_q_points_1d +
359 fe_degree - interp_idx] :
360 constraint_weights[(fe_degree - interp_idx) *
361 n_q_points_1d +
362 fe_degree - i];
363 sum += w * values[real_idx];
364 }
365 }
366 tmp[q_point] = sum;
367 }
368 });
369
370 // The synchronization is done for all the threads in one team with
371 // each team being assigned to one element.
372 team_member.team_barrier();
373
374 Kokkos::parallel_for(
375 Kokkos::TeamThreadRange(team_member, n_q_points),
376 [&](const int &q_point) {
377 const unsigned int x_idx = q_point % n_q_points_1d;
378 const unsigned int y_idx = (q_point / n_q_points_1d) % n_q_points_1d;
379 const unsigned int z_idx = q_point / (n_q_points_1d * n_q_points_1d);
380 const bool constrained_dof =
381 is_constrained_dof_3d<fe_degree, direction>(constraint_mask,
382 x_idx,
383 y_idx,
384 z_idx,
385 face1_type,
386 face2_type,
387 face1,
388 face2,
389 edge);
390 if ((constrained_face != ::internal::MatrixFreeFunctions::
391 ConstraintKinds::unconstrained) &&
392 constrained_dof)
393 values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = tmp[q_point];
394 });
395
396 team_member.team_barrier();
397 }
398
399
400
408 template <int dim,
409 int fe_degree,
410 bool transpose,
411 typename Number,
412 typename ViewType>
413 DEAL_II_HOST_DEVICE void
414 resolve_hanging_nodes(
415 const Kokkos::TeamPolicy<
416 MemorySpace::Default::kokkos_space::execution_space>::member_type
417 &team_member,
418 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
419 constraint_weights,
420 const ::internal::MatrixFreeFunctions::ConstraintKinds
421 constraint_mask,
422 ViewType values)
423 {
424 if constexpr (dim == 2)
425 {
426 interpolate_boundary_2d<fe_degree, 0, transpose>(team_member,
427 constraint_weights,
428 constraint_mask,
429 values);
430
431 interpolate_boundary_2d<fe_degree, 1, transpose>(team_member,
432 constraint_weights,
433 constraint_mask,
434 values);
435 }
436 else if constexpr (dim == 3)
437 {
438 // Interpolate y and z faces (x-direction)
439 interpolate_boundary_3d<fe_degree, 0, transpose>(team_member,
440 constraint_weights,
441 constraint_mask,
442 values);
443 // Interpolate x and z faces (y-direction)
444 interpolate_boundary_3d<fe_degree, 1, transpose>(team_member,
445 constraint_weights,
446 constraint_mask,
447 values);
448 // Interpolate x and y faces (z-direction)
449 interpolate_boundary_3d<fe_degree, 2, transpose>(team_member,
450 constraint_weights,
451 constraint_mask,
452 values);
453 }
454 }
455 } // namespace internal
456} // namespace Portable
457
458DEAL_II_NAMESPACE_CLOSE
459#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:171
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
size
std::size_t size
Definition mpi.cc:733
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:244
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:414
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
Generated by doxygen 1.9.8