Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.1.1
\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
cuda_hanging_nodes_internal.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2019 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>
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
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
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 =
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_Q<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>
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  AssertThrow(
181  (dim == 3) || ((fe_degree % 2) == 1),
182  ExcMessage(
183  "This function is not implemented when dim = 2 and fe_degree is even."));
184 
185  // Set up line-to-cell mapping for edge constraints (only if dim = 3)
187 
188  internal::setup_constraint_weigths<dim>(fe_degree);
189  }
190 
191 
192 
193  template <int dim>
194  void
196  {}
197 
198 
199 
200  template <>
201  void
203  {
204  // In 3D, we can have DoFs on only an edge being constrained (e.g. in a
205  // cartesian 2x2x2 grid, where only the upper left 2 cells are refined).
206  // This sets up a helper data structure in the form of a mapping from
207  // edges (i.e. lines) to neighboring cells.
208 
209  // Mapping from an edge to which children that share that edge.
210  const unsigned int line_to_children[12][2] = {{0, 2},
211  {1, 3},
212  {0, 1},
213  {2, 3},
214  {4, 6},
215  {5, 7},
216  {4, 5},
217  {6, 7},
218  {0, 4},
219  {1, 5},
220  {2, 6},
221  {3, 7}};
222 
223  std::vector<std::vector<std::pair<cell_iterator, unsigned int>>>
224  line_to_inactive_cells(n_raw_lines);
225 
226  // First add active and inactive cells to their lines:
227  for (const auto &cell : dof_handler.cell_iterators())
228  {
229  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
230  ++line)
231  {
232  const unsigned int line_idx = cell->line(line)->index();
233  if (cell->active())
234  line_to_cells[line_idx].push_back(std::make_pair(cell, line));
235  else
236  line_to_inactive_cells[line_idx].push_back(
237  std::make_pair(cell, line));
238  }
239  }
240 
241  // Now, we can access edge-neighboring active cells on same level to also
242  // access of an edge to the edges "children". These are found from looking
243  // at the corresponding edge of children of inactive edge neighbors.
244  for (unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
245  {
246  if ((line_to_cells[line_idx].size() > 0) &&
247  line_to_inactive_cells[line_idx].size() > 0)
248  {
249  // We now have cells to add (active ones) and edges to which they
250  // should be added (inactive cells).
251  const cell_iterator &inactive_cell =
252  line_to_inactive_cells[line_idx][0].first;
253  const unsigned int neighbor_line =
254  line_to_inactive_cells[line_idx][0].second;
255 
256  for (unsigned int c = 0; c < 2; ++c)
257  {
258  const cell_iterator &child =
259  inactive_cell->child(line_to_children[neighbor_line][c]);
260  const unsigned int child_line_idx =
261  child->line(neighbor_line)->index();
262 
263  // Now add all active cells
264  for (const auto cl : line_to_cells[line_idx])
265  line_to_cells[child_line_idx].push_back(cl);
266  }
267  }
268  }
269  }
270 
271 
272 
273  template <int dim>
274  template <typename CellIterator>
275  void
277  std::vector<types::global_dof_index> & dof_indices,
278  const CellIterator & cell,
279  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
280  unsigned int & mask) const
281  {
282  mask = 0;
283  const unsigned int n_dofs_1d = fe_degree + 1;
284  const unsigned int dofs_per_face =
285  Utilities::fixed_power<dim - 1>(n_dofs_1d);
286 
287  std::vector<types::global_dof_index> neighbor_dofs(dofs_per_face);
288 
289  const auto lex_face_mapping =
291  FE_Q<dim - 1>(fe_degree));
292 
293  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
294  ++face)
295  {
296  if ((!cell->at_boundary(face)) &&
297  (cell->neighbor(face)->has_children() == false))
298  {
299  const active_cell_iterator &neighbor = cell->neighbor(face);
300 
301  // Neighbor is coarser than us, i.e., face is constrained
302  if (neighbor->level() < cell->level())
303  {
304  const unsigned int neighbor_face =
305  cell->neighbor_face_no(face);
306 
307  // Find position of face on neighbor
308  unsigned int subface = 0;
309  for (; subface < GeometryInfo<dim>::max_children_per_face;
310  ++subface)
311  if (neighbor->neighbor_child_on_subface(neighbor_face,
312  subface) == cell)
313  break;
314 
315  // Get indices to read
316  neighbor->face(neighbor_face)->get_dof_indices(neighbor_dofs);
317  // If the vector is distributed, we need to transform the
318  // global indices to local ones.
319  if (partitioner)
320  for (auto &index : neighbor_dofs)
321  index = partitioner->global_to_local(index);
322 
323  if (dim == 2)
324  {
325  if (face < 2)
326  {
327  mask |= internal::constr_face_x;
328  if (face == 0)
329  mask |= internal::constr_type_x;
330  if (subface == 0)
331  mask |= internal::constr_type_y;
332  }
333  else
334  {
335  mask |= internal::constr_face_y;
336  if (face == 2)
337  mask |= internal::constr_type_y;
338  if (subface == 0)
339  mask |= internal::constr_type_x;
340  }
341 
342  // Reorder neighbor_dofs and copy into faceth face of
343  // dof_indices
344 
345  // Offset if upper/right face
346  unsigned int offset = (face % 2 == 1) ? fe_degree : 0;
347 
348  for (unsigned int i = 0; i < n_dofs_1d; ++i)
349  {
350  unsigned int idx = 0;
351  // If X-line, i.e., if y = 0 or y = fe_degree
352  if (face > 1)
353  idx = n_dofs_1d * offset + i;
354  // If Y-line, i.e., if x = 0 or x = fe_degree
355  else
356  idx = n_dofs_1d * i + offset;
357 
358  dof_indices[idx] = neighbor_dofs[lex_face_mapping[i]];
359  }
360  }
361  else if (dim == 3)
362  {
363  const bool transpose = !(cell->face_orientation(face));
364 
365  int rotate = 0;
366 
367  if (cell->face_rotation(face))
368  rotate -= 1;
369  if (cell->face_flip(face))
370  rotate -= 2;
371 
372  rotate_face(rotate, n_dofs_1d, neighbor_dofs);
373  rotate_subface_index(rotate, subface);
374 
375  if (transpose)
376  {
377  transpose_face(neighbor_dofs);
378  transpose_subface_index(subface);
379  }
380 
381  // YZ-plane
382  if (face < 2)
383  {
384  mask |= internal::constr_face_x;
385  if (face == 0)
386  mask |= internal::constr_type_x;
387  if (subface % 2 == 0)
388  mask |= internal::constr_type_y;
389  if (subface / 2 == 0)
390  mask |= internal::constr_type_z;
391  }
392  // XZ-plane
393  else if (face < 4)
394  {
395  mask |= internal::constr_face_y;
396  if (face == 2)
397  mask |= internal::constr_type_y;
398  if (subface % 2 == 0)
399  mask |= internal::constr_type_z;
400  if (subface / 2 == 0)
401  mask |= internal::constr_type_x;
402  }
403  // XY-plane
404  else
405  {
406  mask |= internal::constr_face_z;
407  if (face == 4)
408  mask |= internal::constr_type_z;
409  if (subface % 2 == 0)
410  mask |= internal::constr_type_x;
411  if (subface / 2 == 0)
412  mask |= internal::constr_type_y;
413  }
414 
415  // Offset if upper/right/back face
416  unsigned int offset = (face % 2 == 1) ? fe_degree : 0;
417 
418  for (unsigned int i = 0; i < n_dofs_1d; ++i)
419  {
420  for (unsigned int j = 0; j < n_dofs_1d; ++j)
421  {
422  unsigned int idx = 0;
423  // If YZ-plane, i.e., if x = 0 or x = fe_degree,
424  // and orientation standard
425  if (face < 2)
426  idx = n_dofs_1d * n_dofs_1d * i +
427  n_dofs_1d * j + offset;
428  // If XZ-plane, i.e., if y = 0 or y = fe_degree,
429  // and orientation standard
430  else if (face < 4)
431  idx = n_dofs_1d * n_dofs_1d * j +
432  n_dofs_1d * offset + i;
433  // If XY-plane, i.e., if z = 0 or z = fe_degree,
434  // and orientation standard
435  else
436  idx = n_dofs_1d * n_dofs_1d * offset +
437  n_dofs_1d * i + j;
438 
439  dof_indices[idx] =
440  neighbor_dofs[lex_face_mapping[n_dofs_1d * i +
441  j]];
442  }
443  }
444  }
445  else
447  }
448  }
449  }
450 
451  // In 3D we can have a situation where only DoFs on an edge are
452  // constrained. Append these here.
453  if (dim == 3)
454  {
455  // For each line on cell, which faces does it belong to, what is the
456  // edge mask, what is the types of the faces it belong to, and what is
457  // the type along the edge.
458  const unsigned int line_to_edge[12][4] = {
459  {internal::constr_face_x | internal::constr_face_z,
460  internal::constr_edge_zx,
461  internal::constr_type_x | internal::constr_type_z,
462  internal::constr_type_y},
463  {internal::constr_face_x | internal::constr_face_z,
464  internal::constr_edge_zx,
465  internal::constr_type_z,
466  internal::constr_type_y},
467  {internal::constr_face_y | internal::constr_face_z,
468  internal::constr_edge_yz,
469  internal::constr_type_y | internal::constr_type_z,
470  internal::constr_type_x},
471  {internal::constr_face_y | internal::constr_face_z,
472  internal::constr_edge_yz,
473  internal::constr_type_z,
474  internal::constr_type_x},
475  {internal::constr_face_x | internal::constr_face_z,
476  internal::constr_edge_zx,
477  internal::constr_type_x,
478  internal::constr_type_y},
479  {internal::constr_face_x | internal::constr_face_z,
480  internal::constr_edge_zx,
481  0,
482  internal::constr_type_y},
483  {internal::constr_face_y | internal::constr_face_z,
484  internal::constr_edge_yz,
485  internal::constr_type_y,
486  internal::constr_type_x},
487  {internal::constr_face_y | internal::constr_face_z,
488  internal::constr_edge_yz,
489  0,
490  internal::constr_type_x},
491  {internal::constr_face_x | internal::constr_face_y,
492  internal::constr_edge_xy,
493  internal::constr_type_x | internal::constr_type_y,
494  internal::constr_type_z},
495  {internal::constr_face_x | internal::constr_face_y,
496  internal::constr_edge_xy,
497  internal::constr_type_y,
498  internal::constr_type_z},
499  {internal::constr_face_x | internal::constr_face_y,
500  internal::constr_edge_xy,
501  internal::constr_type_x,
502  internal::constr_type_z},
503  {internal::constr_face_x | internal::constr_face_y,
504  internal::constr_edge_xy,
505  0,
506  internal::constr_type_z}};
507 
508  for (unsigned int local_line = 0;
509  local_line < GeometryInfo<dim>::lines_per_cell;
510  ++local_line)
511  {
512  // If we don't already have a constraint for as part of a face
513  if (!(mask & line_to_edge[local_line][0]))
514  {
515  // For each cell which share that edge
516  const unsigned int line = cell->line(local_line)->index();
517  for (const auto edge_neighbor : line_to_cells[line])
518  {
519  // If one of them is coarser than us
520  const cell_iterator neighbor_cell = edge_neighbor.first;
521  if (neighbor_cell->level() < cell->level())
522  {
523  const unsigned int local_line_neighbor =
524  edge_neighbor.second;
525  mask |= line_to_edge[local_line][1] |
526  line_to_edge[local_line][2];
527 
528  bool flipped = false;
529  if (cell->line(local_line)->vertex_index(0) ==
530  neighbor_cell->line(local_line_neighbor)
531  ->vertex_index(0))
532  {
533  // Assuming line directions match axes directions,
534  // we have an unflipped edge of first type
535  mask |= line_to_edge[local_line][3];
536  }
537  else if (cell->line(local_line)->vertex_index(1) ==
538  neighbor_cell->line(local_line_neighbor)
539  ->vertex_index(1))
540  {
541  // We have an unflipped edge of second type
542  }
543  else if (cell->line(local_line)->vertex_index(1) ==
544  neighbor_cell->line(local_line_neighbor)
545  ->vertex_index(0))
546  {
547  // We have a flipped edge of second type
548  flipped = true;
549  }
550  else if (cell->line(local_line)->vertex_index(0) ==
551  neighbor_cell->line(local_line_neighbor)
552  ->vertex_index(1))
553  {
554  // We have a flipped edge of first type
555  mask |= line_to_edge[local_line][3];
556  flipped = true;
557  }
558  else
560 
561  // Copy the unconstrained values
562  neighbor_dofs.resize(n_dofs_1d * n_dofs_1d *
563  n_dofs_1d);
564  neighbor_cell->get_dof_indices(neighbor_dofs);
565  // If the vector is distributed, we need to transform
566  // the global indices to local ones.
567  if (partitioner)
568  for (auto &index : neighbor_dofs)
569  index = partitioner->global_to_local(index);
570 
571  for (unsigned int i = 0; i < n_dofs_1d; ++i)
572  {
573  // Get local dof index along line
574  const unsigned int idx =
575  line_dof_idx(local_line, i, n_dofs_1d);
576  dof_indices[idx] = neighbor_dofs
577  [lexicographic_mapping[line_dof_idx(
578  local_line_neighbor,
579  flipped ? fe_degree - i : i,
580  n_dofs_1d)]];
581  }
582 
583  // Stop looping over edge neighbors
584  break;
585  }
586  }
587  }
588  }
589  }
590  }
591 
592 
593 
594  template <int dim>
595  void
597  unsigned int &subface_index) const
598  {
599  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
600 
601  times = times % 4;
602  times = times < 0 ? times + 4 : times;
603  for (int t = 0; t < times; ++t)
604  subface_index = rot_mapping[subface_index];
605  }
606 
607 
608 
609  template <int dim>
610  void
611  HangingNodes<dim>::rotate_face(
612  int times,
613  unsigned int n_dofs_1d,
614  std::vector<types::global_dof_index> &dofs) const
615  {
616  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
617 
618  times = times % 4;
619  times = times < 0 ? times + 4 : times;
620 
621  std::vector<types::global_dof_index> copy(dofs.size());
622  for (int t = 0; t < times; ++t)
623  {
624  std::swap(copy, dofs);
625 
626  // Vertices
627  for (unsigned int i = 0; i < 4; ++i)
628  dofs[rot_mapping[i]] = copy[i];
629 
630  // Edges
631  const unsigned int n_int = n_dofs_1d - 2;
632  unsigned int offset = 4;
633  for (unsigned int i = 0; i < n_int; ++i)
634  {
635  // Left edge
636  dofs[offset + i] = copy[offset + 2 * n_int + (n_int - 1 - i)];
637  // Right edge
638  dofs[offset + n_int + i] =
639  copy[offset + 3 * n_int + (n_int - 1 - i)];
640  // Bottom edge
641  dofs[offset + 2 * n_int + i] = copy[offset + n_int + i];
642  // Top edge
643  dofs[offset + 3 * n_int + i] = copy[offset + i];
644  }
645 
646  // Interior points
647  offset += 4 * n_int;
648 
649  for (unsigned int i = 0; i < n_int; ++i)
650  for (unsigned int j = 0; j < n_int; ++j)
651  dofs[offset + i * n_int + j] =
652  copy[offset + j * n_int + (n_int - 1 - i)];
653  }
654  }
655 
656 
657 
658  template <int dim>
659  unsigned int
660  HangingNodes<dim>::line_dof_idx(int local_line,
661  unsigned int dof,
662  unsigned int n_dofs_1d) const
663  {
664  unsigned int x, y, z;
665 
666  if (local_line < 8)
667  {
668  x =
669  (local_line % 4 == 0) ? 0 : (local_line % 4 == 1) ? fe_degree : dof;
670  y =
671  (local_line % 4 == 2) ? 0 : (local_line % 4 == 3) ? fe_degree : dof;
672  z = (local_line / 4) * fe_degree;
673  }
674  else
675  {
676  x = ((local_line - 8) % 2) * fe_degree;
677  y = ((local_line - 8) / 2) * fe_degree;
678  z = dof;
679  }
680 
681  return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
682  }
683 
684 
685 
686  template <int dim>
687  void
688  HangingNodes<dim>::transpose_face(
689  std::vector<types::global_dof_index> &dofs) const
690  {
691  const std::vector<types::global_dof_index> copy(dofs);
692 
693  // Vertices
694  dofs[1] = copy[2];
695  dofs[2] = copy[1];
696 
697  // Edges
698  const unsigned int n_int = fe_degree - 1;
699  unsigned int offset = 4;
700  for (unsigned int i = 0; i < n_int; ++i)
701  {
702  // Right edge
703  dofs[offset + i] = copy[offset + 2 * n_int + i];
704  // Left edge
705  dofs[offset + n_int + i] = copy[offset + 3 * n_int + i];
706  // Bottom edge
707  dofs[offset + 2 * n_int + i] = copy[offset + i];
708  // Top edge
709  dofs[offset + 3 * n_int + i] = copy[offset + n_int + i];
710  }
711 
712  // Interior
713  offset += 4 * n_int;
714  for (unsigned int i = 0; i < n_int; ++i)
715  for (unsigned int j = 0; j < n_int; ++j)
716  dofs[offset + i * n_int + j] = copy[offset + j * n_int + i];
717  }
718 
719 
720 
721  template <int dim>
722  void
723  HangingNodes<dim>::transpose_subface_index(unsigned int &subface) const
724  {
725  if (subface == 1)
726  subface = 2;
727  else if (subface == 2)
728  subface = 1;
729  }
730 
731 
732 
733  //-------------------------------------------------------------------------//
734  // Functions for resolving the hanging node constraints //
735  //-------------------------------------------------------------------------//
736  namespace internal
737  {
738  template <unsigned int size>
739  __device__ inline unsigned int
740  index2(unsigned int i, unsigned int j)
741  {
742  return i + size * j;
743  }
744 
745 
746 
747  template <unsigned int size>
748  __device__ inline unsigned int
749  index3(unsigned int i, unsigned int j, unsigned int k)
750  {
751  return i + size * j + size * size * k;
752  }
753 
754 
755 
756  template <unsigned int fe_degree,
757  unsigned int direction,
758  bool transpose,
759  typename Number>
760  __device__ inline void
761  interpolate_boundary_2d(const unsigned int constr, Number *values)
762  {
763  const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
764  const unsigned int y_idx = threadIdx.y;
765 
766  const unsigned int this_type =
767  (direction == 0) ? internal::constr_type_x : internal::constr_type_y;
768 
769  if (constr & (((direction == 0) ? internal::constr_face_y : 0) |
770  ((direction == 1) ? internal::constr_face_x : 0)))
771  {
772  const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
773 
774  __syncthreads();
775 
776  Number t = 0;
777  const bool flag =
778  ((direction == 0) &&
779  ((constr & internal::constr_type_y) ? (y_idx == 0) :
780  (y_idx == fe_degree))) ||
781  ((direction == 1) &&
782  ((constr & internal::constr_type_x) ? (x_idx == 0) :
783  (x_idx == fe_degree)));
784 
785  if (flag)
786  {
787  const bool type = constr & this_type;
788 
789  if (type)
790  {
791  for (unsigned int i = 0; i <= fe_degree; ++i)
792  {
793  const unsigned int real_idx =
794  (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
795  index2<fe_degree + 1>(x_idx, i);
796 
797  const Number w =
798  transpose ?
799  internal::constraint_weights[i * (fe_degree + 1) +
800  interp_idx] :
801  internal::constraint_weights[interp_idx *
802  (fe_degree + 1) +
803  i];
804  t += w * values[real_idx];
805  }
806  }
807  else
808  {
809  for (unsigned int i = 0; i <= fe_degree; ++i)
810  {
811  const unsigned int real_idx =
812  (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
813  index2<fe_degree + 1>(x_idx, i);
814 
815  const Number w =
816  transpose ?
817  internal::constraint_weights[(fe_degree - i) *
818  (fe_degree + 1) +
819  fe_degree -
820  interp_idx] :
821  internal::constraint_weights[(fe_degree -
822  interp_idx) *
823  (fe_degree + 1) +
824  fe_degree - i];
825  t += w * values[real_idx];
826  }
827  }
828  }
829 
830  __syncthreads();
831 
832  if (flag)
833  values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
834  }
835  }
836 
837 
838 
839  template <unsigned int fe_degree,
840  unsigned int direction,
841  bool transpose,
842  typename Number>
843  __device__ inline void
844  interpolate_boundary_3d(const unsigned int constr, 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 
878  if (constr & (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  __syncthreads();
888 
889  Number t = 0;
890  const bool on_face1 = (constr & face1_type) ?
891  (face1_idx == 0) :
892  (face1_idx == fe_degree);
893  const bool on_face2 = (constr & face2_type) ?
894  (face2_idx == 0) :
895  (face2_idx == fe_degree);
896  const bool flag = (((constr & face1) && on_face1) ||
897  ((constr & face2) && on_face2) ||
898  ((constr & edge) && on_face1 && on_face2));
899 
900  if (flag)
901  {
902  const bool type = constr & this_type;
903  if (type)
904  {
905  for (unsigned int i = 0; i <= fe_degree; ++i)
906  {
907  const unsigned int real_idx =
908  (direction == 0) ?
909  index3<fe_degree + 1>(i, y_idx, z_idx) :
910  (direction == 1) ?
911  index3<fe_degree + 1>(x_idx, i, z_idx) :
912  index3<fe_degree + 1>(x_idx, y_idx, i);
913 
914  const Number w =
915  transpose ?
916  internal::constraint_weights[i * (fe_degree + 1) +
917  interp_idx] :
918  internal::constraint_weights[interp_idx *
919  (fe_degree + 1) +
920  i];
921  t += w * values[real_idx];
922  }
923  }
924  else
925  {
926  for (unsigned int i = 0; i <= fe_degree; ++i)
927  {
928  const unsigned int real_idx =
929  (direction == 0) ?
930  index3<fe_degree + 1>(i, y_idx, z_idx) :
931  (direction == 1) ?
932  index3<fe_degree + 1>(x_idx, i, z_idx) :
933  index3<fe_degree + 1>(x_idx, y_idx, i);
934 
935  const Number w =
936  transpose ?
937  internal::constraint_weights[(fe_degree - i) *
938  (fe_degree + 1) +
939  fe_degree -
940  interp_idx] :
941  internal::constraint_weights[(fe_degree -
942  interp_idx) *
943  (fe_degree + 1) +
944  fe_degree - i];
945  t += w * values[real_idx];
946  }
947  }
948  }
949 
950  __syncthreads();
951 
952  if (flag)
953  values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
954  }
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 constr, Number *values)
971  {
972  if (dim == 2)
973  {
974  internal::interpolate_boundary_2d<fe_degree, 0, transpose>(constr,
975  values);
976  internal::interpolate_boundary_2d<fe_degree, 1, transpose>(constr,
977  values);
978  }
979  else if (dim == 3)
980  {
981  // Interpolate y and z faces (x-direction)
982  internal::interpolate_boundary_3d<fe_degree, 0, transpose>(constr,
983  values);
984  // Interpolate x and z faces (y-direction)
985  internal::interpolate_boundary_3d<fe_degree, 1, transpose>(constr,
986  values);
987  // Interpolate x and y faces (z-direction)
988  internal::interpolate_boundary_3d<fe_degree, 2, transpose>(constr,
989  values);
990  }
991  }
992  } // namespace internal
993 } // namespace CUDAWrappers
994 
995 DEAL_II_NAMESPACE_CLOSE
996 
997 #endif
998 
999 #endif
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
T fixed_power(const T t)
Definition: utilities.h:1031
static ::ExceptionBase & ExcMessage(std::string arg1)
Definition: fe_q.h:554
HangingNodes(unsigned int fe_degree, const DoFHandler< dim > &dof_handler, const std::vector< unsigned int > &lexicographic_mapping)
#define AssertCuda(error_code)
Definition: exceptions.h:1722
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
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
void lexicographic_to_hierarchic_numbering(const FiniteElementData< dim > &fe_data, std::vector< unsigned int > &l2h)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInternalError()