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