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