Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_tools_sparsity.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/table.h>
19 #include <deal.II/base/utilities.h>
20 
23 
26 #include <deal.II/dofs/dof_tools.h>
27 
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_hermite.h>
30 #include <deal.II/fe/fe_tools.h>
31 #include <deal.II/fe/fe_values.h>
32 
35 #include <deal.II/grid/tria.h>
37 
39 #include <deal.II/hp/fe_values.h>
41 
44 #include <deal.II/lac/vector.h>
45 
46 #include <algorithm>
47 #include <numeric>
48 
50 
51 
52 
53 namespace DoFTools
54 {
55  template <int dim, int spacedim, typename number>
56  void
58  SparsityPatternBase &sparsity,
59  const AffineConstraints<number> &constraints,
60  const bool keep_constrained_dofs,
62  {
63  const types::global_dof_index n_dofs = dof.n_dofs();
64  (void)n_dofs;
65 
66  Assert(sparsity.n_rows() == n_dofs,
67  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
68  Assert(sparsity.n_cols() == n_dofs,
69  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
70 
71  // If we have a distributed Triangulation only allow locally_owned
72  // subdomain. Not setting a subdomain is also okay, because we skip
73  // ghost cells in the loop below.
74  if (const auto *triangulation = dynamic_cast<
76  &dof.get_triangulation()))
77  {
79  (subdomain_id == triangulation->locally_owned_subdomain()),
80  ExcMessage(
81  "For distributed Triangulation objects and associated "
82  "DoFHandler objects, asking for any subdomain other than the "
83  "locally owned one does not make sense."));
84  }
85 
86  std::vector<types::global_dof_index> dofs_on_this_cell;
87  dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
88 
89  // In case we work with a distributed sparsity pattern of Trilinos
90  // type, we only have to do the work if the current cell is owned by
91  // the calling processor. Otherwise, just continue.
92  for (const auto &cell : dof.active_cell_iterators())
94  (subdomain_id == cell->subdomain_id())) &&
95  cell->is_locally_owned())
96  {
97  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
98  dofs_on_this_cell.resize(dofs_per_cell);
99  cell->get_dof_indices(dofs_on_this_cell);
100 
101  // make sparsity pattern for this cell. if no constraints pattern
102  // was given, then the following call acts as if simply no
103  // constraints existed
104  constraints.add_entries_local_to_global(dofs_on_this_cell,
105  sparsity,
106  keep_constrained_dofs);
107  }
108  }
109 
110 
111 
112  template <int dim, int spacedim, typename number>
113  void
115  const Table<2, Coupling> &couplings,
116  SparsityPatternBase &sparsity,
117  const AffineConstraints<number> &constraints,
118  const bool keep_constrained_dofs,
120  {
121  const types::global_dof_index n_dofs = dof.n_dofs();
122  (void)n_dofs;
123 
124  Assert(sparsity.n_rows() == n_dofs,
125  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
126  Assert(sparsity.n_cols() == n_dofs,
127  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
128  Assert(couplings.n_rows() == dof.get_fe(0).n_components(),
129  ExcDimensionMismatch(couplings.n_rows(),
130  dof.get_fe(0).n_components()));
131  Assert(couplings.n_cols() == dof.get_fe(0).n_components(),
132  ExcDimensionMismatch(couplings.n_cols(),
133  dof.get_fe(0).n_components()));
134 
135  // If we have a distributed Triangulation only allow locally_owned
136  // subdomain. Not setting a subdomain is also okay, because we skip
137  // ghost cells in the loop below.
138  if (const auto *triangulation = dynamic_cast<
140  &dof.get_triangulation()))
141  {
143  (subdomain_id == triangulation->locally_owned_subdomain()),
144  ExcMessage(
145  "For distributed Triangulation objects and associated "
146  "DoFHandler objects, asking for any subdomain other than the "
147  "locally owned one does not make sense."));
148  }
149 
150  const hp::FECollection<dim, spacedim> &fe_collection =
151  dof.get_fe_collection();
152 
153  const std::vector<Table<2, Coupling>> dof_mask //(fe_collection.size())
154  = dof_couplings_from_component_couplings(fe_collection, couplings);
155 
156  // Convert the dof_mask to bool_dof_mask so we can pass it
157  // to constraints.add_entries_local_to_global()
158  std::vector<Table<2, bool>> bool_dof_mask(fe_collection.size());
159  for (unsigned int f = 0; f < fe_collection.size(); ++f)
160  {
161  bool_dof_mask[f].reinit(
162  TableIndices<2>(fe_collection[f].n_dofs_per_cell(),
163  fe_collection[f].n_dofs_per_cell()));
164  bool_dof_mask[f].fill(false);
165  for (unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
166  for (unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
167  if (dof_mask[f](i, j) != none)
168  bool_dof_mask[f](i, j) = true;
169  }
170 
171  std::vector<types::global_dof_index> dofs_on_this_cell(
172  fe_collection.max_dofs_per_cell());
173 
174  // In case we work with a distributed sparsity pattern of Trilinos
175  // type, we only have to do the work if the current cell is owned by
176  // the calling processor. Otherwise, just continue.
177  for (const auto &cell : dof.active_cell_iterators())
179  (subdomain_id == cell->subdomain_id())) &&
180  cell->is_locally_owned())
181  {
182  const types::fe_index fe_index = cell->active_fe_index();
183  const unsigned int dofs_per_cell =
184  fe_collection[fe_index].n_dofs_per_cell();
185 
186  dofs_on_this_cell.resize(dofs_per_cell);
187  cell->get_dof_indices(dofs_on_this_cell);
188 
189 
190  // make sparsity pattern for this cell. if no constraints pattern
191  // was given, then the following call acts as if simply no
192  // constraints existed
193  constraints.add_entries_local_to_global(dofs_on_this_cell,
194  sparsity,
195  keep_constrained_dofs,
196  bool_dof_mask[fe_index]);
197  }
198  }
199 
200 
201 
202  template <int dim, int spacedim>
203  void
205  const DoFHandler<dim, spacedim> &dof_col,
206  SparsityPatternBase &sparsity)
207  {
208  const types::global_dof_index n_dofs_row = dof_row.n_dofs();
209  const types::global_dof_index n_dofs_col = dof_col.n_dofs();
210  (void)n_dofs_row;
211  (void)n_dofs_col;
212 
213  Assert(sparsity.n_rows() == n_dofs_row,
214  ExcDimensionMismatch(sparsity.n_rows(), n_dofs_row));
215  Assert(sparsity.n_cols() == n_dofs_col,
216  ExcDimensionMismatch(sparsity.n_cols(), n_dofs_col));
217 
218  // It doesn't make sense to assemble sparsity patterns when the
219  // Triangulations are both parallel (i.e., different cells are assigned to
220  // different processors) and unequal: no processor will be responsible for
221  // assembling coupling terms between dofs on a cell owned by one processor
222  // and dofs on a cell owned by a different processor.
223  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
224  &dof_row.get_triangulation()) != nullptr ||
225  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
226  &dof_col.get_triangulation()) != nullptr)
227  {
228  Assert(&dof_row.get_triangulation() == &dof_col.get_triangulation(),
229  ExcMessage("This function can only be used with with parallel "
230  "Triangulations when the Triangulations are equal."));
231  }
232 
233  // TODO: Looks like wasteful memory management here
234 
235  using cell_iterator = typename DoFHandler<dim, spacedim>::cell_iterator;
236  std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
237  GridTools::get_finest_common_cells(dof_row, dof_col);
238 
239 #ifdef DEAL_II_WITH_MPI
240  // get_finest_common_cells returns all cells (locally owned and otherwise)
241  // for shared::Tria, but we don't want to assemble on cells that are not
242  // locally owned so remove them
243  if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
244  &dof_row.get_triangulation()) != nullptr ||
246  &dof_col.get_triangulation()) != nullptr)
247  {
248  const types::subdomain_id this_subdomain_id =
250  Assert(this_subdomain_id ==
252  ExcInternalError());
253  cell_list.erase(
254  std::remove_if(
255  cell_list.begin(),
256  cell_list.end(),
257  [=](const std::pair<cell_iterator, cell_iterator> &pair) {
258  return pair.first->subdomain_id() != this_subdomain_id ||
259  pair.second->subdomain_id() != this_subdomain_id;
260  }),
261  cell_list.end());
262  }
263 #endif
264 
265  for (const auto &cell_pair : cell_list)
266  {
267  const cell_iterator cell_row = cell_pair.first;
268  const cell_iterator cell_col = cell_pair.second;
269 
270  if (cell_row->is_active() && cell_col->is_active())
271  {
272  const unsigned int dofs_per_cell_row =
273  cell_row->get_fe().n_dofs_per_cell();
274  const unsigned int dofs_per_cell_col =
275  cell_col->get_fe().n_dofs_per_cell();
276  std::vector<types::global_dof_index> local_dof_indices_row(
277  dofs_per_cell_row);
278  std::vector<types::global_dof_index> local_dof_indices_col(
279  dofs_per_cell_col);
280  cell_row->get_dof_indices(local_dof_indices_row);
281  cell_col->get_dof_indices(local_dof_indices_col);
282  for (const auto &dof : local_dof_indices_row)
283  sparsity.add_row_entries(dof,
284  make_array_view(local_dof_indices_col));
285  }
286  else if (cell_row->has_children())
287  {
288  const std::vector<
290  child_cells =
291  GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
292  cell_row);
293  for (unsigned int i = 0; i < child_cells.size(); ++i)
294  {
296  cell_row_child = child_cells[i];
297  const unsigned int dofs_per_cell_row =
298  cell_row_child->get_fe().n_dofs_per_cell();
299  const unsigned int dofs_per_cell_col =
300  cell_col->get_fe().n_dofs_per_cell();
301  std::vector<types::global_dof_index> local_dof_indices_row(
302  dofs_per_cell_row);
303  std::vector<types::global_dof_index> local_dof_indices_col(
304  dofs_per_cell_col);
305  cell_row_child->get_dof_indices(local_dof_indices_row);
306  cell_col->get_dof_indices(local_dof_indices_col);
307  for (const auto &dof : local_dof_indices_row)
308  sparsity.add_row_entries(
309  dof, make_array_view(local_dof_indices_col));
310  }
311  }
312  else
313  {
314  std::vector<
316  child_cells =
317  GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
318  cell_col);
319  for (unsigned int i = 0; i < child_cells.size(); ++i)
320  {
322  cell_col_child = child_cells[i];
323  const unsigned int dofs_per_cell_row =
324  cell_row->get_fe().n_dofs_per_cell();
325  const unsigned int dofs_per_cell_col =
326  cell_col_child->get_fe().n_dofs_per_cell();
327  std::vector<types::global_dof_index> local_dof_indices_row(
328  dofs_per_cell_row);
329  std::vector<types::global_dof_index> local_dof_indices_col(
330  dofs_per_cell_col);
331  cell_row->get_dof_indices(local_dof_indices_row);
332  cell_col_child->get_dof_indices(local_dof_indices_col);
333  for (const auto &dof : local_dof_indices_row)
334  sparsity.add_row_entries(
335  dof, make_array_view(local_dof_indices_col));
336  }
337  }
338  }
339  }
340 
341 
342 
343  template <int dim, int spacedim>
344  void
346  const DoFHandler<dim, spacedim> &dof,
347  const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
348  SparsityPatternBase &sparsity)
349  {
350  if (dim == 1)
351  {
352  // there are only 2 boundary indicators in 1d, so it is no
353  // performance problem to call the other function
354  std::map<types::boundary_id, const Function<spacedim, double> *>
355  boundary_ids;
356  boundary_ids[0] = nullptr;
357  boundary_ids[1] = nullptr;
358  make_boundary_sparsity_pattern<dim, spacedim>(dof,
359  boundary_ids,
360  dof_to_boundary_mapping,
361  sparsity);
362  return;
363  }
364 
365  const types::global_dof_index n_dofs = dof.n_dofs();
366  (void)n_dofs;
367 
368  AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
369  AssertDimension(sparsity.n_rows(), dof.n_boundary_dofs());
370  AssertDimension(sparsity.n_cols(), dof.n_boundary_dofs());
371 #ifdef DEBUG
372  if (sparsity.n_rows() != 0)
373  {
374  types::global_dof_index max_element = 0;
375  for (const types::global_dof_index index : dof_to_boundary_mapping)
376  if ((index != numbers::invalid_dof_index) && (index > max_element))
377  max_element = index;
378  AssertDimension(max_element, sparsity.n_rows() - 1);
379  }
380 #endif
381 
382  std::vector<types::global_dof_index> dofs_on_this_face;
383  dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
384  std::vector<types::global_dof_index> cols;
385 
386  // loop over all faces to check whether they are at a boundary. note
387  // that we need not take special care of single lines (using
388  // @p{cell->has_boundary_lines}), since we do not support boundaries of
389  // dimension dim-2, and so every boundary line is also part of a
390  // boundary face.
391  for (const auto &cell : dof.active_cell_iterators())
392  for (const unsigned int f : cell->face_indices())
393  if (cell->at_boundary(f))
394  {
395  const unsigned int dofs_per_face =
396  cell->get_fe().n_dofs_per_face(f);
397  dofs_on_this_face.resize(dofs_per_face);
398  cell->face(f)->get_dof_indices(dofs_on_this_face,
399  cell->active_fe_index());
400 
401  // make sparsity pattern for this cell
402  cols.clear();
403  for (const auto &dof : dofs_on_this_face)
404  cols.push_back(dof_to_boundary_mapping[dof]);
405  // We are not guaranteed that the mapping to a second index space
406  // is increasing so sort here to use the faster add_row_entries()
407  // path
408  std::sort(cols.begin(), cols.end());
409  for (const auto &dof : dofs_on_this_face)
410  sparsity.add_row_entries(dof_to_boundary_mapping[dof],
411  make_array_view(cols),
412  true);
413  }
414  }
415 
416 
417 
418  template <int dim, int spacedim, typename number>
419  void
421  const DoFHandler<dim, spacedim> &dof,
422  const std::map<types::boundary_id, const Function<spacedim, number> *>
423  &boundary_ids,
424  const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
425  SparsityPatternBase &sparsity)
426  {
427  if (dim == 1)
428  {
429  // first check left, then right boundary point
430  for (unsigned int direction = 0; direction < 2; ++direction)
431  {
432  // if this boundary is not requested, then go on with next one
433  if (boundary_ids.find(direction) == boundary_ids.end())
434  continue;
435 
436  // find active cell at that boundary: first go to left/right,
437  // then to children
439  dof.begin(0);
440  while (!cell->at_boundary(direction))
441  cell = cell->neighbor(direction);
442  while (!cell->is_active())
443  cell = cell->child(direction);
444 
445  const unsigned int dofs_per_vertex =
446  cell->get_fe().n_dofs_per_vertex();
447  std::vector<types::global_dof_index> boundary_dof_boundary_indices(
448  dofs_per_vertex);
449 
450  // next get boundary mapped dof indices of boundary dofs
451  for (unsigned int i = 0; i < dofs_per_vertex; ++i)
452  boundary_dof_boundary_indices[i] =
453  dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
454 
455  std::sort(boundary_dof_boundary_indices.begin(),
456  boundary_dof_boundary_indices.end());
457  for (const auto &dof : boundary_dof_boundary_indices)
458  sparsity.add_row_entries(
459  dof, make_array_view(boundary_dof_boundary_indices), true);
460  }
461  return;
462  }
463 
464  const types::global_dof_index n_dofs = dof.n_dofs();
465  (void)n_dofs;
466 
467  AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
468  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
469  boundary_ids.end(),
471 
472  const bool fe_is_hermite = (dynamic_cast<const FE_Hermite<dim, spacedim> *>(
473  &(dof.get_fe())) != nullptr);
474 
475  Assert(fe_is_hermite ||
476  sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
477  ExcDimensionMismatch(sparsity.n_rows(),
478  dof.n_boundary_dofs(boundary_ids)));
479  Assert(fe_is_hermite ||
480  sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
481  ExcDimensionMismatch(sparsity.n_cols(),
482  dof.n_boundary_dofs(boundary_ids)));
483  (void)fe_is_hermite;
484 
485 #ifdef DEBUG
486  if (sparsity.n_rows() != 0)
487  {
488  types::global_dof_index max_element = 0;
489  for (const types::global_dof_index index : dof_to_boundary_mapping)
490  if ((index != numbers::invalid_dof_index) && (index > max_element))
491  max_element = index;
492  AssertDimension(max_element, sparsity.n_rows() - 1);
493  }
494 #endif
495 
496  std::vector<types::global_dof_index> dofs_on_this_face;
497  dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
498  std::vector<types::global_dof_index> cols;
499 
500  for (const auto &cell : dof.active_cell_iterators())
501  for (const unsigned int f : cell->face_indices())
502  if (boundary_ids.find(cell->face(f)->boundary_id()) !=
503  boundary_ids.end())
504  {
505  const unsigned int dofs_per_face =
506  cell->get_fe().n_dofs_per_face(f);
507  dofs_on_this_face.resize(dofs_per_face);
508  cell->face(f)->get_dof_indices(dofs_on_this_face,
509  cell->active_fe_index());
510 
511  // make sparsity pattern for this cell
512  cols.clear();
513  for (const auto &dof : dofs_on_this_face)
514  cols.push_back(dof_to_boundary_mapping[dof]);
515 
516  // Like the other one: sort once.
517  std::sort(cols.begin(), cols.end());
518  for (const auto &dof : dofs_on_this_face)
519  sparsity.add_row_entries(dof_to_boundary_mapping[dof],
520  make_array_view(cols),
521  true);
522  }
523  }
524 
525 
526 
527  template <int dim, int spacedim, typename number>
528  void
530  SparsityPatternBase &sparsity,
531  const AffineConstraints<number> &constraints,
532  const bool keep_constrained_dofs,
534 
535  // TODO: QA: reduce the indentation level of this method..., Maier 2012
536 
537  {
538  const types::global_dof_index n_dofs = dof.n_dofs();
539  (void)n_dofs;
540 
541  AssertDimension(sparsity.n_rows(), n_dofs);
542  AssertDimension(sparsity.n_cols(), n_dofs);
543 
544  // If we have a distributed Triangulation only allow locally_owned
545  // subdomain. Not setting a subdomain is also okay, because we skip
546  // ghost cells in the loop below.
547  if (const auto *triangulation = dynamic_cast<
549  &dof.get_triangulation()))
550  {
552  (subdomain_id == triangulation->locally_owned_subdomain()),
553  ExcMessage(
554  "For distributed Triangulation objects and associated "
555  "DoFHandler objects, asking for any subdomain other than the "
556  "locally owned one does not make sense."));
557  }
558 
559  std::vector<types::global_dof_index> dofs_on_this_cell;
560  std::vector<types::global_dof_index> dofs_on_other_cell;
561  dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
562  dofs_on_other_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
563 
564  // TODO: in an old implementation, we used user flags before to tag
565  // faces that were already touched. this way, we could reduce the work
566  // a little bit. now, we instead add only data from one side. this
567  // should be OK, but we need to actually verify it.
568 
569  // In case we work with a distributed sparsity pattern of Trilinos
570  // type, we only have to do the work if the current cell is owned by
571  // the calling processor. Otherwise, just continue.
572  for (const auto &cell : dof.active_cell_iterators())
574  (subdomain_id == cell->subdomain_id())) &&
575  cell->is_locally_owned())
576  {
577  const unsigned int n_dofs_on_this_cell =
578  cell->get_fe().n_dofs_per_cell();
579  dofs_on_this_cell.resize(n_dofs_on_this_cell);
580  cell->get_dof_indices(dofs_on_this_cell);
581 
582  // make sparsity pattern for this cell. if no constraints pattern
583  // was given, then the following call acts as if simply no
584  // constraints existed
585  constraints.add_entries_local_to_global(dofs_on_this_cell,
586  sparsity,
587  keep_constrained_dofs);
588 
589  for (const unsigned int face : cell->face_indices())
590  {
591  typename DoFHandler<dim, spacedim>::face_iterator cell_face =
592  cell->face(face);
593  const bool periodic_neighbor = cell->has_periodic_neighbor(face);
594  if (!cell->at_boundary(face) || periodic_neighbor)
595  {
597  neighbor = cell->neighbor_or_periodic_neighbor(face);
598 
599  // in 1d, we do not need to worry whether the neighbor
600  // might have children and then loop over those children.
601  // rather, we may as well go straight to the cell behind
602  // this particular cell's most terminal child
603  if (dim == 1)
604  while (neighbor->has_children())
605  neighbor = neighbor->child(face == 0 ? 1 : 0);
606 
607  if (neighbor->has_children())
608  {
609  for (unsigned int sub_nr = 0;
610  sub_nr != cell_face->n_active_descendants();
611  ++sub_nr)
612  {
613  const typename DoFHandler<dim, spacedim>::
614  level_cell_iterator sub_neighbor =
615  periodic_neighbor ?
616  cell->periodic_neighbor_child_on_subface(
617  face, sub_nr) :
618  cell->neighbor_child_on_subface(face, sub_nr);
619 
620  const unsigned int n_dofs_on_neighbor =
621  sub_neighbor->get_fe().n_dofs_per_cell();
622  dofs_on_other_cell.resize(n_dofs_on_neighbor);
623  sub_neighbor->get_dof_indices(dofs_on_other_cell);
624 
625  constraints.add_entries_local_to_global(
626  dofs_on_this_cell,
627  dofs_on_other_cell,
628  sparsity,
629  keep_constrained_dofs);
630  constraints.add_entries_local_to_global(
631  dofs_on_other_cell,
632  dofs_on_this_cell,
633  sparsity,
634  keep_constrained_dofs);
635  // only need to add this when the neighbor is not
636  // owned by the current processor, otherwise we add
637  // the entries for the neighbor there
638  if (sub_neighbor->subdomain_id() !=
639  cell->subdomain_id())
640  constraints.add_entries_local_to_global(
641  dofs_on_other_cell,
642  sparsity,
643  keep_constrained_dofs);
644  }
645  }
646  else
647  {
648  // Refinement edges are taken care of by coarser
649  // cells
650  if ((!periodic_neighbor &&
651  cell->neighbor_is_coarser(face)) ||
652  (periodic_neighbor &&
653  cell->periodic_neighbor_is_coarser(face)))
654  if (neighbor->subdomain_id() == cell->subdomain_id())
655  continue;
656 
657  const unsigned int n_dofs_on_neighbor =
658  neighbor->get_fe().n_dofs_per_cell();
659  dofs_on_other_cell.resize(n_dofs_on_neighbor);
660 
661  neighbor->get_dof_indices(dofs_on_other_cell);
662 
663  constraints.add_entries_local_to_global(
664  dofs_on_this_cell,
665  dofs_on_other_cell,
666  sparsity,
667  keep_constrained_dofs);
668 
669  // only need to add these in case the neighbor cell
670  // is not locally owned - otherwise, we touch each
671  // face twice and hence put the indices the other way
672  // around
673  if (!cell->neighbor_or_periodic_neighbor(face)
674  ->is_active() ||
675  (neighbor->subdomain_id() != cell->subdomain_id()))
676  {
677  constraints.add_entries_local_to_global(
678  dofs_on_other_cell,
679  dofs_on_this_cell,
680  sparsity,
681  keep_constrained_dofs);
682  if (neighbor->subdomain_id() != cell->subdomain_id())
683  constraints.add_entries_local_to_global(
684  dofs_on_other_cell,
685  sparsity,
686  keep_constrained_dofs);
687  }
688  }
689  }
690  }
691  }
692  }
693 
694 
695 
696  template <int dim, int spacedim>
697  void
699  SparsityPatternBase &sparsity)
700  {
702  make_flux_sparsity_pattern(dof, sparsity, dummy);
703  }
704 
705  template <int dim, int spacedim>
709  const Table<2, Coupling> &component_couplings)
710  {
711  Assert(component_couplings.n_rows() == fe.n_components(),
712  ExcDimensionMismatch(component_couplings.n_rows(),
713  fe.n_components()));
714  Assert(component_couplings.n_cols() == fe.n_components(),
715  ExcDimensionMismatch(component_couplings.n_cols(),
716  fe.n_components()));
717 
718  const unsigned int n_dofs = fe.n_dofs_per_cell();
719 
720  Table<2, Coupling> dof_couplings(n_dofs, n_dofs);
721 
722  for (unsigned int i = 0; i < n_dofs; ++i)
723  {
724  const unsigned int ii =
725  (fe.is_primitive(i) ?
726  fe.system_to_component_index(i).first :
728  Assert(ii < fe.n_components(), ExcInternalError());
729 
730  for (unsigned int j = 0; j < n_dofs; ++j)
731  {
732  const unsigned int jj =
733  (fe.is_primitive(j) ?
734  fe.system_to_component_index(j).first :
736  Assert(jj < fe.n_components(), ExcInternalError());
737 
738  dof_couplings(i, j) = component_couplings(ii, jj);
739  }
740  }
741  return dof_couplings;
742  }
743 
744 
745 
746  template <int dim, int spacedim>
747  std::vector<Table<2, Coupling>>
750  const Table<2, Coupling> &component_couplings)
751  {
752  std::vector<Table<2, Coupling>> return_value(fe.size());
753  for (unsigned int i = 0; i < fe.size(); ++i)
754  return_value[i] =
755  dof_couplings_from_component_couplings(fe[i], component_couplings);
756 
757  return return_value;
758  }
759 
760 
761 
762  namespace internal
763  {
764  namespace
765  {
766  // helper function
767  template <typename Iterator, typename Iterator2>
768  void
769  add_cell_entries(
770  const Iterator &cell,
771  const unsigned int face_no,
772  const Iterator2 &neighbor,
773  const unsigned int neighbor_face_no,
774  const Table<2, Coupling> &flux_mask,
775  const std::vector<types::global_dof_index> &dofs_on_this_cell,
776  std::vector<types::global_dof_index> &dofs_on_other_cell,
777  std::vector<std::pair<SparsityPatternBase::size_type,
778  SparsityPatternBase::size_type>> &cell_entries)
779  {
780  dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
781  neighbor->get_dof_indices(dofs_on_other_cell);
782 
783  // Keep expensive data structures in separate vectors for inner j loop
784  // in separate vectors
785  boost::container::small_vector<unsigned int, 64>
786  component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
787  boost::container::small_vector<bool, 64> support_on_face_i(
788  neighbor->get_fe().n_dofs_per_cell());
789  boost::container::small_vector<bool, 64> support_on_face_e(
790  neighbor->get_fe().n_dofs_per_cell());
791  for (unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
792  {
793  component_indices_neighbor[j] =
794  (neighbor->get_fe().is_primitive(j) ?
795  neighbor->get_fe().system_to_component_index(j).first :
796  neighbor->get_fe()
797  .get_nonzero_components(j)
798  .first_selected_component());
799  support_on_face_i[j] =
800  neighbor->get_fe().has_support_on_face(j, face_no);
801  support_on_face_e[j] =
802  neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
803  }
804 
805  // For the parallel setting, must include also the diagonal
806  // neighbor-neighbor coupling, otherwise those get included on the
807  // other cell
808  for (int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
809  {
810  const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
811  const auto &dofs_i =
812  (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
813  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
814  {
815  const unsigned int ii =
816  (fe.is_primitive(i) ?
817  fe.system_to_component_index(i).first :
818  fe.get_nonzero_components(i).first_selected_component());
819 
820  Assert(ii < fe.n_components(), ExcInternalError());
821  const bool i_non_zero_i =
822  fe.has_support_on_face(i,
823  (f == 0 ? face_no : neighbor_face_no));
824 
825  for (unsigned int j = 0;
826  j < neighbor->get_fe().n_dofs_per_cell();
827  ++j)
828  {
829  const bool j_non_zero_e = support_on_face_e[j];
830  const unsigned int jj = component_indices_neighbor[j];
831 
832  Assert(jj < neighbor->get_fe().n_components(),
833  ExcInternalError());
834 
835  if ((flux_mask(ii, jj) == always) ||
836  (flux_mask(ii, jj) == nonzero && i_non_zero_i &&
837  j_non_zero_e))
838  cell_entries.emplace_back(dofs_i[i],
839  dofs_on_other_cell[j]);
840  if ((flux_mask(jj, ii) == always) ||
841  (flux_mask(jj, ii) == nonzero && j_non_zero_e &&
842  i_non_zero_i))
843  cell_entries.emplace_back(dofs_on_other_cell[j],
844  dofs_i[i]);
845  }
846  }
847  }
848  }
849 
850 
851 
852  // implementation of the same function in namespace DoFTools for
853  // non-hp-DoFHandlers
854  template <int dim, int spacedim, typename number>
855  void
857  const DoFHandler<dim, spacedim> &dof,
858  SparsityPatternBase &sparsity,
859  const AffineConstraints<number> &constraints,
860  const bool keep_constrained_dofs,
861  const Table<2, Coupling> &int_mask,
862  const Table<2, Coupling> &flux_mask,
864  const std::function<
866  const unsigned int)> &face_has_flux_coupling)
867  {
868  std::vector<types::global_dof_index> rows;
869  std::vector<std::pair<SparsityPatternBase::size_type,
871  cell_entries;
872 
873  const ::hp::FECollection<dim, spacedim> &fe =
874  dof.get_fe_collection();
875 
876  std::vector<types::global_dof_index> dofs_on_this_cell(
878  std::vector<types::global_dof_index> dofs_on_other_cell(
880 
881  const unsigned int n_components = fe.n_components();
882  AssertDimension(int_mask.size(0), n_components);
883  AssertDimension(int_mask.size(1), n_components);
884  AssertDimension(flux_mask.size(0), n_components);
885  AssertDimension(flux_mask.size(1), n_components);
886 
887  // note that we also need to set the respective entries if flux_mask
888  // says so. this is necessary since we need to consider all degrees
889  // of freedom on a cell for interior faces.
890  Table<2, Coupling> int_and_flux_mask(n_components, n_components);
891  for (unsigned int c1 = 0; c1 < n_components; ++c1)
892  for (unsigned int c2 = 0; c2 < n_components; ++c2)
893  if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
894  int_and_flux_mask(c1, c2) = always;
895 
896  // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
897  // to constraints.add_entries_local_to_global()
898  std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
899  dof_couplings_from_component_couplings(fe, int_and_flux_mask);
900 
901  // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
902  // can pass it to constraints.add_entries_local_to_global()
903  std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
904  for (unsigned int f = 0; f < fe.size(); ++f)
905  {
906  bool_int_and_flux_dof_mask[f].reinit(
907  TableIndices<2>(fe[f].n_dofs_per_cell(),
908  fe[f].n_dofs_per_cell()));
909  bool_int_and_flux_dof_mask[f].fill(false);
910  for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
911  for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
912  if (int_and_flux_dof_mask[f](i, j) != none)
913  bool_int_and_flux_dof_mask[f](i, j) = true;
914  }
915 
916 
917  for (const auto &cell : dof.active_cell_iterators())
919  (subdomain_id == cell->subdomain_id())) &&
920  cell->is_locally_owned())
921  {
922  dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
923  cell->get_dof_indices(dofs_on_this_cell);
924 
925  // make sparsity pattern for this cell also taking into
926  // account the couplings due to face contributions on the same
927  // cell
928  constraints.add_entries_local_to_global(
929  dofs_on_this_cell,
930  sparsity,
931  keep_constrained_dofs,
932  bool_int_and_flux_dof_mask[cell->active_fe_index()]);
933 
934  // Loop over interior faces
935  for (const unsigned int face : cell->face_indices())
936  {
937  const bool periodic_neighbor =
938  cell->has_periodic_neighbor(face);
939 
940  if ((!cell->at_boundary(face)) || periodic_neighbor)
941  {
943  neighbor = cell->neighbor_or_periodic_neighbor(face);
944 
945  // If the cells are on the same level (and both are
946  // active, locally-owned cells) then only add to the
947  // sparsity pattern if the current cell is 'greater' in
948  // the total ordering.
949  if (neighbor->level() == cell->level() &&
950  neighbor->index() > cell->index() &&
951  neighbor->is_active() && neighbor->is_locally_owned())
952  continue;
953 
954  // If we are more refined then the neighbor, then we
955  // will automatically find the active neighbor cell when
956  // we call 'neighbor (face)' above. The opposite is not
957  // true; if the neighbor is more refined then the call
958  // 'neighbor (face)' will *not* return an active
959  // cell. Hence, only add things to the sparsity pattern
960  // if (when the levels are different) the neighbor is
961  // coarser than the current cell, except in the case
962  // when the neighbor is not locally owned.
963  if (neighbor->level() != cell->level() &&
964  ((!periodic_neighbor &&
965  !cell->neighbor_is_coarser(face)) ||
966  (periodic_neighbor &&
967  !cell->periodic_neighbor_is_coarser(face))) &&
968  neighbor->is_locally_owned())
969  continue; // (the neighbor is finer)
970 
971  if (!face_has_flux_coupling(cell, face))
972  continue;
973 
974  const unsigned int neighbor_face_no =
975  periodic_neighbor ?
976  cell->periodic_neighbor_face_no(face) :
977  cell->neighbor_face_no(face);
978 
979  // In 1d, go straight to the cell behind this
980  // particular cell's most terminal cell. This makes us
981  // skip the if (neighbor->has_children()) section
982  // below. We need to do this since we otherwise
983  // iterate over the children of the face, which are
984  // always 0 in 1d.
985  if (dim == 1)
986  while (neighbor->has_children())
987  neighbor = neighbor->child(face == 0 ? 1 : 0);
988 
989  if (neighbor->has_children())
990  {
991  for (unsigned int sub_nr = 0;
992  sub_nr != cell->face(face)->n_children();
993  ++sub_nr)
994  {
995  const typename DoFHandler<dim, spacedim>::
996  level_cell_iterator sub_neighbor =
997  periodic_neighbor ?
998  cell->periodic_neighbor_child_on_subface(
999  face, sub_nr) :
1000  cell->neighbor_child_on_subface(face,
1001  sub_nr);
1002  add_cell_entries(cell,
1003  face,
1004  sub_neighbor,
1005  neighbor_face_no,
1006  flux_mask,
1007  dofs_on_this_cell,
1008  dofs_on_other_cell,
1009  cell_entries);
1010  }
1011  }
1012  else
1013  add_cell_entries(cell,
1014  face,
1015  neighbor,
1016  neighbor_face_no,
1017  flux_mask,
1018  dofs_on_this_cell,
1019  dofs_on_other_cell,
1020  cell_entries);
1021  }
1022  }
1023  sparsity.add_entries(make_array_view(cell_entries));
1024  cell_entries.clear();
1025  }
1026  }
1027  } // namespace
1028 
1029  } // namespace internal
1030 
1031 
1032 
1033  template <int dim, int spacedim>
1034  void
1036  SparsityPatternBase &sparsity,
1037  const Table<2, Coupling> &int_mask,
1038  const Table<2, Coupling> &flux_mask,
1040  {
1042 
1043  const bool keep_constrained_dofs = true;
1044 
1046  sparsity,
1047  dummy,
1048  keep_constrained_dofs,
1049  int_mask,
1050  flux_mask,
1051  subdomain_id,
1052  internal::always_couple_on_faces<dim, spacedim>);
1053  }
1054 
1055 
1056 
1057  template <int dim, int spacedim, typename number>
1058  void
1060  const DoFHandler<dim, spacedim> &dof,
1061  SparsityPatternBase &sparsity,
1062  const AffineConstraints<number> &constraints,
1063  const bool keep_constrained_dofs,
1064  const Table<2, Coupling> &int_mask,
1065  const Table<2, Coupling> &flux_mask,
1067  const std::function<
1069  const unsigned int)> &face_has_flux_coupling)
1070  {
1071  // do the error checking and frame code here, and then pass on to more
1072  // specialized functions in the internal namespace
1073  const types::global_dof_index n_dofs = dof.n_dofs();
1074  (void)n_dofs;
1075  const unsigned int n_comp = dof.get_fe(0).n_components();
1076  (void)n_comp;
1077 
1078  Assert(sparsity.n_rows() == n_dofs,
1079  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1080  Assert(sparsity.n_cols() == n_dofs,
1081  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1082  Assert(int_mask.n_rows() == n_comp,
1083  ExcDimensionMismatch(int_mask.n_rows(), n_comp));
1084  Assert(int_mask.n_cols() == n_comp,
1085  ExcDimensionMismatch(int_mask.n_cols(), n_comp));
1086  Assert(flux_mask.n_rows() == n_comp,
1087  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
1088  Assert(flux_mask.n_cols() == n_comp,
1089  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
1090 
1091  // If we have a distributed Triangulation only allow locally_owned
1092  // subdomain. Not setting a subdomain is also okay, because we skip
1093  // ghost cells in the loop below.
1094  if (const auto *triangulation = dynamic_cast<
1096  &dof.get_triangulation()))
1097  {
1099  (subdomain_id == triangulation->locally_owned_subdomain()),
1100  ExcMessage(
1101  "For distributed Triangulation objects and associated "
1102  "DoFHandler objects, asking for any subdomain other than the "
1103  "locally owned one does not make sense."));
1104  }
1105 
1106  Assert(
1107  face_has_flux_coupling,
1108  ExcMessage(
1109  "The function which specifies if a flux coupling occurs over a given "
1110  "face is empty."));
1111 
1113  sparsity,
1114  constraints,
1115  keep_constrained_dofs,
1116  int_mask,
1117  flux_mask,
1118  subdomain_id,
1119  face_has_flux_coupling);
1120  }
1121 
1122 } // end of namespace DoFTools
1123 
1124 
1125 // --------------------------------------------------- explicit instantiations
1126 
1127 #include "dof_tools_sparsity.inst"
1128 
1129 
1130 
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:950
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
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
cell_iterator begin(const unsigned int level=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
types::global_dof_index size_type
size_type n_rows() const
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
size_type n_cols() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int size() const
Definition: collection.h:265
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
Table< 2, Coupling > dof_couplings_from_component_couplings(const FiniteElement< dim, spacedim > &fe, const Table< 2, Coupling > &component_couplings)
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
const types::boundary_id internal_face_boundary_id
Definition: types.h:313
const types::subdomain_id invalid_subdomain_id
Definition: types.h:342
const types::global_dof_index invalid_dof_index
Definition: types.h:253
unsigned int global_dof_index
Definition: types.h:82
unsigned int subdomain_id
Definition: types.h:44
unsigned short int fe_index
Definition: types.h:60
unsigned int boundary_id
Definition: types.h:145
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation