16 #ifndef dealii_hanging_nodes_internal_h
17 #define dealii_hanging_nodes_internal_h
36 namespace MatrixFreeFunctions
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;
108 if (subcell == 0 && face == 0)
115 if (subcell == 0 && face == 0 && edge == 0)
117 else if (0 < face && edge == 0)
119 else if (0 == face && 0 < edge)
121 else if ((face == edge) && (face == 1 || face == 2 || face == 4))
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;
150 return subcell + ((face > 0) << 3) + ((edge > 0) << 4) +
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;
174 subcell + (((flag_0 & 0b01) ? flag_1 : 0) << 3) +
175 (((flag_0 & 0b10) ? flag_1 : 0) << 6));
206 static_cast<std::uint16_t
>(f2));
230 return static_cast<std::uint16_t
>(f1) !=
static_cast<std::uint16_t
>(f2);
242 return static_cast<std::uint16_t
>(f1) <
static_cast<std::uint16_t
>(f2);
254 static_cast<std::uint16_t
>(f2));
283 template <
typename CellIterator>
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,
296 static std::vector<std::vector<bool>>
302 template <
typename CellIterator>
310 template <
typename CellIterator>
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,
318 std::vector<types::global_dof_index> & dof_indices)
const;
332 unsigned int n_dofs_1d,
333 std::vector<types::global_dof_index> &dofs)
const;
338 unsigned int n_dofs_1d)
const;
342 std::vector<types::global_dof_index> &dofs)
const;
347 std::vector<std::vector<
348 std::pair<typename Triangulation<dim>::cell_iterator,
unsigned int>>>
352 {{{{{7, 3}}, {{6, 2}}}},
353 {{{{5, 1}}, {{4, 0}}}},
354 {{{{11, 9}}, {{10, 8}}}}}};
401 return !cell.is_artificial();
405 const unsigned int n_raw_lines =
triangulation.n_raw_lines();
406 this->line_to_cells.resize(n_raw_lines);
414 const unsigned int line_to_children[12][2] = {{0, 2},
427 std::vector<std::vector<
428 std::pair<typename Triangulation<3>::cell_iterator,
unsigned int>>>
429 line_to_inactive_cells(n_raw_lines);
434 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
437 const unsigned int line_idx = cell->
line(line)->index();
439 line_to_cells[line_idx].push_back(std::make_pair(cell, line));
441 line_to_inactive_cells[line_idx].push_back(
442 std::make_pair(cell, line));
449 for (
unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
451 if ((line_to_cells[line_idx].size() > 0) &&
452 line_to_inactive_cells[line_idx].size() > 0)
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;
461 for (
unsigned int c = 0; c < 2; ++c)
464 inactive_cell->child(line_to_children[neighbor_line][c]);
465 const unsigned int child_line_idx =
466 child->line(neighbor_line)->index();
469 for (
const auto &cl : line_to_cells[line_idx])
470 line_to_cells[child_line_idx].push_back(cl);
479 inline std::vector<std::vector<bool>>
481 const ::hp::FECollection<dim> &fe_collection)
483 std::vector<std::vector<bool>> supported_components(
484 fe_collection.size(),
485 std::vector<bool>(fe_collection.n_components(),
false));
487 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
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);
497 &fe_collection[i].base_element(base_element_index)) ==
500 &fe_collection[i].base_element(base_element_index)) ==
502 supported_components[i][comp] =
false;
504 supported_components[i][comp] =
true;
507 return supported_components;
513 template <
typename CellIterator>
516 const CellIterator &cell)
const
519 if ((dim == 3 && line_to_cells.size() == 0) ||
520 (cell->reference_cell().is_hyper_cube() ==
false))
523 if (cell->level() == 0)
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;
532 std::uint16_t face = 0;
533 std::uint16_t edge = 0;
535 for (
unsigned int direction = 0; direction < dim; ++direction)
537 const auto side = (subcell >> direction) & 1U;
538 const auto face_no = direction * 2 + side;
541 if (cell->at_boundary(face_no))
544 const auto &neighbor = cell->neighbor(face_no);
548 if (neighbor->has_children() || neighbor->is_artificial() ||
549 neighbor->level() == cell->level())
553 if (neighbor->get_fe().n_dofs_per_cell() == 0)
556 face |= 1 << direction;
560 for (
unsigned int direction = 0; direction < dim; ++direction)
561 if (face == 0 || face == (1 << direction))
563 const unsigned int line_no =
570 const unsigned int line_index = cell->line(line_no)->index();
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() ==
583 edge_neighbor.first->level() <
585 dof_cell.
get_fe().n_dofs_per_cell() > 0;
588 if (edge_neighbor == line_to_cells[line_index].
end())
591 edge |= 1 << direction;
594 if ((face == 0) && (edge == 0))
597 const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7));
600 inverted_subcell + (face << 3) + (edge << 6));
602 return refinement_configuration;
608 template <
typename CellIterator>
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,
616 std::vector<types::global_dof_index> & dof_indices)
const
618 if (std::find(supported_components[cell->active_fe_index()].begin(),
619 supported_components[cell->active_fe_index()].end(),
621 supported_components[cell->active_fe_index()].end())
624 const auto &fe = cell->get_fe();
628 std::vector<std::vector<unsigned int>>
629 component_to_system_index_face_array(fe.n_components());
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, 0).first]
636 std::vector<unsigned int> idx_offset = {0};
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);
644 idx_offset.push_back(
646 cell->get_fe().base_element(base_element_index).n_dofs_per_cell());
648 std::vector<types::global_dof_index> neighbor_dofs_all(idx_offset.back());
649 std::vector<types::global_dof_index> neighbor_dofs_all_temp(
651 std::vector<types::global_dof_index> neighbor_dofs_face(
652 fe.n_dofs_per_face(0));
655 const auto get_face_idx = [](
const auto n_dofs_1d,
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;
664 return (direction == 0) ? (n_dofs_1d * i + offset) :
665 (n_dofs_1d * offset + i);
670 return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset;
672 return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i;
674 return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j;
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;
693 for (
unsigned int direction = 0; direction < dim; ++direction)
694 if ((face >> direction) & 1U)
696 const auto side = ((subcell >> direction) & 1U) == 0;
697 const auto face_no = direction * 2 + side;
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());
707 for (
auto &
index : neighbor_dofs_face)
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);
717 if (supported_components[cell->active_fe_index()][comp] ==
721 const unsigned int n_dofs_1d =
723 .base_element(base_element_index)
726 const unsigned int dofs_per_face =
728 std::vector<types::global_dof_index> neighbor_dofs(
730 const auto lex_face_mapping =
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]];
744 Assert(cell->face_orientation(face_no),
750 if (cell->face_rotation(face_no))
752 if (cell->face_flip(face_no))
755 rotate_face(
rotate, n_dofs_1d, neighbor_dofs);
757 if (cell->face_orientation(face_no) ==
false)
758 transpose_face(n_dofs_1d - 1, neighbor_dofs);
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);
769 dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) +
771 neighbor_dofs[lex_face_mapping[k]];
776 for (
unsigned int direction = 0; direction < dim; ++direction)
777 if ((edge >> direction) & 1U)
779 const unsigned int line_no =
785 const unsigned int line_index = cell->line(line_no)->index();
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() ==
793 edge_neighbor.first->level() <
797 if (edge_neighbor == line_to_cells[line_index].
end())
800 const auto neighbor_cell = edge_neighbor->first;
801 const auto local_line_neighbor = edge_neighbor->second;
804 &neighbor_cell->get_triangulation(),
805 neighbor_cell->level(),
806 neighbor_cell->index(),
807 &cell->get_dof_handler())
811 for (
auto &
index : neighbor_dofs_all)
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]];
819 cell->line_orientation(line_no) !=
820 neighbor_cell->line_orientation(local_line_neighbor);
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;
827 cell->get_fe().element_multiplicity(base_element_index);
830 if (supported_components[cell->active_fe_index()][comp] ==
834 const unsigned int n_dofs_1d =
836 .base_element(base_element_index)
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,
854 template <
typename CellIterator>
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,
864 const auto supported_components = compute_supported_components(
865 cell->get_dof_handler().get_fe_collection());
867 if ([](
const auto &supported_components) {
868 return std::none_of(supported_components.begin(),
869 supported_components.end(),
871 return *std::max_element(a.begin(), a.end());
873 }(supported_components))
877 const auto refinement_configuration =
878 compute_refinement_configuration(cell);
884 update_dof_indices(cell,
886 lexicographic_mapping,
887 supported_components,
888 refinement_configuration,
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;
904 unsigned int &subface_index)
const
906 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
909 times = times < 0 ? times + 4 : times;
910 for (
int t = 0; t < times; ++t)
911 subface_index = rot_mapping[subface_index];
920 unsigned int n_dofs_1d,
921 std::vector<types::global_dof_index> &dofs)
const
923 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
926 times = times < 0 ? times + 4 : times;
928 std::vector<types::global_dof_index>
copy(dofs.size());
929 for (
int t = 0; t < times; ++t)
934 for (
unsigned int i = 0; i < 4; ++i)
935 dofs[rot_mapping[i]] =
copy[i];
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)
943 dofs[offset + i] =
copy[offset + 2 * n_int + (n_int - 1 - i)];
945 dofs[offset + n_int + i] =
946 copy[offset + 3 * n_int + (n_int - 1 - i)];
948 dofs[offset + 2 * n_int + i] =
copy[offset + n_int + i];
950 dofs[offset + 3 * n_int + i] =
copy[offset + i];
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)];
969 unsigned int n_dofs_1d)
const
971 unsigned int x, y, z;
973 const unsigned int fe_degree = n_dofs_1d - 1;
977 x = (local_line % 4 == 0) ? 0 :
978 (local_line % 4 == 1) ? fe_degree :
980 y = (local_line % 4 == 2) ? 0 :
981 (local_line % 4 == 3) ? fe_degree :
983 z = (local_line / 4) * fe_degree;
987 x = ((local_line - 8) % 2) * fe_degree;
988 y = ((local_line - 8) / 2) * fe_degree;
992 return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
1000 const unsigned int fe_degree,
1001 std::vector<types::global_dof_index> &dofs)
const
1003 const std::vector<types::global_dof_index>
copy(dofs);
1010 const unsigned int n_int = fe_degree - 1;
1011 unsigned int offset = 4;
1012 for (
unsigned int i = 0; i < n_int; ++i)
1015 dofs[offset + i] =
copy[offset + 2 * n_int + i];
1017 dofs[offset + n_int + i] =
copy[offset + 3 * n_int + i];
1019 dofs[offset + 2 * n_int + i] =
copy[offset + i];
1021 dofs[offset + 3 * n_int + i] =
copy[offset + n_int + i];
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];
1039 else if (subface == 2)
const FiniteElement< dimension_, space_dimension_ > & get_fe() const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
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
std::size_t memory_consumption() 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
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
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)
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
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
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation