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