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