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