Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/table.h>
19 #include <deal.II/base/utilities.h>
20 
23 
26 #include <deal.II/dofs/dof_tools.h>
27 
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_tools.h>
30 #include <deal.II/fe/fe_values.h>
31 
35 #include <deal.II/grid/tria.h>
37 
38 #include <deal.II/hp/dof_handler.h>
40 #include <deal.II/hp/fe_values.h>
43 
49 #include <deal.II/lac/vector.h>
50 
51 #include <algorithm>
52 #include <numeric>
53 
55 
56 
57 
58 namespace DoFTools
59 {
60  namespace internal
61  {
76  template <int dim, typename Number = double>
78  {
84  bool
86  const Point<dim, Number> &rhs) const
87  {
88  double downstream_size = 0;
89  double weight = 1.;
90  for (unsigned int d = 0; d < dim; ++d)
91  {
92  downstream_size += (rhs[d] - lhs[d]) * weight;
93  weight *= 1e-5;
94  }
95  if (downstream_size < 0)
96  return false;
97  else if (downstream_size > 0)
98  return true;
99  else
100  {
101  for (unsigned int d = 0; d < dim; ++d)
102  {
103  if (lhs[d] == rhs[d])
104  continue;
105  return lhs[d] < rhs[d];
106  }
107  return false;
108  }
109  }
110  };
111 
112 
113 
114  // return an array that for each dof on the reference cell
115  // lists the corresponding vector component.
116  //
117  // if an element is non-primitive then we assign to each degree of freedom
118  // the following component:
119  // - if the nonzero components that belong to a shape function are not
120  // selected in the component_mask, then the shape function is assigned
121  // to the first nonzero vector component that corresponds to this
122  // shape function
123  // - otherwise, the shape function is assigned the first component selected
124  // in the component_mask that corresponds to this shape function
125  template <int dim, int spacedim>
126  std::vector<unsigned char>
128  const ComponentMask &component_mask)
129  {
130  std::vector<unsigned char> local_component_association(
131  fe.n_dofs_per_cell(), static_cast<unsigned char>(-1));
132 
133  // compute the component each local dof belongs to.
134  // if the shape function is primitive, then this
135  // is simple and we can just associate it with
136  // what system_to_component_index gives us
137  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
138  if (fe.is_primitive(i))
139  local_component_association[i] =
140  fe.system_to_component_index(i).first;
141  else
142  // if the shape function is not primitive, then either use the
143  // component of the first nonzero component corresponding
144  // to this shape function (if the component is not specified
145  // in the component_mask), or the first component of this block
146  // that is listed in the component_mask (if the block this
147  // component corresponds to is indeed specified in the component
148  // mask)
149  {
150  const unsigned int first_comp =
152 
153  if ((fe.get_nonzero_components(i) & component_mask)
154  .n_selected_components(fe.n_components()) == 0)
155  local_component_association[i] = first_comp;
156  else
157  // pick the component selected. we know from the previous 'if'
158  // that within the components that are nonzero for this
159  // shape function there must be at least one for which the
160  // mask is true, so we will for sure run into the break()
161  // at one point
162  for (unsigned int c = first_comp; c < fe.n_components(); ++c)
163  if (component_mask[c] == true)
164  {
165  local_component_association[i] = c;
166  break;
167  }
168  }
169 
170  Assert(std::find(local_component_association.begin(),
171  local_component_association.end(),
172  static_cast<unsigned char>(-1)) ==
173  local_component_association.end(),
174  ExcInternalError());
175 
176  return local_component_association;
177  }
178 
179 
180  // this internal function assigns to each dof the respective component
181  // of the vector system.
182  //
183  // the output array dofs_by_component lists for each dof the
184  // corresponding vector component. if the DoFHandler is based on a
185  // parallel distributed triangulation then the output array is index by
186  // dof.locally_owned_dofs().index_within_set(indices[i])
187  //
188  // if an element is non-primitive then we assign to each degree of
189  // freedom the following component:
190  // - if the nonzero components that belong to a shape function are not
191  // selected in the component_mask, then the shape function is assigned
192  // to the first nonzero vector component that corresponds to this
193  // shape function
194  // - otherwise, the shape function is assigned the first component selected
195  // in the component_mask that corresponds to this shape function
196  template <int dim, int spacedim>
197  void
199  const ComponentMask & component_mask,
200  std::vector<unsigned char> &dofs_by_component)
201  {
202  const ::hp::FECollection<dim, spacedim> &fe_collection =
203  dof.get_fe_collection();
204  Assert(fe_collection.n_components() < 256, ExcNotImplemented());
205  Assert(dofs_by_component.size() == dof.n_locally_owned_dofs(),
206  ExcDimensionMismatch(dofs_by_component.size(),
207  dof.n_locally_owned_dofs()));
208 
209  // next set up a table for the degrees of freedom on each of the
210  // cells (regardless of the fact whether it is listed in the
211  // component_select argument or not)
212  //
213  // for each element 'f' of the FECollection,
214  // local_component_association[f][d] then returns the vector
215  // component that degree of freedom 'd' belongs to
216  std::vector<std::vector<unsigned char>> local_component_association(
217  fe_collection.size());
218  for (unsigned int f = 0; f < fe_collection.size(); ++f)
219  {
220  const FiniteElement<dim, spacedim> &fe = fe_collection[f];
221  local_component_association[f] =
222  get_local_component_association(fe, component_mask);
223  }
224 
225  // then loop over all cells and do the work
226  std::vector<types::global_dof_index> indices;
227  for (const auto &c : 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(),
282  ExcInternalError());
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,
339  fe_collection.component_mask(FEValuesExtractors::Scalar(component)),
340  component_dofs);
341 
342  for (unsigned int i = 0; i < dof_data.size(); ++i)
343  if (component_dofs[i] == static_cast<unsigned char>(component))
344  dof_data(i) = 0;
345  }
346 
347  // count how often we have added a value in the sum for each dof
348  std::vector<unsigned char> touch_count(dof_handler.n_dofs(), 0);
349 
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),
383  ExcInternalError());
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
394  const ComponentMask & component_mask,
395  std::vector<bool> & selected_dofs)
396  {
397  Assert(component_mask.represents_n_components(
398  dof.get_fe_collection().n_components()),
399  ExcMessage(
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
408  if (component_mask.n_selected_components(
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  }
414  else if (component_mask.n_selected_components(
415  dof.get_fe_collection().n_components()) ==
416  dof.get_fe_collection().n_components())
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());
429  internal::get_component_association(dof, component_mask, dofs_by_component);
430 
431  for (types::global_dof_index i = 0; i < dof.n_locally_owned_dofs(); ++i)
432  if (component_mask[dofs_by_component[i]] == true)
433  selected_dofs[i] = true;
434  }
435 
436 
437 
438  template <int dim, int spacedim>
439  IndexSet
441  const ComponentMask & component_mask)
442  {
443  Assert(component_mask.represents_n_components(
444  dof.get_fe_collection().n_components()),
445  ExcMessage(
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
451  if (component_mask.n_selected_components(
452  dof.get_fe_collection().n_components()) == 0)
453  return IndexSet(dof.n_dofs());
454  else if (component_mask.n_selected_components(
455  dof.get_fe_collection().n_components()) ==
456  dof.get_fe_collection().n_components())
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());
462  internal::get_component_association(dof, component_mask, dofs_by_component);
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)
468  if (component_mask[dofs_by_component[i]] == true)
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());
473  result.add_indices(selected_dofs.begin(), selected_dofs.end());
474  return result;
475  }
476 
477 
478 
479  template <int dim, int spacedim>
480  void
482  const BlockMask & block_mask,
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>(
487  dof, dof.get_fe_collection().component_mask(block_mask), selected_dofs);
488  }
489 
490 
491 
492  template <int dim, int spacedim>
493  IndexSet
495  const BlockMask & block_mask)
496  {
497  // simply forward to the function that works based on a component mask
498  return extract_dofs<dim, spacedim>(
499  dof, dof.get_fe_collection().component_mask(block_mask));
500  }
501 
502 
503 
504  template <int dim, int spacedim>
505  std::vector<IndexSet>
507  const ComponentMask &component_mask)
508  {
509  const auto n_comps = dof.get_fe_collection().n_components();
510  Assert(component_mask.represents_n_components(n_comps),
511  ExcMessage(
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());
520  internal::get_component_association(dof, component_mask, dofs_by_component);
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];
527  if (component_mask[comp_i])
528  index_per_comp[comp_i].add_index(
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,
542  const ComponentMask & component_mask,
543  std::vector<bool> & selected_dofs)
544  {
545  const FiniteElement<dim, spacedim> &fe = dof.get_fe();
546 
547  Assert(component_mask.represents_n_components(
548  dof.get_fe_collection().n_components()),
549  ExcMessage(
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
557  if (component_mask.n_selected_components(
558  dof.get_fe_collection().n_components()) == 0)
559  {
560  std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false);
561  return;
562  }
563  else if (component_mask.n_selected_components(
564  dof.get_fe_collection().n_components()) ==
565  dof.get_fe_collection().n_components())
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 =
577  internal::get_local_component_association(fe, component_mask);
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)
580  local_selected_dofs[i] = component_mask[local_component_asssociation[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,
599  const BlockMask & block_mask,
600  std::vector<bool> & selected_dofs)
601  {
602  // simply defer to the other extract_level_dofs() function
603  extract_level_dofs(level,
604  dof,
605  dof.get_fe().component_mask(block_mask),
606  selected_dofs);
607  }
608 
609 
610 
611  template <int dim, int spacedim>
612  void
614  const ComponentMask & component_mask,
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),
621  ExcMessage(
622  "This function can not be used with distributed triangulations. "
623  "See the documentation for more information."));
624 
625  IndexSet indices;
626  extract_boundary_dofs(dof_handler, component_mask, indices, boundary_ids);
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
641  const ComponentMask & component_mask,
642  IndexSet & selected_dofs,
643  const std::set<types::boundary_id> &boundary_ids)
644  {
645  // Simply forward to the other function
646  selected_dofs =
647  extract_boundary_dofs(dof_handler, component_mask, boundary_ids);
648  }
649 
650 
651 
652  template <int dim, int spacedim>
653  IndexSet
655  const ComponentMask & component_mask,
656  const std::set<types::boundary_id> &boundary_ids)
657  {
658  Assert(component_mask.represents_n_components(
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 =
674  ((component_mask.represents_the_all_selected_mask() == false) ||
675  (component_mask.n_selected_components(
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)
718  selected_dofs.add_index(face_dof_indices[i]);
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  {
756  if (component_mask[fe.face_system_to_component_index(
757  i, face)
758  .first] == true)
759  selected_dofs.add_index(face_dof_indices[i]);
760  }
761  else // not primitive
762  {
763  const unsigned int first_nonzero_comp =
764  fe.get_nonzero_components(cell_index)
766  Assert(first_nonzero_comp < fe.n_components(),
767  ExcInternalError());
768 
769  if (component_mask[first_nonzero_comp] == true)
770  selected_dofs.add_index(face_dof_indices[i]);
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,
784  const ComponentMask & component_mask,
785  std::vector<bool> & selected_dofs,
786  const std::set<types::boundary_id> &boundary_ids)
787  {
788  Assert(component_mask.represents_n_components(
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 =
802  (component_mask.represents_the_all_selected_mask() == false);
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]] =
842  (component_mask[fe.system_to_component_index(i)
843  .first] == true);
844  else // not primitive
845  {
846  const unsigned int first_nonzero_comp =
849  Assert(first_nonzero_comp < fe.n_components(),
850  ExcInternalError());
851 
852  selected_dofs[cell_dof_indices[i]] =
853  (component_mask[first_nonzero_comp] == true);
854  }
855  }
856  }
857  }
858  }
859 
860 
861 
862  template <int dim, int spacedim, typename number>
863  IndexSet
865  const DoFHandler<dim, spacedim> &dof_handler,
866  const std::function<
868  & predicate,
869  const AffineConstraints<number> &cm)
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());
953  support_set.add_indices(predicate_dofs.begin(), predicate_dofs.end());
954  support_set.compress();
955 
956  IndexSet halo_set(dof_handler.n_dofs());
957  halo_set.add_indices(dofs_with_support_on_halo_cells.begin(),
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>
974  IndexSet
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>
984  IndexSet
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,
995  // see there for more information
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)
1008  selected_dofs.add_index(
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)
1018  selected_dofs.add_index(
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);
1057  selected_dofs.add_indices(ldi.begin(), ldi.end());
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)
1065  unconstrained_dofs.add_index(
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
1120  // there for more information
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 
1165  dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
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 =
1189  filter_iterators(dof_handler.cell_iterators_on_level(level),
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 
1201  dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
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());
1242  dof_set.add_indices(dofs_on_ghosts.begin(),
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());
1293  dof_set.add_indices(dofs_on_ghosts.begin(),
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
1305  const ComponentMask & component_mask,
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();
1317  Assert(component_mask.represents_n_components(n_components),
1318  ExcDimensionMismatch(n_components, component_mask.size()));
1319 
1320  std::vector<unsigned char> dofs_by_component(
1321  dof_handler.n_locally_owned_dofs());
1323  component_mask,
1324  dofs_by_component);
1325  unsigned int n_selected_dofs = 0;
1326  for (unsigned int i = 0; i < n_components; ++i)
1327  if (component_mask[i] == true)
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)
1337  if (component_mask[dofs_by_component[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)
1358  if (component_mask[data.second[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];
1384  if (component_mask[comp])
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()),
1465  ExcInternalError());
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  {
1482  index_sets[this_subdomain].add_range(i_min, index);
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  {
1491  index_sets[this_subdomain].add_index(i_min);
1492  }
1493 
1494  // otherwise there are at least two different indices
1495  else
1496  {
1497  index_sets[this_subdomain].add_range(i_min,
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 =
1530  locally_owned_dofs_per_subdomain(dof_handler);
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.
1536  for (::types::subdomain_id subdomain_id = 0;
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 
1573  dof_set[subdomain_id].add_indices(
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());
1616  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
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(),
1622  ExcInternalError());
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(),
1674  ExcInternalError());
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.
1704  Assert((dof_handler.get_triangulation().locally_owned_subdomain() ==
1706  (subdomain ==
1707  dof_handler.get_triangulation().locally_owned_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
1747  index_set.add_indices(subdomain_indices.begin(), subdomain_indices.end());
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  {
1815  const FiniteElement<dim, spacedim> &base = fe.base_element(b);
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(),
1865  ExcNotImplemented());
1866  Assert(fe_collection[fe].n_base_elements() ==
1867  fe_collection[0].n_base_elements(),
1868  ExcNotImplemented());
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(),
1873  ExcNotImplemented());
1874  Assert(fe_collection[fe].base_element(b).n_base_elements() ==
1875  fe_collection[0].base_element(b).n_base_elements(),
1876  ExcNotImplemented());
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());
1985  ComponentMask(),
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()),
2006  ExcInternalError());
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
2228  const hp::MappingCollection<dim, spacedim> & mapping,
2229  const DoFHandler<dim, spacedim> & dof_handler,
2230  std::map<types::global_dof_index, Point<spacedim>> &support_points,
2231  const ComponentMask & in_mask)
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(),
2243  q_coll_dummy.push_back(Quadrature<dim>(
2244  fe_collection[fe_index].get_unit_support_points()));
2245  }
2246 
2247  // Take care of components
2248  const ComponentMask mask =
2249  (in_mask.size() == 0 ?
2250  ComponentMask(fe_collection.n_components(), true) :
2251  in_mask);
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 =
2281  fe_values.get_quadrature_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
2289  if (mask[dof_comp])
2290  support_points[local_dof_indices[i]] = points[i];
2291  }
2292  }
2293  }
2294 
2295 
2296  template <int dim, int spacedim>
2297  void
2299  const hp::MappingCollection<dim, spacedim> &mapping,
2300  const DoFHandler<dim, spacedim> & dof_handler,
2301  std::vector<Point<spacedim>> & support_points,
2302  const ComponentMask & mask)
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,
2309  mask);
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(),
2316  ExcInternalError());
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,
2329  const ComponentMask & mask)
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. "
2337  "See the documentation for more information."));
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,
2346  mask);
2347  }
2348 
2349 
2350  template <int dim, int spacedim>
2351  void
2353  const hp::MappingCollection<dim, spacedim> &mapping,
2354  const DoFHandler<dim, spacedim> & dof_handler,
2355  std::vector<Point<spacedim>> & support_points,
2356  const ComponentMask & mask)
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. "
2364  "See the documentation for more information."));
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,
2371  mask);
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,
2381  const ComponentMask & mask)
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,
2392  mask);
2393  }
2394 
2395 
2396  template <int dim, int spacedim>
2397  void
2399  const hp::MappingCollection<dim, spacedim> & mapping,
2400  const DoFHandler<dim, spacedim> & dof_handler,
2401  std::map<types::global_dof_index, Point<spacedim>> &support_points,
2402  const ComponentMask & mask)
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,
2411  mask);
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)
2547  block_list.add(i, indices[j] - offset);
2548  else
2549  {
2550  if (selected_dofs[j])
2551  block_list.add(i, indices[j] - offset);
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])
2593  block_list.add(0, indices[j]);
2594  }
2595  else
2596  {
2597  for (const auto index : indices)
2598  block_list.add(0, index);
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])
2666  block_list.add(block, indices[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();
2684  BlockMask exclude_boundary_dofs = BlockMask(n_blocks, interior_only);
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,
2700  const BlockMask & exclude_boundary_dofs,
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
2728  // about them
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(
2812  fe.face_to_cell_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])
2820  block_list.add(block, indices[j]);
2821  }
2822  else
2823  {
2824  for (const auto index : indices)
2825  block_list.add(block, index);
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)),
2903  ExcInternalError());
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 
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 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
unsigned int n_active_cells() const
Definition: tria.cc:12644
IndexSet dof_indices_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
Definition: dof_tools.cc:1698
void get_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::subdomain_id > &subdomain)
Definition: dof_tools.cc:1585
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
const types::subdomain_id invalid_subdomain_id
Definition: types.h:276
cell_iterator begin(const unsigned int level=0) const
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 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 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
static ::ExceptionBase & ExcIO()
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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
bool operator()(const Point< dim, Number > &lhs, const Point< dim, Number > &rhs) const
Definition: dof_tools.cc:85
cell_iterator end() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
std::vector< IndexSet > locally_relevant_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1509
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1703
bool represents_n_components(const unsigned int n) const
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:86
unsigned int n_dofs_per_vertex() const
bool represents_the_all_selected_mask() const
typename ActiveSelector::CellAccessor cell_accessor
Definition: dof_handler.h:334
void map_dof_to_boundary_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::global_dof_index > &mapping)
Definition: dof_tools.cc:2136
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
Transformed quadrature points.
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
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
IteratorRange< FilteredIterator< BaseIterator > > filter_iterators(IteratorRange< BaseIterator > i, const Predicate &p)
bool is_primitive() const
Definition: fe.h:3302
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
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Definition: point.h:110
std::vector< unsigned char > get_local_component_association(const FiniteElement< dim, spacedim > &fe, const ComponentMask &component_mask)
Definition: dof_tools.cc:127
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1414
unsigned int n_blocks() const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
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 extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1252
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< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
active_cell_iterator begin_active(const unsigned int level=0) const
const std::vector< Point< spacedim > > & get_quadrature_points() const
unsigned int component_to_block_index(const unsigned int component) const
Definition: fe.cc:384
size_type n_constraints() const
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)
std::string compress(const std::string &input)
Definition: utilities.cc:392
unsigned int subdomain_id
Definition: types.h:43
unsigned int n_dofs_per_line() const
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
#define Assert(cond, exc)
Definition: exceptions.h:1465
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:3123
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
types::global_dof_index n_dofs() const
Abstract base class for mapping classes.
Definition: mapping.h:303
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:3282
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
void add(const size_type i, const size_type j)
std::vector< IndexSet > locally_owned_dofs_per_component(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &components=ComponentMask())
Definition: dof_tools.cc:506
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
unsigned int level
Definition: grid_out.cc:4590
VectorType::value_type * end(VectorType &V)
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:502
unsigned int size() const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe.cc:1402
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
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
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_hanging_node_dofs(const DoFHandler< dim, spacedim > &dof_handler, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:1079
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
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
Definition: fe.cc:568
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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 reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda >> &cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition: fe_values.cc:293
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
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Point< 2 > first
Definition: grid_out.cc:4587
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
const IndexSet & locally_owned_dofs() const
void push_back(const Quadrature< dim_in > &new_quadrature)
Definition: q_collection.h:223
unsigned int count_dofs_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
Definition: dof_tools.cc:1681
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
void extract_locally_active_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set, const unsigned int level)
Definition: dof_tools.cc:1174
unsigned int global_dof_index
Definition: types.h:76
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
void extract_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:393
bool has_hp_capabilities() const
void compress() const
Definition: index_set.h:1642
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
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3100
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1746
unsigned int n_components() const
unsigned int n_dofs_per_cell() const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
Definition: fe.h:3254
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
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
VectorType::value_type * begin(VectorType &V)
unsigned int n_selected_blocks(const unsigned int overall_number_of_blocks=numbers::invalid_unsigned_int) const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
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
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
types::global_dof_index n_boundary_dofs() const
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1399
Definition: fe.h:38
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
unsigned int count_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
Definition: dof_tools.cc:2848
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcNotImplemented()
bool is_element(const size_type index) const
Definition: index_set.h:1765
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
const types::boundary_id internal_face_boundary_id
Definition: types.h:255
void extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1210
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
const types::global_dof_index invalid_dof_index
Definition: types.h:211
unsigned int size() const
unsigned int n_base_elements() const
Definition: fe.h:3114
void get_block_association(const DoFHandler< dim, spacedim > &dof, std::vector< unsigned char > &dofs_by_block)
Definition: dof_tools.cc:251
void extract_locally_active_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1139
unsigned int n_components() const
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1880
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
unsigned int cell_index
Definition: grid_tools.cc:1092
types::global_dof_index n_locally_owned_dofs() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
Definition: fe.h:3160
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
Definition: fe.cc:1154
static ::ExceptionBase & ExcInternalError()
const ::FEValues< dim, spacedim > & get_present_fe_values() const
Definition: fe_values.h:693