Reference documentation for deal.II version GIT 85919f3e70 2023-05-28 07:10:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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_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<std::vector<
348  std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>>
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  return MemoryConsumption::memory_consumption(line_to_cells);
376  }
377 
378 
379 
380  template <int dim>
381  inline void
384  {
385  (void)triangulation;
386  }
387 
388 
389 
390  template <>
391  inline void
393  {
394  // Check if we there are no hanging nodes on the current MPI process,
395  // which we do by checking if the second finest level holds no active
396  // non-artificial cell
397  if (triangulation.n_levels() <= 1 ||
398  std::none_of(triangulation.begin_active(triangulation.n_levels() - 2),
399  triangulation.end_active(triangulation.n_levels() - 2),
400  [](const CellAccessor<3, 3> &cell) {
401  return !cell.is_artificial();
402  }))
403  return;
404 
405  const unsigned int n_raw_lines = triangulation.n_raw_lines();
406  this->line_to_cells.resize(n_raw_lines);
407 
408  // In 3d, we can have DoFs on only an edge being constrained (e.g. in a
409  // cartesian 2x2x2 grid, where only the upper left 2 cells are refined).
410  // This sets up a helper data structure in the form of a mapping from
411  // edges (i.e. lines) to neighboring cells.
412 
413  // Mapping from an edge to which children that share that edge.
414  const unsigned int line_to_children[12][2] = {{0, 2},
415  {1, 3},
416  {0, 1},
417  {2, 3},
418  {4, 6},
419  {5, 7},
420  {4, 5},
421  {6, 7},
422  {0, 4},
423  {1, 5},
424  {2, 6},
425  {3, 7}};
426 
427  std::vector<std::vector<
428  std::pair<typename Triangulation<3>::cell_iterator, unsigned int>>>
429  line_to_inactive_cells(n_raw_lines);
430 
431  // First add active and inactive cells to their lines:
432  for (const auto &cell : triangulation.cell_iterators())
433  {
434  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
435  ++line)
436  {
437  const unsigned int line_idx = cell->line(line)->index();
438  if (cell->is_active())
439  line_to_cells[line_idx].push_back(std::make_pair(cell, line));
440  else
441  line_to_inactive_cells[line_idx].push_back(
442  std::make_pair(cell, line));
443  }
444  }
445 
446  // Now, we can access edge-neighboring active cells on same level to also
447  // access of an edge to the edges "children". These are found from looking
448  // at the corresponding edge of children of inactive edge neighbors.
449  for (unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
450  {
451  if ((line_to_cells[line_idx].size() > 0) &&
452  line_to_inactive_cells[line_idx].size() > 0)
453  {
454  // We now have cells to add (active ones) and edges to which they
455  // should be added (inactive cells).
456  const auto &inactive_cell =
457  line_to_inactive_cells[line_idx][0].first;
458  const unsigned int neighbor_line =
459  line_to_inactive_cells[line_idx][0].second;
460 
461  for (unsigned int c = 0; c < 2; ++c)
462  {
463  const auto &child =
464  inactive_cell->child(line_to_children[neighbor_line][c]);
465  const unsigned int child_line_idx =
466  child->line(neighbor_line)->index();
467 
468  // Now add all active cells
469  for (const auto &cl : line_to_cells[line_idx])
470  line_to_cells[child_line_idx].push_back(cl);
471  }
472  }
473  }
474  }
475 
476 
477 
478  template <int dim>
479  inline std::vector<std::vector<bool>>
481  const ::hp::FECollection<dim> &fe_collection)
482  {
483  std::vector<std::vector<bool>> supported_components(
484  fe_collection.size(),
485  std::vector<bool>(fe_collection.n_components(), false));
486 
487  for (unsigned int i = 0; i < fe_collection.size(); ++i)
488  {
489  for (unsigned int base_element_index = 0, comp = 0;
490  base_element_index < fe_collection[i].n_base_elements();
491  ++base_element_index)
492  for (unsigned int c = 0;
493  c < fe_collection[i].element_multiplicity(base_element_index);
494  ++c, ++comp)
495  if (dim == 1 ||
496  (dynamic_cast<const FE_Q<dim> *>(
497  &fe_collection[i].base_element(base_element_index)) ==
498  nullptr &&
499  dynamic_cast<const FE_Q_iso_Q1<dim> *>(
500  &fe_collection[i].base_element(base_element_index)) ==
501  nullptr))
502  supported_components[i][comp] = false;
503  else
504  supported_components[i][comp] = true;
505  }
506 
507  return supported_components;
508  }
509 
510 
511 
512  template <int dim>
513  template <typename CellIterator>
514  inline ConstraintKinds
516  const CellIterator &cell) const
517  {
518  // TODO: for simplex or mixed meshes: nothing to do
519  if ((dim == 3 && line_to_cells.size() == 0) ||
520  (cell->reference_cell().is_hyper_cube() == false))
522 
523  if (cell->level() == 0)
525 
526  const std::uint16_t subcell =
527  cell->parent()->child_iterator_to_index(cell);
528  const std::uint16_t subcell_x = (subcell >> 0) & 1;
529  const std::uint16_t subcell_y = (subcell >> 1) & 1;
530  const std::uint16_t subcell_z = (subcell >> 2) & 1;
531 
532  std::uint16_t face = 0;
533  std::uint16_t edge = 0;
534 
535  for (unsigned int direction = 0; direction < dim; ++direction)
536  {
537  const auto side = (subcell >> direction) & 1U;
538  const auto face_no = direction * 2 + side;
539 
540  // ignore if at boundary
541  if (cell->at_boundary(face_no))
542  continue;
543 
544  const auto &neighbor = cell->neighbor(face_no);
545 
546  // ignore neighbors that are artificial or have the same level or
547  // have children
548  if (neighbor->has_children() || neighbor->is_artificial() ||
549  neighbor->level() == cell->level())
550  continue;
551 
552  // Ignore if the neighbors are FE_Nothing
553  if (neighbor->get_fe().n_dofs_per_cell() == 0)
554  continue;
555 
556  face |= 1 << direction;
557  }
558 
559  if (dim == 3)
560  for (unsigned int direction = 0; direction < dim; ++direction)
561  if (face == 0 || face == (1 << direction))
562  {
563  const unsigned int line_no =
564  direction == 0 ?
565  (local_lines[0][subcell_y == 0][subcell_z == 0]) :
566  (direction == 1 ?
567  (local_lines[1][subcell_x == 0][subcell_z == 0]) :
568  (local_lines[2][subcell_x == 0][subcell_y == 0]));
569 
570  const unsigned int line_index = cell->line(line_no)->index();
571 
572  const auto edge_neighbor =
573  std::find_if(line_to_cells[line_index].begin(),
574  line_to_cells[line_index].end(),
575  [&cell](const auto &edge_neighbor) {
577  &edge_neighbor.first->get_triangulation(),
578  edge_neighbor.first->level(),
579  edge_neighbor.first->index(),
580  &cell->get_dof_handler());
581  return edge_neighbor.first->is_artificial() ==
582  false &&
583  edge_neighbor.first->level() <
584  cell->level() &&
585  dof_cell.get_fe().n_dofs_per_cell() > 0;
586  });
587 
588  if (edge_neighbor == line_to_cells[line_index].end())
589  continue;
590 
591  edge |= 1 << direction;
592  }
593 
594  if ((face == 0) && (edge == 0))
596 
597  const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7));
598 
599  const auto refinement_configuration = static_cast<ConstraintKinds>(
600  inverted_subcell + (face << 3) + (edge << 6));
601  Assert(check(refinement_configuration, dim), ExcInternalError());
602  return refinement_configuration;
603  }
604 
605 
606 
607  template <int dim>
608  template <typename CellIterator>
609  inline void
611  const CellIterator & cell,
612  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
613  const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
614  const std::vector<std::vector<bool>> & supported_components,
615  const ConstraintKinds & refinement_configuration,
616  std::vector<types::global_dof_index> & dof_indices) const
617  {
618  if (std::find(supported_components[cell->active_fe_index()].begin(),
619  supported_components[cell->active_fe_index()].end(),
620  true) ==
621  supported_components[cell->active_fe_index()].end())
622  return;
623 
624  const auto &fe = cell->get_fe();
625 
626  AssertDimension(fe.n_unique_faces(), 1);
627 
628  std::vector<std::vector<unsigned int>>
629  component_to_system_index_face_array(fe.n_components());
630 
631  for (unsigned int i = 0; i < fe.n_dofs_per_face(0); ++i)
632  component_to_system_index_face_array
633  [fe.face_system_to_component_index(i, /*face_no=*/0).first]
634  .push_back(i);
635 
636  std::vector<unsigned int> idx_offset = {0};
637 
638  for (unsigned int base_element_index = 0;
639  base_element_index < cell->get_fe().n_base_elements();
640  ++base_element_index)
641  for (unsigned int c = 0;
642  c < cell->get_fe().element_multiplicity(base_element_index);
643  ++c)
644  idx_offset.push_back(
645  idx_offset.back() +
646  cell->get_fe().base_element(base_element_index).n_dofs_per_cell());
647 
648  std::vector<types::global_dof_index> neighbor_dofs_all(idx_offset.back());
649  std::vector<types::global_dof_index> neighbor_dofs_all_temp(
650  idx_offset.back());
651  std::vector<types::global_dof_index> neighbor_dofs_face(
652  fe.n_dofs_per_face(/*face_no=*/0));
653 
654 
655  const auto get_face_idx = [](const auto n_dofs_1d,
656  const auto face_no,
657  const auto i,
658  const auto j) -> unsigned int {
659  const auto direction = face_no / 2;
660  const auto side = face_no % 2;
661  const auto offset = (side == 1) ? (n_dofs_1d - 1) : 0;
662 
663  if (dim == 2)
664  return (direction == 0) ? (n_dofs_1d * i + offset) :
665  (n_dofs_1d * offset + i);
666  else if (dim == 3)
667  switch (direction)
668  {
669  case 0:
670  return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset;
671  case 1:
672  return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i;
673  case 2:
674  return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j;
675  default:
676  Assert(false, ExcNotImplemented());
677  }
678 
679  Assert(false, ExcNotImplemented());
680 
681  return 0;
682  };
683 
684  const std::uint16_t kind =
685  static_cast<std::uint16_t>(refinement_configuration);
686  const std::uint16_t subcell = (kind >> 0) & 7;
687  const std::uint16_t subcell_x = (subcell >> 0) & 1;
688  const std::uint16_t subcell_y = (subcell >> 1) & 1;
689  const std::uint16_t subcell_z = (subcell >> 2) & 1;
690  const std::uint16_t face = (kind >> 3) & 7;
691  const std::uint16_t edge = (kind >> 6) & 7;
692 
693  for (unsigned int direction = 0; direction < dim; ++direction)
694  if ((face >> direction) & 1U)
695  {
696  const auto side = ((subcell >> direction) & 1U) == 0;
697  const auto face_no = direction * 2 + side;
698 
699  // read DoFs of parent of face, ...
700  cell->neighbor(face_no)
701  ->face(cell->neighbor_face_no(face_no))
702  ->get_dof_indices(neighbor_dofs_face,
703  cell->neighbor(face_no)->active_fe_index());
704 
705  // ... convert the global DoFs to serial ones, and ...
706  if (partitioner)
707  for (auto &index : neighbor_dofs_face)
708  index = partitioner->global_to_local(index);
709 
710  for (unsigned int base_element_index = 0, comp = 0;
711  base_element_index < cell->get_fe().n_base_elements();
712  ++base_element_index)
713  for (unsigned int c = 0;
714  c < cell->get_fe().element_multiplicity(base_element_index);
715  ++c, ++comp)
716  {
717  if (supported_components[cell->active_fe_index()][comp] ==
718  false)
719  continue;
720 
721  const unsigned int n_dofs_1d =
722  cell->get_fe()
723  .base_element(base_element_index)
724  .tensor_degree() +
725  1;
726  const unsigned int dofs_per_face =
727  Utilities::pow(n_dofs_1d, dim - 1);
728  std::vector<types::global_dof_index> neighbor_dofs(
729  dofs_per_face);
730  const auto lex_face_mapping =
732  n_dofs_1d - 1);
733 
734  // ... extract the DoFs of the current component
735  for (unsigned int i = 0; i < dofs_per_face; ++i)
736  neighbor_dofs[i] = neighbor_dofs_face
737  [component_to_system_index_face_array[comp][i]];
738 
739  // fix DoFs depending on orientation, flip, and rotation
740  if (dim == 2)
741  {
742  // TODO: for mixed meshes we need to take care of
743  // orientation here
744  Assert(cell->face_orientation(face_no),
746  }
747  else if (dim == 3)
748  {
749  int rotate = 0; // TODO
750  if (cell->face_rotation(face_no)) //
751  rotate -= 1; //
752  if (cell->face_flip(face_no)) //
753  rotate -= 2; //
754 
755  rotate_face(rotate, n_dofs_1d, neighbor_dofs);
756 
757  if (cell->face_orientation(face_no) == false)
758  transpose_face(n_dofs_1d - 1, neighbor_dofs);
759  }
760  else
761  {
762  Assert(false, ExcNotImplemented());
763  }
764 
765  // update DoF map
766  for (unsigned int i = 0, k = 0; i < n_dofs_1d; ++i)
767  for (unsigned int j = 0; j < (dim == 2 ? 1 : n_dofs_1d);
768  ++j, ++k)
769  dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) +
770  idx_offset[comp]] =
771  neighbor_dofs[lex_face_mapping[k]];
772  }
773  }
774 
775  if (dim == 3)
776  for (unsigned int direction = 0; direction < dim; ++direction)
777  if ((edge >> direction) & 1U)
778  {
779  const unsigned int line_no =
780  direction == 0 ?
781  (local_lines[0][subcell_y][subcell_z]) :
782  (direction == 1 ? (local_lines[1][subcell_x][subcell_z]) :
783  (local_lines[2][subcell_x][subcell_y]));
784 
785  const unsigned int line_index = cell->line(line_no)->index();
786 
787  const auto edge_neighbor =
788  std::find_if(line_to_cells[line_index].begin(),
789  line_to_cells[line_index].end(),
790  [&cell](const auto &edge_neighbor) {
791  return edge_neighbor.first->is_artificial() ==
792  false &&
793  edge_neighbor.first->level() <
794  cell->level();
795  });
796 
797  if (edge_neighbor == line_to_cells[line_index].end())
798  continue;
799 
800  const auto neighbor_cell = edge_neighbor->first;
801  const auto local_line_neighbor = edge_neighbor->second;
802 
804  &neighbor_cell->get_triangulation(),
805  neighbor_cell->level(),
806  neighbor_cell->index(),
807  &cell->get_dof_handler())
808  .get_dof_indices(neighbor_dofs_all);
809 
810  if (partitioner)
811  for (auto &index : neighbor_dofs_all)
812  index = partitioner->global_to_local(index);
813 
814  for (unsigned int i = 0; i < neighbor_dofs_all_temp.size(); ++i)
815  neighbor_dofs_all_temp[i] = neighbor_dofs_all
816  [lexicographic_mapping[cell->active_fe_index()][i]];
817 
818  const bool flipped =
819  cell->line_orientation(line_no) !=
820  neighbor_cell->line_orientation(local_line_neighbor);
821 
822  for (unsigned int base_element_index = 0, comp = 0;
823  base_element_index < cell->get_fe().n_base_elements();
824  ++base_element_index)
825  for (unsigned int c = 0;
826  c <
827  cell->get_fe().element_multiplicity(base_element_index);
828  ++c, ++comp)
829  {
830  if (supported_components[cell->active_fe_index()][comp] ==
831  false)
832  continue;
833 
834  const unsigned int n_dofs_1d =
835  cell->get_fe()
836  .base_element(base_element_index)
837  .tensor_degree() +
838  1;
839 
840  for (unsigned int i = 0; i < n_dofs_1d; ++i)
841  dof_indices[line_dof_idx(line_no, i, n_dofs_1d) +
842  idx_offset[comp]] = neighbor_dofs_all_temp
843  [line_dof_idx(local_line_neighbor,
844  flipped ? (n_dofs_1d - 1 - i) : i,
845  n_dofs_1d) +
846  idx_offset[comp]];
847  }
848  }
849  }
850 
851 
852 
853  template <int dim>
854  template <typename CellIterator>
855  inline bool
857  const CellIterator & cell,
858  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
859  const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
860  std::vector<types::global_dof_index> & dof_indices,
861  const ArrayView<ConstraintKinds> & masks) const
862  {
863  // 1) check if finite elements support fast hanging-node algorithm
864  const auto supported_components = compute_supported_components(
865  cell->get_dof_handler().get_fe_collection());
866 
867  if ([](const auto &supported_components) {
868  return std::none_of(supported_components.begin(),
869  supported_components.end(),
870  [](const auto &a) {
871  return *std::max_element(a.begin(), a.end());
872  });
873  }(supported_components))
874  return false;
875 
876  // 2) determine the refinement configuration of the cell
877  const auto refinement_configuration =
878  compute_refinement_configuration(cell);
879 
880  if (refinement_configuration == ConstraintKinds::unconstrained)
881  return false;
882 
883  // 3) update DoF indices of cell for specified components
884  update_dof_indices(cell,
885  partitioner,
886  lexicographic_mapping,
887  supported_components,
888  refinement_configuration,
889  dof_indices);
890 
891  // 4) TODO: copy refinement configuration to all components
892  for (unsigned int c = 0; c < supported_components[0].size(); ++c)
893  if (supported_components[cell->active_fe_index()][c])
894  masks[c] = refinement_configuration;
895 
896  return true;
897  }
898 
899 
900 
901  template <int dim>
902  inline void
904  unsigned int &subface_index) const
905  {
906  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
907 
908  times = times % 4;
909  times = times < 0 ? times + 4 : times;
910  for (int t = 0; t < times; ++t)
911  subface_index = rot_mapping[subface_index];
912  }
913 
914 
915 
916  template <int dim>
917  inline void
919  int times,
920  unsigned int n_dofs_1d,
921  std::vector<types::global_dof_index> &dofs) const
922  {
923  const unsigned int rot_mapping[4] = {2, 0, 3, 1};
924 
925  times = times % 4;
926  times = times < 0 ? times + 4 : times;
927 
928  std::vector<types::global_dof_index> copy(dofs.size());
929  for (int t = 0; t < times; ++t)
930  {
931  std::swap(copy, dofs);
932 
933  // Vertices
934  for (unsigned int i = 0; i < 4; ++i)
935  dofs[rot_mapping[i]] = copy[i];
936 
937  // Edges
938  const unsigned int n_int = n_dofs_1d - 2;
939  unsigned int offset = 4;
940  for (unsigned int i = 0; i < n_int; ++i)
941  {
942  // Left edge
943  dofs[offset + i] = copy[offset + 2 * n_int + (n_int - 1 - i)];
944  // Right edge
945  dofs[offset + n_int + i] =
946  copy[offset + 3 * n_int + (n_int - 1 - i)];
947  // Bottom edge
948  dofs[offset + 2 * n_int + i] = copy[offset + n_int + i];
949  // Top edge
950  dofs[offset + 3 * n_int + i] = copy[offset + i];
951  }
952 
953  // Interior points
954  offset += 4 * n_int;
955 
956  for (unsigned int i = 0; i < n_int; ++i)
957  for (unsigned int j = 0; j < n_int; ++j)
958  dofs[offset + i * n_int + j] =
959  copy[offset + j * n_int + (n_int - 1 - i)];
960  }
961  }
962 
963 
964 
965  template <int dim>
966  inline unsigned int
968  unsigned int dof,
969  unsigned int n_dofs_1d) const
970  {
971  unsigned int x, y, z;
972 
973  const unsigned int fe_degree = n_dofs_1d - 1;
974 
975  if (local_line < 8)
976  {
977  x = (local_line % 4 == 0) ? 0 :
978  (local_line % 4 == 1) ? fe_degree :
979  dof;
980  y = (local_line % 4 == 2) ? 0 :
981  (local_line % 4 == 3) ? fe_degree :
982  dof;
983  z = (local_line / 4) * fe_degree;
984  }
985  else
986  {
987  x = ((local_line - 8) % 2) * fe_degree;
988  y = ((local_line - 8) / 2) * fe_degree;
989  z = dof;
990  }
991 
992  return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
993  }
994 
995 
996 
997  template <int dim>
998  inline void
1000  const unsigned int fe_degree,
1001  std::vector<types::global_dof_index> &dofs) const
1002  {
1003  const std::vector<types::global_dof_index> copy(dofs);
1004 
1005  // Vertices
1006  dofs[1] = copy[2];
1007  dofs[2] = copy[1];
1008 
1009  // Edges
1010  const unsigned int n_int = fe_degree - 1;
1011  unsigned int offset = 4;
1012  for (unsigned int i = 0; i < n_int; ++i)
1013  {
1014  // Right edge
1015  dofs[offset + i] = copy[offset + 2 * n_int + i];
1016  // Left edge
1017  dofs[offset + n_int + i] = copy[offset + 3 * n_int + i];
1018  // Bottom edge
1019  dofs[offset + 2 * n_int + i] = copy[offset + i];
1020  // Top edge
1021  dofs[offset + 3 * n_int + i] = copy[offset + n_int + i];
1022  }
1023 
1024  // Interior
1025  offset += 4 * n_int;
1026  for (unsigned int i = 0; i < n_int; ++i)
1027  for (unsigned int j = 0; j < n_int; ++j)
1028  dofs[offset + i * n_int + j] = copy[offset + j * n_int + i];
1029  }
1030 
1031 
1032 
1033  template <int dim>
1034  void
1035  HangingNodes<dim>::transpose_subface_index(unsigned int &subface) const
1036  {
1037  if (subface == 1)
1038  subface = 2;
1039  else if (subface == 2)
1040  subface = 1;
1041  }
1042  } // namespace MatrixFreeFunctions
1043 } // namespace internal
1044 
1045 
1047 
1048 #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
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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
void rotate(const double angle, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2141
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
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: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_HOST_DEVICE
Definition: numbers.h:35
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation