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