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