Reference documentation for deal.II version 9.6.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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_hanging_nodes_internal_h
16#define dealii_hanging_nodes_internal_h
17
18#include <deal.II/base/config.h>
19
22
24
25#include <deal.II/fe/fe_q.h>
27#include <deal.II/fe/fe_tools.h>
28
30
32
33namespace 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
83 constexpr compressed_constraint_kind
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
133 inline compressed_constraint_kind
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
201 DEAL_II_HOST_DEVICE inline ConstraintKinds
203 {
204 return static_cast<ConstraintKinds>(static_cast<std::uint16_t>(f1) |
205 static_cast<std::uint16_t>(f2));
206 }
207
208
209
214 DEAL_II_HOST_DEVICE inline ConstraintKinds &
216 {
217 f1 = f1 | f2;
218 return f1;
219 }
220
221
222
226 DEAL_II_HOST_DEVICE 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
249 DEAL_II_HOST_DEVICE inline ConstraintKinds
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 std::size_t
277 memory_consumption() const;
278
282 template <typename CellIterator>
283 bool
285 const CellIterator &cell,
286 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
287 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
288 std::vector<types::global_dof_index> &dof_indices,
289 const ArrayView<ConstraintKinds> &mask) const;
290
295 static std::vector<std::vector<bool>>
296 compute_supported_components(const ::hp::FECollection<dim> &fe);
297
301 template <typename CellIterator>
303 compute_refinement_configuration(const CellIterator &cell) const;
304
309 template <typename CellIterator>
310 void
312 const CellIterator &cell,
313 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
314 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
315 const std::vector<std::vector<bool>> &component_mask,
316 const ConstraintKinds &refinement_configuration,
317 std::vector<types::global_dof_index> &dof_indices) const;
318
319 private:
323 void
325
326 void
327 rotate_subface_index(int times, unsigned int &subface_index) const;
328
329 void
330 rotate_face(int times,
331 unsigned int n_dofs_1d,
332 std::vector<types::global_dof_index> &dofs) const;
333
334 unsigned int
335 line_dof_idx(int local_line,
336 unsigned int dof,
337 unsigned int n_dofs_1d) const;
338
339 void
340 transpose_face(const unsigned int fe_degree,
341 std::vector<types::global_dof_index> &dofs) const;
342
343 void
344 transpose_subface_index(unsigned int &subface) const;
345
346 std::vector<
347 boost::container::small_vector<std::array<unsigned int, 3>, 6>>
349
350 const ::ndarray<unsigned int, 3, 2, 2> local_lines = {
351 {{{{{7, 3}}, {{6, 2}}}},
352 {{{{5, 1}}, {{4, 0}}}},
353 {{{{11, 9}}, {{10, 8}}}}}};
354 };
355
356
357
358 template <int dim>
361 {
362 // Set up line-to-cell mapping for edge constraints (only if dim = 3 and
363 // for pure hex meshes)
365 setup_line_to_cell(triangulation);
366 }
367
368
369
370 template <int dim>
371 inline std::size_t
373 {
374 std::size_t size = 0;
375 for (const auto &a : line_to_cells)
376 size +=
377 (a.capacity() > 6 ? a.capacity() : 0) * sizeof(a[0]) + sizeof(a);
378 return size;
379 }
380
381
382
383 template <int dim>
384 inline void
390
391
392
393 template <>
394 inline void
396 {
397 // Check if we there are no hanging nodes on the current MPI process,
398 // which we do by checking if the second finest level holds no active
399 // non-artificial cell
400 if (triangulation.n_levels() <= 1 ||
403 [](const CellAccessor<3, 3> &cell) {
404 return !cell.is_artificial();
405 }))
406 return;
407
408 const unsigned int n_raw_lines = triangulation.n_raw_lines();
409 this->line_to_cells.resize(n_raw_lines);
410
411 // In 3d, we can have DoFs on only an edge being constrained (e.g. in a
412 // cartesian 2x2x2 grid, where only the upper left 2 cells are refined).
413 // This sets up a helper data structure in the form of a mapping from
414 // edges (i.e. lines) to neighboring cells.
415
416 // Mapping from an edge to which children that share that edge.
417 const unsigned int line_to_children[12][2] = {{0, 2},
418 {1, 3},
419 {0, 1},
420 {2, 3},
421 {4, 6},
422 {5, 7},
423 {4, 5},
424 {6, 7},
425 {0, 4},
426 {1, 5},
427 {2, 6},
428 {3, 7}};
429
430 std::vector<
431 boost::container::small_vector<std::array<unsigned int, 3>, 6>>
432 line_to_inactive_cells(n_raw_lines);
433
434 // First add active and inactive cells to their lines:
435 for (const auto &cell : triangulation.cell_iterators())
436 {
437 const unsigned int cell_level = cell->level();
438 const unsigned int cell_index = cell->index();
439 for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
440 ++line)
441 {
442 const unsigned int line_idx = cell->line_index(line);
443 if (cell->is_active())
444 line_to_cells[line_idx].push_back(
445 {{cell_level, cell_index, line}});
446 else
447 line_to_inactive_cells[line_idx].push_back(
448 {{cell_level, cell_index, line}});
449 }
450 }
451
452 // Now, we can access edge-neighboring active cells on same level to also
453 // access of an edge to the edges "children". These are found from looking
454 // at the corresponding edge of children of inactive edge neighbors.
455 for (unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
456 {
457 if ((line_to_cells[line_idx].size() > 0) &&
458 line_to_inactive_cells[line_idx].size() > 0)
459 {
460 // We now have cells to add (active ones) and edges to which they
461 // should be added (inactive cells).
462 const Triangulation<3>::cell_iterator inactive_cell(
464 line_to_inactive_cells[line_idx][0][0],
465 line_to_inactive_cells[line_idx][0][1]);
466 const unsigned int neighbor_line =
467 line_to_inactive_cells[line_idx][0][2];
468
469 for (unsigned int c = 0; c < 2; ++c)
470 {
471 const auto &child =
472 inactive_cell->child(line_to_children[neighbor_line][c]);
473 const unsigned int child_line_idx =
474 child->line_index(neighbor_line);
475
476 Assert(child->is_active(), ExcInternalError());
477
478 // Now add all active cells
479 for (const auto &cl : line_to_cells[line_idx])
480 line_to_cells[child_line_idx].push_back(cl);
481 }
482 }
483 }
484 }
485
486
487
488 template <int dim>
489 inline std::vector<std::vector<bool>>
491 const ::hp::FECollection<dim> &fe_collection)
492 {
493 std::vector<std::vector<bool>> supported_components(
494 fe_collection.size(),
495 std::vector<bool>(fe_collection.n_components(), false));
496
497 for (unsigned int i = 0; i < fe_collection.size(); ++i)
498 {
499 for (unsigned int base_element_index = 0, comp = 0;
500 base_element_index < fe_collection[i].n_base_elements();
501 ++base_element_index)
502 for (unsigned int c = 0;
503 c < fe_collection[i].element_multiplicity(base_element_index);
504 ++c, ++comp)
505 if (dim == 1 ||
506 (dynamic_cast<const FE_Q<dim> *>(
507 &fe_collection[i].base_element(base_element_index)) ==
508 nullptr &&
509 dynamic_cast<const FE_Q_iso_Q1<dim> *>(
510 &fe_collection[i].base_element(base_element_index)) ==
511 nullptr))
512 supported_components[i][comp] = false;
513 else
514 supported_components[i][comp] = true;
515 }
516
517 return supported_components;
518 }
519
520
521
522 template <int dim>
523 template <typename CellIterator>
524 inline ConstraintKinds
526 const CellIterator &cell) const
527 {
528 // TODO: for simplex or mixed meshes: nothing to do
529 if ((dim == 3 && line_to_cells.empty()) ||
530 (cell->reference_cell().is_hyper_cube() == false))
532
533 if (cell->level() == 0)
535
536 const std::uint16_t subcell =
537 cell->parent()->child_iterator_to_index(cell);
538 const std::uint16_t subcell_x = (subcell >> 0) & 1;
539 const std::uint16_t subcell_y = (subcell >> 1) & 1;
540 const std::uint16_t subcell_z = (subcell >> 2) & 1;
541
542 std::uint16_t face = 0;
543 std::uint16_t edge = 0;
544
545 for (unsigned int direction = 0; direction < dim; ++direction)
546 {
547 const auto side = (subcell >> direction) & 1U;
548 const auto face_no = direction * 2 + side;
549
550 // ignore if at boundary
551 if (cell->at_boundary(face_no))
552 continue;
553
554 const auto &neighbor = cell->neighbor(face_no);
555
556 // ignore neighbors that are artificial or have the same level or
557 // have children
558 if (neighbor->has_children() || neighbor->is_artificial() ||
559 neighbor->level() == cell->level())
560 continue;
561
562 // Ignore if the neighbors are FE_Nothing
563 if (neighbor->get_fe().n_dofs_per_cell() == 0)
564 continue;
565
566 face |= 1 << direction;
567 }
568
569 if (dim == 3)
570 for (unsigned int direction = 0; direction < dim; ++direction)
571 if (face == 0 || face == (1 << direction))
572 {
573 const unsigned int line_no =
574 direction == 0 ?
575 (local_lines[0][subcell_y == 0][subcell_z == 0]) :
576 (direction == 1 ?
577 (local_lines[1][subcell_x == 0][subcell_z == 0]) :
578 (local_lines[2][subcell_x == 0][subcell_y == 0]));
579
580 const unsigned int line_index = cell->line_index(line_no);
581
582 const auto edge_neighbor =
583 std::find_if(line_to_cells[line_index].begin(),
584 line_to_cells[line_index].end(),
585 [&cell](const auto &edge_neighbor) {
587 &cell->get_triangulation(),
588 edge_neighbor[0],
589 edge_neighbor[1],
590 &cell->get_dof_handler());
591 return dof_cell.is_artificial() == false &&
592 dof_cell.level() < cell->level() &&
593 dof_cell.get_fe().n_dofs_per_cell() > 0;
594 });
595
596 if (edge_neighbor == line_to_cells[line_index].end())
597 continue;
598
599 edge |= 1 << direction;
600 }
601
602 if ((face == 0) && (edge == 0))
604
605 const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7));
606
607 const auto refinement_configuration = static_cast<ConstraintKinds>(
608 inverted_subcell + (face << 3) + (edge << 6));
609 Assert(check(refinement_configuration, dim), ExcInternalError());
610 return refinement_configuration;
611 }
612
613
614
615 template <int dim>
616 template <typename CellIterator>
617 inline void
619 const CellIterator &cell,
620 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
621 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
622 const std::vector<std::vector<bool>> &supported_components,
623 const ConstraintKinds &refinement_configuration,
624 std::vector<types::global_dof_index> &dof_indices) const
625 {
626 if (std::find(supported_components[cell->active_fe_index()].begin(),
627 supported_components[cell->active_fe_index()].end(),
628 true) ==
629 supported_components[cell->active_fe_index()].end())
630 return;
631
632 const auto &fe = cell->get_fe();
633
634 AssertDimension(fe.n_unique_faces(), 1);
635
636 std::vector<std::vector<unsigned int>>
637 component_to_system_index_face_array(fe.n_components());
638
639 for (unsigned int i = 0; i < fe.n_dofs_per_face(0); ++i)
640 component_to_system_index_face_array
641 [fe.face_system_to_component_index(i, /*face_no=*/0).first]
642 .push_back(i);
643
644 std::vector<unsigned int> idx_offset = {0};
645
646 for (unsigned int base_element_index = 0;
647 base_element_index < cell->get_fe().n_base_elements();
648 ++base_element_index)
649 for (unsigned int c = 0;
650 c < cell->get_fe().element_multiplicity(base_element_index);
651 ++c)
652 idx_offset.push_back(
653 idx_offset.back() +
654 cell->get_fe().base_element(base_element_index).n_dofs_per_cell());
655
656 std::vector<types::global_dof_index> neighbor_dofs_all(idx_offset.back());
657 std::vector<types::global_dof_index> neighbor_dofs_all_temp(
658 idx_offset.back());
659 std::vector<types::global_dof_index> neighbor_dofs_face(
660 fe.n_dofs_per_face(/*face_no=*/0));
661
662
663 const auto get_face_idx = [](const auto n_dofs_1d,
664 const auto face_no,
665 const auto i,
666 const auto j) -> unsigned int {
667 const auto direction = face_no / 2;
668 const auto side = face_no % 2;
669 const auto offset = (side == 1) ? (n_dofs_1d - 1) : 0;
670
671 if (dim == 2)
672 return (direction == 0) ? (n_dofs_1d * i + offset) :
673 (n_dofs_1d * offset + i);
674 else if (dim == 3)
675 switch (direction)
676 {
677 case 0:
678 return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset;
679 case 1:
680 return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i;
681 case 2:
682 return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j;
683 default:
685 }
686
688
689 return 0;
690 };
691
692 const std::uint16_t kind =
693 static_cast<std::uint16_t>(refinement_configuration);
694 const std::uint16_t subcell = (kind >> 0) & 7;
695 const std::uint16_t subcell_x = (subcell >> 0) & 1;
696 const std::uint16_t subcell_y = (subcell >> 1) & 1;
697 const std::uint16_t subcell_z = (subcell >> 2) & 1;
698 const std::uint16_t face = (kind >> 3) & 7;
699 const std::uint16_t edge = (kind >> 6) & 7;
700
701 for (unsigned int direction = 0; direction < dim; ++direction)
702 if ((face >> direction) & 1U)
703 {
704 const auto side = ((subcell >> direction) & 1U) == 0;
705 const auto face_no = direction * 2 + side;
706
707 // read DoFs of parent of face, ...
708 cell->neighbor(face_no)
709 ->face(cell->neighbor_face_no(face_no))
710 ->get_dof_indices(neighbor_dofs_face,
711 cell->neighbor(face_no)->active_fe_index());
712
713 // ... convert the global DoFs to serial ones, and ...
714 if (partitioner)
715 for (auto &index : neighbor_dofs_face)
716 index = partitioner->global_to_local(index);
717
718 for (unsigned int base_element_index = 0, comp = 0;
719 base_element_index < cell->get_fe().n_base_elements();
720 ++base_element_index)
721 for (unsigned int c = 0;
722 c < cell->get_fe().element_multiplicity(base_element_index);
723 ++c, ++comp)
724 {
725 if (supported_components[cell->active_fe_index()][comp] ==
726 false)
727 continue;
728
729 const unsigned int n_dofs_1d =
730 cell->get_fe()
731 .base_element(base_element_index)
732 .tensor_degree() +
733 1;
734 const unsigned int dofs_per_face =
735 Utilities::pow(n_dofs_1d, dim - 1);
736 std::vector<types::global_dof_index> neighbor_dofs(
737 dofs_per_face);
738 const auto lex_face_mapping =
740 n_dofs_1d - 1);
741
742 // ... extract the DoFs of the current component
743 for (unsigned int i = 0; i < dofs_per_face; ++i)
744 neighbor_dofs[i] = neighbor_dofs_face
745 [component_to_system_index_face_array[comp][i]];
746
747 // fix DoFs depending on orientation, flip, and rotation
748 if (dim == 2)
749 {
750 // TODO: for mixed meshes we need to take care of
751 // orientation here
752 Assert(cell->face_orientation(face_no),
754 }
755 else if (dim == 3)
756 {
757 int rotate = 0; // TODO
758 if (cell->face_rotation(face_no)) //
759 rotate -= 1; //
760 if (cell->face_flip(face_no)) //
761 rotate -= 2; //
762
763 rotate_face(rotate, n_dofs_1d, neighbor_dofs);
764
765 if (cell->face_orientation(face_no) == false)
766 transpose_face(n_dofs_1d - 1, neighbor_dofs);
767 }
768 else
769 {
771 }
772
773 // update DoF map
774 for (unsigned int i = 0, k = 0; i < n_dofs_1d; ++i)
775 for (unsigned int j = 0; j < (dim == 2 ? 1 : n_dofs_1d);
776 ++j, ++k)
777 dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) +
778 idx_offset[comp]] =
779 neighbor_dofs[lex_face_mapping[k]];
780 }
781 }
782
783 if (dim == 3)
784 for (unsigned int direction = 0; direction < dim; ++direction)
785 if ((edge >> direction) & 1U)
786 {
787 const unsigned int line_no =
788 direction == 0 ?
789 (local_lines[0][subcell_y][subcell_z]) :
790 (direction == 1 ? (local_lines[1][subcell_x][subcell_z]) :
791 (local_lines[2][subcell_x][subcell_y]));
792
793 const unsigned int line_index = cell->line(line_no)->index();
794
795 const auto edge_neighbor =
796 std::find_if(line_to_cells[line_index].begin(),
797 line_to_cells[line_index].end(),
798 [&cell](const auto &edge_array) {
800 edge_neighbor(&cell->get_triangulation(),
801 edge_array[0],
802 edge_array[1]);
803 return edge_neighbor->is_artificial() == false &&
804 edge_neighbor->level() < cell->level();
805 });
806
807 if (edge_neighbor == line_to_cells[line_index].end())
808 continue;
809
810 const DoFCellAccessor<dim, dim, false> neighbor_cell(
811 &cell->get_triangulation(),
812 (*edge_neighbor)[0],
813 (*edge_neighbor)[1],
814 &cell->get_dof_handler());
815 const auto local_line_neighbor = (*edge_neighbor)[2];
816
817 neighbor_cell.get_dof_indices(neighbor_dofs_all);
818
819 if (partitioner)
820 for (auto &index : neighbor_dofs_all)
821 index = partitioner->global_to_local(index);
822
823 for (unsigned int i = 0; i < neighbor_dofs_all_temp.size(); ++i)
824 neighbor_dofs_all_temp[i] = neighbor_dofs_all
825 [lexicographic_mapping[cell->active_fe_index()][i]];
826
827 const bool flipped =
828 cell->line_orientation(line_no) !=
829 neighbor_cell.line_orientation(local_line_neighbor);
830
831 for (unsigned int base_element_index = 0, comp = 0;
832 base_element_index < cell->get_fe().n_base_elements();
833 ++base_element_index)
834 for (unsigned int c = 0;
835 c <
836 cell->get_fe().element_multiplicity(base_element_index);
837 ++c, ++comp)
838 {
839 if (supported_components[cell->active_fe_index()][comp] ==
840 false)
841 continue;
842
843 const unsigned int n_dofs_1d =
844 cell->get_fe()
845 .base_element(base_element_index)
846 .tensor_degree() +
847 1;
848
849 for (unsigned int i = 0; i < n_dofs_1d; ++i)
850 dof_indices[line_dof_idx(line_no, i, n_dofs_1d) +
851 idx_offset[comp]] = neighbor_dofs_all_temp
852 [line_dof_idx(local_line_neighbor,
853 flipped ? (n_dofs_1d - 1 - i) : i,
854 n_dofs_1d) +
855 idx_offset[comp]];
856 }
857 }
858 }
859
860
861
862 template <int dim>
863 template <typename CellIterator>
864 inline bool
866 const CellIterator &cell,
867 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
868 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
869 std::vector<types::global_dof_index> &dof_indices,
870 const ArrayView<ConstraintKinds> &masks) const
871 {
872 // 1) check if finite elements support fast hanging-node algorithm
873 const auto supported_components = compute_supported_components(
874 cell->get_dof_handler().get_fe_collection());
875
876 if (std::none_of(supported_components.begin(),
877 supported_components.end(),
878 [](const auto &a) {
879 return *std::max_element(a.begin(), a.end());
880 }))
881 return false;
882
883 // 2) determine the refinement configuration of the cell
884 const auto refinement_configuration =
885 compute_refinement_configuration(cell);
886
887 if (refinement_configuration == ConstraintKinds::unconstrained)
888 return false;
889
890 // 3) update DoF indices of cell for specified components
891 update_dof_indices(cell,
892 partitioner,
893 lexicographic_mapping,
894 supported_components,
895 refinement_configuration,
896 dof_indices);
897
898 // 4) TODO: copy refinement configuration to all components
899 for (unsigned int c = 0; c < supported_components[0].size(); ++c)
900 if (supported_components[cell->active_fe_index()][c])
901 masks[c] = refinement_configuration;
902
903 return true;
904 }
905
906
907
908 template <int dim>
909 inline void
911 unsigned int &subface_index) const
912 {
913 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
914
915 times = times % 4;
916 times = times < 0 ? times + 4 : times;
917 for (int t = 0; t < times; ++t)
918 subface_index = rot_mapping[subface_index];
919 }
920
921
922
923 template <int dim>
924 inline void
926 int times,
927 unsigned int n_dofs_1d,
928 std::vector<types::global_dof_index> &dofs) const
929 {
930 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
931
932 times = times % 4;
933 times = times < 0 ? times + 4 : times;
934
935 std::vector<types::global_dof_index> copy(dofs.size());
936 for (int t = 0; t < times; ++t)
937 {
938 std::swap(copy, dofs);
939
940 // Vertices
941 for (unsigned int i = 0; i < 4; ++i)
942 dofs[rot_mapping[i]] = copy[i];
943
944 // Edges
945 const unsigned int n_int = n_dofs_1d - 2;
946 unsigned int offset = 4;
947 for (unsigned int i = 0; i < n_int; ++i)
948 {
949 // Left edge
950 dofs[offset + i] = copy[offset + 2 * n_int + (n_int - 1 - i)];
951 // Right edge
952 dofs[offset + n_int + i] =
953 copy[offset + 3 * n_int + (n_int - 1 - i)];
954 // Bottom edge
955 dofs[offset + 2 * n_int + i] = copy[offset + n_int + i];
956 // Top edge
957 dofs[offset + 3 * n_int + i] = copy[offset + i];
958 }
959
960 // Interior points
961 offset += 4 * n_int;
962
963 for (unsigned int i = 0; i < n_int; ++i)
964 for (unsigned int j = 0; j < n_int; ++j)
965 dofs[offset + i * n_int + j] =
966 copy[offset + j * n_int + (n_int - 1 - i)];
967 }
968 }
969
970
971
972 template <int dim>
973 inline unsigned int
975 unsigned int dof,
976 unsigned int n_dofs_1d) const
977 {
978 unsigned int x, y, z;
979
980 const unsigned int fe_degree = n_dofs_1d - 1;
981
982 if (local_line < 8)
983 {
984 x = (local_line % 4 == 0) ? 0 :
985 (local_line % 4 == 1) ? fe_degree :
986 dof;
987 y = (local_line % 4 == 2) ? 0 :
988 (local_line % 4 == 3) ? fe_degree :
989 dof;
990 z = (local_line / 4) * fe_degree;
991 }
992 else
993 {
994 x = ((local_line - 8) % 2) * fe_degree;
995 y = ((local_line - 8) / 2) * fe_degree;
996 z = dof;
997 }
998
999 return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
1000 }
1001
1002
1003
1004 template <int dim>
1005 inline void
1007 const unsigned int fe_degree,
1008 std::vector<types::global_dof_index> &dofs) const
1009 {
1010 const std::vector<types::global_dof_index> copy(dofs);
1011
1012 // Vertices
1013 dofs[1] = copy[2];
1014 dofs[2] = copy[1];
1015
1016 // Edges
1017 const unsigned int n_int = fe_degree - 1;
1018 unsigned int offset = 4;
1019 for (unsigned int i = 0; i < n_int; ++i)
1020 {
1021 // Right edge
1022 dofs[offset + i] = copy[offset + 2 * n_int + i];
1023 // Left edge
1024 dofs[offset + n_int + i] = copy[offset + 3 * n_int + i];
1025 // Bottom edge
1026 dofs[offset + 2 * n_int + i] = copy[offset + i];
1027 // Top edge
1028 dofs[offset + 3 * n_int + i] = copy[offset + n_int + i];
1029 }
1030
1031 // Interior
1032 offset += 4 * n_int;
1033 for (unsigned int i = 0; i < n_int; ++i)
1034 for (unsigned int j = 0; j < n_int; ++j)
1035 dofs[offset + i * n_int + j] = copy[offset + j * n_int + i];
1036 }
1037
1038
1039
1040 template <int dim>
1041 void
1043 {
1044 if (subface == 1)
1045 subface = 2;
1046 else if (subface == 2)
1047 subface = 1;
1048 }
1049 } // namespace MatrixFreeFunctions
1050} // namespace internal
1051
1052
1054
1055#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:554
int index() const
int level() const
unsigned int line_index(const unsigned int i) const
bool all_reference_cells_are_hyper_cube() const
unsigned int n_raw_lines() const
unsigned int n_levels() const
active_cell_iterator end_active(const unsigned int level) const
active_cell_iterator begin_active(const unsigned int level=0) 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)
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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
unsigned int cell_index
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
bool operator<(const ConstraintKinds f1, const ConstraintKinds f2)
ConstraintKinds decompress(const compressed_constraint_kind kind_in, const unsigned int dim)
std::uint8_t compressed_constraint_kind
Definition dof_info.h:86
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 &)
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:34
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation