deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
dof_tools_sparsity.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2013 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
14#include <deal.II/base/table.h>
17
20
24
25#include <deal.II/fe/fe.h>
27#include <deal.II/fe/fe_tools.h>
29
32#include <deal.II/grid/tria.h>
34
38
41#include <deal.II/lac/vector.h>
42
43#include <algorithm>
44#include <complex>
45#include <numeric>
46
48
49
50
51namespace DoFTools
52{
53 template <int dim, int spacedim, typename number>
54 void
56 SparsityPatternBase &sparsity,
57 const AffineConstraints<number> &constraints,
58 const bool keep_constrained_dofs,
59 const types::subdomain_id subdomain_id)
60 {
61 const types::global_dof_index n_dofs = dof.n_dofs();
62 (void)n_dofs;
63
64 Assert(sparsity.n_rows() == n_dofs,
65 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
66 Assert(sparsity.n_cols() == n_dofs,
67 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
68
69 // If we have a distributed Triangulation only allow locally_owned
70 // subdomain. Not setting a subdomain is also okay, because we skip
71 // ghost cells in the loop below.
72 if (const auto *triangulation = dynamic_cast<
74 &dof.get_triangulation()))
75 {
76 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
77 (subdomain_id == triangulation->locally_owned_subdomain()),
79 "For distributed Triangulation objects and associated "
80 "DoFHandler objects, asking for any subdomain other than the "
81 "locally owned one does not make sense."));
82 }
83
84 const auto &fe_collection = dof.get_fe_collection();
85 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
86 for (unsigned int f = 0; f < fe_collection.size(); ++f)
87 {
88 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
89 }
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
94 // In case we work with a distributed sparsity pattern of Trilinos
95 // type, we only have to do the work if the current cell is owned by
96 // the calling processor. Otherwise, just continue.
97 for (const auto &cell : dof.active_cell_iterators())
98 if (((subdomain_id == numbers::invalid_subdomain_id) ||
99 (subdomain_id == cell->subdomain_id())) &&
100 cell->is_locally_owned())
101 {
102 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
103 dofs_on_this_cell.resize(dofs_per_cell);
104 cell->get_dof_indices(dofs_on_this_cell);
105
106 // make sparsity pattern for this cell. if no constraints pattern
107 // was given, then the following call acts as if simply no
108 // constraints existed
109 const types::fe_index fe_index = cell->active_fe_index();
110 if (fe_dof_mask[fe_index].empty())
111 constraints.add_entries_local_to_global(dofs_on_this_cell,
112 sparsity,
113 keep_constrained_dofs);
114 else
115 constraints.add_entries_local_to_global(dofs_on_this_cell,
116 sparsity,
117 keep_constrained_dofs,
118 fe_dof_mask[fe_index]);
119 }
120 }
121
122
123
124 template <int dim, int spacedim, typename number>
125 void
127 const Table<2, Coupling> &couplings,
128 SparsityPatternBase &sparsity,
129 const AffineConstraints<number> &constraints,
130 const bool keep_constrained_dofs,
131 const types::subdomain_id subdomain_id)
132 {
133 const types::global_dof_index n_dofs = dof.n_dofs();
134 (void)n_dofs;
135
136 Assert(sparsity.n_rows() == n_dofs,
137 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
138 Assert(sparsity.n_cols() == n_dofs,
139 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
140 Assert(couplings.n_rows() == dof.get_fe(0).n_components(),
141 ExcDimensionMismatch(couplings.n_rows(),
142 dof.get_fe(0).n_components()));
143 Assert(couplings.n_cols() == dof.get_fe(0).n_components(),
144 ExcDimensionMismatch(couplings.n_cols(),
145 dof.get_fe(0).n_components()));
146
147 // If we have a distributed Triangulation only allow locally_owned
148 // subdomain. Not setting a subdomain is also okay, because we skip
149 // ghost cells in the loop below.
150 if (const auto *triangulation = dynamic_cast<
152 &dof.get_triangulation()))
153 {
154 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
155 (subdomain_id == triangulation->locally_owned_subdomain()),
157 "For distributed Triangulation objects and associated "
158 "DoFHandler objects, asking for any subdomain other than the "
159 "locally owned one does not make sense."));
160 }
161
162 const hp::FECollection<dim, spacedim> &fe_collection =
163 dof.get_fe_collection();
164
165 const std::vector<Table<2, Coupling>> dof_mask //(fe_collection.size())
166 = dof_couplings_from_component_couplings(fe_collection, couplings);
167
168 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
169 for (unsigned int f = 0; f < fe_collection.size(); ++f)
170 {
171 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
172 }
173
174 // Convert the dof_mask to bool_dof_mask so we can pass it
175 // to constraints.add_entries_local_to_global()
176 std::vector<Table<2, bool>> bool_dof_mask(fe_collection.size());
177 for (unsigned int f = 0; f < fe_collection.size(); ++f)
178 {
179 bool_dof_mask[f].reinit(
180 TableIndices<2>(fe_collection[f].n_dofs_per_cell(),
181 fe_collection[f].n_dofs_per_cell()));
182 bool_dof_mask[f].fill(false);
183 for (unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
184 for (unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
185 if (dof_mask[f](i, j) != none &&
186 (fe_dof_mask[f].empty() || fe_dof_mask[f](i, j)))
187 bool_dof_mask[f](i, j) = true;
188 }
189
190 std::vector<types::global_dof_index> dofs_on_this_cell(
191 fe_collection.max_dofs_per_cell());
192
193 // In case we work with a distributed sparsity pattern of Trilinos
194 // type, we only have to do the work if the current cell is owned by
195 // the calling processor. Otherwise, just continue.
196 for (const auto &cell : dof.active_cell_iterators())
197 if (((subdomain_id == numbers::invalid_subdomain_id) ||
198 (subdomain_id == cell->subdomain_id())) &&
199 cell->is_locally_owned())
200 {
201 const types::fe_index fe_index = cell->active_fe_index();
202 const unsigned int dofs_per_cell =
203 fe_collection[fe_index].n_dofs_per_cell();
204
205 dofs_on_this_cell.resize(dofs_per_cell);
206 cell->get_dof_indices(dofs_on_this_cell);
207
208
209 // make sparsity pattern for this cell. if no constraints pattern
210 // was given, then the following call acts as if simply no
211 // constraints existed
212 constraints.add_entries_local_to_global(dofs_on_this_cell,
213 sparsity,
214 keep_constrained_dofs,
215 bool_dof_mask[fe_index]);
216 }
217 }
218
219
220
221 template <int dim, int spacedim>
222 void
224 const DoFHandler<dim, spacedim> &dof_col,
225 SparsityPatternBase &sparsity)
226 {
227 const types::global_dof_index n_dofs_row = dof_row.n_dofs();
228 const types::global_dof_index n_dofs_col = dof_col.n_dofs();
229 (void)n_dofs_row;
230 (void)n_dofs_col;
231
232 Assert(sparsity.n_rows() == n_dofs_row,
233 ExcDimensionMismatch(sparsity.n_rows(), n_dofs_row));
234 Assert(sparsity.n_cols() == n_dofs_col,
235 ExcDimensionMismatch(sparsity.n_cols(), n_dofs_col));
236
237 // It doesn't make sense to assemble sparsity patterns when the
238 // Triangulations are both parallel (i.e., different cells are assigned to
239 // different processors) and unequal: no processor will be responsible for
240 // assembling coupling terms between dofs on a cell owned by one processor
241 // and dofs on a cell owned by a different processor.
242 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
243 &dof_row.get_triangulation()) != nullptr ||
244 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
245 &dof_col.get_triangulation()) != nullptr)
246 {
247 Assert(&dof_row.get_triangulation() == &dof_col.get_triangulation(),
248 ExcMessage("This function can only be used with with parallel "
249 "Triangulations when the Triangulations are equal."));
250 }
251
252 // TODO: Looks like wasteful memory management here
253
254 using cell_iterator = typename DoFHandler<dim, spacedim>::cell_iterator;
255 std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
256 GridTools::get_finest_common_cells(dof_row, dof_col);
257
258#ifdef DEAL_II_WITH_MPI
259 // get_finest_common_cells returns all cells (locally owned and otherwise)
260 // for shared::Tria, but we don't want to assemble on cells that are not
261 // locally owned so remove them
262 if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
263 &dof_row.get_triangulation()) != nullptr ||
265 &dof_col.get_triangulation()) != nullptr)
266 {
267 const types::subdomain_id this_subdomain_id =
269 Assert(this_subdomain_id ==
272 cell_list.erase(
273 std::remove_if(
274 cell_list.begin(),
275 cell_list.end(),
276 [=](const std::pair<cell_iterator, cell_iterator> &pair) {
277 return pair.first->subdomain_id() != this_subdomain_id ||
278 pair.second->subdomain_id() != this_subdomain_id;
279 }),
280 cell_list.end());
281 }
282#endif
283
284 for (const auto &cell_pair : cell_list)
285 {
286 const cell_iterator cell_row = cell_pair.first;
287 const cell_iterator cell_col = cell_pair.second;
288
289 if (cell_row->is_active() && cell_col->is_active())
290 {
291 const unsigned int dofs_per_cell_row =
292 cell_row->get_fe().n_dofs_per_cell();
293 const unsigned int dofs_per_cell_col =
294 cell_col->get_fe().n_dofs_per_cell();
295 std::vector<types::global_dof_index> local_dof_indices_row(
296 dofs_per_cell_row);
297 std::vector<types::global_dof_index> local_dof_indices_col(
298 dofs_per_cell_col);
299 cell_row->get_dof_indices(local_dof_indices_row);
300 cell_col->get_dof_indices(local_dof_indices_col);
301 for (const auto &dof : local_dof_indices_row)
302 sparsity.add_row_entries(dof,
303 make_array_view(local_dof_indices_col));
304 }
305 else if (cell_row->has_children())
306 {
307 const std::vector<
309 child_cells =
310 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
311 cell_row);
312 for (unsigned int i = 0; i < child_cells.size(); ++i)
313 {
315 cell_row_child = child_cells[i];
316 const unsigned int dofs_per_cell_row =
317 cell_row_child->get_fe().n_dofs_per_cell();
318 const unsigned int dofs_per_cell_col =
319 cell_col->get_fe().n_dofs_per_cell();
320 std::vector<types::global_dof_index> local_dof_indices_row(
321 dofs_per_cell_row);
322 std::vector<types::global_dof_index> local_dof_indices_col(
323 dofs_per_cell_col);
324 cell_row_child->get_dof_indices(local_dof_indices_row);
325 cell_col->get_dof_indices(local_dof_indices_col);
326 for (const auto &dof : local_dof_indices_row)
327 sparsity.add_row_entries(
328 dof, make_array_view(local_dof_indices_col));
329 }
330 }
331 else
332 {
333 std::vector<
335 child_cells =
336 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
337 cell_col);
338 for (unsigned int i = 0; i < child_cells.size(); ++i)
339 {
341 &cell_col_child = child_cells[i];
342 const unsigned int dofs_per_cell_row =
343 cell_row->get_fe().n_dofs_per_cell();
344 const unsigned int dofs_per_cell_col =
345 cell_col_child->get_fe().n_dofs_per_cell();
346 std::vector<types::global_dof_index> local_dof_indices_row(
347 dofs_per_cell_row);
348 std::vector<types::global_dof_index> local_dof_indices_col(
349 dofs_per_cell_col);
350 cell_row->get_dof_indices(local_dof_indices_row);
351 cell_col_child->get_dof_indices(local_dof_indices_col);
352 for (const auto &dof : local_dof_indices_row)
353 sparsity.add_row_entries(
354 dof, make_array_view(local_dof_indices_col));
355 }
356 }
357 }
358 }
359
360
361
362 template <int dim, int spacedim>
363 void
365 const DoFHandler<dim, spacedim> &dof,
366 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
367 SparsityPatternBase &sparsity)
368 {
369 if (dim == 1)
370 {
371 // there are only 2 boundary indicators in 1d, so it is no
372 // performance problem to call the other function
373 std::map<types::boundary_id, const Function<spacedim, double> *>
374 boundary_ids;
375 boundary_ids[0] = nullptr;
376 boundary_ids[1] = nullptr;
377 make_boundary_sparsity_pattern<dim, spacedim>(dof,
378 boundary_ids,
379 dof_to_boundary_mapping,
380 sparsity);
381 return;
382 }
383
384 const types::global_dof_index n_dofs = dof.n_dofs();
385 (void)n_dofs;
386
387 AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
388 AssertDimension(sparsity.n_rows(), dof.n_boundary_dofs());
389 AssertDimension(sparsity.n_cols(), dof.n_boundary_dofs());
390 if constexpr (running_in_debug_mode())
391 {
392 if (sparsity.n_rows() != 0)
393 {
394 types::global_dof_index max_element = 0;
395 for (const types::global_dof_index index : dof_to_boundary_mapping)
396 if ((index != numbers::invalid_dof_index) &&
397 (index > max_element))
398 max_element = index;
399 AssertDimension(max_element, sparsity.n_rows() - 1);
400 }
401 }
402
403 std::vector<types::global_dof_index> dofs_on_this_face;
404 dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
405 std::vector<types::global_dof_index> cols;
406
407 // loop over all faces to check whether they are at a boundary. note
408 // that we need not take special care of single lines (using
409 // @p{cell->has_boundary_lines}), since we do not support boundaries of
410 // dimension dim-2, and so every boundary line is also part of a
411 // boundary face.
412 for (const auto &cell : dof.active_cell_iterators())
413 for (const unsigned int f : cell->face_indices())
414 if (cell->at_boundary(f))
415 {
416 const unsigned int dofs_per_face =
417 cell->get_fe().n_dofs_per_face(f);
418 dofs_on_this_face.resize(dofs_per_face);
419 cell->face(f)->get_dof_indices(dofs_on_this_face,
420 cell->active_fe_index());
421
422 // make sparsity pattern for this cell
423 cols.clear();
424 for (const auto &dof : dofs_on_this_face)
425 cols.push_back(dof_to_boundary_mapping[dof]);
426 // We are not guaranteed that the mapping to a second index space
427 // is increasing so sort here to use the faster add_row_entries()
428 // path
429 std::sort(cols.begin(), cols.end());
430 for (const auto &dof : dofs_on_this_face)
431 sparsity.add_row_entries(dof_to_boundary_mapping[dof],
432 make_array_view(cols),
433 true);
434 }
435 }
436
437
438
439 template <int dim, int spacedim, typename number>
440 void
442 const DoFHandler<dim, spacedim> &dof,
443 const std::map<types::boundary_id, const Function<spacedim, number> *>
444 &boundary_ids,
445 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
446 SparsityPatternBase &sparsity)
447 {
448 if (dim == 1)
449 {
450 // first check left, then right boundary point
451 for (unsigned int direction = 0; direction < 2; ++direction)
452 {
453 // if this boundary is not requested, then go on with next one
454 if (boundary_ids.find(direction) == boundary_ids.end())
455 continue;
456
457 // find active cell at that boundary: first go to left/right,
458 // then to children
460 dof.begin(0);
461 while (!cell->at_boundary(direction))
462 cell = cell->neighbor(direction);
463 while (!cell->is_active())
464 cell = cell->child(direction);
465
466 const unsigned int dofs_per_vertex =
467 cell->get_fe().n_dofs_per_vertex();
468 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
469 dofs_per_vertex);
470
471 // next get boundary mapped dof indices of boundary dofs
472 for (unsigned int i = 0; i < dofs_per_vertex; ++i)
473 boundary_dof_boundary_indices[i] =
474 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
475
476 std::sort(boundary_dof_boundary_indices.begin(),
477 boundary_dof_boundary_indices.end());
478 for (const auto &dof : boundary_dof_boundary_indices)
479 sparsity.add_row_entries(
480 dof, make_array_view(boundary_dof_boundary_indices), true);
481 }
482 return;
483 }
484
485 const types::global_dof_index n_dofs = dof.n_dofs();
486 (void)n_dofs;
487
488 AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
489 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
490 boundary_ids.end(),
492
493 const bool fe_is_hermite = (dynamic_cast<const FE_Hermite<dim, spacedim> *>(
494 &(dof.get_fe())) != nullptr);
495
496 Assert(fe_is_hermite ||
497 sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
498 ExcDimensionMismatch(sparsity.n_rows(),
499 dof.n_boundary_dofs(boundary_ids)));
500 Assert(fe_is_hermite ||
501 sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
502 ExcDimensionMismatch(sparsity.n_cols(),
503 dof.n_boundary_dofs(boundary_ids)));
504 (void)fe_is_hermite;
505
506 if constexpr (running_in_debug_mode())
507 {
508 if (sparsity.n_rows() != 0)
509 {
510 types::global_dof_index max_element = 0;
511 for (const types::global_dof_index index : dof_to_boundary_mapping)
512 if ((index != numbers::invalid_dof_index) &&
513 (index > max_element))
514 max_element = index;
515 AssertDimension(max_element, sparsity.n_rows() - 1);
516 }
517 }
518
519 std::vector<types::global_dof_index> dofs_on_this_face;
520 dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
521 std::vector<types::global_dof_index> cols;
522
523 for (const auto &cell : dof.active_cell_iterators())
524 for (const unsigned int f : cell->face_indices())
525 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
526 boundary_ids.end())
527 {
528 const unsigned int dofs_per_face =
529 cell->get_fe().n_dofs_per_face(f);
530 dofs_on_this_face.resize(dofs_per_face);
531 cell->face(f)->get_dof_indices(dofs_on_this_face,
532 cell->active_fe_index());
533
534 // make sparsity pattern for this cell
535 cols.clear();
536 for (const auto &dof : dofs_on_this_face)
537 cols.push_back(dof_to_boundary_mapping[dof]);
538
539 // Like the other one: sort once.
540 std::sort(cols.begin(), cols.end());
541 for (const auto &dof : dofs_on_this_face)
542 sparsity.add_row_entries(dof_to_boundary_mapping[dof],
543 make_array_view(cols),
544 true);
545 }
546 }
547
548
549
550 template <int dim, int spacedim, typename number>
551 void
553 SparsityPatternBase &sparsity,
554 const AffineConstraints<number> &constraints,
555 const bool keep_constrained_dofs,
556 const types::subdomain_id subdomain_id)
557
558 // TODO: QA: reduce the indentation level of this method..., Maier 2012
559
560 {
561 const types::global_dof_index n_dofs = dof.n_dofs();
562 (void)n_dofs;
563
564 AssertDimension(sparsity.n_rows(), n_dofs);
565 AssertDimension(sparsity.n_cols(), n_dofs);
566
567 // If we have a distributed Triangulation only allow locally_owned
568 // subdomain. Not setting a subdomain is also okay, because we skip
569 // ghost cells in the loop below.
570 if (const auto *triangulation = dynamic_cast<
572 &dof.get_triangulation()))
573 {
574 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
575 (subdomain_id == triangulation->locally_owned_subdomain()),
577 "For distributed Triangulation objects and associated "
578 "DoFHandler objects, asking for any subdomain other than the "
579 "locally owned one does not make sense."));
580 }
581
582 std::vector<types::global_dof_index> dofs_on_this_cell;
583 std::vector<types::global_dof_index> dofs_on_other_cell;
584 dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
585 dofs_on_other_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
586
587 // TODO: in an old implementation, we used user flags before to tag
588 // faces that were already touched. this way, we could reduce the work
589 // a little bit. now, we instead add only data from one side. this
590 // should be OK, but we need to actually verify it.
591
592 // In case we work with a distributed sparsity pattern of Trilinos
593 // type, we only have to do the work if the current cell is owned by
594 // the calling processor. Otherwise, just continue.
595 for (const auto &cell : dof.active_cell_iterators())
596 if (((subdomain_id == numbers::invalid_subdomain_id) ||
597 (subdomain_id == cell->subdomain_id())) &&
598 cell->is_locally_owned())
599 {
600 const unsigned int n_dofs_on_this_cell =
601 cell->get_fe().n_dofs_per_cell();
602 dofs_on_this_cell.resize(n_dofs_on_this_cell);
603 cell->get_dof_indices(dofs_on_this_cell);
604
605 // make sparsity pattern for this cell. if no constraints pattern
606 // was given, then the following call acts as if simply no
607 // constraints existed
608 constraints.add_entries_local_to_global(dofs_on_this_cell,
609 sparsity,
610 keep_constrained_dofs);
611
612 for (const unsigned int face : cell->face_indices())
613 {
615 cell->face(face);
616 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
617 if (!cell->at_boundary(face) || periodic_neighbor)
618 {
620 neighbor = cell->neighbor_or_periodic_neighbor(face);
621
622 // in 1d, we do not need to worry whether the neighbor
623 // might have children and then loop over those children.
624 // rather, we may as well go straight to the cell behind
625 // this particular cell's most terminal child
626 if (dim == 1)
627 while (neighbor->has_children())
628 neighbor = neighbor->child(face == 0 ? 1 : 0);
629
630 if (neighbor->has_children())
631 {
632 for (unsigned int sub_nr = 0;
633 sub_nr != cell_face->n_active_descendants();
634 ++sub_nr)
635 {
636 const typename DoFHandler<dim, spacedim>::
637 level_cell_iterator sub_neighbor =
638 periodic_neighbor ?
639 cell->periodic_neighbor_child_on_subface(
640 face, sub_nr) :
641 cell->neighbor_child_on_subface(face, sub_nr);
642
643 const unsigned int n_dofs_on_neighbor =
644 sub_neighbor->get_fe().n_dofs_per_cell();
645 dofs_on_other_cell.resize(n_dofs_on_neighbor);
646 sub_neighbor->get_dof_indices(dofs_on_other_cell);
647
648 constraints.add_entries_local_to_global(
649 dofs_on_this_cell,
650 dofs_on_other_cell,
651 sparsity,
652 keep_constrained_dofs);
653 constraints.add_entries_local_to_global(
654 dofs_on_other_cell,
655 dofs_on_this_cell,
656 sparsity,
657 keep_constrained_dofs);
658 // only need to add this when the neighbor is not
659 // owned by the current processor, otherwise we add
660 // the entries for the neighbor there
661 if (sub_neighbor->subdomain_id() !=
662 cell->subdomain_id())
663 constraints.add_entries_local_to_global(
664 dofs_on_other_cell,
665 sparsity,
666 keep_constrained_dofs);
667 }
668 }
669 else
670 {
671 // Refinement edges are taken care of by coarser
672 // cells
673 if ((!periodic_neighbor &&
674 cell->neighbor_is_coarser(face)) ||
675 (periodic_neighbor &&
676 cell->periodic_neighbor_is_coarser(face)))
677 if (neighbor->subdomain_id() == cell->subdomain_id())
678 continue;
679
680 const unsigned int n_dofs_on_neighbor =
681 neighbor->get_fe().n_dofs_per_cell();
682 dofs_on_other_cell.resize(n_dofs_on_neighbor);
683
684 neighbor->get_dof_indices(dofs_on_other_cell);
685
686 constraints.add_entries_local_to_global(
687 dofs_on_this_cell,
688 dofs_on_other_cell,
689 sparsity,
690 keep_constrained_dofs);
691
692 // only need to add these in case the neighbor cell
693 // is not locally owned - otherwise, we touch each
694 // face twice and hence put the indices the other way
695 // around
696 if (!cell->neighbor_or_periodic_neighbor(face)
697 ->is_active() ||
698 (neighbor->subdomain_id() != cell->subdomain_id()))
699 {
700 constraints.add_entries_local_to_global(
701 dofs_on_other_cell,
702 dofs_on_this_cell,
703 sparsity,
704 keep_constrained_dofs);
705 if (neighbor->subdomain_id() != cell->subdomain_id())
706 constraints.add_entries_local_to_global(
707 dofs_on_other_cell,
708 sparsity,
709 keep_constrained_dofs);
710 }
711 }
712 }
713 }
714 }
715 }
716
717
718
719 template <int dim, int spacedim>
720 void
727
728 template <int dim, int spacedim>
732 const Table<2, Coupling> &component_couplings)
733 {
734 Assert(component_couplings.n_rows() == fe.n_components(),
735 ExcDimensionMismatch(component_couplings.n_rows(),
736 fe.n_components()));
737 Assert(component_couplings.n_cols() == fe.n_components(),
738 ExcDimensionMismatch(component_couplings.n_cols(),
739 fe.n_components()));
740
741 const unsigned int n_dofs = fe.n_dofs_per_cell();
742
743 Table<2, Coupling> dof_couplings(n_dofs, n_dofs);
744
745 for (unsigned int i = 0; i < n_dofs; ++i)
746 {
747 const unsigned int ii =
748 (fe.is_primitive(i) ?
749 fe.system_to_component_index(i).first :
752
753 for (unsigned int j = 0; j < n_dofs; ++j)
754 {
755 const unsigned int jj =
756 (fe.is_primitive(j) ?
757 fe.system_to_component_index(j).first :
760
761 dof_couplings(i, j) = component_couplings(ii, jj);
762 }
763 }
764 return dof_couplings;
765 }
766
767
768
769 template <int dim, int spacedim>
770 std::vector<Table<2, Coupling>>
773 const Table<2, Coupling> &component_couplings)
774 {
775 std::vector<Table<2, Coupling>> return_value(fe.size());
776 for (unsigned int i = 0; i < fe.size(); ++i)
777 return_value[i] =
778 dof_couplings_from_component_couplings(fe[i], component_couplings);
779
780 return return_value;
781 }
782
783
784
785 namespace internal
786 {
787 namespace
788 {
789 // helper function
790 template <typename Iterator, typename Iterator2>
791 void
792 add_cell_entries(
793 const Iterator &cell,
794 const unsigned int face_no,
795 const Iterator2 &neighbor,
796 const unsigned int neighbor_face_no,
797 const Table<2, Coupling> &flux_mask,
798 const std::vector<types::global_dof_index> &dofs_on_this_cell,
799 std::vector<types::global_dof_index> &dofs_on_other_cell,
800 std::vector<std::pair<SparsityPatternBase::size_type,
801 SparsityPatternBase::size_type>> &cell_entries)
802 {
803 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
804 neighbor->get_dof_indices(dofs_on_other_cell);
805
806 // Keep expensive data structures in separate vectors for inner j loop
807 // in separate vectors
808 boost::container::small_vector<unsigned int, 64>
809 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
810 boost::container::small_vector<bool, 64> support_on_face_i(
811 neighbor->get_fe().n_dofs_per_cell());
812 boost::container::small_vector<bool, 64> support_on_face_e(
813 neighbor->get_fe().n_dofs_per_cell());
814 for (unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
815 {
816 component_indices_neighbor[j] =
817 (neighbor->get_fe().is_primitive(j) ?
818 neighbor->get_fe().system_to_component_index(j).first :
819 neighbor->get_fe()
820 .get_nonzero_components(j)
821 .first_selected_component());
822 support_on_face_i[j] =
823 neighbor->get_fe().has_support_on_face(j, face_no);
824 support_on_face_e[j] =
825 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
826 }
827
828 // For the parallel setting, must include also the diagonal
829 // neighbor-neighbor coupling, otherwise those get included on the
830 // other cell
831 for (int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
832 {
833 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
834 const auto &dofs_i =
835 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
836 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
837 {
838 const unsigned int ii =
839 (fe.is_primitive(i) ?
840 fe.system_to_component_index(i).first :
841 fe.get_nonzero_components(i).first_selected_component());
842
843 Assert(ii < fe.n_components(), ExcInternalError());
844 const bool i_non_zero_i =
845 fe.has_support_on_face(i,
846 (f == 0 ? face_no : neighbor_face_no));
847
848 for (unsigned int j = 0;
849 j < neighbor->get_fe().n_dofs_per_cell();
850 ++j)
851 {
852 const bool j_non_zero_e = support_on_face_e[j];
853 const unsigned int jj = component_indices_neighbor[j];
854
855 Assert(jj < neighbor->get_fe().n_components(),
857
858 if ((flux_mask(ii, jj) == always) ||
859 (flux_mask(ii, jj) == nonzero && i_non_zero_i &&
860 j_non_zero_e))
861 cell_entries.emplace_back(dofs_i[i],
862 dofs_on_other_cell[j]);
863 if ((flux_mask(jj, ii) == always) ||
864 (flux_mask(jj, ii) == nonzero && j_non_zero_e &&
865 i_non_zero_i))
866 cell_entries.emplace_back(dofs_on_other_cell[j],
867 dofs_i[i]);
868 }
869 }
870 }
871 }
872
873
874
875 // implementation of the same function in namespace DoFTools for
876 // non-hp-DoFHandlers
877 template <int dim, int spacedim, typename number>
878 void
880 const DoFHandler<dim, spacedim> &dof,
881 SparsityPatternBase &sparsity,
882 const AffineConstraints<number> &constraints,
883 const bool keep_constrained_dofs,
884 const Table<2, Coupling> &int_mask,
885 const Table<2, Coupling> &flux_mask,
886 const types::subdomain_id subdomain_id,
887 const std::function<
889 const unsigned int)> &face_has_flux_coupling)
890 {
891 std::vector<types::global_dof_index> rows;
892 std::vector<std::pair<SparsityPatternBase::size_type,
894 cell_entries;
895
896 const ::hp::FECollection<dim, spacedim> &fe =
897 dof.get_fe_collection();
898
899 std::vector<types::global_dof_index> dofs_on_this_cell(
901 std::vector<types::global_dof_index> dofs_on_other_cell(
903
904 const unsigned int n_components = fe.n_components();
905 AssertDimension(int_mask.size(0), n_components);
906 AssertDimension(int_mask.size(1), n_components);
907 AssertDimension(flux_mask.size(0), n_components);
908 AssertDimension(flux_mask.size(1), n_components);
909
910 // note that we also need to set the respective entries if flux_mask
911 // says so. this is necessary since we need to consider all degrees
912 // of freedom on a cell for interior faces.
913 Table<2, Coupling> int_and_flux_mask(n_components, n_components);
914 for (unsigned int c1 = 0; c1 < n_components; ++c1)
915 for (unsigned int c2 = 0; c2 < n_components; ++c2)
916 if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
917 int_and_flux_mask(c1, c2) = always;
918
919 // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
920 // to constraints.add_entries_local_to_global()
921 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
922 dof_couplings_from_component_couplings(fe, int_and_flux_mask);
923
924 // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
925 // can pass it to constraints.add_entries_local_to_global()
926 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
927 for (unsigned int f = 0; f < fe.size(); ++f)
928 {
929 bool_int_and_flux_dof_mask[f].reinit(
930 TableIndices<2>(fe[f].n_dofs_per_cell(),
931 fe[f].n_dofs_per_cell()));
932 bool_int_and_flux_dof_mask[f].fill(false);
933 for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
934 for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
935 if (int_and_flux_dof_mask[f](i, j) != none)
936 bool_int_and_flux_dof_mask[f](i, j) = true;
937 }
938
939
940 for (const auto &cell : dof.active_cell_iterators())
942 (subdomain_id == cell->subdomain_id())) &&
943 cell->is_locally_owned())
944 {
945 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
946 cell->get_dof_indices(dofs_on_this_cell);
947
948 // make sparsity pattern for this cell also taking into
949 // account the couplings due to face contributions on the same
950 // cell
951 constraints.add_entries_local_to_global(
952 dofs_on_this_cell,
953 sparsity,
954 keep_constrained_dofs,
955 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
956
957 // Loop over interior faces
958 for (const unsigned int face : cell->face_indices())
959 {
960 const bool periodic_neighbor =
961 cell->has_periodic_neighbor(face);
962
963 if ((!cell->at_boundary(face)) || periodic_neighbor)
964 {
966 neighbor = cell->neighbor_or_periodic_neighbor(face);
967
968 // If the cells are on the same level (and both are
969 // active, locally-owned cells) then only add to the
970 // sparsity pattern if the current cell is 'greater' in
971 // the total ordering.
972 if (neighbor->level() == cell->level() &&
973 neighbor->index() > cell->index() &&
974 neighbor->is_active() && neighbor->is_locally_owned())
975 continue;
976
977 // If we are more refined then the neighbor, then we
978 // will automatically find the active neighbor cell when
979 // we call 'neighbor (face)' above. The opposite is not
980 // true; if the neighbor is more refined then the call
981 // 'neighbor (face)' will *not* return an active
982 // cell. Hence, only add things to the sparsity pattern
983 // if (when the levels are different) the neighbor is
984 // coarser than the current cell, except in the case
985 // when the neighbor is not locally owned.
986 if (neighbor->level() != cell->level() &&
987 ((!periodic_neighbor &&
988 !cell->neighbor_is_coarser(face)) ||
989 (periodic_neighbor &&
990 !cell->periodic_neighbor_is_coarser(face))) &&
991 neighbor->is_locally_owned())
992 continue; // (the neighbor is finer)
993
994 if (!face_has_flux_coupling(cell, face))
995 continue;
996
997 const unsigned int neighbor_face_no =
998 periodic_neighbor ?
999 cell->periodic_neighbor_face_no(face) :
1000 cell->neighbor_face_no(face);
1001
1002 // In 1d, go straight to the cell behind this
1003 // particular cell's most terminal cell. This makes us
1004 // skip the if (neighbor->has_children()) section
1005 // below. We need to do this since we otherwise
1006 // iterate over the children of the face, which are
1007 // always 0 in 1d.
1008 if (dim == 1)
1009 while (neighbor->has_children())
1010 neighbor = neighbor->child(face == 0 ? 1 : 0);
1011
1012 if (neighbor->has_children())
1013 {
1014 for (unsigned int sub_nr = 0;
1015 sub_nr != cell->face(face)->n_children();
1016 ++sub_nr)
1017 {
1018 const typename DoFHandler<dim, spacedim>::
1019 level_cell_iterator sub_neighbor =
1020 periodic_neighbor ?
1021 cell->periodic_neighbor_child_on_subface(
1022 face, sub_nr) :
1023 cell->neighbor_child_on_subface(face,
1024 sub_nr);
1025 add_cell_entries(cell,
1026 face,
1027 sub_neighbor,
1028 neighbor_face_no,
1029 flux_mask,
1030 dofs_on_this_cell,
1031 dofs_on_other_cell,
1032 cell_entries);
1033 }
1034 }
1035 else
1036 add_cell_entries(cell,
1037 face,
1038 neighbor,
1039 neighbor_face_no,
1040 flux_mask,
1041 dofs_on_this_cell,
1042 dofs_on_other_cell,
1043 cell_entries);
1044 }
1045 }
1046 sparsity.add_entries(make_array_view(cell_entries));
1047 cell_entries.clear();
1048 }
1049 }
1050 } // namespace
1051
1052 } // namespace internal
1053
1054
1055
1056 template <int dim, int spacedim>
1057 void
1059 SparsityPatternBase &sparsity,
1060 const Table<2, Coupling> &int_mask,
1061 const Table<2, Coupling> &flux_mask,
1062 const types::subdomain_id subdomain_id)
1063 {
1065
1066 const bool keep_constrained_dofs = true;
1067
1069 sparsity,
1070 dummy,
1071 keep_constrained_dofs,
1072 int_mask,
1073 flux_mask,
1074 subdomain_id,
1075 internal::always_couple_on_faces<dim, spacedim>);
1076 }
1077
1078
1079
1080 template <int dim, int spacedim, typename number>
1081 void
1083 const DoFHandler<dim, spacedim> &dof,
1084 SparsityPatternBase &sparsity,
1085 const AffineConstraints<number> &constraints,
1086 const bool keep_constrained_dofs,
1087 const Table<2, Coupling> &int_mask,
1088 const Table<2, Coupling> &flux_mask,
1089 const types::subdomain_id subdomain_id,
1090 const std::function<
1092 const unsigned int)> &face_has_flux_coupling)
1093 {
1094 // do the error checking and frame code here, and then pass on to more
1095 // specialized functions in the internal namespace
1096 const types::global_dof_index n_dofs = dof.n_dofs();
1097 (void)n_dofs;
1098 const unsigned int n_comp = dof.get_fe(0).n_components();
1099 (void)n_comp;
1100
1101 Assert(sparsity.n_rows() == n_dofs,
1102 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1103 Assert(sparsity.n_cols() == n_dofs,
1104 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1105 Assert(int_mask.n_rows() == n_comp,
1106 ExcDimensionMismatch(int_mask.n_rows(), n_comp));
1107 Assert(int_mask.n_cols() == n_comp,
1108 ExcDimensionMismatch(int_mask.n_cols(), n_comp));
1109 Assert(flux_mask.n_rows() == n_comp,
1110 ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
1111 Assert(flux_mask.n_cols() == n_comp,
1112 ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
1113
1114 // If we have a distributed Triangulation only allow locally_owned
1115 // subdomain. Not setting a subdomain is also okay, because we skip
1116 // ghost cells in the loop below.
1117 if (const auto *triangulation = dynamic_cast<
1119 &dof.get_triangulation()))
1120 {
1121 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
1122 (subdomain_id == triangulation->locally_owned_subdomain()),
1123 ExcMessage(
1124 "For distributed Triangulation objects and associated "
1125 "DoFHandler objects, asking for any subdomain other than the "
1126 "locally owned one does not make sense."));
1127 }
1128
1129 Assert(
1130 face_has_flux_coupling,
1131 ExcMessage(
1132 "The function which specifies if a flux coupling occurs over a given "
1133 "face is empty."));
1134
1135 internal::make_flux_sparsity_pattern(dof,
1136 sparsity,
1137 constraints,
1138 keep_constrained_dofs,
1139 int_mask,
1140 flux_mask,
1141 subdomain_id,
1142 face_has_flux_coupling);
1143 }
1144
1145} // end of namespace DoFTools
1146
1147
1148// --------------------------------------------------- explicit instantiations
1149
1150#include "dofs/dof_tools_sparsity.inst"
1151
1152
1153
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
typename LevelSelector::cell_iterator level_cell_iterator
cell_iterator begin(const unsigned int level=0) const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
types::global_dof_index size_type
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false)=0
size_type n_cols() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int size() const
Definition collection.h:314
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
constexpr bool running_in_debug_mode()
Definition config.h:76
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
void make_boundary_sparsity_pattern(const DoFHandler< dim, spacedim > &dof, const std::vector< types::global_dof_index > &dof_to_boundary_mapping, SparsityPatternBase &sparsity_pattern)
Table< 2, Coupling > dof_couplings_from_component_couplings(const FiniteElement< dim, spacedim > &fe, const Table< 2, Coupling > &component_couplings)
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
*  *  if(update_pressure &update_flags) *  compute_pressure(constitutive_request
constexpr types::global_dof_index invalid_dof_index
Definition types.h:259
constexpr types::boundary_id internal_face_boundary_id
Definition types.h:319
constexpr types::subdomain_id invalid_subdomain_id
Definition types.h:385
unsigned int subdomain_id
Definition types.h:50
unsigned short int fe_index
Definition types.h:70