16 #include <deal.II/base/quadrature_lib.h> 17 #include <deal.II/base/table.h> 18 #include <deal.II/base/template_constraints.h> 19 #include <deal.II/base/thread_management.h> 20 #include <deal.II/base/utilities.h> 22 #include <deal.II/dofs/dof_accessor.h> 23 #include <deal.II/dofs/dof_handler.h> 24 #include <deal.II/dofs/dof_tools.h> 26 #include <deal.II/fe/fe.h> 27 #include <deal.II/fe/fe_tools.h> 28 #include <deal.II/fe/fe_values.h> 30 #include <deal.II/grid/grid_tools.h> 31 #include <deal.II/grid/intergrid_map.h> 32 #include <deal.II/grid/tria.h> 33 #include <deal.II/grid/tria_iterator.h> 35 #include <deal.II/hp/fe_collection.h> 36 #include <deal.II/hp/fe_values.h> 37 #include <deal.II/hp/q_collection.h> 39 #include <deal.II/lac/affine_constraints.h> 40 #include <deal.II/lac/block_sparsity_pattern.h> 41 #include <deal.II/lac/dynamic_sparsity_pattern.h> 42 #include <deal.II/lac/sparsity_pattern.h> 43 #include <deal.II/lac/trilinos_sparsity_pattern.h> 44 #include <deal.II/lac/vector.h> 46 #include <deal.II/numerics/vector_tools.h> 51 DEAL_II_NAMESPACE_OPEN
57 template <
typename DoFHandlerType,
58 typename SparsityPatternType,
62 SparsityPatternType & sparsity,
64 const bool keep_constrained_dofs,
70 Assert(sparsity.n_rows() == n_dofs,
72 Assert(sparsity.n_cols() == n_dofs,
78 Assert((dof.get_triangulation().locally_owned_subdomain() ==
82 dof.get_triangulation().locally_owned_subdomain()),
84 "For parallel::distributed::Triangulation objects and " 85 "associated DoF handler objects, asking for any subdomain other " 86 "than the locally owned one does not make sense."));
88 std::vector<types::global_dof_index> dofs_on_this_cell;
90 typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
96 for (; cell != endc; ++cell)
98 (subdomain_id == cell->subdomain_id())) &&
99 cell->is_locally_owned())
101 const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
102 dofs_on_this_cell.resize(dofs_per_cell);
103 cell->get_dof_indices(dofs_on_this_cell);
110 keep_constrained_dofs);
116 template <
typename DoFHandlerType,
117 typename SparsityPatternType,
122 SparsityPatternType & sparsity,
124 const bool keep_constrained_dofs,
130 Assert(sparsity.n_rows() == n_dofs,
132 Assert(sparsity.n_cols() == n_dofs,
134 Assert(couplings.n_rows() == dof.get_fe(0).n_components(),
136 dof.get_fe(0).n_components()));
137 Assert(couplings.n_cols() == dof.get_fe(0).n_components(),
139 dof.get_fe(0).n_components()));
144 Assert((dof.get_triangulation().locally_owned_subdomain() ==
148 dof.get_triangulation().locally_owned_subdomain()),
150 "For parallel::distributed::Triangulation objects and " 151 "associated DoF handler objects, asking for any subdomain other " 152 "than the locally owned one does not make sense."));
155 DoFHandlerType::space_dimension> &fe_collection =
156 dof.get_fe_collection();
158 const std::vector<Table<2, Coupling>> dof_mask
163 std::vector<Table<2, bool>> bool_dof_mask(fe_collection.size());
164 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
166 bool_dof_mask[f].reinit(
168 fe_collection[f].dofs_per_cell));
169 bool_dof_mask[f].fill(
false);
170 for (
unsigned int i = 0; i < fe_collection[f].dofs_per_cell; ++i)
171 for (
unsigned int j = 0; j < fe_collection[f].dofs_per_cell; ++j)
172 if (dof_mask[f](i, j) !=
none)
173 bool_dof_mask[f](i, j) =
true;
176 std::vector<types::global_dof_index> dofs_on_this_cell(
177 fe_collection.max_dofs_per_cell());
178 typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
184 for (; cell != endc; ++cell)
186 (subdomain_id == cell->subdomain_id())) &&
187 cell->is_locally_owned())
189 const unsigned int fe_index = cell->active_fe_index();
190 const unsigned int dofs_per_cell =
191 fe_collection[fe_index].dofs_per_cell;
193 dofs_on_this_cell.resize(dofs_per_cell);
194 cell->get_dof_indices(dofs_on_this_cell);
202 keep_constrained_dofs,
203 bool_dof_mask[fe_index]);
209 template <
typename DoFHandlerType,
typename SparsityPatternType>
212 const DoFHandlerType &dof_col,
213 SparsityPatternType & sparsity)
220 Assert(sparsity.n_rows() == n_dofs_row,
222 Assert(sparsity.n_cols() == n_dofs_col,
227 const std::list<std::pair<
typename DoFHandlerType::cell_iterator,
228 typename DoFHandlerType::cell_iterator>>
232 typename std::list<std::pair<
typename DoFHandlerType::cell_iterator,
233 typename DoFHandlerType::cell_iterator>>::
234 const_iterator cell_iter = cell_list.begin();
236 for (; cell_iter != cell_list.end(); ++cell_iter)
238 const typename DoFHandlerType::cell_iterator cell_row =
240 const typename DoFHandlerType::cell_iterator cell_col =
243 if (!cell_row->has_children() && !cell_col->has_children())
245 const unsigned int dofs_per_cell_row =
246 cell_row->get_fe().dofs_per_cell;
247 const unsigned int dofs_per_cell_col =
248 cell_col->get_fe().dofs_per_cell;
249 std::vector<types::global_dof_index> local_dof_indices_row(
251 std::vector<types::global_dof_index> local_dof_indices_col(
253 cell_row->get_dof_indices(local_dof_indices_row);
254 cell_col->get_dof_indices(local_dof_indices_col);
255 for (
unsigned int i = 0; i < dofs_per_cell_row; ++i)
256 sparsity.add_entries(local_dof_indices_row[i],
257 local_dof_indices_col.begin(),
258 local_dof_indices_col.end());
260 else if (cell_row->has_children())
262 const std::vector<typename DoFHandlerType::active_cell_iterator>
264 GridTools::get_active_child_cells<DoFHandlerType>(cell_row);
265 for (
unsigned int i = 0; i < child_cells.size(); i++)
267 const typename DoFHandlerType::cell_iterator cell_row_child =
269 const unsigned int dofs_per_cell_row =
270 cell_row_child->get_fe().dofs_per_cell;
271 const unsigned int dofs_per_cell_col =
272 cell_col->get_fe().dofs_per_cell;
273 std::vector<types::global_dof_index> local_dof_indices_row(
275 std::vector<types::global_dof_index> local_dof_indices_col(
277 cell_row_child->get_dof_indices(local_dof_indices_row);
278 cell_col->get_dof_indices(local_dof_indices_col);
279 for (
unsigned int r = 0; r < dofs_per_cell_row; ++r)
280 sparsity.add_entries(local_dof_indices_row[r],
281 local_dof_indices_col.begin(),
282 local_dof_indices_col.end());
287 std::vector<typename DoFHandlerType::active_cell_iterator>
289 GridTools::get_active_child_cells<DoFHandlerType>(cell_col);
290 for (
unsigned int i = 0; i < child_cells.size(); i++)
292 const typename DoFHandlerType::active_cell_iterator
293 cell_col_child = child_cells[i];
294 const unsigned int dofs_per_cell_row =
295 cell_row->get_fe().dofs_per_cell;
296 const unsigned int dofs_per_cell_col =
297 cell_col_child->get_fe().dofs_per_cell;
298 std::vector<types::global_dof_index> local_dof_indices_row(
300 std::vector<types::global_dof_index> local_dof_indices_col(
302 cell_row->get_dof_indices(local_dof_indices_row);
303 cell_col_child->get_dof_indices(local_dof_indices_col);
304 for (
unsigned int r = 0; r < dofs_per_cell_row; ++r)
305 sparsity.add_entries(local_dof_indices_row[r],
306 local_dof_indices_col.begin(),
307 local_dof_indices_col.end());
315 template <
typename DoFHandlerType,
typename SparsityPatternType>
318 const DoFHandlerType & dof,
319 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
320 SparsityPatternType & sparsity)
322 if (DoFHandlerType::dimension == 1)
329 boundary_ids[0] =
nullptr;
330 boundary_ids[1] =
nullptr;
331 make_boundary_sparsity_pattern<DoFHandlerType, SparsityPatternType>(
332 dof, boundary_ids, dof_to_boundary_mapping, sparsity);
343 if (sparsity.n_rows() != 0)
346 for (std::vector<types::global_dof_index>::const_iterator i =
347 dof_to_boundary_mapping.begin();
348 i != dof_to_boundary_mapping.end();
356 std::vector<types::global_dof_index> dofs_on_this_face;
364 typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
366 for (; cell != endc; ++cell)
367 for (
unsigned int f = 0;
368 f < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
370 if (cell->at_boundary(f))
372 const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
373 dofs_on_this_face.resize(dofs_per_face);
374 cell->face(f)->get_dof_indices(dofs_on_this_face,
375 cell->active_fe_index());
378 for (
unsigned int i = 0; i < dofs_per_face; ++i)
379 for (
unsigned int j = 0; j < dofs_per_face; ++j)
380 sparsity.add(dof_to_boundary_mapping[dofs_on_this_face[i]],
381 dof_to_boundary_mapping[dofs_on_this_face[j]]);
387 template <
typename DoFHandlerType,
388 typename SparsityPatternType,
392 const DoFHandlerType &dof,
396 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
397 SparsityPatternType & sparsity)
399 if (DoFHandlerType::dimension == 1)
402 for (
unsigned int direction = 0; direction < 2; ++direction)
405 if (boundary_ids.find(direction) == boundary_ids.end())
410 typename DoFHandlerType::level_cell_iterator cell = dof.begin(0);
411 while (!cell->at_boundary(direction))
412 cell = cell->neighbor(direction);
413 while (!cell->active())
414 cell = cell->child(direction);
416 const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
417 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
421 for (
unsigned int i = 0; i < dofs_per_vertex; ++i)
422 boundary_dof_boundary_indices[i] =
423 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
425 for (
unsigned int i = 0; i < dofs_per_vertex; ++i)
426 sparsity.add_entries(boundary_dof_boundary_indices[i],
427 boundary_dof_boundary_indices.begin(),
428 boundary_dof_boundary_indices.end());
440 Assert(sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
442 dof.n_boundary_dofs(boundary_ids)));
443 Assert(sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
445 dof.n_boundary_dofs(boundary_ids)));
447 if (sparsity.n_rows() != 0)
450 for (std::vector<types::global_dof_index>::const_iterator i =
451 dof_to_boundary_mapping.begin();
452 i != dof_to_boundary_mapping.end();
460 std::vector<types::global_dof_index> dofs_on_this_face;
462 typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
464 for (; cell != endc; ++cell)
465 for (
unsigned int f = 0;
466 f < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
468 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
471 const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
472 dofs_on_this_face.resize(dofs_per_face);
473 cell->face(f)->get_dof_indices(dofs_on_this_face,
474 cell->active_fe_index());
477 for (
unsigned int i = 0; i < dofs_per_face; ++i)
478 for (
unsigned int j = 0; j < dofs_per_face; ++j)
479 sparsity.add(dof_to_boundary_mapping[dofs_on_this_face[i]],
480 dof_to_boundary_mapping[dofs_on_this_face[j]]);
486 template <
typename DoFHandlerType,
487 typename SparsityPatternType,
491 SparsityPatternType & sparsity,
493 const bool keep_constrained_dofs,
508 Assert((dof.get_triangulation().locally_owned_subdomain() ==
512 dof.get_triangulation().locally_owned_subdomain()),
514 "For parallel::distributed::Triangulation objects and " 515 "associated DoF handler objects, asking for any subdomain other " 516 "than the locally owned one does not make sense."));
518 std::vector<types::global_dof_index> dofs_on_this_cell;
519 std::vector<types::global_dof_index> dofs_on_other_cell;
522 typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
533 for (; cell != endc; ++cell)
535 (subdomain_id == cell->subdomain_id())) &&
536 cell->is_locally_owned())
538 const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell;
539 dofs_on_this_cell.resize(n_dofs_on_this_cell);
540 cell->get_dof_indices(dofs_on_this_cell);
547 keep_constrained_dofs);
549 for (
unsigned int face = 0;
550 face < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
553 typename DoFHandlerType::face_iterator cell_face =
555 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
556 if (!cell->at_boundary(face) || periodic_neighbor)
558 typename DoFHandlerType::level_cell_iterator neighbor =
559 cell->neighbor_or_periodic_neighbor(face);
565 if (DoFHandlerType::dimension == 1)
566 while (neighbor->has_children())
567 neighbor = neighbor->child(face == 0 ? 1 : 0);
569 if (neighbor->has_children())
571 for (
unsigned int sub_nr = 0;
572 sub_nr != cell_face->number_of_children();
575 const typename DoFHandlerType::level_cell_iterator
578 cell->periodic_neighbor_child_on_subface(
580 cell->neighbor_child_on_subface(face, sub_nr);
582 const unsigned int n_dofs_on_neighbor =
583 sub_neighbor->get_fe().dofs_per_cell;
584 dofs_on_other_cell.resize(n_dofs_on_neighbor);
585 sub_neighbor->get_dof_indices(dofs_on_other_cell);
591 keep_constrained_dofs);
596 keep_constrained_dofs);
600 if (sub_neighbor->subdomain_id() !=
601 cell->subdomain_id())
605 keep_constrained_dofs);
612 if ((!periodic_neighbor &&
613 cell->neighbor_is_coarser(face)) ||
614 (periodic_neighbor &&
615 cell->periodic_neighbor_is_coarser(face)))
616 if (neighbor->subdomain_id() == cell->subdomain_id())
619 const unsigned int n_dofs_on_neighbor =
620 neighbor->get_fe().dofs_per_cell;
621 dofs_on_other_cell.resize(n_dofs_on_neighbor);
623 neighbor->get_dof_indices(dofs_on_other_cell);
629 keep_constrained_dofs);
635 if (!cell->neighbor_or_periodic_neighbor(face)
637 (neighbor->subdomain_id() != cell->subdomain_id()))
643 keep_constrained_dofs);
644 if (neighbor->subdomain_id() != cell->subdomain_id())
648 keep_constrained_dofs);
658 template <
typename DoFHandlerType,
typename SparsityPatternType>
661 SparsityPatternType & sparsity)
667 template <
int dim,
int spacedim>
684 for (
unsigned int i = 0; i < n_dofs; ++i)
686 const unsigned int ii =
692 for (
unsigned int j = 0; j < n_dofs; ++j)
694 const unsigned int jj =
700 dof_couplings(i, j) = component_couplings(ii, jj);
703 return dof_couplings;
708 template <
int dim,
int spacedim>
709 std::vector<Table<2, Coupling>>
714 std::vector<Table<2, Coupling>> return_value(fe.
size());
715 for (
unsigned int i = 0; i < fe.
size(); ++i)
730 template <
typename DoFHandlerType,
731 typename SparsityPatternType,
735 SparsityPatternType & sparsity,
737 const bool keep_constrained_dofs,
743 DoFHandlerType::space_dimension> &fe = dof.get_fe();
745 std::vector<types::global_dof_index> dofs_on_this_cell(
747 std::vector<types::global_dof_index> dofs_on_other_cell(
757 for (
unsigned int i = 0; i < fe.dofs_per_cell; ++i)
758 for (
unsigned int f = 0;
759 f < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
761 support_on_face(i, f) = fe.has_support_on_face(i, f);
765 Table<2, bool> bool_int_dof_mask(fe.dofs_per_cell, fe.dofs_per_cell);
766 bool_int_dof_mask.
fill(
false);
767 for (
unsigned int i = 0; i < fe.dofs_per_cell; ++i)
768 for (
unsigned int j = 0; j < fe.dofs_per_cell; ++j)
769 if (int_dof_mask(i, j) !=
none)
770 bool_int_dof_mask(i, j) =
true;
772 typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
774 for (; cell != endc; ++cell)
776 (subdomain_id == cell->subdomain_id())) &&
777 cell->is_locally_owned())
779 cell->get_dof_indices(dofs_on_this_cell);
783 keep_constrained_dofs,
786 for (
unsigned int face_n = 0;
791 const typename DoFHandlerType::face_iterator cell_face =
794 const bool periodic_neighbor =
795 cell->has_periodic_neighbor(face_n);
797 if (cell->at_boundary(face_n) && (!periodic_neighbor))
799 for (
unsigned int i = 0; i < fe.dofs_per_cell; ++i)
801 const bool i_non_zero_i = support_on_face(i, face_n);
802 for (
unsigned int j = 0; j < fe.dofs_per_cell; ++j)
804 const bool j_non_zero_i =
805 support_on_face(j, face_n);
807 if (flux_dof_mask(i, j) ==
always ||
808 (flux_dof_mask(i, j) ==
nonzero &&
809 i_non_zero_i && j_non_zero_i))
810 sparsity.add(dofs_on_this_cell[i],
811 dofs_on_this_cell[j]);
817 typename DoFHandlerType::level_cell_iterator neighbor =
818 cell->neighbor_or_periodic_neighbor(face_n);
823 if (neighbor->level() == cell->level() &&
824 neighbor->index() > cell->index() &&
825 neighbor->active() && neighbor->is_locally_owned())
838 if (neighbor->level() != cell->level() &&
839 ((!periodic_neighbor &&
840 !cell->neighbor_is_coarser(face_n)) ||
841 (periodic_neighbor &&
842 !cell->periodic_neighbor_is_coarser(face_n))) &&
843 neighbor->is_locally_owned())
846 const unsigned int neighbor_face_n =
848 cell->periodic_neighbor_face_no(face_n) :
849 cell->neighbor_face_no(face_n);
851 if (neighbor->has_children())
853 for (
unsigned int sub_nr = 0;
854 sub_nr != cell_face->n_children();
857 const typename DoFHandlerType::level_cell_iterator
860 cell->periodic_neighbor_child_on_subface(
862 cell->neighbor_child_on_subface(face_n,
865 sub_neighbor->get_dof_indices(dofs_on_other_cell);
866 for (
unsigned int i = 0; i < fe.dofs_per_cell;
869 const bool i_non_zero_i =
870 support_on_face(i, face_n);
871 const bool i_non_zero_e =
872 support_on_face(i, neighbor_face_n);
873 for (
unsigned int j = 0; j < fe.dofs_per_cell;
876 const bool j_non_zero_i =
877 support_on_face(j, face_n);
878 const bool j_non_zero_e =
879 support_on_face(j, neighbor_face_n);
881 if (flux_dof_mask(i, j) ==
always)
883 sparsity.add(dofs_on_this_cell[i],
884 dofs_on_other_cell[j]);
885 sparsity.add(dofs_on_other_cell[i],
886 dofs_on_this_cell[j]);
887 sparsity.add(dofs_on_this_cell[i],
888 dofs_on_this_cell[j]);
889 sparsity.add(dofs_on_other_cell[i],
890 dofs_on_other_cell[j]);
892 else if (flux_dof_mask(i, j) ==
nonzero)
894 if (i_non_zero_i && j_non_zero_e)
895 sparsity.add(dofs_on_this_cell[i],
896 dofs_on_other_cell[j]);
897 if (i_non_zero_e && j_non_zero_i)
898 sparsity.add(dofs_on_other_cell[i],
899 dofs_on_this_cell[j]);
900 if (i_non_zero_i && j_non_zero_i)
901 sparsity.add(dofs_on_this_cell[i],
902 dofs_on_this_cell[j]);
903 if (i_non_zero_e && j_non_zero_e)
904 sparsity.add(dofs_on_other_cell[i],
905 dofs_on_other_cell[j]);
908 if (flux_dof_mask(j, i) ==
always)
910 sparsity.add(dofs_on_this_cell[j],
911 dofs_on_other_cell[i]);
912 sparsity.add(dofs_on_other_cell[j],
913 dofs_on_this_cell[i]);
914 sparsity.add(dofs_on_this_cell[j],
915 dofs_on_this_cell[i]);
916 sparsity.add(dofs_on_other_cell[j],
917 dofs_on_other_cell[i]);
919 else if (flux_dof_mask(j, i) ==
nonzero)
921 if (j_non_zero_i && i_non_zero_e)
922 sparsity.add(dofs_on_this_cell[j],
923 dofs_on_other_cell[i]);
924 if (j_non_zero_e && i_non_zero_i)
925 sparsity.add(dofs_on_other_cell[j],
926 dofs_on_this_cell[i]);
927 if (j_non_zero_i && i_non_zero_i)
928 sparsity.add(dofs_on_this_cell[j],
929 dofs_on_this_cell[i]);
930 if (j_non_zero_e && i_non_zero_e)
931 sparsity.add(dofs_on_other_cell[j],
932 dofs_on_other_cell[i]);
940 neighbor->get_dof_indices(dofs_on_other_cell);
941 for (
unsigned int i = 0; i < fe.dofs_per_cell; ++i)
943 const bool i_non_zero_i =
944 support_on_face(i, face_n);
945 const bool i_non_zero_e =
946 support_on_face(i, neighbor_face_n);
947 for (
unsigned int j = 0; j < fe.dofs_per_cell;
950 const bool j_non_zero_i =
951 support_on_face(j, face_n);
952 const bool j_non_zero_e =
953 support_on_face(j, neighbor_face_n);
954 if (flux_dof_mask(i, j) ==
always)
956 sparsity.add(dofs_on_this_cell[i],
957 dofs_on_other_cell[j]);
958 sparsity.add(dofs_on_other_cell[i],
959 dofs_on_this_cell[j]);
960 sparsity.add(dofs_on_this_cell[i],
961 dofs_on_this_cell[j]);
962 sparsity.add(dofs_on_other_cell[i],
963 dofs_on_other_cell[j]);
965 if (flux_dof_mask(i, j) ==
nonzero)
967 if (i_non_zero_i && j_non_zero_e)
968 sparsity.add(dofs_on_this_cell[i],
969 dofs_on_other_cell[j]);
970 if (i_non_zero_e && j_non_zero_i)
971 sparsity.add(dofs_on_other_cell[i],
972 dofs_on_this_cell[j]);
973 if (i_non_zero_i && j_non_zero_i)
974 sparsity.add(dofs_on_this_cell[i],
975 dofs_on_this_cell[j]);
976 if (i_non_zero_e && j_non_zero_e)
977 sparsity.add(dofs_on_other_cell[i],
978 dofs_on_other_cell[j]);
981 if (flux_dof_mask(j, i) ==
always)
983 sparsity.add(dofs_on_this_cell[j],
984 dofs_on_other_cell[i]);
985 sparsity.add(dofs_on_other_cell[j],
986 dofs_on_this_cell[i]);
987 sparsity.add(dofs_on_this_cell[j],
988 dofs_on_this_cell[i]);
989 sparsity.add(dofs_on_other_cell[j],
990 dofs_on_other_cell[i]);
992 if (flux_dof_mask(j, i) ==
nonzero)
994 if (j_non_zero_i && i_non_zero_e)
995 sparsity.add(dofs_on_this_cell[j],
996 dofs_on_other_cell[i]);
997 if (j_non_zero_e && i_non_zero_i)
998 sparsity.add(dofs_on_other_cell[j],
999 dofs_on_this_cell[i]);
1000 if (j_non_zero_i && i_non_zero_i)
1001 sparsity.add(dofs_on_this_cell[j],
1002 dofs_on_this_cell[i]);
1003 if (j_non_zero_e && i_non_zero_e)
1004 sparsity.add(dofs_on_other_cell[j],
1005 dofs_on_other_cell[i]);
1020 typename SparsityPatternType,
1024 const ::hp::DoFHandler<dim, spacedim> &dof,
1025 SparsityPatternType & sparsity,
1027 const bool keep_constrained_dofs,
1039 const ::hp::FECollection<dim, spacedim> &fe =
1040 dof.get_fe_collection();
1042 std::vector<types::global_dof_index> dofs_on_this_cell(
1044 std::vector<types::global_dof_index> dofs_on_other_cell(
1059 if (int_mask(c1, c2) !=
none || flux_mask(c1, c2) !=
none)
1060 int_and_flux_mask(c1, c2) =
always;
1062 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
1067 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
1068 for (
unsigned int f = 0; f < fe.size(); ++f)
1070 bool_int_and_flux_dof_mask[f].reinit(
1072 bool_int_and_flux_dof_mask[f].fill(
false);
1073 for (
unsigned int i = 0; i < fe[f].dofs_per_cell; ++i)
1074 for (
unsigned int j = 0; j < fe[f].dofs_per_cell; ++j)
1075 if (int_and_flux_dof_mask[f](i, j) !=
none)
1076 bool_int_and_flux_dof_mask[f](i, j) =
true;
1080 typename ::hp::DoFHandler<dim, spacedim>::active_cell_iterator
1081 cell = dof.begin_active(),
1083 for (; cell != endc; ++cell)
1085 (subdomain_id == cell->subdomain_id())) &&
1086 cell->is_locally_owned())
1088 dofs_on_this_cell.resize(cell->get_fe().dofs_per_cell);
1089 cell->get_dof_indices(dofs_on_this_cell);
1096 keep_constrained_dofs,
1097 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
1100 for (
unsigned int face = 0;
1101 face < GeometryInfo<dim>::faces_per_cell;
1104 const typename ::hp::DoFHandler<dim,
1105 spacedim>::face_iterator
1106 cell_face = cell->face(face);
1108 const bool periodic_neighbor =
1109 cell->has_periodic_neighbor(face);
1111 if ((!cell->at_boundary(face)) || periodic_neighbor)
1113 typename ::hp::DoFHandler<dim, spacedim>::
1114 level_cell_iterator neighbor =
1115 cell->neighbor_or_periodic_neighbor(face);
1121 if (neighbor->level() == cell->level() &&
1122 neighbor->index() > cell->index() &&
1123 neighbor->active() && neighbor->is_locally_owned())
1137 if (neighbor->level() != cell->level() &&
1138 ((!periodic_neighbor &&
1139 !cell->neighbor_is_coarser(face)) ||
1140 (periodic_neighbor &&
1141 !cell->periodic_neighbor_is_coarser(face))) &&
1142 neighbor->is_locally_owned())
1145 if (neighbor->has_children())
1147 for (
unsigned int sub_nr = 0;
1148 sub_nr != cell_face->n_children();
1151 const typename ::hp::DoFHandler<
1153 spacedim>::level_cell_iterator sub_neighbor =
1155 cell->periodic_neighbor_child_on_subface(
1157 cell->neighbor_child_on_subface(face, sub_nr);
1159 dofs_on_other_cell.resize(
1160 sub_neighbor->get_fe().dofs_per_cell);
1161 sub_neighbor->get_dof_indices(dofs_on_other_cell);
1162 for (
unsigned int i = 0;
1163 i < cell->get_fe().dofs_per_cell;
1166 const unsigned int ii =
1167 (cell->get_fe().is_primitive(i) ?
1169 .system_to_component_index(i)
1172 .get_nonzero_components(i)
1173 .first_selected_component());
1178 for (
unsigned int j = 0;
1179 j < sub_neighbor->get_fe().dofs_per_cell;
1182 const unsigned int jj =
1183 (sub_neighbor->get_fe().is_primitive(
1185 sub_neighbor->get_fe()
1186 .system_to_component_index(j)
1188 sub_neighbor->get_fe()
1189 .get_nonzero_components(j)
1190 .first_selected_component());
1192 Assert(jj < sub_neighbor->get_fe()
1196 if ((flux_mask(ii, jj) ==
always) ||
1197 (flux_mask(ii, jj) ==
nonzero))
1199 sparsity.add(dofs_on_this_cell[i],
1200 dofs_on_other_cell[j]);
1203 if ((flux_mask(jj, ii) ==
always) ||
1204 (flux_mask(jj, ii) ==
nonzero))
1206 sparsity.add(dofs_on_other_cell[j],
1207 dofs_on_this_cell[i]);
1215 dofs_on_other_cell.resize(
1216 neighbor->get_fe().dofs_per_cell);
1217 neighbor->get_dof_indices(dofs_on_other_cell);
1218 for (
unsigned int i = 0;
1219 i < cell->get_fe().dofs_per_cell;
1222 const unsigned int ii =
1223 (cell->get_fe().is_primitive(i) ?
1225 .system_to_component_index(i)
1228 .get_nonzero_components(i)
1229 .first_selected_component());
1234 for (
unsigned int j = 0;
1235 j < neighbor->get_fe().dofs_per_cell;
1238 const unsigned int jj =
1239 (neighbor->get_fe().is_primitive(j) ?
1241 .system_to_component_index(j)
1244 .get_nonzero_components(j)
1245 .first_selected_component());
1250 if ((flux_mask(ii, jj) ==
always) ||
1251 (flux_mask(ii, jj) ==
nonzero))
1253 sparsity.add(dofs_on_this_cell[i],
1254 dofs_on_other_cell[j]);
1257 if ((flux_mask(jj, ii) ==
always) ||
1258 (flux_mask(jj, ii) ==
nonzero))
1260 sparsity.add(dofs_on_other_cell[j],
1261 dofs_on_this_cell[i]);
1276 template <
typename DoFHandlerType,
typename SparsityPatternType>
1279 SparsityPatternType & sparsity,
1286 const bool keep_constrained_dofs =
true;
1291 keep_constrained_dofs,
1297 template <
typename DoFHandlerType,
1298 typename SparsityPatternType,
1302 SparsityPatternType & sparsity,
1304 const bool keep_constrained_dofs,
1313 const unsigned int n_comp = dof.get_fe(0).n_components();
1316 Assert(sparsity.n_rows() == n_dofs,
1318 Assert(sparsity.n_cols() == n_dofs,
1320 Assert(int_mask.n_rows() == n_comp,
1322 Assert(int_mask.n_cols() == n_comp,
1324 Assert(flux_mask.n_rows() == n_comp,
1326 Assert(flux_mask.n_cols() == n_comp,
1332 Assert((dof.get_triangulation().locally_owned_subdomain() ==
1336 dof.get_triangulation().locally_owned_subdomain()),
1338 "For parallel::distributed::Triangulation objects and " 1339 "associated DoF handler objects, asking for any subdomain other " 1340 "than the locally owned one does not make sense."));
1345 keep_constrained_dofs,
1357 #include "dof_tools_sparsity.inst" 1361 DEAL_II_NAMESPACE_CLOSE
#define AssertDimension(dim1, dim2)
const types::subdomain_id invalid_subdomain_id
void make_boundary_sparsity_pattern(const DoFHandlerType &dof, const std::vector< types::global_dof_index > &dof_to_boundary_mapping, SparsityPatternType &sparsity_pattern)
unsigned int size() const
bool is_primitive() const
unsigned long long int global_dof_index
static::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
size_type size(const unsigned int i) const
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void make_flux_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern)
const ComponentMask & get_nonzero_components(const unsigned int i) const
unsigned int n_components() const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const unsigned int dofs_per_cell
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
static::ExceptionBase & ExcInvalidBoundaryIndicator()
const types::boundary_id internal_face_boundary_id
const types::global_dof_index invalid_dof_index
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void fill(InputIterator entries, const bool C_style_indexing=true)
static::ExceptionBase & ExcInternalError()