Reference documentation for deal.II version Git 6d836ad036 2021-06-16 20:57:45 +0200
\(\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 - 2021 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>
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>
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  // Initialize best_cell.first to the end iterator
454  best_cell.first = mesh.end();
455 
456  // Find closest vertex and determine
457  // all adjacent cells
458  std::vector<active_cell_iterator> adjacent_cells_tmp =
460  mesh, find_closest_vertex(mapping, mesh, p, marked_vertices));
461 
462  // Make sure that we have found
463  // at least one cell adjacent to vertex.
464  Assert(adjacent_cells_tmp.size() > 0, ExcInternalError());
465 
466  // Copy all the cells into a std::set
467  std::set<active_cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(),
468  adjacent_cells_tmp.end());
469  std::set<active_cell_iterator> searched_cells;
470 
471  // Determine the maximal number of cells
472  // in the grid.
473  // As long as we have not found
474  // the cell and have not searched
475  // every cell in the triangulation,
476  // we keep on looking.
477  const unsigned int n_active_cells =
478  mesh.get_triangulation().n_active_cells();
479  bool found = false;
480  unsigned int cells_searched = 0;
481  while (!found && cells_searched < n_active_cells)
482  {
483  typename std::set<active_cell_iterator>::const_iterator
484  cell = adjacent_cells.begin(),
485  endc = adjacent_cells.end();
486  for (; cell != endc; ++cell)
487  {
488  if ((*cell)->is_artificial() == false)
489  {
490  try
491  {
492  const Point<dim> p_cell =
493  mapping.transform_real_to_unit_cell(*cell, p);
494 
495  // calculate the infinity norm of
496  // the distance vector to the unit cell.
497  const double dist =
499 
500  // We compare if the point is inside the
501  // unit cell (or at least not too far
502  // outside). If it is, it is also checked
503  // that the cell has a more refined state
504  if ((dist < best_distance) ||
505  ((dist == best_distance) &&
506  ((*cell)->level() > best_level)))
507  {
508  found = true;
509  best_distance = dist;
510  best_level = (*cell)->level();
511  best_cell = std::make_pair(*cell, p_cell);
512  }
513  }
514  catch (
515  typename MappingQGeneric<dim,
516  spacedim>::ExcTransformationFailed &)
517  {
518  // ok, the transformation
519  // failed presumably
520  // because the point we
521  // are looking for lies
522  // outside the current
523  // cell. this means that
524  // the current cell can't
525  // be the cell around the
526  // point, so just ignore
527  // this cell and move on
528  // to the next
529  }
530  }
531  }
532 
533  // update the number of cells searched
534  cells_searched += adjacent_cells.size();
535 
536  // if the user provided a custom mask for vertices,
537  // terminate the search without trying to expand the search
538  // to all cells of the triangulation, as done below.
539  if (marked_vertices.size() > 0)
540  cells_searched = n_active_cells;
541 
542  // if we have not found the cell in
543  // question and have not yet searched every
544  // cell, we expand our search to
545  // all the not already searched neighbors of
546  // the cells in adjacent_cells. This is
547  // what find_active_cell_around_point_internal
548  // is for.
549  if (!found && cells_searched < n_active_cells)
550  {
551  find_active_cell_around_point_internal<dim, MeshType, spacedim>(
552  mesh, searched_cells, adjacent_cells);
553  }
554  }
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  const auto cell_and_point = find_active_cell_around_point(
578  mapping, mesh, p, marked_vertices, tolerance);
579 
580  if (cell_and_point.first == mesh.end())
581  return {};
582 
584  mapping, mesh, p, tolerance, cell_and_point);
585  }
586 
587 
588 
589  template <int dim, template <int, int> class MeshType, int spacedim>
590 #ifndef _MSC_VER
591  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
592  Point<dim>>>
593 #else
594  std::vector<std::pair<
595  typename ::internal::
596  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
597  Point<dim>>>
598 #endif
600  const Mapping<dim, spacedim> & mapping,
601  const MeshType<dim, spacedim> &mesh,
602  const Point<spacedim> & p,
603  const double tolerance,
604  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
605  Point<dim>> & first_cell)
606  {
607  std::vector<
608  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
609  Point<dim>>>
610  cells_and_points;
611 
612  // insert the fist cell and point into the vector
613  cells_and_points.push_back(first_cell);
614 
615  // check if the given point is on the surface of the unit cell. If yes,
616  // need to find all neighbors
617  const Point<dim> unit_point = cells_and_points.front().second;
618  const auto my_cell = cells_and_points.front().first;
619  Tensor<1, dim> distance_to_center;
620  unsigned int n_dirs_at_threshold = 0;
621  unsigned int last_point_at_threshold = numbers::invalid_unsigned_int;
622  for (unsigned int d = 0; d < dim; ++d)
623  {
624  distance_to_center[d] = std::abs(unit_point[d] - 0.5);
625  if (distance_to_center[d] > 0.5 - tolerance)
626  {
627  ++n_dirs_at_threshold;
628  last_point_at_threshold = d;
629  }
630  }
631 
632  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
633  cells_to_add;
634  // point is within face -> only need neighbor
635  if (n_dirs_at_threshold == 1)
636  {
637  unsigned int neighbor_index =
638  2 * last_point_at_threshold +
639  (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
640  if (!my_cell->at_boundary(neighbor_index))
641  cells_to_add.push_back(my_cell->neighbor(neighbor_index));
642  }
643  // corner point -> use all neighbors
644  else if (n_dirs_at_threshold == dim)
645  {
646  unsigned int local_vertex_index = 0;
647  for (unsigned int d = 0; d < dim; ++d)
648  local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
649  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
651  mesh, my_cell->vertex_index(local_vertex_index));
652  for (const auto &cell : cells)
653  if (cell != my_cell)
654  cells_to_add.push_back(cell);
655  }
656  // point on line in 3D: We cannot simply take the intersection between
657  // the two vertices of cells because of hanging nodes. So instead we
658  // list the vertices around both points and then select the
659  // appropriate cells according to the result of read_to_unit_cell
660  // below.
661  else if (n_dirs_at_threshold == 2)
662  {
663  std::pair<unsigned int, unsigned int> vertex_indices[3];
664  unsigned int count_vertex_indices = 0;
665  unsigned int free_direction = numbers::invalid_unsigned_int;
666  for (unsigned int d = 0; d < dim; ++d)
667  {
668  if (distance_to_center[d] > 0.5 - tolerance)
669  {
670  vertex_indices[count_vertex_indices].first = d;
671  vertex_indices[count_vertex_indices].second =
672  unit_point[d] > 0.5 ? 1 : 0;
673  count_vertex_indices++;
674  }
675  else
676  free_direction = d;
677  }
678 
679  AssertDimension(count_vertex_indices, 2);
680  Assert(free_direction != numbers::invalid_unsigned_int,
681  ExcInternalError());
682 
683  const unsigned int first_vertex =
684  (vertex_indices[0].second << vertex_indices[0].first) +
685  (vertex_indices[1].second << vertex_indices[1].first);
686  for (unsigned int d = 0; d < 2; ++d)
687  {
688  auto tentative_cells = find_cells_adjacent_to_vertex(
689  mesh,
690  my_cell->vertex_index(first_vertex + (d << free_direction)));
691  for (const auto &cell : tentative_cells)
692  {
693  bool cell_not_yet_present = true;
694  for (const auto &other_cell : cells_to_add)
695  if (cell == other_cell)
696  {
697  cell_not_yet_present = false;
698  break;
699  }
700  if (cell_not_yet_present)
701  cells_to_add.push_back(cell);
702  }
703  }
704  }
705 
706  const double original_distance_to_unit_cell =
708  for (const auto &cell : cells_to_add)
709  {
710  if (cell != my_cell)
711  try
712  {
713  const Point<dim> p_unit =
714  mapping.transform_real_to_unit_cell(cell, p);
716  original_distance_to_unit_cell + tolerance)
717  cells_and_points.emplace_back(cell, p_unit);
718  }
719  catch (typename Mapping<dim>::ExcTransformationFailed &)
720  {}
721  }
722 
723  std::sort(
724  cells_and_points.begin(),
725  cells_and_points.end(),
726  [](const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
727  Point<dim>> &a,
728  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
729  Point<dim>> &b) { return a.first < b.first; });
730 
731  return cells_and_points;
732  }
733 
734 
735 
736  template <class MeshType>
737  std::vector<typename MeshType::active_cell_iterator>
739  const MeshType &mesh,
740  const std::function<bool(const typename MeshType::active_cell_iterator &)>
741  &predicate)
742  {
743  std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
744  std::vector<bool> locally_active_vertices_on_subdomain(
745  mesh.get_triangulation().n_vertices(), false);
746 
747  // Find the cells for which the predicate is true
748  // These are the cells around which we wish to construct
749  // the halo layer
750  for (const auto &cell : mesh.active_cell_iterators())
751  if (predicate(cell)) // True predicate --> Part of subdomain
752  for (const auto v : cell->vertex_indices())
753  locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true;
754 
755  // Find the cells that do not conform to the predicate
756  // but share a vertex with the selected subdomain
757  // These comprise the halo layer
758  for (const auto &cell : mesh.active_cell_iterators())
759  if (!predicate(cell)) // False predicate --> Potential halo cell
760  for (const auto v : cell->vertex_indices())
761  if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
762  true)
763  {
764  active_halo_layer.push_back(cell);
765  break;
766  }
767 
768  return active_halo_layer;
769  }
770 
771 
772 
773  template <class MeshType>
774  std::vector<typename MeshType::cell_iterator>
776  const MeshType &mesh,
777  const std::function<bool(const typename MeshType::cell_iterator &)>
778  & predicate,
779  const unsigned int level)
780  {
781  std::vector<typename MeshType::cell_iterator> level_halo_layer;
782  std::vector<bool> locally_active_vertices_on_level_subdomain(
783  mesh.get_triangulation().n_vertices(), false);
784 
785  // Find the cells for which the predicate is true
786  // These are the cells around which we wish to construct
787  // the halo layer
788  for (typename MeshType::cell_iterator cell = mesh.begin(level);
789  cell != mesh.end(level);
790  ++cell)
791  if (predicate(cell)) // True predicate --> Part of subdomain
792  for (const unsigned int v : cell->vertex_indices())
793  locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
794  true;
795 
796  // Find the cells that do not conform to the predicate
797  // but share a vertex with the selected subdomain on that level
798  // These comprise the halo layer
799  for (typename MeshType::cell_iterator cell = mesh.begin(level);
800  cell != mesh.end(level);
801  ++cell)
802  if (!predicate(cell)) // False predicate --> Potential halo cell
803  for (const unsigned int v : cell->vertex_indices())
804  if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
805  v)] == true)
806  {
807  level_halo_layer.push_back(cell);
808  break;
809  }
810 
811  return level_halo_layer;
812  }
813 
814 
815  namespace
816  {
817  template <class MeshType>
818  bool
819  contains_locally_owned_cells(
820  const std::vector<typename MeshType::active_cell_iterator> &cells)
821  {
822  for (typename std::vector<
823  typename MeshType::active_cell_iterator>::const_iterator it =
824  cells.begin();
825  it != cells.end();
826  ++it)
827  {
828  if ((*it)->is_locally_owned())
829  return true;
830  }
831  return false;
832  }
833 
834  template <class MeshType>
835  bool
836  contains_artificial_cells(
837  const std::vector<typename MeshType::active_cell_iterator> &cells)
838  {
839  for (typename std::vector<
840  typename MeshType::active_cell_iterator>::const_iterator it =
841  cells.begin();
842  it != cells.end();
843  ++it)
844  {
845  if ((*it)->is_artificial())
846  return true;
847  }
848  return false;
849  }
850  } // namespace
851 
852 
853 
854  template <class MeshType>
855  std::vector<typename MeshType::active_cell_iterator>
856  compute_ghost_cell_halo_layer(const MeshType &mesh)
857  {
858  std::function<bool(const typename MeshType::active_cell_iterator &)>
859  predicate = IteratorFilters::LocallyOwnedCell();
860 
861  const std::vector<typename MeshType::active_cell_iterator>
862  active_halo_layer = compute_active_cell_halo_layer(mesh, predicate);
863 
864  // Check that we never return locally owned or artificial cells
865  // What is left should only be the ghost cells
866  Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) == false,
867  ExcMessage("Halo layer contains locally owned cells"));
868  Assert(contains_artificial_cells<MeshType>(active_halo_layer) == false,
869  ExcMessage("Halo layer contains artificial cells"));
870 
871  return active_halo_layer;
872  }
873 
874 
875 
876  template <class MeshType>
877  std::vector<typename MeshType::active_cell_iterator>
879  const MeshType &mesh,
880  const std::function<bool(const typename MeshType::active_cell_iterator &)>
881  & predicate,
882  const double layer_thickness)
883  {
884  std::vector<typename MeshType::active_cell_iterator>
885  subdomain_boundary_cells, active_cell_layer_within_distance;
886  std::vector<bool> vertices_outside_subdomain(
887  mesh.get_triangulation().n_vertices(), false);
888 
889  const unsigned int spacedim = MeshType::space_dimension;
890 
891  unsigned int n_non_predicate_cells = 0; // Number of non predicate cells
892 
893  // Find the layer of cells for which predicate is true and that
894  // are on the boundary with other cells. These are
895  // subdomain boundary cells.
896 
897  // Find the cells for which the predicate is false
898  // These are the cells which are around the predicate subdomain
899  for (const auto &cell : mesh.active_cell_iterators())
900  if (!predicate(cell)) // Negation of predicate --> Not Part of subdomain
901  {
902  for (const unsigned int v : cell->vertex_indices())
903  vertices_outside_subdomain[cell->vertex_index(v)] = true;
904  n_non_predicate_cells++;
905  }
906 
907  // If all the active cells conform to the predicate
908  // or if none of the active cells conform to the predicate
909  // there is no active cell layer around the predicate
910  // subdomain (within any distance)
911  if (n_non_predicate_cells == 0 ||
912  n_non_predicate_cells == mesh.get_triangulation().n_active_cells())
913  return std::vector<typename MeshType::active_cell_iterator>();
914 
915  // Find the cells that conform to the predicate
916  // but share a vertex with the cell not in the predicate subdomain
917  for (const auto &cell : mesh.active_cell_iterators())
918  if (predicate(cell)) // True predicate --> Potential boundary cell of the
919  // subdomain
920  for (const unsigned int v : cell->vertex_indices())
921  if (vertices_outside_subdomain[cell->vertex_index(v)] == true)
922  {
923  subdomain_boundary_cells.push_back(cell);
924  break; // No need to go through remaining vertices
925  }
926 
927  // To cheaply filter out some cells located far away from the predicate
928  // subdomain, get the bounding box of the predicate subdomain.
929  std::pair<Point<spacedim>, Point<spacedim>> bounding_box =
930  compute_bounding_box(mesh, predicate);
931 
932  // DOUBLE_EPSILON to compare really close double values
933  const double DOUBLE_EPSILON = 100. * std::numeric_limits<double>::epsilon();
934 
935  // Add layer_thickness to the bounding box
936  for (unsigned int d = 0; d < spacedim; ++d)
937  {
938  bounding_box.first[d] -= (layer_thickness + DOUBLE_EPSILON); // minp
939  bounding_box.second[d] += (layer_thickness + DOUBLE_EPSILON); // maxp
940  }
941 
942  std::vector<Point<spacedim>>
943  subdomain_boundary_cells_centers; // cache all the subdomain boundary
944  // cells centers here
945  std::vector<double>
946  subdomain_boundary_cells_radii; // cache all the subdomain boundary cells
947  // radii
948  subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
949  subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
950  // compute cell radius for each boundary cell of the predicate subdomain
951  for (typename std::vector<typename MeshType::active_cell_iterator>::
952  const_iterator subdomain_boundary_cell_iterator =
953  subdomain_boundary_cells.begin();
954  subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
955  ++subdomain_boundary_cell_iterator)
956  {
957  const std::pair<Point<spacedim>, double>
958  &subdomain_boundary_cell_enclosing_ball =
959  (*subdomain_boundary_cell_iterator)->enclosing_ball();
960 
961  subdomain_boundary_cells_centers.push_back(
962  subdomain_boundary_cell_enclosing_ball.first);
963  subdomain_boundary_cells_radii.push_back(
964  subdomain_boundary_cell_enclosing_ball.second);
965  }
966  AssertThrow(subdomain_boundary_cells_radii.size() ==
967  subdomain_boundary_cells_centers.size(),
968  ExcInternalError());
969 
970  // Find the cells that are within layer_thickness of predicate subdomain
971  // boundary distance but are inside the extended bounding box. Most cells
972  // might be outside the extended bounding box, so we could skip them. Those
973  // cells that are inside the extended bounding box but are not part of the
974  // predicate subdomain are possible candidates to be within the distance to
975  // the boundary cells of the predicate subdomain.
976  for (const auto &cell : mesh.active_cell_iterators())
977  {
978  // Ignore all the cells that are in the predicate subdomain
979  if (predicate(cell))
980  continue;
981 
982  const std::pair<Point<spacedim>, double> &cell_enclosing_ball =
983  cell->enclosing_ball();
984 
985  const Point<spacedim> cell_enclosing_ball_center =
986  cell_enclosing_ball.first;
987  const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
988 
989  bool cell_inside = true; // reset for each cell
990 
991  for (unsigned int d = 0; d < spacedim; ++d)
992  cell_inside &=
993  (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius >
994  bounding_box.first[d]) &&
995  (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius <
996  bounding_box.second[d]);
997  // cell_inside is true if its enclosing ball intersects the extended
998  // bounding box
999 
1000  // Ignore all the cells that are outside the extended bounding box
1001  if (cell_inside)
1002  for (unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1003  ++i)
1004  if (cell_enclosing_ball_center.distance_square(
1005  subdomain_boundary_cells_centers[i]) <
1006  Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1007  subdomain_boundary_cells_radii[i] +
1008  layer_thickness + DOUBLE_EPSILON))
1009  {
1010  active_cell_layer_within_distance.push_back(cell);
1011  break; // Exit the loop checking all the remaining subdomain
1012  // boundary cells
1013  }
1014  }
1015  return active_cell_layer_within_distance;
1016  }
1017 
1018 
1019 
1020  template <class MeshType>
1021  std::vector<typename MeshType::active_cell_iterator>
1023  const double layer_thickness)
1024  {
1025  IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate;
1026  std::function<bool(const typename MeshType::active_cell_iterator &)>
1027  predicate(locally_owned_cell_predicate);
1028 
1029  const std::vector<typename MeshType::active_cell_iterator>
1030  ghost_cell_layer_within_distance =
1032  predicate,
1033  layer_thickness);
1034 
1035  // Check that we never return locally owned or artificial cells
1036  // What is left should only be the ghost cells
1037  Assert(
1038  contains_locally_owned_cells<MeshType>(
1039  ghost_cell_layer_within_distance) == false,
1040  ExcMessage(
1041  "Ghost cells within layer_thickness contains locally owned cells."));
1042  Assert(
1043  contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1044  false,
1045  ExcMessage(
1046  "Ghost cells within layer_thickness contains artificial cells. "
1047  "The function compute_ghost_cell_layer_within_distance "
1048  "is probably called while using parallel::distributed::Triangulation. "
1049  "In such case please refer to the description of this function."));
1050 
1051  return ghost_cell_layer_within_distance;
1052  }
1053 
1054 
1055 
1056  template <class MeshType>
1057  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1059  const MeshType &mesh,
1060  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1061  &predicate)
1062  {
1063  std::vector<bool> locally_active_vertices_on_subdomain(
1064  mesh.get_triangulation().n_vertices(), false);
1065 
1066  const unsigned int spacedim = MeshType::space_dimension;
1067 
1068  // Two extreme points can define the bounding box
1069  // around the active cells that conform to the given predicate.
1071 
1072  // initialize minp and maxp with the first predicate cell center
1073  for (const auto &cell : mesh.active_cell_iterators())
1074  if (predicate(cell))
1075  {
1076  minp = cell->center();
1077  maxp = cell->center();
1078  break;
1079  }
1080 
1081  // Run through all the cells to check if it belongs to predicate domain,
1082  // if it belongs to the predicate domain, extend the bounding box.
1083  for (const auto &cell : mesh.active_cell_iterators())
1084  if (predicate(cell)) // True predicate --> Part of subdomain
1085  for (const unsigned int v : cell->vertex_indices())
1086  if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1087  false)
1088  {
1089  locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1090  true;
1091  for (unsigned int d = 0; d < spacedim; ++d)
1092  {
1093  minp[d] = std::min(minp[d], cell->vertex(v)[d]);
1094  maxp[d] = std::max(maxp[d], cell->vertex(v)[d]);
1095  }
1096  }
1097 
1098  return std::make_pair(minp, maxp);
1099  }
1100 
1101 
1102 
1103  template <typename MeshType>
1104  std::list<std::pair<typename MeshType::cell_iterator,
1105  typename MeshType::cell_iterator>>
1106  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
1107  {
1108  Assert(have_same_coarse_mesh(mesh_1, mesh_2),
1109  ExcMessage("The two meshes must be represent triangulations that "
1110  "have the same coarse meshes"));
1111  // We will allow the output to contain ghost cells when we have shared
1112  // Triangulations (i.e., so that each processor will get exactly the same
1113  // list of cell pairs), but not when we have two distributed
1114  // Triangulations (so that all active cells are partitioned by processor).
1115  // Non-parallel Triangulations have no ghost or artificial cells, so they
1116  // work the same way as shared Triangulations here.
1117  bool remove_ghost_cells = false;
1118 #ifdef DEAL_II_WITH_MPI
1119  {
1120  constexpr int dim = MeshType::dimension;
1121  constexpr int spacedim = MeshType::space_dimension;
1123  *>(&mesh_1.get_triangulation()) != nullptr ||
1125  *>(&mesh_2.get_triangulation()) != nullptr)
1126  {
1127  Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1128  ExcMessage("This function can only be used with meshes "
1129  "corresponding to distributed Triangulations when "
1130  "both Triangulations are equal."));
1131  remove_ghost_cells = true;
1132  }
1133  }
1134 #endif
1135 
1136  // the algorithm goes as follows: first, we fill a list with pairs of
1137  // iterators common to the two meshes on the coarsest level. then we
1138  // traverse the list; each time, we find a pair of iterators for which
1139  // both correspond to non-active cells, we delete this item and push the
1140  // pairs of iterators to their children to the back. if these again both
1141  // correspond to non-active cells, we will get to the later on for further
1142  // consideration
1143  using CellList = std::list<std::pair<typename MeshType::cell_iterator,
1144  typename MeshType::cell_iterator>>;
1145  CellList cell_list;
1146 
1147  // first push the coarse level cells
1148  typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1149  cell_2 = mesh_2.begin(0);
1150  for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1151  cell_list.emplace_back(cell_1, cell_2);
1152 
1153  // then traverse list as described above
1154  typename CellList::iterator cell_pair = cell_list.begin();
1155  while (cell_pair != cell_list.end())
1156  {
1157  // if both cells in this pair have children, then erase this element
1158  // and push their children instead
1159  if (cell_pair->first->has_children() &&
1160  cell_pair->second->has_children())
1161  {
1162  Assert(cell_pair->first->refinement_case() ==
1163  cell_pair->second->refinement_case(),
1164  ExcNotImplemented());
1165  for (unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1166  cell_list.emplace_back(cell_pair->first->child(c),
1167  cell_pair->second->child(c));
1168 
1169  // erasing an iterator keeps other iterators valid, so already
1170  // advance the present iterator by one and then delete the element
1171  // we've visited before
1172  const auto previous_cell_pair = cell_pair;
1173  ++cell_pair;
1174  cell_list.erase(previous_cell_pair);
1175  }
1176  else
1177  {
1178  // at least one cell is active
1179  if (remove_ghost_cells &&
1180  ((cell_pair->first->is_active() &&
1181  !cell_pair->first->is_locally_owned()) ||
1182  (cell_pair->second->is_active() &&
1183  !cell_pair->second->is_locally_owned())))
1184  {
1185  // we only exclude ghost cells for distributed Triangulations
1186  const auto previous_cell_pair = cell_pair;
1187  ++cell_pair;
1188  cell_list.erase(previous_cell_pair);
1189  }
1190  else
1191  ++cell_pair;
1192  }
1193  }
1194 
1195  // just to make sure everything is ok, validate that all pairs have at
1196  // least one active iterator or have different refinement_cases
1197  for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1198  ++cell_pair)
1199  Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1200  (cell_pair->first->refinement_case() !=
1201  cell_pair->second->refinement_case()),
1202  ExcInternalError());
1203 
1204  return cell_list;
1205  }
1206 
1207 
1208 
1209  template <int dim, int spacedim>
1210  bool
1212  const Triangulation<dim, spacedim> &mesh_2)
1213  {
1214  // make sure the two meshes have
1215  // the same number of coarse cells
1216  if (mesh_1.n_cells(0) != mesh_2.n_cells(0))
1217  return false;
1218 
1219  // if so, also make sure they have
1220  // the same vertices on the cells
1221  // of the coarse mesh
1223  mesh_1.begin(0),
1224  cell_2 =
1225  mesh_2.begin(0),
1226  endc = mesh_1.end(0);
1227  for (; cell_1 != endc; ++cell_1, ++cell_2)
1228  {
1229  if (cell_1->n_vertices() != cell_2->n_vertices())
1230  return false;
1231  for (const unsigned int v : cell_1->vertex_indices())
1232  if (cell_1->vertex(v) != cell_2->vertex(v))
1233  return false;
1234  }
1235 
1236  // if we've gotten through all
1237  // this, then the meshes really
1238  // seem to have a common coarse
1239  // mesh
1240  return true;
1241  }
1242 
1243 
1244 
1245  template <typename MeshType>
1246  bool
1247  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2)
1248  {
1249  return have_same_coarse_mesh(mesh_1.get_triangulation(),
1250  mesh_2.get_triangulation());
1251  }
1252 
1253 
1254 
1255  template <int dim, int spacedim>
1256  std::pair<typename ::DoFHandler<dim, spacedim>::active_cell_iterator,
1257  Point<dim>>
1259  const hp::MappingCollection<dim, spacedim> &mapping,
1260  const DoFHandler<dim, spacedim> & mesh,
1261  const Point<spacedim> & p,
1262  const double tolerance)
1263  {
1264  Assert((mapping.size() == 1) ||
1265  (mapping.size() == mesh.get_fe_collection().size()),
1266  ExcMessage("Mapping collection needs to have either size 1 "
1267  "or size equal to the number of elements in "
1268  "the FECollection."));
1269 
1270  using cell_iterator =
1271  typename ::DoFHandler<dim, spacedim>::active_cell_iterator;
1272 
1273  std::pair<cell_iterator, Point<dim>> best_cell;
1274  // If we have only one element in the MappingCollection,
1275  // we use find_active_cell_around_point using only one
1276  // mapping.
1277  if (mapping.size() == 1)
1278  {
1279  const std::vector<bool> marked_vertices = {};
1280  best_cell = find_active_cell_around_point(
1281  mapping[0], mesh, p, marked_vertices, tolerance);
1282  }
1283  else
1284  {
1285  // The best distance is set to the
1286  // maximum allowable distance from
1287  // the unit cell
1288  double best_distance = tolerance;
1289  int best_level = -1;
1290 
1291 
1292  // Find closest vertex and determine
1293  // all adjacent cells
1294  unsigned int vertex = find_closest_vertex(mesh, p);
1295 
1296  std::vector<cell_iterator> adjacent_cells_tmp =
1297  find_cells_adjacent_to_vertex(mesh, vertex);
1298 
1299  // Make sure that we have found
1300  // at least one cell adjacent to vertex.
1301  Assert(adjacent_cells_tmp.size() > 0, ExcInternalError());
1302 
1303  // Copy all the cells into a std::set
1304  std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(),
1305  adjacent_cells_tmp.end());
1306  std::set<cell_iterator> searched_cells;
1307 
1308  // Determine the maximal number of cells
1309  // in the grid.
1310  // As long as we have not found
1311  // the cell and have not searched
1312  // every cell in the triangulation,
1313  // we keep on looking.
1314  const unsigned int n_cells = mesh.get_triangulation().n_cells();
1315  bool found = false;
1316  unsigned int cells_searched = 0;
1317  while (!found && cells_searched < n_cells)
1318  {
1319  typename std::set<cell_iterator>::const_iterator
1320  cell = adjacent_cells.begin(),
1321  endc = adjacent_cells.end();
1322  for (; cell != endc; ++cell)
1323  {
1324  try
1325  {
1326  const Point<dim> p_cell =
1327  mapping[(*cell)->active_fe_index()]
1328  .transform_real_to_unit_cell(*cell, p);
1329 
1330 
1331  // calculate the infinity norm of
1332  // the distance vector to the unit cell.
1333  const double dist =
1335 
1336  // We compare if the point is inside the
1337  // unit cell (or at least not too far
1338  // outside). If it is, it is also checked
1339  // that the cell has a more refined state
1340  if (dist < best_distance || (dist == best_distance &&
1341  (*cell)->level() > best_level))
1342  {
1343  found = true;
1344  best_distance = dist;
1345  best_level = (*cell)->level();
1346  best_cell = std::make_pair(*cell, p_cell);
1347  }
1348  }
1349  catch (
1350  typename MappingQGeneric<dim,
1351  spacedim>::ExcTransformationFailed &)
1352  {
1353  // ok, the transformation
1354  // failed presumably
1355  // because the point we
1356  // are looking for lies
1357  // outside the current
1358  // cell. this means that
1359  // the current cell can't
1360  // be the cell around the
1361  // point, so just ignore
1362  // this cell and move on
1363  // to the next
1364  }
1365  }
1366  // update the number of cells searched
1367  cells_searched += adjacent_cells.size();
1368  // if we have not found the cell in
1369  // question and have not yet searched every
1370  // cell, we expand our search to
1371  // all the not already searched neighbors of
1372  // the cells in adjacent_cells.
1373  if (!found && cells_searched < n_cells)
1374  {
1375  find_active_cell_around_point_internal<dim,
1376  ::DoFHandler,
1377  spacedim>(
1378  mesh, searched_cells, adjacent_cells);
1379  }
1380  }
1381  }
1382 
1383  return best_cell;
1384  }
1385 
1386 
1387  template <class MeshType>
1388  std::vector<typename MeshType::active_cell_iterator>
1389  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
1390  {
1391  Assert(cell->is_locally_owned(),
1392  ExcMessage("This function only makes sense if the cell for "
1393  "which you are asking for a patch, is locally "
1394  "owned."));
1395 
1396  std::vector<typename MeshType::active_cell_iterator> patch;
1397  patch.push_back(cell);
1398  for (const unsigned int face_number : cell->face_indices())
1399  if (cell->face(face_number)->at_boundary() == false)
1400  {
1401  if (cell->neighbor(face_number)->has_children() == false)
1402  patch.push_back(cell->neighbor(face_number));
1403  else
1404  // the neighbor is refined. in 2d/3d, we can simply ask for the
1405  // children of the neighbor because they can not be further refined
1406  // and, consequently, the children is active
1407  if (MeshType::dimension > 1)
1408  {
1409  for (unsigned int subface = 0;
1410  subface < cell->face(face_number)->n_children();
1411  ++subface)
1412  patch.push_back(
1413  cell->neighbor_child_on_subface(face_number, subface));
1414  }
1415  else
1416  {
1417  // in 1d, we need to work a bit harder: iterate until we find
1418  // the child by going from cell to child to child etc
1419  typename MeshType::cell_iterator neighbor =
1420  cell->neighbor(face_number);
1421  while (neighbor->has_children())
1422  neighbor = neighbor->child(1 - face_number);
1423 
1424  Assert(neighbor->neighbor(1 - face_number) == cell,
1425  ExcInternalError());
1426  patch.push_back(neighbor);
1427  }
1428  }
1429  return patch;
1430  }
1431 
1432 
1433 
1434  template <class Container>
1435  std::vector<typename Container::cell_iterator>
1437  const std::vector<typename Container::active_cell_iterator> &patch)
1438  {
1439  Assert(patch.size() > 0,
1440  ExcMessage(
1441  "Vector containing patch cells should not be an empty vector!"));
1442  // In order to extract the set of cells with the coarsest common level from
1443  // the give vector of cells: First it finds the number associated with the
1444  // minimum level of refinement, namely "min_level"
1445  int min_level = patch[0]->level();
1446 
1447  for (unsigned int i = 0; i < patch.size(); ++i)
1448  min_level = std::min(min_level, patch[i]->level());
1449  std::set<typename Container::cell_iterator> uniform_cells;
1450  typename std::vector<
1451  typename Container::active_cell_iterator>::const_iterator patch_cell;
1452  // it loops through all cells of the input vector
1453  for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1454  {
1455  // If the refinement level of each cell i the loop be equal to the
1456  // min_level, so that that cell inserted into the set of uniform_cells,
1457  // as the set of cells with the coarsest common refinement level
1458  if ((*patch_cell)->level() == min_level)
1459  uniform_cells.insert(*patch_cell);
1460  else
1461  // If not, it asks for the parent of the cell, until it finds the
1462  // parent cell with the refinement level equal to the min_level and
1463  // inserts that parent cell into the the set of uniform_cells, as the
1464  // set of cells with the coarsest common refinement level.
1465  {
1466  typename Container::cell_iterator parent = *patch_cell;
1467 
1468  while (parent->level() > min_level)
1469  parent = parent->parent();
1470  uniform_cells.insert(parent);
1471  }
1472  }
1473 
1474  return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1475  uniform_cells.end());
1476  }
1477 
1478 
1479 
1480  template <class Container>
1481  void
1483  const std::vector<typename Container::active_cell_iterator> &patch,
1485  &local_triangulation,
1486  std::map<
1487  typename Triangulation<Container::dimension,
1488  Container::space_dimension>::active_cell_iterator,
1489  typename Container::active_cell_iterator> &patch_to_global_tria_map)
1490 
1491  {
1492  const std::vector<typename Container::cell_iterator> uniform_cells =
1493  get_cells_at_coarsest_common_level<Container>(patch);
1494  // First it creates triangulation from the vector of "uniform_cells"
1495  local_triangulation.clear();
1496  std::vector<Point<Container::space_dimension>> vertices;
1497  const unsigned int n_uniform_cells = uniform_cells.size();
1498  std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1499  unsigned int k = 0; // for enumerating cells
1500  unsigned int i = 0; // for enumerating vertices
1501  typename std::vector<typename Container::cell_iterator>::const_iterator
1502  uniform_cell;
1503  for (uniform_cell = uniform_cells.begin();
1504  uniform_cell != uniform_cells.end();
1505  ++uniform_cell)
1506  {
1507  for (const unsigned int v : (*uniform_cell)->vertex_indices())
1508  {
1510  (*uniform_cell)->vertex(v);
1511  bool repeat_vertex = false;
1512 
1513  for (unsigned int m = 0; m < i; ++m)
1514  {
1515  if (position == vertices[m])
1516  {
1517  repeat_vertex = true;
1518  cells[k].vertices[v] = m;
1519  break;
1520  }
1521  }
1522  if (repeat_vertex == false)
1523  {
1524  vertices.push_back(position);
1525  cells[k].vertices[v] = i;
1526  i = i + 1;
1527  }
1528 
1529  } // for vertices_per_cell
1530  k = k + 1;
1531  }
1532  local_triangulation.create_triangulation(vertices, cells, SubCellData());
1533  Assert(local_triangulation.n_active_cells() == uniform_cells.size(),
1534  ExcInternalError());
1535  local_triangulation.clear_user_flags();
1536  unsigned int index = 0;
1537  // Create a map between cells of class DoFHandler into class Triangulation
1538  std::map<typename Triangulation<Container::dimension,
1539  Container::space_dimension>::cell_iterator,
1540  typename Container::cell_iterator>
1541  patch_to_global_tria_map_tmp;
1542  for (typename Triangulation<Container::dimension,
1543  Container::space_dimension>::cell_iterator
1544  coarse_cell = local_triangulation.begin();
1545  coarse_cell != local_triangulation.end();
1546  ++coarse_cell, ++index)
1547  {
1548  patch_to_global_tria_map_tmp.insert(
1549  std::make_pair(coarse_cell, uniform_cells[index]));
1550  // To ensure that the cells with the same coordinates (here, we compare
1551  // their centers) are mapped into each other.
1552 
1553  Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1554  1e-15 * coarse_cell->diameter(),
1555  ExcInternalError());
1556  }
1557  bool refinement_necessary;
1558  // In this loop we start to do refinement on the above coarse triangulation
1559  // to reach to the same level of refinement as the patch cells are really on
1560  do
1561  {
1562  refinement_necessary = false;
1563  for (const auto &active_tria_cell :
1564  local_triangulation.active_cell_iterators())
1565  {
1566  if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1567  {
1568  active_tria_cell->set_refine_flag();
1569  refinement_necessary = true;
1570  }
1571  else
1572  for (unsigned int i = 0; i < patch.size(); ++i)
1573  {
1574  // Even though vertices may not be exactly the same, the
1575  // appropriate cells will match since == for TriAccessors
1576  // checks only cell level and index.
1577  if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1578  patch[i])
1579  {
1580  // adjust the cell vertices of the local_triangulation to
1581  // match cell vertices of the global triangulation
1582  for (const unsigned int v :
1583  active_tria_cell->vertex_indices())
1584  active_tria_cell->vertex(v) = patch[i]->vertex(v);
1585 
1586  Assert(active_tria_cell->center().distance(
1587  patch_to_global_tria_map_tmp[active_tria_cell]
1588  ->center()) <=
1589  1e-15 * active_tria_cell->diameter(),
1590  ExcInternalError());
1591 
1592  active_tria_cell->set_user_flag();
1593  break;
1594  }
1595  }
1596  }
1597 
1598  if (refinement_necessary)
1599  {
1600  local_triangulation.execute_coarsening_and_refinement();
1601 
1602  for (typename Triangulation<
1603  Container::dimension,
1604  Container::space_dimension>::cell_iterator cell =
1605  local_triangulation.begin();
1606  cell != local_triangulation.end();
1607  ++cell)
1608  {
1609  if (patch_to_global_tria_map_tmp.find(cell) !=
1610  patch_to_global_tria_map_tmp.end())
1611  {
1612  if (cell->has_children())
1613  {
1614  // Note: Since the cell got children, then it should not
1615  // be in the map anymore children may be added into the
1616  // map, instead
1617 
1618  // these children may not yet be in the map
1619  for (unsigned int c = 0; c < cell->n_children(); ++c)
1620  {
1621  if (patch_to_global_tria_map_tmp.find(cell->child(
1622  c)) == patch_to_global_tria_map_tmp.end())
1623  {
1624  patch_to_global_tria_map_tmp.insert(
1625  std::make_pair(
1626  cell->child(c),
1627  patch_to_global_tria_map_tmp[cell]->child(
1628  c)));
1629 
1630  // One might be tempted to assert that the cell
1631  // being added here has the same center as the
1632  // equivalent cell in the global triangulation,
1633  // but it may not be the case. For
1634  // triangulations that have been perturbed or
1635  // smoothed, the cell indices and levels may be
1636  // the same, but the vertex locations may not.
1637  // We adjust the vertices of the cells that have
1638  // no children (ie the active cells) to be
1639  // consistent with the global triangulation
1640  // later on and add assertions at that time
1641  // to guarantee the cells in the
1642  // local_triangulation are physically at the
1643  // same locations of the cells in the patch of
1644  // the global triangulation.
1645  }
1646  }
1647  // The parent cell whose children were added
1648  // into the map should be deleted from the map
1649  patch_to_global_tria_map_tmp.erase(cell);
1650  }
1651  }
1652  }
1653  }
1654  }
1655  while (refinement_necessary);
1656 
1657 
1658  // Last assertion check to make sure we have the right cells and centers
1659  // in the map, and hence the correct vertices of the triangulation
1660  for (typename Triangulation<Container::dimension,
1661  Container::space_dimension>::cell_iterator
1662  cell = local_triangulation.begin();
1663  cell != local_triangulation.end();
1664  ++cell)
1665  {
1666  if (cell->user_flag_set())
1667  {
1668  Assert(patch_to_global_tria_map_tmp.find(cell) !=
1669  patch_to_global_tria_map_tmp.end(),
1670  ExcInternalError());
1671 
1672  Assert(cell->center().distance(
1673  patch_to_global_tria_map_tmp[cell]->center()) <=
1674  1e-15 * cell->diameter(),
1675  ExcInternalError());
1676  }
1677  }
1678 
1679 
1680  typename std::map<
1681  typename Triangulation<Container::dimension,
1682  Container::space_dimension>::cell_iterator,
1683  typename Container::cell_iterator>::iterator
1684  map_tmp_it = patch_to_global_tria_map_tmp.begin(),
1685  map_tmp_end = patch_to_global_tria_map_tmp.end();
1686  // Now we just need to take the temporary map of pairs of type cell_iterator
1687  // "patch_to_global_tria_map_tmp" making pair of active_cell_iterators so
1688  // that filling out the final map "patch_to_global_tria_map"
1689  for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1690  patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1691  }
1692 
1693 
1694 
1695  template <int dim, int spacedim>
1696  std::map<
1698  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1700  {
1701  // This is the map from global_dof_index to
1702  // a set of cells on patch. We first map into
1703  // a set because it is very likely that we
1704  // will attempt to add a cell more than once
1705  // to a particular patch and we want to preserve
1706  // uniqueness of cell iterators. std::set does this
1707  // automatically for us. Later after it is all
1708  // constructed, we will copy to a map of vectors
1709  // since that is the preferred output for other
1710  // functions.
1711  std::map<types::global_dof_index,
1712  std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1713  dof_to_set_of_cells_map;
1714 
1715  std::vector<types::global_dof_index> local_dof_indices;
1716  std::vector<types::global_dof_index> local_face_dof_indices;
1717  std::vector<types::global_dof_index> local_line_dof_indices;
1718 
1719  // a place to save the dof_handler user flags and restore them later
1720  // to maintain const of dof_handler.
1721  std::vector<bool> user_flags;
1722 
1723 
1724  // in 3d, we need pointers from active lines to the
1725  // active parent lines, so we construct it as needed.
1726  std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1728  lines_to_parent_lines_map;
1729  if (dim == 3)
1730  {
1731  // save user flags as they will be modified and then later restored
1732  dof_handler.get_triangulation().save_user_flags(user_flags);
1733  const_cast<::Triangulation<dim, spacedim> &>(
1734  dof_handler.get_triangulation())
1735  .clear_user_flags();
1736 
1737 
1739  cell = dof_handler.begin_active(),
1740  endc = dof_handler.end();
1741  for (; cell != endc; ++cell)
1742  {
1743  // We only want lines that are locally_relevant
1744  // although it doesn't hurt to have lines that
1745  // are children of ghost cells since there are
1746  // few and we don't have to use them.
1747  if (cell->is_artificial() == false)
1748  {
1749  for (unsigned int l = 0; l < cell->n_lines(); ++l)
1750  if (cell->line(l)->has_children())
1751  for (unsigned int c = 0; c < cell->line(l)->n_children();
1752  ++c)
1753  {
1754  lines_to_parent_lines_map[cell->line(l)->child(c)] =
1755  cell->line(l);
1756  // set flags to know that child
1757  // line has an active parent.
1758  cell->line(l)->child(c)->set_user_flag();
1759  }
1760  }
1761  }
1762  }
1763 
1764 
1765  // We loop through all cells and add cell to the
1766  // map for the dofs that it immediately touches
1767  // and then account for all the other dofs of
1768  // which it is a part, mainly the ones that must
1769  // be added on account of adaptivity hanging node
1770  // constraints.
1772  cell = dof_handler.begin_active(),
1773  endc = dof_handler.end();
1774  for (; cell != endc; ++cell)
1775  {
1776  // Need to loop through all cells that could
1777  // be in the patch of dofs on locally_owned
1778  // cells including ghost cells
1779  if (cell->is_artificial() == false)
1780  {
1781  const unsigned int n_dofs_per_cell =
1782  cell->get_fe().n_dofs_per_cell();
1783  local_dof_indices.resize(n_dofs_per_cell);
1784 
1785  // Take care of adding cell pointer to each
1786  // dofs that exists on cell.
1787  cell->get_dof_indices(local_dof_indices);
1788  for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
1789  dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell);
1790 
1791  // In the case of the adjacent cell (over
1792  // faces or edges) being more refined, we
1793  // want to add all of the children to the
1794  // patch since the support function at that
1795  // dof could be non-zero along that entire
1796  // face (or line).
1797 
1798  // Take care of dofs on neighbor faces
1799  for (const unsigned int f : cell->face_indices())
1800  {
1801  if (cell->face(f)->has_children())
1802  {
1803  for (unsigned int c = 0; c < cell->face(f)->n_children();
1804  ++c)
1805  {
1806  // Add cell to dofs of all subfaces
1807  //
1808  // *-------------------*----------*---------*
1809  // | | add cell | |
1810  // | |<- to dofs| |
1811  // | |of subface| |
1812  // | cell *----------*---------*
1813  // | | add cell | |
1814  // | |<- to dofs| |
1815  // | |of subface| |
1816  // *-------------------*----------*---------*
1817  //
1818  Assert(cell->face(f)->child(c)->has_children() == false,
1819  ExcInternalError());
1820 
1821  const unsigned int n_dofs_per_face =
1822  cell->get_fe().n_dofs_per_face(f, c);
1823  local_face_dof_indices.resize(n_dofs_per_face);
1824 
1825  cell->face(f)->child(c)->get_dof_indices(
1826  local_face_dof_indices);
1827  for (unsigned int i = 0; i < n_dofs_per_face; ++i)
1828  dof_to_set_of_cells_map[local_face_dof_indices[i]]
1829  .insert(cell);
1830  }
1831  }
1832  else if ((cell->face(f)->at_boundary() == false) &&
1833  (cell->neighbor_is_coarser(f)))
1834  {
1835  // Add cell to dofs of parent face and all
1836  // child faces of parent face
1837  //
1838  // *-------------------*----------*---------*
1839  // | | | |
1840  // | | cell | |
1841  // | add cell | | |
1842  // | to dofs -> *----------*---------*
1843  // | of parent | add cell | |
1844  // | face |<- to dofs| |
1845  // | |of subface| |
1846  // *-------------------*----------*---------*
1847  //
1848 
1849  // Add cell to all dofs of parent face
1850  std::pair<unsigned int, unsigned int>
1851  neighbor_face_no_subface_no =
1852  cell->neighbor_of_coarser_neighbor(f);
1853  unsigned int face_no = neighbor_face_no_subface_no.first;
1854  unsigned int subface = neighbor_face_no_subface_no.second;
1855 
1856  const unsigned int n_dofs_per_face =
1857  cell->get_fe().n_dofs_per_face(face_no);
1858  local_face_dof_indices.resize(n_dofs_per_face);
1859 
1860  cell->neighbor(f)->face(face_no)->get_dof_indices(
1861  local_face_dof_indices);
1862  for (unsigned int i = 0; i < n_dofs_per_face; ++i)
1863  dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(
1864  cell);
1865 
1866  // Add cell to all dofs of children of
1867  // parent face
1868  for (unsigned int c = 0;
1869  c < cell->neighbor(f)->face(face_no)->n_children();
1870  ++c)
1871  {
1872  if (c != subface) // don't repeat work on dofs of
1873  // original cell
1874  {
1875  const unsigned int n_dofs_per_face =
1876  cell->get_fe().n_dofs_per_face(face_no, c);
1877  local_face_dof_indices.resize(n_dofs_per_face);
1878 
1879  Assert(cell->neighbor(f)
1880  ->face(face_no)
1881  ->child(c)
1882  ->has_children() == false,
1883  ExcInternalError());
1884  cell->neighbor(f)
1885  ->face(face_no)
1886  ->child(c)
1887  ->get_dof_indices(local_face_dof_indices);
1888  for (unsigned int i = 0; i < n_dofs_per_face; ++i)
1889  dof_to_set_of_cells_map[local_face_dof_indices[i]]
1890  .insert(cell);
1891  }
1892  }
1893  }
1894  }
1895 
1896 
1897  // If 3d, take care of dofs on lines in the
1898  // same pattern as faces above. That is, if
1899  // a cell's line has children, distribute
1900  // cell to dofs of children of line, and
1901  // if cell's line has an active parent, then
1902  // distribute cell to dofs on parent line
1903  // and dofs on all children of parent line.
1904  if (dim == 3)
1905  {
1906  for (unsigned int l = 0; l < cell->n_lines(); ++l)
1907  {
1908  if (cell->line(l)->has_children())
1909  {
1910  for (unsigned int c = 0;
1911  c < cell->line(l)->n_children();
1912  ++c)
1913  {
1914  Assert(cell->line(l)->child(c)->has_children() ==
1915  false,
1916  ExcInternalError());
1917 
1918  // dofs_per_line returns number of dofs
1919  // on line not including the vertices of the line.
1920  const unsigned int n_dofs_per_line =
1921  2 * cell->get_fe().n_dofs_per_vertex() +
1922  cell->get_fe().n_dofs_per_line();
1923  local_line_dof_indices.resize(n_dofs_per_line);
1924 
1925  cell->line(l)->child(c)->get_dof_indices(
1926  local_line_dof_indices);
1927  for (unsigned int i = 0; i < n_dofs_per_line; ++i)
1928  dof_to_set_of_cells_map[local_line_dof_indices[i]]
1929  .insert(cell);
1930  }
1931  }
1932  // user flag was set above to denote that
1933  // an active parent line exists so add
1934  // cell to dofs of parent and all it's
1935  // children
1936  else if (cell->line(l)->user_flag_set() == true)
1937  {
1939  parent_line =
1940  lines_to_parent_lines_map[cell->line(l)];
1941  Assert(parent_line->has_children(), ExcInternalError());
1942 
1943  // dofs_per_line returns number of dofs
1944  // on line not including the vertices of the line.
1945  const unsigned int n_dofs_per_line =
1946  2 * cell->get_fe().n_dofs_per_vertex() +
1947  cell->get_fe().n_dofs_per_line();
1948  local_line_dof_indices.resize(n_dofs_per_line);
1949 
1950  parent_line->get_dof_indices(local_line_dof_indices);
1951  for (unsigned int i = 0; i < n_dofs_per_line; ++i)
1952  dof_to_set_of_cells_map[local_line_dof_indices[i]]
1953  .insert(cell);
1954 
1955  for (unsigned int c = 0; c < parent_line->n_children();
1956  ++c)
1957  {
1958  Assert(parent_line->child(c)->has_children() ==
1959  false,
1960  ExcInternalError());
1961 
1962  const unsigned int n_dofs_per_line =
1963  2 * cell->get_fe().n_dofs_per_vertex() +
1964  cell->get_fe().n_dofs_per_line();
1965  local_line_dof_indices.resize(n_dofs_per_line);
1966 
1967  parent_line->child(c)->get_dof_indices(
1968  local_line_dof_indices);
1969  for (unsigned int i = 0; i < n_dofs_per_line; ++i)
1970  dof_to_set_of_cells_map[local_line_dof_indices[i]]
1971  .insert(cell);
1972  }
1973  }
1974  } // for lines l
1975  } // if dim == 3
1976  } // if cell->is_locally_owned()
1977  } // for cells
1978 
1979 
1980  if (dim == 3)
1981  {
1982  // finally, restore user flags that were changed above
1983  // to when we constructed the pointers to parent of lines
1984  // Since dof_handler is const, we must leave it unchanged.
1985  const_cast<::Triangulation<dim, spacedim> &>(
1986  dof_handler.get_triangulation())
1987  .load_user_flags(user_flags);
1988  }
1989 
1990  // Finally, we copy map of sets to
1991  // map of vectors using the std::vector::assign() function
1992  std::map<
1994  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1995  dof_to_cell_patches;
1996 
1997  typename std::map<
1999  std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
2000  iterator it = dof_to_set_of_cells_map.begin(),
2001  it_end = dof_to_set_of_cells_map.end();
2002  for (; it != it_end; ++it)
2003  dof_to_cell_patches[it->first].assign(it->second.begin(),
2004  it->second.end());
2005 
2006  return dof_to_cell_patches;
2007  }
2008 
2009  /*
2010  * Internally used in collect_periodic_faces
2011  */
2012  template <typename CellIterator>
2013  void
2015  std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2016  std::set<std::pair<typename identity<CellIterator>::type, unsigned int>>
2017  & pairs2,
2018  const int direction,
2019  std::vector<PeriodicFacePair<CellIterator>> &matched_pairs,
2020  const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2021  & offset,
2022  const FullMatrix<double> &matrix)
2023  {
2024  static const int space_dim = CellIterator::AccessorType::space_dimension;
2025  (void)space_dim;
2026  AssertIndexRange(direction, space_dim);
2027 
2028 #ifdef DEBUG
2029  {
2030  constexpr int dim = CellIterator::AccessorType::dimension;
2031  constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2032  // For parallel::fullydistributed::Triangulation there might be unmatched
2033  // faces on periodic boundaries on the coarse grid, which results that
2034  // this assert is not fulfilled (not a bug!). See also the discussion in
2035  // the method collect_periodic_faces.
2036  if (!(((pairs1.size() > 0) &&
2037  (dynamic_cast<const parallel::fullydistributed::
2038  Triangulation<dim, spacedim> *>(
2039  &pairs1.begin()->first->get_triangulation()) != nullptr)) ||
2040  ((pairs2.size() > 0) &&
2041  (dynamic_cast<
2043  *>(&pairs2.begin()->first->get_triangulation()) != nullptr))))
2044  Assert(pairs1.size() == pairs2.size(),
2045  ExcMessage("Unmatched faces on periodic boundaries"));
2046  }
2047 #endif
2048 
2049  unsigned int n_matches = 0;
2050 
2051  // Match with a complexity of O(n^2). This could be improved...
2052  std::bitset<3> orientation;
2053  using PairIterator =
2054  typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2055  for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2056  {
2057  for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2058  {
2059  const CellIterator cell1 = it1->first;
2060  const CellIterator cell2 = it2->first;
2061  const unsigned int face_idx1 = it1->second;
2062  const unsigned int face_idx2 = it2->second;
2063  if (GridTools::orthogonal_equality(orientation,
2064  cell1->face(face_idx1),
2065  cell2->face(face_idx2),
2066  direction,
2067  offset,
2068  matrix))
2069  {
2070  // We have a match, so insert the matching pairs and
2071  // remove the matched cell in pairs2 to speed up the
2072  // matching:
2073  const PeriodicFacePair<CellIterator> matched_face = {
2074  {cell1, cell2}, {face_idx1, face_idx2}, orientation, matrix};
2075  matched_pairs.push_back(matched_face);
2076  pairs2.erase(it2);
2077  ++n_matches;
2078  break;
2079  }
2080  }
2081  }
2082 
2083  // Assure that all faces are matched if not
2084  // parallel::fullydistributed::Triangulation is used. This is related to the
2085  // fact that the faces might not be successfully matched on the coarse grid
2086  // (not a bug!). See also the comment above and in the
2087  // method collect_periodic_faces.
2088  {
2089  constexpr int dim = CellIterator::AccessorType::dimension;
2090  constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2091  if (!(((pairs1.size() > 0) &&
2092  (dynamic_cast<const parallel::fullydistributed::
2093  Triangulation<dim, spacedim> *>(
2094  &pairs1.begin()->first->get_triangulation()) != nullptr)) ||
2095  ((pairs2.size() > 0) &&
2096  (dynamic_cast<
2098  *>(&pairs2.begin()->first->get_triangulation()) != nullptr))))
2099  AssertThrow(n_matches == pairs1.size() && pairs2.size() == 0,
2100  ExcMessage("Unmatched faces on periodic boundaries"));
2101  }
2102  }
2103 
2104 
2105 
2106  template <typename MeshType>
2107  void
2109  const MeshType & mesh,
2110  const types::boundary_id b_id,
2111  const int direction,
2113  & matched_pairs,
2115  const FullMatrix<double> & matrix)
2116  {
2117  static const int dim = MeshType::dimension;
2118  static const int space_dim = MeshType::space_dimension;
2119  (void)dim;
2120  (void)space_dim;
2121  AssertIndexRange(direction, space_dim);
2122 
2123  Assert(dim == space_dim, ExcNotImplemented());
2124 
2125  // Loop over all cells on the highest level and collect all boundary
2126  // faces 2*direction and 2*direction*1:
2127 
2128  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2129  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2130 
2131  for (typename MeshType::cell_iterator cell = mesh.begin(0);
2132  cell != mesh.end(0);
2133  ++cell)
2134  {
2135  const typename MeshType::face_iterator face_1 =
2136  cell->face(2 * direction);
2137  const typename MeshType::face_iterator face_2 =
2138  cell->face(2 * direction + 1);
2139 
2140  if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2141  {
2142  const std::pair<typename MeshType::cell_iterator, unsigned int>
2143  pair1 = std::make_pair(cell, 2 * direction);
2144  pairs1.insert(pair1);
2145  }
2146 
2147  if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2148  {
2149  const std::pair<typename MeshType::cell_iterator, unsigned int>
2150  pair2 = std::make_pair(cell, 2 * direction + 1);
2151  pairs2.insert(pair2);
2152  }
2153  }
2154 
2155  Assert(pairs1.size() == pairs2.size(),
2156  ExcMessage("Unmatched faces on periodic boundaries"));
2157 
2158  Assert(pairs1.size() > 0,
2159  ExcMessage("No new periodic face pairs have been found. "
2160  "Are you sure that you've selected the correct boundary "
2161  "id's and that the coarsest level mesh is colorized?"));
2162 
2163 #ifdef DEBUG
2164  const unsigned int size_old = matched_pairs.size();
2165 #endif
2166 
2167  // and call match_periodic_face_pairs that does the actual matching:
2169  pairs1, pairs2, direction, matched_pairs, offset, matrix);
2170 
2171 #ifdef DEBUG
2172  // check for standard orientation
2173  const unsigned int size_new = matched_pairs.size();
2174  for (unsigned int i = size_old; i < size_new; ++i)
2175  {
2176  Assert(matched_pairs[i].orientation == 1,
2177  ExcMessage(
2178  "Found a face match with non standard orientation. "
2179  "This function is only suitable for meshes with cells "
2180  "in default orientation"));
2181  }
2182 #endif
2183  }
2184 
2185 
2186 
2187  template <typename MeshType>
2188  void
2190  const MeshType & mesh,
2191  const types::boundary_id b_id1,
2192  const types::boundary_id b_id2,
2193  const int direction,
2195  & matched_pairs,
2197  const FullMatrix<double> & matrix)
2198  {
2199  static const int dim = MeshType::dimension;
2200  static const int space_dim = MeshType::space_dimension;
2201  (void)dim;
2202  (void)space_dim;
2203  AssertIndexRange(direction, space_dim);
2204 
2205  // Loop over all cells on the highest level and collect all boundary
2206  // faces belonging to b_id1 and b_id2:
2207 
2208  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2209  std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2210 
2211  for (typename MeshType::cell_iterator cell = mesh.begin(0);
2212  cell != mesh.end(0);
2213  ++cell)
2214  {
2215  for (unsigned int i : cell->face_indices())
2216  {
2217  const typename MeshType::face_iterator face = cell->face(i);
2218  if (face->at_boundary() && face->boundary_id() == b_id1)
2219  {
2220  const std::pair<typename MeshType::cell_iterator, unsigned int>
2221  pair1 = std::make_pair(cell, i);
2222  pairs1.insert(pair1);
2223  }
2224 
2225  if (face->at_boundary() && face->boundary_id() == b_id2)
2226  {
2227  const std::pair<typename MeshType::cell_iterator, unsigned int>
2228  pair2 = std::make_pair(cell, i);
2229  pairs2.insert(pair2);
2230  }
2231  }
2232  }
2233 
2234  // Assure that all faces are matched on the coare grid. This requirement
2235  // can only fulfilled if a process owns the complete coarse grid. This is
2236  // not the case for a parallel::fullydistributed::Triangulation, i.e. this
2237  // requirement has not to be met neither (consider faces on the outside of a
2238  // ghost cell that are periodic but for which the ghost neighbor doesn't
2239  // exist).
2240  if (!(((pairs1.size() > 0) &&
2241  (dynamic_cast<
2243  *>(&pairs1.begin()->first->get_triangulation()) != nullptr)) ||
2244  ((pairs2.size() > 0) &&
2245  (dynamic_cast<
2247  *>(&pairs2.begin()->first->get_triangulation()) != nullptr))))
2248  Assert(pairs1.size() == pairs2.size(),
2249  ExcMessage("Unmatched faces on periodic boundaries"));
2250 
2251  Assert(
2252  (pairs1.size() > 0 ||
2253  (dynamic_cast<
2255  &mesh.begin()->get_triangulation()) != nullptr)),
2256  ExcMessage("No new periodic face pairs have been found. "
2257  "Are you sure that you've selected the correct boundary "
2258  "id's and that the coarsest level mesh is colorized?"));
2259 
2260  // and call match_periodic_face_pairs that does the actual matching:
2262  pairs1, pairs2, direction, matched_pairs, offset, matrix);
2263  }
2264 
2265 
2266 
2267  /*
2268  * Internally used in orthogonal_equality
2269  *
2270  * An orthogonal equality test for points:
2271  *
2272  * point1 and point2 are considered equal, if
2273  * matrix.point1 + offset - point2
2274  * is parallel to the unit vector in <direction>
2275  */
2276  template <int spacedim>
2277  inline bool
2279  const Point<spacedim> & point2,
2280  const int direction,
2281  const Tensor<1, spacedim> &offset,
2282  const FullMatrix<double> & matrix)
2283  {
2284  AssertIndexRange(direction, spacedim);
2285 
2286  Assert(matrix.m() == matrix.n(), ExcInternalError());
2287 
2288  Point<spacedim> distance;
2289 
2290  if (matrix.m() == spacedim)
2291  for (int i = 0; i < spacedim; ++i)
2292  for (int j = 0; j < spacedim; ++j)
2293  distance(i) += matrix(i, j) * point1(j);
2294  else
2295  distance = point1;
2296 
2297  distance += offset - point2;
2298 
2299  for (int i = 0; i < spacedim; ++i)
2300  {
2301  // Only compare coordinate-components != direction:
2302  if (i == direction)
2303  continue;
2304 
2305  if (std::abs(distance(i)) > 1.e-10)
2306  return false;
2307  }
2308 
2309  return true;
2310  }
2311 
2312 
2313  /*
2314  * Internally used in orthogonal_equality
2315  *
2316  * A lookup table to transform vertex matchings to orientation flags of
2317  * the form (face_orientation, face_flip, face_rotation)
2318  *
2319  * See the comment on the next function as well as the detailed
2320  * documentation of make_periodicity_constraints and
2321  * collect_periodic_faces for details
2322  */
2323  template <int dim>
2325  {};
2326 
2327  template <>
2329  {
2330  using MATCH_T =
2331  std::array<unsigned int, GeometryInfo<1>::vertices_per_face>;
2332  static inline std::bitset<3>
2333  lookup(const MATCH_T &)
2334  {
2335  // The 1D case is trivial
2336  return 1; // [true ,false,false]
2337  }
2338  };
2339 
2340  template <>
2342  {
2343  using MATCH_T =
2344  std::array<unsigned int, GeometryInfo<2>::vertices_per_face>;
2345  static inline std::bitset<3>
2346  lookup(const MATCH_T &matching)
2347  {
2348  // In 2D matching faces (=lines) results in two cases: Either
2349  // they are aligned or flipped. We store this "line_flip"
2350  // property somewhat sloppy as "face_flip"
2351  // (always: face_orientation = true, face_rotation = false)
2352 
2353  static const MATCH_T m_tff = {{0, 1}};
2354  if (matching == m_tff)
2355  return 1; // [true ,false,false]
2356  static const MATCH_T m_ttf = {{1, 0}};
2357  if (matching == m_ttf)
2358  return 3; // [true ,true ,false]
2359  Assert(false, ExcInternalError());
2360  // what follows is dead code, but it avoids warnings about the lack
2361  // of a return value
2362  return 0;
2363  }
2364  };
2365 
2366  template <>
2368  {
2369  using MATCH_T =
2370  std::array<unsigned int, GeometryInfo<3>::vertices_per_face>;
2371 
2372  static inline std::bitset<3>
2373  lookup(const MATCH_T &matching)
2374  {
2375  // The full fledged 3D case. *Yay*
2376  // See the documentation in include/deal.II/base/geometry_info.h
2377  // as well as the actual implementation in source/grid/tria.cc
2378  // for more details...
2379 
2380  static const MATCH_T m_tff = {{0, 1, 2, 3}};
2381  if (matching == m_tff)
2382  return 1; // [true ,false,false]
2383  static const MATCH_T m_tft = {{1, 3, 0, 2}};
2384  if (matching == m_tft)
2385  return 5; // [true ,false,true ]
2386  static const MATCH_T m_ttf = {{3, 2, 1, 0}};
2387  if (matching == m_ttf)
2388  return 3; // [true ,true ,false]
2389  static const MATCH_T m_ttt = {{2, 0, 3, 1}};
2390  if (matching == m_ttt)
2391  return 7; // [true ,true ,true ]
2392  static const MATCH_T m_fff = {{0, 2, 1, 3}};
2393  if (matching == m_fff)
2394  return 0; // [false,false,false]
2395  static const MATCH_T m_fft = {{2, 3, 0, 1}};
2396  if (matching == m_fft)
2397  return 4; // [false,false,true ]
2398  static const MATCH_T m_ftf = {{3, 1, 2, 0}};
2399  if (matching == m_ftf)
2400  return 2; // [false,true ,false]
2401  static const MATCH_T m_ftt = {{1, 0, 3, 2}};
2402  if (matching == m_ftt)
2403  return 6; // [false,true ,true ]
2404  Assert(false, ExcInternalError());
2405  // what follows is dead code, but it avoids warnings about the lack
2406  // of a return value
2407  return 0;
2408  }
2409  };
2410 
2411 
2412 
2413  template <typename FaceIterator>
2414  inline bool orthogonal_equality(
2415  std::bitset<3> & orientation,
2416  const FaceIterator & face1,
2417  const FaceIterator & face2,
2418  const int direction,
2420  const FullMatrix<double> & matrix)
2421  {
2422  Assert(matrix.m() == matrix.n(),
2423  ExcMessage("The supplied matrix must be a square matrix"));
2424 
2425  static const int dim = FaceIterator::AccessorType::dimension;
2426 
2427  // Do a full matching of the face vertices:
2428 
2429  std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
2430 
2431  AssertDimension(face1->n_vertices(), face2->n_vertices());
2432 
2433  std::set<unsigned int> face2_vertices;
2434  for (unsigned int i = 0; i < face1->n_vertices(); ++i)
2435  face2_vertices.insert(i);
2436 
2437  for (unsigned int i = 0; i < face1->n_vertices(); ++i)
2438  {
2439  for (std::set<unsigned int>::iterator it = face2_vertices.begin();
2440  it != face2_vertices.end();
2441  ++it)
2442  {
2443  if (orthogonal_equality(face1->vertex(i),
2444  face2->vertex(*it),
2445  direction,
2446  offset,
2447  matrix))
2448  {
2449  matching[i] = *it;
2450  face2_vertices.erase(it);
2451  break; // jump out of the innermost loop
2452  }
2453  }
2454  }
2455 
2456  // And finally, a lookup to determine the ordering bitmask:
2457  if (face2_vertices.empty())
2458  orientation = OrientationLookupTable<dim>::lookup(matching);
2459 
2460  return face2_vertices.empty();
2461  }
2462 
2463 
2464 
2465  template <typename FaceIterator>
2466  inline bool
2468  const FaceIterator & face1,
2469  const FaceIterator & face2,
2470  const int direction,
2472  const FullMatrix<double> & matrix)
2473  {
2474  // Call the function above with a dummy orientation array
2475  std::bitset<3> dummy;
2476  return orthogonal_equality(dummy, face1, face2, direction, offset, matrix);
2477  }
2478 
2479 
2480 
2481 } // namespace GridTools
2482 
2483 
2484 #include "grid_tools_dof_handlers.inst"
2485 
2486 
size_type m() const
unsigned int n_active_cells() const
Definition: tria.cc:13821
static std::bitset< 3 > lookup(const MATCH_T &matching)
const Triangulation< dim, spacedim > & get_triangulation() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1656
void clear_user_flags()
Definition: tria.cc:12295
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:13813
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:394
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:13331
cell_iterator end() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1721
std::array< unsigned int, GeometryInfo< 3 >::vertices_per_face > MATCH_T
#define AssertThrow(cond, exc)
Definition: exceptions.h:1588
Point< 2 > second
Definition: grid_out.cc:4588
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)
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:13114
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:6187
static std::bitset< 3 > lookup(const MATCH_T &)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
AdjacentCell adjacent_cells[2]
Definition: grid_tools.cc:1176
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:13225
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:14488
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:11622
#define Assert(cond, exc)
Definition: exceptions.h:1478
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:303
::DoFHandler< dim, spacedim > DoFHandler
Definition: dof_handler.h:124
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1256
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:399
unsigned int level
Definition: grid_out.cc:4590
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)
const hp::FECollection< dim, spacedim > & get_fe_collection() const
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13771
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:13764
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Point< 2 > first
Definition: grid_out.cc:4587
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:14365
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:398
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:2568
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()
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)
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:6206
virtual void clear()
Definition: tria.cc:11285
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:14441