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