Reference documentation for deal.II version 9.2.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\}}\)
cuda_hanging_nodes_internal.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_cuda_hanging_nodes_internal_h
17 #define dealii_cuda_hanging_nodes_internal_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_COMPILER_CUDA_AWARE
22 
23 # include <deal.II/base/utilities.h>
24 
25 # include <deal.II/dofs/dof_accessor.h>
26 # include <deal.II/dofs/dof_handler.h>
27 
28 # include <deal.II/fe/fe_q.h>
29 # include <deal.II/fe/fe_tools.h>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 namespace CUDAWrappers
34 {
35  namespace internal
36  {
45  template <int dim>
46  class HangingNodes
47  {
48  public:
52  HangingNodes(unsigned int fe_degree,
53  const DoFHandler<dim> & dof_handler,
54  const std::vector<unsigned int> &lexicographic_mapping);
55 
59  template <typename CellIterator>
60  void
61  setup_constraints(
62  std::vector<types::global_dof_index> & dof_indices,
63  const CellIterator & cell,
64  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
65  unsigned int & mask) const;
66 
67  private:
71  void
72  setup_line_to_cell();
73 
74  void
75  rotate_subface_index(int times, unsigned int &subface_index) const;
76 
77  void
78  rotate_face(int times,
79  unsigned int n_dofs_1d,
80  std::vector<types::global_dof_index> &dofs) const;
81 
82  unsigned int
83  line_dof_idx(int local_line,
84  unsigned int dof,
85  unsigned int n_dofs_1d) const;
86 
87  void
88  transpose_face(std::vector<types::global_dof_index> &dofs) const;
89 
90  void
91  transpose_subface_index(unsigned int &subface) const;
92 
93  using cell_iterator = typename DoFHandler<dim>::cell_iterator;
94  using active_cell_iterator =
95  typename DoFHandler<dim>::active_cell_iterator;
96  const unsigned int n_raw_lines;
97  std::vector<std::vector<std::pair<cell_iterator, unsigned int>>>
98  line_to_cells;
99  const std::vector<unsigned int> &lexicographic_mapping;
100  const unsigned int fe_degree;
101  const DoFHandler<dim> & dof_handler;
102  };
103 
104  namespace internal
105  {
106  // TODO: use a template parameter instead of a macro
107 # define DEAL_II_MAX_ELEM_DEGREE 10
108  __constant__ double constraint_weights[(DEAL_II_MAX_ELEM_DEGREE + 1) *
109  (DEAL_II_MAX_ELEM_DEGREE + 1)];
110 
111  // Here is the system for how we store constraint types in a binary mask.
112  // This is not a complete contradiction-free system, i.e., there are
113  // invalid states that we just assume that we never get.
114 
115  // If the mask is zero, there are no constraints. Then, there are three
116  // different fields with one bit per dimension. The first field determines
117  // the type, or the position of an element along each direction. The
118  // second field determines if there is a constrained face with that
119  // direction as normal. The last field determines if there is a
120  // constrained edge of a given pair of coordinate planes, but where
121  // neither of the corresponding faces are constrained (only valid in 3D).
122 
123  // The element is placed in the 'first position' along *-axis. These also
124  // determine which face is constrained. For example, in 2D, if
125  // constr_face_x and constr_type are set, then x = 0 is constrained.
126  constexpr unsigned int constr_type_x = 1 << 0;
127  constexpr unsigned int constr_type_y = 1 << 1;
128  constexpr unsigned int constr_type_z = 1 << 2;
129 
130  // Element has as a constraint at * = 0 or * = fe_degree face
131  constexpr unsigned int constr_face_x = 1 << 3;
132  constexpr unsigned int constr_face_y = 1 << 4;
133  constexpr unsigned int constr_face_z = 1 << 5;
134 
135  // Element has as a constraint at * = 0 or * = fe_degree edge
136  constexpr unsigned int constr_edge_xy = 1 << 6;
137  constexpr unsigned int constr_edge_yz = 1 << 7;
138  constexpr unsigned int constr_edge_zx = 1 << 8;
139 
140  template <int dim>
141  void
142  setup_constraint_weigths(unsigned int fe_degree)
143  {
144  FE_Q<2> fe_q(fe_degree);
145  FullMatrix<double> interpolation_matrix(fe_q.dofs_per_face,
146  fe_q.dofs_per_face);
147  fe_q.get_subface_interpolation_matrix(fe_q, 0, interpolation_matrix);
148 
149  std::vector<unsigned int> mapping =
150  FETools::lexicographic_to_hierarchic_numbering<1>(fe_degree);
151 
152  FullMatrix<double> mapped_matrix(fe_q.dofs_per_face,
153  fe_q.dofs_per_face);
154  for (unsigned int i = 0; i < fe_q.dofs_per_face; ++i)
155  for (unsigned int j = 0; j < fe_q.dofs_per_face; ++j)
156  mapped_matrix(i, j) = interpolation_matrix(mapping[i], mapping[j]);
157 
158  cudaError_t error_code =
159  cudaMemcpyToSymbol(internal::constraint_weights,
160  &mapped_matrix[0][0],
161  sizeof(double) * fe_q.dofs_per_face *
162  fe_q.dofs_per_face);
163  AssertCuda(error_code);
164  }
165  } // namespace internal
166 
167 
168 
169  template <int dim>
170  inline HangingNodes<dim>::HangingNodes(
171  unsigned int fe_degree,
172  const DoFHandler<dim> & dof_handler,
173  const std::vector<unsigned int> &lexicographic_mapping)
174  : n_raw_lines(dof_handler.get_triangulation().n_raw_lines())
175  , line_to_cells(dim == 3 ? n_raw_lines : 0)
176  , lexicographic_mapping(lexicographic_mapping)
177  , fe_degree(fe_degree)
178  , dof_handler(dof_handler)
179  {
180  // Set up line-to-cell mapping for edge constraints (only if dim = 3)
181  setup_line_to_cell();
182 
183  internal::setup_constraint_weigths<dim>(fe_degree);
184  }
185 
186 
187 
188  template <int dim>
189  inline void
190  HangingNodes<dim>::setup_line_to_cell()
191  {}
192 
193 
194 
195  template <>
196  inline void
197  HangingNodes<3>::setup_line_to_cell()
198  {
199  // In 3D, we can have DoFs on only an edge being constrained (e.g. in a
200  // cartesian 2x2x2 grid, where only the upper left 2 cells are refined).
201  // This sets up a helper data structure in the form of a mapping from
202  // edges (i.e. lines) to neighboring cells.
203 
204  // Mapping from an edge to which children that share that edge.
205  const unsigned int line_to_children[12][2] = {{0, 2},
206  {1, 3},
207  {0, 1},
208  {2, 3},
209  {4, 6},
210  {5, 7},
211  {4, 5},
212  {6, 7},
213  {0, 4},
214  {1, 5},
215  {2, 6},
216  {3, 7}};
217 
218  std::vector<std::vector<std::pair<cell_iterator, unsigned int>>>
219  line_to_inactive_cells(n_raw_lines);
220 
221  // First add active and inactive cells to their lines:
222  for (const auto &cell : dof_handler.cell_iterators())
223  {
224  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
225  ++line)
226  {
227  const unsigned int line_idx = cell->line(line)->index();
228  if (cell->is_active())
229  line_to_cells[line_idx].push_back(std::make_pair(cell, line));
230  else
231  line_to_inactive_cells[line_idx].push_back(
232  std::make_pair(cell, line));
233  }
234  }
235 
236  // Now, we can access edge-neighboring active cells on same level to also
237  // access of an edge to the edges "children". These are found from looking
238  // at the corresponding edge of children of inactive edge neighbors.
239  for (unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
240  {
241  if ((line_to_cells[line_idx].size() > 0) &&
242  line_to_inactive_cells[line_idx].size() > 0)
243  {
244  // We now have cells to add (active ones) and edges to which they
245  // should be added (inactive cells).
246  const cell_iterator &inactive_cell =
247  line_to_inactive_cells[line_idx][0].first;
248  const unsigned int neighbor_line =
249  line_to_inactive_cells[line_idx][0].second;
250 
251  for (unsigned int c = 0; c < 2; ++c)
252  {
253  const cell_iterator &child =
254  inactive_cell->child(line_to_children[neighbor_line][c]);
255  const unsigned int child_line_idx =
256  child->line(neighbor_line)->index();
257 
258  // Now add all active cells
259  for (const auto cl : line_to_cells[line_idx])
260  line_to_cells[child_line_idx].push_back(cl);
261  }
262  }
263  }
264  }
265 
266 
267 
268  template <int dim>
269  template <typename CellIterator>
270  inline void
271  HangingNodes<dim>::setup_constraints(
272  std::vector<types::global_dof_index> & dof_indices,
273  const CellIterator & cell,
274  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
275  unsigned int & mask) const
276  {
277  mask = 0;
278  const unsigned int n_dofs_1d = fe_degree + 1;
279  const unsigned int dofs_per_face =
280  Utilities::fixed_power<dim - 1>(n_dofs_1d);
281 
282  std::vector<types::global_dof_index> neighbor_dofs(dofs_per_face);
283 
284  const auto lex_face_mapping =
285  FETools::lexicographic_to_hierarchic_numbering<dim - 1>(fe_degree);
286 
287  for (const unsigned int face : GeometryInfo<dim>::face_indices())
288  {
289  if ((!cell->at_boundary(face)) &&
290  (cell->neighbor(face)->has_children() == false))
291  {
292  const active_cell_iterator &neighbor = cell->neighbor(face);
293 
294  // Neighbor is coarser than us, i.e., face is constrained
295  if (neighbor->level() < cell->level())
296  {
297  const unsigned int neighbor_face =
298  cell->neighbor_face_no(face);
299 
300  // Find position of face on neighbor
301  unsigned int subface = 0;
302  for (; subface < GeometryInfo<dim>::max_children_per_face;
303  ++subface)
304  if (neighbor->neighbor_child_on_subface(neighbor_face,
305  subface) == cell)
306  break;
307 
308  // Get indices to read
309  neighbor->face(neighbor_face)->get_dof_indices(neighbor_dofs);
310  // If the vector is distributed, we need to transform the
311  // global indices to local ones.
312  if (partitioner)
313  for (auto &index : neighbor_dofs)
314  index = partitioner->global_to_local(index);
315 
316  if (dim == 2)
317  {
318  if (face < 2)
319  {
320  mask |= internal::constr_face_x;
321  if (face == 0)
322  mask |= internal::constr_type_x;
323  if (subface == 0)
324  mask |= internal::constr_type_y;
325  }
326  else
327  {
328  mask |= internal::constr_face_y;
329  if (face == 2)
330  mask |= internal::constr_type_y;
331  if (subface == 0)
332  mask |= internal::constr_type_x;
333  }
334 
335  // Reorder neighbor_dofs and copy into faceth face of
336  // dof_indices
337 
338  // Offset if upper/right face
339  unsigned int offset = (face % 2 == 1) ? fe_degree : 0;
340 
341  for (unsigned int i = 0; i < n_dofs_1d; ++i)
342  {
343  unsigned int idx = 0;
344  // If X-line, i.e., if y = 0 or y = fe_degree
345  if (face > 1)
346  idx = n_dofs_1d * offset + i;
347  // If Y-line, i.e., if x = 0 or x = fe_degree
348  else
349  idx = n_dofs_1d * i + offset;
350 
351  dof_indices[idx] = neighbor_dofs[lex_face_mapping[i]];
352  }
353  }
354  else if (dim == 3)
355  {
356  const bool transpose = !(cell->face_orientation(face));
357 
358  int rotate = 0;
359 
360  if (cell->face_rotation(face))
361  rotate -= 1;
362  if (cell->face_flip(face))
363  rotate -= 2;
364 
365  rotate_face(rotate, n_dofs_1d, neighbor_dofs);
366  rotate_subface_index(rotate, subface);
367 
368  if (transpose)
369  {
370  transpose_face(neighbor_dofs);
371  transpose_subface_index(subface);
372  }
373 
374  // YZ-plane
375  if (face < 2)
376  {
377  mask |= internal::constr_face_x;
378  if (face == 0)
379  mask |= internal::constr_type_x;
380  if (subface % 2 == 0)
381  mask |= internal::constr_type_y;
382  if (subface / 2 == 0)
383  mask |= internal::constr_type_z;
384  }
385  // XZ-plane
386  else if (face < 4)
387  {
388  mask |= internal::constr_face_y;
389  if (face == 2)
390  mask |= internal::constr_type_y;
391  if (subface % 2 == 0)
392  mask |= internal::constr_type_z;
393  if (subface / 2 == 0)
394  mask |= internal::constr_type_x;
395  }
396  // XY-plane
397  else
398  {
399  mask |= internal::constr_face_z;
400  if (face == 4)
401  mask |= internal::constr_type_z;
402  if (subface % 2 == 0)
403  mask |= internal::constr_type_x;
404  if (subface / 2 == 0)
405  mask |= internal::constr_type_y;
406  }
407 
408  // Offset if upper/right/back face
409  unsigned int offset = (face % 2 == 1) ? fe_degree : 0;
410 
411  for (unsigned int i = 0; i < n_dofs_1d; ++i)
412  {
413  for (unsigned int j = 0; j < n_dofs_1d; ++j)
414  {
415  unsigned int idx = 0;
416  // If YZ-plane, i.e., if x = 0 or x = fe_degree,
417  // and orientation standard
418  if (face < 2)
419  idx = n_dofs_1d * n_dofs_1d * i +
420  n_dofs_1d * j + offset;
421  // If XZ-plane, i.e., if y = 0 or y = fe_degree,
422  // and orientation standard
423  else if (face < 4)
424  idx = n_dofs_1d * n_dofs_1d * j +
425  n_dofs_1d * offset + i;
426  // If XY-plane, i.e., if z = 0 or z = fe_degree,
427  // and orientation standard
428  else
429  idx = n_dofs_1d * n_dofs_1d * offset +
430  n_dofs_1d * i + j;
431 
432  dof_indices[idx] =
433  neighbor_dofs[lex_face_mapping[n_dofs_1d * i +
434  j]];
435  }
436  }
437  }
438  else
439  ExcNotImplemented();
440  }
441  }
442  }
443 
444  // In 3D we can have a situation where only DoFs on an edge are
445  // constrained. Append these here.
446  if (dim == 3)
447  {
448  // For each line on cell, which faces does it belong to, what is the
449  // edge mask, what is the types of the faces it belong to, and what is
450  // the type along the edge.
451  const unsigned int line_to_edge[12][4] = {
452  {internal::constr_face_x | internal::constr_face_z,
453  internal::constr_edge_zx,
454  internal::constr_type_x | internal::constr_type_z,
455  internal::constr_type_y},
456  {internal::constr_face_x | internal::constr_face_z,
457  internal::constr_edge_zx,
458  internal::constr_type_z,
459  internal::constr_type_y},
460  {internal::constr_face_y | internal::constr_face_z,
461  internal::constr_edge_yz,
462  internal::constr_type_y | internal::constr_type_z,
463  internal::constr_type_x},
464  {internal::constr_face_y | internal::constr_face_z,
465  internal::constr_edge_yz,
466  internal::constr_type_z,
467  internal::constr_type_x},
468  {internal::constr_face_x | internal::constr_face_z,
469  internal::constr_edge_zx,
470  internal::constr_type_x,
471  internal::constr_type_y},
472  {internal::constr_face_x | internal::constr_face_z,
473  internal::constr_edge_zx,
474  0,
475  internal::constr_type_y},
476  {internal::constr_face_y | internal::constr_face_z,
477  internal::constr_edge_yz,
478  internal::constr_type_y,
479  internal::constr_type_x},
480  {internal::constr_face_y | internal::constr_face_z,
481  internal::constr_edge_yz,
482  0,
483  internal::constr_type_x},
484  {internal::constr_face_x | internal::constr_face_y,
485  internal::constr_edge_xy,
486  internal::constr_type_x | internal::constr_type_y,
487  internal::constr_type_z},
488  {internal::constr_face_x | internal::constr_face_y,
489  internal::constr_edge_xy,
490  internal::constr_type_y,
491  internal::constr_type_z},
492  {internal::constr_face_x | internal::constr_face_y,
493  internal::constr_edge_xy,
494  internal::constr_type_x,
495  internal::constr_type_z},
496  {internal::constr_face_x | internal::constr_face_y,
497  internal::constr_edge_xy,
498  0,
499  internal::constr_type_z}};
500 
501  for (unsigned int local_line = 0;
502  local_line < GeometryInfo<dim>::lines_per_cell;
503  ++local_line)
504  {
505  // If we don't already have a constraint for as part of a face
506  if (!(mask & line_to_edge[local_line][0]))
507  {
508  // For each cell which share that edge
509  const unsigned int line = cell->line(local_line)->index();
510  for (const auto edge_neighbor : line_to_cells[line])
511  {
512  // If one of them is coarser than us
513  const cell_iterator neighbor_cell = edge_neighbor.first;
514  if (neighbor_cell->level() < cell->level())
515  {
516  const unsigned int local_line_neighbor =
517  edge_neighbor.second;
518  mask |= line_to_edge[local_line][1] |
519  line_to_edge[local_line][2];
520 
521  bool flipped = false;
522  if (cell->line(local_line)->vertex_index(0) ==
523  neighbor_cell->line(local_line_neighbor)
524  ->vertex_index(0))
525  {
526  // Assuming line directions match axes directions,
527  // we have an unflipped edge of first type
528  mask |= line_to_edge[local_line][3];
529  }
530  else if (cell->line(local_line)->vertex_index(1) ==
531  neighbor_cell->line(local_line_neighbor)
532  ->vertex_index(1))
533  {
534  // We have an unflipped edge of second type
535  }
536  else if (cell->line(local_line)->vertex_index(1) ==
537  neighbor_cell->line(local_line_neighbor)
538  ->vertex_index(0))
539  {
540  // We have a flipped edge of second type
541  flipped = true;
542  }
543  else if (cell->line(local_line)->vertex_index(0) ==
544  neighbor_cell->line(local_line_neighbor)
545  ->vertex_index(1))
546  {
547  // We have a flipped edge of first type
548  mask |= line_to_edge[local_line][3];
549  flipped = true;
550  }
551  else
552  ExcInternalError();
553 
554  // Copy the unconstrained values
555  neighbor_dofs.resize(n_dofs_1d * n_dofs_1d *
556  n_dofs_1d);
557  neighbor_cell->get_dof_indices(neighbor_dofs);
558  // If the vector is distributed, we need to transform
559  // the global indices to local ones.
560  if (partitioner)
561  for (auto &index : neighbor_dofs)
562  index = partitioner->global_to_local(index);
563 
564  for (unsigned int i = 0; i < n_dofs_1d; ++i)
565  {
566  // Get local dof index along line
567  const unsigned int idx =
568  line_dof_idx(local_line, i, n_dofs_1d);
569  dof_indices[idx] = neighbor_dofs
570  [lexicographic_mapping[line_dof_idx(
571  local_line_neighbor,
572  flipped ? fe_degree - i : i,
573  n_dofs_1d)]];
574  }
575 
576  // Stop looping over edge neighbors
577  break;
578  }
579  }
580  }
581  }
582  }
583  }
584 
585 
586 
587  template <int dim>
588  inline void
589  HangingNodes<dim>::rotate_subface_index(int times,
590  unsigned int &subface_index) const
591  {
592  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
593 
594  times = times % 4;
595  times = times < 0 ? times + 4 : times;
596  for (int t = 0; t < times; ++t)
597  subface_index = rot_mapping[subface_index];
598  }
599 
600 
601 
602  template <int dim>
603  inline void
604  HangingNodes<dim>::rotate_face(
605  int times,
606  unsigned int n_dofs_1d,
607  std::vector<types::global_dof_index> &dofs) const
608  {
609  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
610 
611  times = times % 4;
612  times = times < 0 ? times + 4 : times;
613 
614  std::vector<types::global_dof_index> copy(dofs.size());
615  for (int t = 0; t < times; ++t)
616  {
617  std::swap(copy, dofs);
618 
619  // Vertices
620  for (unsigned int i = 0; i < 4; ++i)
621  dofs[rot_mapping[i]] = copy[i];
622 
623  // Edges
624  const unsigned int n_int = n_dofs_1d - 2;
625  unsigned int offset = 4;
626  for (unsigned int i = 0; i < n_int; ++i)
627  {
628  // Left edge
629  dofs[offset + i] = copy[offset + 2 * n_int + (n_int - 1 - i)];
630  // Right edge
631  dofs[offset + n_int + i] =
632  copy[offset + 3 * n_int + (n_int - 1 - i)];
633  // Bottom edge
634  dofs[offset + 2 * n_int + i] = copy[offset + n_int + i];
635  // Top edge
636  dofs[offset + 3 * n_int + i] = copy[offset + i];
637  }
638 
639  // Interior points
640  offset += 4 * n_int;
641 
642  for (unsigned int i = 0; i < n_int; ++i)
643  for (unsigned int j = 0; j < n_int; ++j)
644  dofs[offset + i * n_int + j] =
645  copy[offset + j * n_int + (n_int - 1 - i)];
646  }
647  }
648 
649 
650 
651  template <int dim>
652  inline unsigned int
653  HangingNodes<dim>::line_dof_idx(int local_line,
654  unsigned int dof,
655  unsigned int n_dofs_1d) const
656  {
657  unsigned int x, y, z;
658 
659  if (local_line < 8)
660  {
661  x =
662  (local_line % 4 == 0) ? 0 : (local_line % 4 == 1) ? fe_degree : dof;
663  y =
664  (local_line % 4 == 2) ? 0 : (local_line % 4 == 3) ? fe_degree : dof;
665  z = (local_line / 4) * fe_degree;
666  }
667  else
668  {
669  x = ((local_line - 8) % 2) * fe_degree;
670  y = ((local_line - 8) / 2) * fe_degree;
671  z = dof;
672  }
673 
674  return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
675  }
676 
677 
678 
679  template <int dim>
680  inline void
681  HangingNodes<dim>::transpose_face(
682  std::vector<types::global_dof_index> &dofs) const
683  {
684  const std::vector<types::global_dof_index> copy(dofs);
685 
686  // Vertices
687  dofs[1] = copy[2];
688  dofs[2] = copy[1];
689 
690  // Edges
691  const unsigned int n_int = fe_degree - 1;
692  unsigned int offset = 4;
693  for (unsigned int i = 0; i < n_int; ++i)
694  {
695  // Right edge
696  dofs[offset + i] = copy[offset + 2 * n_int + i];
697  // Left edge
698  dofs[offset + n_int + i] = copy[offset + 3 * n_int + i];
699  // Bottom edge
700  dofs[offset + 2 * n_int + i] = copy[offset + i];
701  // Top edge
702  dofs[offset + 3 * n_int + i] = copy[offset + n_int + i];
703  }
704 
705  // Interior
706  offset += 4 * n_int;
707  for (unsigned int i = 0; i < n_int; ++i)
708  for (unsigned int j = 0; j < n_int; ++j)
709  dofs[offset + i * n_int + j] = copy[offset + j * n_int + i];
710  }
711 
712 
713 
714  template <int dim>
715  void
716  HangingNodes<dim>::transpose_subface_index(unsigned int &subface) const
717  {
718  if (subface == 1)
719  subface = 2;
720  else if (subface == 2)
721  subface = 1;
722  }
723 
724 
725 
726  //-------------------------------------------------------------------------//
727  // Functions for resolving the hanging node constraints //
728  //-------------------------------------------------------------------------//
729  namespace internal
730  {
731  template <unsigned int size>
732  __device__ inline unsigned int
733  index2(unsigned int i, unsigned int j)
734  {
735  return i + size * j;
736  }
737 
738 
739 
740  template <unsigned int size>
741  __device__ inline unsigned int
742  index3(unsigned int i, unsigned int j, unsigned int k)
743  {
744  return i + size * j + size * size * k;
745  }
746 
747 
748 
749  template <unsigned int fe_degree,
750  unsigned int direction,
751  bool transpose,
752  typename Number>
753  __device__ inline void
754  interpolate_boundary_2d(const unsigned int constraint_mask,
755  Number * values)
756  {
757  const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
758  const unsigned int y_idx = threadIdx.y;
759 
760  const unsigned int this_type =
761  (direction == 0) ? internal::constr_type_x : internal::constr_type_y;
762 
763  const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
764 
765  Number t = 0;
766  // Flag is true if dof is constrained for the given direction and the
767  // given face.
768  const bool constrained_face =
769  (constraint_mask &
770  (((direction == 0) ? internal::constr_face_y : 0) |
771  ((direction == 1) ? internal::constr_face_x : 0)));
772 
773  // Flag is true if for the given direction, the dof is constrained with
774  // the right type and is on the correct side (left (= 0) or right (=
775  // fe_degree))
776  const bool constrained_dof =
777  ((direction == 0) && ((constraint_mask & internal::constr_type_y) ?
778  (y_idx == 0) :
779  (y_idx == fe_degree))) ||
780  ((direction == 1) && ((constraint_mask & internal::constr_type_x) ?
781  (x_idx == 0) :
782  (x_idx == fe_degree)));
783 
784  if (constrained_face && constrained_dof)
785  {
786  const bool type = constraint_mask & this_type;
787 
788  if (type)
789  {
790  for (unsigned int i = 0; i <= fe_degree; ++i)
791  {
792  const unsigned int real_idx =
793  (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
794  index2<fe_degree + 1>(x_idx, i);
795 
796  const Number w =
797  transpose ?
798  internal::constraint_weights[i * (fe_degree + 1) +
799  interp_idx] :
800  internal::constraint_weights[interp_idx *
801  (fe_degree + 1) +
802  i];
803  t += w * values[real_idx];
804  }
805  }
806  else
807  {
808  for (unsigned int i = 0; i <= fe_degree; ++i)
809  {
810  const unsigned int real_idx =
811  (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
812  index2<fe_degree + 1>(x_idx, i);
813 
814  const Number w =
815  transpose ?
816  internal::constraint_weights[(fe_degree - i) *
817  (fe_degree + 1) +
818  fe_degree - interp_idx] :
819  internal::constraint_weights[(fe_degree - interp_idx) *
820  (fe_degree + 1) +
821  fe_degree - i];
822  t += w * values[real_idx];
823  }
824  }
825  }
826 
827  // The synchronization is done for all the threads in one block with
828  // each block being assigned to one element.
829  __syncthreads();
830  if (constrained_face && constrained_dof)
831  values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
832 
833  __syncthreads();
834  }
835 
836 
837 
838  template <unsigned int fe_degree,
839  unsigned int direction,
840  bool transpose,
841  typename Number>
842  __device__ inline void
843  interpolate_boundary_3d(const unsigned int constraint_mask,
844  Number * values)
845  {
846  const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
847  const unsigned int y_idx = threadIdx.y;
848  const unsigned int z_idx = threadIdx.z;
849 
850  const unsigned int this_type =
851  (direction == 0) ? internal::constr_type_x :
852  (direction == 1) ? internal::constr_type_y :
853  internal::constr_type_z;
854  const unsigned int face1_type =
855  (direction == 0) ? internal::constr_type_y :
856  (direction == 1) ? internal::constr_type_z :
857  internal::constr_type_x;
858  const unsigned int face2_type =
859  (direction == 0) ? internal::constr_type_z :
860  (direction == 1) ? internal::constr_type_x :
861  internal::constr_type_y;
862 
863  // If computing in x-direction, need to match against constr_face_y or
864  // constr_face_z
865  const unsigned int face1 = (direction == 0) ? internal::constr_face_y :
866  (direction == 1) ?
867  internal::constr_face_z :
868  internal::constr_face_x;
869  const unsigned int face2 = (direction == 0) ? internal::constr_face_z :
870  (direction == 1) ?
871  internal::constr_face_x :
872  internal::constr_face_y;
873  const unsigned int edge = (direction == 0) ? internal::constr_edge_yz :
874  (direction == 1) ?
875  internal::constr_edge_zx :
876  internal::constr_edge_xy;
877  const unsigned int constrained_face =
878  constraint_mask & (face1 | face2 | edge);
879 
880  const unsigned int interp_idx =
881  (direction == 0) ? x_idx : (direction == 1) ? y_idx : z_idx;
882  const unsigned int face1_idx =
883  (direction == 0) ? y_idx : (direction == 1) ? z_idx : x_idx;
884  const unsigned int face2_idx =
885  (direction == 0) ? z_idx : (direction == 1) ? x_idx : y_idx;
886 
887  Number t = 0;
888  const bool on_face1 = (constraint_mask & face1_type) ?
889  (face1_idx == 0) :
890  (face1_idx == fe_degree);
891  const bool on_face2 = (constraint_mask & face2_type) ?
892  (face2_idx == 0) :
893  (face2_idx == fe_degree);
894  const bool constrained_dof =
895  (((constraint_mask & face1) && on_face1) ||
896  ((constraint_mask & face2) && on_face2) ||
897  ((constraint_mask & edge) && on_face1 && on_face2));
898 
899  if (constrained_face && constrained_dof)
900  {
901  const bool type = constraint_mask & this_type;
902  if (type)
903  {
904  for (unsigned int i = 0; i <= fe_degree; ++i)
905  {
906  const unsigned int real_idx =
907  (direction == 0) ?
908  index3<fe_degree + 1>(i, y_idx, z_idx) :
909  (direction == 1) ?
910  index3<fe_degree + 1>(x_idx, i, z_idx) :
911  index3<fe_degree + 1>(x_idx, y_idx, i);
912 
913  const Number w =
914  transpose ?
915  internal::constraint_weights[i * (fe_degree + 1) +
916  interp_idx] :
917  internal::constraint_weights[interp_idx *
918  (fe_degree + 1) +
919  i];
920  t += w * values[real_idx];
921  }
922  }
923  else
924  {
925  for (unsigned int i = 0; i <= fe_degree; ++i)
926  {
927  const unsigned int real_idx =
928  (direction == 0) ?
929  index3<fe_degree + 1>(i, y_idx, z_idx) :
930  (direction == 1) ?
931  index3<fe_degree + 1>(x_idx, i, z_idx) :
932  index3<fe_degree + 1>(x_idx, y_idx, i);
933 
934  const Number w =
935  transpose ?
936  internal::constraint_weights[(fe_degree - i) *
937  (fe_degree + 1) +
938  fe_degree - interp_idx] :
939  internal::constraint_weights[(fe_degree - interp_idx) *
940  (fe_degree + 1) +
941  fe_degree - i];
942  t += w * values[real_idx];
943  }
944  }
945  }
946 
947  // The synchronization is done for all the threads in one block with
948  // each block being assigned to one element.
949  __syncthreads();
950 
951  if (constrained_face && constrained_dof)
952  values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
953 
954  __syncthreads();
955  }
956  } // namespace internal
957 
958 
959 
968  template <int dim, int fe_degree, bool transpose, typename Number>
969  __device__ void
970  resolve_hanging_nodes(const unsigned int constraint_mask, Number *values)
971  {
972  if (dim == 2)
973  {
974  internal::interpolate_boundary_2d<fe_degree, 0, transpose>(
975  constraint_mask, values);
976 
977  internal::interpolate_boundary_2d<fe_degree, 1, transpose>(
978  constraint_mask, values);
979  }
980  else if (dim == 3)
981  {
982  // Interpolate y and z faces (x-direction)
983  internal::interpolate_boundary_3d<fe_degree, 0, transpose>(
984  constraint_mask, values);
985  // Interpolate x and z faces (y-direction)
986  internal::interpolate_boundary_3d<fe_degree, 1, transpose>(
987  constraint_mask, values);
988  // Interpolate x and y faces (z-direction)
989  internal::interpolate_boundary_3d<fe_degree, 2, transpose>(
990  constraint_mask, values);
991  }
992  }
993  } // namespace internal
994 } // namespace CUDAWrappers
995 
996 DEAL_II_NAMESPACE_CLOSE
997 
998 #endif
999 
1000 #endif
CUDAWrappers::internal::HangingNodes::HangingNodes
HangingNodes(unsigned int fe_degree, const DoFHandler< dim > &dof_handler, const std::vector< unsigned int > &lexicographic_mapping)
Definition: cuda_hanging_nodes_internal.h:170
CUDAWrappers::internal::internal::constr_type_x
constexpr unsigned int constr_type_x
Definition: cuda_hanging_nodes_internal.h:126
DEAL_II_MAX_ELEM_DEGREE
#define DEAL_II_MAX_ELEM_DEGREE
Definition: cuda_hanging_nodes_internal.h:107
CUDAWrappers::internal::HangingNodes::active_cell_iterator
typename DoFHandler< dim >::active_cell_iterator active_cell_iterator
Definition: cuda_hanging_nodes_internal.h:95
CUDAWrappers::internal::internal::index2
unsigned int index2(unsigned int i, unsigned int j)
Definition: cuda_hanging_nodes_internal.h:733
CUDAWrappers::internal::internal::constr_edge_yz
constexpr unsigned int constr_edge_yz
Definition: cuda_hanging_nodes_internal.h:137
CUDAWrappers::internal::HangingNodes::rotate_face
void rotate_face(int times, unsigned int n_dofs_1d, std::vector< types::global_dof_index > &dofs) const
Definition: cuda_hanging_nodes_internal.h:604
FE_Q
Definition: fe_q.h:554
CUDAWrappers::internal::HangingNodes::setup_line_to_cell
void setup_line_to_cell()
Definition: cuda_hanging_nodes_internal.h:190
CUDAWrappers::internal::internal::constr_type_y
constexpr unsigned int constr_type_y
Definition: cuda_hanging_nodes_internal.h:127
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
utilities.h
CUDAWrappers::internal::HangingNodes::n_raw_lines
const unsigned int n_raw_lines
Definition: cuda_hanging_nodes_internal.h:96
CUDAWrappers::internal::internal::constr_face_z
constexpr unsigned int constr_face_z
Definition: cuda_hanging_nodes_internal.h:133
GeometryInfo
Definition: geometry_info.h:1224
CUDAWrappers::internal::internal::constr_face_y
constexpr unsigned int constr_face_y
Definition: cuda_hanging_nodes_internal.h:132
CUDAWrappers::internal::HangingNodes::line_to_cells
std::vector< std::vector< std::pair< cell_iterator, unsigned int > > > line_to_cells
Definition: cuda_hanging_nodes_internal.h:98
DoFHandler< dim >
AssertCuda
#define AssertCuda(error_code)
Definition: exceptions.h:1734
CUDAWrappers::internal::HangingNodes::setup_constraints
void setup_constraints(std::vector< types::global_dof_index > &dof_indices, const CellIterator &cell, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner, unsigned int &mask) const
Definition: cuda_hanging_nodes_internal.h:271
Physics::Elasticity::Kinematics::w
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
CUDAWrappers::internal::internal::interpolate_boundary_2d
void interpolate_boundary_2d(const unsigned int constraint_mask, Number *values)
Definition: cuda_hanging_nodes_internal.h:754
CUDAWrappers::internal::HangingNodes::dof_handler
const DoFHandler< dim > & dof_handler
Definition: cuda_hanging_nodes_internal.h:101
CUDAWrappers
Definition: cuda_size.h:23
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1072
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
CUDAWrappers::internal::internal::index3
unsigned int index3(unsigned int i, unsigned int j, unsigned int k)
Definition: cuda_hanging_nodes_internal.h:742
CUDAWrappers::internal::internal::constraint_weights
__constant__ double constraint_weights[(DEAL_II_MAX_ELEM_DEGREE+1) *(DEAL_II_MAX_ELEM_DEGREE+1)]
Definition: cuda_hanging_nodes_internal.h:109
MemorySpace::swap
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
Definition: memory_space.h:103
fe_q.h
CUDAWrappers::internal::HangingNodes::line_dof_idx
unsigned int line_dof_idx(int local_line, unsigned int dof, unsigned int n_dofs_1d) const
Definition: cuda_hanging_nodes_internal.h:653
Utilities::MPI::Partitioner::global_to_local
unsigned int global_to_local(const types::global_dof_index global_index) const
CUDAWrappers::internal::internal::setup_constraint_weigths
void setup_constraint_weigths(unsigned int fe_degree)
Definition: cuda_hanging_nodes_internal.h:142
FE_Q_Base< TensorProductPolynomials< dim >, dim, dim >::get_subface_interpolation_matrix
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
Definition: fe_q_base.cc:583
fe_tools.h
FiniteElementData::dofs_per_face
const unsigned int dofs_per_face
Definition: fe_base.h:275
CUDAWrappers::internal::HangingNodes::transpose_face
void transpose_face(std::vector< types::global_dof_index > &dofs) const
Definition: cuda_hanging_nodes_internal.h:681
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
CUDAWrappers::internal::HangingNodes::cell_iterator
typename DoFHandler< dim >::cell_iterator cell_iterator
Definition: cuda_hanging_nodes_internal.h:93
CUDAWrappers::internal::internal::constr_edge_zx
constexpr unsigned int constr_edge_zx
Definition: cuda_hanging_nodes_internal.h:138
CUDAWrappers::internal::HangingNodes::fe_degree
const unsigned int fe_degree
Definition: cuda_hanging_nodes_internal.h:100
dof_handler.h
CUDAWrappers::internal::internal::constr_edge_xy
constexpr unsigned int constr_edge_xy
Definition: cuda_hanging_nodes_internal.h:136
dof_accessor.h
CUDAWrappers::internal::HangingNodes::lexicographic_mapping
const std::vector< unsigned int > & lexicographic_mapping
Definition: cuda_hanging_nodes_internal.h:99
CUDAWrappers::internal::HangingNodes::rotate_subface_index
void rotate_subface_index(int times, unsigned int &subface_index) const
Definition: cuda_hanging_nodes_internal.h:589
config.h
FullMatrix< double >
internal
Definition: aligned_vector.h:369
GridTools::rotate
void rotate(const double angle, Triangulation< dim > &triangulation)
CUDAWrappers::internal::HangingNodes::transpose_subface_index
void transpose_subface_index(unsigned int &subface) const
Definition: cuda_hanging_nodes_internal.h:716
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
CUDAWrappers::internal::resolve_hanging_nodes
void resolve_hanging_nodes(const unsigned int constraint_mask, Number *values)
Definition: cuda_hanging_nodes_internal.h:970
CUDAWrappers::internal::internal::constr_type_z
constexpr unsigned int constr_type_z
Definition: cuda_hanging_nodes_internal.h:128
CUDAWrappers::internal::internal::interpolate_boundary_3d
void interpolate_boundary_3d(const unsigned int constraint_mask, Number *values)
Definition: cuda_hanging_nodes_internal.h:843
CUDAWrappers::internal::internal::constr_face_x
constexpr unsigned int constr_face_x
Definition: cuda_hanging_nodes_internal.h:131
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
CUDAWrappers::internal::HangingNodes
Definition: cuda_hanging_nodes_internal.h:46
FETools::lexicographic_to_hierarchic_numbering
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
Generated by   doxygen 1.8.17