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