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