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