Reference documentation for deal.II version GIT relicensing-878-g25c11e0287 2024-06-21 14:50:02+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.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 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
15#include <deal.II/base/mpi.h>
17#include <deal.II/base/table.h>
19
22
26
27#include <deal.II/fe/fe.h>
28#include <deal.II/fe/fe_tools.h>
30
33#include <deal.II/grid/tria.h>
35
40
44#include <deal.II/lac/vector.h>
45
46#include <algorithm>
47#include <numeric>
48
50
51
52
53namespace DoFTools
54{
55 namespace internal
56 {
71 template <int dim, typename Number = double>
73 {
79 bool
81 const Point<dim, Number> &rhs) const
82 {
83 double downstream_size = 0;
84 double weight = 1.;
85 for (unsigned int d = 0; d < dim; ++d)
86 {
87 downstream_size += (rhs[d] - lhs[d]) * weight;
88 weight *= 1e-5;
89 }
90 if (downstream_size < 0)
91 return false;
92 else if (downstream_size > 0)
93 return true;
94 else
95 {
96 for (unsigned int d = 0; d < dim; ++d)
97 {
98 if (lhs[d] == rhs[d])
99 continue;
100 return lhs[d] < rhs[d];
101 }
102 return false;
103 }
104 }
105 };
106
107
108
109 // return an array that for each dof on the reference cell
110 // lists the corresponding vector component.
111 //
112 // if an element is non-primitive then we assign to each degree of freedom
113 // the following component:
114 // - if the nonzero components that belong to a shape function are not
115 // selected in the component_mask, then the shape function is assigned
116 // to the first nonzero vector component that corresponds to this
117 // shape function
118 // - otherwise, the shape function is assigned the first component selected
119 // in the component_mask that corresponds to this shape function
120 template <int dim, int spacedim>
121 std::vector<unsigned char>
123 const ComponentMask &component_mask)
124 {
125 std::vector<unsigned char> local_component_association(
126 fe.n_dofs_per_cell(), static_cast<unsigned char>(-1));
127
128 // compute the component each local dof belongs to.
129 // if the shape function is primitive, then this
130 // is simple and we can just associate it with
131 // what system_to_component_index gives us
132 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
133 if (fe.is_primitive(i))
134 local_component_association[i] =
135 fe.system_to_component_index(i).first;
136 else
137 // if the shape function is not primitive, then either use the
138 // component of the first nonzero component corresponding
139 // to this shape function (if the component is not specified
140 // in the component_mask), or the first component of this block
141 // that is listed in the component_mask (if the block this
142 // component corresponds to is indeed specified in the component
143 // mask)
144 {
145 const unsigned int first_comp =
147
148 if ((fe.get_nonzero_components(i) & component_mask)
149 .n_selected_components(fe.n_components()) == 0)
150 local_component_association[i] = first_comp;
151 else
152 // pick the component selected. we know from the previous 'if'
153 // that within the components that are nonzero for this
154 // shape function there must be at least one for which the
155 // mask is true, so we will for sure run into the break()
156 // at one point
157 for (unsigned int c = first_comp; c < fe.n_components(); ++c)
158 if (component_mask[c] == true)
159 {
160 local_component_association[i] = c;
161 break;
162 }
163 }
164
165 Assert(std::find(local_component_association.begin(),
166 local_component_association.end(),
167 static_cast<unsigned char>(-1)) ==
168 local_component_association.end(),
170
171 return local_component_association;
172 }
173
174
175 // this internal function assigns to each dof the respective component
176 // of the vector system.
177 //
178 // the output array dofs_by_component lists for each dof the
179 // corresponding vector component. if the DoFHandler is based on a
180 // parallel distributed triangulation then the output array is index by
181 // dof_handler.locally_owned_dofs().index_within_set(indices[i])
182 //
183 // if an element is non-primitive then we assign to each degree of
184 // freedom the following component:
185 // - if the nonzero components that belong to a shape function are not
186 // selected in the component_mask, then the shape function is assigned
187 // to the first nonzero vector component that corresponds to this
188 // shape function
189 // - otherwise, the shape function is assigned the first component selected
190 // in the component_mask that corresponds to this shape function
191 template <int dim, int spacedim>
192 void
194 const DoFHandler<dim, spacedim> &dof_handler,
195 const ComponentMask &component_mask,
196 std::vector<unsigned char> &dofs_by_component,
197 const unsigned int mg_level = numbers::invalid_unsigned_int)
198 {
199 const ::hp::FECollection<dim, spacedim> &fe_collection =
200 dof_handler.get_fe_collection();
201 Assert(fe_collection.n_components() < 256, ExcNotImplemented());
202
203 const auto &locally_owned_dofs =
204 (mg_level == numbers::invalid_unsigned_int) ?
205 dof_handler.locally_owned_dofs() :
206 dof_handler.locally_owned_mg_dofs(mg_level);
207
208 AssertDimension(dofs_by_component.size(),
209 locally_owned_dofs.n_elements());
210
211 // next set up a table for the degrees of freedom on each of the
212 // cells (regardless of the fact whether it is listed in the
213 // component_select argument or not)
214 //
215 // for each element 'f' of the FECollection,
216 // local_component_association[f][d] then returns the vector
217 // component that degree of freedom 'd' belongs to
218 std::vector<std::vector<unsigned char>> local_component_association(
219 fe_collection.size());
220 for (unsigned int f = 0; f < fe_collection.size(); ++f)
221 {
222 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
223 local_component_association[f] =
224 get_local_component_association(fe, component_mask);
225 }
226
227 // then loop over all cells and do the work
228 std::vector<types::global_dof_index> indices;
229
230 const auto runner = [&](const auto &task) {
231 if (mg_level == numbers::invalid_unsigned_int)
232 {
233 for (const auto &cell : dof_handler.active_cell_iterators() |
235 {
236 indices.resize(cell->get_fe().n_dofs_per_cell());
237 cell->get_dof_indices(indices);
238
239 task(cell);
240 }
241 }
242 else
243 {
244 for (const auto &cell :
245 dof_handler.cell_iterators_on_level(mg_level))
246 if (cell->is_locally_owned_on_level())
247 {
248 indices.resize(cell->get_fe().n_dofs_per_cell());
249 cell->get_mg_dof_indices(indices);
250
251 task(cell);
252 }
253 }
254 };
255
256 runner([&](const auto &cell) {
257 const types::fe_index fe_index = cell->active_fe_index();
258 for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
259 if (locally_owned_dofs.is_element(indices[i]))
260 dofs_by_component[locally_owned_dofs.index_within_set(indices[i])] =
261 local_component_association[fe_index][i];
262 });
263 }
264
265
266 // this is the function corresponding to the one above but working on
267 // blocks instead of components.
268 //
269 // the output array dofs_by_block lists for each dof the corresponding
270 // vector block. if the DoFHandler is based on a parallel distributed
271 // triangulation then the output array is index by
272 // dof.locally_owned_dofs().index_within_set(indices[i])
273 template <int dim, int spacedim>
274 inline void
276 std::vector<unsigned char> &dofs_by_block)
277 {
278 const ::hp::FECollection<dim, spacedim> &fe_collection =
279 dof.get_fe_collection();
280 Assert(fe_collection.n_components() < 256, ExcNotImplemented());
281 Assert(dofs_by_block.size() == dof.n_locally_owned_dofs(),
282 ExcDimensionMismatch(dofs_by_block.size(),
283 dof.n_locally_owned_dofs()));
284
285 // next set up a table for the degrees of freedom on each of the
286 // cells (regardless of the fact whether it is listed in the
287 // component_select argument or not)
288 //
289 // for each element 'f' of the FECollection,
290 // local_block_association[f][d] then returns the vector block that
291 // degree of freedom 'd' belongs to
292 std::vector<std::vector<unsigned char>> local_block_association(
293 fe_collection.size());
294 for (unsigned int f = 0; f < fe_collection.size(); ++f)
295 {
296 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
297 local_block_association[f].resize(fe.n_dofs_per_cell(),
298 static_cast<unsigned char>(-1));
299 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
300 local_block_association[f][i] = fe.system_to_block_index(i).first;
301
302 Assert(std::find(local_block_association[f].begin(),
303 local_block_association[f].end(),
304 static_cast<unsigned char>(-1)) ==
305 local_block_association[f].end(),
307 }
308
309 // then loop over all cells and do the work
310 std::vector<types::global_dof_index> indices;
311 for (const auto &cell : dof.active_cell_iterators())
312 if (cell->is_locally_owned())
313 {
314 const types::fe_index fe_index = cell->active_fe_index();
315 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
316 indices.resize(dofs_per_cell);
317 cell->get_dof_indices(indices);
318 for (unsigned int i = 0; i < dofs_per_cell; ++i)
319 if (dof.locally_owned_dofs().is_element(indices[i]))
320 dofs_by_block[dof.locally_owned_dofs().index_within_set(
321 indices[i])] = local_block_association[fe_index][i];
322 }
323 }
324 } // namespace internal
325
326
327
328 template <int dim, int spacedim, typename Number>
329 void
331 const Vector<Number> &cell_data,
332 Vector<double> &dof_data,
333 const unsigned int component)
334 {
335 const Triangulation<dim, spacedim> &tria = dof_handler.get_triangulation();
336 (void)tria;
337
338 AssertDimension(cell_data.size(), tria.n_active_cells());
339 AssertDimension(dof_data.size(), dof_handler.n_dofs());
340 const auto &fe_collection = dof_handler.get_fe_collection();
341 AssertIndexRange(component, fe_collection.n_components());
342 for (unsigned int i = 0; i < fe_collection.size(); ++i)
343 {
344 Assert(fe_collection[i].is_primitive() == true,
346 }
347
348 // store a flag whether we should care about different components. this
349 // is just a simplification, we could ask for this at every single
350 // place equally well
351 const bool consider_components =
352 (dof_handler.get_fe_collection().n_components() != 1);
353
354 // zero out the components that we will touch
355 if (consider_components == false)
356 dof_data = 0;
357 else
358 {
359 std::vector<unsigned char> component_dofs(
360 dof_handler.n_locally_owned_dofs());
362 dof_handler,
363 fe_collection.component_mask(FEValuesExtractors::Scalar(component)),
364 component_dofs);
365
366 for (unsigned int i = 0; i < dof_data.size(); ++i)
367 if (component_dofs[i] == static_cast<unsigned char>(component))
368 dof_data(i) = 0;
369 }
370
371 // count how often we have added a value in the sum for each dof
372 std::vector<unsigned char> touch_count(dof_handler.n_dofs(), 0);
373
374 std::vector<types::global_dof_index> dof_indices;
375 dof_indices.reserve(fe_collection.max_dofs_per_cell());
376
377 for (const auto &cell : dof_handler.active_cell_iterators())
378 {
379 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
380 dof_indices.resize(dofs_per_cell);
381 cell->get_dof_indices(dof_indices);
382
383 for (unsigned int i = 0; i < dofs_per_cell; ++i)
384 // consider this dof only if it is the right component. if there
385 // is only one component, short cut the test
386 if (!consider_components ||
387 (cell->get_fe().system_to_component_index(i).first == component))
388 {
389 // sum up contribution of the present_cell to this dof
390 dof_data(dof_indices[i]) += cell_data(cell->active_cell_index());
391
392 // note that we added another summand
393 ++touch_count[dof_indices[i]];
394 }
395 }
396
397 // compute the mean value on all the dofs by dividing with the number
398 // of summands.
399 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
400 {
401 // assert that each dof was used at least once. this needs not be
402 // the case if the vector has more than one component
403 Assert(consider_components || (touch_count[i] != 0),
405 if (touch_count[i] != 0)
406 dof_data(i) /= touch_count[i];
407 }
408 }
409
410
411
412 template <int dim, int spacedim>
415 const ComponentMask &component_mask)
416 {
417 Assert(component_mask.represents_n_components(
420 "The given component mask is not sized correctly to represent the "
421 "components of the given finite element."));
422
423 // Two special cases: no component is selected, and all components are
424 // selected; both rather stupid, but easy to catch
425 if (component_mask.n_selected_components(
426 dof.get_fe_collection().n_components()) == 0)
427 return IndexSet(dof.n_dofs());
428 else if (component_mask.n_selected_components(
431 return dof.locally_owned_dofs();
432
433 // get the component association of each DoF and then select the ones
434 // that match the given set of components
435 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs());
436 internal::get_component_association(dof, component_mask, dofs_by_component);
437
438 // fill the selected components in a vector
439 std::vector<types::global_dof_index> selected_dofs;
440 selected_dofs.reserve(dof.n_locally_owned_dofs());
441 for (types::global_dof_index i = 0; i < dofs_by_component.size(); ++i)
442 if (component_mask[dofs_by_component[i]] == true)
443 selected_dofs.push_back(dof.locally_owned_dofs().nth_index_in_set(i));
444
445 // fill vector of indices to return argument
446 IndexSet result(dof.n_dofs());
447 result.add_indices(selected_dofs.begin(), selected_dofs.end());
448 return result;
449 }
450
451
452
453 template <int dim, int spacedim>
456 const BlockMask &block_mask)
457 {
458 // simply forward to the function that works based on a component mask
459 return extract_dofs<dim, spacedim>(
460 dof, dof.get_fe_collection().component_mask(block_mask));
461 }
462
463
464
465 template <int dim, int spacedim>
466 std::vector<IndexSet>
468 const ComponentMask &component_mask)
469 {
470 const auto n_comps = dof.get_fe_collection().n_components();
471 Assert(component_mask.represents_n_components(n_comps),
473 "The given component mask is not sized correctly to represent the "
474 "components of the given finite element."));
475
476 const auto &locally_owned_dofs = dof.locally_owned_dofs();
477
478 // get the component association of each DoF and then select the ones
479 // that match the given set of components
480 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs());
481 internal::get_component_association(dof, component_mask, dofs_by_component);
482
483 std::vector<IndexSet> index_per_comp(n_comps, IndexSet(dof.n_dofs()));
484
485 for (types::global_dof_index i = 0; i < dof.n_locally_owned_dofs(); ++i)
486 {
487 const auto &comp_i = dofs_by_component[i];
488 if (component_mask[comp_i])
489 index_per_comp[comp_i].add_index(
490 locally_owned_dofs.nth_index_in_set(i));
491 }
492 for (auto &c : index_per_comp)
493 c.compress();
494 return index_per_comp;
495 }
496
497
498
499 template <int dim, int spacedim>
500 void
501 extract_level_dofs(const unsigned int level,
502 const DoFHandler<dim, spacedim> &dof,
503 const ComponentMask &component_mask,
504 std::vector<bool> &selected_dofs)
505 {
506 const FiniteElement<dim, spacedim> &fe = dof.get_fe();
507
508 Assert(component_mask.represents_n_components(
511 "The given component mask is not sized correctly to represent the "
512 "components of the given finite element."));
513 Assert(selected_dofs.size() == dof.n_dofs(level),
514 ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
515
516 // two special cases: no component is selected, and all components are
517 // selected, both rather stupid, but easy to catch
518 if (component_mask.n_selected_components(
519 dof.get_fe_collection().n_components()) == 0)
520 {
521 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false);
522 return;
523 }
524 else if (component_mask.n_selected_components(
527 {
528 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), true);
529 return;
530 }
531
532 // preset all values by false
533 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false);
534
535 // next set up a table for the degrees of freedom on each of the cells
536 // whether it is something interesting or not
537 std::vector<unsigned char> local_component_association =
539 std::vector<bool> local_selected_dofs(fe.n_dofs_per_cell());
540 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
541 local_selected_dofs[i] = component_mask[local_component_association[i]];
542
543 // then loop over all cells and do work
544 std::vector<types::global_dof_index> indices(fe.n_dofs_per_cell());
545 for (const auto &cell : dof.cell_iterators_on_level(level))
546 {
547 cell->get_mg_dof_indices(indices);
548 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
549 selected_dofs[indices[i]] = local_selected_dofs[i];
550 }
551 }
552
553
554
555 template <int dim, int spacedim>
556 void
557 extract_level_dofs(const unsigned int level,
558 const DoFHandler<dim, spacedim> &dof,
559 const BlockMask &block_mask,
560 std::vector<bool> &selected_dofs)
561 {
562 // simply defer to the other extract_level_dofs() function
564 dof,
565 dof.get_fe().component_mask(block_mask),
566 selected_dofs);
567 }
568
569
570
571 template <int dim, int spacedim>
572 void
574 const ComponentMask &component_mask,
575 std::vector<bool> &selected_dofs,
576 const std::set<types::boundary_id> &boundary_ids)
577 {
578 Assert((dynamic_cast<
580 &dof_handler.get_triangulation()) == nullptr),
582 "This function can not be used with distributed triangulations. "
583 "See the documentation for more information."));
584
585 IndexSet indices =
586 extract_boundary_dofs(dof_handler, component_mask, boundary_ids);
587
588 // clear and reset array by default values
589 selected_dofs.clear();
590 selected_dofs.resize(dof_handler.n_dofs(), false);
591
592 // then convert the values computed above to the binary vector
593 indices.fill_binary_vector(selected_dofs);
594 }
595
596
597
598 template <int dim, int spacedim>
599 void
601 const ComponentMask &component_mask,
602 IndexSet &selected_dofs,
603 const std::set<types::boundary_id> &boundary_ids)
604 {
605 // Simply forward to the other function
606 selected_dofs =
607 extract_boundary_dofs(dof_handler, component_mask, boundary_ids);
608 }
609
610
611
612 template <int dim, int spacedim>
615 const ComponentMask &component_mask,
616 const std::set<types::boundary_id> &boundary_ids)
617 {
618 Assert(component_mask.represents_n_components(
619 dof_handler.get_fe_collection().n_components()),
620 ExcMessage("Component mask has invalid size."));
621 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
622 boundary_ids.end(),
624
625 IndexSet selected_dofs(dof_handler.n_dofs());
626
627 // let's see whether we have to check for certain boundary indicators
628 // or whether we can accept all
629 const bool check_boundary_id = (boundary_ids.size() != 0);
630
631 // also see whether we have to check whether a certain vector component
632 // is selected, or all
633 const bool check_vector_component =
634 ((component_mask.represents_the_all_selected_mask() == false) ||
635 (component_mask.n_selected_components(
636 dof_handler.get_fe_collection().n_components()) !=
637 dof_handler.get_fe_collection().n_components()));
638
639 std::vector<types::global_dof_index> face_dof_indices;
640 face_dof_indices.reserve(
641 dof_handler.get_fe_collection().max_dofs_per_face());
642
643 // now loop over all cells and check whether their faces are at the
644 // boundary. note that we need not take special care of single lines
645 // being at the boundary (using @p{cell->has_boundary_lines}), since we
646 // do not support boundaries of dimension dim-2, and so every isolated
647 // boundary line is also part of a boundary face which we will be
648 // visiting sooner or later
649 for (const auto &cell : dof_handler.active_cell_iterators())
650 // only work on cells that are either locally owned or at least ghost
651 // cells
652 if (cell->is_artificial() == false)
653 for (const unsigned int face : cell->face_indices())
654 if (cell->at_boundary(face))
655 if (!check_boundary_id ||
656 (boundary_ids.find(cell->face(face)->boundary_id()) !=
657 boundary_ids.end()))
658 {
659 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
660
661 const auto reference_cell = cell->reference_cell();
662
663 const unsigned int n_vertices_per_cell =
664 reference_cell.n_vertices();
665 const unsigned int n_lines_per_cell = reference_cell.n_lines();
666 const unsigned int n_vertices_per_face =
667 reference_cell.face_reference_cell(face).n_vertices();
668 const unsigned int n_lines_per_face =
669 reference_cell.face_reference_cell(face).n_lines();
670
671 const unsigned int dofs_per_face = fe.n_dofs_per_face(face);
672 face_dof_indices.resize(dofs_per_face);
673 cell->face(face)->get_dof_indices(face_dof_indices,
674 cell->active_fe_index());
675
676 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
677 if (!check_vector_component)
678 selected_dofs.add_index(face_dof_indices[i]);
679 else
680 // check for component is required. somewhat tricky as
681 // usual for the case that the shape function is
682 // non-primitive, but use usual convention (see docs)
683 {
684 // first get at the cell-global number of a face dof,
685 // to ask the FE certain questions
686 const unsigned int cell_index =
687 (dim == 1 ?
688 i :
689 (dim == 2 ?
690 (i < 2 * fe.n_dofs_per_vertex() ?
691 i :
692 i + 2 * fe.n_dofs_per_vertex()) :
693 (dim == 3 ? (i < n_vertices_per_face *
694 fe.n_dofs_per_vertex() ?
695 i :
696 (i < n_vertices_per_face *
697 fe.n_dofs_per_vertex() +
698 n_lines_per_face *
699 fe.n_dofs_per_line() ?
700 (i - n_vertices_per_face *
701 fe.n_dofs_per_vertex()) +
702 n_vertices_per_cell *
703 fe.n_dofs_per_vertex() :
704 (i -
705 n_vertices_per_face *
706 fe.n_dofs_per_vertex() -
707 n_lines_per_face *
708 fe.n_dofs_per_line()) +
709 n_vertices_per_cell *
710 fe.n_dofs_per_vertex() +
711 n_lines_per_cell *
712 fe.n_dofs_per_line())) :
714 if (fe.is_primitive(cell_index))
715 {
716 if (component_mask[fe.face_system_to_component_index(
717 i, face)
718 .first] == true)
719 selected_dofs.add_index(face_dof_indices[i]);
720 }
721 else // not primitive
722 {
723 const unsigned int first_nonzero_comp =
726 Assert(first_nonzero_comp < fe.n_components(),
728
729 if (component_mask[first_nonzero_comp] == true)
730 selected_dofs.add_index(face_dof_indices[i]);
731 }
732 }
733 }
734
735 return selected_dofs;
736 }
737
738
739
740 template <int dim, int spacedim>
741 void
743 const DoFHandler<dim, spacedim> &dof_handler,
744 const ComponentMask &component_mask,
745 std::vector<bool> &selected_dofs,
746 const std::set<types::boundary_id> &boundary_ids)
747 {
748 Assert(component_mask.represents_n_components(
749 dof_handler.get_fe_collection().n_components()),
750 ExcMessage("This component mask has the wrong size."));
751 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
752 boundary_ids.end(),
754
755 // let's see whether we have to check for certain boundary indicators
756 // or whether we can accept all
757 const bool check_boundary_id = (boundary_ids.size() != 0);
758
759 // also see whether we have to check whether a certain vector component
760 // is selected, or all
761 const bool check_vector_component =
762 (component_mask.represents_the_all_selected_mask() == false);
763
764 // clear and reset array by default values
765 selected_dofs.clear();
766 selected_dofs.resize(dof_handler.n_dofs(), false);
767 std::vector<types::global_dof_index> cell_dof_indices;
768 cell_dof_indices.reserve(
769 dof_handler.get_fe_collection().max_dofs_per_cell());
770
771 // now loop over all cells and check whether their faces are at the
772 // boundary. note that we need not take special care of single lines
773 // being at the boundary (using @p{cell->has_boundary_lines}), since we
774 // do not support boundaries of dimension dim-2, and so every isolated
775 // boundary line is also part of a boundary face which we will be
776 // visiting sooner or later
777 for (const auto &cell : dof_handler.active_cell_iterators())
778 for (const unsigned int face : cell->face_indices())
779 if (cell->at_boundary(face))
780 if (!check_boundary_id ||
781 (boundary_ids.find(cell->face(face)->boundary_id()) !=
782 boundary_ids.end()))
783 {
784 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
785
786 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
787 cell_dof_indices.resize(dofs_per_cell);
788 cell->get_dof_indices(cell_dof_indices);
789
790 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
791 if (fe.has_support_on_face(i, face))
792 {
793 if (!check_vector_component)
794 selected_dofs[cell_dof_indices[i]] = true;
795 else
796 // check for component is required. somewhat tricky
797 // as usual for the case that the shape function is
798 // non-primitive, but use usual convention (see docs)
799 {
800 if (fe.is_primitive(i))
801 selected_dofs[cell_dof_indices[i]] =
802 (component_mask[fe.system_to_component_index(i)
803 .first] == true);
804 else // not primitive
805 {
806 const unsigned int first_nonzero_comp =
809 Assert(first_nonzero_comp < fe.n_components(),
811
812 selected_dofs[cell_dof_indices[i]] =
813 (component_mask[first_nonzero_comp] == true);
814 }
815 }
816 }
817 }
818 }
819
820
821
822 template <int dim, int spacedim, typename number>
825 const DoFHandler<dim, spacedim> &dof_handler,
826 const std::function<
828 &predicate,
830 {
831 const std::function<bool(
833 predicate_local =
834 [=](
836 -> bool { return cell->is_locally_owned() && predicate(cell); };
837
838 std::vector<types::global_dof_index> local_dof_indices;
839 local_dof_indices.reserve(
840 dof_handler.get_fe_collection().max_dofs_per_cell());
841
842 // Get all the dofs that live on the subdomain:
843 std::set<types::global_dof_index> predicate_dofs;
844
845 for (const auto &cell : dof_handler.active_cell_iterators())
846 if (!cell->is_artificial() && predicate(cell))
847 {
848 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
849 cell->get_dof_indices(local_dof_indices);
850 predicate_dofs.insert(local_dof_indices.begin(),
851 local_dof_indices.end());
852 }
853
854 // Get halo layer and accumulate its DoFs
855 std::set<types::global_dof_index> dofs_with_support_on_halo_cells;
856
857 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
858 halo_cells =
859 GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local);
860 for (typename std::vector<typename DoFHandler<dim, spacedim>::
861 active_cell_iterator>::const_iterator it =
862 halo_cells.begin();
863 it != halo_cells.end();
864 ++it)
865 {
866 // skip halo cells that satisfy the given predicate and thereby
867 // shall be a part of the index set on another MPI core.
868 // Those could only be ghost cells with p::d::Tria.
869 if (predicate(*it))
870 {
871 Assert((*it)->is_ghost(), ExcInternalError());
872 continue;
873 }
874
875 const unsigned int dofs_per_cell = (*it)->get_fe().n_dofs_per_cell();
876 local_dof_indices.resize(dofs_per_cell);
877 (*it)->get_dof_indices(local_dof_indices);
878 dofs_with_support_on_halo_cells.insert(local_dof_indices.begin(),
879 local_dof_indices.end());
880 }
881
882 // A complication coming from the fact that DoFs living on locally
883 // owned cells which satisfy predicate may participate in constraints for
884 // DoFs outside of this region.
885 if (cm.n_constraints() > 0)
886 {
887 // Remove DoFs that are part of constraints for DoFs outside
888 // of predicate. Since we will subtract halo_dofs from predicate_dofs,
889 // achieve this by extending halo_dofs with DoFs to which
890 // halo_dofs are constrained.
891 std::set<types::global_dof_index> extra_halo;
892 for (const auto dof : dofs_with_support_on_halo_cells)
893 // if halo DoF is constrained, add all DoFs to which it's constrained
894 // because after resolving constraints, the support of the DoFs that
895 // constrain the current DoF will extend to the halo cells.
896 if (const auto *line_ptr = cm.get_constraint_entries(dof))
897 {
898 const unsigned int line_size = line_ptr->size();
899 for (unsigned int j = 0; j < line_size; ++j)
900 extra_halo.insert((*line_ptr)[j].first);
901 }
902
903 dofs_with_support_on_halo_cells.insert(extra_halo.begin(),
904 extra_halo.end());
905 }
906
907 // Rework our std::set's into IndexSet and subtract halo DoFs from the
908 // predicate_dofs:
909 IndexSet support_set(dof_handler.n_dofs());
910 support_set.add_indices(predicate_dofs.begin(), predicate_dofs.end());
911 support_set.compress();
912
913 IndexSet halo_set(dof_handler.n_dofs());
914 halo_set.add_indices(dofs_with_support_on_halo_cells.begin(),
915 dofs_with_support_on_halo_cells.end());
916 halo_set.compress();
917
918 support_set.subtract_set(halo_set);
919
920 // we intentionally do not want to limit the output to locally owned DoFs.
921 return support_set;
922 }
923
924
925
926 namespace internal
927 {
928 namespace
929 {
930 template <int spacedim>
933 {
934 // there are no hanging nodes in 1d
935 return IndexSet(dof_handler.n_dofs());
936 }
937
938
939 template <int spacedim>
942 {
943 const unsigned int dim = 2;
944
945 IndexSet selected_dofs(dof_handler.n_dofs());
946
947 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
948
949 // this function is similar to the make_sparsity_pattern function,
950 // see there for more information
951 for (const auto &cell : dof_handler.active_cell_iterators())
952 if (!cell->is_artificial())
953 {
954 for (const unsigned int face : cell->face_indices())
955 if (cell->face(face)->has_children())
956 {
958 line = cell->face(face);
959
960 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
961 ++dof)
962 selected_dofs.add_index(
963 line->child(0)->vertex_dof_index(1, dof));
964
965 for (unsigned int child = 0; child < 2; ++child)
966 {
967 if (cell->neighbor_child_on_subface(face, child)
968 ->is_artificial())
969 continue;
970 for (unsigned int dof = 0; dof != fe.n_dofs_per_line();
971 ++dof)
972 selected_dofs.add_index(
973 line->child(child)->dof_index(dof));
974 }
975 }
976 }
977
978 selected_dofs.compress();
979 return selected_dofs;
980 }
981
982
983 template <int spacedim>
986 {
987 const unsigned int dim = 3;
988
989 IndexSet selected_dofs(dof_handler.n_dofs());
990 IndexSet unconstrained_dofs(dof_handler.n_dofs());
991
992 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
993
994 for (const auto &cell : dof_handler.active_cell_iterators())
995 if (!cell->is_artificial())
996 for (auto f : cell->face_indices())
997 {
998 const typename DoFHandler<dim, spacedim>::face_iterator face =
999 cell->face(f);
1000 if (cell->face(f)->has_children())
1001 {
1002 for (unsigned int child = 0; child < 4; ++child)
1003 if (!cell->neighbor_child_on_subface(f, child)
1004 ->is_artificial())
1005 {
1006 // simply take all DoFs that live on this subface
1007 std::vector<types::global_dof_index> ldi(
1008 fe.n_dofs_per_face(f, child));
1009 face->child(child)->get_dof_indices(ldi);
1010 selected_dofs.add_indices(ldi.begin(), ldi.end());
1011 }
1012
1013 // and subtract (in the end) all the indices which a shared
1014 // between this face and its subfaces
1015 for (unsigned int vertex = 0; vertex < 4; ++vertex)
1016 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
1017 ++dof)
1018 unconstrained_dofs.add_index(
1019 face->vertex_dof_index(vertex, dof));
1020 }
1021 }
1022 selected_dofs.subtract_set(unconstrained_dofs);
1023 return selected_dofs;
1024 }
1025 } // namespace
1026 } // namespace internal
1027
1028
1029
1030 template <int dim, int spacedim>
1031 IndexSet
1033 {
1034 return internal::extract_hanging_node_dofs(dof_handler);
1035 }
1036
1037
1038
1039 template <int dim, int spacedim>
1040 void
1042 const types::subdomain_id subdomain_id,
1043 std::vector<bool> &selected_dofs)
1044 {
1045 Assert(selected_dofs.size() == dof_handler.n_dofs(),
1046 ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
1047
1048 // preset all values by false
1049 std::fill_n(selected_dofs.begin(), dof_handler.n_dofs(), false);
1050
1051 std::vector<types::global_dof_index> local_dof_indices;
1052 local_dof_indices.reserve(
1053 dof_handler.get_fe_collection().max_dofs_per_cell());
1054
1055 // this function is similar to the make_sparsity_pattern function, see
1056 // there for more information
1057 for (const auto &cell : dof_handler.active_cell_iterators())
1058 if (cell->subdomain_id() == subdomain_id)
1059 {
1060 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1061 local_dof_indices.resize(dofs_per_cell);
1062 cell->get_dof_indices(local_dof_indices);
1063 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1064 selected_dofs[local_dof_indices[i]] = true;
1065 }
1066 }
1067
1068
1069
1070 template <int dim, int spacedim>
1071 IndexSet
1073 {
1074 // collect all the locally owned dofs
1075 IndexSet dof_set = dof_handler.locally_owned_dofs();
1076
1077 // add the DoF on the adjacent ghost cells to the IndexSet, cache them
1078 // in a set. need to check each dof manually because we can't be sure
1079 // that the dof range of locally_owned_dofs is really contiguous.
1080 std::vector<types::global_dof_index> dof_indices;
1081 std::set<types::global_dof_index> global_dof_indices;
1082
1083 for (const auto &cell : dof_handler.active_cell_iterators() |
1085 {
1086 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1087 cell->get_dof_indices(dof_indices);
1088
1089 for (const types::global_dof_index dof_index : dof_indices)
1090 if (!dof_set.is_element(dof_index))
1091 global_dof_indices.insert(dof_index);
1092 }
1093
1094 dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
1095
1096 dof_set.compress();
1097
1098 return dof_set;
1099 }
1100
1101
1102
1103 template <int dim, int spacedim>
1104 void
1106 IndexSet &dof_set)
1107 {
1108 dof_set = extract_locally_active_dofs(dof_handler);
1109 }
1110
1111
1112
1113 template <int dim, int spacedim>
1114 IndexSet
1116 const DoFHandler<dim, spacedim> &dof_handler,
1117 const unsigned int level)
1118 {
1119 // collect all the locally owned dofs
1120 IndexSet dof_set = dof_handler.locally_owned_mg_dofs(level);
1121
1122 // add the DoF on the adjacent ghost cells to the IndexSet, cache them
1123 // in a set. need to check each dof manually because we can't be sure
1124 // that the dof range of locally_owned_dofs is really contiguous.
1125 std::vector<types::global_dof_index> dof_indices;
1126 std::set<types::global_dof_index> global_dof_indices;
1127
1128 const auto filtered_iterators_range =
1131 for (const auto &cell : filtered_iterators_range)
1132 {
1133 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1134 cell->get_mg_dof_indices(dof_indices);
1135
1136 for (const types::global_dof_index dof_index : dof_indices)
1137 if (!dof_set.is_element(dof_index))
1138 global_dof_indices.insert(dof_index);
1139 }
1140
1141 dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
1142
1143 dof_set.compress();
1144
1145 return dof_set;
1146 }
1147
1148
1149
1150 template <int dim, int spacedim>
1151 void
1153 const DoFHandler<dim, spacedim> &dof_handler,
1154 IndexSet &dof_set,
1155 const unsigned int level)
1156 {
1157 dof_set = extract_locally_active_level_dofs(dof_handler, level);
1158 }
1159
1160
1161
1162 template <int dim, int spacedim>
1163 IndexSet
1165 {
1166 // collect all the locally owned dofs
1167 IndexSet dof_set = dof_handler.locally_owned_dofs();
1168
1169 // now add the DoF on the adjacent ghost cells to the IndexSet
1170
1171 // Note: For certain meshes (in particular in 3d and with many
1172 // processors), it is really necessary to cache intermediate data. After
1173 // trying several objects such as std::set, a vector that is always kept
1174 // sorted, and a vector that is initially unsorted and sorted once at the
1175 // end, the latter has been identified to provide the best performance.
1176 // Martin Kronbichler
1177 std::vector<types::global_dof_index> dof_indices;
1178 std::vector<types::global_dof_index> dofs_on_ghosts;
1179
1180 for (const auto &cell : dof_handler.active_cell_iterators())
1181 if (cell->is_ghost())
1182 {
1183 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1184 cell->get_dof_indices(dof_indices);
1185 for (const auto dof_index : dof_indices)
1186 if (!dof_set.is_element(dof_index))
1187 dofs_on_ghosts.push_back(dof_index);
1188 }
1189
1190 // sort and put into an index set
1191 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1192 dof_set.add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1193 dof_set.compress();
1194
1195 return dof_set;
1196 }
1197
1198
1199
1200 template <int dim, int spacedim>
1201 void
1203 IndexSet &dof_set)
1204 {
1205 dof_set = extract_locally_relevant_dofs(dof_handler);
1206 }
1207
1208
1209
1210 template <int dim, int spacedim>
1211 IndexSet
1213 const DoFHandler<dim, spacedim> &dof_handler,
1214 const unsigned int level)
1215 {
1216 // collect all the locally owned dofs
1217 IndexSet dof_set = dof_handler.locally_owned_mg_dofs(level);
1218
1219 // add the DoF on the adjacent ghost cells to the IndexSet
1220
1221 // Note: For certain meshes (in particular in 3d and with many
1222 // processors), it is really necessary to cache intermediate data. After
1223 // trying several objects such as std::set, a vector that is always kept
1224 // sorted, and a vector that is initially unsorted and sorted once at the
1225 // end, the latter has been identified to provide the best performance.
1226 // Martin Kronbichler
1227 std::vector<types::global_dof_index> dof_indices;
1228 std::vector<types::global_dof_index> dofs_on_ghosts;
1229
1230 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
1231 {
1232 const types::subdomain_id id = cell->level_subdomain_id();
1233
1234 // skip artificial and own cells (only look at ghost cells)
1235 if (id == dof_handler.get_triangulation().locally_owned_subdomain() ||
1237 continue;
1238
1239 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1240 cell->get_mg_dof_indices(dof_indices);
1241 for (const auto dof_index : dof_indices)
1242 if (!dof_set.is_element(dof_index))
1243 dofs_on_ghosts.push_back(dof_index);
1244 }
1245
1246 // sort and fill into an index set
1247 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1248 dof_set.add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1249 dof_set.compress();
1250
1251 return dof_set;
1252 }
1253
1254
1255
1256 template <int dim, int spacedim>
1257 void
1259 const DoFHandler<dim, spacedim> &dof_handler,
1260 const unsigned int level,
1261 IndexSet &dof_set)
1262 {
1263 dof_set = extract_locally_relevant_level_dofs(dof_handler, level);
1264 }
1265
1266
1267 namespace internal
1268 {
1269 template <int dim, int spacedim>
1270 void
1272 const ComponentMask &component_mask,
1273 const unsigned int mg_level,
1274 std::vector<std::vector<bool>> &constant_modes)
1275 {
1276 const auto &locally_owned_dofs =
1277 (mg_level == numbers::invalid_unsigned_int) ?
1278 dof_handler.locally_owned_dofs() :
1279 dof_handler.locally_owned_mg_dofs(mg_level);
1280
1281 // If there are no locally owned DoFs, return with an empty
1282 // constant_modes object:
1283 if (locally_owned_dofs.n_elements() == 0)
1284 {
1285 constant_modes = std::vector<std::vector<bool>>(0);
1286 return;
1287 }
1288
1289 const unsigned int n_components = dof_handler.get_fe(0).n_components();
1290 Assert(component_mask.represents_n_components(n_components),
1291 ExcDimensionMismatch(n_components, component_mask.size()));
1292
1293 std::vector<unsigned char> dofs_by_component(
1294 locally_owned_dofs.n_elements());
1296 component_mask,
1297 dofs_by_component,
1298 mg_level);
1299 unsigned int n_selected_dofs = 0;
1300 for (unsigned int i = 0; i < n_components; ++i)
1301 if (component_mask[i] == true)
1302 n_selected_dofs +=
1303 std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
1304
1305 // Find local numbering within the selected components
1306 std::vector<unsigned int> component_numbering(
1307 locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int);
1308 for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements();
1309 ++i)
1310 if (component_mask[dofs_by_component[i]])
1311 component_numbering[i] = count++;
1312
1313 // get the element constant modes and find a translation table between
1314 // index in the constant modes and the components.
1315 //
1316 // TODO: We might be able to extend this also for elements which do not
1317 // have the same constant modes, but that is messy...
1318 const ::hp::FECollection<dim, spacedim> &fe_collection =
1319 dof_handler.get_fe_collection();
1320 std::vector<Table<2, bool>> element_constant_modes;
1321 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
1322 constant_mode_to_component_translation(n_components);
1323 {
1324 unsigned int n_constant_modes = 0;
1325 int first_non_empty_constant_mode = -1;
1326 for (unsigned int f = 0; f < fe_collection.size(); ++f)
1327 {
1328 std::pair<Table<2, bool>, std::vector<unsigned int>> data =
1329 fe_collection[f].get_constant_modes();
1330
1331 // Store the index of the current element if it is the first that
1332 // has non-empty constant modes.
1333 if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
1334 {
1335 first_non_empty_constant_mode = f;
1336 // This is the first non-empty constant mode, so we figure out
1337 // the translation between index in the constant modes and the
1338 // components
1339 for (unsigned int i = 0; i < data.second.size(); ++i)
1340 if (component_mask[data.second[i]])
1341 constant_mode_to_component_translation[data.second[i]]
1342 .emplace_back(n_constant_modes++, i);
1343 }
1344
1345 // Add the constant modes of this element to the list and assert
1346 // that there are as many constant modes as for the other elements
1347 // (or zero constant modes).
1348 element_constant_modes.push_back(data.first);
1349 Assert(element_constant_modes.back().n_rows() == 0 ||
1350 element_constant_modes.back().n_rows() ==
1351 element_constant_modes[first_non_empty_constant_mode]
1352 .n_rows(),
1354 }
1355 AssertIndexRange(first_non_empty_constant_mode, fe_collection.size());
1356
1357 // Now we know the number of constant modes and resize the return vector
1358 // accordingly
1359 constant_modes.clear();
1360 constant_modes.resize(n_constant_modes,
1361 std::vector<bool>(n_selected_dofs, false));
1362 }
1363
1364 // Loop over all owned cells and ask the element for the constant modes
1365 std::vector<types::global_dof_index> dof_indices;
1366
1367 const auto runner = [&](const auto &task) {
1368 if (mg_level == numbers::invalid_unsigned_int)
1369 {
1370 for (const auto &cell : dof_handler.active_cell_iterators() |
1372 {
1373 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1374 cell->get_dof_indices(dof_indices);
1375
1376 task(cell);
1377 }
1378 }
1379 else
1380 {
1381 for (const auto &cell :
1382 dof_handler.cell_iterators_on_level(mg_level))
1383 if (cell->is_locally_owned_on_level())
1384 {
1385 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1386 cell->get_mg_dof_indices(dof_indices);
1387
1388 task(cell);
1389 }
1390 }
1391 };
1392
1393 runner([&](const auto &cell) {
1394 for (unsigned int i = 0; i < dof_indices.size(); ++i)
1395 if (locally_owned_dofs.is_element(dof_indices[i]))
1396 {
1397 const unsigned int loc_index =
1398 locally_owned_dofs.index_within_set(dof_indices[i]);
1399 const unsigned int comp = dofs_by_component[loc_index];
1400 if (component_mask[comp])
1401 for (auto &indices :
1402 constant_mode_to_component_translation[comp])
1403 constant_modes[indices
1404 .first][component_numbering[loc_index]] =
1405 element_constant_modes[cell->active_fe_index()](
1406 indices.second, i);
1407 }
1408 });
1409 }
1410 } // namespace internal
1411
1412
1413
1414 template <int dim, int spacedim>
1415 void
1417 const ComponentMask &component_mask,
1418 std::vector<std::vector<bool>> &constant_modes)
1419 {
1421 component_mask,
1423 constant_modes);
1424 }
1425
1426
1427
1428 template <int dim, int spacedim>
1429 void
1431 const DoFHandler<dim, spacedim> &dof_handler,
1432 const ComponentMask &component_mask,
1433 std::vector<std::vector<bool>> &constant_modes)
1434 {
1436 component_mask,
1437 level,
1438 constant_modes);
1439 }
1440
1441
1442
1443 template <int dim, int spacedim>
1444 std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1445 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1446 unsigned int>>
1450 &c1_to_c0_face,
1451 const DoFHandler<dim, spacedim> &c0_dh,
1452 const DoFHandler<dim - 1, spacedim> &c1_dh)
1453 {
1454 // This is the returned object: a map of codimension-1 active dof cell
1455 // iterators to codimension-0 cells and face indices
1456 std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1457 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1458 unsigned int>>
1459 c1_to_c0_cell_and_face;
1460
1461 // Shortcut if there are no faces to check
1462 if (c1_to_c0_face.empty())
1463 return c1_to_c0_cell_and_face;
1464
1465 // This is the partial inverse of the map passed as input, for dh
1466 std::map<typename Triangulation<dim, spacedim>::face_iterator,
1467 typename DoFHandler<dim - 1, spacedim>::active_cell_iterator>
1468 c0_to_c1;
1469
1470 // map volume mesh face -> codimension 1 dof cell
1471 for (const auto &[c1_cell, c0_cell] : c1_to_c0_face)
1472 if (!c1_cell->has_children())
1473 c0_to_c1[c0_cell] = c1_cell->as_dof_handler_iterator(c1_dh);
1474
1475 // generate a mapping that maps codimension-1 cells
1476 // to codimension-0 cells and faces
1477 for (const auto &cell :
1478 c0_dh.active_cell_iterators()) // disp_dof.active_cell_iterators())
1479 for (const auto f : cell->face_indices())
1480 if (cell->face(f)->at_boundary())
1481 {
1482 const auto &it = c0_to_c1.find(cell->face(f));
1483 if (it != c0_to_c1.end())
1484 {
1485 const auto &c1_cell = it->second;
1486 c1_to_c0_cell_and_face[c1_cell] = {cell, f};
1487 c0_to_c1.erase(it);
1488 }
1489 }
1490 // Check the dimensions: make sure all active cells we had have been mapped.
1491 AssertDimension(c0_to_c1.size(), 0);
1492 return c1_to_c0_cell_and_face;
1493 }
1494
1495
1496
1497 template <int dim, int spacedim>
1498 void
1500 std::vector<unsigned int> &active_fe_indices)
1501 {
1502 AssertDimension(active_fe_indices.size(),
1503 dof_handler.get_triangulation().n_active_cells());
1504
1505 std::vector<types::fe_index> indices = dof_handler.get_active_fe_indices();
1506
1507 active_fe_indices.assign(indices.begin(), indices.end());
1508 }
1509
1510 template <int dim, int spacedim>
1511 std::vector<IndexSet>
1513 {
1514 Assert(dof_handler.n_dofs() > 0,
1515 ExcMessage("The given DoFHandler has no DoFs."));
1516
1517 // If the Triangulation is distributed, the only thing we can usefully
1518 // ask is for its locally owned subdomain
1519 Assert((dynamic_cast<
1521 &dof_handler.get_triangulation()) == nullptr),
1522 ExcMessage(
1523 "For parallel::distributed::Triangulation objects and "
1524 "associated DoF handler objects, asking for any information "
1525 "related to a subdomain other than the locally owned one does "
1526 "not make sense."));
1527
1528 // The following is a random process (flip of a coin), thus should be called
1529 // once only.
1530 std::vector<::types::subdomain_id> subdomain_association(
1531 dof_handler.n_dofs());
1533 subdomain_association);
1534
1535 // Figure out how many subdomain ids there are.
1536 //
1537 // if this is a parallel triangulation, then we can just ask the
1538 // triangulation for this. if this is a sequential triangulation, we loop
1539 // over all cells and take the largest subdomain_id value we find; the
1540 // number of subdomains is then the largest found value plus one. (we here
1541 // assume that all subdomain ids up to the largest are actually used; this
1542 // may not be true for a sequential triangulation where these values have
1543 // been set by hand and not in accordance with some MPI communicator; but
1544 // the function returns an array indexed starting at zero, so we need to
1545 // collect information for each subdomain index anyway, not just for the
1546 // used one.)
1547 const unsigned int n_subdomains =
1548 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1549 &dof_handler.get_triangulation()) == nullptr ?
1550 [&dof_handler]() {
1551 unsigned int max_subdomain_id = 0;
1552 for (const auto &cell : dof_handler.active_cell_iterators())
1553 max_subdomain_id =
1554 std::max(max_subdomain_id, cell->subdomain_id());
1555 return max_subdomain_id + 1;
1556 }() :
1558 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1559 &dof_handler.get_triangulation())
1560 ->get_communicator()));
1561 Assert(n_subdomains > *std::max_element(subdomain_association.begin(),
1562 subdomain_association.end()),
1564
1565 std::vector<::IndexSet> index_sets(
1566 n_subdomains, ::IndexSet(dof_handler.n_dofs()));
1567
1568 // loop over subdomain_association and populate IndexSet when a
1569 // change in subdomain ID is found
1570 ::types::global_dof_index i_min = 0;
1571 ::types::global_dof_index this_subdomain = subdomain_association[0];
1572
1573 for (::types::global_dof_index index = 1;
1574 index < subdomain_association.size();
1575 ++index)
1576 {
1577 // found index different from the current one
1578 if (subdomain_association[index] != this_subdomain)
1579 {
1580 index_sets[this_subdomain].add_range(i_min, index);
1581 i_min = index;
1582 this_subdomain = subdomain_association[index];
1583 }
1584 }
1585
1586 // the very last element is of different index
1587 if (i_min == subdomain_association.size() - 1)
1588 {
1589 index_sets[this_subdomain].add_index(i_min);
1590 }
1591
1592 // otherwise there are at least two different indices
1593 else
1594 {
1595 index_sets[this_subdomain].add_range(i_min,
1596 subdomain_association.size());
1597 }
1598
1599 for (unsigned int i = 0; i < n_subdomains; ++i)
1600 index_sets[i].compress();
1601
1602 return index_sets;
1603 }
1604
1605 template <int dim, int spacedim>
1606 std::vector<IndexSet>
1608 const DoFHandler<dim, spacedim> &dof_handler)
1609 {
1610 // If the Triangulation is distributed, the only thing we can usefully
1611 // ask is for its locally owned subdomain
1612 Assert((dynamic_cast<
1614 &dof_handler.get_triangulation()) == nullptr),
1615 ExcMessage(
1616 "For parallel::distributed::Triangulation objects and "
1617 "associated DoF handler objects, asking for any information "
1618 "related to a subdomain other than the locally owned one does "
1619 "not make sense."));
1620
1621 // Collect all the locally owned DoFs
1622 // Note: Even though the distribution of DoFs by the
1623 // locally_owned_dofs_per_subdomain function is pseudo-random, we will
1624 // collect all the DoFs on the subdomain and its layer cell. Therefore, the
1625 // random nature of this function does not play a role in the extraction of
1626 // the locally relevant DoFs
1627 std::vector<IndexSet> dof_set =
1629 const ::types::subdomain_id n_subdomains = dof_set.size();
1630
1631 // Add the DoFs on the adjacent (equivalent ghost) cells to the IndexSet,
1632 // cache them in a set. Need to check each DoF manually because we can't
1633 // be sure that the DoF range of locally_owned_dofs is really contiguous.
1634 for (::types::subdomain_id subdomain_id = 0;
1635 subdomain_id < n_subdomains;
1636 ++subdomain_id)
1637 {
1638 // Extract the layer of cells around this subdomain
1639 std::function<bool(
1641 predicate = IteratorFilters::SubdomainEqualTo(subdomain_id);
1642 const std::vector<
1644 active_halo_layer =
1645 GridTools::compute_active_cell_halo_layer(dof_handler, predicate);
1646
1647 // Extract DoFs associated with halo layer
1648 std::vector<types::global_dof_index> local_dof_indices;
1649 std::set<types::global_dof_index> subdomain_halo_global_dof_indices;
1650 for (typename std::vector<
1652 const_iterator it_cell = active_halo_layer.begin();
1653 it_cell != active_halo_layer.end();
1654 ++it_cell)
1655 {
1657 &cell = *it_cell;
1658 Assert(
1659 cell->subdomain_id() != subdomain_id,
1660 ExcMessage(
1661 "The subdomain ID of the halo cell should not match that of the vector entry."));
1662
1663 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1664 cell->get_dof_indices(local_dof_indices);
1665
1666 for (const types::global_dof_index local_dof_index :
1667 local_dof_indices)
1668 subdomain_halo_global_dof_indices.insert(local_dof_index);
1669 }
1670
1671 dof_set[subdomain_id].add_indices(
1672 subdomain_halo_global_dof_indices.begin(),
1673 subdomain_halo_global_dof_indices.end());
1674
1675 dof_set[subdomain_id].compress();
1676 }
1677
1678 return dof_set;
1679 }
1680
1681 template <int dim, int spacedim>
1682 void
1684 const DoFHandler<dim, spacedim> &dof_handler,
1685 std::vector<types::subdomain_id> &subdomain_association)
1686 {
1687 // if the Triangulation is distributed, the only thing we can usefully
1688 // ask is for its locally owned subdomain
1689 Assert((dynamic_cast<
1691 &dof_handler.get_triangulation()) == nullptr),
1692 ExcMessage(
1693 "For parallel::distributed::Triangulation objects and "
1694 "associated DoF handler objects, asking for any subdomain other "
1695 "than the locally owned one does not make sense."));
1696
1697 Assert(subdomain_association.size() == dof_handler.n_dofs(),
1698 ExcDimensionMismatch(subdomain_association.size(),
1699 dof_handler.n_dofs()));
1700
1701 // catch an error that happened in some versions of the shared tria
1702 // distribute_dofs() function where we were trying to call this
1703 // function at a point in time when not all internal DoFHandler
1704 // structures were quite set up yet.
1705 Assert(dof_handler.n_dofs() > 0, ExcInternalError());
1706
1707 // In case this function is executed with parallel::shared::Triangulation
1708 // with possibly artificial cells, we need to take "true" subdomain IDs
1709 // (i.e. without artificial cells). Otherwise we are good to use
1710 // subdomain_id as stored in cell->subdomain_id().
1711 std::vector<types::subdomain_id> cell_owners(
1712 dof_handler.get_triangulation().n_active_cells());
1715 &dof_handler.get_triangulation())))
1716 {
1717 cell_owners = tr->get_true_subdomain_ids_of_cells();
1718 Assert(tr->get_true_subdomain_ids_of_cells().size() ==
1719 tr->n_active_cells(),
1721 }
1722 else
1723 {
1724 for (const auto &cell : dof_handler.active_cell_iterators() |
1726 cell_owners[cell->active_cell_index()] = cell->subdomain_id();
1727 }
1728
1729 // preset all values by an invalid value
1730 std::fill_n(subdomain_association.begin(),
1731 dof_handler.n_dofs(),
1733
1734 std::vector<types::global_dof_index> local_dof_indices;
1735 local_dof_indices.reserve(
1736 dof_handler.get_fe_collection().max_dofs_per_cell());
1737
1738 // loop over all cells and record which subdomain a DoF belongs to.
1739 // give to the smaller subdomain_id in case it is on an interface
1740 for (const auto &cell : dof_handler.active_cell_iterators())
1741 {
1742 const types::subdomain_id subdomain_id =
1743 cell_owners[cell->active_cell_index()];
1744 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1745 local_dof_indices.resize(dofs_per_cell);
1746 cell->get_dof_indices(local_dof_indices);
1747
1748 // set subdomain ids. if dofs already have their values set then
1749 // they must be on partition interfaces. in that case assign them
1750 // to either the previous association or the current processor
1751 // with the smaller subdomain id.
1752 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1753 if (subdomain_association[local_dof_indices[i]] ==
1755 subdomain_association[local_dof_indices[i]] = subdomain_id;
1756 else if (subdomain_association[local_dof_indices[i]] > subdomain_id)
1757 {
1758 subdomain_association[local_dof_indices[i]] = subdomain_id;
1759 }
1760 }
1761
1762 Assert(std::find(subdomain_association.begin(),
1763 subdomain_association.end(),
1765 subdomain_association.end(),
1767 }
1768
1769
1770
1771 template <int dim, int spacedim>
1772 unsigned int
1774 const DoFHandler<dim, spacedim> &dof_handler,
1775 const types::subdomain_id subdomain)
1776 {
1777 std::vector<types::subdomain_id> subdomain_association(
1778 dof_handler.n_dofs());
1779 get_subdomain_association(dof_handler, subdomain_association);
1780
1781 return std::count(subdomain_association.begin(),
1782 subdomain_association.end(),
1783 subdomain);
1784 }
1785
1786
1787
1788 template <int dim, int spacedim>
1789 IndexSet
1791 const DoFHandler<dim, spacedim> &dof_handler,
1792 const types::subdomain_id subdomain)
1793 {
1794 // If we have a distributed::Triangulation only allow locally_owned
1795 // subdomain.
1798 (subdomain ==
1800 ExcMessage(
1801 "For parallel::distributed::Triangulation objects and "
1802 "associated DoF handler objects, asking for any subdomain other "
1803 "than the locally owned one does not make sense."));
1804
1805 IndexSet index_set(dof_handler.n_dofs());
1806
1807 std::vector<types::global_dof_index> local_dof_indices;
1808 local_dof_indices.reserve(
1809 dof_handler.get_fe_collection().max_dofs_per_cell());
1810
1811 // first generate an unsorted list of all indices which we fill from
1812 // the back. could also insert them directly into the IndexSet, but
1813 // that inserts indices in the middle, which is an O(n^2) algorithm and
1814 // hence too expensive. Could also use std::set, but that is in general
1815 // more expensive than a vector
1816 std::vector<types::global_dof_index> subdomain_indices;
1817
1818 for (const auto &cell : dof_handler.active_cell_iterators())
1819 if ((cell->is_artificial() == false) &&
1820 (cell->subdomain_id() == subdomain))
1821 {
1822 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1823 local_dof_indices.resize(dofs_per_cell);
1824 cell->get_dof_indices(local_dof_indices);
1825 subdomain_indices.insert(subdomain_indices.end(),
1826 local_dof_indices.begin(),
1827 local_dof_indices.end());
1828 }
1829 // sort indices and put into an index set:
1830 std::sort(subdomain_indices.begin(), subdomain_indices.end());
1831 index_set.add_indices(subdomain_indices.begin(), subdomain_indices.end());
1832 index_set.compress();
1833
1834 return index_set;
1835 }
1836
1837
1838
1839 template <int dim, int spacedim>
1840 void
1842 const DoFHandler<dim, spacedim> &dof_handler,
1843 const types::subdomain_id subdomain,
1844 std::vector<unsigned int> &n_dofs_on_subdomain)
1845 {
1846 Assert(n_dofs_on_subdomain.size() == dof_handler.get_fe(0).n_components(),
1847 ExcDimensionMismatch(n_dofs_on_subdomain.size(),
1848 dof_handler.get_fe(0).n_components()));
1849 std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
1850
1851 // Make sure there are at least some cells with this subdomain id
1852 Assert(std::any_of(
1853 dof_handler.begin_active(),
1855 dof_handler.end()},
1856 [subdomain](
1857 const typename DoFHandler<dim, spacedim>::cell_accessor &cell) {
1858 return cell.subdomain_id() == subdomain;
1859 }),
1860 ExcMessage("There are no cells for the given subdomain!"));
1861
1862 std::vector<types::subdomain_id> subdomain_association(
1863 dof_handler.n_dofs());
1864 get_subdomain_association(dof_handler, subdomain_association);
1865
1866 std::vector<unsigned char> component_association(dof_handler.n_dofs());
1868 ComponentMask(std::vector<bool>()),
1869 component_association);
1870
1871 for (unsigned int c = 0; c < dof_handler.get_fe(0).n_components(); ++c)
1872 {
1873 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
1874 if ((subdomain_association[i] == subdomain) &&
1875 (component_association[i] == static_cast<unsigned char>(c)))
1876 ++n_dofs_on_subdomain[c];
1877 }
1878 }
1879
1880
1881
1882 namespace internal
1883 {
1884 // TODO: why is this function so complicated? It would be nice to have
1885 // comments that explain why we can't just loop over all components and
1886 // count the entries in dofs_by_component that have this component's
1887 // index
1888 template <int dim, int spacedim>
1889 void
1891 const std::vector<unsigned char> &dofs_by_component,
1892 const std::vector<unsigned int> &target_component,
1893 const bool only_once,
1894 std::vector<types::global_dof_index> &dofs_per_component,
1895 unsigned int &component)
1896 {
1897 for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
1898 {
1899 const FiniteElement<dim, spacedim> &base = fe.base_element(b);
1900 // Dimension of base element
1901 unsigned int d = base.n_components();
1902
1903 for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
1904 {
1905 if (base.n_base_elements() > 1)
1906 resolve_components(base,
1907 dofs_by_component,
1908 target_component,
1909 only_once,
1910 dofs_per_component,
1911 component);
1912 else
1913 {
1914 for (unsigned int dd = 0; dd < d; ++dd, ++component)
1915 dofs_per_component[target_component[component]] +=
1916 std::count(dofs_by_component.begin(),
1917 dofs_by_component.end(),
1918 component);
1919
1920 // if we have non-primitive FEs and want all components
1921 // to show the number of dofs, need to copy the result to
1922 // those components
1923 if (!base.is_primitive() && !only_once)
1924 for (unsigned int dd = 1; dd < d; ++dd)
1925 dofs_per_component[target_component[component - d + dd]] =
1926 dofs_per_component[target_component[component - d]];
1927 }
1928 }
1929 }
1930 }
1931
1932
1933 template <int dim, int spacedim>
1934 void
1936 const std::vector<unsigned char> &dofs_by_component,
1937 const std::vector<unsigned int> &target_component,
1938 const bool only_once,
1939 std::vector<types::global_dof_index> &dofs_per_component,
1940 unsigned int &component)
1941 {
1942 // assert that all elements in the collection have the same structure
1943 // (base elements and multiplicity, components per base element) and
1944 // then simply call the function above
1945 for (unsigned int fe = 1; fe < fe_collection.size(); ++fe)
1946 {
1947 Assert(fe_collection[fe].n_components() ==
1948 fe_collection[0].n_components(),
1950 Assert(fe_collection[fe].n_base_elements() ==
1951 fe_collection[0].n_base_elements(),
1953 for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b)
1954 {
1955 Assert(fe_collection[fe].base_element(b).n_components() ==
1956 fe_collection[0].base_element(b).n_components(),
1958 Assert(fe_collection[fe].base_element(b).n_base_elements() ==
1959 fe_collection[0].base_element(b).n_base_elements(),
1961 }
1962 }
1963
1964 resolve_components(fe_collection[0],
1965 dofs_by_component,
1966 target_component,
1967 only_once,
1968 dofs_per_component,
1969 component);
1970 }
1971 } // namespace internal
1972
1973
1974
1975 namespace internal
1976 {
1977 namespace
1978 {
1982 template <int dim, int spacedim>
1983 bool
1984 all_elements_are_primitive(const FiniteElement<dim, spacedim> &fe)
1985 {
1986 return fe.is_primitive();
1987 }
1988
1989
1994 template <int dim, int spacedim>
1995 bool
1996 all_elements_are_primitive(
1997 const ::hp::FECollection<dim, spacedim> &fe_collection)
1998 {
1999 for (unsigned int i = 0; i < fe_collection.size(); ++i)
2000 if (fe_collection[i].is_primitive() == false)
2001 return false;
2002
2003 return true;
2004 }
2005 } // namespace
2006 } // namespace internal
2007
2008
2009
2010 template <int dim, int spacedim>
2011 std::vector<types::global_dof_index>
2013 const DoFHandler<dim, spacedim> &dof_handler,
2014 const bool only_once,
2015 const std::vector<unsigned int> &target_component_)
2016 {
2017 const unsigned int n_components = dof_handler.get_fe(0).n_components();
2018
2019 // If the empty vector was given as default argument, set up this
2020 // vector as identity.
2021 std::vector<unsigned int> target_component = target_component_;
2022 if (target_component.empty())
2023 {
2024 target_component.resize(n_components);
2025 for (unsigned int i = 0; i < n_components; ++i)
2026 target_component[i] = i;
2027 }
2028 else
2029 Assert(target_component.size() == n_components,
2030 ExcDimensionMismatch(target_component.size(), n_components));
2031
2032
2033 const unsigned int max_component =
2034 *std::max_element(target_component.begin(), target_component.end());
2035 const unsigned int n_target_components = max_component + 1;
2036
2037 std::vector<types::global_dof_index> dofs_per_component(
2038 n_target_components, types::global_dof_index(0));
2039
2040 // special case for only one component. treat this first since it does
2041 // not require any computations
2042 if (n_components == 1)
2043 {
2044 dofs_per_component[0] = dof_handler.n_locally_owned_dofs();
2045 return dofs_per_component;
2046 }
2047
2048
2049 // otherwise determine the number of dofs in each component separately.
2050 // do so in parallel
2051 std::vector<unsigned char> dofs_by_component(
2052 dof_handler.n_locally_owned_dofs());
2054 ComponentMask(),
2055 dofs_by_component);
2056
2057 // next count what we got
2058 unsigned int component = 0;
2060 dofs_by_component,
2061 target_component,
2062 only_once,
2063 dofs_per_component,
2064 component);
2065 Assert(n_components == component, ExcInternalError());
2066
2067 // finally sanity check. this is only valid if the finite element is
2068 // actually primitive, so exclude other elements from this
2069 Assert((internal::all_elements_are_primitive(
2070 dof_handler.get_fe_collection()) == false) ||
2071 (std::accumulate(dofs_per_component.begin(),
2072 dofs_per_component.end(),
2074 dof_handler.n_locally_owned_dofs()),
2076
2077 // reduce information from all CPUs
2078#ifdef DEAL_II_WITH_MPI
2079
2081 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2082 &dof_handler.get_triangulation())))
2083 {
2084 std::vector<types::global_dof_index> local_dof_count =
2085 dofs_per_component;
2086
2087 const int ierr = MPI_Allreduce(local_dof_count.data(),
2088 dofs_per_component.data(),
2089 n_target_components,
2091 MPI_SUM,
2092 tria->get_communicator());
2093 AssertThrowMPI(ierr);
2094 }
2095#endif
2096
2097 return dofs_per_component;
2098 }
2099
2100
2101
2102 template <int dim, int spacedim>
2103 std::vector<types::global_dof_index>
2105 const std::vector<unsigned int> &target_block_)
2106 {
2107 const ::hp::FECollection<dim, spacedim> &fe_collection =
2108 dof_handler.get_fe_collection();
2109 Assert(fe_collection.size() < 256, ExcNotImplemented());
2110
2111 // If the empty vector for target_block(e.g., as default argument), then
2112 // set up this vector as identity. We do this set up with the first
2113 // element of the collection, but the whole thing can only work if
2114 // all elements have the same number of blocks anyway -- so check
2115 // that right after
2116 const unsigned int n_blocks = fe_collection[0].n_blocks();
2117
2118 std::vector<unsigned int> target_block = target_block_;
2119 if (target_block.empty())
2120 {
2121 target_block.resize(fe_collection[0].n_blocks());
2122 for (unsigned int i = 0; i < n_blocks; ++i)
2123 target_block[i] = i;
2124 }
2125 else
2126 Assert(target_block.size() == n_blocks,
2127 ExcDimensionMismatch(target_block.size(), n_blocks));
2128 for (unsigned int f = 1; f < fe_collection.size(); ++f)
2129 Assert(fe_collection[0].n_blocks() == fe_collection[f].n_blocks(),
2130 ExcMessage("This function can only work if all elements in a "
2131 "collection have the same number of blocks."));
2132
2133 // special case for only one block. treat this first since it does
2134 // not require any computations
2135 if (n_blocks == 1)
2136 {
2137 std::vector<types::global_dof_index> dofs_per_block(1);
2138 dofs_per_block[0] = dof_handler.n_dofs();
2139 return dofs_per_block;
2140 }
2141
2142 // Otherwise set up the right-sized object and start working
2143 const unsigned int max_block =
2144 *std::max_element(target_block.begin(), target_block.end());
2145 const unsigned int n_target_blocks = max_block + 1;
2146
2147 std::vector<types::global_dof_index> dofs_per_block(n_target_blocks);
2148
2149 // Loop over the elements of the collection, but really only consider
2150 // the last element (see #9271)
2151 for (unsigned int this_fe = fe_collection.size() - 1;
2152 this_fe < fe_collection.size();
2153 ++this_fe)
2154 {
2155 const FiniteElement<dim, spacedim> &fe = fe_collection[this_fe];
2156
2157 std::vector<unsigned char> dofs_by_block(
2158 dof_handler.n_locally_owned_dofs());
2159 internal::get_block_association(dof_handler, dofs_by_block);
2160
2161 // next count what we got
2162 for (unsigned int block = 0; block < fe.n_blocks(); ++block)
2163 dofs_per_block[target_block[block]] +=
2164 std::count(dofs_by_block.begin(), dofs_by_block.end(), block);
2165
2166#ifdef DEAL_II_WITH_MPI
2167 // if we are working on a parallel mesh, we now need to collect
2168 // this information from all processors
2170 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2171 &dof_handler.get_triangulation())))
2172 {
2173 std::vector<types::global_dof_index> local_dof_count =
2174 dofs_per_block;
2175 const int ierr = MPI_Allreduce(local_dof_count.data(),
2176 dofs_per_block.data(),
2177 n_target_blocks,
2179 MPI_SUM,
2180 tria->get_communicator());
2181 AssertThrowMPI(ierr);
2182 }
2183#endif
2184 }
2185
2186 return dofs_per_block;
2187 }
2188
2189
2190
2191 template <int dim, int spacedim>
2192 void
2194 std::vector<types::global_dof_index> &mapping)
2195 {
2196 mapping.clear();
2197 mapping.insert(mapping.end(),
2198 dof_handler.n_dofs(),
2200
2201 std::vector<types::global_dof_index> dofs_on_face;
2202 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2203 types::global_dof_index next_boundary_index = 0;
2204
2205 // now loop over all cells and check whether their faces are at the
2206 // boundary. note that we need not take special care of single lines
2207 // being at the boundary (using @p{cell->has_boundary_lines}), since we
2208 // do not support boundaries of dimension dim-2, and so every isolated
2209 // boundary line is also part of a boundary face which we will be
2210 // visiting sooner or later
2211 for (const auto &cell : dof_handler.active_cell_iterators())
2212 for (const unsigned int f : cell->face_indices())
2213 if (cell->at_boundary(f))
2214 {
2215 const unsigned int dofs_per_face =
2216 cell->get_fe().n_dofs_per_face(f);
2217 dofs_on_face.resize(dofs_per_face);
2218 cell->face(f)->get_dof_indices(dofs_on_face,
2219 cell->active_fe_index());
2220 for (unsigned int i = 0; i < dofs_per_face; ++i)
2221 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2222 mapping[dofs_on_face[i]] = next_boundary_index++;
2223 }
2224
2225 AssertDimension(next_boundary_index, dof_handler.n_boundary_dofs());
2226 }
2227
2228
2229
2230 template <int dim, int spacedim>
2231 void
2233 const std::set<types::boundary_id> &boundary_ids,
2234 std::vector<types::global_dof_index> &mapping)
2235 {
2236 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2237 boundary_ids.end(),
2239
2240 mapping.clear();
2241 mapping.insert(mapping.end(),
2242 dof_handler.n_dofs(),
2244
2245 // return if there is nothing to do
2246 if (boundary_ids.empty())
2247 return;
2248
2249 std::vector<types::global_dof_index> dofs_on_face;
2250 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2251 types::global_dof_index next_boundary_index = 0;
2252
2253 for (const auto &cell : dof_handler.active_cell_iterators())
2254 for (const unsigned int f : cell->face_indices())
2255 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
2256 boundary_ids.end())
2257 {
2258 const unsigned int dofs_per_face =
2259 cell->get_fe().n_dofs_per_face(f);
2260 dofs_on_face.resize(dofs_per_face);
2261 cell->face(f)->get_dof_indices(dofs_on_face,
2262 cell->active_fe_index());
2263 for (unsigned int i = 0; i < dofs_per_face; ++i)
2264 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2265 mapping[dofs_on_face[i]] = next_boundary_index++;
2266 }
2267
2268 AssertDimension(next_boundary_index,
2269 dof_handler.n_boundary_dofs(boundary_ids));
2270 }
2271
2272 namespace internal
2273 {
2274 namespace
2275 {
2276 template <int dim, int spacedim>
2277 std::map<types::global_dof_index, Point<spacedim>>
2280 const DoFHandler<dim, spacedim> &dof_handler,
2281 const ComponentMask &in_mask)
2282 {
2283 std::map<types::global_dof_index, Point<spacedim>> support_points;
2284
2285 const hp::FECollection<dim, spacedim> &fe_collection =
2286 dof_handler.get_fe_collection();
2287 hp::QCollection<dim> q_coll_dummy;
2288
2289 for (unsigned int fe_index = 0; fe_index < fe_collection.size();
2290 ++fe_index)
2291 {
2292 // check whether every FE in the collection has support points
2293 Assert((fe_collection[fe_index].n_dofs_per_cell() == 0) ||
2294 (fe_collection[fe_index].has_support_points()),
2296 q_coll_dummy.push_back(Quadrature<dim>(
2297 fe_collection[fe_index].get_unit_support_points()));
2298 }
2299
2300 // Take care of components
2301 const ComponentMask mask =
2302 (in_mask.size() == 0 ?
2303 ComponentMask(fe_collection.n_components(), true) :
2304 in_mask);
2305
2306 // Now loop over all cells and enquire the support points on each
2307 // of these. we use dummy quadrature formulas where the quadrature
2308 // points are located at the unit support points to enquire the
2309 // location of the support points in real space.
2310 //
2311 // The weights of the quadrature rule have been set to invalid
2312 // values by the used constructor.
2313 hp::FEValues<dim, spacedim> hp_fe_values(mapping,
2314 fe_collection,
2315 q_coll_dummy,
2317
2318 std::vector<types::global_dof_index> local_dof_indices;
2319 for (const auto &cell : dof_handler.active_cell_iterators())
2320 // only work on locally relevant cells
2321 if (cell->is_artificial() == false)
2322 {
2323 hp_fe_values.reinit(cell);
2324 const FEValues<dim, spacedim> &fe_values =
2325 hp_fe_values.get_present_fe_values();
2326
2327 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2328 cell->get_dof_indices(local_dof_indices);
2329
2330 const std::vector<Point<spacedim>> &points =
2331 fe_values.get_quadrature_points();
2332 for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell();
2333 ++i)
2334 {
2335 const unsigned int dof_comp =
2336 cell->get_fe().system_to_component_index(i).first;
2337
2338 // insert the values into the map if it is a valid component
2339 if (mask[dof_comp])
2340 support_points[local_dof_indices[i]] = points[i];
2341 }
2342 }
2343
2344 return support_points;
2345 }
2346
2347
2348 template <int dim, int spacedim>
2349 std::vector<Point<spacedim>>
2350 map_dofs_to_support_points_vector(
2352 const DoFHandler<dim, spacedim> &dof_handler,
2353 const ComponentMask &mask)
2354 {
2355 std::vector<Point<spacedim>> support_points(dof_handler.n_dofs());
2356
2357 // get the data in the form of the map as above
2358 const std::map<types::global_dof_index, Point<spacedim>>
2359 x_support_points =
2360 map_dofs_to_support_points(mapping, dof_handler, mask);
2361
2362 // now convert from the map to the linear vector. make sure every
2363 // entry really appeared in the map
2364 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
2365 {
2366 Assert(x_support_points.find(i) != x_support_points.end(),
2368
2369 support_points[i] = x_support_points.find(i)->second;
2370 }
2371
2372 return support_points;
2373 }
2374 } // namespace
2375 } // namespace internal
2376
2377
2378 template <int dim, int spacedim>
2379 void
2381 const DoFHandler<dim, spacedim> &dof_handler,
2382 std::vector<Point<spacedim>> &support_points,
2383 const ComponentMask &mask)
2384 {
2385 AssertDimension(support_points.size(), dof_handler.n_dofs());
2386 Assert((dynamic_cast<
2388 &dof_handler.get_triangulation()) == nullptr),
2389 ExcMessage(
2390 "This function can not be used with distributed triangulations. "
2391 "See the documentation for more information."));
2392
2393 // Let the internal function do all the work, just make sure that it
2394 // gets a MappingCollection
2395 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2396
2397 support_points =
2398 internal::map_dofs_to_support_points_vector(mapping_collection,
2399 dof_handler,
2400 mask);
2401 }
2402
2403
2404 template <int dim, int spacedim>
2405 void
2408 const DoFHandler<dim, spacedim> &dof_handler,
2409 std::vector<Point<spacedim>> &support_points,
2410 const ComponentMask &mask)
2411 {
2412 AssertDimension(support_points.size(), dof_handler.n_dofs());
2413 Assert((dynamic_cast<
2415 &dof_handler.get_triangulation()) == nullptr),
2416 ExcMessage(
2417 "This function can not be used with distributed triangulations. "
2418 "See the documentation for more information."));
2419
2420 // Let the internal function do all the work, just make sure that it
2421 // gets a MappingCollection
2422 support_points =
2423 internal::map_dofs_to_support_points_vector(mapping, dof_handler, mask);
2424 }
2425
2426
2427 // This function is deprecated:
2428 template <int dim, int spacedim>
2431 const Mapping<dim, spacedim> &mapping,
2432 const DoFHandler<dim, spacedim> &dof_handler,
2433 std::map<types::global_dof_index, Point<spacedim>> &support_points,
2434 const ComponentMask &mask)
2435 {
2436 support_points.clear();
2437
2438 // Let the internal function do all the work, just make sure that it
2439 // gets a MappingCollection
2440 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2441
2442 support_points = internal::map_dofs_to_support_points(mapping_collection,
2443 dof_handler,
2444 mask);
2445 }
2446
2447
2448 // This function is deprecated:
2449 template <int dim, int spacedim>
2453 const DoFHandler<dim, spacedim> &dof_handler,
2454 std::map<types::global_dof_index, Point<spacedim>> &support_points,
2455 const ComponentMask &mask)
2456 {
2457 support_points.clear();
2458
2459 // Let the internal function do all the work, just make sure that it
2460 // gets a MappingCollection
2461 support_points =
2462 internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2463 }
2464
2465
2466 template <int dim, int spacedim>
2467 std::map<types::global_dof_index, Point<spacedim>>
2469 const DoFHandler<dim, spacedim> &dof_handler,
2470 const ComponentMask &mask)
2471 {
2472 // Let the internal function do all the work, just make sure that it
2473 // gets a MappingCollection
2474 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2475
2476 return internal::map_dofs_to_support_points(mapping_collection,
2477 dof_handler,
2478 mask);
2479 }
2480
2481
2482 template <int dim, int spacedim>
2483 std::map<types::global_dof_index, Point<spacedim>>
2486 const DoFHandler<dim, spacedim> &dof_handler,
2487 const ComponentMask &mask)
2488 {
2489 return internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2490 }
2491
2492
2493 template <int spacedim>
2494 void
2496 std::ostream &out,
2497 const std::map<types::global_dof_index, Point<spacedim>> &support_points)
2498 {
2499 AssertThrow(out.fail() == false, ExcIO());
2500
2501 std::map<Point<spacedim>,
2502 std::vector<types::global_dof_index>,
2504 point_map;
2505
2506 // convert to map point -> list of DoFs
2507 for (const auto &it : support_points)
2508 {
2509 std::vector<types::global_dof_index> &v = point_map[it.second];
2510 v.push_back(it.first);
2511 }
2512
2513 // print the newly created map:
2514 for (const auto &it : point_map)
2515 {
2516 out << it.first << " \"";
2517 const std::vector<types::global_dof_index> &v = it.second;
2518 for (unsigned int i = 0; i < v.size(); ++i)
2519 {
2520 if (i > 0)
2521 out << ", ";
2522 out << v[i];
2523 }
2524
2525 out << "\"\n";
2526 }
2527
2528 out << std::flush;
2529 }
2530
2531
2532 template <int dim, int spacedim>
2533 void
2535 const Table<2, Coupling> &table,
2536 std::vector<Table<2, Coupling>> &tables_by_block)
2537 {
2538 if (dof_handler.has_hp_capabilities() == false)
2539 {
2540 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
2541 const unsigned int nb = fe.n_blocks();
2542
2543 tables_by_block.resize(1);
2544 tables_by_block[0].reinit(nb, nb);
2545 tables_by_block[0].fill(none);
2546
2547 for (unsigned int i = 0; i < fe.n_components(); ++i)
2548 {
2549 const unsigned int ib = fe.component_to_block_index(i);
2550 for (unsigned int j = 0; j < fe.n_components(); ++j)
2551 {
2552 const unsigned int jb = fe.component_to_block_index(j);
2553 tables_by_block[0](ib, jb) |= table(i, j);
2554 }
2555 }
2556 }
2557 else
2558 {
2559 const hp::FECollection<dim> &fe_collection =
2560 dof_handler.get_fe_collection();
2561 tables_by_block.resize(fe_collection.size());
2562
2563 for (unsigned int f = 0; f < fe_collection.size(); ++f)
2564 {
2565 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
2566
2567 const unsigned int nb = fe.n_blocks();
2568 tables_by_block[f].reinit(nb, nb);
2569 tables_by_block[f].fill(none);
2570 for (unsigned int i = 0; i < fe.n_components(); ++i)
2571 {
2572 const unsigned int ib = fe.component_to_block_index(i);
2573 for (unsigned int j = 0; j < fe.n_components(); ++j)
2574 {
2575 const unsigned int jb = fe.component_to_block_index(j);
2576 tables_by_block[f](ib, jb) |= table(i, j);
2577 }
2578 }
2579 }
2580 }
2581 }
2582
2583
2584
2585 template <int dim, int spacedim>
2586 void
2588 const DoFHandler<dim, spacedim> &dof_handler,
2589 const unsigned int level,
2590 const std::vector<bool> &selected_dofs,
2591 const types::global_dof_index offset)
2592 {
2593 std::vector<types::global_dof_index> indices;
2594
2595 unsigned int i = 0;
2596
2597 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2598 if (cell->is_locally_owned_on_level())
2599 ++i;
2600 block_list.reinit(i,
2601 dof_handler.n_dofs(),
2602 dof_handler.get_fe().n_dofs_per_cell());
2603 i = 0;
2604 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2605 if (cell->is_locally_owned_on_level())
2606 {
2607 indices.resize(cell->get_fe().n_dofs_per_cell());
2608 cell->get_mg_dof_indices(indices);
2609
2610 if (selected_dofs.size() != 0)
2611 AssertDimension(indices.size(), selected_dofs.size());
2612
2613 for (types::global_dof_index j = 0; j < indices.size(); ++j)
2614 {
2615 if (selected_dofs.empty())
2616 block_list.add(i, indices[j] - offset);
2617 else
2618 {
2619 if (selected_dofs[j])
2620 block_list.add(i, indices[j] - offset);
2621 }
2622 }
2623 ++i;
2624 }
2625 }
2626
2627
2628 template <int dim, int spacedim>
2629 void
2631 const DoFHandler<dim, spacedim> &dof_handler,
2632 const unsigned int level,
2633 const bool interior_only)
2634 {
2635 const FiniteElement<dim> &fe = dof_handler.get_fe();
2636 block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level));
2637
2638 std::vector<types::global_dof_index> indices;
2639 std::vector<bool> exclude;
2640
2641 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2642 {
2643 indices.resize(cell->get_fe().n_dofs_per_cell());
2644 cell->get_mg_dof_indices(indices);
2645
2646 if (interior_only)
2647 {
2648 // Exclude degrees of freedom on faces opposite to the vertex
2649 exclude.resize(fe.n_dofs_per_cell());
2650 std::fill(exclude.begin(), exclude.end(), false);
2651
2652 for (const unsigned int face : cell->face_indices())
2653 if (cell->at_boundary(face) ||
2654 cell->neighbor(face)->level() != cell->level())
2655 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2656 exclude[fe.face_to_cell_index(i, face)] = true;
2657 for (types::global_dof_index j = 0; j < indices.size(); ++j)
2658 if (!exclude[j])
2659 block_list.add(0, indices[j]);
2660 }
2661 else
2662 {
2663 for (const auto index : indices)
2664 block_list.add(0, index);
2665 }
2666 }
2667 }
2668
2669
2670 template <int dim, int spacedim>
2671 void
2673 const DoFHandler<dim, spacedim> &dof_handler,
2674 const unsigned int level,
2675 const bool interior_dofs_only,
2676 const bool boundary_dofs)
2677 {
2678 Assert(level > 0 && level < dof_handler.get_triangulation().n_levels(),
2679 ExcIndexRange(level, 1, dof_handler.get_triangulation().n_levels()));
2680
2681 std::vector<types::global_dof_index> indices;
2682 std::vector<bool> exclude;
2683
2684 unsigned int block = 0;
2685 for (const auto &pcell : dof_handler.cell_iterators_on_level(level - 1))
2686 {
2687 if (pcell->is_active())
2688 continue;
2689
2690 for (unsigned int child = 0; child < pcell->n_children(); ++child)
2691 {
2692 const auto cell = pcell->child(child);
2693
2694 // For hp, only this line here would have to be replaced.
2695 const FiniteElement<dim> &fe = dof_handler.get_fe();
2696 const unsigned int n_dofs = fe.n_dofs_per_cell();
2697 indices.resize(n_dofs);
2698 exclude.resize(n_dofs);
2699 std::fill(exclude.begin(), exclude.end(), false);
2700 cell->get_mg_dof_indices(indices);
2701
2702 if (interior_dofs_only)
2703 {
2704 // Eliminate dofs on faces of the child which are on faces
2705 // of the parent
2706 for (unsigned int d = 0; d < dim; ++d)
2707 {
2708 const unsigned int face =
2710 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2711 exclude[fe.face_to_cell_index(i, face)] = true;
2712 }
2713
2714 // Now remove all degrees of freedom on the domain boundary
2715 // from the exclusion list
2716 if (boundary_dofs)
2717 for (const unsigned int face :
2719 if (cell->at_boundary(face))
2720 for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
2721 ++i)
2722 exclude[fe.face_to_cell_index(i, face)] = false;
2723 }
2724
2725 for (unsigned int i = 0; i < n_dofs; ++i)
2726 if (!exclude[i])
2727 block_list.add(block, indices[i]);
2728 }
2729 ++block;
2730 }
2731 }
2732
2733 template <int dim, int spacedim>
2734 std::vector<unsigned int>
2736 const DoFHandler<dim, spacedim> &dof_handler,
2737 const unsigned int level,
2738 const bool interior_only,
2739 const bool boundary_patches,
2740 const bool level_boundary_patches,
2741 const bool single_cell_patches,
2742 const bool invert_vertex_mapping)
2743 {
2744 const unsigned int n_blocks = dof_handler.get_fe().n_blocks();
2745 BlockMask exclude_boundary_dofs = BlockMask(n_blocks, interior_only);
2746 return make_vertex_patches(block_list,
2747 dof_handler,
2748 level,
2749 exclude_boundary_dofs,
2750 boundary_patches,
2751 level_boundary_patches,
2752 single_cell_patches,
2753 invert_vertex_mapping);
2754 }
2755
2756 template <int dim, int spacedim>
2757 std::vector<unsigned int>
2759 const DoFHandler<dim, spacedim> &dof_handler,
2760 const unsigned int level,
2761 const BlockMask &exclude_boundary_dofs,
2762 const bool boundary_patches,
2763 const bool level_boundary_patches,
2764 const bool single_cell_patches,
2765 const bool invert_vertex_mapping)
2766 {
2767 // Vector mapping from vertex index in the triangulation to consecutive
2768 // block indices on this level The number of cells at a vertex
2769 std::vector<unsigned int> vertex_cell_count(
2770 dof_handler.get_triangulation().n_vertices(), 0);
2771
2772 // Is a vertex at the boundary?
2773 std::vector<bool> vertex_boundary(
2774 dof_handler.get_triangulation().n_vertices(), false);
2775
2776 std::vector<unsigned int> vertex_mapping(
2777 dof_handler.get_triangulation().n_vertices(),
2779
2780 // Estimate for the number of dofs at this point
2781 std::vector<unsigned int> vertex_dof_count(
2782 dof_handler.get_triangulation().n_vertices(), 0);
2783
2784 // Identify all vertices active on this level and remember some data
2785 // about them
2786 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2787 for (const unsigned int v : cell->vertex_indices())
2788 {
2789 const unsigned int vg = cell->vertex_index(v);
2790 vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell();
2791 ++vertex_cell_count[vg];
2792 for (unsigned int d = 0; d < dim; ++d)
2793 {
2794 const unsigned int face = GeometryInfo<dim>::vertex_to_face[v][d];
2795 if (cell->at_boundary(face))
2796 vertex_boundary[vg] = true;
2797 else if ((!level_boundary_patches) &&
2798 (cell->neighbor(face)->level() !=
2799 static_cast<int>(level)))
2800 vertex_boundary[vg] = true;
2801 }
2802 }
2803 // From now on, only vertices with positive dof count are "in".
2804
2805 // Remove vertices at boundaries or in corners
2806 for (unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg)
2807 if ((!single_cell_patches && vertex_cell_count[vg] < 2) ||
2808 (!boundary_patches && vertex_boundary[vg]))
2809 vertex_dof_count[vg] = 0;
2810
2811 // Create a mapping from all vertices to the ones used here
2812 unsigned int n_vertex_count = 0;
2813 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2814 if (vertex_dof_count[vg] != 0)
2815 vertex_mapping[vg] = n_vertex_count++;
2816
2817 // Compactify dof count
2818 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2819 if (vertex_dof_count[vg] != 0)
2820 vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
2821
2822 // Now that we have all the data, we reduce it to the part we actually
2823 // want
2824 vertex_dof_count.resize(n_vertex_count);
2825
2826 // At this point, the list of patches is ready. Now we enter the dofs
2827 // into the sparsity pattern.
2828 block_list.reinit(vertex_dof_count.size(),
2829 dof_handler.n_dofs(level),
2830 vertex_dof_count);
2831
2832 std::vector<types::global_dof_index> indices;
2833 std::vector<bool> exclude;
2834
2835 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2836 {
2837 const FiniteElement<dim> &fe = cell->get_fe();
2838 indices.resize(fe.n_dofs_per_cell());
2839 cell->get_mg_dof_indices(indices);
2840
2841 for (const unsigned int v : cell->vertex_indices())
2842 {
2843 const unsigned int vg = cell->vertex_index(v);
2844 const unsigned int block = vertex_mapping[vg];
2845 if (block == numbers::invalid_unsigned_int)
2846 continue;
2847
2848 // Collect excluded dofs for some block(s) if boundary dofs
2849 // for a block are decided to be excluded
2850 if (exclude_boundary_dofs.size() == 0 ||
2851 exclude_boundary_dofs.n_selected_blocks() != 0)
2852 {
2853 // Exclude degrees of freedom on faces opposite to the
2854 // vertex
2855 exclude.resize(fe.n_dofs_per_cell());
2856 std::fill(exclude.begin(), exclude.end(), false);
2857
2858 for (unsigned int d = 0; d < dim; ++d)
2859 {
2860 const unsigned int a_face =
2862 const unsigned int face =
2864 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2865 {
2866 // For each dof, get the block it is in and decide to
2867 // exclude it or not
2868 if (exclude_boundary_dofs[fe.system_to_block_index(
2870 i, face))
2871 .first] == true)
2872 exclude[fe.face_to_cell_index(i, face)] = true;
2873 }
2874 }
2875 for (unsigned int j = 0; j < indices.size(); ++j)
2876 if (!exclude[j])
2877 block_list.add(block, indices[j]);
2878 }
2879 else
2880 {
2881 for (const auto index : indices)
2882 block_list.add(block, index);
2883 }
2884 }
2885 }
2886
2887 if (invert_vertex_mapping)
2888 {
2889 // Compress vertex mapping
2890 unsigned int n_vertex_count = 0;
2891 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2892 if (vertex_mapping[vg] != numbers::invalid_unsigned_int)
2893 vertex_mapping[n_vertex_count++] = vg;
2894
2895 // Now we reduce it to the part we actually want
2896 vertex_mapping.resize(n_vertex_count);
2897 }
2898
2899 return vertex_mapping;
2900 }
2901
2902
2903 template <int dim, int spacedim>
2904 unsigned int
2906 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2907 &patch)
2908 {
2909 std::set<types::global_dof_index> dofs_on_patch;
2910 std::vector<types::global_dof_index> local_dof_indices;
2911
2912 // loop over the cells in the patch and get the DoFs on each.
2913 // add all of them to a std::set which automatically makes sure
2914 // all duplicates are ignored
2915 for (unsigned int i = 0; i < patch.size(); ++i)
2916 {
2918 patch[i];
2919 Assert(cell->is_artificial() == false,
2920 ExcMessage("This function can not be called with cells that are "
2921 "not either locally owned or ghost cells."));
2922 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2923 cell->get_dof_indices(local_dof_indices);
2924 dofs_on_patch.insert(local_dof_indices.begin(),
2925 local_dof_indices.end());
2926 }
2927
2928 // now return the number of DoFs (duplicates were ignored)
2929 return dofs_on_patch.size();
2930 }
2931
2932
2933
2934 template <int dim, int spacedim>
2935 std::vector<types::global_dof_index>
2937 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2938 &patch)
2939 {
2940 std::set<types::global_dof_index> dofs_on_patch;
2941 std::vector<types::global_dof_index> local_dof_indices;
2942
2943 // loop over the cells in the patch and get the DoFs on each.
2944 // add all of them to a std::set which automatically makes sure
2945 // all duplicates are ignored
2946 for (unsigned int i = 0; i < patch.size(); ++i)
2947 {
2949 patch[i];
2950 Assert(cell->is_artificial() == false,
2951 ExcMessage("This function can not be called with cells that are "
2952 "not either locally owned or ghost cells."));
2953 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2954 cell->get_dof_indices(local_dof_indices);
2955 dofs_on_patch.insert(local_dof_indices.begin(),
2956 local_dof_indices.end());
2957 }
2958
2959 Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)),
2961
2962 // return a vector with the content of the set above. copying
2963 // also ensures that we retain sortedness as promised in the
2964 // documentation and as necessary to retain the block structure
2965 // also on the local system
2966 return std::vector<types::global_dof_index>(dofs_on_patch.begin(),
2967 dofs_on_patch.end());
2968 }
2969
2970
2971} // end of namespace DoFTools
2972
2973
2974
2975// explicit instantiations
2976
2977#include "dof_tools.inst"
2978
2979
2980
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
size_type n_constraints() const
unsigned int size() const
unsigned int n_selected_blocks(const unsigned int overall_number_of_blocks=numbers::invalid_unsigned_int) const
bool represents_n_components(const unsigned int n) const
bool represents_the_all_selected_mask() const
unsigned int size() const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) 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
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
const IndexSet & locally_owned_dofs() const
active_cell_iterator begin_active(const unsigned int level=0) const
std::vector< types::fe_index > get_active_fe_indices() const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
types::global_dof_index n_locally_owned_dofs() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
const FEValues< dim, spacedim > & get_present_fe_values() const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_blocks() const
unsigned int n_components() const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
const ComponentMask & get_nonzero_components(const unsigned int i) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const unsigned char combined_orientation=ReferenceCell::default_combined_face_orientation()) const
bool is_primitive() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
unsigned int component_to_block_index(const unsigned int component) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
unsigned int element_multiplicity(const unsigned int index) const
unsigned int n_base_elements() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1990
bool is_element(const size_type index) const
Definition index_set.h:1883
void add_index(const size_type index)
Definition index_set.h:1784
void fill_binary_vector(VectorType &vector) const
void subtract_set(const IndexSet &other)
Definition index_set.cc:470
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1971
void compress() const
Definition index_set.h:1773
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1820
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
void add(const size_type i, const size_type j)
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
unsigned int n_levels() const
unsigned int n_vertices() const
virtual size_type size() const override
unsigned int size() const
Definition collection.h:264
unsigned int max_dofs_per_face() const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
unsigned int n_components() const
unsigned int max_dofs_per_cell() const
void push_back(const Quadrature< dim_in > &new_quadrature)
virtual MPI_Comm get_communicator() const override
Definition tria_base.cc:160
#define DEAL_II_DEPRECATED
Definition config.h:207
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
unsigned int level
Definition grid_out.cc:4626
unsigned int cell_index
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< FilteredIterator< BaseIterator > > filter_iterators(IteratorRange< BaseIterator > i, const Predicate &p)
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::line_iterator line_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
typename ActiveSelector::CellAccessor cell_accessor
@ update_quadrature_points
Transformed quadrature points.
void get_block_association(const DoFHandler< dim, spacedim > &dof, std::vector< unsigned char > &dofs_by_block)
Definition dof_tools.cc:275
void extract_constant_modes(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, const unsigned int mg_level, std::vector< std::vector< bool > > &constant_modes)
std::vector< unsigned char > get_local_component_association(const FiniteElement< dim, spacedim > &fe, const ComponentMask &component_mask)
Definition dof_tools.cc:122
void get_component_association(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< unsigned char > &dofs_by_component, const unsigned int mg_level=numbers::invalid_unsigned_int)
Definition dof_tools.cc:193
void resolve_components(const FiniteElement< dim, spacedim > &fe, const std::vector< unsigned char > &dofs_by_component, const std::vector< unsigned int > &target_component, const bool only_once, std::vector< types::global_dof_index > &dofs_per_component, unsigned int &component)
void get_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::subdomain_id > &subdomain)
IndexSet dof_indices_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
IndexSet extract_dofs_with_support_contained_within(const DoFHandler< dim, spacedim > &dof_handler, const std::function< bool(const typename DoFHandler< dim, spacedim >::active_cell_iterator &)> &predicate, const AffineConstraints< number > &constraints={})
Definition dof_tools.cc:824
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_boundary_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask={}, const std::set< types::boundary_id > &boundary_ids={})
Definition dof_tools.cc:614
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
std::map< typename DoFHandler< dim - 1, spacedim >::active_cell_iterator, std::pair< typename DoFHandler< dim, spacedim >::active_cell_iterator, unsigned int > > map_boundary_to_bulk_dof_iterators(const std::map< typename Triangulation< dim - 1, spacedim >::cell_iterator, typename Triangulation< dim, spacedim >::face_iterator > &c1_to_c0, const DoFHandler< dim, spacedim > &c0_dh, const DoFHandler< dim - 1, spacedim > &c1_dh)
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
void make_cell_patches(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const std::vector< bool > &selected_dofs={}, const types::global_dof_index offset=0)
IndexSet extract_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask)
Definition dof_tools.cc:414
IndexSet extract_locally_active_dofs(const DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
void map_dofs_to_support_points(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, std::vector< Point< spacedim > > &support_points, const ComponentMask &mask={})
void write_gnuplot_dof_support_point_info(std::ostream &out, const std::map< types::global_dof_index, Point< spacedim > > &support_points)
void extract_subdomain_dofs(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain_id, std::vector< bool > &selected_dofs)
IndexSet extract_locally_active_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
std::vector< unsigned int > make_vertex_patches(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool interior_dofs_only, const bool boundary_patches=false, const bool level_boundary_patches=false, const bool single_cell_patches=false, const bool invert_vertex_mapping=false)
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
void extract_level_dofs(const unsigned int level, const DoFHandler< dim, spacedim > &dof, const ComponentMask &component_mask, std::vector< bool > &selected_dofs)
Definition dof_tools.cc:501
void extract_dofs_with_support_on_boundary(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
Definition dof_tools.cc:742
void distribute_cell_to_dof_vector(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &cell_data, Vector< double > &dof_data, const unsigned int component=0)
Definition dof_tools.cc:330
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
void make_child_patches(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool interior_dofs_only, const bool boundary_dofs=false)
void map_dof_to_boundary_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::global_dof_index > &mapping)
std::vector< IndexSet > locally_relevant_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
void make_single_patch(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool interior_dofs_only=false)
std::vector< types::global_dof_index > get_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
std::vector< IndexSet > locally_owned_dofs_per_component(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &components={})
Definition dof_tools.cc:467
unsigned int count_dofs_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
unsigned int count_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
void convert_couplings_to_blocks(const DoFHandler< dim, spacedim > &dof_handler, const Table< 2, Coupling > &table_by_component, std::vector< Table< 2, Coupling > > &tables_by_block)
void extract_constant_modes(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< std::vector< bool > > &constant_modes)
IndexSet extract_hanging_node_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void extract_level_constant_modes(const unsigned int level, const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< std::vector< bool > > &constant_modes)
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
const types::boundary_id internal_face_boundary_id
Definition types.h:312
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
const types::subdomain_id invalid_subdomain_id
Definition types.h:341
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::global_dof_index invalid_dof_index
Definition types.h:252
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:59
bool operator()(const Point< dim, Number > &lhs, const Point< dim, Number > &rhs) const
Definition dof_tools.cc:80
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition types.h:102