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