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