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