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