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