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