Reference documentation for deal.II version GIT 85a786e179 2022-11-27 23:15: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 - 2022 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_tools.h>
30 #include <deal.II/fe/fe_values.h>
31 
34 #include <deal.II/grid/tria.h>
36 
38 #include <deal.II/hp/fe_values.h>
40 
43 #include <deal.II/lac/vector.h>
44 
45 #include <algorithm>
46 #include <numeric>
47 
49 
50 
51 
52 namespace DoFTools
53 {
54  template <int dim, int spacedim, typename number>
55  void
57  SparsityPatternBase & sparsity,
58  const AffineConstraints<number> &constraints,
59  const bool keep_constrained_dofs,
61  {
62  const types::global_dof_index n_dofs = dof.n_dofs();
63  (void)n_dofs;
64 
65  Assert(sparsity.n_rows() == n_dofs,
66  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
67  Assert(sparsity.n_cols() == n_dofs,
68  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
69 
70  // If we have a distributed Triangulation only allow locally_owned
71  // subdomain. Not setting a subdomain is also okay, because we skip
72  // ghost cells in the loop below.
73  if (const auto *triangulation = dynamic_cast<
75  &dof.get_triangulation()))
76  {
78  (subdomain_id == triangulation->locally_owned_subdomain()),
79  ExcMessage(
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."));
83  }
84 
85  std::vector<types::global_dof_index> dofs_on_this_cell;
86  dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
87 
88  // In case we work with a distributed sparsity pattern of Trilinos
89  // type, we only have to do the work if the current cell is owned by
90  // the calling processor. Otherwise, just continue.
91  for (const auto &cell : dof.active_cell_iterators())
93  (subdomain_id == cell->subdomain_id())) &&
94  cell->is_locally_owned())
95  {
96  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
97  dofs_on_this_cell.resize(dofs_per_cell);
98  cell->get_dof_indices(dofs_on_this_cell);
99 
100  // make sparsity pattern for this cell. if no constraints pattern
101  // was given, then the following call acts as if simply no
102  // constraints existed
103  constraints.add_entries_local_to_global(dofs_on_this_cell,
104  sparsity,
105  keep_constrained_dofs);
106  }
107  }
108 
109 
110 
111  template <int dim, int spacedim, typename number>
112  void
114  const Table<2, Coupling> & couplings,
115  SparsityPatternBase & sparsity,
116  const AffineConstraints<number> &constraints,
117  const bool keep_constrained_dofs,
119  {
120  const types::global_dof_index n_dofs = dof.n_dofs();
121  (void)n_dofs;
122 
123  Assert(sparsity.n_rows() == n_dofs,
124  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
125  Assert(sparsity.n_cols() == n_dofs,
126  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
127  Assert(couplings.n_rows() == dof.get_fe(0).n_components(),
128  ExcDimensionMismatch(couplings.n_rows(),
129  dof.get_fe(0).n_components()));
130  Assert(couplings.n_cols() == dof.get_fe(0).n_components(),
131  ExcDimensionMismatch(couplings.n_cols(),
132  dof.get_fe(0).n_components()));
133 
134  // If we have a distributed Triangulation only allow locally_owned
135  // subdomain. Not setting a subdomain is also okay, because we skip
136  // ghost cells in the loop below.
137  if (const auto *triangulation = dynamic_cast<
139  &dof.get_triangulation()))
140  {
142  (subdomain_id == triangulation->locally_owned_subdomain()),
143  ExcMessage(
144  "For distributed Triangulation objects and associated "
145  "DoFHandler objects, asking for any subdomain other than the "
146  "locally owned one does not make sense."));
147  }
148 
149  const hp::FECollection<dim, spacedim> &fe_collection =
150  dof.get_fe_collection();
151 
152  const std::vector<Table<2, Coupling>> dof_mask //(fe_collection.size())
153  = dof_couplings_from_component_couplings(fe_collection, couplings);
154 
155  // Convert the dof_mask to bool_dof_mask so we can pass it
156  // to constraints.add_entries_local_to_global()
157  std::vector<Table<2, bool>> bool_dof_mask(fe_collection.size());
158  for (unsigned int f = 0; f < fe_collection.size(); ++f)
159  {
160  bool_dof_mask[f].reinit(
161  TableIndices<2>(fe_collection[f].n_dofs_per_cell(),
162  fe_collection[f].n_dofs_per_cell()));
163  bool_dof_mask[f].fill(false);
164  for (unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
165  for (unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
166  if (dof_mask[f](i, j) != none)
167  bool_dof_mask[f](i, j) = true;
168  }
169 
170  std::vector<types::global_dof_index> dofs_on_this_cell(
171  fe_collection.max_dofs_per_cell());
172 
173  // In case we work with a distributed sparsity pattern of Trilinos
174  // type, we only have to do the work if the current cell is owned by
175  // the calling processor. Otherwise, just continue.
176  for (const auto &cell : dof.active_cell_iterators())
178  (subdomain_id == cell->subdomain_id())) &&
179  cell->is_locally_owned())
180  {
181  const types::fe_index fe_index = cell->active_fe_index();
182  const unsigned int dofs_per_cell =
183  fe_collection[fe_index].n_dofs_per_cell();
184 
185  dofs_on_this_cell.resize(dofs_per_cell);
186  cell->get_dof_indices(dofs_on_this_cell);
187 
188 
189  // make sparsity pattern for this cell. if no constraints pattern
190  // was given, then the following call acts as if simply no
191  // constraints existed
192  constraints.add_entries_local_to_global(dofs_on_this_cell,
193  sparsity,
194  keep_constrained_dofs,
195  bool_dof_mask[fe_index]);
196  }
197  }
198 
199 
200 
201  template <int dim, int spacedim>
202  void
204  const DoFHandler<dim, spacedim> &dof_col,
205  SparsityPatternBase & sparsity)
206  {
207  const types::global_dof_index n_dofs_row = dof_row.n_dofs();
208  const types::global_dof_index n_dofs_col = dof_col.n_dofs();
209  (void)n_dofs_row;
210  (void)n_dofs_col;
211 
212  Assert(sparsity.n_rows() == n_dofs_row,
213  ExcDimensionMismatch(sparsity.n_rows(), n_dofs_row));
214  Assert(sparsity.n_cols() == n_dofs_col,
215  ExcDimensionMismatch(sparsity.n_cols(), n_dofs_col));
216 
217  // It doesn't make sense to assemble sparsity patterns when the
218  // Triangulations are both parallel (i.e., different cells are assigned to
219  // different processors) and unequal: no processor will be responsible for
220  // assembling coupling terms between dofs on a cell owned by one processor
221  // and dofs on a cell owned by a different processor.
222  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
223  &dof_row.get_triangulation()) != nullptr ||
224  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
225  &dof_col.get_triangulation()) != nullptr)
226  {
227  Assert(&dof_row.get_triangulation() == &dof_col.get_triangulation(),
228  ExcMessage("This function can only be used with with parallel "
229  "Triangulations when the Triangulations are equal."));
230  }
231 
232  // TODO: Looks like wasteful memory management here
233 
234  using cell_iterator = typename DoFHandler<dim, spacedim>::cell_iterator;
235  std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
236  GridTools::get_finest_common_cells(dof_row, dof_col);
237 
238 #ifdef DEAL_II_WITH_MPI
239  // get_finest_common_cells returns all cells (locally owned and otherwise)
240  // for shared::Tria, but we don't want to assemble on cells that are not
241  // locally owned so remove them
242  if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
243  &dof_row.get_triangulation()) != nullptr ||
245  &dof_col.get_triangulation()) != nullptr)
246  {
247  const types::subdomain_id this_subdomain_id =
249  Assert(this_subdomain_id ==
251  ExcInternalError());
252  cell_list.erase(
253  std::remove_if(
254  cell_list.begin(),
255  cell_list.end(),
256  [=](const std::pair<cell_iterator, cell_iterator> &pair) {
257  return pair.first->subdomain_id() != this_subdomain_id ||
258  pair.second->subdomain_id() != this_subdomain_id;
259  }),
260  cell_list.end());
261  }
262 #endif
263 
264  for (const auto &cell_pair : cell_list)
265  {
266  const cell_iterator cell_row = cell_pair.first;
267  const cell_iterator cell_col = cell_pair.second;
268 
269  if (cell_row->is_active() && cell_col->is_active())
270  {
271  const unsigned int dofs_per_cell_row =
272  cell_row->get_fe().n_dofs_per_cell();
273  const unsigned int dofs_per_cell_col =
274  cell_col->get_fe().n_dofs_per_cell();
275  std::vector<types::global_dof_index> local_dof_indices_row(
276  dofs_per_cell_row);
277  std::vector<types::global_dof_index> local_dof_indices_col(
278  dofs_per_cell_col);
279  cell_row->get_dof_indices(local_dof_indices_row);
280  cell_col->get_dof_indices(local_dof_indices_col);
281  for (const auto &dof : local_dof_indices_row)
282  sparsity.add_row_entries(dof,
283  make_array_view(local_dof_indices_col));
284  }
285  else if (cell_row->has_children())
286  {
287  const std::vector<
289  child_cells =
290  GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
291  cell_row);
292  for (unsigned int i = 0; i < child_cells.size(); ++i)
293  {
295  cell_row_child = child_cells[i];
296  const unsigned int dofs_per_cell_row =
297  cell_row_child->get_fe().n_dofs_per_cell();
298  const unsigned int dofs_per_cell_col =
299  cell_col->get_fe().n_dofs_per_cell();
300  std::vector<types::global_dof_index> local_dof_indices_row(
301  dofs_per_cell_row);
302  std::vector<types::global_dof_index> local_dof_indices_col(
303  dofs_per_cell_col);
304  cell_row_child->get_dof_indices(local_dof_indices_row);
305  cell_col->get_dof_indices(local_dof_indices_col);
306  for (const auto &dof : local_dof_indices_row)
307  sparsity.add_row_entries(
308  dof, make_array_view(local_dof_indices_col));
309  }
310  }
311  else
312  {
313  std::vector<
315  child_cells =
316  GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
317  cell_col);
318  for (unsigned int i = 0; i < child_cells.size(); ++i)
319  {
321  cell_col_child = child_cells[i];
322  const unsigned int dofs_per_cell_row =
323  cell_row->get_fe().n_dofs_per_cell();
324  const unsigned int dofs_per_cell_col =
325  cell_col_child->get_fe().n_dofs_per_cell();
326  std::vector<types::global_dof_index> local_dof_indices_row(
327  dofs_per_cell_row);
328  std::vector<types::global_dof_index> local_dof_indices_col(
329  dofs_per_cell_col);
330  cell_row->get_dof_indices(local_dof_indices_row);
331  cell_col_child->get_dof_indices(local_dof_indices_col);
332  for (const auto &dof : local_dof_indices_row)
333  sparsity.add_row_entries(
334  dof, make_array_view(local_dof_indices_col));
335  }
336  }
337  }
338  }
339 
340 
341 
342  template <int dim, int spacedim>
343  void
345  const DoFHandler<dim, spacedim> & dof,
346  const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
347  SparsityPatternBase & sparsity)
348  {
349  if (dim == 1)
350  {
351  // there are only 2 boundary indicators in 1d, so it is no
352  // performance problem to call the other function
353  std::map<types::boundary_id, const Function<spacedim, double> *>
354  boundary_ids;
355  boundary_ids[0] = nullptr;
356  boundary_ids[1] = nullptr;
357  make_boundary_sparsity_pattern<dim, spacedim>(dof,
358  boundary_ids,
359  dof_to_boundary_mapping,
360  sparsity);
361  return;
362  }
363 
364  const types::global_dof_index n_dofs = dof.n_dofs();
365  (void)n_dofs;
366 
367  AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
368  AssertDimension(sparsity.n_rows(), dof.n_boundary_dofs());
369  AssertDimension(sparsity.n_cols(), dof.n_boundary_dofs());
370 #ifdef DEBUG
371  if (sparsity.n_rows() != 0)
372  {
373  types::global_dof_index max_element = 0;
374  for (const types::global_dof_index index : dof_to_boundary_mapping)
375  if ((index != numbers::invalid_dof_index) && (index > max_element))
376  max_element = index;
377  AssertDimension(max_element, sparsity.n_rows() - 1);
378  }
379 #endif
380 
381  std::vector<types::global_dof_index> dofs_on_this_face;
382  dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
383  std::vector<types::global_dof_index> cols;
384 
385  // loop over all faces to check whether they are at a boundary. note
386  // that we need not take special care of single lines (using
387  // @p{cell->has_boundary_lines}), since we do not support boundaries of
388  // dimension dim-2, and so every boundary line is also part of a
389  // boundary face.
390  for (const auto &cell : dof.active_cell_iterators())
391  for (const unsigned int f : cell->face_indices())
392  if (cell->at_boundary(f))
393  {
394  const unsigned int dofs_per_face =
395  cell->get_fe().n_dofs_per_face(f);
396  dofs_on_this_face.resize(dofs_per_face);
397  cell->face(f)->get_dof_indices(dofs_on_this_face,
398  cell->active_fe_index());
399 
400  // make sparsity pattern for this cell
401  cols.clear();
402  for (const auto &dof : dofs_on_this_face)
403  cols.push_back(dof_to_boundary_mapping[dof]);
404  // We are not guaranteed that the mapping to a second index space
405  // is increasing so sort here to use the faster add_row_entries()
406  // path
407  std::sort(cols.begin(), cols.end());
408  for (const auto &dof : dofs_on_this_face)
409  sparsity.add_row_entries(dof_to_boundary_mapping[dof],
410  make_array_view(cols),
411  true);
412  }
413  }
414 
415 
416 
417  template <int dim, int spacedim, typename number>
418  void
420  const DoFHandler<dim, spacedim> &dof,
421  const std::map<types::boundary_id, const Function<spacedim, number> *>
422  & boundary_ids,
423  const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
424  SparsityPatternBase & sparsity)
425  {
426  if (dim == 1)
427  {
428  // first check left, then right boundary point
429  for (unsigned int direction = 0; direction < 2; ++direction)
430  {
431  // if this boundary is not requested, then go on with next one
432  if (boundary_ids.find(direction) == boundary_ids.end())
433  continue;
434 
435  // find active cell at that boundary: first go to left/right,
436  // then to children
438  dof.begin(0);
439  while (!cell->at_boundary(direction))
440  cell = cell->neighbor(direction);
441  while (!cell->is_active())
442  cell = cell->child(direction);
443 
444  const unsigned int dofs_per_vertex =
445  cell->get_fe().n_dofs_per_vertex();
446  std::vector<types::global_dof_index> boundary_dof_boundary_indices(
447  dofs_per_vertex);
448 
449  // next get boundary mapped dof indices of boundary dofs
450  for (unsigned int i = 0; i < dofs_per_vertex; ++i)
451  boundary_dof_boundary_indices[i] =
452  dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
453 
454  std::sort(boundary_dof_boundary_indices.begin(),
455  boundary_dof_boundary_indices.end());
456  for (const auto &dof : boundary_dof_boundary_indices)
457  sparsity.add_row_entries(
458  dof, make_array_view(boundary_dof_boundary_indices), true);
459  }
460  return;
461  }
462 
463  const types::global_dof_index n_dofs = dof.n_dofs();
464  (void)n_dofs;
465 
466  AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
467  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
468  boundary_ids.end(),
470  Assert(sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
471  ExcDimensionMismatch(sparsity.n_rows(),
472  dof.n_boundary_dofs(boundary_ids)));
473  Assert(sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
474  ExcDimensionMismatch(sparsity.n_cols(),
475  dof.n_boundary_dofs(boundary_ids)));
476 #ifdef DEBUG
477  if (sparsity.n_rows() != 0)
478  {
479  types::global_dof_index max_element = 0;
480  for (const types::global_dof_index index : dof_to_boundary_mapping)
481  if ((index != numbers::invalid_dof_index) && (index > max_element))
482  max_element = index;
483  AssertDimension(max_element, sparsity.n_rows() - 1);
484  }
485 #endif
486 
487  std::vector<types::global_dof_index> dofs_on_this_face;
488  dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
489  std::vector<types::global_dof_index> cols;
490 
491  for (const auto &cell : dof.active_cell_iterators())
492  for (const unsigned int f : cell->face_indices())
493  if (boundary_ids.find(cell->face(f)->boundary_id()) !=
494  boundary_ids.end())
495  {
496  const unsigned int dofs_per_face =
497  cell->get_fe().n_dofs_per_face(f);
498  dofs_on_this_face.resize(dofs_per_face);
499  cell->face(f)->get_dof_indices(dofs_on_this_face,
500  cell->active_fe_index());
501 
502  // make sparsity pattern for this cell
503  cols.clear();
504  for (const auto &dof : dofs_on_this_face)
505  cols.push_back(dof_to_boundary_mapping[dof]);
506 
507  // Like the other one: sort once.
508  std::sort(cols.begin(), cols.end());
509  for (const auto &dof : dofs_on_this_face)
510  sparsity.add_row_entries(dof_to_boundary_mapping[dof],
511  make_array_view(cols),
512  true);
513  }
514  }
515 
516 
517 
518  template <int dim, int spacedim, typename number>
519  void
521  SparsityPatternBase & sparsity,
522  const AffineConstraints<number> &constraints,
523  const bool keep_constrained_dofs,
525 
526  // TODO: QA: reduce the indentation level of this method..., Maier 2012
527 
528  {
529  const types::global_dof_index n_dofs = dof.n_dofs();
530  (void)n_dofs;
531 
532  AssertDimension(sparsity.n_rows(), n_dofs);
533  AssertDimension(sparsity.n_cols(), n_dofs);
534 
535  // If we have a distributed Triangulation only allow locally_owned
536  // subdomain. Not setting a subdomain is also okay, because we skip
537  // ghost cells in the loop below.
538  if (const auto *triangulation = dynamic_cast<
540  &dof.get_triangulation()))
541  {
543  (subdomain_id == triangulation->locally_owned_subdomain()),
544  ExcMessage(
545  "For distributed Triangulation objects and associated "
546  "DoFHandler objects, asking for any subdomain other than the "
547  "locally owned one does not make sense."));
548  }
549 
550  std::vector<types::global_dof_index> dofs_on_this_cell;
551  std::vector<types::global_dof_index> dofs_on_other_cell;
552  dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
553  dofs_on_other_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
554 
555  // TODO: in an old implementation, we used user flags before to tag
556  // faces that were already touched. this way, we could reduce the work
557  // a little bit. now, we instead add only data from one side. this
558  // should be OK, but we need to actually verify it.
559 
560  // In case we work with a distributed sparsity pattern of Trilinos
561  // type, we only have to do the work if the current cell is owned by
562  // the calling processor. Otherwise, just continue.
563  for (const auto &cell : dof.active_cell_iterators())
565  (subdomain_id == cell->subdomain_id())) &&
566  cell->is_locally_owned())
567  {
568  const unsigned int n_dofs_on_this_cell =
569  cell->get_fe().n_dofs_per_cell();
570  dofs_on_this_cell.resize(n_dofs_on_this_cell);
571  cell->get_dof_indices(dofs_on_this_cell);
572 
573  // make sparsity pattern for this cell. if no constraints pattern
574  // was given, then the following call acts as if simply no
575  // constraints existed
576  constraints.add_entries_local_to_global(dofs_on_this_cell,
577  sparsity,
578  keep_constrained_dofs);
579 
580  for (const unsigned int face : cell->face_indices())
581  {
582  typename DoFHandler<dim, spacedim>::face_iterator cell_face =
583  cell->face(face);
584  const bool periodic_neighbor = cell->has_periodic_neighbor(face);
585  if (!cell->at_boundary(face) || periodic_neighbor)
586  {
588  neighbor = cell->neighbor_or_periodic_neighbor(face);
589 
590  // in 1d, we do not need to worry whether the neighbor
591  // might have children and then loop over those children.
592  // rather, we may as well go straight to the cell behind
593  // this particular cell's most terminal child
594  if (dim == 1)
595  while (neighbor->has_children())
596  neighbor = neighbor->child(face == 0 ? 1 : 0);
597 
598  if (neighbor->has_children())
599  {
600  for (unsigned int sub_nr = 0;
601  sub_nr != cell_face->n_active_descendants();
602  ++sub_nr)
603  {
604  const typename DoFHandler<dim, spacedim>::
605  level_cell_iterator sub_neighbor =
606  periodic_neighbor ?
607  cell->periodic_neighbor_child_on_subface(
608  face, sub_nr) :
609  cell->neighbor_child_on_subface(face, sub_nr);
610 
611  const unsigned int n_dofs_on_neighbor =
612  sub_neighbor->get_fe().n_dofs_per_cell();
613  dofs_on_other_cell.resize(n_dofs_on_neighbor);
614  sub_neighbor->get_dof_indices(dofs_on_other_cell);
615 
616  constraints.add_entries_local_to_global(
617  dofs_on_this_cell,
618  dofs_on_other_cell,
619  sparsity,
620  keep_constrained_dofs);
621  constraints.add_entries_local_to_global(
622  dofs_on_other_cell,
623  dofs_on_this_cell,
624  sparsity,
625  keep_constrained_dofs);
626  // only need to add this when the neighbor is not
627  // owned by the current processor, otherwise we add
628  // the entries for the neighbor there
629  if (sub_neighbor->subdomain_id() !=
630  cell->subdomain_id())
631  constraints.add_entries_local_to_global(
632  dofs_on_other_cell,
633  sparsity,
634  keep_constrained_dofs);
635  }
636  }
637  else
638  {
639  // Refinement edges are taken care of by coarser
640  // cells
641  if ((!periodic_neighbor &&
642  cell->neighbor_is_coarser(face)) ||
643  (periodic_neighbor &&
644  cell->periodic_neighbor_is_coarser(face)))
645  if (neighbor->subdomain_id() == cell->subdomain_id())
646  continue;
647 
648  const unsigned int n_dofs_on_neighbor =
649  neighbor->get_fe().n_dofs_per_cell();
650  dofs_on_other_cell.resize(n_dofs_on_neighbor);
651 
652  neighbor->get_dof_indices(dofs_on_other_cell);
653 
654  constraints.add_entries_local_to_global(
655  dofs_on_this_cell,
656  dofs_on_other_cell,
657  sparsity,
658  keep_constrained_dofs);
659 
660  // only need to add these in case the neighbor cell
661  // is not locally owned - otherwise, we touch each
662  // face twice and hence put the indices the other way
663  // around
664  if (!cell->neighbor_or_periodic_neighbor(face)
665  ->is_active() ||
666  (neighbor->subdomain_id() != cell->subdomain_id()))
667  {
668  constraints.add_entries_local_to_global(
669  dofs_on_other_cell,
670  dofs_on_this_cell,
671  sparsity,
672  keep_constrained_dofs);
673  if (neighbor->subdomain_id() != cell->subdomain_id())
674  constraints.add_entries_local_to_global(
675  dofs_on_other_cell,
676  sparsity,
677  keep_constrained_dofs);
678  }
679  }
680  }
681  }
682  }
683  }
684 
685 
686 
687  template <int dim, int spacedim>
688  void
690  SparsityPatternBase & sparsity)
691  {
693  make_flux_sparsity_pattern(dof, sparsity, dummy);
694  }
695 
696  template <int dim, int spacedim>
700  const Table<2, Coupling> & component_couplings)
701  {
702  Assert(component_couplings.n_rows() == fe.n_components(),
703  ExcDimensionMismatch(component_couplings.n_rows(),
704  fe.n_components()));
705  Assert(component_couplings.n_cols() == fe.n_components(),
706  ExcDimensionMismatch(component_couplings.n_cols(),
707  fe.n_components()));
708 
709  const unsigned int n_dofs = fe.n_dofs_per_cell();
710 
711  Table<2, Coupling> dof_couplings(n_dofs, n_dofs);
712 
713  for (unsigned int i = 0; i < n_dofs; ++i)
714  {
715  const unsigned int ii =
716  (fe.is_primitive(i) ?
717  fe.system_to_component_index(i).first :
719  Assert(ii < fe.n_components(), ExcInternalError());
720 
721  for (unsigned int j = 0; j < n_dofs; ++j)
722  {
723  const unsigned int jj =
724  (fe.is_primitive(j) ?
725  fe.system_to_component_index(j).first :
727  Assert(jj < fe.n_components(), ExcInternalError());
728 
729  dof_couplings(i, j) = component_couplings(ii, jj);
730  }
731  }
732  return dof_couplings;
733  }
734 
735 
736 
737  template <int dim, int spacedim>
738  std::vector<Table<2, Coupling>>
741  const Table<2, Coupling> & component_couplings)
742  {
743  std::vector<Table<2, Coupling>> return_value(fe.size());
744  for (unsigned int i = 0; i < fe.size(); ++i)
745  return_value[i] =
746  dof_couplings_from_component_couplings(fe[i], component_couplings);
747 
748  return return_value;
749  }
750 
751 
752 
753  namespace internal
754  {
755  namespace
756  {
757  // implementation of the same function in namespace DoFTools for
758  // non-hp-DoFHandlers
759  template <int dim, int spacedim, typename number>
760  void
762  const DoFHandler<dim, spacedim> &dof,
763  SparsityPatternBase & sparsity,
764  const AffineConstraints<number> &constraints,
765  const bool keep_constrained_dofs,
766  const Table<2, Coupling> & int_mask,
767  const Table<2, Coupling> & flux_mask,
769  const std::function<
771  const unsigned int)> &face_has_flux_coupling)
772  {
773  std::vector<types::global_dof_index> rows;
774  std::vector<std::pair<SparsityPatternBase::size_type,
776  cell_entries;
777 
778  if (dof.has_hp_capabilities() == false)
779  {
780  const FiniteElement<dim, spacedim> &fe = dof.get_fe();
781 
782  std::vector<types::global_dof_index> dofs_on_this_cell(
783  fe.n_dofs_per_cell());
784  std::vector<types::global_dof_index> dofs_on_other_cell(
785  fe.n_dofs_per_cell());
786 
787  const Table<2, Coupling>
788  int_dof_mask =
790  flux_dof_mask =
792 
793  Table<2, bool> support_on_face(fe.n_dofs_per_cell(),
795  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
796  for (const unsigned int f : GeometryInfo<dim>::face_indices())
797  support_on_face(i, f) = fe.has_support_on_face(i, f);
798 
799  // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
800  // to constraints.add_entries_local_to_global()
801  Table<2, bool> bool_int_dof_mask(fe.n_dofs_per_cell(),
802  fe.n_dofs_per_cell());
803  bool_int_dof_mask.fill(false);
804  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
805  for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
806  if (int_dof_mask(i, j) != none)
807  bool_int_dof_mask(i, j) = true;
808 
809  for (const auto &cell : dof.active_cell_iterators())
811  (subdomain_id == cell->subdomain_id())) &&
812  cell->is_locally_owned())
813  {
814  cell->get_dof_indices(dofs_on_this_cell);
815  // make sparsity pattern for this cell
816  constraints.add_entries_local_to_global(dofs_on_this_cell,
817  sparsity,
818  keep_constrained_dofs,
819  bool_int_dof_mask);
820  // Loop over all interior neighbors
821  for (const unsigned int face_n : cell->face_indices())
822  {
824  cell_face = cell->face(face_n);
825 
826  const bool periodic_neighbor =
827  cell->has_periodic_neighbor(face_n);
828 
829  if (cell->at_boundary(face_n) && (!periodic_neighbor))
830  {
831  for (unsigned int i = 0; i < fe.n_dofs_per_cell();
832  ++i)
833  {
834  const bool i_non_zero_i =
835  support_on_face(i, face_n);
836  for (unsigned int j = 0; j < fe.n_dofs_per_cell();
837  ++j)
838  {
839  const bool j_non_zero_i =
840  support_on_face(j, face_n);
841 
842  if (flux_dof_mask(i, j) == always ||
843  (flux_dof_mask(i, j) == nonzero &&
844  i_non_zero_i && j_non_zero_i))
845  cell_entries.emplace_back(
846  dofs_on_this_cell[i],
847  dofs_on_this_cell[j]);
848  }
849  }
850  sparsity.add_entries(make_array_view(cell_entries));
851  cell_entries.clear();
852  }
853  else
854  {
855  if (!face_has_flux_coupling(cell, face_n))
856  continue;
857 
858  typename DoFHandler<dim,
859  spacedim>::level_cell_iterator
860  neighbor =
861  cell->neighbor_or_periodic_neighbor(face_n);
862  // If the cells are on the same level (and both are
863  // active, locally-owned cells) then only add to the
864  // sparsity pattern if the current cell is 'greater'
865  // in the total ordering.
866  if (neighbor->level() == cell->level() &&
867  neighbor->index() > cell->index() &&
868  neighbor->is_active() &&
869  neighbor->is_locally_owned())
870  continue;
871  // If we are more refined then the neighbor, then we
872  // will automatically find the active neighbor cell
873  // when we call 'neighbor (face_n)' above. The
874  // opposite is not true; if the neighbor is more
875  // refined then the call 'neighbor (face_n)' will
876  // *not* return an active cell. Hence, only add things
877  // to the sparsity pattern if (when the levels are
878  // different) the neighbor is coarser than the current
879  // cell.
880  //
881  // Like above, do not use this optimization if the
882  // neighbor is not locally owned.
883  if (neighbor->level() != cell->level() &&
884  ((!periodic_neighbor &&
885  !cell->neighbor_is_coarser(face_n)) ||
886  (periodic_neighbor &&
887  !cell->periodic_neighbor_is_coarser(face_n))) &&
888  neighbor->is_locally_owned())
889  continue; // (the neighbor is finer)
890 
891  const unsigned int neighbor_face_n =
892  periodic_neighbor ?
893  cell->periodic_neighbor_face_no(face_n) :
894  cell->neighbor_face_no(face_n);
895 
896 
897  // In 1D, go straight to the cell behind this
898  // particular cell's most terminal cell. This makes us
899  // skip the if (neighbor->has_children()) section
900  // below. We need to do this since we otherwise
901  // iterate over the children of the face, which are
902  // always 0 in 1D.
903  if (dim == 1)
904  while (neighbor->has_children())
905  neighbor = neighbor->child(face_n == 0 ? 1 : 0);
906 
907  if (neighbor->has_children())
908  {
909  for (unsigned int sub_nr = 0;
910  sub_nr != cell_face->n_children();
911  ++sub_nr)
912  {
913  const typename DoFHandler<dim, spacedim>::
914  level_cell_iterator sub_neighbor =
915  periodic_neighbor ?
916  cell
917  ->periodic_neighbor_child_on_subface(
918  face_n, sub_nr) :
919  cell->neighbor_child_on_subface(face_n,
920  sub_nr);
921 
922  sub_neighbor->get_dof_indices(
923  dofs_on_other_cell);
924  for (unsigned int i = 0;
925  i < fe.n_dofs_per_cell();
926  ++i)
927  {
928  const bool i_non_zero_i =
929  support_on_face(i, face_n);
930  const bool i_non_zero_e =
931  support_on_face(i, neighbor_face_n);
932  for (unsigned int j = 0;
933  j < fe.n_dofs_per_cell();
934  ++j)
935  {
936  const bool j_non_zero_i =
937  support_on_face(j, face_n);
938  const bool j_non_zero_e =
939  support_on_face(j, neighbor_face_n);
940 
941  if (flux_dof_mask(i, j) == always)
942  {
943  cell_entries.emplace_back(
944  dofs_on_this_cell[i],
945  dofs_on_other_cell[j]);
946  cell_entries.emplace_back(
947  dofs_on_other_cell[i],
948  dofs_on_this_cell[j]);
949  cell_entries.emplace_back(
950  dofs_on_this_cell[i],
951  dofs_on_this_cell[j]);
952  cell_entries.emplace_back(
953  dofs_on_other_cell[i],
954  dofs_on_other_cell[j]);
955  }
956  else if (flux_dof_mask(i, j) ==
957  nonzero)
958  {
959  if (i_non_zero_i && j_non_zero_e)
960  cell_entries.emplace_back(
961  dofs_on_this_cell[i],
962  dofs_on_other_cell[j]);
963  if (i_non_zero_e && j_non_zero_i)
964  cell_entries.emplace_back(
965  dofs_on_other_cell[i],
966  dofs_on_this_cell[j]);
967  if (i_non_zero_i && j_non_zero_i)
968  cell_entries.emplace_back(
969  dofs_on_this_cell[i],
970  dofs_on_this_cell[j]);
971  if (i_non_zero_e && j_non_zero_e)
972  cell_entries.emplace_back(
973  dofs_on_other_cell[i],
974  dofs_on_other_cell[j]);
975  }
976 
977  if (flux_dof_mask(j, i) == always)
978  {
979  cell_entries.emplace_back(
980  dofs_on_this_cell[j],
981  dofs_on_other_cell[i]);
982  cell_entries.emplace_back(
983  dofs_on_other_cell[j],
984  dofs_on_this_cell[i]);
985  cell_entries.emplace_back(
986  dofs_on_this_cell[j],
987  dofs_on_this_cell[i]);
988  cell_entries.emplace_back(
989  dofs_on_other_cell[j],
990  dofs_on_other_cell[i]);
991  }
992  else if (flux_dof_mask(j, i) ==
993  nonzero)
994  {
995  if (j_non_zero_i && i_non_zero_e)
996  cell_entries.emplace_back(
997  dofs_on_this_cell[j],
998  dofs_on_other_cell[i]);
999  if (j_non_zero_e && i_non_zero_i)
1000  cell_entries.emplace_back(
1001  dofs_on_other_cell[j],
1002  dofs_on_this_cell[i]);
1003  if (j_non_zero_i && i_non_zero_i)
1004  cell_entries.emplace_back(
1005  dofs_on_this_cell[j],
1006  dofs_on_this_cell[i]);
1007  if (j_non_zero_e && i_non_zero_e)
1008  cell_entries.emplace_back(
1009  dofs_on_other_cell[j],
1010  dofs_on_other_cell[i]);
1011  }
1012  }
1013  }
1014  }
1015  sparsity.add_entries(
1016  make_array_view(cell_entries));
1017  cell_entries.clear();
1018  }
1019  else
1020  {
1021  neighbor->get_dof_indices(dofs_on_other_cell);
1022  for (unsigned int i = 0; i < fe.n_dofs_per_cell();
1023  ++i)
1024  {
1025  const bool i_non_zero_i =
1026  support_on_face(i, face_n);
1027  const bool i_non_zero_e =
1028  support_on_face(i, neighbor_face_n);
1029  for (unsigned int j = 0;
1030  j < fe.n_dofs_per_cell();
1031  ++j)
1032  {
1033  const bool j_non_zero_i =
1034  support_on_face(j, face_n);
1035  const bool j_non_zero_e =
1036  support_on_face(j, neighbor_face_n);
1037  if (flux_dof_mask(i, j) == always)
1038  {
1039  cell_entries.emplace_back(
1040  dofs_on_this_cell[i],
1041  dofs_on_other_cell[j]);
1042  cell_entries.emplace_back(
1043  dofs_on_other_cell[i],
1044  dofs_on_this_cell[j]);
1045  cell_entries.emplace_back(
1046  dofs_on_this_cell[i],
1047  dofs_on_this_cell[j]);
1048  cell_entries.emplace_back(
1049  dofs_on_other_cell[i],
1050  dofs_on_other_cell[j]);
1051  }
1052  if (flux_dof_mask(i, j) == nonzero)
1053  {
1054  if (i_non_zero_i && j_non_zero_e)
1055  cell_entries.emplace_back(
1056  dofs_on_this_cell[i],
1057  dofs_on_other_cell[j]);
1058  if (i_non_zero_e && j_non_zero_i)
1059  cell_entries.emplace_back(
1060  dofs_on_other_cell[i],
1061  dofs_on_this_cell[j]);
1062  if (i_non_zero_i && j_non_zero_i)
1063  cell_entries.emplace_back(
1064  dofs_on_this_cell[i],
1065  dofs_on_this_cell[j]);
1066  if (i_non_zero_e && j_non_zero_e)
1067  cell_entries.emplace_back(
1068  dofs_on_other_cell[i],
1069  dofs_on_other_cell[j]);
1070  }
1071 
1072  if (flux_dof_mask(j, i) == always)
1073  {
1074  cell_entries.emplace_back(
1075  dofs_on_this_cell[j],
1076  dofs_on_other_cell[i]);
1077  cell_entries.emplace_back(
1078  dofs_on_other_cell[j],
1079  dofs_on_this_cell[i]);
1080  cell_entries.emplace_back(
1081  dofs_on_this_cell[j],
1082  dofs_on_this_cell[i]);
1083  cell_entries.emplace_back(
1084  dofs_on_other_cell[j],
1085  dofs_on_other_cell[i]);
1086  }
1087  if (flux_dof_mask(j, i) == nonzero)
1088  {
1089  if (j_non_zero_i && i_non_zero_e)
1090  cell_entries.emplace_back(
1091  dofs_on_this_cell[j],
1092  dofs_on_other_cell[i]);
1093  if (j_non_zero_e && i_non_zero_i)
1094  cell_entries.emplace_back(
1095  dofs_on_other_cell[j],
1096  dofs_on_this_cell[i]);
1097  if (j_non_zero_i && i_non_zero_i)
1098  cell_entries.emplace_back(
1099  dofs_on_this_cell[j],
1100  dofs_on_this_cell[i]);
1101  if (j_non_zero_e && i_non_zero_e)
1102  cell_entries.emplace_back(
1103  dofs_on_other_cell[j],
1104  dofs_on_other_cell[i]);
1105  }
1106  }
1107  }
1108  sparsity.add_entries(
1109  make_array_view(cell_entries));
1110  cell_entries.clear();
1111  }
1112  }
1113  }
1114  }
1115  }
1116  else
1117  {
1118  // while the implementation above is quite optimized and caches a
1119  // lot of data (see e.g. the int/flux_dof_mask tables), this is no
1120  // longer practical for the hp-version since we would have to have
1121  // it for all combinations of elements in the hp::FECollection.
1122  // consequently, the implementation here is simpler and probably
1123  // less efficient but at least readable...
1124 
1125  const ::hp::FECollection<dim, spacedim> &fe =
1126  dof.get_fe_collection();
1127 
1128  std::vector<types::global_dof_index> dofs_on_this_cell(
1130  std::vector<types::global_dof_index> dofs_on_other_cell(
1132 
1133  const unsigned int n_components = fe.n_components();
1134  AssertDimension(int_mask.size(0), n_components);
1135  AssertDimension(int_mask.size(1), n_components);
1136  AssertDimension(flux_mask.size(0), n_components);
1137  AssertDimension(flux_mask.size(1), n_components);
1138 
1139  // note that we also need to set the respective entries if flux_mask
1140  // says so. this is necessary since we need to consider all degrees
1141  // of freedom on a cell for interior faces.
1142  Table<2, Coupling> int_and_flux_mask(n_components, n_components);
1143  for (unsigned int c1 = 0; c1 < n_components; ++c1)
1144  for (unsigned int c2 = 0; c2 < n_components; ++c2)
1145  if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
1146  int_and_flux_mask(c1, c2) = always;
1147 
1148  std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
1149  dof_couplings_from_component_couplings(fe, int_and_flux_mask);
1150 
1151  // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
1152  // can pass it to constraints.add_entries_local_to_global()
1153  std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
1154  for (unsigned int f = 0; f < fe.size(); ++f)
1155  {
1156  bool_int_and_flux_dof_mask[f].reinit(
1157  TableIndices<2>(fe[f].n_dofs_per_cell(),
1158  fe[f].n_dofs_per_cell()));
1159  bool_int_and_flux_dof_mask[f].fill(false);
1160  for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
1161  for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
1162  if (int_and_flux_dof_mask[f](i, j) != none)
1163  bool_int_and_flux_dof_mask[f](i, j) = true;
1164  }
1165 
1166 
1167  for (const auto &cell : dof.active_cell_iterators())
1169  (subdomain_id == cell->subdomain_id())) &&
1170  cell->is_locally_owned())
1171  {
1172  dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
1173  cell->get_dof_indices(dofs_on_this_cell);
1174 
1175  // make sparsity pattern for this cell also taking into
1176  // account the couplings due to face contributions on the same
1177  // cell
1178  constraints.add_entries_local_to_global(
1179  dofs_on_this_cell,
1180  sparsity,
1181  keep_constrained_dofs,
1182  bool_int_and_flux_dof_mask[cell->active_fe_index()]);
1183 
1184  // Loop over interior faces
1185  for (const unsigned int face : cell->face_indices())
1186  {
1187  const typename ::DoFHandler<dim,
1188  spacedim>::face_iterator
1189  cell_face = cell->face(face);
1190 
1191  const bool periodic_neighbor =
1192  cell->has_periodic_neighbor(face);
1193 
1194  if ((!cell->at_boundary(face)) || periodic_neighbor)
1195  {
1196  typename ::DoFHandler<dim, spacedim>::
1197  level_cell_iterator neighbor =
1198  cell->neighbor_or_periodic_neighbor(face);
1199 
1200  if (!face_has_flux_coupling(cell, face))
1201  continue;
1202 
1203  // Like the non-hp-case: If the cells are on the same
1204  // level (and both are active, locally-owned cells)
1205  // then only add to the sparsity pattern if the
1206  // current cell is 'greater' in the total ordering.
1207  if (neighbor->level() == cell->level() &&
1208  neighbor->index() > cell->index() &&
1209  neighbor->is_active() &&
1210  neighbor->is_locally_owned())
1211  continue;
1212  // Again, like the non-hp-case: If we are more refined
1213  // then the neighbor, then we will automatically find
1214  // the active neighbor cell when we call 'neighbor
1215  // (face)' above. The opposite is not true; if the
1216  // neighbor is more refined then the call 'neighbor
1217  // (face)' will *not* return an active cell. Hence,
1218  // only add things to the sparsity pattern if (when
1219  // the levels are different) the neighbor is coarser
1220  // than the current cell.
1221  //
1222  // Like above, do not use this optimization if the
1223  // neighbor is not locally owned.
1224  if (neighbor->level() != cell->level() &&
1225  ((!periodic_neighbor &&
1226  !cell->neighbor_is_coarser(face)) ||
1227  (periodic_neighbor &&
1228  !cell->periodic_neighbor_is_coarser(face))) &&
1229  neighbor->is_locally_owned())
1230  continue; // (the neighbor is finer)
1231 
1232  // In 1D, go straight to the cell behind this
1233  // particular cell's most terminal cell. This makes us
1234  // skip the if (neighbor->has_children()) section
1235  // below. We need to do this since we otherwise
1236  // iterate over the children of the face, which are
1237  // always 0 in 1D.
1238  if (dim == 1)
1239  while (neighbor->has_children())
1240  neighbor = neighbor->child(face == 0 ? 1 : 0);
1241 
1242  if (neighbor->has_children())
1243  {
1244  for (unsigned int sub_nr = 0;
1245  sub_nr != cell_face->n_children();
1246  ++sub_nr)
1247  {
1248  const typename ::DoFHandler<dim,
1249  spacedim>::
1250  level_cell_iterator sub_neighbor =
1251  periodic_neighbor ?
1252  cell
1253  ->periodic_neighbor_child_on_subface(
1254  face, sub_nr) :
1255  cell->neighbor_child_on_subface(face,
1256  sub_nr);
1257 
1258  dofs_on_other_cell.resize(
1259  sub_neighbor->get_fe().n_dofs_per_cell());
1260  sub_neighbor->get_dof_indices(
1261  dofs_on_other_cell);
1262  for (unsigned int i = 0;
1263  i < cell->get_fe().n_dofs_per_cell();
1264  ++i)
1265  {
1266  const unsigned int ii =
1267  (cell->get_fe().is_primitive(i) ?
1268  cell->get_fe()
1269  .system_to_component_index(i)
1270  .first :
1271  cell->get_fe()
1272  .get_nonzero_components(i)
1273  .first_selected_component());
1274 
1275  Assert(ii < cell->get_fe().n_components(),
1276  ExcInternalError());
1277 
1278  for (unsigned int j = 0;
1279  j < sub_neighbor->get_fe()
1280  .n_dofs_per_cell();
1281  ++j)
1282  {
1283  const unsigned int jj =
1284  (sub_neighbor->get_fe()
1285  .is_primitive(j) ?
1286  sub_neighbor->get_fe()
1287  .system_to_component_index(j)
1288  .first :
1289  sub_neighbor->get_fe()
1290  .get_nonzero_components(j)
1291  .first_selected_component());
1292 
1293  Assert(jj < sub_neighbor->get_fe()
1294  .n_components(),
1295  ExcInternalError());
1296 
1297  if ((flux_mask(ii, jj) == always) ||
1298  (flux_mask(ii, jj) == nonzero))
1299  {
1300  cell_entries.emplace_back(
1301  dofs_on_this_cell[i],
1302  dofs_on_other_cell[j]);
1303  }
1304 
1305  if ((flux_mask(jj, ii) == always) ||
1306  (flux_mask(jj, ii) == nonzero))
1307  {
1308  cell_entries.emplace_back(
1309  dofs_on_other_cell[j],
1310  dofs_on_this_cell[i]);
1311  }
1312  }
1313  }
1314  }
1315  }
1316  else
1317  {
1318  dofs_on_other_cell.resize(
1319  neighbor->get_fe().n_dofs_per_cell());
1320  neighbor->get_dof_indices(dofs_on_other_cell);
1321  for (unsigned int i = 0;
1322  i < cell->get_fe().n_dofs_per_cell();
1323  ++i)
1324  {
1325  const unsigned int ii =
1326  (cell->get_fe().is_primitive(i) ?
1327  cell->get_fe()
1328  .system_to_component_index(i)
1329  .first :
1330  cell->get_fe()
1331  .get_nonzero_components(i)
1332  .first_selected_component());
1333 
1334  Assert(ii < cell->get_fe().n_components(),
1335  ExcInternalError());
1336 
1337  for (unsigned int j = 0;
1338  j < neighbor->get_fe().n_dofs_per_cell();
1339  ++j)
1340  {
1341  const unsigned int jj =
1342  (neighbor->get_fe().is_primitive(j) ?
1343  neighbor->get_fe()
1344  .system_to_component_index(j)
1345  .first :
1346  neighbor->get_fe()
1347  .get_nonzero_components(j)
1348  .first_selected_component());
1349 
1350  Assert(
1351  jj < neighbor->get_fe().n_components(),
1352  ExcInternalError());
1353 
1354  if ((flux_mask(ii, jj) == always) ||
1355  (flux_mask(ii, jj) == nonzero))
1356  {
1357  cell_entries.emplace_back(
1358  dofs_on_this_cell[i],
1359  dofs_on_other_cell[j]);
1360  }
1361 
1362  if ((flux_mask(jj, ii) == always) ||
1363  (flux_mask(jj, ii) == nonzero))
1364  {
1365  cell_entries.emplace_back(
1366  dofs_on_other_cell[j],
1367  dofs_on_this_cell[i]);
1368  }
1369  }
1370  }
1371  }
1372  }
1373  }
1374  sparsity.add_entries(make_array_view(cell_entries));
1375  cell_entries.clear();
1376  }
1377  }
1378  }
1379  } // namespace
1380 
1381  } // namespace internal
1382 
1383 
1384 
1385  template <int dim, int spacedim>
1386  void
1388  SparsityPatternBase & sparsity,
1389  const Table<2, Coupling> & int_mask,
1390  const Table<2, Coupling> & flux_mask,
1392  {
1394 
1395  const bool keep_constrained_dofs = true;
1396 
1398  sparsity,
1399  dummy,
1400  keep_constrained_dofs,
1401  int_mask,
1402  flux_mask,
1403  subdomain_id,
1404  internal::always_couple_on_faces<dim, spacedim>);
1405  }
1406 
1407 
1408 
1409  template <int dim, int spacedim, typename number>
1410  void
1412  const DoFHandler<dim, spacedim> &dof,
1413  SparsityPatternBase & sparsity,
1414  const AffineConstraints<number> &constraints,
1415  const bool keep_constrained_dofs,
1416  const Table<2, Coupling> & int_mask,
1417  const Table<2, Coupling> & flux_mask,
1419  const std::function<
1421  const unsigned int)> &face_has_flux_coupling)
1422  {
1423  // do the error checking and frame code here, and then pass on to more
1424  // specialized functions in the internal namespace
1425  const types::global_dof_index n_dofs = dof.n_dofs();
1426  (void)n_dofs;
1427  const unsigned int n_comp = dof.get_fe(0).n_components();
1428  (void)n_comp;
1429 
1430  Assert(sparsity.n_rows() == n_dofs,
1431  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1432  Assert(sparsity.n_cols() == n_dofs,
1433  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1434  Assert(int_mask.n_rows() == n_comp,
1435  ExcDimensionMismatch(int_mask.n_rows(), n_comp));
1436  Assert(int_mask.n_cols() == n_comp,
1437  ExcDimensionMismatch(int_mask.n_cols(), n_comp));
1438  Assert(flux_mask.n_rows() == n_comp,
1439  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
1440  Assert(flux_mask.n_cols() == n_comp,
1441  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
1442 
1443  // If we have a distributed Triangulation only allow locally_owned
1444  // subdomain. Not setting a subdomain is also okay, because we skip
1445  // ghost cells in the loop below.
1446  if (const auto *triangulation = dynamic_cast<
1448  &dof.get_triangulation()))
1449  {
1451  (subdomain_id == triangulation->locally_owned_subdomain()),
1452  ExcMessage(
1453  "For distributed Triangulation objects and associated "
1454  "DoFHandler objects, asking for any subdomain other than the "
1455  "locally owned one does not make sense."));
1456  }
1457 
1458  Assert(
1459  face_has_flux_coupling,
1460  ExcMessage(
1461  "The function which specifies if a flux coupling occurs over a given "
1462  "face is empty."));
1463 
1465  sparsity,
1466  constraints,
1467  keep_constrained_dofs,
1468  int_mask,
1469  flux_mask,
1470  subdomain_id,
1471  face_has_flux_coupling);
1472  }
1473 
1474 } // end of namespace DoFTools
1475 
1476 
1477 // --------------------------------------------------- explicit instantiations
1478 
1479 #include "dof_tools_sparsity.inst"
1480 
1481 
1482 
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:699
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
bool has_hp_capabilities() 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
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) 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:264
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
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:1501
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
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:270
const types::subdomain_id invalid_subdomain_id
Definition: types.h:291
const types::global_dof_index invalid_dof_index
Definition: types.h:226
unsigned int global_dof_index
Definition: types.h:81
unsigned int subdomain_id
Definition: types.h:43
unsigned short int fe_index
Definition: types.h:59
unsigned int boundary_id
Definition: types.h:134
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation