Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18:55:02+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\}}\)
hanging_nodes_internal.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2022 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_hanging_nodes_internal_h
17 #define dealii_hanging_nodes_internal_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/ndarray.h>
22 #include <deal.II/base/utilities.h>
23 
25 
26 #include <deal.II/fe/fe_q.h>
27 #include <deal.II/fe/fe_tools.h>
28 
30 
32 
33 namespace internal
34 {
35  namespace MatrixFreeFunctions
36  {
50  enum class ConstraintKinds : std::uint16_t
51  {
52  // default: unconstrained cell
53  unconstrained = 0,
54 
55  // subcell
56  subcell_x = 1 << 0,
57  subcell_y = 1 << 1,
58  subcell_z = 1 << 2,
59 
60  // face is constrained
61  face_x = 1 << 3,
62  face_y = 1 << 4,
63  face_z = 1 << 5,
64 
65  // edge is constrained
66  edge_x = 1 << 6,
67  edge_y = 1 << 7,
68  edge_z = 1 << 8
69  };
70 
71 
72 
76  using compressed_constraint_kind = std::uint8_t;
77 
78 
79 
85 
86 
87 
91  inline bool
92  check(const ConstraintKinds kind_in, const unsigned int dim)
93  {
94  const std::uint16_t kind = static_cast<std::uint16_t>(kind_in);
95  const std::uint16_t subcell = (kind >> 0) & 7;
96  const std::uint16_t face = (kind >> 3) & 7;
97  const std::uint16_t edge = (kind >> 6) & 7;
98 
99  if ((kind >> 9) > 0)
100  return false;
101 
102  if (dim == 2)
103  {
104  if (edge > 0)
105  return false; // in 2D there are no edge constraints
106 
107  if (subcell == 0 && face == 0)
108  return true; // no constraints
109  else if (0 < face)
110  return true; // at least one face is constrained
111  }
112  else if (dim == 3)
113  {
114  if (subcell == 0 && face == 0 && edge == 0)
115  return true; // no constraints
116  else if (0 < face && edge == 0)
117  return true; // at least one face is constrained
118  else if (0 == face && 0 < edge)
119  return true; // at least one edge is constrained
120  else if ((face == edge) && (face == 1 || face == 2 || face == 4))
121  return true; // one face and its orthogonal edge is constrained
122  }
123 
124  return false;
125  }
126 
127 
128 
134  compress(const ConstraintKinds kind_in, const unsigned int dim)
135  {
136  Assert(check(kind_in, dim), ExcInternalError());
137 
138  if (dim == 2)
139  return static_cast<compressed_constraint_kind>(kind_in);
140 
141  if (kind_in == ConstraintKinds::unconstrained)
143 
144  const std::uint16_t kind = static_cast<std::uint16_t>(kind_in);
145  const std::uint16_t subcell = (kind >> 0) & 7;
146  const std::uint16_t face = (kind >> 3) & 7;
147  const std::uint16_t edge = (kind >> 6) & 7;
148 
149  return subcell + ((face > 0) << 3) + ((edge > 0) << 4) +
150  (std::max(face, edge) << 5);
151  }
152 
153 
154 
159  inline ConstraintKinds
160  decompress(const compressed_constraint_kind kind_in, const unsigned int dim)
161  {
162  if (dim == 2)
163  return static_cast<ConstraintKinds>(kind_in);
164 
167 
168  const std::uint16_t subcell = (kind_in >> 0) & 7;
169  const std::uint16_t flag_0 = (kind_in >> 3) & 3;
170  const std::uint16_t flag_1 = (kind_in >> 5) & 7;
171 
172  const auto result = static_cast<ConstraintKinds>(
173  subcell + (((flag_0 & 0b01) ? flag_1 : 0) << 3) +
174  (((flag_0 & 0b10) ? flag_1 : 0) << 6));
175 
176  Assert(check(result, dim), ExcInternalError());
177 
178  return result;
179  }
180 
181 
182 
186  inline std::size_t
188  {
189  return sizeof(ConstraintKinds);
190  }
191 
192 
193 
203  {
204  return static_cast<ConstraintKinds>(static_cast<std::uint16_t>(f1) |
205  static_cast<std::uint16_t>(f2));
206  }
207 
208 
209 
216  {
217  f1 = f1 | f2;
218  return f1;
219  }
220 
221 
222 
226  DEAL_II_CUDA_HOST_DEV inline bool
228  {
229  return static_cast<std::uint16_t>(f1) != static_cast<std::uint16_t>(f2);
230  }
231 
232 
233 
239  operator<(const ConstraintKinds f1, const ConstraintKinds f2)
240  {
241  return static_cast<std::uint16_t>(f1) < static_cast<std::uint16_t>(f2);
242  }
243 
244 
245 
251  {
252  return static_cast<ConstraintKinds>(static_cast<std::uint16_t>(f1) &
253  static_cast<std::uint16_t>(f2));
254  }
255 
256 
257 
264  template <int dim>
266  {
267  public:
271  HangingNodes(const Triangulation<dim> &triangualtion);
272 
276  template <typename CellIterator>
277  bool
279  const CellIterator & cell,
280  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
281  const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
282  std::vector<types::global_dof_index> & dof_indices,
283  const ArrayView<ConstraintKinds> & mask) const;
284 
289  static std::vector<std::vector<bool>>
290  compute_supported_components(const ::hp::FECollection<dim> &fe);
291 
295  template <typename CellIterator>
297  compute_refinement_configuration(const CellIterator &cell) const;
298 
303  template <typename CellIterator>
304  void
306  const CellIterator & cell,
307  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
308  const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
309  const std::vector<std::vector<bool>> & component_mask,
310  const ConstraintKinds & refinement_configuration,
311  std::vector<types::global_dof_index> & dof_indices) const;
312 
313  private:
317  void
319 
320  void
321  rotate_subface_index(int times, unsigned int &subface_index) const;
322 
323  void
324  rotate_face(int times,
325  unsigned int n_dofs_1d,
326  std::vector<types::global_dof_index> &dofs) const;
327 
328  unsigned int
329  line_dof_idx(int local_line,
330  unsigned int dof,
331  unsigned int n_dofs_1d) const;
332 
333  void
334  transpose_face(const unsigned int fe_degree,
335  std::vector<types::global_dof_index> &dofs) const;
336 
337  void
338  transpose_subface_index(unsigned int &subface) const;
339 
340  std::vector<std::vector<
341  std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>>
343 
344  const ::ndarray<unsigned int, 3, 2, 2> local_lines = {
345  {{{{{7, 3}}, {{6, 2}}}},
346  {{{{5, 1}}, {{4, 0}}}},
347  {{{{11, 9}}, {{10, 8}}}}}};
348  };
349 
350 
351 
352  template <int dim>
355  {
356  // Set up line-to-cell mapping for edge constraints (only if dim = 3 and
357  // for pure hex meshes)
358  if (triangulation.all_reference_cells_are_hyper_cube())
359  setup_line_to_cell(triangulation);
360  }
361 
362 
363 
364  template <int dim>
365  inline void
368  {
369  (void)triangulation;
370  }
371 
372 
373 
374  template <>
375  inline void
377  {
378  // Check if we there are no hanging nodes on the current MPI process,
379  // which we do by checking if the second finest level holds no active
380  // non-artificial cell
381  if (triangulation.n_levels() <= 1 ||
382  std::none_of(triangulation.begin_active(triangulation.n_levels() - 2),
383  triangulation.end_active(triangulation.n_levels() - 2),
384  [](const CellAccessor<3, 3> &cell) {
385  return !cell.is_artificial();
386  }))
387  return;
388 
389  const unsigned int n_raw_lines = triangulation.n_raw_lines();
390  this->line_to_cells.resize(n_raw_lines);
391 
392  // In 3D, we can have DoFs on only an edge being constrained (e.g. in a
393  // cartesian 2x2x2 grid, where only the upper left 2 cells are refined).
394  // This sets up a helper data structure in the form of a mapping from
395  // edges (i.e. lines) to neighboring cells.
396 
397  // Mapping from an edge to which children that share that edge.
398  const unsigned int line_to_children[12][2] = {{0, 2},
399  {1, 3},
400  {0, 1},
401  {2, 3},
402  {4, 6},
403  {5, 7},
404  {4, 5},
405  {6, 7},
406  {0, 4},
407  {1, 5},
408  {2, 6},
409  {3, 7}};
410 
411  std::vector<std::vector<
412  std::pair<typename Triangulation<3>::cell_iterator, unsigned int>>>
413  line_to_inactive_cells(n_raw_lines);
414 
415  // First add active and inactive cells to their lines:
416  for (const auto &cell : triangulation.cell_iterators())
417  {
418  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
419  ++line)
420  {
421  const unsigned int line_idx = cell->line(line)->index();
422  if (cell->is_active())
423  line_to_cells[line_idx].push_back(std::make_pair(cell, line));
424  else
425  line_to_inactive_cells[line_idx].push_back(
426  std::make_pair(cell, line));
427  }
428  }
429 
430  // Now, we can access edge-neighboring active cells on same level to also
431  // access of an edge to the edges "children". These are found from looking
432  // at the corresponding edge of children of inactive edge neighbors.
433  for (unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
434  {
435  if ((line_to_cells[line_idx].size() > 0) &&
436  line_to_inactive_cells[line_idx].size() > 0)
437  {
438  // We now have cells to add (active ones) and edges to which they
439  // should be added (inactive cells).
440  const auto &inactive_cell =
441  line_to_inactive_cells[line_idx][0].first;
442  const unsigned int neighbor_line =
443  line_to_inactive_cells[line_idx][0].second;
444 
445  for (unsigned int c = 0; c < 2; ++c)
446  {
447  const auto &child =
448  inactive_cell->child(line_to_children[neighbor_line][c]);
449  const unsigned int child_line_idx =
450  child->line(neighbor_line)->index();
451 
452  // Now add all active cells
453  for (const auto &cl : line_to_cells[line_idx])
454  line_to_cells[child_line_idx].push_back(cl);
455  }
456  }
457  }
458  }
459 
460 
461 
462  template <int dim>
463  inline std::vector<std::vector<bool>>
465  const ::hp::FECollection<dim> &fe_collection)
466  {
467  std::vector<std::vector<bool>> supported_components(
468  fe_collection.size(),
469  std::vector<bool>(fe_collection.n_components(), false));
470 
471  for (unsigned int i = 0; i < fe_collection.size(); ++i)
472  {
473  for (unsigned int base_element_index = 0, comp = 0;
474  base_element_index < fe_collection[i].n_base_elements();
475  ++base_element_index)
476  for (unsigned int c = 0;
477  c < fe_collection[i].element_multiplicity(base_element_index);
478  ++c, ++comp)
479  if (dim == 1 || dynamic_cast<const FE_Q<dim> *>(
480  &fe_collection[i].base_element(
481  base_element_index)) == nullptr)
482  supported_components[i][comp] = false;
483  else
484  supported_components[i][comp] = true;
485  }
486 
487  return supported_components;
488  }
489 
490 
491 
492  template <int dim>
493  template <typename CellIterator>
494  inline ConstraintKinds
496  const CellIterator &cell) const
497  {
498  // TODO: for simplex or mixed meshes: nothing to do
499  if ((dim == 3 && line_to_cells.size() == 0) ||
500  (cell->reference_cell().is_hyper_cube() == false))
502 
503  if (cell->level() == 0)
505 
506  const std::uint16_t subcell =
507  cell->parent()->child_iterator_to_index(cell);
508  const std::uint16_t subcell_x = (subcell >> 0) & 1;
509  const std::uint16_t subcell_y = (subcell >> 1) & 1;
510  const std::uint16_t subcell_z = (subcell >> 2) & 1;
511 
512  std::uint16_t face = 0;
513  std::uint16_t edge = 0;
514 
515  for (unsigned int direction = 0; direction < dim; ++direction)
516  {
517  const auto side = (subcell >> direction) & 1U;
518  const auto face_no = direction * 2 + side;
519 
520  // ignore if at boundary
521  if (cell->at_boundary(face_no))
522  continue;
523 
524  const auto &neighbor = cell->neighbor(face_no);
525 
526  // ignore neighbors that are artificial or have the same level or
527  // have children
528  if (neighbor->has_children() || neighbor->is_artificial() ||
529  neighbor->level() == cell->level())
530  continue;
531 
532  face |= 1 << direction;
533  }
534 
535  if (dim == 3)
536  for (unsigned int direction = 0; direction < dim; ++direction)
537  if (face == 0 || face == (1 << direction))
538  {
539  const unsigned int line_no =
540  direction == 0 ?
541  (local_lines[0][subcell_y == 0][subcell_z == 0]) :
542  (direction == 1 ?
543  (local_lines[1][subcell_x == 0][subcell_z == 0]) :
544  (local_lines[2][subcell_x == 0][subcell_y == 0]));
545 
546  const unsigned int line_index = cell->line(line_no)->index();
547 
548  const auto edge_neighbor =
549  std::find_if(line_to_cells[line_index].begin(),
550  line_to_cells[line_index].end(),
551  [&cell](const auto &edge_neighbor) {
552  return edge_neighbor.first->is_artificial() ==
553  false &&
554  edge_neighbor.first->level() <
555  cell->level();
556  });
557 
558  if (edge_neighbor == line_to_cells[line_index].end())
559  continue;
560 
561  edge |= 1 << direction;
562  }
563 
564  if ((face == 0) && (edge == 0))
566 
567  const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7));
568 
569  const auto refinement_configuration = static_cast<ConstraintKinds>(
570  inverted_subcell + (face << 3) + (edge << 6));
571  Assert(check(refinement_configuration, dim), ExcInternalError());
572  return refinement_configuration;
573  }
574 
575 
576 
577  template <int dim>
578  template <typename CellIterator>
579  inline void
581  const CellIterator & cell,
582  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
583  const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
584  const std::vector<std::vector<bool>> & supported_components,
585  const ConstraintKinds & refinement_configuration,
586  std::vector<types::global_dof_index> & dof_indices) const
587  {
588  if (std::find(supported_components[cell->active_fe_index()].begin(),
589  supported_components[cell->active_fe_index()].end(),
590  true) ==
591  supported_components[cell->active_fe_index()].end())
592  return;
593 
594  const auto &fe = cell->get_fe();
595 
596  AssertDimension(fe.n_unique_faces(), 1);
597 
598  std::vector<std::vector<unsigned int>>
599  component_to_system_index_face_array(fe.n_components());
600 
601  for (unsigned int i = 0; i < fe.n_dofs_per_face(0); ++i)
602  component_to_system_index_face_array
603  [fe.face_system_to_component_index(i, /*face_no=*/0).first]
604  .push_back(i);
605 
606  std::vector<unsigned int> idx_offset = {0};
607 
608  for (unsigned int base_element_index = 0;
609  base_element_index < cell->get_fe().n_base_elements();
610  ++base_element_index)
611  for (unsigned int c = 0;
612  c < cell->get_fe().element_multiplicity(base_element_index);
613  ++c)
614  idx_offset.push_back(
615  idx_offset.back() +
616  cell->get_fe().base_element(base_element_index).n_dofs_per_cell());
617 
618  std::vector<types::global_dof_index> neighbor_dofs_all(idx_offset.back());
619  std::vector<types::global_dof_index> neighbor_dofs_all_temp(
620  idx_offset.back());
621  std::vector<types::global_dof_index> neighbor_dofs_face(
622  fe.n_dofs_per_face(/*face_no=*/0));
623 
624 
625  const auto get_face_idx = [](const auto n_dofs_1d,
626  const auto face_no,
627  const auto i,
628  const auto j) -> unsigned int {
629  const auto direction = face_no / 2;
630  const auto side = face_no % 2;
631  const auto offset = (side == 1) ? (n_dofs_1d - 1) : 0;
632 
633  if (dim == 2)
634  return (direction == 0) ? (n_dofs_1d * i + offset) :
635  (n_dofs_1d * offset + i);
636  else if (dim == 3)
637  switch (direction)
638  {
639  case 0:
640  return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset;
641  case 1:
642  return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i;
643  case 2:
644  return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j;
645  default:
646  Assert(false, ExcNotImplemented());
647  }
648 
649  Assert(false, ExcNotImplemented());
650 
651  return 0;
652  };
653 
654  const std::uint16_t kind =
655  static_cast<std::uint16_t>(refinement_configuration);
656  const std::uint16_t subcell = (kind >> 0) & 7;
657  const std::uint16_t subcell_x = (subcell >> 0) & 1;
658  const std::uint16_t subcell_y = (subcell >> 1) & 1;
659  const std::uint16_t subcell_z = (subcell >> 2) & 1;
660  const std::uint16_t face = (kind >> 3) & 7;
661  const std::uint16_t edge = (kind >> 6) & 7;
662 
663  for (unsigned int direction = 0; direction < dim; ++direction)
664  if ((face >> direction) & 1U)
665  {
666  const auto side = ((subcell >> direction) & 1U) == 0;
667  const auto face_no = direction * 2 + side;
668 
669  // read DoFs of parent of face, ...
670  cell->neighbor(face_no)
671  ->face(cell->neighbor_face_no(face_no))
672  ->get_dof_indices(neighbor_dofs_face,
673  cell->neighbor(face_no)->active_fe_index());
674 
675  // ... convert the global DoFs to serial ones, and ...
676  if (partitioner)
677  for (auto &index : neighbor_dofs_face)
678  index = partitioner->global_to_local(index);
679 
680  for (unsigned int base_element_index = 0, comp = 0;
681  base_element_index < cell->get_fe().n_base_elements();
682  ++base_element_index)
683  for (unsigned int c = 0;
684  c < cell->get_fe().element_multiplicity(base_element_index);
685  ++c, ++comp)
686  {
687  if (supported_components[cell->active_fe_index()][comp] ==
688  false)
689  continue;
690 
691  const unsigned int n_dofs_1d =
692  cell->get_fe()
693  .base_element(base_element_index)
694  .tensor_degree() +
695  1;
696  const unsigned int dofs_per_face =
697  Utilities::pow(n_dofs_1d, dim - 1);
698  std::vector<types::global_dof_index> neighbor_dofs(
699  dofs_per_face);
700  const auto lex_face_mapping =
702  n_dofs_1d - 1);
703 
704  // ... extract the DoFs of the current component
705  for (unsigned int i = 0; i < dofs_per_face; ++i)
706  neighbor_dofs[i] = neighbor_dofs_face
707  [component_to_system_index_face_array[comp][i]];
708 
709  // fix DoFs depending on orientation, flip, and rotation
710  if (dim == 2)
711  {
712  // TODO: for mixed meshes we need to take care of
713  // orientation here
714  Assert(cell->face_orientation(face_no),
716  }
717  else if (dim == 3)
718  {
719  int rotate = 0; // TODO
720  if (cell->face_rotation(face_no)) //
721  rotate -= 1; //
722  if (cell->face_flip(face_no)) //
723  rotate -= 2; //
724 
725  rotate_face(rotate, n_dofs_1d, neighbor_dofs);
726 
727  if (cell->face_orientation(face_no) == false)
728  transpose_face(n_dofs_1d - 1, neighbor_dofs);
729  }
730  else
731  {
732  Assert(false, ExcNotImplemented());
733  }
734 
735  // update DoF map
736  for (unsigned int i = 0, k = 0; i < n_dofs_1d; ++i)
737  for (unsigned int j = 0; j < (dim == 2 ? 1 : n_dofs_1d);
738  ++j, ++k)
739  dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) +
740  idx_offset[comp]] =
741  neighbor_dofs[lex_face_mapping[k]];
742  }
743  }
744 
745  if (dim == 3)
746  for (unsigned int direction = 0; direction < dim; ++direction)
747  if ((edge >> direction) & 1U)
748  {
749  const unsigned int line_no =
750  direction == 0 ?
751  (local_lines[0][subcell_y][subcell_z]) :
752  (direction == 1 ? (local_lines[1][subcell_x][subcell_z]) :
753  (local_lines[2][subcell_x][subcell_y]));
754 
755  const unsigned int line_index = cell->line(line_no)->index();
756 
757  const auto edge_neighbor =
758  std::find_if(line_to_cells[line_index].begin(),
759  line_to_cells[line_index].end(),
760  [&cell](const auto &edge_neighbor) {
761  return edge_neighbor.first->is_artificial() ==
762  false &&
763  edge_neighbor.first->level() <
764  cell->level();
765  });
766 
767  if (edge_neighbor == line_to_cells[line_index].end())
768  continue;
769 
770  const auto neighbor_cell = edge_neighbor->first;
771  const auto local_line_neighbor = edge_neighbor->second;
772 
774  &neighbor_cell->get_triangulation(),
775  neighbor_cell->level(),
776  neighbor_cell->index(),
777  &cell->get_dof_handler())
778  .get_dof_indices(neighbor_dofs_all);
779 
780  if (partitioner)
781  for (auto &index : neighbor_dofs_all)
782  index = partitioner->global_to_local(index);
783 
784  for (unsigned int i = 0; i < neighbor_dofs_all_temp.size(); ++i)
785  neighbor_dofs_all_temp[i] = neighbor_dofs_all
786  [lexicographic_mapping[cell->active_fe_index()][i]];
787 
788  const bool flipped =
789  cell->line_orientation(line_no) !=
790  neighbor_cell->line_orientation(local_line_neighbor);
791 
792  for (unsigned int base_element_index = 0, comp = 0;
793  base_element_index < cell->get_fe().n_base_elements();
794  ++base_element_index)
795  for (unsigned int c = 0;
796  c <
797  cell->get_fe().element_multiplicity(base_element_index);
798  ++c, ++comp)
799  {
800  if (supported_components[cell->active_fe_index()][comp] ==
801  false)
802  continue;
803 
804  const unsigned int n_dofs_1d =
805  cell->get_fe()
806  .base_element(base_element_index)
807  .tensor_degree() +
808  1;
809 
810  for (unsigned int i = 0; i < n_dofs_1d; ++i)
811  dof_indices[line_dof_idx(line_no, i, n_dofs_1d) +
812  idx_offset[comp]] = neighbor_dofs_all_temp
813  [line_dof_idx(local_line_neighbor,
814  flipped ? (n_dofs_1d - 1 - i) : i,
815  n_dofs_1d) +
816  idx_offset[comp]];
817  }
818  }
819  }
820 
821 
822 
823  template <int dim>
824  template <typename CellIterator>
825  inline bool
827  const CellIterator & cell,
828  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
829  const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
830  std::vector<types::global_dof_index> & dof_indices,
831  const ArrayView<ConstraintKinds> & masks) const
832  {
833  // 1) check if finite elements support fast hanging-node algorithm
834  const auto supported_components = compute_supported_components(
835  cell->get_dof_handler().get_fe_collection());
836 
837  if ([](const auto &supported_components) {
838  return std::none_of(supported_components.begin(),
839  supported_components.end(),
840  [](const auto &a) {
841  return *std::max_element(a.begin(), a.end());
842  });
843  }(supported_components))
844  return false;
845 
846  // 2) determine the refinement configuration of the cell
847  const auto refinement_configuration =
848  compute_refinement_configuration(cell);
849 
850  if (refinement_configuration == ConstraintKinds::unconstrained)
851  return false;
852 
853  // 3) update DoF indices of cell for specified components
854  update_dof_indices(cell,
855  partitioner,
856  lexicographic_mapping,
857  supported_components,
858  refinement_configuration,
859  dof_indices);
860 
861  // 4) TODO: copy refinement configuration to all components
862  for (unsigned int c = 0; c < supported_components[0].size(); ++c)
863  if (supported_components[cell->active_fe_index()][c])
864  masks[c] = refinement_configuration;
865 
866  return true;
867  }
868 
869 
870 
871  template <int dim>
872  inline void
874  unsigned int &subface_index) const
875  {
876  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
877 
878  times = times % 4;
879  times = times < 0 ? times + 4 : times;
880  for (int t = 0; t < times; ++t)
881  subface_index = rot_mapping[subface_index];
882  }
883 
884 
885 
886  template <int dim>
887  inline void
889  int times,
890  unsigned int n_dofs_1d,
891  std::vector<types::global_dof_index> &dofs) const
892  {
893  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
894 
895  times = times % 4;
896  times = times < 0 ? times + 4 : times;
897 
898  std::vector<types::global_dof_index> copy(dofs.size());
899  for (int t = 0; t < times; ++t)
900  {
901  std::swap(copy, dofs);
902 
903  // Vertices
904  for (unsigned int i = 0; i < 4; ++i)
905  dofs[rot_mapping[i]] = copy[i];
906 
907  // Edges
908  const unsigned int n_int = n_dofs_1d - 2;
909  unsigned int offset = 4;
910  for (unsigned int i = 0; i < n_int; ++i)
911  {
912  // Left edge
913  dofs[offset + i] = copy[offset + 2 * n_int + (n_int - 1 - i)];
914  // Right edge
915  dofs[offset + n_int + i] =
916  copy[offset + 3 * n_int + (n_int - 1 - i)];
917  // Bottom edge
918  dofs[offset + 2 * n_int + i] = copy[offset + n_int + i];
919  // Top edge
920  dofs[offset + 3 * n_int + i] = copy[offset + i];
921  }
922 
923  // Interior points
924  offset += 4 * n_int;
925 
926  for (unsigned int i = 0; i < n_int; ++i)
927  for (unsigned int j = 0; j < n_int; ++j)
928  dofs[offset + i * n_int + j] =
929  copy[offset + j * n_int + (n_int - 1 - i)];
930  }
931  }
932 
933 
934 
935  template <int dim>
936  inline unsigned int
938  unsigned int dof,
939  unsigned int n_dofs_1d) const
940  {
941  unsigned int x, y, z;
942 
943  const unsigned int fe_degree = n_dofs_1d - 1;
944 
945  if (local_line < 8)
946  {
947  x = (local_line % 4 == 0) ? 0 :
948  (local_line % 4 == 1) ? fe_degree :
949  dof;
950  y = (local_line % 4 == 2) ? 0 :
951  (local_line % 4 == 3) ? fe_degree :
952  dof;
953  z = (local_line / 4) * fe_degree;
954  }
955  else
956  {
957  x = ((local_line - 8) % 2) * fe_degree;
958  y = ((local_line - 8) / 2) * fe_degree;
959  z = dof;
960  }
961 
962  return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
963  }
964 
965 
966 
967  template <int dim>
968  inline void
970  const unsigned int fe_degree,
971  std::vector<types::global_dof_index> &dofs) const
972  {
973  const std::vector<types::global_dof_index> copy(dofs);
974 
975  // Vertices
976  dofs[1] = copy[2];
977  dofs[2] = copy[1];
978 
979  // Edges
980  const unsigned int n_int = fe_degree - 1;
981  unsigned int offset = 4;
982  for (unsigned int i = 0; i < n_int; ++i)
983  {
984  // Right edge
985  dofs[offset + i] = copy[offset + 2 * n_int + i];
986  // Left edge
987  dofs[offset + n_int + i] = copy[offset + 3 * n_int + i];
988  // Bottom edge
989  dofs[offset + 2 * n_int + i] = copy[offset + i];
990  // Top edge
991  dofs[offset + 3 * n_int + i] = copy[offset + n_int + i];
992  }
993 
994  // Interior
995  offset += 4 * n_int;
996  for (unsigned int i = 0; i < n_int; ++i)
997  for (unsigned int j = 0; j < n_int; ++j)
998  dofs[offset + i * n_int + j] = copy[offset + j * n_int + i];
999  }
1000 
1001 
1002 
1003  template <int dim>
1004  void
1005  HangingNodes<dim>::transpose_subface_index(unsigned int &subface) const
1006  {
1007  if (subface == 1)
1008  subface = 2;
1009  else if (subface == 2)
1010  subface = 1;
1011  }
1012  } // namespace MatrixFreeFunctions
1013 } // namespace internal
1014 
1015 
1017 
1018 #endif
bool is_active() const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
Definition: fe_q.h:549
typename ::internal::TriangulationImplementation::Iterators< dim, spacedim >::line_iterator line(const unsigned int i) const
ConstraintKinds compute_refinement_configuration(const CellIterator &cell) const
void setup_line_to_cell(const Triangulation< dim > &triangulation)
void transpose_face(const unsigned int fe_degree, std::vector< types::global_dof_index > &dofs) const
static std::vector< std::vector< bool > > compute_supported_components(const ::hp::FECollection< dim > &fe)
std::vector< std::vector< std::pair< typename Triangulation< dim >::cell_iterator, unsigned int > > > line_to_cells
HangingNodes(const Triangulation< dim > &triangualtion)
void transpose_subface_index(unsigned int &subface) const
void rotate_face(int times, unsigned int n_dofs_1d, std::vector< types::global_dof_index > &dofs) const
unsigned int line_dof_idx(int local_line, unsigned int dof, unsigned int n_dofs_1d) const
bool setup_constraints(const CellIterator &cell, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner, const std::vector< std::vector< unsigned int >> &lexicographic_mapping, std::vector< types::global_dof_index > &dof_indices, const ArrayView< ConstraintKinds > &mask) const
const ::ndarray< unsigned int, 3, 2, 2 > local_lines
void update_dof_indices(const CellIterator &cell, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner, const std::vector< std::vector< unsigned int >> &lexicographic_mapping, const std::vector< std::vector< bool >> &component_mask, const ConstraintKinds &refinement_configuration, std::vector< types::global_dof_index > &dof_indices) const
void rotate_subface_index(int times, unsigned int &subface_index) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
void rotate(const double angle, Triangulation< dim > &triangulation)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
bool operator<(const ConstraintKinds f1, const ConstraintKinds f2)
ConstraintKinds decompress(const compressed_constraint_kind kind_in, const unsigned int dim)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
compressed_constraint_kind compress(const ConstraintKinds kind_in, const unsigned int dim)
std::size_t memory_consumption(const ConstraintKinds &)
std::uint8_t compressed_constraint_kind
Definition: dof_info.h:83
ConstraintKinds operator&(const ConstraintKinds f1, const ConstraintKinds f2)
bool operator!=(const ConstraintKinds f1, const ConstraintKinds f2)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
ConstraintKinds & operator|=(ConstraintKinds &f1, const ConstraintKinds f2)
ConstraintKinds operator|(const ConstraintKinds f1, const ConstraintKinds f2)
void copy(const T *begin, const T *end, U *dest)
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation