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