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