Reference documentation for deal.II version Git 99efdf013c 2021-02-27 20:52:41 -0500
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools_dof_handlers.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/point.h>
18 #include <deal.II/base/tensor.h>
19 
23 
26 
27 #include <deal.II/fe/mapping_q1.h>
29 
32 #include <deal.II/grid/tria.h>
35 
36 #include <deal.II/hp/dof_handler.h>
38 
40 
41 #include <algorithm>
42 #include <array>
43 #include <cmath>
44 #include <list>
45 #include <map>
46 #include <numeric>
47 #include <set>
48 #include <vector>
49 
50 
52 
53 namespace GridTools
54 {
55  template <int dim, template <int, int> class MeshType, int spacedim>
56  unsigned int
57  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
58  const Point<spacedim> & p,
59  const std::vector<bool> & marked_vertices)
60  {
61  // first get the underlying
62  // triangulation from the
63  // mesh and determine vertices
64  // and used vertices
66 
67  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
68 
69  Assert(tria.get_vertices().size() == marked_vertices.size() ||
70  marked_vertices.size() == 0,
71  ExcDimensionMismatch(tria.get_vertices().size(),
72  marked_vertices.size()));
73 
74  // If p is an element of marked_vertices,
75  // and q is that of used_Vertices,
76  // the vector marked_vertices does NOT
77  // contain unused vertices if p implies q.
78  // I.e., if p is true q must be true
79  // (if p is false, q could be false or true).
80  // p implies q logic is encapsulated in ~p|q.
81  Assert(
82  marked_vertices.size() == 0 ||
83  std::equal(marked_vertices.begin(),
84  marked_vertices.end(),
85  tria.get_used_vertices().begin(),
86  [](bool p, bool q) { return !p || q; }),
87  ExcMessage(
88  "marked_vertices should be a subset of used vertices in the triangulation "
89  "but marked_vertices contains one or more vertices that are not used vertices!"));
90 
91  // In addition, if a vector bools
92  // is specified (marked_vertices)
93  // marking all the vertices which
94  // could be the potentially closest
95  // vertex to the point, use it instead
96  // of used vertices
97  const std::vector<bool> &used = (marked_vertices.size() == 0) ?
98  tria.get_used_vertices() :
99  marked_vertices;
100 
101  // At the beginning, the first
102  // used vertex is the closest one
103  std::vector<bool>::const_iterator first =
104  std::find(used.begin(), used.end(), true);
105 
106  // Assert that at least one vertex
107  // is actually used
108  Assert(first != used.end(), ExcInternalError());
109 
110  unsigned int best_vertex = std::distance(used.begin(), first);
111  double best_dist = (p - vertices[best_vertex]).norm_square();
112 
113  // For all remaining vertices, test
114  // whether they are any closer
115  for (unsigned int j = best_vertex + 1; j < vertices.size(); j++)
116  if (used[j])
117  {
118  double dist = (p - vertices[j]).norm_square();
119  if (dist < best_dist)
120  {
121  best_vertex = j;
122  best_dist = dist;
123  }
124  }
125 
126  return best_vertex;
127  }
128 
129 
130 
131  template <int dim, template <int, int> class MeshType, int spacedim>
132  unsigned int
134  const MeshType<dim, spacedim> &mesh,
135  const Point<spacedim> & p,
136  const std::vector<bool> & marked_vertices)
137  {
138  // Take a shortcut in the simple case.
139  if (mapping.preserves_vertex_locations() == true)
140  return find_closest_vertex(mesh, p, marked_vertices);
141 
142  // first get the underlying
143  // triangulation from the
144  // mesh and determine vertices
145  // and used vertices
146  const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();
147 
148  auto vertices = extract_used_vertices(tria, mapping);
149 
150  Assert(tria.get_vertices().size() == marked_vertices.size() ||
151  marked_vertices.size() == 0,
152  ExcDimensionMismatch(tria.get_vertices().size(),
153  marked_vertices.size()));
154 
155  // If p is an element of marked_vertices,
156  // and q is that of used_Vertices,
157  // the vector marked_vertices does NOT
158  // contain unused vertices if p implies q.
159  // I.e., if p is true q must be true
160  // (if p is false, q could be false or true).
161  // p implies q logic is encapsulated in ~p|q.
162  Assert(
163  marked_vertices.size() == 0 ||
164  std::equal(marked_vertices.begin(),
165  marked_vertices.end(),
166  tria.get_used_vertices().begin(),
167  [](bool p, bool q) { return !p || q; }),
168  ExcMessage(
169  "marked_vertices should be a subset of used vertices in the triangulation "
170  "but marked_vertices contains one or more vertices that are not used vertices!"));
171 
172  // Remove from the map unwanted elements.
173  if (marked_vertices.size())
174  for (auto it = vertices.begin(); it != vertices.end();)
175  {
176  if (marked_vertices[it->first] == false)
177  {
178  vertices.erase(it++);
179  }
180  else
181  {
182  ++it;
183  }
184  }
185 
186  return find_closest_vertex(vertices, p);
187  }
188 
189 
190 
191  template <int dim, template <int, int> class MeshType, int spacedim>
192 #ifndef _MSC_VER
193  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
194 #else
195  std::vector<
196  typename ::internal::
197  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
198 #endif
199  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &mesh,
200  const unsigned int vertex)
201  {
202  // make sure that the given vertex is
203  // an active vertex of the underlying
204  // triangulation
205  AssertIndexRange(vertex, mesh.get_triangulation().n_vertices());
206  Assert(mesh.get_triangulation().get_used_vertices()[vertex],
207  ExcVertexNotUsed(vertex));
208 
209  // use a set instead of a vector
210  // to ensure that cells are inserted only
211  // once
212  std::set<typename ::internal::
213  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
214  adjacent_cells;
215 
216  typename ::internal::
217  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
218  cell = mesh.begin_active(),
219  endc = mesh.end();
220 
221  // go through all active cells and look if the vertex is part of that cell
222  //
223  // in 1d, this is all we need to care about. in 2d/3d we also need to worry
224  // that the vertex might be a hanging node on a face or edge of a cell; in
225  // this case, we would want to add those cells as well on whose faces the
226  // vertex is located but for which it is not a vertex itself.
227  //
228  // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging
229  // node can only be in the middle of a face and we can query the neighboring
230  // cell from the current cell. on the other hand, in 3d a hanging node
231  // vertex can also be on an edge but there can be many other cells on
232  // this edge and we can not access them from the cell we are currently
233  // on.
234  //
235  // so, in the 3d case, if we run the algorithm as in 2d, we catch all
236  // those cells for which the vertex we seek is on a *subface*, but we
237  // miss the case of cells for which the vertex we seek is on a
238  // sub-edge for which there is no corresponding sub-face (because the
239  // immediate neighbor behind this face is not refined), see for example
240  // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
241  // haven't yet found the vertex for the current cell we also need to
242  // look at the mid-points of edges
243  //
244  // as a final note, deciding whether a neighbor is actually coarser is
245  // simple in the case of isotropic refinement (we just need to look at
246  // the level of the current and the neighboring cell). however, this
247  // isn't so simple if we have used anisotropic refinement since then
248  // the level of a cell is not indicative of whether it is coarser or
249  // not than the current cell. ultimately, we want to add all cells on
250  // which the vertex is, independent of whether they are coarser or
251  // finer and so in the 2d case below we simply add *any* *active* neighbor.
252  // in the worst case, we add cells multiple times to the adjacent_cells
253  // list, but std::set throws out those cells already entered
254  for (; cell != endc; ++cell)
255  {
256  for (const unsigned int v : cell->vertex_indices())
257  if (cell->vertex_index(v) == vertex)
258  {
259  // OK, we found a cell that contains
260  // the given vertex. We add it
261  // to the list.
262  adjacent_cells.insert(cell);
263 
264  // as explained above, in 2+d we need to check whether
265  // this vertex is on a face behind which there is a
266  // (possibly) coarser neighbor. if this is the case,
267  // then we need to also add this neighbor
268  if (dim >= 2)
269  for (const auto face :
270  cell->reference_cell().faces_for_given_vertex(v))
271  if (!cell->at_boundary(face) &&
272  cell->neighbor(face)->is_active())
273  {
274  // there is a (possibly) coarser cell behind a
275  // face to which the vertex belongs. the
276  // vertex we are looking at is then either a
277  // vertex of that coarser neighbor, or it is a
278  // hanging node on one of the faces of that
279  // cell. in either case, it is adjacent to the
280  // vertex, so add it to the list as well (if
281  // the cell was already in the list then the
282  // std::set makes sure that we get it only
283  // once)
284  adjacent_cells.insert(cell->neighbor(face));
285  }
286 
287  // in any case, we have found a cell, so go to the next cell
288  goto next_cell;
289  }
290 
291  // in 3d also loop over the edges
292  if (dim >= 3)
293  {
294  for (unsigned int e = 0; e < cell->n_lines(); ++e)
295  if (cell->line(e)->has_children())
296  // the only place where this vertex could have been
297  // hiding is on the mid-edge point of the edge we
298  // are looking at
299  if (cell->line(e)->child(0)->vertex_index(1) == vertex)
300  {
301  adjacent_cells.insert(cell);
302 
303  // jump out of this tangle of nested loops
304  goto next_cell;
305  }
306  }
307 
308  // in more than 3d we would probably have to do the same as
309  // above also for even lower-dimensional objects
310  Assert(dim <= 3, ExcNotImplemented());
311 
312  // move on to the next cell if we have found the
313  // vertex on the current one
314  next_cell:;
315  }
316 
317  // if this was an active vertex then there needs to have been
318  // at least one cell to which it is adjacent!
319  Assert(adjacent_cells.size() > 0, ExcInternalError());
320 
321  // return the result as a vector, rather than the set we built above
322  return std::vector<
323  typename ::internal::
324  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
325  adjacent_cells.begin(), adjacent_cells.end());
326  }
327 
328 
329 
330  namespace
331  {
332  template <int dim, template <int, int> class MeshType, int spacedim>
333  void
334  find_active_cell_around_point_internal(
335  const MeshType<dim, spacedim> &mesh,
336 #ifndef _MSC_VER
337  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>
338  &searched_cells,
339  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>
340  &adjacent_cells)
341 #else
342  std::set<
343  typename ::internal::
344  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
345  &searched_cells,
346  std::set<
347  typename ::internal::
348  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
349  &adjacent_cells)
350 #endif
351  {
352 #ifndef _MSC_VER
353  using cell_iterator =
354  typename MeshType<dim, spacedim>::active_cell_iterator;
355 #else
356  using cell_iterator = typename ::internal::
357  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
358 #endif
359 
360  // update the searched cells
361  searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end());
362  // now we to collect all neighbors
363  // of the cells in adjacent_cells we
364  // have not yet searched.
365  std::set<cell_iterator> adjacent_cells_new;
366 
367  typename std::set<cell_iterator>::const_iterator cell =
368  adjacent_cells.begin(),
369  endc =
370  adjacent_cells.end();
371  for (; cell != endc; ++cell)
372  {
373  std::vector<cell_iterator> active_neighbors;
374  get_active_neighbors<MeshType<dim, spacedim>>(*cell,
375  active_neighbors);
376  for (unsigned int i = 0; i < active_neighbors.size(); ++i)
377  if (searched_cells.find(active_neighbors[i]) ==
378  searched_cells.end())
379  adjacent_cells_new.insert(active_neighbors[i]);
380  }
381  adjacent_cells.clear();
382  adjacent_cells.insert(adjacent_cells_new.begin(),
383  adjacent_cells_new.end());
384  if (adjacent_cells.size() == 0)
385  {
386  // we haven't found any other cell that would be a
387  // neighbor of a previously found cell, but we know
388  // that we haven't checked all cells yet. that means
389  // that the domain is disconnected. in that case,
390  // choose the first previously untouched cell we
391  // can find
392  cell_iterator it = mesh.begin_active();
393  for (; it != mesh.end(); ++it)
394  if (searched_cells.find(it) == searched_cells.end())
395  {
396  adjacent_cells.insert(it);
397  break;
398  }
399  }
400  }
401  } // namespace
402 
403 
404 
405  template <int dim, template <int, int> class MeshType, int spacedim>
406 #ifndef _MSC_VER
407  typename MeshType<dim, spacedim>::active_cell_iterator
408 #else
409  typename ::internal::
410  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
411 #endif
412  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
413  const Point<spacedim> & p,
414  const std::vector<bool> & marked_vertices,
415  const double tolerance)
416  {
417  return find_active_cell_around_point<dim, MeshType, spacedim>(
418  get_default_linear_mapping(mesh.get_triangulation()),
419  mesh,
420  p,
421  marked_vertices,
422  tolerance)
423  .first;
424  }
425 
426 
427 
428  template <int dim, template <int, int> class MeshType, int spacedim>
429 #ifndef _MSC_VER
430  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
431 #else
432  std::pair<typename ::internal::
433  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
434  Point<dim>>
435 #endif
437  const MeshType<dim, spacedim> &mesh,
438  const Point<spacedim> & p,
439  const std::vector<bool> & marked_vertices,
440  const double tolerance)
441  {
442  using active_cell_iterator = typename ::internal::
443  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
444 
445  // The best distance is set to the
446  // maximum allowable distance from
447  // the unit cell; we assume a
448  // max. deviation of the given tolerance
449  double best_distance = tolerance;
450  int best_level = -1;
451  std::pair<active_cell_iterator, Point<dim>> best_cell;
452 
453  // Find closest vertex and determine
454  // all adjacent cells
455  std::vector<active_cell_iterator> adjacent_cells_tmp =
457  mesh, find_closest_vertex(mapping, mesh, p, marked_vertices));
458 
459  // Make sure that we have found
460  // at least one cell adjacent to vertex.
461  Assert(adjacent_cells_tmp.size() > 0, ExcInternalError());
462 
463  // Copy all the cells into a std::set
464  std::set<active_cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(),
465  adjacent_cells_tmp.end());
466  std::set<active_cell_iterator> searched_cells;
467 
468  // Determine the maximal number of cells
469  // in the grid.
470  // As long as we have not found
471  // the cell and have not searched
472  // every cell in the triangulation,
473  // we keep on looking.
474  const unsigned int n_active_cells =
475  mesh.get_triangulation().n_active_cells();
476  bool found = false;
477  unsigned int cells_searched = 0;
478  while (!found && cells_searched < n_active_cells)
479  {
480  typename std::set<active_cell_iterator>::const_iterator
481  cell = adjacent_cells.begin(),
482  endc = adjacent_cells.end();
483  for (; cell != endc; ++cell)
484  {
485  if ((*cell)->is_artificial() == false)
486  {
487  try
488  {
489  const Point<dim> p_cell =
490  mapping.transform_real_to_unit_cell(*cell, p);
491 
492  // calculate the infinity norm of
493  // the distance vector to the unit cell.
494  const double dist =
496 
497  // We compare if the point is inside the
498  // unit cell (or at least not too far
499  // outside). If it is, it is also checked
500  // that the cell has a more refined state
501  if ((dist < best_distance) ||
502  ((dist == best_distance) &&
503  ((*cell)->level() > best_level)))
504  {
505  found = true;
506  best_distance = dist;
507  best_level = (*cell)->level();
508  best_cell = std::make_pair(*cell, p_cell);
509  }
510  }
511  catch (
512  typename MappingQGeneric<dim,
513  spacedim>::ExcTransformationFailed &)
514  {
515  // ok, the transformation
516  // failed presumably
517  // because the point we
518  // are looking for lies
519  // outside the current
520  // cell. this means that
521  // the current cell can't
522  // be the cell around the
523  // point, so just ignore
524  // this cell and move on
525  // to the next
526  }
527  }
528  }
529 
530  // update the number of cells searched
531  cells_searched += adjacent_cells.size();
532 
533  // if the user provided a custom mask for vertices,
534  // terminate the search without trying to expand the search
535  // to all cells of the triangulation, as done below.
536  if (marked_vertices.size() > 0)
537  cells_searched = n_active_cells;
538 
539  // if we have not found the cell in
540  // question and have not yet searched every
541  // cell, we expand our search to
542  // all the not already searched neighbors of
543  // the cells in adjacent_cells. This is
544  // what find_active_cell_around_point_internal
545  // is for.
546  if (!found && cells_searched < n_active_cells)
547  {
548  find_active_cell_around_point_internal<dim, MeshType, spacedim>(
549  mesh, searched_cells, adjacent_cells);
550  }
551  }
552 
553  AssertThrow(best_cell.first.state() == IteratorState::valid,
554  ExcPointNotFound<spacedim>(p));
555 
556  return best_cell;
557  }
558 
559 
560 
561  template <int dim, template <int, int> class MeshType, int spacedim>
562 #ifndef _MSC_VER
563  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
564  Point<dim>>>
565 #else
566  std::vector<std::pair<
567  typename ::internal::
568  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
569  Point<dim>>>
570 #endif
572  const MeshType<dim, spacedim> &mesh,
573  const Point<spacedim> & p,
574  const double tolerance,
575  const std::vector<bool> &marked_vertices)
576  {
577  try
578  {
579  const auto cell_and_point = find_active_cell_around_point(
580  mapping, mesh, p, marked_vertices, tolerance);
581 
583  mapping, mesh, p, tolerance, cell_and_point);
584  }
585  catch (ExcPointNotFound<spacedim> &)
586  {}
587 
588  return {};
589  }
590 
591 
592 
593  template <int dim, template <int, int> class MeshType, int spacedim>
594 #ifndef _MSC_VER
595  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
596  Point<dim>>>
597 #else
598  std::vector<std::pair<
599  typename ::internal::
600  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
601  Point<dim>>>
602 #endif
604  const Mapping<dim, spacedim> & mapping,
605  const MeshType<dim, spacedim> &mesh,
606  const Point<spacedim> & p,
607  const double tolerance,
608  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
609  Point<dim>> & first_cell)
610  {
611  std::vector<
612  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
613  Point<dim>>>
614  cells_and_points;
615 
616  // insert the fist cell and point into the vector
617  cells_and_points.push_back(first_cell);
618 
619  // check if the given point is on the surface of the unit cell. if yes,
620  // need to find all neighbors
621  const Point<dim> unit_point = cells_and_points.front().second;
622  const auto my_cell = cells_and_points.front().first;
623  Tensor<1, dim> distance_to_center;
624  unsigned int n_dirs_at_threshold = 0;
625  unsigned int last_point_at_threshold = numbers::invalid_unsigned_int;
626  for (unsigned int d = 0; d < dim; ++d)
627  {
628  distance_to_center[d] = std::abs(unit_point[d] - 0.5);
629  if (distance_to_center[d] > 0.5 - tolerance)
630  {
631  ++n_dirs_at_threshold;
632  last_point_at_threshold = d;
633  }
634  }
635 
636  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
637  cells_to_add;
638  // point is within face -> only need neighbor
639  if (n_dirs_at_threshold == 1)
640  {
641  unsigned int neighbor_index =
642  2 * last_point_at_threshold +
643  (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
644  if (!my_cell->at_boundary(neighbor_index))
645  cells_to_add.push_back(my_cell->neighbor(neighbor_index));
646  }
647  // corner point -> use all neighbors
648  else if (n_dirs_at_threshold == dim)
649  {
650  unsigned int local_vertex_index = 0;
651  for (unsigned int d = 0; d < dim; ++d)
652  local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
653  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
655  mesh, my_cell->vertex_index(local_vertex_index));
656  for (const auto &cell : cells)
657  if (cell != my_cell)
658  cells_to_add.push_back(cell);
659  }
660  // point on line in 3D: We cannot simply take the intersection between
661  // the two vertices of cells because of hanging nodes. So instead we
662  // list the vertices around both points and then select the
663  // appropriate cells according to the result of read_to_unit_cell
664  // below.
665  else if (n_dirs_at_threshold == 2)
666  {
667  std::pair<unsigned int, unsigned int> vertex_indices[3];
668  unsigned int count_vertex_indices = 0;
669  unsigned int free_direction = numbers::invalid_unsigned_int;
670  for (unsigned int d = 0; d < dim; ++d)
671  {
672  if (distance_to_center[d] > 0.5 - tolerance)
673  {
674  vertex_indices[count_vertex_indices].first = d;
675  vertex_indices[count_vertex_indices].second =
676  unit_point[d] > 0.5 ? 1 : 0;
677  count_vertex_indices++;
678  }
679  else
680  free_direction = d;
681  }
682 
683  AssertDimension(count_vertex_indices, 2);
684  Assert(free_direction != numbers::invalid_unsigned_int,
685  ExcInternalError());
686 
687  const unsigned int first_vertex =
688  (vertex_indices[0].second << vertex_indices[0].first) +
689  (vertex_indices[1].second << vertex_indices[1].first);
690  for (unsigned int d = 0; d < 2; ++d)
691  {
692  auto tentative_cells = find_cells_adjacent_to_vertex(
693  mesh,
694  my_cell->vertex_index(first_vertex + (d << free_direction)));
695  for (const auto &cell : tentative_cells)
696  {
697  bool cell_not_yet_present = true;
698  for (const auto &other_cell : cells_to_add)
699  if (cell == other_cell)
700  {
701  cell_not_yet_present = false;
702  break;
703  }
704  if (cell_not_yet_present)
705  cells_to_add.push_back(cell);
706  }
707  }
708  }
709 
710  const double original_distance_to_unit_cell =
712  for (const auto &cell : cells_to_add)
713  {
714  if (cell != my_cell)
715  try
716  {
717  const Point<dim> p_unit =
718  mapping.transform_real_to_unit_cell(cell, p);
720  original_distance_to_unit_cell + tolerance)
721  cells_and_points.emplace_back(cell, p_unit);
722  }
723  catch (typename Mapping<dim>::ExcTransformationFailed &)
724  {}
725  }
726 
727  std::sort(
728  cells_and_points.begin(),
729  cells_and_points.end(),
730  [](const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
731  Point<dim>> &a,
732  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
733  Point<dim>> &b) { return a.first < b.first; });
734 
735  return cells_and_points;
736  }
737 
738 
739 
740  template <class MeshType>
741  std::vector<typename MeshType::active_cell_iterator>
743  const MeshType &mesh,
744  const std::function<bool(const typename MeshType::active_cell_iterator &)>
745  &predicate)
746  {
747  std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
748  std::vector<bool> locally_active_vertices_on_subdomain(
749  mesh.get_triangulation().n_vertices(), false);
750 
751  // Find the cells for which the predicate is true
752  // These are the cells around which we wish to construct
753  // the halo layer
754  for (const auto &cell : mesh.active_cell_iterators())
755  if (predicate(cell)) // True predicate --> Part of subdomain
756  for (const auto v : cell->vertex_indices())
757  locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true;
758 
759  // Find the cells that do not conform to the predicate
760  // but share a vertex with the selected subdomain
761  // These comprise the halo layer
762  for (const auto &cell : mesh.active_cell_iterators())
763  if (!predicate(cell)) // False predicate --> Potential halo cell
764  for (const auto v : cell->vertex_indices())
765  if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
766  true)
767  {
768  active_halo_layer.push_back(cell);
769  break;
770  }
771 
772  return active_halo_layer;
773  }
774 
775 
776 
777  template <class MeshType>
778  std::vector<typename MeshType::cell_iterator>
780  const MeshType &mesh,
781  const std::function<bool(const typename MeshType::cell_iterator &)>
782  & predicate,
783  const unsigned int level)
784  {
785  std::vector<typename MeshType::cell_iterator> level_halo_layer;
786  std::vector<bool> locally_active_vertices_on_level_subdomain(
787  mesh.get_triangulation().n_vertices(), false);
788 
789  // Find the cells for which the predicate is true
790  // These are the cells around which we wish to construct
791  // the halo layer
792  for (typename MeshType::cell_iterator cell = mesh.begin(level);
793  cell != mesh.end(level);
794  ++cell)
795  if (predicate(cell)) // True predicate --> Part of subdomain
796  for (const unsigned int v : cell->vertex_indices())
797  locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
798  true;
799 
800  // Find the cells that do not conform to the predicate
801  // but share a vertex with the selected subdomain on that level
802  // These comprise the halo layer
803  for (typename MeshType::cell_iterator cell = mesh.begin(level);
804  cell != mesh.end(level);
805  ++cell)
806  if (!predicate(cell)) // False predicate --> Potential halo cell
807  for (const unsigned int v : cell->vertex_indices())
808  if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
809  v)] == true)
810  {
811  level_halo_layer.push_back(cell);
812  break;
813  }
814 
815  return level_halo_layer;
816  }
817 
818 
819  namespace
820  {
821  template <class MeshType>
822  bool
823  contains_locally_owned_cells(
824  const std::vector<typename MeshType::active_cell_iterator> &cells)
825  {
826  for (typename std::vector<
827  typename MeshType::active_cell_iterator>::const_iterator it =
828  cells.begin();
829  it != cells.end();
830  ++it)
831  {
832  if ((*it)->is_locally_owned())
833  return true;
834  }
835  return false;
836  }
837 
838  template <class MeshType>
839  bool
840  contains_artificial_cells(
841  const std::vector<typename MeshType::active_cell_iterator> &cells)
842  {
843  for (typename std::vector<
844  typename MeshType::active_cell_iterator>::const_iterator it =
845  cells.begin();
846  it != cells.end();
847  ++it)
848  {
849  if ((*it)->is_artificial())
850  return true;
851  }
852  return false;
853  }
854  } // namespace
855 
856 
857 
858  template <class MeshType>
859  std::vector<typename MeshType::active_cell_iterator>
860  compute_ghost_cell_halo_layer(const MeshType &mesh)
861  {
862  std::function<bool(const typename MeshType::active_cell_iterator &)>
863  predicate = IteratorFilters::LocallyOwnedCell();
864 
865  const std::vector<typename MeshType::active_cell_iterator>
866  active_halo_layer = compute_active_cell_halo_layer(mesh, predicate);
867 
868  // Check that we never return locally owned or artificial cells
869  // What is left should only be the ghost cells
870  Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) == false,
871  ExcMessage("Halo layer contains locally owned cells"));
872  Assert(contains_artificial_cells<MeshType>(active_halo_layer) == false,
873  ExcMessage("Halo layer contains artificial cells"));
874 
875  return active_halo_layer;
876  }
877 
878 
879 
880  template <class MeshType>
881  std::vector<typename MeshType::active_cell_iterator>
883  const MeshType &mesh,
884  const std::function<bool(const typename MeshType::active_cell_iterator &)>
885  & predicate,
886  const double layer_thickness)
887  {
888  std::vector<typename MeshType::active_cell_iterator>
889  subdomain_boundary_cells, active_cell_layer_within_distance;
890  std::vector<bool> vertices_outside_subdomain(
891  mesh.get_triangulation().n_vertices(), false);
892 
893  const unsigned int spacedim = MeshType::space_dimension;
894 
895  unsigned int n_non_predicate_cells = 0; // Number of non predicate cells
896 
897  // Find the layer of cells for which predicate is true and that
898  // are on the boundary with other cells. These are
899  // subdomain boundary cells.
900 
901  // Find the cells for which the predicate is false
902  // These are the cells which are around the predicate subdomain
903  for (const auto &cell : mesh.active_cell_iterators())
904  if (!predicate(cell)) // Negation of predicate --> Not Part of subdomain
905  {
906  for (const unsigned int v : cell->vertex_indices())
907  vertices_outside_subdomain[cell->vertex_index(v)] = true;
908  n_non_predicate_cells++;
909  }
910 
911  // If all the active cells conform to the predicate
912  // or if none of the active cells conform to the predicate
913  // there is no active cell layer around the predicate
914  // subdomain (within any distance)
915  if (n_non_predicate_cells == 0 ||
916  n_non_predicate_cells == mesh.get_triangulation().n_active_cells())
917  return std::vector<typename MeshType::active_cell_iterator>();
918 
919  // Find the cells that conform to the predicate
920  // but share a vertex with the cell not in the predicate subdomain
921  for (const auto &cell : mesh.active_cell_iterators())
922  if (predicate(cell)) // True predicate --> Potential boundary cell of the
923  // subdomain
924  for (const unsigned int v : cell->vertex_indices())
925  if (vertices_outside_subdomain[cell->vertex_index(v)] == true)
926  {
927  subdomain_boundary_cells.push_back(cell);
928  break; // No need to go through remaining vertices
929  }
930 
931  // To cheaply filter out some cells located far away from the predicate
932  // subdomain, get the bounding box of the predicate subdomain.
933  std::pair<Point<spacedim>, Point<spacedim>> bounding_box =
934  compute_bounding_box(mesh, predicate);
935 
936  // DOUBLE_EPSILON to compare really close double values
937  const double DOUBLE_EPSILON = 100. * std::numeric_limits<double>::epsilon();
938 
939  // Add layer_thickness to the bounding box
940  for (unsigned int d = 0; d < spacedim; ++d)
941  {
942  bounding_box.first[d] -= (layer_thickness + DOUBLE_EPSILON); // minp
943  bounding_box.second[d] += (layer_thickness + DOUBLE_EPSILON); // maxp
944  }
945 
946  std::vector<Point<spacedim>>
947  subdomain_boundary_cells_centers; // cache all the subdomain boundary
948  // cells centers here
949  std::vector<double>
950  subdomain_boundary_cells_radii; // cache all the subdomain boundary cells
951  // radii
952  subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
953  subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
954  // compute cell radius for each boundary cell of the predicate subdomain
955  for (typename std::vector<typename MeshType::active_cell_iterator>::
956  const_iterator subdomain_boundary_cell_iterator =
957  subdomain_boundary_cells.begin();
958  subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
959  ++subdomain_boundary_cell_iterator)
960  {
961  const std::pair<Point<spacedim>, double>
962  &subdomain_boundary_cell_enclosing_ball =
963  (*subdomain_boundary_cell_iterator)->enclosing_ball();
964 
965  subdomain_boundary_cells_centers.push_back(
966  subdomain_boundary_cell_enclosing_ball.first);
967  subdomain_boundary_cells_radii.push_back(
968  subdomain_boundary_cell_enclosing_ball.second);
969  }
970  AssertThrow(subdomain_boundary_cells_radii.size() ==
971  subdomain_boundary_cells_centers.size(),
972  ExcInternalError());
973 
974  // Find the cells that are within layer_thickness of predicate subdomain
975  // boundary distance but are inside the extended bounding box. Most cells
976  // might be outside the extended bounding box, so we could skip them. Those
977  // cells that are inside the extended bounding box but are not part of the
978  // predicate subdomain are possible candidates to be within the distance to
979  // the boundary cells of the predicate subdomain.
980  for (const auto &cell : mesh.active_cell_iterators())
981  {
982  // Ignore all the cells that are in the predicate subdomain
983  if (predicate(cell))
984  continue;
985 
986  const std::pair<Point<spacedim>, double> &cell_enclosing_ball =
987  cell->enclosing_ball();
988 
989  const Point<spacedim> cell_enclosing_ball_center =
990  cell_enclosing_ball.first;
991  const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
992 
993  bool cell_inside = true; // reset for each cell
994 
995  for (unsigned int d = 0; d < spacedim; ++d)
996  cell_inside &=
997  (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius >
998  bounding_box.first[d]) &&
999  (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius <
1000  bounding_box.second[d]);
1001  // cell_inside is true if its enclosing ball intersects the extended
1002  // bounding box
1003 
1004  // Ignore all the cells that are outside the extended bounding box
1005  if (cell_inside)
1006  for (unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1007  ++i)
1008  if (cell_enclosing_ball_center.distance_square(
1009  subdomain_boundary_cells_centers[i]) <
1010  Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1011  subdomain_boundary_cells_radii[i] +
1012  layer_thickness + DOUBLE_EPSILON))
1013  {
1014  active_cell_layer_within_distance.push_back(cell);
1015  break; // Exit the loop checking all the remaining subdomain
1016  // boundary cells
1017  }
1018  }
1019  return active_cell_layer_within_distance;
1020  }
1021 
1022 
1023 
1024  template <class MeshType>
1025  std::vector<typename MeshType::active_cell_iterator>
1027  const double layer_thickness)
1028  {
1029  IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate;
1030  std::function<bool(const typename MeshType::active_cell_iterator &)>
1031  predicate(locally_owned_cell_predicate);
1032 
1033  const std::vector<typename MeshType::active_cell_iterator>
1034  ghost_cell_layer_within_distance =
1036  predicate,
1037  layer_thickness);
1038 
1039  // Check that we never return locally owned or artificial cells
1040  // What is left should only be the ghost cells
1041  Assert(
1042  contains_locally_owned_cells<MeshType>(
1043  ghost_cell_layer_within_distance) == false,
1044  ExcMessage(
1045  "Ghost cells within layer_thickness contains locally owned cells."));
1046  Assert(
1047  contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1048  false,
1049  ExcMessage(
1050  "Ghost cells within layer_thickness contains artificial cells. "
1051  "The function compute_ghost_cell_layer_within_distance "
1052  "is probably called while using parallel::distributed::Triangulation. "
1053  "In such case please refer to the description of this function."));
1054 
1055  return ghost_cell_layer_within_distance;
1056  }
1057 
1058 
1059 
1060  template <class MeshType>
1061  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1063  const MeshType &mesh,
1064  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1065  &predicate)
1066  {
1067  std::vector<bool> locally_active_vertices_on_subdomain(
1068  mesh.get_triangulation().n_vertices(), false);
1069 
1070  const unsigned int spacedim = MeshType::space_dimension;
1071 
1072  // Two extreme points can define the bounding box
1073  // around the active cells that conform to the given predicate.
1075 
1076  // initialize minp and maxp with the first predicate cell center
1077  for (const auto &cell : mesh.active_cell_iterators())
1078  if (predicate(cell))
1079  {
1080  minp = cell->center();
1081  maxp = cell->center();
1082  break;
1083  }
1084 
1085  // Run through all the cells to check if it belongs to predicate domain,
1086  // if it belongs to the predicate domain, extend the bounding box.
1087  for (const auto &cell : mesh.active_cell_iterators())
1088  if (predicate(cell)) // True predicate --> Part of subdomain
1089  for (const unsigned int v : cell->vertex_indices())
1090  if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1091  false)
1092  {
1093  locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1094  true;
1095  for (unsigned int d = 0; d < spacedim; ++d)
1096  {
1097  minp[d] = std::min(minp[d], cell->vertex(v)[d]);
1098  maxp[d] = std::max(maxp[d], cell->vertex(v)[d]);
1099  }
1100  }
1101 
1102  return std::make_pair(minp, maxp);
1103  }
1104 
1105 
1106 
1107  template <typename MeshType>
1108  std::list<std::pair<typename MeshType::cell_iterator,
1109  typename MeshType::cell_iterator>>
1110  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
1111  {
1112  Assert(have_same_coarse_mesh(mesh_1, mesh_2),
1113  ExcMessage("The two meshes must be represent triangulations that "
1114  "have the same coarse meshes"));
1115  // We will allow the output to contain ghost cells when we have shared
1116  // Triangulations (i.e., so that each processor will get exactly the same
1117  // list of cell pairs), but not when we have two distributed
1118  // Triangulations (so that all active cells are partitioned by processor).
1119  // Non-parallel Triangulations have no ghost or artificial cells, so they
1120  // work the same way as shared Triangulations here.
1121  bool remove_ghost_cells = false;
1122 #ifdef DEAL_II_WITH_MPI
1123  {
1124  constexpr int dim = MeshType::dimension;
1125  constexpr int spacedim = MeshType::space_dimension;
1127  *>(&mesh_1.get_triangulation()) != nullptr ||
1129  *>(&mesh_2.get_triangulation()) != nullptr)
1130  {
1131  Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1132  ExcMessage("This function can only be used with meshes "
1133  "corresponding to distributed Triangulations when "
1134  "both Triangulations are equal."));
1135  remove_ghost_cells = true;
1136  }
1137  }
1138 #endif
1139 
1140  // the algorithm goes as follows: first, we fill a list with pairs of
1141  // iterators common to the two meshes on the coarsest level. then we
1142  // traverse the list; each time, we find a pair of iterators for which
1143  // both correspond to non-active cells, we delete this item and push the
1144  // pairs of iterators to their children to the back. if these again both
1145  // correspond to non-active cells, we will get to the later on for further
1146  // consideration
1147  using CellList = std::list<std::pair<typename MeshType::cell_iterator,
1148  typename MeshType::cell_iterator>>;
1149  CellList cell_list;
1150 
1151  // first push the coarse level cells
1152  typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1153  cell_2 = mesh_2.begin(0);
1154  for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1155  cell_list.emplace_back(cell_1, cell_2);
1156 
1157  // then traverse list as described above
1158  typename CellList::iterator cell_pair = cell_list.begin();
1159  while (cell_pair != cell_list.end())
1160  {
1161  // if both cells in this pair have children, then erase this element
1162  // and push their children instead
1163  if (cell_pair->first->has_children() &&
1164  cell_pair->second->has_children())
1165  {
1166  Assert(cell_pair->first->refinement_case() ==
1167  cell_pair->second->refinement_case(),
1168  ExcNotImplemented());
1169  for (unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1170  cell_list.emplace_back(cell_pair->first->child(c),
1171  cell_pair->second->child(c));
1172 
1173  // erasing an iterator keeps other iterators valid, so already
1174  // advance the present iterator by one and then delete the element
1175  // we've visited before
1176  const auto previous_cell_pair = cell_pair;
1177  ++cell_pair;
1178  cell_list.erase(previous_cell_pair);
1179  }
1180  else
1181  {
1182  // at least one cell is active
1183  if (remove_ghost_cells &&
1184  ((cell_pair->first->is_active() &&
1185  !cell_pair->first->is_locally_owned()) ||
1186  (cell_pair->second->is_active() &&
1187  !cell_pair->second->is_locally_owned())))
1188  {
1189  // we only exclude ghost cells for distributed Triangulations
1190  const auto previous_cell_pair = cell_pair;
1191  ++cell_pair;
1192  cell_list.erase(previous_cell_pair);
1193  }
1194  else
1195  ++cell_pair;
1196  }
1197  }
1198 
1199  // just to make sure everything is ok, validate that all pairs have at
1200  // least one active iterator or have different refinement_cases
1201  for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1202  ++cell_pair)
1203  Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1204  (cell_pair->first->refinement_case() !=
1205  cell_pair->second->refinement_case()),
1206  ExcInternalError());
1207 
1208  return cell_list;
1209  }
1210 
1211 
1212 
1213  template <int dim, int spacedim>
1214  bool
1216  const Triangulation<dim, spacedim> &mesh_2)
1217  {
1218  // make sure the two meshes have
1219  // the same number of coarse cells
1220  if (mesh_1.n_cells(0) != mesh_2.n_cells(0))
1221  return false;
1222 
1223  // if so, also make sure they have
1224  // the same vertices on the cells
1225  // of the coarse mesh
1227  mesh_1.begin(0),
1228  cell_2 =
1229  mesh_2.begin(0),
1230  endc = mesh_1.end(0);
1231  for (; cell_1 != endc; ++cell_1, ++cell_2)
1232  {
1233  if (cell_1->n_vertices() != cell_2->n_vertices())
1234  return false;
1235  for (const unsigned int v : cell_1->vertex_indices())
1236  if (cell_1->vertex(v) != cell_2->vertex(v))
1237  return false;
1238  }
1239 
1240  // if we've gotten through all
1241  // this, then the meshes really
1242  // seem to have a common coarse
1243  // mesh
1244  return true;
1245  }
1246 
1247 
1248 
1249  template <typename MeshType>
1250  bool
1251  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2)
1252  {
1253  return have_same_coarse_mesh(mesh_1.get_triangulation(),
1254  mesh_2.get_triangulation());
1255  }
1256 
1257 
1258 
1259  template <int dim, int spacedim>
1260  std::pair<typename ::DoFHandler<dim, spacedim>::active_cell_iterator,
1261  Point<dim>>
1263  const hp::MappingCollection<dim, spacedim> &mapping,
1264  const DoFHandler<dim, spacedim> & mesh,
1265  const Point<spacedim> & p,
1266  const double tolerance)
1267  {
1268  Assert((mapping.size() == 1) ||
1269  (mapping.size() == mesh.get_fe_collection().size()),
1270  ExcMessage("Mapping collection needs to have either size 1 "
1271  "or size equal to the number of elements in "
1272  "the FECollection."));
1273 
1274  using cell_iterator =
1275  typename ::DoFHandler<dim, spacedim>::active_cell_iterator;
1276 
1277  std::pair<cell_iterator, Point<dim>> best_cell;
1278  // If we have only one element in the MappingCollection,
1279  // we use find_active_cell_around_point using only one
1280  // mapping.
1281  if (mapping.size() == 1)
1282  {
1283  const std::vector<bool> marked_vertices = {};
1284  best_cell = find_active_cell_around_point(
1285  mapping[0], mesh, p, marked_vertices, tolerance);
1286  }
1287  else
1288  {
1289  // The best distance is set to the
1290  // maximum allowable distance from
1291  // the unit cell
1292  double best_distance = tolerance;
1293  int best_level = -1;
1294 
1295 
1296  // Find closest vertex and determine
1297  // all adjacent cells
1298  unsigned int vertex = find_closest_vertex(mesh, p);
1299 
1300  std::vector<cell_iterator> adjacent_cells_tmp =
1301  find_cells_adjacent_to_vertex(mesh, vertex);
1302 
1303  // Make sure that we have found
1304  // at least one cell adjacent to vertex.
1305  Assert(adjacent_cells_tmp.size() > 0, ExcInternalError());
1306 
1307  // Copy all the cells into a std::set
1308  std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(),
1309  adjacent_cells_tmp.end());
1310  std::set<cell_iterator> searched_cells;
1311 
1312  // Determine the maximal number of cells
1313  // in the grid.
1314  // As long as we have not found
1315  // the cell and have not searched
1316  // every cell in the triangulation,
1317  // we keep on looking.
1318  const unsigned int n_cells = mesh.get_triangulation().n_cells();
1319  bool found = false;
1320  unsigned int cells_searched = 0;
1321  while (!found && cells_searched < n_cells)
1322  {
1323  typename std::set<cell_iterator>::const_iterator
1324  cell = adjacent_cells.begin(),
1325  endc = adjacent_cells.end();
1326  for (; cell != endc; ++cell)
1327  {
1328  try
1329  {
1330  const Point<dim> p_cell =
1331  mapping[(*cell)->active_fe_index()]
1332  .transform_real_to_unit_cell(*cell, p);
1333 
1334 
1335  // calculate the infinity norm of
1336  // the distance vector to the unit cell.
1337  const double dist =
1339 
1340  // We compare if the point is inside the
1341  // unit cell (or at least not too far
1342  // outside). If it is, it is also checked
1343  // that the cell has a more refined state
1344  if (dist < best_distance || (dist == best_distance &&
1345  (*cell)->level() > best_level))
1346  {
1347  found = true;
1348  best_distance = dist;
1349  best_level = (*cell)->level();
1350  best_cell = std::make_pair(*cell, p_cell);
1351  }
1352  }
1353  catch (
1354  typename MappingQGeneric<dim,
1355  spacedim>::ExcTransformationFailed &)
1356  {
1357  // ok, the transformation
1358  // failed presumably
1359  // because the point we
1360  // are looking for lies
1361  // outside the current
1362  // cell. this means that
1363  // the current cell can't
1364  // be the cell around the
1365  // point, so just ignore
1366  // this cell and move on
1367  // to the next
1368  }
1369  }
1370  // update the number of cells searched
1371  cells_searched += adjacent_cells.size();
1372  // if we have not found the cell in
1373  // question and have not yet searched every
1374  // cell, we expand our search to
1375  // all the not already searched neighbors of
1376  // the cells in adjacent_cells.
1377  if (!found && cells_searched < n_cells)
1378  {
1379  find_active_cell_around_point_internal<dim,
1380  ::DoFHandler,
1381  spacedim>(
1382  mesh, searched_cells, adjacent_cells);
1383  }
1384  }
1385  }
1386 
1387  AssertThrow(best_cell.first.state() == IteratorState::valid,
1388  ExcPointNotFound<spacedim>(p));
1389 
1390  return best_cell;
1391  }
1392 
1393 
1394  template <class MeshType>
1395  std::vector<typename MeshType::active_cell_iterator>
1396  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
1397  {
1398  Assert(cell->is_locally_owned(),
1399  ExcMessage("This function only makes sense if the cell for "
1400  "which you are asking for a patch, is locally "
1401  "owned."));
1402 
1403  std::vector<typename MeshType::active_cell_iterator> patch;
1404  patch.push_back(cell);
1405  for (const unsigned int face_number : cell->face_indices())
1406  if (cell->face(face_number)->at_boundary() == false)
1407  {
1408  if (cell->neighbor(face_number)->has_children() == false)
1409  patch.push_back(cell->neighbor(face_number));
1410  else
1411  // the neighbor is refined. in 2d/3d, we can simply ask for the
1412  // children of the neighbor because they can not be further refined
1413  // and, consequently, the children is active
1414  if (MeshType::dimension > 1)
1415  {
1416  for (unsigned int subface = 0;
1417  subface < cell->face(face_number)->n_children();
1418  ++subface)
1419  patch.push_back(
1420  cell->neighbor_child_on_subface(face_number, subface));
1421  }
1422  else
1423  {
1424  // in 1d, we need to work a bit harder: iterate until we find
1425  // the child by going from cell to child to child etc
1426  typename MeshType::cell_iterator neighbor =
1427  cell->neighbor(face_number);
1428  while (neighbor->has_children())
1429  neighbor = neighbor->child(1 - face_number);
1430 
1431  Assert(neighbor->neighbor(1 - face_number) == cell,
1432  ExcInternalError());
1433  patch.push_back(neighbor);
1434  }
1435  }
1436  return patch;
1437  }
1438 
1439 
1440 
1441  template <class Container>
1442  std::vector<typename Container::cell_iterator>
1444  const std::vector<typename Container::active_cell_iterator> &patch)
1445  {
1446  Assert(patch.size() > 0,
1447  ExcMessage(
1448  "Vector containing patch cells should not be an empty vector!"));
1449  // In order to extract the set of cells with the coarsest common level from
1450  // the give vector of cells: First it finds the number associated with the
1451  // minimum level of refinement, namely "min_level"
1452  int min_level = patch[0]->level();
1453 
1454  for (unsigned int i = 0; i < patch.size(); ++i)
1455  min_level = std::min(min_level, patch[i]->level());
1456  std::set<typename Container::cell_iterator> uniform_cells;
1457  typename std::vector<
1458  typename Container::active_cell_iterator>::const_iterator patch_cell;
1459  // it loops through all cells of the input vector
1460  for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1461  {
1462  // If the refinement level of each cell i the loop be equal to the
1463  // min_level, so that that cell inserted into the set of uniform_cells,
1464  // as the set of cells with the coarsest common refinement level
1465  if ((*patch_cell)->level() == min_level)
1466  uniform_cells.insert(*patch_cell);
1467  else
1468  // If not, it asks for the parent of the cell, until it finds the
1469  // parent cell with the refinement level equal to the min_level and
1470  // inserts that parent cell into the the set of uniform_cells, as the
1471  // set of cells with the coarsest common refinement level.
1472  {
1473  typename Container::cell_iterator parent = *patch_cell;
1474 
1475  while (parent->level() > min_level)
1476  parent = parent->parent();
1477  uniform_cells.insert(parent);
1478  }
1479  }
1480 
1481  return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1482  uniform_cells.end());
1483  }
1484 
1485 
1486 
1487  template <class Container>
1488  void
1490  const std::vector<typename Container::active_cell_iterator> &patch,
1492  &local_triangulation,
1493  std::map<
1494  typename Triangulation<Container::dimension,
1495  Container::space_dimension>::active_cell_iterator,
1496  typename Container::active_cell_iterator> &patch_to_global_tria_map)
1497 
1498  {
1499  const std::vector<typename Container::cell_iterator> uniform_cells =
1500  get_cells_at_coarsest_common_level<Container>(patch);
1501  // First it creates triangulation from the vector of "uniform_cells"
1502  local_triangulation.clear();
1503  std::vector<Point<Container::space_dimension>> vertices;
1504  const unsigned int n_uniform_cells = uniform_cells.size();
1505  std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1506  unsigned int k = 0; // for enumerating cells
1507  unsigned int i = 0; // for enumerating vertices
1508  typename std::vector<typename Container::cell_iterator>::const_iterator
1509  uniform_cell;
1510  for (uniform_cell = uniform_cells.begin();
1511  uniform_cell != uniform_cells.end();
1512  ++uniform_cell)
1513  {
1514  for (const unsigned int v : (*uniform_cell)->vertex_indices())
1515  {
1517  (*uniform_cell)->vertex(v);
1518  bool repeat_vertex = false;
1519 
1520  for (unsigned int m = 0; m < i; ++m)
1521  {
1522  if (position == vertices[m])
1523  {
1524  repeat_vertex = true;
1525  cells[k].vertices[v] = m;
1526  break;
1527  }
1528  }
1529  if (repeat_vertex == false)
1530  {
1531  vertices.push_back(position);
1532  cells[k].vertices[v] = i;
1533  i = i + 1;
1534  }
1535 
1536  } // for vertices_per_cell
1537  k = k + 1;
1538  }
1539  local_triangulation.create_triangulation(vertices, cells, SubCellData());
1540  Assert(local_triangulation.n_active_cells() == uniform_cells.size(),
1541  ExcInternalError());
1542  local_triangulation.clear_user_flags();
1543  unsigned int index = 0;
1544  // Create a map between cells of class DoFHandler into class Triangulation
1545  std::map<typename Triangulation<Container::dimension,
1546  Container::space_dimension>::cell_iterator,
1547  typename Container::cell_iterator>
1548  patch_to_global_tria_map_tmp;
1549  for (typename Triangulation<Container::dimension,
1550  Container::space_dimension>::cell_iterator
1551  coarse_cell = local_triangulation.begin();
1552  coarse_cell != local_triangulation.end();
1553  ++coarse_cell, ++index)
1554  {
1555  patch_to_global_tria_map_tmp.insert(
1556  std::make_pair(coarse_cell, uniform_cells[index]));
1557  // To ensure that the cells with the same coordinates (here, we compare
1558  // their centers) are mapped into each other.
1559 
1560  Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1561  1e-15 * coarse_cell->diameter(),
1562  ExcInternalError());
1563  }
1564  bool refinement_necessary;
1565  // In this loop we start to do refinement on the above coarse triangulation
1566  // to reach to the same level of refinement as the patch cells are really on
1567  do
1568  {
1569  refinement_necessary = false;
1570  for (const auto &active_tria_cell :
1571  local_triangulation.active_cell_iterators())
1572  {
1573  if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1574  {
1575  active_tria_cell->set_refine_flag();
1576  refinement_necessary = true;
1577  }
1578  else
1579  for (unsigned int i = 0; i < patch.size(); ++i)
1580  {
1581  // Even though vertices may not be exactly the same, the
1582  // appropriate cells will match since == for TriAccessors
1583  // checks only cell level and index.
1584  if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1585  patch[i])
1586  {
1587  // adjust the cell vertices of the local_triangulation to
1588  // match cell vertices of the global triangulation
1589  for (const unsigned int v :
1590  active_tria_cell->vertex_indices())
1591  active_tria_cell->vertex(v) = patch[i]->vertex(v);
1592 
1593  Assert(active_tria_cell->center().distance(
1594  patch_to_global_tria_map_tmp[active_tria_cell]
1595  ->center()) <=
1596  1e-15 * active_tria_cell->diameter(),
1597  ExcInternalError());
1598 
1599  active_tria_cell->set_user_flag();
1600  break;
1601  }
1602  }
1603  }
1604 
1605  if (refinement_necessary)
1606  {
1607  local_triangulation.execute_coarsening_and_refinement();
1608 
1609  for (typename Triangulation<
1610  Container::dimension,
1611  Container::space_dimension>::cell_iterator cell =
1612  local_triangulation.begin();
1613  cell != local_triangulation.end();
1614  ++cell)
1615  {
1616  if (patch_to_global_tria_map_tmp.find(cell) !=
1617  patch_to_global_tria_map_tmp.end())
1618  {
1619  if (cell->has_children())
1620  {
1621  // Note: Since the cell got children, then it should not
1622  // be in the map anymore children may be added into the
1623  // map, instead
1624 
1625  // these children may not yet be in the map
1626  for (unsigned int c = 0; c < cell->n_children(); ++c)
1627  {
1628  if (patch_to_global_tria_map_tmp.find(cell->child(
1629  c)) == patch_to_global_tria_map_tmp.end())
1630  {
1631  patch_to_global_tria_map_tmp.insert(
1632  std::make_pair(
1633  cell->child(c),
1634  patch_to_global_tria_map_tmp[cell]->child(
1635  c)));
1636 
1637  // One might be tempted to assert that the cell
1638  // being added here has the same center as the
1639  // equivalent cell in the global triangulation,
1640  // but it may not be the case. For
1641  // triangulations that have been perturbed or
1642  // smoothed, the cell indices and levels may be
1643  // the same, but the vertex locations may not.
1644  // We adjust the vertices of the cells that have
1645  // no children (ie the active cells) to be
1646  // consistent with the global triangulation
1647  // later on and add assertions at that time
1648  // to guarantee the cells in the
1649  // local_triangulation are physically at the
1650  // same locations of the cells in the patch of
1651  // the global triangulation.
1652  }
1653  }
1654  // The parent cell whose children were added
1655  // into the map should be deleted from the map
1656  patch_to_global_tria_map_tmp.erase(cell);
1657  }
1658  }
1659  }
1660  }
1661  }
1662  while (refinement_necessary);
1663 
1664 
1665  // Last assertion check to make sure we have the right cells and centers
1666  // in the map, and hence the correct vertices of the triangulation
1667  for (typename Triangulation<Container::dimension,
1668  Container::space_dimension>::cell_iterator
1669  cell = local_triangulation.begin();
1670  cell != local_triangulation.end();
1671  ++cell)
1672  {
1673  if (cell->user_flag_set())
1674  {
1675  Assert(patch_to_global_tria_map_tmp.find(cell) !=
1676  patch_to_global_tria_map_tmp.end(),
1677  ExcInternalError());
1678 
1679  Assert(cell->center().distance(
1680  patch_to_global_tria_map_tmp[cell]->center()) <=
1681  1e-15 * cell->diameter(),
1682  ExcInternalError());
1683  }
1684  }
1685 
1686 
1687  typename std::map<
1688  typename Triangulation<Container::dimension,
1689  Container::space_dimension>::cell_iterator,
1690  typename Container::cell_iterator>::iterator
1691  map_tmp_it = patch_to_global_tria_map_tmp.begin(),
1692  map_tmp_end = patch_to_global_tria_map_tmp.end();
1693  // Now we just need to take the temporary map of pairs of type cell_iterator
1694  // "patch_to_global_tria_map_tmp" making pair of active_cell_iterators so
1695  // that filling out the final map "patch_to_global_tria_map"
1696  for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1697  patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1698  }
1699 
1700 
1701 
1702  template <int dim, int spacedim>
1703  std::map<
1705  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1707  {
1708  // This is the map from global_dof_index to
1709  // a set of cells on patch. We first map into
1710  // a set because it is very likely that we
1711  // will attempt to add a cell more than once
1712  // to a particular patch and we want to preserve
1713  // uniqueness of cell iterators. std::set does this
1714  // automatically for us. Later after it is all
1715  // constructed, we will copy to a map of vectors
1716  // since that is the preferred output for other
1717  // functions.
1718  std::map<types::global_dof_index,
1719  std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1720  dof_to_set_of_cells_map;
1721 
1722  std::vector<types::global_dof_index> local_dof_indices;
1723  std::vector<types::global_dof_index> local_face_dof_indices;
1724  std::vector<types::global_dof_index> local_line_dof_indices;
1725 
1726  // a place to save the dof_handler user flags and restore them later
1727  // to maintain const of dof_handler.
1728  std::vector<bool> user_flags;
1729 
1730 
1731  // in 3d, we need pointers from active lines to the
1732  // active parent lines, so we construct it as needed.
1733  std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1735  lines_to_parent_lines_map;
1736  if (dim == 3)
1737  {
1738  // save user flags as they will be modified and then later restored
1739  dof_handler.get_triangulation().save_user_flags(user_flags);
1740  const_cast<::Triangulation<dim, spacedim> &>(
1741  dof_handler.get_triangulation())
1742  .clear_user_flags();
1743 
1744 
1746  cell = dof_handler.begin_active(),
1747  endc = dof_handler.end();
1748  for (; cell != endc; ++cell)
1749  {
1750  // We only want lines that are locally_relevant
1751  // although it doesn't hurt to have lines that
1752  // are children of ghost cells since there are
1753  // few and we don't have to use them.
1754  if (cell->is_artificial() == false)
1755  {
1756  for (unsigned int l = 0; l < cell->n_lines(); ++l)
1757  if (cell->line(l)->has_children())
1758  for (unsigned int c = 0; c < cell->line(l)->n_children();
1759  ++c)
1760  {
1761  lines_to_parent_lines_map[cell->line(l)->child(c)] =
1762  cell->line(l);
1763  // set flags to know that child
1764  // line has an active parent.
1765  cell->line(l)->child(c)->set_user_flag();
1766  }
1767  }
1768  }
1769  }
1770 
1771 
1772  // We loop through all cells and add cell to the
1773  // map for the dofs that it immediately touches
1774  // and then account for all the other dofs of
1775  // which it is a part, mainly the ones that must
1776  // be added on account of adaptivity hanging node
1777  // constraints.
1779  cell = dof_handler.begin_active(),
1780  endc = dof_handler.end();
1781  for (; cell != endc; ++cell)
1782  {
1783  // Need to loop through all cells that could
1784  // be in the patch of dofs on locally_owned
1785  // cells including ghost cells
1786  if (cell->is_artificial() == false)
1787  {
1788  const unsigned int n_dofs_per_cell =
1789  cell->get_fe().n_dofs_per_cell();
1790  local_dof_indices.resize(n_dofs_per_cell);
1791 
1792  // Take care of adding cell pointer to each
1793  // dofs that exists on cell.
1794  cell->get_dof_indices(local_dof_indices);
1795  for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
1796  dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell);
1797 
1798  // In the case of the adjacent cell (over
1799  // faces or edges) being more refined, we
1800  // want to add all of the children to the
1801  // patch since the support function at that
1802  // dof could be non-zero along that entire
1803  // face (or line).
1804 
1805  // Take care of dofs on neighbor faces
1806  for (const unsigned int f : cell->face_indices())
1807  {
1808  if (cell->face(f)->has_children())
1809  {
1810  for (unsigned int c = 0; c < cell->face(f)->n_children();
1811  ++c)
1812  {
1813  // Add cell to dofs of all subfaces
1814  //
1815  // *-------------------*----------*---------*
1816  // | | add cell | |
1817  // | |<- to dofs| |
1818  // | |of subface| |
1819  // | cell *----------*---------*
1820  // | | add cell | |
1821  // | |<- to dofs| |
1822  // | |of subface| |
1823  // *-------------------*----------*---------*
1824  //
1825  Assert(cell->face(f)->child(c)->has_children() == false,
1826  ExcInternalError());
1827 
1828  const unsigned int n_dofs_per_face =
1829  cell->get_fe().n_dofs_per_face(f, c);
1830  local_face_dof_indices.resize(n_dofs_per_face);
1831 
1832  cell->face(f)->child(c)->get_dof_indices(
1833  local_face_dof_indices);
1834  for (unsigned int i = 0; i < n_dofs_per_face; ++i)
1835  dof_to_set_of_cells_map[local_face_dof_indices[i]]
1836  .insert(cell);
1837  }
1838  }
1839  else if ((cell->face(f)->at_boundary() == false) &&
1840  (cell->neighbor_is_coarser(f)))
1841  {
1842  // Add cell to dofs of parent face and all
1843  // child faces of parent face
1844  //
1845  // *-------------------*----------*---------*
1846  // | | | |
1847  // | | cell | |
1848  // | add cell | | |
1849  // | to dofs -> *----------*---------*
1850  // | of parent | add cell | |
1851  // | face |<- to dofs| |
1852  // | |of subface| |
1853  // *-------------------*----------*---------*
1854  //
1855 
1856  // Add cell to all dofs of parent face
1857  std::pair<unsigned int, unsigned int>
1858  neighbor_face_no_subface_no =
1859  cell->neighbor_of_coarser_neighbor(f);
1860  unsigned int face_no = neighbor_face_no_subface_no.first;
1861  unsigned int subface = neighbor_face_no_subface_no.second;
1862 
1863  const unsigned int n_dofs_per_face =
1864  cell->get_fe().n_dofs_per_face(face_no);
1865  local_face_dof_indices.resize(n_dofs_per_face);
1866 
1867  cell->neighbor(f)->face(face_no)->get_dof_indices(
1868  local_face_dof_indices);
1869  for (unsigned int i = 0; i < n_dofs_per_face; ++i)
1870  dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(
1871  cell);
1872 
1873  // Add cell to all dofs of children of
1874  // parent face
1875  for (unsigned int c = 0;
1876  c < cell->neighbor(f)->face(face_no)->n_children();
1877  ++c)
1878  {
1879  if (c != subface) // don't repeat work on dofs of
1880  // original cell
1881  {
1882  const unsigned int n_dofs_per_face =
1883  cell->get_fe().n_dofs_per_face(face_no, c);
1884  local_face_dof_indices.resize(n_dofs_per_face);
1885 
1886  Assert(cell->neighbor(f)
1887  ->face(face_no)
1888  ->child(c)
1889  ->has_children() == false,
1890  ExcInternalError());
1891  cell->neighbor(f)
1892  ->face(face_no)
1893  ->child(c)
1894  ->get_dof_indices(local_face_dof_indices);
1895  for (unsigned int i = 0; i < n_dofs_per_face; ++i)
1896  dof_to_set_of_cells_map[local_face_dof_indices[i]]
1897  .insert(cell);
1898  }
1899  }
1900  }
1901  }
1902 
1903 
1904  // If 3d, take care of dofs on lines in the
1905  // same pattern as faces above. That is, if
1906  // a cell's line has children, distribute
1907  // cell to dofs of children of line, and
1908  // if cell's line has an active parent, then
1909  // distribute cell to dofs on parent line
1910  // and dofs on all children of parent line.
1911  if (dim == 3)
1912  {
1913  for (unsigned int l = 0; l < cell->n_lines(); ++l)
1914  {
1915  if (cell->line(l)->has_children())
1916  {
1917  for (unsigned int c = 0;
1918  c < cell->line(l)->n_children();
1919  ++c)
1920  {
1921  Assert(cell->line(l)->child(c)->has_children() ==
1922  false,
1923  ExcInternalError());
1924 
1925  // dofs_per_line returns number of dofs
1926  // on line not including the vertices of the line.
1927  const unsigned int n_dofs_per_line =
1928  2 * cell->get_fe().n_dofs_per_vertex() +
1929  cell->get_fe().n_dofs_per_line();
1930  local_line_dof_indices.resize(n_dofs_per_line);
1931 
1932  cell->line(l)->child(c)->get_dof_indices(
1933  local_line_dof_indices);
1934  for (unsigned int i = 0; i < n_dofs_per_line; ++i)
1935  dof_to_set_of_cells_map[local_line_dof_indices[i]]
1936  .insert(cell);
1937  }
1938  }
1939  // user flag was set above to denote that
1940  // an active parent line exists so add
1941  // cell to dofs of parent and all it's
1942  // children
1943  else if (cell->line(l)->user_flag_set() == true)
1944  {
1946  parent_line =
1947  lines_to_parent_lines_map[cell->line(l)];
1948  Assert(parent_line->has_children(), ExcInternalError());
1949 
1950  // dofs_per_line returns number of dofs
1951  // on line not including the vertices of the line.
1952  const unsigned int n_dofs_per_line =
1953  2 * cell->get_fe().n_dofs_per_vertex() +
1954  cell->get_fe().n_dofs_per_line();
1955  local_line_dof_indices.resize(n_dofs_per_line);
1956 
1957  parent_line->get_dof_indices(local_line_dof_indices);
1958  for (unsigned int i = 0; i < n_dofs_per_line; ++i)
1959  dof_to_set_of_cells_map[local_line_dof_indices[i]]
1960  .insert(cell);
1961 
1962  for (unsigned int c = 0; c < parent_line->n_children();
1963  ++c)
1964  {
1965  Assert(parent_line->child(c)->has_children() ==
1966  false,
1967  ExcInternalError());
1968 
1969  const unsigned int n_dofs_per_line =
1970  2 * cell->get_fe().n_dofs_per_vertex() +
1971  cell->get_fe().n_dofs_per_line();
1972  local_line_dof_indices.resize(n_dofs_per_line);
1973 
1974  parent_line->child(c)->get_dof_indices(
1975  local_line_dof_indices);
1976  for (unsigned int i = 0; i < n_dofs_per_line; ++i)
1977  dof_to_set_of_cells_map[local_line_dof_indices[i]]
1978  .insert(cell);
1979  }
1980  }
1981  } // for lines l
1982  } // if dim == 3
1983  } // if cell->is_locally_owned()
1984  } // for cells
1985 
1986 
1987  if (dim == 3)
1988  {
1989  // finally, restore user flags that were changed above
1990  // to when we constructed the pointers to parent of lines
1991  // Since dof_handler is const, we must leave it unchanged.
1992  const_cast<::Triangulation<dim, spacedim> &>(
1993  dof_handler.get_triangulation())
1994  .load_user_flags(user_flags);
1995  }
1996 
1997  // Finally, we copy map of sets to
1998  // map of vectors using the std::vector::assign() function
1999  std::map<
2001  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2002  dof_to_cell_patches;
2003 
2004  typename std::map<
2006  std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
2007  iterator it = dof_to_set_of_cells_map.begin(),
2008  it_end = dof_to_set_of_cells_map.end();
2009  for (; it != it_end; ++it)
2010  dof_to_cell_patches[it->first].assign(it->second.begin(),
2011  it->second.end());
2012 
2013  return dof_to_cell_patches;
2014  }
2015 
2016  /*
2017  * Internally used in collect_periodic_faces
2018  */
2019  template <typename CellIterator>
2020  void
2022  std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2023  std::set<std::pair<typename identity<CellIterator>::type, unsigned int>>
2024  & pairs2,
2025  const int direction,
2026  std::vector<PeriodicFacePair<CellIterator>> &matched_pairs,
2027  const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2028  & offset,
2029  const FullMatrix<double> &matrix)
2030  {
2031  static const int space_dim = CellIterator::AccessorType::space_dimension;
2032  (void)space_dim;
2033  AssertIndexRange(direction, space_dim);
2034 
2035 #ifdef DEBUG
2036  {
2037  constexpr int dim = CellIterator::AccessorType::dimension;
2038  constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2039  // For parallel::fullydistributed::Triangulation there might be unmatched
2040  // faces on periodic boundaries on the coarse grid, which results that
2041  // this assert is not fulfilled (not a bug!). See also the discussion in
2042  // the method collect_periodic_faces.
2043  if (!(((pairs1.size() > 0) &&
2044  (dynamic_cast<const parallel::fullydistributed::
2045  Triangulation<dim, spacedim> *>(
2046  &pairs1.begin()->first->get_triangulation()) != nullptr)) ||
2047  ((pairs2.size() > 0) &&
2048  (dynamic_cast<
2050  *>(&pairs2.begin()->first->get_triangulation()) != nullptr))))
2051  Assert(pairs1.size() == pairs2.size(),
2052  ExcMessage("Unmatched faces on periodic boundaries"));
2053  }
2054 #endif
2055 
2056  unsigned int n_matches = 0;
2057 
2058  // Match with a complexity of O(n^2). This could be improved...
2059  std::bitset<3> orientation;
2060  using PairIterator =
2061  typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2062  for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2063  {
2064  for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2065  {
2066  const CellIterator cell1 = it1->first;
2067  const CellIterator cell2 = it2->first;
2068  const unsigned int face_idx1 = it1->second;
2069  const unsigned int face_idx2 = it2->second;
2070  if (GridTools::orthogonal_equality(orientation,
2071  cell1->face(face_idx1),
2072  cell2->face(face_idx2),
2073  direction,
2074  offset,
2075  matrix))
2076  {
2077  // We have a match, so insert the matching pairs and
2078  // remove the matched cell in pairs2 to speed up the
2079  // matching:
2080  const PeriodicFacePair<CellIterator> matched_face = {
2081  {cell1, cell2}, {face_idx1, face_idx2}, orientation, matrix};
2082  matched_pairs.push_back(matched_face);
2083  pairs2.erase(it2);
2084  ++n_matches;
2085  break;
2086  }
2087  }
2088  }
2089 
2090  // Assure that all faces are matched if not
2091  // parallel::fullydistributed::Triangulation is used. This is related to the
2092  // fact that the faces might not be successfully matched on the coarse grid
2093  // (not a bug!). See also the comment above and in the
2094  // method collect_periodic_faces.
2095  {
2096  constexpr int dim = CellIterator::AccessorType::dimension;
2097  constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2098  if (!(((pairs1.size() > 0) &&
2099  (dynamic_cast<const parallel::fullydistributed::
2100  Triangulation<dim, spacedim> *>(
2101  &pairs1.begin()->first->get_triangulation()) != nullptr)) ||
2102  ((pairs2.size() > 0) &&
2103  (dynamic_cast<
2105  *>(&pairs2.begin()->first->get_triangulation()) != nullptr))))
2106  AssertThrow(n_matches == pairs1.size() && pairs2.size() == 0,
2107  ExcMessage("Unmatched faces on periodic boundaries"));
2108  }
2109  }
2110 
2111 
2112 
2113  template <typename MeshType>
2114  void
2116  const MeshType & mesh,
2117  const types::boundary_id b_id,
2118  const int direction,
2120  & matched_pairs,
2122  const FullMatrix<double> & matrix)
2123  {
2124  static const int dim = MeshType::dimension;
2125  static const int space_dim = MeshType::space_dimension;
2126  (void)dim;
2127  (void)space_dim;
2128  AssertIndexRange(direction, space_dim);
2129 
2130  Assert(dim == space_dim, ExcNotImplemented());
2131 
2132  // Loop over all cells on the highest level and collect all boundary
2133  // faces 2*direction and 2*direction*1:
2134 
2135  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2136  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2137 
2138  for (typename MeshType::cell_iterator cell = mesh.begin(0);
2139  cell != mesh.end(0);
2140  ++cell)
2141  {
2142  const typename MeshType::face_iterator face_1 =
2143  cell->face(2 * direction);
2144  const typename MeshType::face_iterator face_2 =
2145  cell->face(2 * direction + 1);
2146 
2147  if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2148  {
2149  const std::pair<typename MeshType::cell_iterator, unsigned int>
2150  pair1 = std::make_pair(cell, 2 * direction);
2151  pairs1.insert(pair1);
2152  }
2153 
2154  if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2155  {
2156  const std::pair<typename MeshType::cell_iterator, unsigned int>
2157  pair2 = std::make_pair(cell, 2 * direction + 1);
2158  pairs2.insert(pair2);
2159  }
2160  }
2161 
2162  Assert(pairs1.size() == pairs2.size(),
2163  ExcMessage("Unmatched faces on periodic boundaries"));
2164 
2165  Assert(pairs1.size() > 0,
2166  ExcMessage("No new periodic face pairs have been found. "
2167  "Are you sure that you've selected the correct boundary "
2168  "id's and that the coarsest level mesh is colorized?"));
2169 
2170 #ifdef DEBUG
2171  const unsigned int size_old = matched_pairs.size();
2172 #endif
2173 
2174  // and call match_periodic_face_pairs that does the actual matching:
2176  pairs1, pairs2, direction, matched_pairs, offset, matrix);
2177 
2178 #ifdef DEBUG
2179  // check for standard orientation
2180  const unsigned int size_new = matched_pairs.size();
2181  for (unsigned int i = size_old; i < size_new; ++i)
2182  {
2183  Assert(matched_pairs[i].orientation == 1,
2184  ExcMessage(
2185  "Found a face match with non standard orientation. "
2186  "This function is only suitable for meshes with cells "
2187  "in default orientation"));
2188  }
2189 #endif
2190  }
2191 
2192 
2193 
2194  template <typename MeshType>
2195  void
2197  const MeshType & mesh,
2198  const types::boundary_id b_id1,
2199  const types::boundary_id b_id2,
2200  const int direction,
2202  & matched_pairs,
2204  const FullMatrix<double> & matrix)
2205  {
2206  static const int dim = MeshType::dimension;
2207  static const int space_dim = MeshType::space_dimension;
2208  (void)dim;
2209  (void)space_dim;
2210  AssertIndexRange(direction, space_dim);
2211 
2212  // Loop over all cells on the highest level and collect all boundary
2213  // faces belonging to b_id1 and b_id2:
2214 
2215  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2216  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2217 
2218  for (typename MeshType::cell_iterator cell = mesh.begin(0);
2219  cell != mesh.end(0);
2220  ++cell)
2221  {
2222  for (unsigned int i : cell->face_indices())
2223  {
2224  const typename MeshType::face_iterator face = cell->face(i);
2225  if (face->at_boundary() && face->boundary_id() == b_id1)
2226  {
2227  const std::pair<typename MeshType::cell_iterator, unsigned int>
2228  pair1 = std::make_pair(cell, i);
2229  pairs1.insert(pair1);
2230  }
2231 
2232  if (face->at_boundary() && face->boundary_id() == b_id2)
2233  {
2234  const std::pair<typename MeshType::cell_iterator, unsigned int>
2235  pair2 = std::make_pair(cell, i);
2236  pairs2.insert(pair2);
2237  }
2238  }
2239  }
2240 
2241  // Assure that all faces are matched on the coare grid. This requirement
2242  // can only fulfilled if a process owns the complete coarse grid. This is
2243  // not the case for a parallel::fullydistributed::Triangulation, i.e. this
2244  // requirement has not to be met neither (consider faces on the outside of a
2245  // ghost cell that are periodic but for which the ghost neighbor doesn't
2246  // exist).
2247  if (!(((pairs1.size() > 0) &&
2248  (dynamic_cast<
2250  *>(&pairs1.begin()->first->get_triangulation()) != nullptr)) ||
2251  ((pairs2.size() > 0) &&
2252  (dynamic_cast<
2254  *>(&pairs2.begin()->first->get_triangulation()) != nullptr))))
2255  Assert(pairs1.size() == pairs2.size(),
2256  ExcMessage("Unmatched faces on periodic boundaries"));
2257 
2258  Assert(
2259  (pairs1.size() > 0 ||
2260  (dynamic_cast<
2262  &mesh.begin()->get_triangulation()) != nullptr)),
2263  ExcMessage("No new periodic face pairs have been found. "
2264  "Are you sure that you've selected the correct boundary "
2265  "id's and that the coarsest level mesh is colorized?"));
2266 
2267  // and call match_periodic_face_pairs that does the actual matching:
2269  pairs1, pairs2, direction, matched_pairs, offset, matrix);
2270  }
2271 
2272 
2273 
2274  /*
2275  * Internally used in orthogonal_equality
2276  *
2277  * An orthogonal equality test for points:
2278  *
2279  * point1 and point2 are considered equal, if
2280  * matrix.point1 + offset - point2
2281  * is parallel to the unit vector in <direction>
2282  */
2283  template <int spacedim>
2284  inline bool
2286  const Point<spacedim> & point2,
2287  const int direction,
2288  const Tensor<1, spacedim> &offset,
2289  const FullMatrix<double> & matrix)
2290  {
2291  AssertIndexRange(direction, spacedim);
2292 
2293  Assert(matrix.m() == matrix.n(), ExcInternalError());
2294 
2295  Point<spacedim> distance;
2296 
2297  if (matrix.m() == spacedim)
2298  for (int i = 0; i < spacedim; ++i)
2299  for (int j = 0; j < spacedim; ++j)
2300  distance(i) += matrix(i, j) * point1(j);
2301  else
2302  distance = point1;
2303 
2304  distance += offset - point2;
2305 
2306  for (int i = 0; i < spacedim; ++i)
2307  {
2308  // Only compare coordinate-components != direction:
2309  if (i == direction)
2310  continue;
2311 
2312  if (std::abs(distance(i)) > 1.e-10)
2313  return false;
2314  }
2315 
2316  return true;
2317  }
2318 
2319 
2320  /*
2321  * Internally used in orthogonal_equality
2322  *
2323  * A lookup table to transform vertex matchings to orientation flags of
2324  * the form (face_orientation, face_flip, face_rotation)
2325  *
2326  * See the comment on the next function as well as the detailed
2327  * documentation of make_periodicity_constraints and
2328  * collect_periodic_faces for details
2329  */
2330  template <int dim>
2332  {};
2333 
2334  template <>
2336  {
2337  using MATCH_T =
2338  std::array<unsigned int, GeometryInfo<1>::vertices_per_face>;
2339  static inline std::bitset<3>
2340  lookup(const MATCH_T &)
2341  {
2342  // The 1D case is trivial
2343  return 1; // [true ,false,false]
2344  }
2345  };
2346 
2347  template <>
2349  {
2350  using MATCH_T =
2351  std::array<unsigned int, GeometryInfo<2>::vertices_per_face>;
2352  static inline std::bitset<3>
2353  lookup(const MATCH_T &matching)
2354  {
2355  // In 2D matching faces (=lines) results in two cases: Either
2356  // they are aligned or flipped. We store this "line_flip"
2357  // property somewhat sloppy as "face_flip"
2358  // (always: face_orientation = true, face_rotation = false)
2359 
2360  static const MATCH_T m_tff = {{0, 1}};
2361  if (matching == m_tff)
2362  return 1; // [true ,false,false]
2363  static const MATCH_T m_ttf = {{1, 0}};
2364  if (matching == m_ttf)
2365  return 3; // [true ,true ,false]
2366  Assert(false, ExcInternalError());
2367  // what follows is dead code, but it avoids warnings about the lack
2368  // of a return value
2369  return 0;
2370  }
2371  };
2372 
2373  template <>
2375  {
2376  using MATCH_T =
2377  std::array<unsigned int, GeometryInfo<3>::vertices_per_face>;
2378 
2379  static inline std::bitset<3>
2380  lookup(const MATCH_T &matching)
2381  {
2382  // The full fledged 3D case. *Yay*
2383  // See the documentation in include/deal.II/base/geometry_info.h
2384  // as well as the actual implementation in source/grid/tria.cc
2385  // for more details...
2386 
2387  static const MATCH_T m_tff = {{0, 1, 2, 3}};
2388  if (matching == m_tff)
2389  return 1; // [true ,false,false]
2390  static const MATCH_T m_tft = {{1, 3, 0, 2}};
2391  if (matching == m_tft)
2392  return 5; // [true ,false,true ]
2393  static const MATCH_T m_ttf = {{3, 2, 1, 0}};
2394  if (matching == m_ttf)
2395  return 3; // [true ,true ,false]
2396  static const MATCH_T m_ttt = {{2, 0, 3, 1}};
2397  if (matching == m_ttt)
2398  return 7; // [true ,true ,true ]
2399  static const MATCH_T m_fff = {{0, 2, 1, 3}};
2400  if (matching == m_fff)
2401  return 0; // [false,false,false]
2402  static const MATCH_T m_fft = {{2, 3, 0, 1}};
2403  if (matching == m_fft)
2404  return 4; // [false,false,true ]
2405  static const MATCH_T m_ftf = {{3, 1, 2, 0}};
2406  if (matching == m_ftf)
2407  return 2; // [false,true ,false]
2408  static const MATCH_T m_ftt = {{1, 0, 3, 2}};
2409  if (matching == m_ftt)
2410  return 6; // [false,true ,true ]
2411  Assert(false, ExcInternalError());
2412  // what follows is dead code, but it avoids warnings about the lack
2413  // of a return value
2414  return 0;
2415  }
2416  };
2417 
2418 
2419 
2420  template <typename FaceIterator>
2421  inline bool orthogonal_equality(
2422  std::bitset<3> & orientation,
2423  const FaceIterator & face1,
2424  const FaceIterator & face2,
2425  const int direction,
2427  const FullMatrix<double> & matrix)
2428  {
2429  Assert(matrix.m() == matrix.n(),
2430  ExcMessage("The supplied matrix must be a square matrix"));
2431 
2432  static const int dim = FaceIterator::AccessorType::dimension;
2433 
2434  // Do a full matching of the face vertices:
2435 
2436  std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
2437 
2438  AssertDimension(face1->n_vertices(), face2->n_vertices());
2439 
2440  std::set<unsigned int> face2_vertices;
2441  for (unsigned int i = 0; i < face1->n_vertices(); ++i)
2442  face2_vertices.insert(i);
2443 
2444  for (unsigned int i = 0; i < face1->n_vertices(); ++i)
2445  {
2446  for (std::set<unsigned int>::iterator it = face2_vertices.begin();
2447  it != face2_vertices.end();
2448  ++it)
2449  {
2450  if (orthogonal_equality(face1->vertex(i),
2451  face2->vertex(*it),
2452  direction,
2453  offset,
2454  matrix))
2455  {
2456  matching[i] = *it;
2457  face2_vertices.erase(it);
2458  break; // jump out of the innermost loop
2459  }
2460  }
2461  }
2462 
2463  // And finally, a lookup to determine the ordering bitmask:
2464  if (face2_vertices.empty())
2465  orientation = OrientationLookupTable<dim>::lookup(matching);
2466 
2467  return face2_vertices.empty();
2468  }
2469 
2470 
2471 
2472  template <typename FaceIterator>
2473  inline bool
2475  const FaceIterator & face1,
2476  const FaceIterator & face2,
2477  const int direction,
2479  const FullMatrix<double> & matrix)
2480  {
2481  // Call the function above with a dummy orientation array
2482  std::bitset<3> dummy;
2483  return orthogonal_equality(dummy, face1, face2, direction, offset, matrix);
2484  }
2485 
2486 
2487 
2488 } // namespace GridTools
2489 
2490 
2491 #include "grid_tools_dof_handlers.inst"
2492 
2493 
size_type m() const
unsigned int n_active_cells() const
Definition: tria.cc:12697
static std::bitset< 3 > lookup(const MATCH_T &matching)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
void clear_user_flags()
Definition: tria.cc:11171
Contents is actually a matrix.
bool have_same_coarse_mesh(const Triangulation< dim, spacedim > &mesh_1, const Triangulation< dim, spacedim > &mesh_2)
unsigned int n_cells() const
Definition: tria.cc:12689
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:393
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:240
bool orthogonal_equality(std::bitset< 3 > &orientation, const FaceIterator &face1, const FaceIterator &face2, const int direction, const Tensor< 1, FaceIterator::AccessorType::space_dimension > &offset=Tensor< 1, FaceIterator::AccessorType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
std::array< unsigned int, GeometryInfo< 2 >::vertices_per_face > MATCH_T
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12207
cell_iterator end() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
std::array< unsigned int, GeometryInfo< 3 >::vertices_per_face > MATCH_T
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
Point< 2 > second
Definition: grid_out.cc:4576
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
std::vector< typename MeshType::cell_iterator > compute_cell_halo_layer_on_level(const MeshType &mesh, const std::function< bool(const typename MeshType::cell_iterator &)> &predicate, const unsigned int level)
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11990
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:5362
static std::bitset< 3 > lookup(const MATCH_T &)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
static double distance_to_unit_cell(const Point< dim > &p)
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)
cell_iterator end() const
Definition: tria.cc:12101
active_cell_iterator begin_active(const unsigned int level=0) const
void match_periodic_face_pairs(std::set< std::pair< CellIterator, unsigned int >> &pairs1, std::set< std::pair< typename identity< CellIterator >::type, unsigned int >> &pairs2, const int direction, std::vector< PeriodicFacePair< CellIterator >> &matched_pairs, const ::Tensor< 1, CellIterator::AccessorType::space_dimension > &offset, const FullMatrix< double > &matrix)
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:13364
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
std::map< types::global_dof_index, std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > > get_dof_to_support_patch_map(DoFHandler< dim, spacedim > &dof_handler)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:10485
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:303
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:394
unsigned int level
Definition: grid_out.cc:4578
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={}, const double tolerance=1.e-10)
Point< 3 > vertices[4]
void build_triangulation_from_patch(const std::vector< typename Container::active_cell_iterator > &patch, Triangulation< Container::dimension, Container::space_dimension > &local_triangulation, std::map< typename Triangulation< Container::dimension, Container::space_dimension >::active_cell_iterator, typename Container::active_cell_iterator > &patch_to_global_tria_map)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12647
std::array< unsigned int, GeometryInfo< 1 >::vertices_per_face > MATCH_T
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12640
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Point< 2 > first
Definition: grid_out.cc:4575
unsigned int global_dof_index
Definition: types.h:76
__global__ void set(Number *val, const Number s, const size_type N)
typename ActiveSelector::line_iterator line_iterator
Definition: dof_handler.h:357
std::vector< std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > > find_all_active_cells_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const double tolerance, const std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim >> &first_cell)
const std::vector< bool > & get_used_vertices() const
Definition: tria.cc:13241
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:393
T min(const T &t, const MPI_Comm &mpi_communicator)
std::vector< typename MeshType< dim, spacedim >::active_cell_iterator > find_cells_adjacent_to_vertex(const MeshType< dim, spacedim > &container, const unsigned int vertex_index)
Definition: grid_tools.cc:1479
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator >> &matched_pairs, const Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
Point< 3 > center
virtual bool preserves_vertex_locations() const =0
static ::ExceptionBase & ExcNotImplemented()
Iterator points to a valid object.
static std::bitset< 3 > lookup(const MATCH_T &matching)
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
unsigned int size() const
T max(const T &t, const MPI_Comm &mpi_communicator)
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_layer_within_distance(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const double layer_thickness)
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:5381
virtual void clear()
Definition: tria.cc:10180
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
Triangulation< dim, spacedim > & get_triangulation()
Definition: tria.cc:13317