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