53 template <
int dim,
int spacedim,
typename number>
58 const bool keep_constrained_dofs,
72 if (
const auto *triangulation =
dynamic_cast<
77 (subdomain_id == triangulation->locally_owned_subdomain()),
79 "For distributed Triangulation objects and associated "
80 "DoFHandler objects, asking for any subdomain other than the "
81 "locally owned one does not make sense."));
85 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
86 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
88 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
91 std::vector<types::global_dof_index> dofs_on_this_cell;
99 (subdomain_id == cell->subdomain_id())) &&
100 cell->is_locally_owned())
102 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
103 dofs_on_this_cell.resize(dofs_per_cell);
104 cell->get_dof_indices(dofs_on_this_cell);
110 if (fe_dof_mask[fe_index].empty())
113 keep_constrained_dofs);
117 keep_constrained_dofs,
118 fe_dof_mask[fe_index]);
124 template <
int dim,
int spacedim,
typename number>
130 const bool keep_constrained_dofs,
150 if (
const auto *triangulation =
dynamic_cast<
155 (subdomain_id == triangulation->locally_owned_subdomain()),
157 "For distributed Triangulation objects and associated "
158 "DoFHandler objects, asking for any subdomain other than the "
159 "locally owned one does not make sense."));
165 const std::vector<Table<2, Coupling>> dof_mask
168 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.
size());
169 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
171 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
176 std::vector<Table<2, bool>> bool_dof_mask(fe_collection.
size());
177 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
179 bool_dof_mask[f].reinit(
181 fe_collection[f].n_dofs_per_cell()));
182 bool_dof_mask[f].fill(
false);
183 for (
unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
184 for (
unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
185 if (dof_mask[f](i, j) !=
none &&
186 (fe_dof_mask[f].empty() || fe_dof_mask[f](i, j)))
187 bool_dof_mask[f](i, j) =
true;
190 std::vector<types::global_dof_index> dofs_on_this_cell(
198 (subdomain_id == cell->subdomain_id())) &&
199 cell->is_locally_owned())
202 const unsigned int dofs_per_cell =
203 fe_collection[fe_index].n_dofs_per_cell();
205 dofs_on_this_cell.resize(dofs_per_cell);
206 cell->get_dof_indices(dofs_on_this_cell);
214 keep_constrained_dofs,
215 bool_dof_mask[fe_index]);
221 template <
int dim,
int spacedim>
248 ExcMessage(
"This function can only be used with with parallel "
249 "Triangulations when the Triangulations are equal."));
255 std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
258#ifdef DEAL_II_WITH_MPI
269 Assert(this_subdomain_id ==
276 [=](
const std::pair<cell_iterator, cell_iterator> &pair) {
277 return pair.first->subdomain_id() != this_subdomain_id ||
278 pair.second->subdomain_id() != this_subdomain_id;
284 for (
const auto &cell_pair : cell_list)
286 const cell_iterator cell_row = cell_pair.first;
287 const cell_iterator cell_col = cell_pair.second;
289 if (cell_row->is_active() && cell_col->is_active())
291 const unsigned int dofs_per_cell_row =
292 cell_row->get_fe().n_dofs_per_cell();
293 const unsigned int dofs_per_cell_col =
294 cell_col->get_fe().n_dofs_per_cell();
295 std::vector<types::global_dof_index> local_dof_indices_row(
297 std::vector<types::global_dof_index> local_dof_indices_col(
299 cell_row->get_dof_indices(local_dof_indices_row);
300 cell_col->get_dof_indices(local_dof_indices_col);
301 for (
const auto &dof : local_dof_indices_row)
305 else if (cell_row->has_children())
310 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
312 for (
unsigned int i = 0; i < child_cells.size(); ++i)
315 cell_row_child = child_cells[i];
316 const unsigned int dofs_per_cell_row =
318 const unsigned int dofs_per_cell_col =
319 cell_col->get_fe().n_dofs_per_cell();
320 std::vector<types::global_dof_index> local_dof_indices_row(
322 std::vector<types::global_dof_index> local_dof_indices_col(
324 cell_row_child->get_dof_indices(local_dof_indices_row);
325 cell_col->get_dof_indices(local_dof_indices_col);
326 for (
const auto &dof : local_dof_indices_row)
336 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
338 for (
unsigned int i = 0; i < child_cells.size(); ++i)
341 &cell_col_child = child_cells[i];
342 const unsigned int dofs_per_cell_row =
343 cell_row->get_fe().n_dofs_per_cell();
344 const unsigned int dofs_per_cell_col =
346 std::vector<types::global_dof_index> local_dof_indices_row(
348 std::vector<types::global_dof_index> local_dof_indices_col(
350 cell_row->get_dof_indices(local_dof_indices_row);
351 cell_col_child->get_dof_indices(local_dof_indices_col);
352 for (
const auto &dof : local_dof_indices_row)
362 template <
int dim,
int spacedim>
366 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
373 std::map<types::boundary_id, const Function<spacedim, double> *>
375 boundary_ids[0] =
nullptr;
376 boundary_ids[1] =
nullptr;
377 make_boundary_sparsity_pattern<dim, spacedim>(dof,
379 dof_to_boundary_mapping,
392 if (sparsity.
n_rows() != 0)
397 (index > max_element))
403 std::vector<types::global_dof_index> dofs_on_this_face;
405 std::vector<types::global_dof_index> cols;
413 for (
const unsigned int f : cell->face_indices())
414 if (cell->at_boundary(f))
416 const unsigned int dofs_per_face =
417 cell->get_fe().n_dofs_per_face(f);
418 dofs_on_this_face.resize(dofs_per_face);
419 cell->face(f)->get_dof_indices(dofs_on_this_face,
420 cell->active_fe_index());
424 for (
const auto &dof : dofs_on_this_face)
425 cols.push_back(dof_to_boundary_mapping[dof]);
429 std::sort(cols.begin(), cols.end());
430 for (
const auto &dof : dofs_on_this_face)
439 template <
int dim,
int spacedim,
typename number>
445 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
451 for (
unsigned int direction = 0; direction < 2; ++direction)
454 if (boundary_ids.find(direction) == boundary_ids.end())
461 while (!cell->at_boundary(direction))
462 cell = cell->neighbor(direction);
463 while (!cell->is_active())
464 cell = cell->child(direction);
466 const unsigned int dofs_per_vertex =
468 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
472 for (
unsigned int i = 0; i < dofs_per_vertex; ++i)
473 boundary_dof_boundary_indices[i] =
474 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
476 std::sort(boundary_dof_boundary_indices.begin(),
477 boundary_dof_boundary_indices.end());
478 for (
const auto &dof : boundary_dof_boundary_indices)
494 &(dof.
get_fe())) !=
nullptr);
508 if (sparsity.
n_rows() != 0)
513 (index > max_element))
519 std::vector<types::global_dof_index> dofs_on_this_face;
521 std::vector<types::global_dof_index> cols;
524 for (
const unsigned int f : cell->face_indices())
525 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
528 const unsigned int dofs_per_face =
529 cell->get_fe().n_dofs_per_face(f);
530 dofs_on_this_face.resize(dofs_per_face);
531 cell->face(f)->get_dof_indices(dofs_on_this_face,
532 cell->active_fe_index());
536 for (
const auto &dof : dofs_on_this_face)
537 cols.push_back(dof_to_boundary_mapping[dof]);
540 std::sort(cols.begin(), cols.end());
541 for (
const auto &dof : dofs_on_this_face)
550 template <
int dim,
int spacedim,
typename number>
555 const bool keep_constrained_dofs,
570 if (
const auto *triangulation =
dynamic_cast<
575 (subdomain_id == triangulation->locally_owned_subdomain()),
577 "For distributed Triangulation objects and associated "
578 "DoFHandler objects, asking for any subdomain other than the "
579 "locally owned one does not make sense."));
582 std::vector<types::global_dof_index> dofs_on_this_cell;
583 std::vector<types::global_dof_index> dofs_on_other_cell;
597 (subdomain_id == cell->subdomain_id())) &&
598 cell->is_locally_owned())
600 const unsigned int n_dofs_on_this_cell =
601 cell->get_fe().n_dofs_per_cell();
602 dofs_on_this_cell.resize(n_dofs_on_this_cell);
603 cell->get_dof_indices(dofs_on_this_cell);
610 keep_constrained_dofs);
612 for (
const unsigned int face : cell->face_indices())
616 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
617 if (!cell->at_boundary(face) || periodic_neighbor)
620 neighbor = cell->neighbor_or_periodic_neighbor(face);
627 while (neighbor->has_children())
628 neighbor = neighbor->child(face == 0 ? 1 : 0);
630 if (neighbor->has_children())
632 for (
unsigned int sub_nr = 0;
633 sub_nr != cell_face->n_active_descendants();
637 level_cell_iterator sub_neighbor =
639 cell->periodic_neighbor_child_on_subface(
641 cell->neighbor_child_on_subface(face, sub_nr);
643 const unsigned int n_dofs_on_neighbor =
645 dofs_on_other_cell.resize(n_dofs_on_neighbor);
646 sub_neighbor->get_dof_indices(dofs_on_other_cell);
652 keep_constrained_dofs);
657 keep_constrained_dofs);
661 if (sub_neighbor->subdomain_id() !=
662 cell->subdomain_id())
666 keep_constrained_dofs);
673 if ((!periodic_neighbor &&
674 cell->neighbor_is_coarser(face)) ||
675 (periodic_neighbor &&
676 cell->periodic_neighbor_is_coarser(face)))
677 if (neighbor->subdomain_id() == cell->subdomain_id())
680 const unsigned int n_dofs_on_neighbor =
682 dofs_on_other_cell.resize(n_dofs_on_neighbor);
684 neighbor->get_dof_indices(dofs_on_other_cell);
690 keep_constrained_dofs);
696 if (!cell->neighbor_or_periodic_neighbor(face)
698 (neighbor->subdomain_id() != cell->subdomain_id()))
704 keep_constrained_dofs);
705 if (neighbor->subdomain_id() != cell->subdomain_id())
709 keep_constrained_dofs);
719 template <
int dim,
int spacedim>
728 template <
int dim,
int spacedim>
745 for (
unsigned int i = 0; i < n_dofs; ++i)
747 const unsigned int ii =
753 for (
unsigned int j = 0; j < n_dofs; ++j)
755 const unsigned int jj =
761 dof_couplings(i, j) = component_couplings(ii, jj);
764 return dof_couplings;
769 template <
int dim,
int spacedim>
770 std::vector<Table<2, Coupling>>
775 std::vector<Table<2, Coupling>> return_value(fe.
size());
776 for (
unsigned int i = 0; i < fe.
size(); ++i)
790 template <
typename Iterator,
typename Iterator2>
793 const Iterator &cell,
794 const unsigned int face_no,
795 const Iterator2 &neighbor,
796 const unsigned int neighbor_face_no,
798 const std::vector<types::global_dof_index> &dofs_on_this_cell,
799 std::vector<types::global_dof_index> &dofs_on_other_cell,
803 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
804 neighbor->get_dof_indices(dofs_on_other_cell);
808 boost::container::small_vector<unsigned int, 64>
809 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
810 boost::container::small_vector<bool, 64> support_on_face_i(
811 neighbor->get_fe().n_dofs_per_cell());
812 boost::container::small_vector<bool, 64> support_on_face_e(
813 neighbor->get_fe().n_dofs_per_cell());
814 for (
unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
816 component_indices_neighbor[j] =
817 (neighbor->get_fe().is_primitive(j) ?
818 neighbor->get_fe().system_to_component_index(j).first :
820 .get_nonzero_components(j)
821 .first_selected_component());
822 support_on_face_i[j] =
823 neighbor->get_fe().has_support_on_face(j, face_no);
824 support_on_face_e[j] =
825 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
831 for (
int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
833 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
835 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
836 for (
unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
838 const unsigned int ii =
839 (fe.is_primitive(i) ?
840 fe.system_to_component_index(i).first :
841 fe.get_nonzero_components(i).first_selected_component());
844 const bool i_non_zero_i =
845 fe.has_support_on_face(i,
846 (f == 0 ? face_no : neighbor_face_no));
848 for (
unsigned int j = 0;
849 j < neighbor->get_fe().n_dofs_per_cell();
852 const bool j_non_zero_e = support_on_face_e[j];
853 const unsigned int jj = component_indices_neighbor[j];
855 Assert(jj < neighbor->get_fe().n_components(),
858 if ((flux_mask(ii, jj) ==
always) ||
859 (flux_mask(ii, jj) ==
nonzero && i_non_zero_i &&
861 cell_entries.emplace_back(dofs_i[i],
862 dofs_on_other_cell[j]);
863 if ((flux_mask(jj, ii) ==
always) ||
864 (flux_mask(jj, ii) ==
nonzero && j_non_zero_e &&
866 cell_entries.emplace_back(dofs_on_other_cell[j],
877 template <
int dim,
int spacedim,
typename number>
883 const bool keep_constrained_dofs,
889 const unsigned int)> &face_has_flux_coupling)
891 std::vector<types::global_dof_index> rows;
896 const ::hp::FECollection<dim, spacedim> &fe =
899 std::vector<types::global_dof_index> dofs_on_this_cell(
901 std::vector<types::global_dof_index> dofs_on_other_cell(
904 const unsigned int n_components = fe.n_components();
914 for (
unsigned int c1 = 0; c1 < n_components; ++c1)
915 for (
unsigned int c2 = 0; c2 < n_components; ++c2)
916 if (int_mask(c1, c2) !=
none || flux_mask(c1, c2) !=
none)
917 int_and_flux_mask(c1, c2) =
always;
921 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
926 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
927 for (
unsigned int f = 0; f < fe.size(); ++f)
929 bool_int_and_flux_dof_mask[f].reinit(
931 fe[f].n_dofs_per_cell()));
932 bool_int_and_flux_dof_mask[f].fill(
false);
933 for (
unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
934 for (
unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
935 if (int_and_flux_dof_mask[f](i, j) !=
none)
936 bool_int_and_flux_dof_mask[f](i, j) =
true;
940 for (
const auto &cell : dof.active_cell_iterators())
943 cell->is_locally_owned())
945 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
946 cell->get_dof_indices(dofs_on_this_cell);
954 keep_constrained_dofs,
955 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
958 for (
const unsigned int face : cell->face_indices())
960 const bool periodic_neighbor =
961 cell->has_periodic_neighbor(face);
963 if ((!cell->at_boundary(face)) || periodic_neighbor)
966 neighbor = cell->neighbor_or_periodic_neighbor(face);
972 if (neighbor->level() == cell->level() &&
973 neighbor->index() > cell->index() &&
974 neighbor->is_active() && neighbor->is_locally_owned())
986 if (neighbor->level() != cell->level() &&
987 ((!periodic_neighbor &&
988 !cell->neighbor_is_coarser(face)) ||
989 (periodic_neighbor &&
990 !cell->periodic_neighbor_is_coarser(face))) &&
991 neighbor->is_locally_owned())
994 if (!face_has_flux_coupling(cell, face))
997 const unsigned int neighbor_face_no =
999 cell->periodic_neighbor_face_no(face) :
1000 cell->neighbor_face_no(face);
1009 while (neighbor->has_children())
1010 neighbor = neighbor->child(face == 0 ? 1 : 0);
1012 if (neighbor->has_children())
1014 for (
unsigned int sub_nr = 0;
1015 sub_nr != cell->face(face)->n_children();
1019 level_cell_iterator sub_neighbor =
1021 cell->periodic_neighbor_child_on_subface(
1023 cell->neighbor_child_on_subface(face,
1025 add_cell_entries(cell,
1036 add_cell_entries(cell,
1047 cell_entries.clear();
1056 template <
int dim,
int spacedim>
1066 const bool keep_constrained_dofs =
true;
1071 keep_constrained_dofs,
1075 internal::always_couple_on_faces<dim, spacedim>);
1080 template <
int dim,
int spacedim,
typename number>
1086 const bool keep_constrained_dofs,
1090 const std::function<
1092 const unsigned int)> &face_has_flux_coupling)
1105 Assert(int_mask.n_rows() == n_comp,
1107 Assert(int_mask.n_cols() == n_comp,
1109 Assert(flux_mask.n_rows() == n_comp,
1111 Assert(flux_mask.n_cols() == n_comp,
1117 if (
const auto *triangulation =
dynamic_cast<
1122 (subdomain_id == triangulation->locally_owned_subdomain()),
1124 "For distributed Triangulation objects and associated "
1125 "DoFHandler objects, asking for any subdomain other than the "
1126 "locally owned one does not make sense."));
1130 face_has_flux_coupling,
1132 "The function which specifies if a flux coupling occurs over a given "
1135 internal::make_flux_sparsity_pattern(dof,
1138 keep_constrained_dofs,
1142 face_has_flux_coupling);
1150#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)
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
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