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