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