Reference documentation for deal.II version 9.3.3
dof_tools.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2021 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
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>
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
149 {
150 const unsigned int first_comp =
152
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)
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
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] =
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 : dof.active_cell_iterators())
228 if (c->is_locally_owned())
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 {
311 const Triangulation<dim, spacedim> &tria = dof_handler.get_triangulation();
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,
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
351 cell = dof_handler.begin_active(),
352 endc = dof_handler.end();
353 std::vector<types::global_dof_index> dof_indices;
354 dof_indices.reserve(fe_collection.max_dofs_per_cell());
355
356 for (unsigned int present_cell = 0; cell != endc; ++cell, ++present_cell)
357 {
358 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
359 dof_indices.resize(dofs_per_cell);
360 cell->get_dof_indices(dof_indices);
361
362 for (unsigned int i = 0; i < dofs_per_cell; ++i)
363 // consider this dof only if it is the right component. if there
364 // is only one component, short cut the test
365 if (!consider_components ||
366 (cell->get_fe().system_to_component_index(i).first == component))
367 {
368 // sum up contribution of the present_cell to this dof
369 dof_data(dof_indices[i]) += cell_data(present_cell);
370
371 // note that we added another summand
372 ++touch_count[dof_indices[i]];
373 }
374 }
375
376 // compute the mean value on all the dofs by dividing with the number
377 // of summands.
378 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
379 {
380 // assert that each dof was used at least once. this needs not be
381 // the case if the vector has more than one component
382 Assert(consider_components || (touch_count[i] != 0),
384 if (touch_count[i] != 0)
385 dof_data(i) /= touch_count[i];
386 }
387 }
388
389
390
391 template <int dim, int spacedim>
392 void
395 std::vector<bool> & selected_dofs)
396 {
400 "The given component mask is not sized correctly to represent the "
401 "components of the given finite element."));
402 Assert(selected_dofs.size() == dof.n_locally_owned_dofs(),
403 ExcDimensionMismatch(selected_dofs.size(),
404 dof.n_locally_owned_dofs()));
405
406 // two special cases: no component is selected, and all components are
407 // selected; both rather stupid, but easy to catch
409 dof.get_fe_collection().n_components()) == 0)
410 {
411 std::fill_n(selected_dofs.begin(), dof.n_locally_owned_dofs(), false);
412 return;
413 }
417 {
418 std::fill_n(selected_dofs.begin(), dof.n_locally_owned_dofs(), true);
419 return;
420 }
421
422
423 // preset all values by false
424 std::fill_n(selected_dofs.begin(), dof.n_locally_owned_dofs(), false);
425
426 // get the component association of each DoF and then select the ones
427 // that match the given set of components
428 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs());
430
431 for (types::global_dof_index i = 0; i < dof.n_locally_owned_dofs(); ++i)
433 selected_dofs[i] = true;
434 }
435
436
437
438 template <int dim, int spacedim>
442 {
446 "The given component mask is not sized correctly to represent the "
447 "components of the given finite element."));
448
449 // Two special cases: no component is selected, and all components are
450 // selected; both rather stupid, but easy to catch
452 dof.get_fe_collection().n_components()) == 0)
453 return IndexSet(dof.n_dofs());
457 return dof.locally_owned_dofs();
458
459 // get the component association of each DoF and then select the ones
460 // that match the given set of components
461 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs());
463
464 // fill the selected components in a vector
465 std::vector<types::global_dof_index> selected_dofs;
466 selected_dofs.reserve(dof.n_locally_owned_dofs());
467 for (types::global_dof_index i = 0; i < dofs_by_component.size(); ++i)
469 selected_dofs.push_back(dof.locally_owned_dofs().nth_index_in_set(i));
470
471 // fill vector of indices to return argument
472 IndexSet result(dof.n_dofs());
474 return result;
475 }
476
477
478
479 template <int dim, int spacedim>
480 void
483 std::vector<bool> & selected_dofs)
484 {
485 // simply forward to the function that works based on a component mask
486 extract_dofs<dim, spacedim>(
488 }
489
490
491
492 template <int dim, int spacedim>
496 {
497 // simply forward to the function that works based on a component mask
498 return extract_dofs<dim, spacedim>(
500 }
501
502
503
504 template <int dim, int spacedim>
505 std::vector<IndexSet>
508 {
509 const auto n_comps = dof.get_fe_collection().n_components();
512 "The given component mask is not sized correctly to represent the "
513 "components of the given finite element."));
514
515 const auto &locally_owned_dofs = dof.locally_owned_dofs();
516
517 // get the component association of each DoF and then select the ones
518 // that match the given set of components
519 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs());
521
522 std::vector<IndexSet> index_per_comp(n_comps, IndexSet(dof.n_dofs()));
523
524 for (types::global_dof_index i = 0; i < dof.n_locally_owned_dofs(); ++i)
525 {
526 const auto &comp_i = dofs_by_component[i];
529 locally_owned_dofs.nth_index_in_set(i));
530 }
531 for (auto &c : index_per_comp)
532 c.compress();
533 return index_per_comp;
534 }
535
536
537
538 template <int dim, int spacedim>
539 void
540 extract_level_dofs(const unsigned int level,
541 const DoFHandler<dim, spacedim> &dof,
543 std::vector<bool> & selected_dofs)
544 {
545 const FiniteElement<dim, spacedim> &fe = dof.get_fe();
546
550 "The given component mask is not sized correctly to represent the "
551 "components of the given finite element."));
552 Assert(selected_dofs.size() == dof.n_dofs(level),
553 ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
554
555 // two special cases: no component is selected, and all components are
556 // selected, both rather stupid, but easy to catch
558 dof.get_fe_collection().n_components()) == 0)
559 {
560 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false);
561 return;
562 }
566 {
567 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), true);
568 return;
569 }
570
571 // preset all values by false
572 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false);
573
574 // next set up a table for the degrees of freedom on each of the cells
575 // whether it is something interesting or not
576 std::vector<unsigned char> local_component_asssociation =
578 std::vector<bool> local_selected_dofs(fe.n_dofs_per_cell());
579 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
581
582 // then loop over all cells and do work
583 std::vector<types::global_dof_index> indices(fe.n_dofs_per_cell());
585 for (c = dof.begin(level); c != dof.end(level); ++c)
586 {
587 c->get_mg_dof_indices(indices);
588 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
589 selected_dofs[indices[i]] = local_selected_dofs[i];
590 }
591 }
592
593
594
595 template <int dim, int spacedim>
596 void
597 extract_level_dofs(const unsigned int level,
598 const DoFHandler<dim, spacedim> &dof,
600 std::vector<bool> & selected_dofs)
601 {
602 // simply defer to the other extract_level_dofs() function
604 dof,
606 selected_dofs);
607 }
608
609
610
611 template <int dim, int spacedim>
612 void
615 std::vector<bool> & selected_dofs,
616 const std::set<types::boundary_id> &boundary_ids)
617 {
618 Assert((dynamic_cast<
620 &dof_handler.get_triangulation()) == nullptr),
622 "This function can not be used with distributed triangulations. "
624
625 IndexSet indices;
627
628 // clear and reset array by default values
629 selected_dofs.clear();
630 selected_dofs.resize(dof_handler.n_dofs(), false);
631
632 // then convert the values computed above to the binary vector
633 indices.fill_binary_vector(selected_dofs);
634 }
635
636
637
638 template <int dim, int spacedim>
639 void
642 IndexSet & selected_dofs,
643 const std::set<types::boundary_id> &boundary_ids)
644 {
645 // Simply forward to the other function
646 selected_dofs =
648 }
649
650
651
652 template <int dim, int spacedim>
656 const std::set<types::boundary_id> &boundary_ids)
657 {
659 dof_handler.get_fe_collection().n_components()),
660 ExcMessage("Component mask has invalid size."));
661 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
662 boundary_ids.end(),
664
665 IndexSet selected_dofs(dof_handler.n_dofs());
666
667 // let's see whether we have to check for certain boundary indicators
668 // or whether we can accept all
669 const bool check_boundary_id = (boundary_ids.size() != 0);
670
671 // also see whether we have to check whether a certain vector component
672 // is selected, or all
673 const bool check_vector_component =
676 dof_handler.get_fe_collection().n_components()) !=
677 dof_handler.get_fe_collection().n_components()));
678
679 std::vector<types::global_dof_index> face_dof_indices;
680 face_dof_indices.reserve(
681 dof_handler.get_fe_collection().max_dofs_per_face());
682
683 // now loop over all cells and check whether their faces are at the
684 // boundary. note that we need not take special care of single lines
685 // being at the boundary (using @p{cell->has_boundary_lines}), since we
686 // do not support boundaries of dimension dim-2, and so every isolated
687 // boundary line is also part of a boundary face which we will be
688 // visiting sooner or later
689 for (const auto &cell : dof_handler.active_cell_iterators())
690 // only work on cells that are either locally owned or at least ghost
691 // cells
692 if (cell->is_artificial() == false)
693 for (const unsigned int face : cell->face_indices())
694 if (cell->at_boundary(face))
695 if (!check_boundary_id ||
696 (boundary_ids.find(cell->face(face)->boundary_id()) !=
697 boundary_ids.end()))
698 {
699 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
700
701 const auto reference_cell = cell->reference_cell();
702
703 const unsigned int n_vertices_per_cell =
704 reference_cell.n_vertices();
705 const unsigned int n_lines_per_cell = reference_cell.n_lines();
706 const unsigned int n_vertices_per_face =
707 reference_cell.face_reference_cell(face).n_vertices();
708 const unsigned int n_lines_per_face =
709 reference_cell.face_reference_cell(face).n_lines();
710
711 const unsigned int dofs_per_face = fe.n_dofs_per_face(face);
712 face_dof_indices.resize(dofs_per_face);
713 cell->face(face)->get_dof_indices(face_dof_indices,
714 cell->active_fe_index());
715
716 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
717 if (!check_vector_component)
719 else
720 // check for component is required. somewhat tricky as
721 // usual for the case that the shape function is
722 // non-primitive, but use usual convention (see docs)
723 {
724 // first get at the cell-global number of a face dof,
725 // to ask the FE certain questions
726 const unsigned int cell_index =
727 (dim == 1 ?
728 i :
729 (dim == 2 ?
730 (i < 2 * fe.n_dofs_per_vertex() ?
731 i :
732 i + 2 * fe.n_dofs_per_vertex()) :
733 (dim == 3 ? (i < n_vertices_per_face *
734 fe.n_dofs_per_vertex() ?
735 i :
736 (i < n_vertices_per_face *
737 fe.n_dofs_per_vertex() +
738 n_lines_per_face *
739 fe.n_dofs_per_line() ?
740 (i - n_vertices_per_face *
741 fe.n_dofs_per_vertex()) +
742 n_vertices_per_cell *
743 fe.n_dofs_per_vertex() :
744 (i -
745 n_vertices_per_face *
746 fe.n_dofs_per_vertex() -
747 n_lines_per_face *
748 fe.n_dofs_per_line()) +
749 n_vertices_per_cell *
750 fe.n_dofs_per_vertex() +
751 n_lines_per_cell *
752 fe.n_dofs_per_line())) :
754 if (fe.is_primitive(cell_index))
755 {
757 i, face)
758 .first] == true)
760 }
761 else // not primitive
762 {
763 const unsigned int first_nonzero_comp =
766 Assert(first_nonzero_comp < fe.n_components(),
768
771 }
772 }
773 }
774
775 return selected_dofs;
776 }
777
778
779
780 template <int dim, int spacedim>
781 void
783 const DoFHandler<dim, spacedim> & dof_handler,
785 std::vector<bool> & selected_dofs,
786 const std::set<types::boundary_id> &boundary_ids)
787 {
789 dof_handler.get_fe_collection().n_components()),
790 ExcMessage("This component mask has the wrong size."));
791 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
792 boundary_ids.end(),
794
795 // let's see whether we have to check for certain boundary indicators
796 // or whether we can accept all
797 const bool check_boundary_id = (boundary_ids.size() != 0);
798
799 // also see whether we have to check whether a certain vector component
800 // is selected, or all
801 const bool check_vector_component =
803
804 // clear and reset array by default values
805 selected_dofs.clear();
806 selected_dofs.resize(dof_handler.n_dofs(), false);
807 std::vector<types::global_dof_index> cell_dof_indices;
808 cell_dof_indices.reserve(
809 dof_handler.get_fe_collection().max_dofs_per_cell());
810
811 // now loop over all cells and check whether their faces are at the
812 // boundary. note that we need not take special care of single lines
813 // being at the boundary (using @p{cell->has_boundary_lines}), since we
814 // do not support boundaries of dimension dim-2, and so every isolated
815 // boundary line is also part of a boundary face which we will be
816 // visiting sooner or later
817 for (const auto &cell : dof_handler.active_cell_iterators())
818 for (const unsigned int face : cell->face_indices())
819 if (cell->at_boundary(face))
820 if (!check_boundary_id ||
821 (boundary_ids.find(cell->face(face)->boundary_id()) !=
822 boundary_ids.end()))
823 {
824 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
825
826 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
827 cell_dof_indices.resize(dofs_per_cell);
828 cell->get_dof_indices(cell_dof_indices);
829
830 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
831 if (fe.has_support_on_face(i, face))
832 {
833 if (!check_vector_component)
834 selected_dofs[cell_dof_indices[i]] = true;
835 else
836 // check for component is required. somewhat tricky
837 // as usual for the case that the shape function is
838 // non-primitive, but use usual convention (see docs)
839 {
840 if (fe.is_primitive(i))
841 selected_dofs[cell_dof_indices[i]] =
843 .first] == true);
844 else // not primitive
845 {
846 const unsigned int first_nonzero_comp =
849 Assert(first_nonzero_comp < fe.n_components(),
851
852 selected_dofs[cell_dof_indices[i]] =
854 }
855 }
856 }
857 }
858 }
859
860
861
862 template <int dim, int spacedim, typename number>
865 const DoFHandler<dim, spacedim> &dof_handler,
866 const std::function<
868 & predicate,
870 {
871 const std::function<bool(
873 predicate_local =
874 [=](
876 -> bool { return cell->is_locally_owned() && predicate(cell); };
877
878 std::vector<types::global_dof_index> local_dof_indices;
879 local_dof_indices.reserve(
880 dof_handler.get_fe_collection().max_dofs_per_cell());
881
882 // Get all the dofs that live on the subdomain:
883 std::set<types::global_dof_index> predicate_dofs;
884
886 cell = dof_handler.begin_active(),
887 endc = dof_handler.end();
888 for (; cell != endc; ++cell)
889 if (!cell->is_artificial() && predicate(cell))
890 {
891 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
892 cell->get_dof_indices(local_dof_indices);
893 predicate_dofs.insert(local_dof_indices.begin(),
894 local_dof_indices.end());
895 }
896
897 // Get halo layer and accumulate its DoFs
898 std::set<types::global_dof_index> dofs_with_support_on_halo_cells;
899
900 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
901 halo_cells =
902 GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local);
903 for (typename std::vector<typename DoFHandler<dim, spacedim>::
904 active_cell_iterator>::const_iterator it =
905 halo_cells.begin();
906 it != halo_cells.end();
907 ++it)
908 {
909 // skip halo cells that satisfy the given predicate and thereby
910 // shall be a part of the index set on another MPI core.
911 // Those could only be ghost cells with p::d::Tria.
912 if (predicate(*it))
913 {
914 Assert((*it)->is_ghost(), ExcInternalError());
915 continue;
916 }
917
918 const unsigned int dofs_per_cell = (*it)->get_fe().n_dofs_per_cell();
919 local_dof_indices.resize(dofs_per_cell);
920 (*it)->get_dof_indices(local_dof_indices);
921 dofs_with_support_on_halo_cells.insert(local_dof_indices.begin(),
922 local_dof_indices.end());
923 }
924
925 // A complication coming from the fact that DoFs living on locally
926 // owned cells which satisfy predicate may participate in constraints for
927 // DoFs outside of this region.
928 if (cm.n_constraints() > 0)
929 {
930 // Remove DoFs that are part of constraints for DoFs outside
931 // of predicate. Since we will subtract halo_dofs from predicate_dofs,
932 // achieve this by extending halo_dofs with DoFs to which
933 // halo_dofs are constrained.
934 std::set<types::global_dof_index> extra_halo;
935 for (const auto dof : dofs_with_support_on_halo_cells)
936 // if halo DoF is constrained, add all DoFs to which it's constrained
937 // because after resolving constraints, the support of the DoFs that
938 // constrain the current DoF will extend to the halo cells.
939 if (const auto *line_ptr = cm.get_constraint_entries(dof))
940 {
941 const unsigned int line_size = line_ptr->size();
942 for (unsigned int j = 0; j < line_size; ++j)
943 extra_halo.insert((*line_ptr)[j].first);
944 }
945
946 dofs_with_support_on_halo_cells.insert(extra_halo.begin(),
947 extra_halo.end());
948 }
949
950 // Rework our std::set's into IndexSet and subtract halo DoFs from the
951 // predicate_dofs:
952 IndexSet support_set(dof_handler.n_dofs());
954 support_set.compress();
955
956 IndexSet halo_set(dof_handler.n_dofs());
958 dofs_with_support_on_halo_cells.end());
959 halo_set.compress();
960
961 support_set.subtract_set(halo_set);
962
963 // we intentionally do not want to limit the output to locally owned DoFs.
964 return support_set;
965 }
966
967
968
969 namespace internal
970 {
971 namespace
972 {
973 template <int spacedim>
976 const ::DoFHandler<1, spacedim> &dof_handler)
977 {
978 // there are no hanging nodes in 1d
979 return IndexSet(dof_handler.n_dofs());
980 }
981
982
983 template <int spacedim>
986 const ::DoFHandler<2, spacedim> &dof_handler)
987 {
988 const unsigned int dim = 2;
989
990 IndexSet selected_dofs(dof_handler.n_dofs());
991
992 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
993
994 // this function is similar to the make_sparsity_pattern function,
996 for (const auto &cell : dof_handler.active_cell_iterators())
997 if (!cell->is_artificial())
998 {
999 for (const unsigned int face : cell->face_indices())
1000 if (cell->face(face)->has_children())
1001 {
1002 const typename ::DoFHandler<dim,
1003 spacedim>::line_iterator
1004 line = cell->face(face);
1005
1006 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
1007 ++dof)
1009 line->child(0)->vertex_dof_index(1, dof));
1010
1011 for (unsigned int child = 0; child < 2; ++child)
1012 {
1013 if (cell->neighbor_child_on_subface(face, child)
1014 ->is_artificial())
1015 continue;
1016 for (unsigned int dof = 0; dof != fe.n_dofs_per_line();
1017 ++dof)
1019 line->child(child)->dof_index(dof));
1020 }
1021 }
1022 }
1023
1024 selected_dofs.compress();
1025 return selected_dofs;
1026 }
1027
1028
1029 template <int spacedim>
1030 IndexSet
1032 const ::DoFHandler<3, spacedim> &dof_handler)
1033 {
1034 const unsigned int dim = 3;
1035
1036 IndexSet selected_dofs(dof_handler.n_dofs());
1037 IndexSet unconstrained_dofs(dof_handler.n_dofs());
1038
1039 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
1040
1041 for (const auto &cell : dof_handler.active_cell_iterators())
1042 if (!cell->is_artificial())
1043 for (auto f : cell->face_indices())
1044 {
1045 const typename ::DoFHandler<dim, spacedim>::face_iterator
1046 face = cell->face(f);
1047 if (cell->face(f)->has_children())
1048 {
1049 for (unsigned int child = 0; child < 4; ++child)
1050 if (!cell->neighbor_child_on_subface(f, child)
1051 ->is_artificial())
1052 {
1053 // simply take all DoFs that live on this subface
1054 std::vector<types::global_dof_index> ldi(
1055 fe.n_dofs_per_face(f, child));
1056 face->child(child)->get_dof_indices(ldi);
1058 }
1059
1060 // and subtract (in the end) all the indices which a shared
1061 // between this face and its subfaces
1062 for (unsigned int vertex = 0; vertex < 4; ++vertex)
1063 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
1064 ++dof)
1066 face->vertex_dof_index(vertex, dof));
1067 }
1068 }
1069 selected_dofs.subtract_set(unconstrained_dofs);
1070 return selected_dofs;
1071 }
1072 } // namespace
1073 } // namespace internal
1074
1075
1076
1077 template <int dim, int spacedim>
1078 void
1080 std::vector<bool> & selected_dofs)
1081 {
1082 const IndexSet selected_dofs_as_index_set =
1083 extract_hanging_node_dofs(dof_handler);
1084 Assert(selected_dofs.size() == dof_handler.n_dofs(),
1085 ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
1086 // preset all values by false
1087 std::fill(selected_dofs.begin(), selected_dofs.end(), false);
1088 for (const auto index : selected_dofs_as_index_set)
1089 selected_dofs[index] = true;
1090 }
1091
1092
1093
1094 template <int dim, int spacedim>
1095 IndexSet
1097 {
1098 return internal::extract_hanging_node_dofs(dof_handler);
1099 }
1100
1101
1102
1103 template <int dim, int spacedim>
1104 void
1107 std::vector<bool> & selected_dofs)
1108 {
1109 Assert(selected_dofs.size() == dof_handler.n_dofs(),
1110 ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
1111
1112 // preset all values by false
1113 std::fill_n(selected_dofs.begin(), dof_handler.n_dofs(), false);
1114
1115 std::vector<types::global_dof_index> local_dof_indices;
1116 local_dof_indices.reserve(
1117 dof_handler.get_fe_collection().max_dofs_per_cell());
1118
1119 // this function is similar to the make_sparsity_pattern function, see
1122 cell = dof_handler.begin_active(),
1123 endc = dof_handler.end();
1124 for (; cell != endc; ++cell)
1125 if (cell->subdomain_id() == subdomain_id)
1126 {
1127 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1128 local_dof_indices.resize(dofs_per_cell);
1129 cell->get_dof_indices(local_dof_indices);
1130 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1131 selected_dofs[local_dof_indices[i]] = true;
1132 }
1133 }
1134
1135
1136
1137 template <int dim, int spacedim>
1138 void
1140 IndexSet & dof_set)
1141 {
1142 // collect all the locally owned dofs
1143 dof_set = dof_handler.locally_owned_dofs();
1144
1145 // add the DoF on the adjacent ghost cells to the IndexSet, cache them
1146 // in a set. need to check each dof manually because we can't be sure
1147 // that the dof range of locally_owned_dofs is really contiguous.
1148 std::vector<types::global_dof_index> dof_indices;
1149 std::set<types::global_dof_index> global_dof_indices;
1150
1152 cell = dof_handler.begin_active(),
1153 endc = dof_handler.end();
1154 for (; cell != endc; ++cell)
1155 if (cell->is_locally_owned())
1156 {
1157 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1158 cell->get_dof_indices(dof_indices);
1159
1160 for (const types::global_dof_index dof_index : dof_indices)
1161 if (!dof_set.is_element(dof_index))
1162 global_dof_indices.insert(dof_index);
1163 }
1164
1166
1167 dof_set.compress();
1168 }
1169
1170
1171
1172 template <int dim, int spacedim>
1173 void
1175 const DoFHandler<dim, spacedim> &dof_handler,
1176 IndexSet & dof_set,
1177 const unsigned int level)
1178 {
1179 // collect all the locally owned dofs
1180 dof_set = dof_handler.locally_owned_mg_dofs(level);
1181
1182 // add the DoF on the adjacent ghost cells to the IndexSet, cache them
1183 // in a set. need to check each dof manually because we can't be sure
1184 // that the dof range of locally_owned_dofs is really contiguous.
1185 std::vector<types::global_dof_index> dof_indices;
1186 std::set<types::global_dof_index> global_dof_indices;
1187
1188 const auto filtered_iterators_range =
1191 for (const auto &cell : filtered_iterators_range)
1192 {
1193 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1194 cell->get_mg_dof_indices(dof_indices);
1195
1196 for (const types::global_dof_index dof_index : dof_indices)
1197 if (!dof_set.is_element(dof_index))
1198 global_dof_indices.insert(dof_index);
1199 }
1200
1202
1203 dof_set.compress();
1204 }
1205
1206
1207
1208 template <int dim, int spacedim>
1209 void
1211 IndexSet & dof_set)
1212 {
1213 // collect all the locally owned dofs
1214 dof_set = dof_handler.locally_owned_dofs();
1215
1216 // now add the DoF on the adjacent ghost cells to the IndexSet
1217
1218 // Note: For certain meshes (in particular in 3D and with many
1219 // processors), it is really necessary to cache intermediate data. After
1220 // trying several objects such as std::set, a vector that is always kept
1221 // sorted, and a vector that is initially unsorted and sorted once at the
1222 // end, the latter has been identified to provide the best performance.
1223 // Martin Kronbichler
1224 std::vector<types::global_dof_index> dof_indices;
1225 std::vector<types::global_dof_index> dofs_on_ghosts;
1226
1228 cell = dof_handler.begin_active(),
1229 endc = dof_handler.end();
1230 for (; cell != endc; ++cell)
1231 if (cell->is_ghost())
1232 {
1233 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1234 cell->get_dof_indices(dof_indices);
1235 for (const auto dof_index : dof_indices)
1236 if (!dof_set.is_element(dof_index))
1237 dofs_on_ghosts.push_back(dof_index);
1238 }
1239
1240 // sort, compress out duplicates, fill into index set
1241 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1243 std::unique(dofs_on_ghosts.begin(),
1244 dofs_on_ghosts.end()));
1245 dof_set.compress();
1246 }
1247
1248
1249
1250 template <int dim, int spacedim>
1251 void
1253 const DoFHandler<dim, spacedim> &dof_handler,
1254 const unsigned int level,
1255 IndexSet & dof_set)
1256 {
1257 // collect all the locally owned dofs
1258 dof_set = dof_handler.locally_owned_mg_dofs(level);
1259
1260 // add the DoF on the adjacent ghost cells to the IndexSet
1261
1262 // Note: For certain meshes (in particular in 3D and with many
1263 // processors), it is really necessary to cache intermediate data. After
1264 // trying several objects such as std::set, a vector that is always kept
1265 // sorted, and a vector that is initially unsorted and sorted once at the
1266 // end, the latter has been identified to provide the best performance.
1267 // Martin Kronbichler
1268 std::vector<types::global_dof_index> dof_indices;
1269 std::vector<types::global_dof_index> dofs_on_ghosts;
1270
1271 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof_handler.begin(
1272 level),
1273 endc =
1274 dof_handler.end(level);
1275 for (; cell != endc; ++cell)
1276 {
1277 const types::subdomain_id id = cell->level_subdomain_id();
1278
1279 // skip artificial and own cells (only look at ghost cells)
1280 if (id == dof_handler.get_triangulation().locally_owned_subdomain() ||
1282 continue;
1283
1284 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1285 cell->get_mg_dof_indices(dof_indices);
1286 for (const auto dof_index : dof_indices)
1287 if (!dof_set.is_element(dof_index))
1288 dofs_on_ghosts.push_back(dof_index);
1289 }
1290
1291 // sort, compress out duplicates, fill into index set
1292 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1294 std::unique(dofs_on_ghosts.begin(),
1295 dofs_on_ghosts.end()));
1296
1297 dof_set.compress();
1298 }
1299
1300
1301
1302 template <int dim, int spacedim>
1303 void
1306 std::vector<std::vector<bool>> & constant_modes)
1307 {
1308 // If there are no locally owned DoFs, return with an empty
1309 // constant_modes object:
1310 if (dof_handler.n_locally_owned_dofs() == 0)
1311 {
1312 constant_modes = std::vector<std::vector<bool>>(0);
1313 return;
1314 }
1315
1316 const unsigned int n_components = dof_handler.get_fe(0).n_components();
1319
1320 std::vector<unsigned char> dofs_by_component(
1321 dof_handler.n_locally_owned_dofs());
1324 dofs_by_component);
1325 unsigned int n_selected_dofs = 0;
1326 for (unsigned int i = 0; i < n_components; ++i)
1328 n_selected_dofs +=
1329 std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
1330
1331 // Find local numbering within the selected components
1332 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
1333 std::vector<unsigned int> component_numbering(
1334 locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int);
1335 for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements();
1336 ++i)
1338 component_numbering[i] = count++;
1339
1340 // get the element constant modes and find a translation table between
1341 // index in the constant modes and the components.
1342 //
1343 // TODO: We might be able to extend this also for elements which do not
1344 // have the same constant modes, but that is messy...
1345 const ::hp::FECollection<dim, spacedim> &fe_collection =
1346 dof_handler.get_fe_collection();
1347 std::vector<Table<2, bool>> element_constant_modes;
1348 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
1349 constant_mode_to_component_translation(n_components);
1350 unsigned int n_constant_modes = 0;
1351 for (unsigned int f = 0; f < fe_collection.size(); ++f)
1352 {
1353 std::pair<Table<2, bool>, std::vector<unsigned int>> data =
1354 fe_collection[f].get_constant_modes();
1355 element_constant_modes.push_back(data.first);
1356 if (f == 0)
1357 for (unsigned int i = 0; i < data.second.size(); ++i)
1359 constant_mode_to_component_translation[data.second[i]]
1360 .emplace_back(n_constant_modes++, i);
1361 AssertDimension(element_constant_modes.back().n_rows(),
1362 element_constant_modes[0].n_rows());
1363 }
1364
1365 // First count the number of dofs in the current component.
1366 constant_modes.clear();
1367 constant_modes.resize(n_constant_modes,
1368 std::vector<bool>(n_selected_dofs, false));
1369
1370 // Loop over all owned cells and ask the element for the constant modes
1371 std::vector<types::global_dof_index> dof_indices;
1372 for (const auto &cell : dof_handler.active_cell_iterators())
1373 if (cell->is_locally_owned())
1374 {
1375 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1376 cell->get_dof_indices(dof_indices);
1377
1378 for (unsigned int i = 0; i < dof_indices.size(); ++i)
1379 if (locally_owned_dofs.is_element(dof_indices[i]))
1380 {
1381 const unsigned int loc_index =
1382 locally_owned_dofs.index_within_set(dof_indices[i]);
1383 const unsigned int comp = dofs_by_component[loc_index];
1385 for (auto &indices :
1386 constant_mode_to_component_translation[comp])
1387 constant_modes[indices
1388 .first][component_numbering[loc_index]] =
1389 element_constant_modes[cell->active_fe_index()](
1390 indices.second, i);
1391 }
1392 }
1393 }
1394
1395
1396
1397 template <int dim, int spacedim>
1398 void
1400 std::vector<unsigned int> & active_fe_indices)
1401 {
1402 AssertDimension(active_fe_indices.size(),
1403 dof_handler.get_triangulation().n_active_cells());
1404
1406 cell = dof_handler.begin_active(),
1407 endc = dof_handler.end();
1408 for (; cell != endc; ++cell)
1409 active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
1410 }
1411
1412 template <int dim, int spacedim>
1413 std::vector<IndexSet>
1415 {
1416 Assert(dof_handler.n_dofs() > 0,
1417 ExcMessage("The given DoFHandler has no DoFs."));
1418
1419 // If the Triangulation is distributed, the only thing we can usefully
1420 // ask is for its locally owned subdomain
1421 Assert((dynamic_cast<
1423 &dof_handler.get_triangulation()) == nullptr),
1424 ExcMessage(
1425 "For parallel::distributed::Triangulation objects and "
1426 "associated DoF handler objects, asking for any information "
1427 "related to a subdomain other than the locally owned one does "
1428 "not make sense."));
1429
1430 // The following is a random process (flip of a coin), thus should be called
1431 // once only.
1432 std::vector<::types::subdomain_id> subdomain_association(
1433 dof_handler.n_dofs());
1435 subdomain_association);
1436
1437 // Figure out how many subdomain ids there are.
1438 //
1439 // if this is a parallel triangulation, then we can just ask the
1440 // triangulation for this. if this is a sequential triangulation, we loop
1441 // over all cells and take the largest subdomain_id value we find; the
1442 // number of subdomains is then the largest found value plus one. (we here
1443 // assume that all subdomain ids up to the largest are actually used; this
1444 // may not be true for a sequential triangulation where these values have
1445 // been set by hand and not in accordance with some MPI communicator; but
1446 // the function returns an array indexed starting at zero, so we need to
1447 // collect information for each subdomain index anyway, not just for the
1448 // used one.)
1449 const unsigned int n_subdomains =
1450 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1451 &dof_handler.get_triangulation()) == nullptr ?
1452 [&dof_handler]() {
1453 unsigned int max_subdomain_id = 0;
1454 for (const auto &cell : dof_handler.active_cell_iterators())
1455 max_subdomain_id =
1456 std::max(max_subdomain_id, cell->subdomain_id());
1457 return max_subdomain_id + 1;
1458 }() :
1460 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1461 &dof_handler.get_triangulation())
1462 ->get_communicator()));
1463 Assert(n_subdomains > *std::max_element(subdomain_association.begin(),
1464 subdomain_association.end()),
1466
1467 std::vector<::IndexSet> index_sets(
1468 n_subdomains, ::IndexSet(dof_handler.n_dofs()));
1469
1470 // loop over subdomain_association and populate IndexSet when a
1471 // change in subdomain ID is found
1472 ::types::global_dof_index i_min = 0;
1473 ::types::global_dof_index this_subdomain = subdomain_association[0];
1474
1475 for (::types::global_dof_index index = 1;
1476 index < subdomain_association.size();
1477 ++index)
1478 {
1479 // found index different from the current one
1480 if (subdomain_association[index] != this_subdomain)
1481 {
1483 i_min = index;
1484 this_subdomain = subdomain_association[index];
1485 }
1486 }
1487
1488 // the very last element is of different index
1489 if (i_min == subdomain_association.size() - 1)
1490 {
1492 }
1493
1494 // otherwise there are at least two different indices
1495 else
1496 {
1498 subdomain_association.size());
1499 }
1500
1501 for (unsigned int i = 0; i < n_subdomains; i++)
1502 index_sets[i].compress();
1503
1504 return index_sets;
1505 }
1506
1507 template <int dim, int spacedim>
1508 std::vector<IndexSet>
1510 const DoFHandler<dim, spacedim> &dof_handler)
1511 {
1512 // If the Triangulation is distributed, the only thing we can usefully
1513 // ask is for its locally owned subdomain
1514 Assert((dynamic_cast<
1516 &dof_handler.get_triangulation()) == nullptr),
1517 ExcMessage(
1518 "For parallel::distributed::Triangulation objects and "
1519 "associated DoF handler objects, asking for any information "
1520 "related to a subdomain other than the locally owned one does "
1521 "not make sense."));
1522
1523 // Collect all the locally owned DoFs
1524 // Note: Even though the distribution of DoFs by the
1525 // locally_owned_dofs_per_subdomain function is pseudo-random, we will
1526 // collect all the DoFs on the subdomain and its layer cell. Therefore, the
1527 // random nature of this function does not play a role in the extraction of
1528 // the locally relevant DoFs
1529 std::vector<IndexSet> dof_set =
1531 const ::types::subdomain_id n_subdomains = dof_set.size();
1532
1533 // Add the DoFs on the adjacent (equivalent ghost) cells to the IndexSet,
1534 // cache them in a set. Need to check each DoF manually because we can't
1535 // be sure that the DoF range of locally_owned_dofs is really contiguous.
1537 subdomain_id < n_subdomains;
1538 ++subdomain_id)
1539 {
1540 // Extract the layer of cells around this subdomain
1541 std::function<bool(
1544 const std::vector<
1546 active_halo_layer =
1547 GridTools::compute_active_cell_halo_layer(dof_handler, predicate);
1548
1549 // Extract DoFs associated with halo layer
1550 std::vector<types::global_dof_index> local_dof_indices;
1551 std::set<types::global_dof_index> subdomain_halo_global_dof_indices;
1552 for (typename std::vector<
1554 const_iterator it_cell = active_halo_layer.begin();
1555 it_cell != active_halo_layer.end();
1556 ++it_cell)
1557 {
1559 &cell = *it_cell;
1560 Assert(
1561 cell->subdomain_id() != subdomain_id,
1562 ExcMessage(
1563 "The subdomain ID of the halo cell should not match that of the vector entry."));
1564
1565 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1566 cell->get_dof_indices(local_dof_indices);
1567
1568 for (const types::global_dof_index local_dof_index :
1569 local_dof_indices)
1570 subdomain_halo_global_dof_indices.insert(local_dof_index);
1571 }
1572
1574 subdomain_halo_global_dof_indices.begin(),
1575 subdomain_halo_global_dof_indices.end());
1576
1577 dof_set[subdomain_id].compress();
1578 }
1579
1580 return dof_set;
1581 }
1582
1583 template <int dim, int spacedim>
1584 void
1586 const DoFHandler<dim, spacedim> & dof_handler,
1587 std::vector<types::subdomain_id> &subdomain_association)
1588 {
1589 // if the Triangulation is distributed, the only thing we can usefully
1590 // ask is for its locally owned subdomain
1591 Assert((dynamic_cast<
1593 &dof_handler.get_triangulation()) == nullptr),
1594 ExcMessage(
1595 "For parallel::distributed::Triangulation objects and "
1596 "associated DoF handler objects, asking for any subdomain other "
1597 "than the locally owned one does not make sense."));
1598
1599 Assert(subdomain_association.size() == dof_handler.n_dofs(),
1600 ExcDimensionMismatch(subdomain_association.size(),
1601 dof_handler.n_dofs()));
1602
1603 // catch an error that happened in some versions of the shared tria
1604 // distribute_dofs() function where we were trying to call this
1605 // function at a point in time when not all internal DoFHandler
1606 // structures were quite set up yet.
1607 Assert(dof_handler.n_dofs() > 0, ExcInternalError());
1608
1609 // In case this function is executed with parallel::shared::Triangulation
1610 // with possibly artificial cells, we need to take "true" subdomain IDs
1611 // (i.e. without artificial cells). Otherwise we are good to use
1612 // subdomain_id as stored in cell->subdomain_id().
1613 std::vector<types::subdomain_id> cell_owners(
1614 dof_handler.get_triangulation().n_active_cells());
1617 &dof_handler.get_triangulation())))
1618 {
1619 cell_owners = tr->get_true_subdomain_ids_of_cells();
1620 Assert(tr->get_true_subdomain_ids_of_cells().size() ==
1621 tr->n_active_cells(),
1623 }
1624 else
1625 {
1627 dof_handler.begin_active();
1628 cell != dof_handler.end();
1629 cell++)
1630 if (cell->is_locally_owned())
1631 cell_owners[cell->active_cell_index()] = cell->subdomain_id();
1632 }
1633
1634 // preset all values by an invalid value
1635 std::fill_n(subdomain_association.begin(),
1636 dof_handler.n_dofs(),
1638
1639 std::vector<types::global_dof_index> local_dof_indices;
1640 local_dof_indices.reserve(
1641 dof_handler.get_fe_collection().max_dofs_per_cell());
1642
1643 // loop over all cells and record which subdomain a DoF belongs to.
1644 // give to the smaller subdomain_id in case it is on an interface
1646 cell = dof_handler.begin_active(),
1647 endc = dof_handler.end();
1648 for (; cell != endc; ++cell)
1649 {
1651 cell_owners[cell->active_cell_index()];
1652 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1653 local_dof_indices.resize(dofs_per_cell);
1654 cell->get_dof_indices(local_dof_indices);
1655
1656 // set subdomain ids. if dofs already have their values set then
1657 // they must be on partition interfaces. in that case assign them
1658 // to either the previous association or the current processor
1659 // with the smaller subdomain id.
1660 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1661 if (subdomain_association[local_dof_indices[i]] ==
1663 subdomain_association[local_dof_indices[i]] = subdomain_id;
1664 else if (subdomain_association[local_dof_indices[i]] > subdomain_id)
1665 {
1666 subdomain_association[local_dof_indices[i]] = subdomain_id;
1667 }
1668 }
1669
1670 Assert(std::find(subdomain_association.begin(),
1671 subdomain_association.end(),
1673 subdomain_association.end(),
1675 }
1676
1677
1678
1679 template <int dim, int spacedim>
1680 unsigned int
1682 const DoFHandler<dim, spacedim> &dof_handler,
1683 const types::subdomain_id subdomain)
1684 {
1685 std::vector<types::subdomain_id> subdomain_association(
1686 dof_handler.n_dofs());
1687 get_subdomain_association(dof_handler, subdomain_association);
1688
1689 return std::count(subdomain_association.begin(),
1690 subdomain_association.end(),
1691 subdomain);
1692 }
1693
1694
1695
1696 template <int dim, int spacedim>
1697 IndexSet
1699 const DoFHandler<dim, spacedim> &dof_handler,
1700 const types::subdomain_id subdomain)
1701 {
1702 // If we have a distributed::Triangulation only allow locally_owned
1703 // subdomain.
1706 (subdomain ==
1708 ExcMessage(
1709 "For parallel::distributed::Triangulation objects and "
1710 "associated DoF handler objects, asking for any subdomain other "
1711 "than the locally owned one does not make sense."));
1712
1713 IndexSet index_set(dof_handler.n_dofs());
1714
1715 std::vector<types::global_dof_index> local_dof_indices;
1716 local_dof_indices.reserve(
1717 dof_handler.get_fe_collection().max_dofs_per_cell());
1718
1719 // first generate an unsorted list of all indices which we fill from
1720 // the back. could also insert them directly into the IndexSet, but
1721 // that inserts indices in the middle, which is an O(n^2) algorithm and
1722 // hence too expensive. Could also use std::set, but that is in general
1723 // more expensive than a vector
1724 std::vector<types::global_dof_index> subdomain_indices;
1725
1727 cell = dof_handler.begin_active(),
1728 endc = dof_handler.end();
1729 for (; cell != endc; ++cell)
1730 if ((cell->is_artificial() == false) &&
1731 (cell->subdomain_id() == subdomain))
1732 {
1733 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1734 local_dof_indices.resize(dofs_per_cell);
1735 cell->get_dof_indices(local_dof_indices);
1736 subdomain_indices.insert(subdomain_indices.end(),
1737 local_dof_indices.begin(),
1738 local_dof_indices.end());
1739 }
1740 // sort indices and remove duplicates
1741 std::sort(subdomain_indices.begin(), subdomain_indices.end());
1742 subdomain_indices.erase(std::unique(subdomain_indices.begin(),
1743 subdomain_indices.end()),
1744 subdomain_indices.end());
1745
1746 // insert into IndexSet
1748 index_set.compress();
1749
1750 return index_set;
1751 }
1752
1753
1754
1755 template <int dim, int spacedim>
1756 void
1758 const DoFHandler<dim, spacedim> &dof_handler,
1759 const types::subdomain_id subdomain,
1760 std::vector<unsigned int> & n_dofs_on_subdomain)
1761 {
1762 Assert(n_dofs_on_subdomain.size() == dof_handler.get_fe(0).n_components(),
1763 ExcDimensionMismatch(n_dofs_on_subdomain.size(),
1764 dof_handler.get_fe(0).n_components()));
1765 std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
1766
1767 // Make sure there are at least some cells with this subdomain id
1768 Assert(std::any_of(
1769 dof_handler.begin_active(),
1771 dof_handler.end()},
1772 [subdomain](
1773 const typename DoFHandler<dim, spacedim>::cell_accessor &cell) {
1774 return cell.subdomain_id() == subdomain;
1775 }),
1776 ExcMessage("There are no cells for the given subdomain!"));
1777
1778 std::vector<types::subdomain_id> subdomain_association(
1779 dof_handler.n_dofs());
1780 get_subdomain_association(dof_handler, subdomain_association);
1781
1782 std::vector<unsigned char> component_association(dof_handler.n_dofs());
1784 std::vector<bool>(),
1785 component_association);
1786
1787 for (unsigned int c = 0; c < dof_handler.get_fe(0).n_components(); ++c)
1788 {
1789 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
1790 if ((subdomain_association[i] == subdomain) &&
1791 (component_association[i] == static_cast<unsigned char>(c)))
1792 ++n_dofs_on_subdomain[c];
1793 }
1794 }
1795
1796
1797
1798 namespace internal
1799 {
1800 // TODO: why is this function so complicated? It would be nice to have
1801 // comments that explain why we can't just loop over all components and
1802 // count the entries in dofs_by_component that have this component's
1803 // index
1804 template <int dim, int spacedim>
1805 void
1807 const std::vector<unsigned char> & dofs_by_component,
1808 const std::vector<unsigned int> & target_component,
1809 const bool only_once,
1810 std::vector<types::global_dof_index> &dofs_per_component,
1811 unsigned int & component)
1812 {
1813 for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
1814 {
1816 // Dimension of base element
1817 unsigned int d = base.n_components();
1818
1819 for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
1820 {
1821 if (base.n_base_elements() > 1)
1822 resolve_components(base,
1823 dofs_by_component,
1824 target_component,
1825 only_once,
1826 dofs_per_component,
1827 component);
1828 else
1829 {
1830 for (unsigned int dd = 0; dd < d; ++dd, ++component)
1831 dofs_per_component[target_component[component]] +=
1832 std::count(dofs_by_component.begin(),
1833 dofs_by_component.end(),
1834 component);
1835
1836 // if we have non-primitive FEs and want all components
1837 // to show the number of dofs, need to copy the result to
1838 // those components
1839 if (!base.is_primitive() && !only_once)
1840 for (unsigned int dd = 1; dd < d; ++dd)
1841 dofs_per_component[target_component[component - d + dd]] =
1842 dofs_per_component[target_component[component - d]];
1843 }
1844 }
1845 }
1846 }
1847
1848
1849 template <int dim, int spacedim>
1850 void
1852 const std::vector<unsigned char> & dofs_by_component,
1853 const std::vector<unsigned int> & target_component,
1854 const bool only_once,
1855 std::vector<types::global_dof_index> &dofs_per_component,
1856 unsigned int & component)
1857 {
1858 // assert that all elements in the collection have the same structure
1859 // (base elements and multiplicity, components per base element) and
1860 // then simply call the function above
1861 for (unsigned int fe = 1; fe < fe_collection.size(); ++fe)
1862 {
1863 Assert(fe_collection[fe].n_components() ==
1864 fe_collection[0].n_components(),
1866 Assert(fe_collection[fe].n_base_elements() ==
1867 fe_collection[0].n_base_elements(),
1869 for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b)
1870 {
1871 Assert(fe_collection[fe].base_element(b).n_components() ==
1872 fe_collection[0].base_element(b).n_components(),
1874 Assert(fe_collection[fe].base_element(b).n_base_elements() ==
1875 fe_collection[0].base_element(b).n_base_elements(),
1877 }
1878 }
1879
1880 resolve_components(fe_collection[0],
1881 dofs_by_component,
1882 target_component,
1883 only_once,
1884 dofs_per_component,
1885 component);
1886 }
1887 } // namespace internal
1888
1889
1890
1891 namespace internal
1892 {
1893 namespace
1894 {
1898 template <int dim, int spacedim>
1899 bool
1900 all_elements_are_primitive(const FiniteElement<dim, spacedim> &fe)
1901 {
1902 return fe.is_primitive();
1903 }
1904
1905
1910 template <int dim, int spacedim>
1911 bool
1912 all_elements_are_primitive(
1913 const ::hp::FECollection<dim, spacedim> &fe_collection)
1914 {
1915 for (unsigned int i = 0; i < fe_collection.size(); ++i)
1916 if (fe_collection[i].is_primitive() == false)
1917 return false;
1918
1919 return true;
1920 }
1921 } // namespace
1922 } // namespace internal
1923
1924
1925
1926 // deprecated function
1927 template <int dim, int spacedim>
1928 void
1930 const DoFHandler<dim, spacedim> & dof_handler,
1931 std::vector<types::global_dof_index> &dofs_per_component,
1932 const bool only_once,
1933 const std::vector<unsigned int> & target_component)
1934 {
1935 dofs_per_component =
1936 count_dofs_per_fe_component(dof_handler, only_once, target_component);
1937 }
1938
1939
1940
1941 template <int dim, int spacedim>
1942 std::vector<types::global_dof_index>
1944 const DoFHandler<dim, spacedim> &dof_handler,
1945 const bool only_once,
1946 const std::vector<unsigned int> &target_component_)
1947 {
1948 const unsigned int n_components = dof_handler.get_fe(0).n_components();
1949
1950 // If the empty vector was given as default argument, set up this
1951 // vector as identity.
1952 std::vector<unsigned int> target_component = target_component_;
1953 if (target_component.size() == 0)
1954 {
1955 target_component.resize(n_components);
1956 for (unsigned int i = 0; i < n_components; ++i)
1957 target_component[i] = i;
1958 }
1959 else
1960 Assert(target_component.size() == n_components,
1961 ExcDimensionMismatch(target_component.size(), n_components));
1962
1963
1964 const unsigned int max_component =
1965 *std::max_element(target_component.begin(), target_component.end());
1966 const unsigned int n_target_components = max_component + 1;
1967
1968 std::vector<types::global_dof_index> dofs_per_component(
1969 n_target_components, types::global_dof_index(0));
1970
1971 // special case for only one component. treat this first since it does
1972 // not require any computations
1973 if (n_components == 1)
1974 {
1975 dofs_per_component[0] = dof_handler.n_locally_owned_dofs();
1976 return dofs_per_component;
1977 }
1978
1979
1980 // otherwise determine the number of dofs in each component separately.
1981 // do so in parallel
1982 std::vector<unsigned char> dofs_by_component(
1983 dof_handler.n_locally_owned_dofs());
1986 dofs_by_component);
1987
1988 // next count what we got
1989 unsigned int component = 0;
1991 dofs_by_component,
1992 target_component,
1993 only_once,
1994 dofs_per_component,
1995 component);
1996 Assert(n_components == component, ExcInternalError());
1997
1998 // finally sanity check. this is only valid if the finite element is
1999 // actually primitive, so exclude other elements from this
2000 Assert((internal::all_elements_are_primitive(
2001 dof_handler.get_fe_collection()) == false) ||
2002 (std::accumulate(dofs_per_component.begin(),
2003 dofs_per_component.end(),
2005 dof_handler.n_locally_owned_dofs()),
2007
2008 // reduce information from all CPUs
2009#ifdef DEAL_II_WITH_MPI
2010
2012 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2013 &dof_handler.get_triangulation())))
2014 {
2015 std::vector<types::global_dof_index> local_dof_count =
2016 dofs_per_component;
2017
2018 const int ierr = MPI_Allreduce(local_dof_count.data(),
2019 dofs_per_component.data(),
2020 n_target_components,
2022 MPI_SUM,
2023 tria->get_communicator());
2024 AssertThrowMPI(ierr);
2025 }
2026#endif
2027
2028 return dofs_per_component;
2029 }
2030
2031
2032
2033 // deprecated function
2034 template <int dim, int spacedim>
2035 void
2037 std::vector<types::global_dof_index> &dofs_per_block,
2038 const std::vector<unsigned int> & target_block)
2039 {
2040 dofs_per_block = count_dofs_per_fe_block(dof_handler, target_block);
2041 }
2042
2043
2044
2045 template <int dim, int spacedim>
2046 std::vector<types::global_dof_index>
2048 const std::vector<unsigned int> &target_block_)
2049 {
2050 const ::hp::FECollection<dim, spacedim> &fe_collection =
2051 dof_handler.get_fe_collection();
2052 Assert(fe_collection.size() < 256, ExcNotImplemented());
2053
2054 // If the empty vector for target_block(e.g., as default argument), then
2055 // set up this vector as identity. We do this set up with the first
2056 // element of the collection, but the whole thing can only work if
2057 // all elements have the same number of blocks anyway -- so check
2058 // that right after
2059 const unsigned int n_blocks = fe_collection[0].n_blocks();
2060
2061 std::vector<unsigned int> target_block = target_block_;
2062 if (target_block.size() == 0)
2063 {
2064 target_block.resize(fe_collection[0].n_blocks());
2065 for (unsigned int i = 0; i < n_blocks; ++i)
2066 target_block[i] = i;
2067 }
2068 else
2069 Assert(target_block.size() == n_blocks,
2070 ExcDimensionMismatch(target_block.size(), n_blocks));
2071 for (unsigned int f = 1; f < fe_collection.size(); ++f)
2072 Assert(fe_collection[0].n_blocks() == fe_collection[f].n_blocks(),
2073 ExcMessage("This function can only work if all elements in a "
2074 "collection have the same number of blocks."));
2075
2076 // special case for only one block. treat this first since it does
2077 // not require any computations
2078 if (n_blocks == 1)
2079 {
2080 std::vector<types::global_dof_index> dofs_per_block(1);
2081 dofs_per_block[0] = dof_handler.n_dofs();
2082 return dofs_per_block;
2083 }
2084
2085 // Otherwise set up the right-sized object and start working
2086 const unsigned int max_block =
2087 *std::max_element(target_block.begin(), target_block.end());
2088 const unsigned int n_target_blocks = max_block + 1;
2089
2090 std::vector<types::global_dof_index> dofs_per_block(n_target_blocks);
2091
2092 // Loop over the elements of the collection, but really only consider
2093 // the last element (see #9271)
2094 for (unsigned int this_fe = fe_collection.size() - 1;
2095 this_fe < fe_collection.size();
2096 ++this_fe)
2097 {
2098 const FiniteElement<dim, spacedim> &fe = fe_collection[this_fe];
2099
2100 std::vector<unsigned char> dofs_by_block(
2101 dof_handler.n_locally_owned_dofs());
2102 internal::get_block_association(dof_handler, dofs_by_block);
2103
2104 // next count what we got
2105 for (unsigned int block = 0; block < fe.n_blocks(); ++block)
2106 dofs_per_block[target_block[block]] +=
2107 std::count(dofs_by_block.begin(), dofs_by_block.end(), block);
2108
2109#ifdef DEAL_II_WITH_MPI
2110 // if we are working on a parallel mesh, we now need to collect
2111 // this information from all processors
2113 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2114 &dof_handler.get_triangulation())))
2115 {
2116 std::vector<types::global_dof_index> local_dof_count =
2117 dofs_per_block;
2118 const int ierr = MPI_Allreduce(local_dof_count.data(),
2119 dofs_per_block.data(),
2120 n_target_blocks,
2122 MPI_SUM,
2123 tria->get_communicator());
2124 AssertThrowMPI(ierr);
2125 }
2126#endif
2127 }
2128
2129 return dofs_per_block;
2130 }
2131
2132
2133
2134 template <int dim, int spacedim>
2135 void
2137 std::vector<types::global_dof_index> &mapping)
2138 {
2139 mapping.clear();
2140 mapping.insert(mapping.end(),
2141 dof_handler.n_dofs(),
2143
2144 std::vector<types::global_dof_index> dofs_on_face;
2145 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2146 types::global_dof_index next_boundary_index = 0;
2147
2148 // now loop over all cells and check whether their faces are at the
2149 // boundary. note that we need not take special care of single lines
2150 // being at the boundary (using @p{cell->has_boundary_lines}), since we
2151 // do not support boundaries of dimension dim-2, and so every isolated
2152 // boundary line is also part of a boundary face which we will be
2153 // visiting sooner or later
2155 cell = dof_handler.begin_active(),
2156 endc = dof_handler.end();
2157 for (; cell != endc; ++cell)
2158 for (const unsigned int f : cell->face_indices())
2159 if (cell->at_boundary(f))
2160 {
2161 const unsigned int dofs_per_face =
2162 cell->get_fe().n_dofs_per_face(f);
2163 dofs_on_face.resize(dofs_per_face);
2164 cell->face(f)->get_dof_indices(dofs_on_face,
2165 cell->active_fe_index());
2166 for (unsigned int i = 0; i < dofs_per_face; ++i)
2167 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2168 mapping[dofs_on_face[i]] = next_boundary_index++;
2169 }
2170
2171 AssertDimension(next_boundary_index, dof_handler.n_boundary_dofs());
2172 }
2173
2174
2175
2176 template <int dim, int spacedim>
2177 void
2179 const std::set<types::boundary_id> &boundary_ids,
2180 std::vector<types::global_dof_index> &mapping)
2181 {
2182 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2183 boundary_ids.end(),
2185
2186 mapping.clear();
2187 mapping.insert(mapping.end(),
2188 dof_handler.n_dofs(),
2190
2191 // return if there is nothing to do
2192 if (boundary_ids.size() == 0)
2193 return;
2194
2195 std::vector<types::global_dof_index> dofs_on_face;
2196 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2197 types::global_dof_index next_boundary_index = 0;
2198
2200 cell = dof_handler.begin_active(),
2201 endc = dof_handler.end();
2202 for (; cell != endc; ++cell)
2203 for (const unsigned int f : cell->face_indices())
2204 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
2205 boundary_ids.end())
2206 {
2207 const unsigned int dofs_per_face =
2208 cell->get_fe().n_dofs_per_face(f);
2209 dofs_on_face.resize(dofs_per_face);
2210 cell->face(f)->get_dof_indices(dofs_on_face,
2211 cell->active_fe_index());
2212 for (unsigned int i = 0; i < dofs_per_face; ++i)
2213 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2214 mapping[dofs_on_face[i]] = next_boundary_index++;
2215 }
2216
2217 AssertDimension(next_boundary_index,
2218 dof_handler.n_boundary_dofs(boundary_ids));
2219 }
2220
2221 namespace internal
2222 {
2223 namespace
2224 {
2225 template <int dim, int spacedim>
2226 void
2229 const DoFHandler<dim, spacedim> & dof_handler,
2230 std::map<types::global_dof_index, Point<spacedim>> &support_points,
2232 {
2233 const hp::FECollection<dim, spacedim> &fe_collection =
2234 dof_handler.get_fe_collection();
2235 hp::QCollection<dim> q_coll_dummy;
2236
2237 for (unsigned int fe_index = 0; fe_index < fe_collection.size();
2238 ++fe_index)
2239 {
2240 // check whether every FE in the collection has support points
2241 Assert(fe_collection[fe_index].has_support_points(),
2244 fe_collection[fe_index].get_unit_support_points()));
2245 }
2246
2247 // Take care of components
2252
2253 // Now loop over all cells and enquire the support points on each
2254 // of these. we use dummy quadrature formulas where the quadrature
2255 // points are located at the unit support points to enquire the
2256 // location of the support points in real space.
2257 //
2258 // The weights of the quadrature rule have been set to invalid
2259 // values by the used constructor.
2260 hp::FEValues<dim, spacedim> hp_fe_values(mapping,
2261 fe_collection,
2262 q_coll_dummy,
2265 cell = dof_handler.begin_active(),
2266 endc = dof_handler.end();
2267
2268 std::vector<types::global_dof_index> local_dof_indices;
2269 for (; cell != endc; ++cell)
2270 // only work on locally relevant cells
2271 if (cell->is_artificial() == false)
2272 {
2273 hp_fe_values.reinit(cell);
2274 const FEValues<dim, spacedim> &fe_values =
2275 hp_fe_values.get_present_fe_values();
2276
2277 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2278 cell->get_dof_indices(local_dof_indices);
2279
2280 const std::vector<Point<spacedim>> &points =
2282 for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell();
2283 ++i)
2284 {
2285 const unsigned int dof_comp =
2286 cell->get_fe().system_to_component_index(i).first;
2287
2288 // insert the values into the map if it is a valid component
2290 support_points[local_dof_indices[i]] = points[i];
2291 }
2292 }
2293 }
2294
2295
2296 template <int dim, int spacedim>
2297 void
2300 const DoFHandler<dim, spacedim> & dof_handler,
2301 std::vector<Point<spacedim>> & support_points,
2303 {
2304 // get the data in the form of the map as above
2305 std::map<types::global_dof_index, Point<spacedim>> x_support_points;
2307 dof_handler,
2308 x_support_points,
2310
2311 // now convert from the map to the linear vector. make sure every
2312 // entry really appeared in the map
2313 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
2314 {
2315 Assert(x_support_points.find(i) != x_support_points.end(),
2317
2318 support_points[i] = x_support_points[i];
2319 }
2320 }
2321 } // namespace
2322 } // namespace internal
2323
2324 template <int dim, int spacedim>
2325 void
2327 const DoFHandler<dim, spacedim> &dof_handler,
2328 std::vector<Point<spacedim>> & support_points,
2330 {
2331 AssertDimension(support_points.size(), dof_handler.n_dofs());
2332 Assert((dynamic_cast<
2334 &dof_handler.get_triangulation()) == nullptr),
2335 ExcMessage(
2336 "This function can not be used with distributed triangulations. "
2338
2339 // Let the internal function do all the work, just make sure that it
2340 // gets a MappingCollection
2341 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2342
2343 internal::map_dofs_to_support_points(mapping_collection,
2344 dof_handler,
2345 support_points,
2347 }
2348
2349
2350 template <int dim, int spacedim>
2351 void
2354 const DoFHandler<dim, spacedim> & dof_handler,
2355 std::vector<Point<spacedim>> & support_points,
2357 {
2358 AssertDimension(support_points.size(), dof_handler.n_dofs());
2359 Assert((dynamic_cast<
2361 &dof_handler.get_triangulation()) == nullptr),
2362 ExcMessage(
2363 "This function can not be used with distributed triangulations. "
2365
2366 // Let the internal function do all the work, just make sure that it
2367 // gets a MappingCollection
2369 dof_handler,
2370 support_points,
2372 }
2373
2374
2375 template <int dim, int spacedim>
2376 void
2378 const Mapping<dim, spacedim> & mapping,
2379 const DoFHandler<dim, spacedim> & dof_handler,
2380 std::map<types::global_dof_index, Point<spacedim>> &support_points,
2382 {
2383 support_points.clear();
2384
2385 // Let the internal function do all the work, just make sure that it
2386 // gets a MappingCollection
2387 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2388
2389 internal::map_dofs_to_support_points(mapping_collection,
2390 dof_handler,
2391 support_points,
2393 }
2394
2395
2396 template <int dim, int spacedim>
2397 void
2400 const DoFHandler<dim, spacedim> & dof_handler,
2401 std::map<types::global_dof_index, Point<spacedim>> &support_points,
2403 {
2404 support_points.clear();
2405
2406 // Let the internal function do all the work, just make sure that it
2407 // gets a MappingCollection
2409 dof_handler,
2410 support_points,
2412 }
2413
2414 template <int spacedim>
2415 void
2417 std::ostream & out,
2418 const std::map<types::global_dof_index, Point<spacedim>> &support_points)
2419 {
2420 AssertThrow(out, ExcIO());
2421
2422 using dof_map_t = std::map<types::global_dof_index, Point<spacedim>>;
2423
2424 using point_map_t = std::map<Point<spacedim>,
2425 std::vector<types::global_dof_index>,
2427
2428 point_map_t point_map;
2429
2430 // convert to map point -> list of DoFs
2431 for (typename dof_map_t::const_iterator it = support_points.begin();
2432 it != support_points.end();
2433 ++it)
2434 {
2435 std::vector<types::global_dof_index> &v = point_map[it->second];
2436 v.push_back(it->first);
2437 }
2438
2439 // print the newly created map:
2440 for (typename point_map_t::iterator it = point_map.begin();
2441 it != point_map.end();
2442 ++it)
2443 {
2444 out << it->first << " \"";
2445 const std::vector<types::global_dof_index> &v = it->second;
2446 for (unsigned int i = 0; i < v.size(); ++i)
2447 {
2448 if (i > 0)
2449 out << ", ";
2450 out << v[i];
2451 }
2452
2453 out << "\"\n";
2454 }
2455
2456 out << std::flush;
2457 }
2458
2459
2460 template <int dim, int spacedim>
2461 void
2463 const Table<2, Coupling> & table,
2464 std::vector<Table<2, Coupling>> &tables_by_block)
2465 {
2466 if (dof_handler.has_hp_capabilities() == false)
2467 {
2468 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
2469 const unsigned int nb = fe.n_blocks();
2470
2471 tables_by_block.resize(1);
2472 tables_by_block[0].reinit(nb, nb);
2473 tables_by_block[0].fill(none);
2474
2475 for (unsigned int i = 0; i < fe.n_components(); ++i)
2476 {
2477 const unsigned int ib = fe.component_to_block_index(i);
2478 for (unsigned int j = 0; j < fe.n_components(); ++j)
2479 {
2480 const unsigned int jb = fe.component_to_block_index(j);
2481 tables_by_block[0](ib, jb) |= table(i, j);
2482 }
2483 }
2484 }
2485 else
2486 {
2487 const hp::FECollection<dim> &fe_collection =
2488 dof_handler.get_fe_collection();
2489 tables_by_block.resize(fe_collection.size());
2490
2491 for (unsigned int f = 0; f < fe_collection.size(); ++f)
2492 {
2493 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
2494
2495 const unsigned int nb = fe.n_blocks();
2496 tables_by_block[f].reinit(nb, nb);
2497 tables_by_block[f].fill(none);
2498 for (unsigned int i = 0; i < fe.n_components(); ++i)
2499 {
2500 const unsigned int ib = fe.component_to_block_index(i);
2501 for (unsigned int j = 0; j < fe.n_components(); ++j)
2502 {
2503 const unsigned int jb = fe.component_to_block_index(j);
2504 tables_by_block[f](ib, jb) |= table(i, j);
2505 }
2506 }
2507 }
2508 }
2509 }
2510
2511
2512
2513 template <int dim, int spacedim>
2514 void
2516 const DoFHandler<dim, spacedim> &dof_handler,
2517 const unsigned int level,
2518 const std::vector<bool> & selected_dofs,
2519 const types::global_dof_index offset)
2520 {
2523 dof_handler.end(level);
2524 std::vector<types::global_dof_index> indices;
2525
2526 unsigned int i = 0;
2527
2528 for (cell = dof_handler.begin(level); cell != endc; ++cell)
2529 if (cell->is_locally_owned_on_level())
2530 ++i;
2531 block_list.reinit(i,
2532 dof_handler.n_dofs(),
2533 dof_handler.get_fe().n_dofs_per_cell());
2534 i = 0;
2535 for (cell = dof_handler.begin(level); cell != endc; ++cell)
2536 if (cell->is_locally_owned_on_level())
2537 {
2538 indices.resize(cell->get_fe().n_dofs_per_cell());
2539 cell->get_mg_dof_indices(indices);
2540
2541 if (selected_dofs.size() != 0)
2542 AssertDimension(indices.size(), selected_dofs.size());
2543
2544 for (types::global_dof_index j = 0; j < indices.size(); ++j)
2545 {
2546 if (selected_dofs.size() == 0)
2548 else
2549 {
2550 if (selected_dofs[j])
2552 }
2553 }
2554 ++i;
2555 }
2556 }
2557
2558
2559 template <int dim, int spacedim>
2560 void
2562 const DoFHandler<dim, spacedim> &dof_handler,
2563 const unsigned int level,
2564 const bool interior_only)
2565 {
2566 const FiniteElement<dim> &fe = dof_handler.get_fe();
2567 block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level));
2570 dof_handler.end(level);
2571
2572 std::vector<types::global_dof_index> indices;
2573 std::vector<bool> exclude;
2574
2575 for (cell = dof_handler.begin(level); cell != endc; ++cell)
2576 {
2577 indices.resize(cell->get_fe().n_dofs_per_cell());
2578 cell->get_mg_dof_indices(indices);
2579
2580 if (interior_only)
2581 {
2582 // Exclude degrees of freedom on faces opposite to the vertex
2583 exclude.resize(fe.n_dofs_per_cell());
2584 std::fill(exclude.begin(), exclude.end(), false);
2585
2586 for (const unsigned int face : cell->face_indices())
2587 if (cell->at_boundary(face) ||
2588 cell->neighbor(face)->level() != cell->level())
2589 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2590 exclude[fe.face_to_cell_index(i, face)] = true;
2591 for (types::global_dof_index j = 0; j < indices.size(); ++j)
2592 if (!exclude[j])
2594 }
2595 else
2596 {
2597 for (const auto index : indices)
2599 }
2600 }
2601 }
2602
2603
2604 template <int dim, int spacedim>
2605 void
2607 const DoFHandler<dim, spacedim> &dof_handler,
2608 const unsigned int level,
2609 const bool interior_dofs_only,
2610 const bool boundary_dofs)
2611 {
2612 Assert(level > 0 && level < dof_handler.get_triangulation().n_levels(),
2613 ExcIndexRange(level, 1, dof_handler.get_triangulation().n_levels()));
2614
2616 dof_handler.begin(level - 1);
2618 dof_handler.end(level - 1);
2619
2620 std::vector<types::global_dof_index> indices;
2621 std::vector<bool> exclude;
2622
2623 for (unsigned int block = 0; pcell != endc; ++pcell)
2624 {
2625 if (pcell->is_active())
2626 continue;
2627
2628 for (unsigned int child = 0; child < pcell->n_children(); ++child)
2629 {
2631 pcell->child(child);
2632
2633 // For hp, only this line here would have to be replaced.
2634 const FiniteElement<dim> &fe = dof_handler.get_fe();
2635 const unsigned int n_dofs = fe.n_dofs_per_cell();
2636 indices.resize(n_dofs);
2637 exclude.resize(n_dofs);
2638 std::fill(exclude.begin(), exclude.end(), false);
2639 cell->get_mg_dof_indices(indices);
2640
2641 if (interior_dofs_only)
2642 {
2643 // Eliminate dofs on faces of the child which are on faces
2644 // of the parent
2645 for (unsigned int d = 0; d < dim; ++d)
2646 {
2647 const unsigned int face =
2649 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2650 exclude[fe.face_to_cell_index(i, face)] = true;
2651 }
2652
2653 // Now remove all degrees of freedom on the domain boundary
2654 // from the exclusion list
2655 if (boundary_dofs)
2656 for (const unsigned int face :
2658 if (cell->at_boundary(face))
2659 for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
2660 ++i)
2661 exclude[fe.face_to_cell_index(i, face)] = false;
2662 }
2663
2664 for (unsigned int i = 0; i < n_dofs; ++i)
2665 if (!exclude[i])
2667 }
2668 ++block;
2669 }
2670 }
2671
2672 template <int dim, int spacedim>
2673 std::vector<unsigned int>
2675 const DoFHandler<dim, spacedim> &dof_handler,
2676 const unsigned int level,
2677 const bool interior_only,
2678 const bool boundary_patches,
2679 const bool level_boundary_patches,
2680 const bool single_cell_patches,
2681 const bool invert_vertex_mapping)
2682 {
2683 const unsigned int n_blocks = dof_handler.get_fe().n_blocks();
2685 return make_vertex_patches(block_list,
2686 dof_handler,
2687 level,
2688 exclude_boundary_dofs,
2689 boundary_patches,
2690 level_boundary_patches,
2691 single_cell_patches,
2692 invert_vertex_mapping);
2693 }
2694
2695 template <int dim, int spacedim>
2696 std::vector<unsigned int>
2698 const DoFHandler<dim, spacedim> &dof_handler,
2699 const unsigned int level,
2701 const bool boundary_patches,
2702 const bool level_boundary_patches,
2703 const bool single_cell_patches,
2704 const bool invert_vertex_mapping)
2705 {
2708 dof_handler.end(level);
2709
2710 // Vector mapping from vertex index in the triangulation to consecutive
2711 // block indices on this level The number of cells at a vertex
2712 std::vector<unsigned int> vertex_cell_count(
2713 dof_handler.get_triangulation().n_vertices(), 0);
2714
2715 // Is a vertex at the boundary?
2716 std::vector<bool> vertex_boundary(
2717 dof_handler.get_triangulation().n_vertices(), false);
2718
2719 std::vector<unsigned int> vertex_mapping(
2720 dof_handler.get_triangulation().n_vertices(),
2722
2723 // Estimate for the number of dofs at this point
2724 std::vector<unsigned int> vertex_dof_count(
2725 dof_handler.get_triangulation().n_vertices(), 0);
2726
2727 // Identify all vertices active on this level and remember some data
2729 for (cell = dof_handler.begin(level); cell != endc; ++cell)
2730 for (const unsigned int v : cell->vertex_indices())
2731 {
2732 const unsigned int vg = cell->vertex_index(v);
2733 vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell();
2734 ++vertex_cell_count[vg];
2735 for (unsigned int d = 0; d < dim; ++d)
2736 {
2737 const unsigned int face = GeometryInfo<dim>::vertex_to_face[v][d];
2738 if (cell->at_boundary(face))
2739 vertex_boundary[vg] = true;
2740 else if ((!level_boundary_patches) &&
2741 (cell->neighbor(face)->level() !=
2742 static_cast<int>(level)))
2743 vertex_boundary[vg] = true;
2744 }
2745 }
2746 // From now on, only vertices with positive dof count are "in".
2747
2748 // Remove vertices at boundaries or in corners
2749 for (unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg)
2750 if ((!single_cell_patches && vertex_cell_count[vg] < 2) ||
2751 (!boundary_patches && vertex_boundary[vg]))
2752 vertex_dof_count[vg] = 0;
2753
2754 // Create a mapping from all vertices to the ones used here
2755 unsigned int n_vertex_count = 0;
2756 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2757 if (vertex_dof_count[vg] != 0)
2758 vertex_mapping[vg] = n_vertex_count++;
2759
2760 // Compactify dof count
2761 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2762 if (vertex_dof_count[vg] != 0)
2763 vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
2764
2765 // Now that we have all the data, we reduce it to the part we actually
2766 // want
2767 vertex_dof_count.resize(n_vertex_count);
2768
2769 // At this point, the list of patches is ready. Now we enter the dofs
2770 // into the sparsity pattern.
2771 block_list.reinit(vertex_dof_count.size(),
2772 dof_handler.n_dofs(level),
2773 vertex_dof_count);
2774
2775 std::vector<types::global_dof_index> indices;
2776 std::vector<bool> exclude;
2777
2778 for (cell = dof_handler.begin(level); cell != endc; ++cell)
2779 {
2780 const FiniteElement<dim> &fe = cell->get_fe();
2781 indices.resize(fe.n_dofs_per_cell());
2782 cell->get_mg_dof_indices(indices);
2783
2784 for (const unsigned int v : cell->vertex_indices())
2785 {
2786 const unsigned int vg = cell->vertex_index(v);
2787 const unsigned int block = vertex_mapping[vg];
2788 if (block == numbers::invalid_unsigned_int)
2789 continue;
2790
2791 // Collect excluded dofs for some block(s) if boundary dofs
2792 // for a block are decided to be excluded
2793 if (exclude_boundary_dofs.size() == 0 ||
2794 exclude_boundary_dofs.n_selected_blocks() != 0)
2795 {
2796 // Exclude degrees of freedom on faces opposite to the
2797 // vertex
2798 exclude.resize(fe.n_dofs_per_cell());
2799 std::fill(exclude.begin(), exclude.end(), false);
2800
2801 for (unsigned int d = 0; d < dim; ++d)
2802 {
2803 const unsigned int a_face =
2805 const unsigned int face =
2807 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2808 {
2809 // For each dof, get the block it is in and decide to
2810 // exclude it or not
2811 if (exclude_boundary_dofs[fe.system_to_block_index(
2813 i, face))
2814 .first] == true)
2815 exclude[fe.face_to_cell_index(i, face)] = true;
2816 }
2817 }
2818 for (unsigned int j = 0; j < indices.size(); ++j)
2819 if (!exclude[j])
2821 }
2822 else
2823 {
2824 for (const auto index : indices)
2826 }
2827 }
2828 }
2829
2830 if (invert_vertex_mapping)
2831 {
2832 // Compress vertex mapping
2833 unsigned int n_vertex_count = 0;
2834 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2835 if (vertex_mapping[vg] != numbers::invalid_unsigned_int)
2836 vertex_mapping[n_vertex_count++] = vg;
2837
2838 // Now we reduce it to the part we actually want
2839 vertex_mapping.resize(n_vertex_count);
2840 }
2841
2842 return vertex_mapping;
2843 }
2844
2845
2846 template <int dim, int spacedim>
2847 unsigned int
2849 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2850 &patch)
2851 {
2852 std::set<types::global_dof_index> dofs_on_patch;
2853 std::vector<types::global_dof_index> local_dof_indices;
2854
2855 // loop over the cells in the patch and get the DoFs on each.
2856 // add all of them to a std::set which automatically makes sure
2857 // all duplicates are ignored
2858 for (unsigned int i = 0; i < patch.size(); ++i)
2859 {
2861 patch[i];
2862 Assert(cell->is_artificial() == false,
2863 ExcMessage("This function can not be called with cells that are "
2864 "not either locally owned or ghost cells."));
2865 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2866 cell->get_dof_indices(local_dof_indices);
2867 dofs_on_patch.insert(local_dof_indices.begin(),
2868 local_dof_indices.end());
2869 }
2870
2871 // now return the number of DoFs (duplicates were ignored)
2872 return dofs_on_patch.size();
2873 }
2874
2875
2876
2877 template <int dim, int spacedim>
2878 std::vector<types::global_dof_index>
2880 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2881 &patch)
2882 {
2883 std::set<types::global_dof_index> dofs_on_patch;
2884 std::vector<types::global_dof_index> local_dof_indices;
2885
2886 // loop over the cells in the patch and get the DoFs on each.
2887 // add all of them to a std::set which automatically makes sure
2888 // all duplicates are ignored
2889 for (unsigned int i = 0; i < patch.size(); ++i)
2890 {
2892 patch[i];
2893 Assert(cell->is_artificial() == false,
2894 ExcMessage("This function can not be called with cells that are "
2895 "not either locally owned or ghost cells."));
2896 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2897 cell->get_dof_indices(local_dof_indices);
2898 dofs_on_patch.insert(local_dof_indices.begin(),
2899 local_dof_indices.end());
2900 }
2901
2902 Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)),
2904
2905 // return a vector with the content of the set above. copying
2906 // also ensures that we retain sortedness as promised in the
2907 // documentation and as necessary to retain the block structure
2908 // also on the local system
2909 return std::vector<types::global_dof_index>(dofs_on_patch.begin(),
2910 dofs_on_patch.end());
2911 }
2912
2913
2914} // end of namespace DoFTools
2915
2916
2917
2918// explicit instantiations
2919
2920#include "dof_tools.inst"
2921
2922
2923
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
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
cell_iterator end() 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
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:1921
size_type n_elements() const
Definition: index_set.h:1832
bool is_element(const size_type index) const
Definition: index_set.h:1765
Definition: index_set.h:1653
void fill_binary_vector(VectorType &vector) const
void subtract_set(const IndexSet &other)
Definition: index_set.cc:258
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1880
void compress() const
Definition: index_set.h:1642
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1703
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:110
unsigned int size() const
Definition: collection.h:109
unsigned int max_dofs_per_face() const
unsigned int n_components() const
unsigned int max_dofs_per_cell() const
Definition: q_collection.h:223
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
unsigned int level
Definition: grid_out.cc:4590
unsigned int cell_index
Definition: grid_tools.cc:1092
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()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1746
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths) override
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
void add(const size_type i, const size_type j)
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:1806
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:1585
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:613
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:864
IndexSet dof_indices_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
Definition: dof_tools.cc:1698
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1414
void count_dofs_per_component(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::global_dof_index > &dofs_per_component, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
Definition: dof_tools.cc:1929
void extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1252
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:2326
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1399
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:2515
void extract_locally_active_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1139
std::vector< IndexSet > locally_owned_dofs_per_component(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &components=ComponentMask())
Definition: dof_tools.cc:506
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:2416
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:1105
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:2674
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:2047
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:540
void extract_hanging_node_dofs(const DoFHandler< dim, spacedim > &dof_handler, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:1079
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:782
void map_dofs_to_support_points(const hp::MappingCollection< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, std::map< types::global_dof_index, Point< spacedim > > &support_points, const ComponentMask &mask)
Definition: dof_tools.cc:2398
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:1943
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:2606
void map_dof_to_boundary_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::global_dof_index > &mapping)
Definition: dof_tools.cc:2136
std::vector< IndexSet > locally_relevant_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1509
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:2561
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:2879
void extract_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:393
unsigned int count_dofs_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
Definition: dof_tools.cc:1681
void extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1210
unsigned int count_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
Definition: dof_tools.cc:2848
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:2462
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:1304
void count_dofs_per_block(const DoFHandler< dim, spacedim > &dof, std::vector< types::global_dof_index > &dofs_per_block, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition: dof_tools.cc:2036
void extract_locally_active_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set, const unsigned int level)
Definition: dof_tools.cc:1174
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
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)
std::enable_if< IsBlockVector< VectorType >::value, unsignedint >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
std::string compress(const std::string &input)
Definition: utilities.cc:392
const types::boundary_id internal_face_boundary_id
Definition: types.h:255
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
const types::subdomain_id invalid_subdomain_id
Definition: types.h:276
static const unsigned int invalid_unsigned_int
Definition: types.h:196
const types::global_dof_index invalid_dof_index
Definition: types.h:211
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned int subdomain_id
Definition: types.h:43
bool operator()(const Point< dim, Number > &lhs, const Point< dim, Number > &rhs) const
Definition: dof_tools.cc:85
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:86