deal.II version GIT relicensing-2289-g1e5549a87a 2024-12-21 21:30: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#ifdef DEBUG
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) && (index > max_element))
397 max_element = index;
398 AssertDimension(max_element, sparsity.n_rows() - 1);
399 }
400#endif
401
402 std::vector<types::global_dof_index> dofs_on_this_face;
403 dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
404 std::vector<types::global_dof_index> cols;
405
406 // loop over all faces to check whether they are at a boundary. note
407 // that we need not take special care of single lines (using
408 // @p{cell->has_boundary_lines}), since we do not support boundaries of
409 // dimension dim-2, and so every boundary line is also part of a
410 // boundary face.
411 for (const auto &cell : dof.active_cell_iterators())
412 for (const unsigned int f : cell->face_indices())
413 if (cell->at_boundary(f))
414 {
415 const unsigned int dofs_per_face =
416 cell->get_fe().n_dofs_per_face(f);
417 dofs_on_this_face.resize(dofs_per_face);
418 cell->face(f)->get_dof_indices(dofs_on_this_face,
419 cell->active_fe_index());
420
421 // make sparsity pattern for this cell
422 cols.clear();
423 for (const auto &dof : dofs_on_this_face)
424 cols.push_back(dof_to_boundary_mapping[dof]);
425 // We are not guaranteed that the mapping to a second index space
426 // is increasing so sort here to use the faster add_row_entries()
427 // path
428 std::sort(cols.begin(), cols.end());
429 for (const auto &dof : dofs_on_this_face)
430 sparsity.add_row_entries(dof_to_boundary_mapping[dof],
431 make_array_view(cols),
432 true);
433 }
434 }
435
436
437
438 template <int dim, int spacedim, typename number>
439 void
441 const DoFHandler<dim, spacedim> &dof,
442 const std::map<types::boundary_id, const Function<spacedim, number> *>
443 &boundary_ids,
444 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
445 SparsityPatternBase &sparsity)
446 {
447 if (dim == 1)
448 {
449 // first check left, then right boundary point
450 for (unsigned int direction = 0; direction < 2; ++direction)
451 {
452 // if this boundary is not requested, then go on with next one
453 if (boundary_ids.find(direction) == boundary_ids.end())
454 continue;
455
456 // find active cell at that boundary: first go to left/right,
457 // then to children
459 dof.begin(0);
460 while (!cell->at_boundary(direction))
461 cell = cell->neighbor(direction);
462 while (!cell->is_active())
463 cell = cell->child(direction);
464
465 const unsigned int dofs_per_vertex =
466 cell->get_fe().n_dofs_per_vertex();
467 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
468 dofs_per_vertex);
469
470 // next get boundary mapped dof indices of boundary dofs
471 for (unsigned int i = 0; i < dofs_per_vertex; ++i)
472 boundary_dof_boundary_indices[i] =
473 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
474
475 std::sort(boundary_dof_boundary_indices.begin(),
476 boundary_dof_boundary_indices.end());
477 for (const auto &dof : boundary_dof_boundary_indices)
478 sparsity.add_row_entries(
479 dof, make_array_view(boundary_dof_boundary_indices), true);
480 }
481 return;
482 }
483
484 const types::global_dof_index n_dofs = dof.n_dofs();
485 (void)n_dofs;
486
487 AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
488 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
489 boundary_ids.end(),
491
492 const bool fe_is_hermite = (dynamic_cast<const FE_Hermite<dim, spacedim> *>(
493 &(dof.get_fe())) != nullptr);
494
495 Assert(fe_is_hermite ||
496 sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
497 ExcDimensionMismatch(sparsity.n_rows(),
498 dof.n_boundary_dofs(boundary_ids)));
499 Assert(fe_is_hermite ||
500 sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
501 ExcDimensionMismatch(sparsity.n_cols(),
502 dof.n_boundary_dofs(boundary_ids)));
503 (void)fe_is_hermite;
504
505#ifdef DEBUG
506 if (sparsity.n_rows() != 0)
507 {
508 types::global_dof_index max_element = 0;
509 for (const types::global_dof_index index : dof_to_boundary_mapping)
510 if ((index != numbers::invalid_dof_index) && (index > max_element))
511 max_element = index;
512 AssertDimension(max_element, sparsity.n_rows() - 1);
513 }
514#endif
515
516 std::vector<types::global_dof_index> dofs_on_this_face;
517 dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
518 std::vector<types::global_dof_index> cols;
519
520 for (const auto &cell : dof.active_cell_iterators())
521 for (const unsigned int f : cell->face_indices())
522 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
523 boundary_ids.end())
524 {
525 const unsigned int dofs_per_face =
526 cell->get_fe().n_dofs_per_face(f);
527 dofs_on_this_face.resize(dofs_per_face);
528 cell->face(f)->get_dof_indices(dofs_on_this_face,
529 cell->active_fe_index());
530
531 // make sparsity pattern for this cell
532 cols.clear();
533 for (const auto &dof : dofs_on_this_face)
534 cols.push_back(dof_to_boundary_mapping[dof]);
535
536 // Like the other one: sort once.
537 std::sort(cols.begin(), cols.end());
538 for (const auto &dof : dofs_on_this_face)
539 sparsity.add_row_entries(dof_to_boundary_mapping[dof],
540 make_array_view(cols),
541 true);
542 }
543 }
544
545
546
547 template <int dim, int spacedim, typename number>
548 void
550 SparsityPatternBase &sparsity,
551 const AffineConstraints<number> &constraints,
552 const bool keep_constrained_dofs,
553 const types::subdomain_id subdomain_id)
554
555 // TODO: QA: reduce the indentation level of this method..., Maier 2012
556
557 {
558 const types::global_dof_index n_dofs = dof.n_dofs();
559 (void)n_dofs;
560
561 AssertDimension(sparsity.n_rows(), n_dofs);
562 AssertDimension(sparsity.n_cols(), n_dofs);
563
564 // If we have a distributed Triangulation only allow locally_owned
565 // subdomain. Not setting a subdomain is also okay, because we skip
566 // ghost cells in the loop below.
567 if (const auto *triangulation = dynamic_cast<
569 &dof.get_triangulation()))
570 {
571 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
572 (subdomain_id == triangulation->locally_owned_subdomain()),
574 "For distributed Triangulation objects and associated "
575 "DoFHandler objects, asking for any subdomain other than the "
576 "locally owned one does not make sense."));
577 }
578
579 std::vector<types::global_dof_index> dofs_on_this_cell;
580 std::vector<types::global_dof_index> dofs_on_other_cell;
581 dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
582 dofs_on_other_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
583
584 // TODO: in an old implementation, we used user flags before to tag
585 // faces that were already touched. this way, we could reduce the work
586 // a little bit. now, we instead add only data from one side. this
587 // should be OK, but we need to actually verify it.
588
589 // In case we work with a distributed sparsity pattern of Trilinos
590 // type, we only have to do the work if the current cell is owned by
591 // the calling processor. Otherwise, just continue.
592 for (const auto &cell : dof.active_cell_iterators())
593 if (((subdomain_id == numbers::invalid_subdomain_id) ||
594 (subdomain_id == cell->subdomain_id())) &&
595 cell->is_locally_owned())
596 {
597 const unsigned int n_dofs_on_this_cell =
598 cell->get_fe().n_dofs_per_cell();
599 dofs_on_this_cell.resize(n_dofs_on_this_cell);
600 cell->get_dof_indices(dofs_on_this_cell);
601
602 // make sparsity pattern for this cell. if no constraints pattern
603 // was given, then the following call acts as if simply no
604 // constraints existed
605 constraints.add_entries_local_to_global(dofs_on_this_cell,
606 sparsity,
607 keep_constrained_dofs);
608
609 for (const unsigned int face : cell->face_indices())
610 {
612 cell->face(face);
613 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
614 if (!cell->at_boundary(face) || periodic_neighbor)
615 {
617 neighbor = cell->neighbor_or_periodic_neighbor(face);
618
619 // in 1d, we do not need to worry whether the neighbor
620 // might have children and then loop over those children.
621 // rather, we may as well go straight to the cell behind
622 // this particular cell's most terminal child
623 if (dim == 1)
624 while (neighbor->has_children())
625 neighbor = neighbor->child(face == 0 ? 1 : 0);
626
627 if (neighbor->has_children())
628 {
629 for (unsigned int sub_nr = 0;
630 sub_nr != cell_face->n_active_descendants();
631 ++sub_nr)
632 {
633 const typename DoFHandler<dim, spacedim>::
634 level_cell_iterator sub_neighbor =
635 periodic_neighbor ?
636 cell->periodic_neighbor_child_on_subface(
637 face, sub_nr) :
638 cell->neighbor_child_on_subface(face, sub_nr);
639
640 const unsigned int n_dofs_on_neighbor =
641 sub_neighbor->get_fe().n_dofs_per_cell();
642 dofs_on_other_cell.resize(n_dofs_on_neighbor);
643 sub_neighbor->get_dof_indices(dofs_on_other_cell);
644
645 constraints.add_entries_local_to_global(
646 dofs_on_this_cell,
647 dofs_on_other_cell,
648 sparsity,
649 keep_constrained_dofs);
650 constraints.add_entries_local_to_global(
651 dofs_on_other_cell,
652 dofs_on_this_cell,
653 sparsity,
654 keep_constrained_dofs);
655 // only need to add this when the neighbor is not
656 // owned by the current processor, otherwise we add
657 // the entries for the neighbor there
658 if (sub_neighbor->subdomain_id() !=
659 cell->subdomain_id())
660 constraints.add_entries_local_to_global(
661 dofs_on_other_cell,
662 sparsity,
663 keep_constrained_dofs);
664 }
665 }
666 else
667 {
668 // Refinement edges are taken care of by coarser
669 // cells
670 if ((!periodic_neighbor &&
671 cell->neighbor_is_coarser(face)) ||
672 (periodic_neighbor &&
673 cell->periodic_neighbor_is_coarser(face)))
674 if (neighbor->subdomain_id() == cell->subdomain_id())
675 continue;
676
677 const unsigned int n_dofs_on_neighbor =
678 neighbor->get_fe().n_dofs_per_cell();
679 dofs_on_other_cell.resize(n_dofs_on_neighbor);
680
681 neighbor->get_dof_indices(dofs_on_other_cell);
682
683 constraints.add_entries_local_to_global(
684 dofs_on_this_cell,
685 dofs_on_other_cell,
686 sparsity,
687 keep_constrained_dofs);
688
689 // only need to add these in case the neighbor cell
690 // is not locally owned - otherwise, we touch each
691 // face twice and hence put the indices the other way
692 // around
693 if (!cell->neighbor_or_periodic_neighbor(face)
694 ->is_active() ||
695 (neighbor->subdomain_id() != cell->subdomain_id()))
696 {
697 constraints.add_entries_local_to_global(
698 dofs_on_other_cell,
699 dofs_on_this_cell,
700 sparsity,
701 keep_constrained_dofs);
702 if (neighbor->subdomain_id() != cell->subdomain_id())
703 constraints.add_entries_local_to_global(
704 dofs_on_other_cell,
705 sparsity,
706 keep_constrained_dofs);
707 }
708 }
709 }
710 }
711 }
712 }
713
714
715
716 template <int dim, int spacedim>
717 void
724
725 template <int dim, int spacedim>
729 const Table<2, Coupling> &component_couplings)
730 {
731 Assert(component_couplings.n_rows() == fe.n_components(),
732 ExcDimensionMismatch(component_couplings.n_rows(),
733 fe.n_components()));
734 Assert(component_couplings.n_cols() == fe.n_components(),
735 ExcDimensionMismatch(component_couplings.n_cols(),
736 fe.n_components()));
737
738 const unsigned int n_dofs = fe.n_dofs_per_cell();
739
740 Table<2, Coupling> dof_couplings(n_dofs, n_dofs);
741
742 for (unsigned int i = 0; i < n_dofs; ++i)
743 {
744 const unsigned int ii =
745 (fe.is_primitive(i) ?
746 fe.system_to_component_index(i).first :
749
750 for (unsigned int j = 0; j < n_dofs; ++j)
751 {
752 const unsigned int jj =
753 (fe.is_primitive(j) ?
754 fe.system_to_component_index(j).first :
757
758 dof_couplings(i, j) = component_couplings(ii, jj);
759 }
760 }
761 return dof_couplings;
762 }
763
764
765
766 template <int dim, int spacedim>
767 std::vector<Table<2, Coupling>>
770 const Table<2, Coupling> &component_couplings)
771 {
772 std::vector<Table<2, Coupling>> return_value(fe.size());
773 for (unsigned int i = 0; i < fe.size(); ++i)
774 return_value[i] =
775 dof_couplings_from_component_couplings(fe[i], component_couplings);
776
777 return return_value;
778 }
779
780
781
782 namespace internal
783 {
784 namespace
785 {
786 // helper function
787 template <typename Iterator, typename Iterator2>
788 void
789 add_cell_entries(
790 const Iterator &cell,
791 const unsigned int face_no,
792 const Iterator2 &neighbor,
793 const unsigned int neighbor_face_no,
794 const Table<2, Coupling> &flux_mask,
795 const std::vector<types::global_dof_index> &dofs_on_this_cell,
796 std::vector<types::global_dof_index> &dofs_on_other_cell,
797 std::vector<std::pair<SparsityPatternBase::size_type,
798 SparsityPatternBase::size_type>> &cell_entries)
799 {
800 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
801 neighbor->get_dof_indices(dofs_on_other_cell);
802
803 // Keep expensive data structures in separate vectors for inner j loop
804 // in separate vectors
805 boost::container::small_vector<unsigned int, 64>
806 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
807 boost::container::small_vector<bool, 64> support_on_face_i(
808 neighbor->get_fe().n_dofs_per_cell());
809 boost::container::small_vector<bool, 64> support_on_face_e(
810 neighbor->get_fe().n_dofs_per_cell());
811 for (unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
812 {
813 component_indices_neighbor[j] =
814 (neighbor->get_fe().is_primitive(j) ?
815 neighbor->get_fe().system_to_component_index(j).first :
816 neighbor->get_fe()
817 .get_nonzero_components(j)
818 .first_selected_component());
819 support_on_face_i[j] =
820 neighbor->get_fe().has_support_on_face(j, face_no);
821 support_on_face_e[j] =
822 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
823 }
824
825 // For the parallel setting, must include also the diagonal
826 // neighbor-neighbor coupling, otherwise those get included on the
827 // other cell
828 for (int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
829 {
830 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
831 const auto &dofs_i =
832 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
833 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
834 {
835 const unsigned int ii =
836 (fe.is_primitive(i) ?
837 fe.system_to_component_index(i).first :
838 fe.get_nonzero_components(i).first_selected_component());
839
840 Assert(ii < fe.n_components(), ExcInternalError());
841 const bool i_non_zero_i =
842 fe.has_support_on_face(i,
843 (f == 0 ? face_no : neighbor_face_no));
844
845 for (unsigned int j = 0;
846 j < neighbor->get_fe().n_dofs_per_cell();
847 ++j)
848 {
849 const bool j_non_zero_e = support_on_face_e[j];
850 const unsigned int jj = component_indices_neighbor[j];
851
852 Assert(jj < neighbor->get_fe().n_components(),
854
855 if ((flux_mask(ii, jj) == always) ||
856 (flux_mask(ii, jj) == nonzero && i_non_zero_i &&
857 j_non_zero_e))
858 cell_entries.emplace_back(dofs_i[i],
859 dofs_on_other_cell[j]);
860 if ((flux_mask(jj, ii) == always) ||
861 (flux_mask(jj, ii) == nonzero && j_non_zero_e &&
862 i_non_zero_i))
863 cell_entries.emplace_back(dofs_on_other_cell[j],
864 dofs_i[i]);
865 }
866 }
867 }
868 }
869
870
871
872 // implementation of the same function in namespace DoFTools for
873 // non-hp-DoFHandlers
874 template <int dim, int spacedim, typename number>
875 void
877 const DoFHandler<dim, spacedim> &dof,
878 SparsityPatternBase &sparsity,
879 const AffineConstraints<number> &constraints,
880 const bool keep_constrained_dofs,
881 const Table<2, Coupling> &int_mask,
882 const Table<2, Coupling> &flux_mask,
883 const types::subdomain_id subdomain_id,
884 const std::function<
886 const unsigned int)> &face_has_flux_coupling)
887 {
888 std::vector<types::global_dof_index> rows;
889 std::vector<std::pair<SparsityPatternBase::size_type,
891 cell_entries;
892
893 const ::hp::FECollection<dim, spacedim> &fe =
894 dof.get_fe_collection();
895
896 std::vector<types::global_dof_index> dofs_on_this_cell(
898 std::vector<types::global_dof_index> dofs_on_other_cell(
900
901 const unsigned int n_components = fe.n_components();
902 AssertDimension(int_mask.size(0), n_components);
903 AssertDimension(int_mask.size(1), n_components);
904 AssertDimension(flux_mask.size(0), n_components);
905 AssertDimension(flux_mask.size(1), n_components);
906
907 // note that we also need to set the respective entries if flux_mask
908 // says so. this is necessary since we need to consider all degrees
909 // of freedom on a cell for interior faces.
910 Table<2, Coupling> int_and_flux_mask(n_components, n_components);
911 for (unsigned int c1 = 0; c1 < n_components; ++c1)
912 for (unsigned int c2 = 0; c2 < n_components; ++c2)
913 if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
914 int_and_flux_mask(c1, c2) = always;
915
916 // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
917 // to constraints.add_entries_local_to_global()
918 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
919 dof_couplings_from_component_couplings(fe, int_and_flux_mask);
920
921 // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
922 // can pass it to constraints.add_entries_local_to_global()
923 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
924 for (unsigned int f = 0; f < fe.size(); ++f)
925 {
926 bool_int_and_flux_dof_mask[f].reinit(
927 TableIndices<2>(fe[f].n_dofs_per_cell(),
928 fe[f].n_dofs_per_cell()));
929 bool_int_and_flux_dof_mask[f].fill(false);
930 for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
931 for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
932 if (int_and_flux_dof_mask[f](i, j) != none)
933 bool_int_and_flux_dof_mask[f](i, j) = true;
934 }
935
936
937 for (const auto &cell : dof.active_cell_iterators())
939 (subdomain_id == cell->subdomain_id())) &&
940 cell->is_locally_owned())
941 {
942 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
943 cell->get_dof_indices(dofs_on_this_cell);
944
945 // make sparsity pattern for this cell also taking into
946 // account the couplings due to face contributions on the same
947 // cell
948 constraints.add_entries_local_to_global(
949 dofs_on_this_cell,
950 sparsity,
951 keep_constrained_dofs,
952 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
953
954 // Loop over interior faces
955 for (const unsigned int face : cell->face_indices())
956 {
957 const bool periodic_neighbor =
958 cell->has_periodic_neighbor(face);
959
960 if ((!cell->at_boundary(face)) || periodic_neighbor)
961 {
963 neighbor = cell->neighbor_or_periodic_neighbor(face);
964
965 // If the cells are on the same level (and both are
966 // active, locally-owned cells) then only add to the
967 // sparsity pattern if the current cell is 'greater' in
968 // the total ordering.
969 if (neighbor->level() == cell->level() &&
970 neighbor->index() > cell->index() &&
971 neighbor->is_active() && neighbor->is_locally_owned())
972 continue;
973
974 // If we are more refined then the neighbor, then we
975 // will automatically find the active neighbor cell when
976 // we call 'neighbor (face)' above. The opposite is not
977 // true; if the neighbor is more refined then the call
978 // 'neighbor (face)' will *not* return an active
979 // cell. Hence, only add things to the sparsity pattern
980 // if (when the levels are different) the neighbor is
981 // coarser than the current cell, except in the case
982 // when the neighbor is not locally owned.
983 if (neighbor->level() != cell->level() &&
984 ((!periodic_neighbor &&
985 !cell->neighbor_is_coarser(face)) ||
986 (periodic_neighbor &&
987 !cell->periodic_neighbor_is_coarser(face))) &&
988 neighbor->is_locally_owned())
989 continue; // (the neighbor is finer)
990
991 if (!face_has_flux_coupling(cell, face))
992 continue;
993
994 const unsigned int neighbor_face_no =
995 periodic_neighbor ?
996 cell->periodic_neighbor_face_no(face) :
997 cell->neighbor_face_no(face);
998
999 // In 1d, go straight to the cell behind this
1000 // particular cell's most terminal cell. This makes us
1001 // skip the if (neighbor->has_children()) section
1002 // below. We need to do this since we otherwise
1003 // iterate over the children of the face, which are
1004 // always 0 in 1d.
1005 if (dim == 1)
1006 while (neighbor->has_children())
1007 neighbor = neighbor->child(face == 0 ? 1 : 0);
1008
1009 if (neighbor->has_children())
1010 {
1011 for (unsigned int sub_nr = 0;
1012 sub_nr != cell->face(face)->n_children();
1013 ++sub_nr)
1014 {
1015 const typename DoFHandler<dim, spacedim>::
1016 level_cell_iterator sub_neighbor =
1017 periodic_neighbor ?
1018 cell->periodic_neighbor_child_on_subface(
1019 face, sub_nr) :
1020 cell->neighbor_child_on_subface(face,
1021 sub_nr);
1022 add_cell_entries(cell,
1023 face,
1024 sub_neighbor,
1025 neighbor_face_no,
1026 flux_mask,
1027 dofs_on_this_cell,
1028 dofs_on_other_cell,
1029 cell_entries);
1030 }
1031 }
1032 else
1033 add_cell_entries(cell,
1034 face,
1035 neighbor,
1036 neighbor_face_no,
1037 flux_mask,
1038 dofs_on_this_cell,
1039 dofs_on_other_cell,
1040 cell_entries);
1041 }
1042 }
1043 sparsity.add_entries(make_array_view(cell_entries));
1044 cell_entries.clear();
1045 }
1046 }
1047 } // namespace
1048
1049 } // namespace internal
1050
1051
1052
1053 template <int dim, int spacedim>
1054 void
1056 SparsityPatternBase &sparsity,
1057 const Table<2, Coupling> &int_mask,
1058 const Table<2, Coupling> &flux_mask,
1059 const types::subdomain_id subdomain_id)
1060 {
1062
1063 const bool keep_constrained_dofs = true;
1064
1066 sparsity,
1067 dummy,
1068 keep_constrained_dofs,
1069 int_mask,
1070 flux_mask,
1071 subdomain_id,
1072 internal::always_couple_on_faces<dim, spacedim>);
1073 }
1074
1075
1076
1077 template <int dim, int spacedim, typename number>
1078 void
1080 const DoFHandler<dim, spacedim> &dof,
1081 SparsityPatternBase &sparsity,
1082 const AffineConstraints<number> &constraints,
1083 const bool keep_constrained_dofs,
1084 const Table<2, Coupling> &int_mask,
1085 const Table<2, Coupling> &flux_mask,
1086 const types::subdomain_id subdomain_id,
1087 const std::function<
1089 const unsigned int)> &face_has_flux_coupling)
1090 {
1091 // do the error checking and frame code here, and then pass on to more
1092 // specialized functions in the internal namespace
1093 const types::global_dof_index n_dofs = dof.n_dofs();
1094 (void)n_dofs;
1095 const unsigned int n_comp = dof.get_fe(0).n_components();
1096 (void)n_comp;
1097
1098 Assert(sparsity.n_rows() == n_dofs,
1099 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1100 Assert(sparsity.n_cols() == n_dofs,
1101 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1102 Assert(int_mask.n_rows() == n_comp,
1103 ExcDimensionMismatch(int_mask.n_rows(), n_comp));
1104 Assert(int_mask.n_cols() == n_comp,
1105 ExcDimensionMismatch(int_mask.n_cols(), n_comp));
1106 Assert(flux_mask.n_rows() == n_comp,
1107 ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
1108 Assert(flux_mask.n_cols() == n_comp,
1109 ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
1110
1111 // If we have a distributed Triangulation only allow locally_owned
1112 // subdomain. Not setting a subdomain is also okay, because we skip
1113 // ghost cells in the loop below.
1114 if (const auto *triangulation = dynamic_cast<
1116 &dof.get_triangulation()))
1117 {
1118 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
1119 (subdomain_id == triangulation->locally_owned_subdomain()),
1120 ExcMessage(
1121 "For distributed Triangulation objects and associated "
1122 "DoFHandler objects, asking for any subdomain other than the "
1123 "locally owned one does not make sense."));
1124 }
1125
1126 Assert(
1127 face_has_flux_coupling,
1128 ExcMessage(
1129 "The function which specifies if a flux coupling occurs over a given "
1130 "face is empty."));
1131
1132 internal::make_flux_sparsity_pattern(dof,
1133 sparsity,
1134 constraints,
1135 keep_constrained_dofs,
1136 int_mask,
1137 flux_mask,
1138 subdomain_id,
1139 face_has_flux_coupling);
1140 }
1141
1142} // end of namespace DoFTools
1143
1144
1145// --------------------------------------------------- explicit instantiations
1146
1147#include "dof_tools_sparsity.inst"
1148
1149
1150
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:949
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:315
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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)
const types::boundary_id internal_face_boundary_id
Definition types.h:312
const types::subdomain_id invalid_subdomain_id
Definition types.h:341
const types::global_dof_index invalid_dof_index
Definition types.h:252
unsigned int subdomain_id
Definition types.h:43
unsigned short int fe_index
Definition types.h:59
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation