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