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