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