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