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