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