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