54 template <
int dim,
int spacedim,
typename number>
59 const bool keep_constrained_dofs,
80 "For distributed Triangulation objects and associated "
81 "DoFHandler objects, asking for any subdomain other than the "
82 "locally owned one does not make sense."));
86 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
87 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
89 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
92 std::vector<types::global_dof_index> dofs_on_this_cell;
100 (subdomain_id == cell->subdomain_id())) &&
101 cell->is_locally_owned())
103 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
104 dofs_on_this_cell.resize(dofs_per_cell);
105 cell->get_dof_indices(dofs_on_this_cell);
111 if (fe_dof_mask[fe_index].empty())
114 keep_constrained_dofs);
118 keep_constrained_dofs,
119 fe_dof_mask[fe_index]);
125 template <
int dim,
int spacedim,
typename number>
131 const bool keep_constrained_dofs,
158 "For distributed Triangulation objects and associated "
159 "DoFHandler objects, asking for any subdomain other than the "
160 "locally owned one does not make sense."));
166 const std::vector<Table<2, Coupling>> dof_mask
169 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.
size());
170 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
172 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
177 std::vector<Table<2, bool>> bool_dof_mask(fe_collection.
size());
178 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
180 bool_dof_mask[f].reinit(
182 fe_collection[f].n_dofs_per_cell()));
183 bool_dof_mask[f].fill(
false);
184 for (
unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
185 for (
unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
186 if (dof_mask[f](i, j) !=
none &&
187 (fe_dof_mask[f].empty() || fe_dof_mask[f](i, j)))
188 bool_dof_mask[f](i, j) =
true;
191 std::vector<types::global_dof_index> dofs_on_this_cell(
199 (subdomain_id == cell->subdomain_id())) &&
200 cell->is_locally_owned())
203 const unsigned int dofs_per_cell =
204 fe_collection[fe_index].n_dofs_per_cell();
206 dofs_on_this_cell.resize(dofs_per_cell);
207 cell->get_dof_indices(dofs_on_this_cell);
215 keep_constrained_dofs,
216 bool_dof_mask[fe_index]);
222 template <
int dim,
int spacedim>
249 ExcMessage(
"This function can only be used with with parallel "
250 "Triangulations when the Triangulations are equal."));
256 std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
259#ifdef DEAL_II_WITH_MPI
270 Assert(this_subdomain_id ==
277 [=](
const std::pair<cell_iterator, cell_iterator> &pair) {
278 return pair.first->subdomain_id() != this_subdomain_id ||
279 pair.second->subdomain_id() != this_subdomain_id;
285 for (
const auto &cell_pair : cell_list)
287 const cell_iterator cell_row = cell_pair.first;
288 const cell_iterator cell_col = cell_pair.second;
290 if (cell_row->is_active() && cell_col->is_active())
292 const unsigned int dofs_per_cell_row =
293 cell_row->get_fe().n_dofs_per_cell();
294 const unsigned int dofs_per_cell_col =
295 cell_col->get_fe().n_dofs_per_cell();
296 std::vector<types::global_dof_index> local_dof_indices_row(
298 std::vector<types::global_dof_index> local_dof_indices_col(
300 cell_row->get_dof_indices(local_dof_indices_row);
301 cell_col->get_dof_indices(local_dof_indices_col);
302 for (
const auto &dof : local_dof_indices_row)
306 else if (cell_row->has_children())
311 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
313 for (
unsigned int i = 0; i < child_cells.size(); ++i)
316 cell_row_child = child_cells[i];
317 const unsigned int dofs_per_cell_row =
319 const unsigned int dofs_per_cell_col =
320 cell_col->get_fe().n_dofs_per_cell();
321 std::vector<types::global_dof_index> local_dof_indices_row(
323 std::vector<types::global_dof_index> local_dof_indices_col(
325 cell_row_child->get_dof_indices(local_dof_indices_row);
326 cell_col->get_dof_indices(local_dof_indices_col);
327 for (
const auto &dof : local_dof_indices_row)
337 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
339 for (
unsigned int i = 0; i < child_cells.size(); ++i)
342 &cell_col_child = child_cells[i];
343 const unsigned int dofs_per_cell_row =
344 cell_row->get_fe().n_dofs_per_cell();
345 const unsigned int dofs_per_cell_col =
347 std::vector<types::global_dof_index> local_dof_indices_row(
349 std::vector<types::global_dof_index> local_dof_indices_col(
351 cell_row->get_dof_indices(local_dof_indices_row);
352 cell_col_child->get_dof_indices(local_dof_indices_col);
353 for (
const auto &dof : local_dof_indices_row)
363 template <
int dim,
int spacedim>
367 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
374 std::map<types::boundary_id, const Function<spacedim, double> *>
376 boundary_ids[0] =
nullptr;
377 boundary_ids[1] =
nullptr;
378 make_boundary_sparsity_pattern<dim, spacedim>(dof,
380 dof_to_boundary_mapping,
393 if (sparsity.
n_rows() != 0)
398 (index > max_element))
404 std::vector<types::global_dof_index> dofs_on_this_face;
406 std::vector<types::global_dof_index> cols;
414 for (
const unsigned int f : cell->face_indices())
415 if (cell->at_boundary(f))
417 const unsigned int dofs_per_face =
418 cell->get_fe().n_dofs_per_face(f);
419 dofs_on_this_face.resize(dofs_per_face);
420 cell->face(f)->get_dof_indices(dofs_on_this_face,
421 cell->active_fe_index());
425 for (
const auto &dof : dofs_on_this_face)
426 cols.push_back(dof_to_boundary_mapping[dof]);
430 std::sort(cols.begin(), cols.end());
431 for (
const auto &dof : dofs_on_this_face)
440 template <
int dim,
int spacedim,
typename number>
446 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
452 for (
unsigned int direction = 0; direction < 2; ++direction)
455 if (boundary_ids.find(direction) == boundary_ids.end())
462 while (!cell->at_boundary(direction))
463 cell = cell->neighbor(direction);
464 while (!cell->is_active())
465 cell = cell->child(direction);
467 const unsigned int dofs_per_vertex =
469 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
473 for (
unsigned int i = 0; i < dofs_per_vertex; ++i)
474 boundary_dof_boundary_indices[i] =
475 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
477 std::sort(boundary_dof_boundary_indices.begin(),
478 boundary_dof_boundary_indices.end());
479 for (
const auto &dof : boundary_dof_boundary_indices)
495 &(dof.
get_fe())) !=
nullptr);
509 if (sparsity.
n_rows() != 0)
514 (index > max_element))
520 std::vector<types::global_dof_index> dofs_on_this_face;
522 std::vector<types::global_dof_index> cols;
525 for (
const unsigned int f : cell->face_indices())
526 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
529 const unsigned int dofs_per_face =
530 cell->get_fe().n_dofs_per_face(f);
531 dofs_on_this_face.resize(dofs_per_face);
532 cell->face(f)->get_dof_indices(dofs_on_this_face,
533 cell->active_fe_index());
537 for (
const auto &dof : dofs_on_this_face)
538 cols.push_back(dof_to_boundary_mapping[dof]);
541 std::sort(cols.begin(), cols.end());
542 for (
const auto &dof : dofs_on_this_face)
551 template <
int dim,
int spacedim,
typename number>
556 const bool keep_constrained_dofs,
578 "For distributed Triangulation objects and associated "
579 "DoFHandler objects, asking for any subdomain other than the "
580 "locally owned one does not make sense."));
583 std::vector<types::global_dof_index> dofs_on_this_cell;
584 std::vector<types::global_dof_index> dofs_on_other_cell;
598 (subdomain_id == cell->subdomain_id())) &&
599 cell->is_locally_owned())
601 const unsigned int n_dofs_on_this_cell =
602 cell->get_fe().n_dofs_per_cell();
603 dofs_on_this_cell.resize(n_dofs_on_this_cell);
604 cell->get_dof_indices(dofs_on_this_cell);
611 keep_constrained_dofs);
613 for (
const unsigned int face : cell->face_indices())
617 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
618 if (!cell->at_boundary(face) || periodic_neighbor)
621 neighbor = cell->neighbor_or_periodic_neighbor(face);
628 while (neighbor->has_children())
629 neighbor = neighbor->child(face == 0 ? 1 : 0);
631 if (neighbor->has_children())
633 for (
unsigned int sub_nr = 0;
634 sub_nr != cell_face->n_active_descendants();
638 level_cell_iterator sub_neighbor =
640 cell->periodic_neighbor_child_on_subface(
642 cell->neighbor_child_on_subface(face, sub_nr);
644 const unsigned int n_dofs_on_neighbor =
646 dofs_on_other_cell.resize(n_dofs_on_neighbor);
647 sub_neighbor->get_dof_indices(dofs_on_other_cell);
653 keep_constrained_dofs);
658 keep_constrained_dofs);
662 if (sub_neighbor->subdomain_id() !=
663 cell->subdomain_id())
667 keep_constrained_dofs);
674 if ((!periodic_neighbor &&
675 cell->neighbor_is_coarser(face)) ||
676 (periodic_neighbor &&
677 cell->periodic_neighbor_is_coarser(face)))
678 if (neighbor->subdomain_id() == cell->subdomain_id())
681 const unsigned int n_dofs_on_neighbor =
683 dofs_on_other_cell.resize(n_dofs_on_neighbor);
685 neighbor->get_dof_indices(dofs_on_other_cell);
691 keep_constrained_dofs);
697 if (!cell->neighbor_or_periodic_neighbor(face)
699 (neighbor->subdomain_id() != cell->subdomain_id()))
705 keep_constrained_dofs);
706 if (neighbor->subdomain_id() != cell->subdomain_id())
710 keep_constrained_dofs);
720 template <
int dim,
int spacedim>
729 template <
int dim,
int spacedim>
746 for (
unsigned int i = 0; i < n_dofs; ++i)
748 const unsigned int ii =
754 for (
unsigned int j = 0; j < n_dofs; ++j)
756 const unsigned int jj =
762 dof_couplings(i, j) = component_couplings(ii, jj);
765 return dof_couplings;
770 template <
int dim,
int spacedim>
771 std::vector<Table<2, Coupling>>
776 std::vector<Table<2, Coupling>> return_value(fe.
size());
777 for (
unsigned int i = 0; i < fe.
size(); ++i)
791 template <
typename Iterator,
typename Iterator2>
794 const Iterator &cell,
795 const unsigned int face_no,
796 const Iterator2 &neighbor,
797 const unsigned int neighbor_face_no,
799 const std::vector<types::global_dof_index> &dofs_on_this_cell,
800 std::vector<types::global_dof_index> &dofs_on_other_cell,
804 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
805 neighbor->get_dof_indices(dofs_on_other_cell);
809 boost::container::small_vector<unsigned int, 64>
810 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
811 boost::container::small_vector<bool, 64> support_on_face_i(
812 neighbor->get_fe().n_dofs_per_cell());
813 boost::container::small_vector<bool, 64> support_on_face_e(
814 neighbor->get_fe().n_dofs_per_cell());
815 for (
unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
817 component_indices_neighbor[j] =
818 (neighbor->get_fe().is_primitive(j) ?
819 neighbor->get_fe().system_to_component_index(j).first :
821 .get_nonzero_components(j)
822 .first_selected_component());
823 support_on_face_i[j] =
824 neighbor->get_fe().has_support_on_face(j, face_no);
825 support_on_face_e[j] =
826 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
832 for (
int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
834 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
836 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
837 for (
unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
839 const unsigned int ii =
840 (fe.is_primitive(i) ?
841 fe.system_to_component_index(i).first :
842 fe.get_nonzero_components(i).first_selected_component());
845 const bool i_non_zero_i =
846 fe.has_support_on_face(i,
847 (f == 0 ? face_no : neighbor_face_no));
849 for (
unsigned int j = 0;
850 j < neighbor->get_fe().n_dofs_per_cell();
853 const bool j_non_zero_e = support_on_face_e[j];
854 const unsigned int jj = component_indices_neighbor[j];
856 Assert(jj < neighbor->get_fe().n_components(),
859 if ((flux_mask(ii, jj) ==
always) ||
860 (flux_mask(ii, jj) ==
nonzero && i_non_zero_i &&
862 cell_entries.emplace_back(dofs_i[i],
863 dofs_on_other_cell[j]);
864 if ((flux_mask(jj, ii) ==
always) ||
865 (flux_mask(jj, ii) ==
nonzero && j_non_zero_e &&
867 cell_entries.emplace_back(dofs_on_other_cell[j],
878 template <
int dim,
int spacedim,
typename number>
884 const bool keep_constrained_dofs,
890 const unsigned int)> &face_has_flux_coupling)
892 std::vector<types::global_dof_index> rows;
897 const ::hp::FECollection<dim, spacedim> &fe =
900 std::vector<types::global_dof_index> dofs_on_this_cell(
902 std::vector<types::global_dof_index> dofs_on_other_cell(
905 const unsigned int n_components = fe.n_components();
915 for (
unsigned int c1 = 0; c1 < n_components; ++c1)
916 for (
unsigned int c2 = 0; c2 < n_components; ++c2)
917 if (int_mask(c1, c2) !=
none || flux_mask(c1, c2) !=
none)
918 int_and_flux_mask(c1, c2) =
always;
922 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
927 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
928 for (
unsigned int f = 0; f < fe.size(); ++f)
930 bool_int_and_flux_dof_mask[f].reinit(
932 fe[f].n_dofs_per_cell()));
933 bool_int_and_flux_dof_mask[f].fill(
false);
934 for (
unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
935 for (
unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
936 if (int_and_flux_dof_mask[f](i, j) !=
none)
937 bool_int_and_flux_dof_mask[f](i, j) =
true;
941 for (
const auto &cell : dof.active_cell_iterators())
944 cell->is_locally_owned())
946 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
947 cell->get_dof_indices(dofs_on_this_cell);
955 keep_constrained_dofs,
956 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
959 for (
const unsigned int face : cell->face_indices())
961 const bool periodic_neighbor =
962 cell->has_periodic_neighbor(face);
964 if ((!cell->at_boundary(face)) || periodic_neighbor)
967 neighbor = cell->neighbor_or_periodic_neighbor(face);
973 if (neighbor->level() == cell->level() &&
974 neighbor->index() > cell->index() &&
975 neighbor->is_active() && neighbor->is_locally_owned())
987 if (neighbor->level() != cell->level() &&
988 ((!periodic_neighbor &&
989 !cell->neighbor_is_coarser(face)) ||
990 (periodic_neighbor &&
991 !cell->periodic_neighbor_is_coarser(face))) &&
992 neighbor->is_locally_owned())
995 if (!face_has_flux_coupling(cell, face))
998 const unsigned int neighbor_face_no =
1000 cell->periodic_neighbor_face_no(face) :
1001 cell->neighbor_face_no(face);
1010 while (neighbor->has_children())
1011 neighbor = neighbor->child(face == 0 ? 1 : 0);
1013 if (neighbor->has_children())
1015 for (
unsigned int sub_nr = 0;
1016 sub_nr != cell->face(face)->n_children();
1020 level_cell_iterator sub_neighbor =
1022 cell->periodic_neighbor_child_on_subface(
1024 cell->neighbor_child_on_subface(face,
1026 add_cell_entries(cell,
1037 add_cell_entries(cell,
1048 cell_entries.clear();
1057 template <
int dim,
int spacedim>
1067 const bool keep_constrained_dofs =
true;
1072 keep_constrained_dofs,
1076 internal::always_couple_on_faces<dim, spacedim>);
1081 template <
int dim,
int spacedim,
typename number>
1087 const bool keep_constrained_dofs,
1091 const std::function<
1093 const unsigned int)> &face_has_flux_coupling)
1106 Assert(int_mask.n_rows() == n_comp,
1108 Assert(int_mask.n_cols() == n_comp,
1110 Assert(flux_mask.n_rows() == n_comp,
1112 Assert(flux_mask.n_cols() == n_comp,
1125 "For distributed Triangulation objects and associated "
1126 "DoFHandler objects, asking for any subdomain other than the "
1127 "locally owned one does not make sense."));
1131 face_has_flux_coupling,
1133 "The function which specifies if a flux coupling occurs over a given "
1136 internal::make_flux_sparsity_pattern(dof,
1139 keep_constrained_dofs,
1143 face_has_flux_coupling);
1151#include "dofs/dof_tools_sparsity.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
typename LevelSelector::cell_iterator level_cell_iterator
cell_iterator begin(const unsigned int level=0) const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
types::global_dof_index size_type
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false)=0
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int size() const
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
void make_boundary_sparsity_pattern(const DoFHandler< dim, spacedim > &dof, const std::vector< types::global_dof_index > &dof_to_boundary_mapping, SparsityPatternBase &sparsity_pattern)
constexpr types::global_dof_index invalid_dof_index
constexpr types::boundary_id internal_face_boundary_id
constexpr types::subdomain_id invalid_subdomain_id
unsigned int subdomain_id
unsigned short int fe_index
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation