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;
348 boost::container::small_vector<std::array<unsigned int, 3>, 6>>
352 {{{{{7, 3}}, {{6, 2}}}},
353 {{{{5, 1}}, {{4, 0}}}},
354 {{{{11, 9}}, {{10, 8}}}}}};
375 std::size_t size = 0;
376 for (
const auto &a : line_to_cells)
378 (a.capacity() > 6 ? a.capacity() : 0) *
sizeof(a[0]) +
sizeof(a);
405 return !cell.is_artificial();
409 const unsigned int n_raw_lines =
triangulation.n_raw_lines();
410 this->line_to_cells.resize(n_raw_lines);
418 const unsigned int line_to_children[12][2] = {{0, 2},
432 boost::container::small_vector<std::array<unsigned int, 3>, 6>>
433 line_to_inactive_cells(n_raw_lines);
438 const unsigned int cell_level = cell->
level();
440 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
443 const unsigned int line_idx = cell->
line_index(line);
445 line_to_cells[line_idx].push_back(
448 line_to_inactive_cells[line_idx].push_back(
456 for (
unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
458 if ((line_to_cells[line_idx].size() > 0) &&
459 line_to_inactive_cells[line_idx].size() > 0)
465 line_to_inactive_cells[line_idx][0][0],
466 line_to_inactive_cells[line_idx][0][1]);
467 const unsigned int neighbor_line =
468 line_to_inactive_cells[line_idx][0][2];
470 for (
unsigned int c = 0; c < 2; ++c)
473 inactive_cell->child(line_to_children[neighbor_line][c]);
474 const unsigned int child_line_idx =
475 child->line_index(neighbor_line);
480 for (
const auto &cl : line_to_cells[line_idx])
481 line_to_cells[child_line_idx].push_back(cl);
490 inline std::vector<std::vector<bool>>
492 const ::hp::FECollection<dim> &fe_collection)
494 std::vector<std::vector<bool>> supported_components(
495 fe_collection.size(),
496 std::vector<bool>(fe_collection.n_components(),
false));
498 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
500 for (
unsigned int base_element_index = 0, comp = 0;
501 base_element_index < fe_collection[i].n_base_elements();
502 ++base_element_index)
503 for (
unsigned int c = 0;
504 c < fe_collection[i].element_multiplicity(base_element_index);
508 &fe_collection[i].base_element(base_element_index)) ==
511 &fe_collection[i].base_element(base_element_index)) ==
513 supported_components[i][comp] =
false;
515 supported_components[i][comp] =
true;
518 return supported_components;
524 template <
typename CellIterator>
527 const CellIterator &cell)
const
530 if ((dim == 3 && line_to_cells.empty()) ||
531 (cell->reference_cell().is_hyper_cube() ==
false))
534 if (cell->level() == 0)
537 const std::uint16_t subcell =
538 cell->parent()->child_iterator_to_index(cell);
539 const std::uint16_t
subcell_x = (subcell >> 0) & 1;
540 const std::uint16_t
subcell_y = (subcell >> 1) & 1;
541 const std::uint16_t
subcell_z = (subcell >> 2) & 1;
543 std::uint16_t face = 0;
544 std::uint16_t edge = 0;
546 for (
unsigned int direction = 0; direction < dim; ++direction)
548 const auto side = (subcell >> direction) & 1U;
549 const auto face_no = direction * 2 + side;
552 if (cell->at_boundary(face_no))
555 const auto &neighbor = cell->neighbor(face_no);
559 if (neighbor->has_children() || neighbor->is_artificial() ||
560 neighbor->level() == cell->level())
564 if (neighbor->get_fe().n_dofs_per_cell() == 0)
567 face |= 1 << direction;
571 for (
unsigned int direction = 0; direction < dim; ++direction)
572 if (face == 0 || face == (1 << direction))
574 const unsigned int line_no =
581 const unsigned int line_index = cell->line_index(line_no);
583 const auto edge_neighbor =
584 std::find_if(line_to_cells[line_index].
begin(),
585 line_to_cells[line_index].
end(),
586 [&cell](
const auto &edge_neighbor) {
588 &cell->get_triangulation(),
591 &cell->get_dof_handler());
592 return dof_cell.is_artificial() ==
false &&
593 dof_cell.level() < cell->level() &&
594 dof_cell.
get_fe().n_dofs_per_cell() > 0;
597 if (edge_neighbor == line_to_cells[line_index].
end())
600 edge |= 1 << direction;
603 if ((face == 0) && (edge == 0))
606 const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7));
609 inverted_subcell + (face << 3) + (edge << 6));
611 return refinement_configuration;
617 template <
typename CellIterator>
620 const CellIterator &cell,
621 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
622 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
623 const std::vector<std::vector<bool>> &supported_components,
625 std::vector<types::global_dof_index> &dof_indices)
const
627 if (std::find(supported_components[cell->active_fe_index()].begin(),
628 supported_components[cell->active_fe_index()].end(),
630 supported_components[cell->active_fe_index()].end())
633 const auto &fe = cell->get_fe();
637 std::vector<std::vector<unsigned int>>
638 component_to_system_index_face_array(fe.n_components());
640 for (
unsigned int i = 0; i < fe.n_dofs_per_face(0); ++i)
641 component_to_system_index_face_array
642 [fe.face_system_to_component_index(i, 0).first]
645 std::vector<unsigned int> idx_offset = {0};
647 for (
unsigned int base_element_index = 0;
648 base_element_index < cell->get_fe().n_base_elements();
649 ++base_element_index)
650 for (
unsigned int c = 0;
651 c < cell->get_fe().element_multiplicity(base_element_index);
653 idx_offset.push_back(
655 cell->get_fe().base_element(base_element_index).n_dofs_per_cell());
657 std::vector<types::global_dof_index> neighbor_dofs_all(idx_offset.back());
658 std::vector<types::global_dof_index> neighbor_dofs_all_temp(
660 std::vector<types::global_dof_index> neighbor_dofs_face(
661 fe.n_dofs_per_face(0));
664 const auto get_face_idx = [](
const auto n_dofs_1d,
667 const auto j) ->
unsigned int {
668 const auto direction = face_no / 2;
669 const auto side = face_no % 2;
670 const auto offset = (side == 1) ? (n_dofs_1d - 1) : 0;
673 return (direction == 0) ? (n_dofs_1d * i + offset) :
674 (n_dofs_1d * offset + i);
679 return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset;
681 return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i;
683 return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j;
693 const std::uint16_t kind =
694 static_cast<std::uint16_t
>(refinement_configuration);
695 const std::uint16_t subcell = (kind >> 0) & 7;
696 const std::uint16_t
subcell_x = (subcell >> 0) & 1;
697 const std::uint16_t
subcell_y = (subcell >> 1) & 1;
698 const std::uint16_t
subcell_z = (subcell >> 2) & 1;
699 const std::uint16_t face = (kind >> 3) & 7;
700 const std::uint16_t edge = (kind >> 6) & 7;
702 for (
unsigned int direction = 0; direction < dim; ++direction)
703 if ((face >> direction) & 1U)
705 const auto side = ((subcell >> direction) & 1U) == 0;
706 const auto face_no = direction * 2 + side;
709 cell->neighbor(face_no)
710 ->face(cell->neighbor_face_no(face_no))
711 ->get_dof_indices(neighbor_dofs_face,
712 cell->neighbor(face_no)->active_fe_index());
716 for (
auto &
index : neighbor_dofs_face)
719 for (
unsigned int base_element_index = 0, comp = 0;
720 base_element_index < cell->get_fe().n_base_elements();
721 ++base_element_index)
722 for (
unsigned int c = 0;
723 c < cell->get_fe().element_multiplicity(base_element_index);
726 if (supported_components[cell->active_fe_index()][comp] ==
730 const unsigned int n_dofs_1d =
732 .base_element(base_element_index)
735 const unsigned int dofs_per_face =
737 std::vector<types::global_dof_index> neighbor_dofs(
739 const auto lex_face_mapping =
744 for (
unsigned int i = 0; i < dofs_per_face; ++i)
745 neighbor_dofs[i] = neighbor_dofs_face
746 [component_to_system_index_face_array[comp][i]];
753 Assert(cell->face_orientation(face_no),
759 if (cell->face_rotation(face_no))
761 if (cell->face_flip(face_no))
764 rotate_face(
rotate, n_dofs_1d, neighbor_dofs);
766 if (cell->face_orientation(face_no) ==
false)
767 transpose_face(n_dofs_1d - 1, neighbor_dofs);
775 for (
unsigned int i = 0, k = 0; i < n_dofs_1d; ++i)
776 for (
unsigned int j = 0; j < (dim == 2 ? 1 : n_dofs_1d);
778 dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) +
780 neighbor_dofs[lex_face_mapping[k]];
785 for (
unsigned int direction = 0; direction < dim; ++direction)
786 if ((edge >> direction) & 1U)
788 const unsigned int line_no =
794 const unsigned int line_index = cell->line(line_no)->index();
796 const auto edge_neighbor =
797 std::find_if(line_to_cells[line_index].
begin(),
798 line_to_cells[line_index].
end(),
799 [&cell](
const auto &edge_array) {
801 edge_neighbor(&cell->get_triangulation(),
804 return edge_neighbor->is_artificial() ==
false &&
805 edge_neighbor->level() < cell->level();
808 if (edge_neighbor == line_to_cells[line_index].
end())
812 &cell->get_triangulation(),
815 &cell->get_dof_handler());
816 const auto local_line_neighbor = (*edge_neighbor)[2];
821 for (
auto &
index : neighbor_dofs_all)
824 for (
unsigned int i = 0; i < neighbor_dofs_all_temp.size(); ++i)
825 neighbor_dofs_all_temp[i] = neighbor_dofs_all
826 [lexicographic_mapping[cell->active_fe_index()][i]];
829 cell->line_orientation(line_no) !=
830 neighbor_cell.line_orientation(local_line_neighbor);
832 for (
unsigned int base_element_index = 0, comp = 0;
833 base_element_index < cell->get_fe().n_base_elements();
834 ++base_element_index)
835 for (
unsigned int c = 0;
837 cell->get_fe().element_multiplicity(base_element_index);
840 if (supported_components[cell->active_fe_index()][comp] ==
844 const unsigned int n_dofs_1d =
846 .base_element(base_element_index)
850 for (
unsigned int i = 0; i < n_dofs_1d; ++i)
851 dof_indices[line_dof_idx(line_no, i, n_dofs_1d) +
852 idx_offset[comp]] = neighbor_dofs_all_temp
853 [line_dof_idx(local_line_neighbor,
854 flipped ? (n_dofs_1d - 1 - i) : i,
864 template <
typename CellIterator>
867 const CellIterator &cell,
868 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
869 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
870 std::vector<types::global_dof_index> &dof_indices,
874 const auto supported_components = compute_supported_components(
875 cell->get_dof_handler().get_fe_collection());
877 if (std::none_of(supported_components.begin(),
878 supported_components.end(),
880 return *std::max_element(a.begin(), a.end());
885 const auto refinement_configuration =
886 compute_refinement_configuration(cell);
892 update_dof_indices(cell,
894 lexicographic_mapping,
895 supported_components,
896 refinement_configuration,
900 for (
unsigned int c = 0; c < supported_components[0].size(); ++c)
901 if (supported_components[cell->active_fe_index()][c])
902 masks[c] = refinement_configuration;
912 unsigned int &subface_index)
const
914 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
917 times = times < 0 ? times + 4 : times;
918 for (
int t = 0; t < times; ++t)
919 subface_index = rot_mapping[subface_index];
928 unsigned int n_dofs_1d,
929 std::vector<types::global_dof_index> &dofs)
const
931 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
934 times = times < 0 ? times + 4 : times;
936 std::vector<types::global_dof_index>
copy(dofs.size());
937 for (
int t = 0; t < times; ++t)
942 for (
unsigned int i = 0; i < 4; ++i)
943 dofs[rot_mapping[i]] =
copy[i];
946 const unsigned int n_int = n_dofs_1d - 2;
947 unsigned int offset = 4;
948 for (
unsigned int i = 0; i < n_int; ++i)
951 dofs[offset + i] =
copy[offset + 2 * n_int + (n_int - 1 - i)];
953 dofs[offset + n_int + i] =
954 copy[offset + 3 * n_int + (n_int - 1 - i)];
956 dofs[offset + 2 * n_int + i] =
copy[offset + n_int + i];
958 dofs[offset + 3 * n_int + i] =
copy[offset + i];
964 for (
unsigned int i = 0; i < n_int; ++i)
965 for (
unsigned int j = 0; j < n_int; ++j)
966 dofs[offset + i * n_int + j] =
967 copy[offset + j * n_int + (n_int - 1 - i)];
977 unsigned int n_dofs_1d)
const
979 unsigned int x, y, z;
981 const unsigned int fe_degree = n_dofs_1d - 1;
985 x = (local_line % 4 == 0) ? 0 :
986 (local_line % 4 == 1) ? fe_degree :
988 y = (local_line % 4 == 2) ? 0 :
989 (local_line % 4 == 3) ? fe_degree :
991 z = (local_line / 4) * fe_degree;
995 x = ((local_line - 8) % 2) * fe_degree;
996 y = ((local_line - 8) / 2) * fe_degree;
1000 return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
1008 const unsigned int fe_degree,
1009 std::vector<types::global_dof_index> &dofs)
const
1011 const std::vector<types::global_dof_index>
copy(dofs);
1018 const unsigned int n_int = fe_degree - 1;
1019 unsigned int offset = 4;
1020 for (
unsigned int i = 0; i < n_int; ++i)
1023 dofs[offset + i] =
copy[offset + 2 * n_int + i];
1025 dofs[offset + n_int + i] =
copy[offset + 3 * n_int + i];
1027 dofs[offset + 2 * n_int + i] =
copy[offset + i];
1029 dofs[offset + 3 * n_int + i] =
copy[offset + n_int + i];
1033 offset += 4 * n_int;
1034 for (
unsigned int i = 0; i < n_int; ++i)
1035 for (
unsigned int j = 0; j < n_int; ++j)
1036 dofs[offset + i * n_int + j] =
copy[offset + j * n_int + i];
1047 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
unsigned int line_index(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
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)
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)
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