deal.II version GIT relicensing-2848-g5241f990fb 2025-03-16 19:30:00+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_topology.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
19
22
23#include <boost/container/small_vector.hpp>
24
25#include <algorithm>
26#include <map>
27#include <numeric>
28#include <set>
29#include <vector>
30
32
33namespace GridTools
34{
35 // Generic functions for appending face data in 2d or 3d. TODO: we can
36 // remove these once we have 'if constexpr'.
37 namespace internal
38 {
39 inline void
40 append_face_data(const CellData<1> &face_data, SubCellData &subcell_data)
41 {
42 subcell_data.boundary_lines.push_back(face_data);
43 }
44
45
46
47 inline void
48 append_face_data(const CellData<2> &face_data, SubCellData &subcell_data)
49 {
50 subcell_data.boundary_quads.push_back(face_data);
51 }
52
53
54
55 // Lexical comparison for sorting CellData objects.
56 template <int structdim>
58 {
59 bool
61 const CellData<structdim> &b) const
62 {
63 // Check vertices:
64 if (std::lexicographical_compare(std::begin(a.vertices),
65 std::end(a.vertices),
66 std::begin(b.vertices),
67 std::end(b.vertices)))
68 return true;
69 // it should never be necessary to check the material or manifold
70 // ids as a 'tiebreaker' (since they must be equal if the vertex
71 // indices are equal). Assert it anyway:
72 if constexpr (running_in_debug_mode())
73 {
74 if (std::equal(std::begin(a.vertices),
75 std::end(a.vertices),
76 std::begin(b.vertices)))
77 {
78 Assert(a.material_id == b.material_id &&
79 a.manifold_id == b.manifold_id,
81 "Two CellData objects with equal vertices must "
82 "have the same material/boundary ids and manifold "
83 "ids."));
84 }
85 }
86 return false;
87 }
88 };
89
90
100 template <int dim>
102 {
103 public:
107 template <typename FaceIteratorType>
108 void
109 insert_face_data(const FaceIteratorType &face)
110 {
111 CellData<dim - 1> face_cell_data(face->n_vertices());
112 for (unsigned int vertex_n = 0; vertex_n < face->n_vertices();
113 ++vertex_n)
114 face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
115 face_cell_data.boundary_id = face->boundary_id();
116 face_cell_data.manifold_id = face->manifold_id();
117
118 face_data.insert(std::move(face_cell_data));
119 }
120
126 {
127 SubCellData subcell_data;
128
129 for (const CellData<dim - 1> &face_cell_data : face_data)
130 internal::append_face_data(face_cell_data, subcell_data);
131 return subcell_data;
132 }
133
134
135 private:
136 std::set<CellData<dim - 1>, internal::CellDataComparator<dim - 1>>
138 };
139
140
141 // Do nothing for dim=1:
142 template <>
144 {
145 public:
146 template <typename FaceIteratorType>
147 void
148 insert_face_data(const FaceIteratorType &)
149 {}
150
153 {
154 return SubCellData();
155 }
156 };
157 } // namespace internal
158
159
160
161 template <int dim, int spacedim>
162 std::
163 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
165 {
166 Assert(tria.n_levels() >= 1,
167 ExcMessage("The input triangulation must be non-empty."));
168
169 std::vector<Point<spacedim>> vertices = tria.get_vertices();
170 std::vector<CellData<dim>> cells;
171
173 std::set<CellData<1>, internal::CellDataComparator<1>>
174 line_data; // only used in 3d
175
176 for (const auto &cell : tria.cell_iterators_on_level(0))
177 {
178 // Save cell data
179 CellData<dim> cell_data(cell->n_vertices());
180 for (const unsigned int cell_vertex_n : cell->vertex_indices())
181 {
182 Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
184 cell_data.vertices[cell_vertex_n] =
185 cell->vertex_index(cell_vertex_n);
186 }
187 cell_data.material_id = cell->material_id();
188 cell_data.manifold_id = cell->manifold_id();
189 cells.emplace_back(std::move(cell_data));
190
191 // Save face data
192 if (dim > 1)
193 {
194 for (const unsigned int face_n : cell->face_indices())
195 // We don't need to insert anything if we have default values
196 {
197 const auto face = cell->face(face_n);
198 if (face->boundary_id() != numbers::internal_face_boundary_id ||
199 face->manifold_id() != numbers::flat_manifold_id)
200 face_data.insert_face_data(face);
201 }
202 }
203 // Save line data
204 if (dim == 3)
205 {
206 for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
207 {
208 const auto line = cell->line(line_n);
209 // We don't need to insert anything if we have default values
210 if (line->boundary_id() != numbers::internal_face_boundary_id ||
211 line->manifold_id() != numbers::flat_manifold_id)
212 {
213 CellData<1> line_cell_data(line->n_vertices());
214 for (const unsigned int vertex_n : line->vertex_indices())
215 line_cell_data.vertices[vertex_n] =
216 line->vertex_index(vertex_n);
217 line_cell_data.boundary_id = line->boundary_id();
218 line_cell_data.manifold_id = line->manifold_id();
219 line_data.insert(std::move(line_cell_data));
220 }
221 }
222 }
223 }
224
225 SubCellData subcell_data = face_data.get();
226
227 if (dim == 3)
228 for (const CellData<1> &face_line_data : line_data)
229 subcell_data.boundary_lines.push_back(face_line_data);
230
231 // We end up with a 'vertices' array that uses some of the entries,
232 // but not all -- specifically, all vertices referenced by level-0
233 // cells. We can compress the array:
234 GridTools::delete_unused_vertices(vertices, cells, subcell_data);
235
236 return std::tuple<std::vector<Point<spacedim>>,
237 std::vector<CellData<dim>>,
238 SubCellData>(std::move(vertices),
239 std::move(cells),
240 std::move(subcell_data));
241 }
242
243
244
245 template <int dim, int spacedim>
246 void
248 std::vector<CellData<dim>> &cells,
249 SubCellData &subcelldata)
250 {
251 Assert(
252 subcelldata.check_consistency(dim),
254 "Invalid SubCellData supplied according to ::check_consistency(). "
255 "This is caused by data containing objects for the wrong dimension."));
256
257 // first check which vertices are actually used
258 std::vector<bool> vertex_used(vertices.size(), false);
259 for (unsigned int c = 0; c < cells.size(); ++c)
260 for (unsigned int v = 0; v < cells[c].vertices.size(); ++v)
261 {
262 Assert(cells[c].vertices[v] < vertices.size(),
263 ExcMessage("Invalid vertex index encountered! cells[" +
264 Utilities::int_to_string(c) + "].vertices[" +
265 Utilities::int_to_string(v) + "]=" +
266 Utilities::int_to_string(cells[c].vertices[v]) +
267 " is invalid, because only " +
268 Utilities::int_to_string(vertices.size()) +
269 " vertices were supplied."));
270 vertex_used[cells[c].vertices[v]] = true;
271 }
272
273
274 // then renumber the vertices that are actually used in the same order as
275 // they were beforehand
276 const unsigned int invalid_vertex = numbers::invalid_unsigned_int;
277 std::vector<unsigned int> new_vertex_numbers(vertices.size(),
278 invalid_vertex);
279 unsigned int next_free_number = 0;
280 for (unsigned int i = 0; i < vertices.size(); ++i)
281 if (vertex_used[i] == true)
282 {
283 new_vertex_numbers[i] = next_free_number;
284 ++next_free_number;
285 }
286
287 // next replace old vertex numbers by the new ones
288 for (unsigned int c = 0; c < cells.size(); ++c)
289 for (auto &v : cells[c].vertices)
290 v = new_vertex_numbers[v];
291
292 // same for boundary data
293 for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT
294 ++c)
295 for (unsigned int v = 0;
296 v < subcelldata.boundary_lines[c].vertices.size();
297 ++v)
298 {
299 Assert(subcelldata.boundary_lines[c].vertices[v] <
300 new_vertex_numbers.size(),
302 "Invalid vertex index in subcelldata.boundary_lines. "
303 "subcelldata.boundary_lines[" +
304 Utilities::int_to_string(c) + "].vertices[" +
305 Utilities::int_to_string(v) + "]=" +
307 subcelldata.boundary_lines[c].vertices[v]) +
308 " is invalid, because only " +
309 Utilities::int_to_string(vertices.size()) +
310 " vertices were supplied."));
311 subcelldata.boundary_lines[c].vertices[v] =
312 new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
313 }
314
315 for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT
316 ++c)
317 for (unsigned int v = 0;
318 v < subcelldata.boundary_quads[c].vertices.size();
319 ++v)
320 {
321 Assert(subcelldata.boundary_quads[c].vertices[v] <
322 new_vertex_numbers.size(),
324 "Invalid vertex index in subcelldata.boundary_quads. "
325 "subcelldata.boundary_quads[" +
326 Utilities::int_to_string(c) + "].vertices[" +
327 Utilities::int_to_string(v) + "]=" +
329 subcelldata.boundary_quads[c].vertices[v]) +
330 " is invalid, because only " +
331 Utilities::int_to_string(vertices.size()) +
332 " vertices were supplied."));
333
334 subcelldata.boundary_quads[c].vertices[v] =
335 new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
336 }
337
338 // finally copy over the vertices which we really need to a new array and
339 // replace the old one by the new one
340 std::vector<Point<spacedim>> tmp;
341 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(), true));
342 for (unsigned int v = 0; v < vertices.size(); ++v)
343 if (vertex_used[v] == true)
344 tmp.push_back(vertices[v]);
345 swap(vertices, tmp);
346 }
347
348
349
350 template <int dim, int spacedim>
351 void
353 std::vector<CellData<dim>> &cells,
354 SubCellData &subcelldata,
355 std::vector<unsigned int> &considered_vertices,
356 const double tol)
357 {
358 if (tol == 0.0)
359 return; // nothing to do per definition
360
361 AssertIndexRange(2, vertices.size());
362 std::vector<unsigned int> new_vertex_numbers(vertices.size());
363 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
364
365 // if the considered_vertices vector is empty, consider all vertices
366 if (considered_vertices.empty())
367 considered_vertices = new_vertex_numbers;
368 Assert(considered_vertices.size() <= vertices.size(), ExcInternalError());
369
370 // The algorithm below improves upon the naive O(n^2) algorithm by first
371 // sorting vertices by their value in one component and then only
372 // comparing vertices for equality which are nearly equal in that
373 // component. For example, if @p vertices form a cube, then we will only
374 // compare points that have the same x coordinate when we try to find
375 // duplicated vertices.
376
377 // Start by finding the longest coordinate direction. This minimizes the
378 // number of points that need to be compared against each-other in a
379 // single set for typical geometries.
380 const BoundingBox<spacedim> bbox(vertices);
381
382 unsigned int longest_coordinate_direction = 0;
383 double longest_coordinate_length = bbox.side_length(0);
384 for (unsigned int d = 1; d < spacedim; ++d)
385 {
386 const double coordinate_length = bbox.side_length(d);
387 if (longest_coordinate_length < coordinate_length)
388 {
389 longest_coordinate_length = coordinate_length;
390 longest_coordinate_direction = d;
391 }
392 }
393
394 // Sort vertices (while preserving their vertex numbers) along that
395 // coordinate direction:
396 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
397 sorted_vertices.reserve(vertices.size());
398 for (const unsigned int vertex_n : considered_vertices)
399 {
400 AssertIndexRange(vertex_n, vertices.size());
401 sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]);
402 }
403 std::sort(sorted_vertices.begin(),
404 sorted_vertices.end(),
405 [&](const std::pair<unsigned int, Point<spacedim>> &a,
406 const std::pair<unsigned int, Point<spacedim>> &b) {
407 return a.second[longest_coordinate_direction] <
408 b.second[longest_coordinate_direction];
409 });
410
411 auto within_tolerance = [=](const Point<spacedim> &a,
412 const Point<spacedim> &b) {
413 for (unsigned int d = 0; d < spacedim; ++d)
414 if (std::abs(a[d] - b[d]) > tol)
415 return false;
416 return true;
417 };
418
419 // Find a range of numbers that have the same component in the longest
420 // coordinate direction:
421 auto range_start = sorted_vertices.begin();
422 while (range_start != sorted_vertices.end())
423 {
424 auto range_end = range_start + 1;
425 while (range_end != sorted_vertices.end() &&
426 std::abs(range_end->second[longest_coordinate_direction] -
427 range_start->second[longest_coordinate_direction]) <
428 tol)
429 ++range_end;
430
431 // preserve behavior with older versions of this function by replacing
432 // higher vertex numbers by lower vertex numbers
433 std::sort(range_start,
434 range_end,
435 [](const std::pair<unsigned int, Point<spacedim>> &a,
436 const std::pair<unsigned int, Point<spacedim>> &b) {
437 return a.first < b.first;
438 });
439
440 // Now de-duplicate [range_start, range_end)
441 //
442 // We have identified all points that are within a strip of width 'tol'
443 // in one coordinate direction. Now we need to figure out which of these
444 // are also close in other coordinate directions. If two are close, we
445 // can mark the second one for deletion.
446 for (auto reference = range_start; reference != range_end; ++reference)
447 {
448 if (reference->first != numbers::invalid_unsigned_int)
449 for (auto it = reference + 1; it != range_end; ++it)
450 {
451 if (within_tolerance(reference->second, it->second))
452 {
453 new_vertex_numbers[it->first] = reference->first;
454 // skip the replaced vertex in the future
456 }
457 }
458 }
459 range_start = range_end;
460 }
461
462 // now we got a renumbering list. simply renumber all vertices
463 // (non-duplicate vertices get renumbered to themselves, so nothing bad
464 // happens). after that, the duplicate vertices will be unused, so call
465 // delete_unused_vertices() to do that part of the job.
466 for (auto &cell : cells)
467 for (auto &vertex_index : cell.vertices)
468 vertex_index = new_vertex_numbers[vertex_index];
469 for (auto &quad : subcelldata.boundary_quads)
470 for (auto &vertex_index : quad.vertices)
471 vertex_index = new_vertex_numbers[vertex_index];
472 for (auto &line : subcelldata.boundary_lines)
473 for (auto &vertex_index : line.vertices)
474 vertex_index = new_vertex_numbers[vertex_index];
475
476 delete_unused_vertices(vertices, cells, subcelldata);
477 }
478
479
480
481 template <int dim>
482 void
484 const double tol)
485 {
486 if (vertices.empty())
487 return;
488
489 // 1) map point to local vertex index
490 std::map<Point<dim>, unsigned int, FloatingPointComparator<double>>
491 map_point_to_local_vertex_index{FloatingPointComparator<double>(tol)};
492
493 // 2) initialize map with existing points uniquely
494 for (unsigned int i = 0; i < vertices.size(); ++i)
495 map_point_to_local_vertex_index[vertices[i]] = i;
496
497 // no duplicate points are found
498 if (map_point_to_local_vertex_index.size() == vertices.size())
499 return;
500
501 // 3) remove duplicate entries from vertices
502 vertices.resize(map_point_to_local_vertex_index.size());
503 {
504 unsigned int j = 0;
505 for (const auto &p : map_point_to_local_vertex_index)
506 vertices[j++] = p.first;
507 }
508 }
509
510
511
512 template <int dim, int spacedim>
513 std::size_t
515 const std::vector<Point<spacedim>> &all_vertices,
516 std::vector<CellData<dim>> &cells)
517 {
518 // This function is presently only implemented for volumetric (codimension
519 // 0) elements.
520
521 if (dim == 1)
522 return 0;
523 if (dim == 2 && spacedim == 3)
525
526 std::size_t n_negative_cells = 0;
527 std::size_t cell_no = 0;
528 for (auto &cell : cells)
529 {
530 const ArrayView<const unsigned int> vertices(cell.vertices);
531 // Some pathologically twisted cells can have exactly zero measure but
532 // we can still fix them
533 if (GridTools::cell_measure(all_vertices, vertices) <= 0)
534 {
535 ++n_negative_cells;
536 const auto reference_cell =
538
539 if (reference_cell.is_hyper_cube())
540 {
541 if (dim == 2)
542 {
543 // flip the cell across the y = x line in 2d
544 std::swap(cell.vertices[1], cell.vertices[2]);
545 }
546 else if (dim == 3)
547 {
548 // swap the front and back faces in 3d
549 std::swap(cell.vertices[0], cell.vertices[2]);
550 std::swap(cell.vertices[1], cell.vertices[3]);
551 std::swap(cell.vertices[4], cell.vertices[6]);
552 std::swap(cell.vertices[5], cell.vertices[7]);
553 }
554 }
555 else if (reference_cell.is_simplex())
556 {
557 // By basic rules for computing determinants we can just swap
558 // two vertices to fix a negative volume. Arbitrarily pick the
559 // last two.
560 std::swap(cell.vertices[cell.vertices.size() - 2],
561 cell.vertices[cell.vertices.size() - 1]);
562 }
563 else if (reference_cell == ReferenceCells::Wedge)
564 {
565 // swap the two triangular faces
566 std::swap(cell.vertices[0], cell.vertices[3]);
567 std::swap(cell.vertices[1], cell.vertices[4]);
568 std::swap(cell.vertices[2], cell.vertices[5]);
569 }
570 else if (reference_cell == ReferenceCells::Pyramid)
571 {
572 // Try swapping two vertices in the base - perhaps things were
573 // read in the UCD (counter-clockwise) order instead of lexical
574 std::swap(cell.vertices[2], cell.vertices[3]);
575 }
576 else
577 {
579 }
580 // Check whether the resulting cell is now ok.
581 // If not, then the grid is seriously broken and
582 // we just give up.
583 AssertThrow(GridTools::cell_measure(all_vertices, vertices) > 0,
584 ExcGridHasInvalidCell(cell_no));
585 }
586 ++cell_no;
587 }
588 return n_negative_cells;
589 }
590
591
592 template <int dim, int spacedim>
593 void
595 const std::vector<Point<spacedim>> &all_vertices,
596 std::vector<CellData<dim>> &cells)
597 {
598 const std::size_t n_negative_cells =
599 invert_cells_with_negative_measure(all_vertices, cells);
600
601 // We assume that all cells of a grid have
602 // either positive or negative volumes but
603 // not both mixed. Although above reordering
604 // might work also on single cells, grids
605 // with both kind of cells are very likely to
606 // be broken. Check for this here.
607 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
609 std::string(
610 "This function assumes that either all cells have positive "
611 "volume, or that all cells have been specified in an "
612 "inverted vertex order so that their volume is negative. "
613 "(In the latter case, this class automatically inverts "
614 "every cell.) However, the mesh you have specified "
615 "appears to have both cells with positive and cells with "
616 "negative volume. You need to check your mesh which "
617 "cells these are and how they got there.\n"
618 "As a hint, of the total ") +
619 std::to_string(cells.size()) + " cells in the mesh, " +
620 std::to_string(n_negative_cells) +
621 " appear to have a negative volume."));
622 }
623
624
625
626 // Functions and classes for consistently_order_cells
627 namespace
628 {
634 struct CheapEdge
635 {
639 CheapEdge(const unsigned int v0, const unsigned int v1)
640 : v0(v0)
641 , v1(v1)
642 {}
643
648 bool
649 operator<(const CheapEdge &e) const
650 {
651 return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
652 }
653
654 private:
658 const unsigned int v0, v1;
659 };
660
661
670 template <int dim>
671 bool
672 is_consistent(const std::vector<CellData<dim>> &cells)
673 {
674 std::set<CheapEdge> edges;
675
676 for (typename std::vector<CellData<dim>>::const_iterator c =
677 cells.begin();
678 c != cells.end();
679 ++c)
680 {
681 // construct the edges in reverse order. for each of them,
682 // ensure that the reverse edge is not yet in the list of
683 // edges (return false if the reverse edge already *is* in
684 // the list) and then add the actual edge to it; std::set
685 // eliminates duplicates automatically
686 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
687 {
688 const CheapEdge reverse_edge(
691 if (edges.find(reverse_edge) != edges.end())
692 return false;
693
694
695 // ok, not. insert edge in correct order
696 const CheapEdge correct_edge(
699 edges.insert(correct_edge);
700 }
701 }
702
703 // no conflicts found, so return true
704 return true;
705 }
706
707
714 template <int dim>
715 struct ParallelEdges
716 {
722 static const unsigned int starter_edges[dim];
723
728 static const unsigned int n_other_parallel_edges = (1 << (dim - 1)) - 1;
729 static const unsigned int
732 };
733
734 template <>
735 const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2};
736
737 template <>
738 const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1},
739 {0},
740 {3},
741 {2}};
742
743 template <>
744 const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8};
745
746 template <>
747 const unsigned int ParallelEdges<3>::parallel_edges[12][3] = {
748 {1, 4, 5}, // line 0
749 {0, 4, 5}, // line 1
750 {3, 6, 7}, // line 2
751 {2, 6, 7}, // line 3
752 {0, 1, 5}, // line 4
753 {0, 1, 4}, // line 5
754 {2, 3, 7}, // line 6
755 {2, 3, 6}, // line 7
756 {9, 10, 11}, // line 8
757 {8, 10, 11}, // line 9
758 {8, 9, 11}, // line 10
759 {8, 9, 10} // line 11
760 };
761
762
767 struct AdjacentCell
768 {
772 AdjacentCell()
775 {}
776
780 AdjacentCell(const unsigned int cell_index,
781 const unsigned int edge_within_cell)
784 {}
785
786
787 unsigned int cell_index;
788 unsigned int edge_within_cell;
789 };
790
791
792
793 template <int dim>
794 class AdjacentCells;
795
801 template <>
802 class AdjacentCells<2>
803 {
804 public:
809 using const_iterator = const AdjacentCell *;
810
819 void
820 push_back(const AdjacentCell &adjacent_cell)
821 {
823 adjacent_cells[0] = adjacent_cell;
824 else
825 {
829 adjacent_cells[1] = adjacent_cell;
830 }
831 }
832
833
838 const_iterator
839 begin() const
840 {
841 return adjacent_cells;
842 }
843
844
850 const_iterator
851 end() const
852 {
853 // check whether the current object stores zero, one, or two
854 // adjacent cells, and use this to point to the element past the
855 // last valid one
857 return adjacent_cells;
859 return adjacent_cells + 1;
860 else
861 return adjacent_cells + 2;
862 }
863
864 private:
871 AdjacentCell adjacent_cells[2];
872 };
873
874
875
883 template <>
884 class AdjacentCells<3> : public std::vector<AdjacentCell>
885 {};
886
887
897 template <int dim>
898 class Edge
899 {
900 public:
906 Edge(const CellData<dim> &cell, const unsigned int edge_number)
907 : orientation_status(not_oriented)
908 {
911
912 // copy vertices for this particular line
913 vertex_indices[0] =
914 cell
916 vertex_indices[1] =
917 cell
919
920 // bring them into standard orientation
921 if (vertex_indices[0] > vertex_indices[1])
922 std::swap(vertex_indices[0], vertex_indices[1]);
923 }
924
929 bool
930 operator<(const Edge<dim> &e) const
931 {
932 return ((vertex_indices[0] < e.vertex_indices[0]) ||
933 ((vertex_indices[0] == e.vertex_indices[0]) &&
934 (vertex_indices[1] < e.vertex_indices[1])));
935 }
936
940 bool
941 operator==(const Edge<dim> &e) const
942 {
943 return ((vertex_indices[0] == e.vertex_indices[0]) &&
944 (vertex_indices[1] == e.vertex_indices[1]));
945 }
946
951 unsigned int vertex_indices[2];
952
957 enum OrientationStatus
958 {
959 not_oriented,
960 forward,
961 backward
962 };
963
964 OrientationStatus orientation_status;
965
970 AdjacentCells<dim> adjacent_cells;
971 };
972
973
974
979 template <int dim>
980 struct Cell
981 {
987 Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
988 {
989 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
990 vertex_indices[i] = c.vertices[i];
991
992 // now for each of the edges of this cell, find the location inside the
993 // given edge_list array and store than index
994 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
995 {
996 const Edge<dim> e(c, l);
997 edge_indices[l] =
998 (std::lower_bound(edge_list.begin(), edge_list.end(), e) -
999 edge_list.begin());
1000 Assert(edge_indices[l] < edge_list.size(), ExcInternalError());
1001 Assert(edge_list[edge_indices[l]] == e, ExcInternalError());
1002 }
1003 }
1004
1009
1015 };
1016
1017
1018
1019 template <int dim>
1020 class EdgeDeltaSet;
1021
1031 template <>
1032 class EdgeDeltaSet<2>
1033 {
1034 public:
1038 using const_iterator = const unsigned int *;
1039
1044 EdgeDeltaSet()
1045 {
1047 }
1048
1049
1053 void
1054 clear()
1055 {
1057 }
1058
1063 void
1064 insert(const unsigned int edge_index)
1065 {
1067 edge_indices[0] = edge_index;
1068 else
1069 {
1072 edge_indices[1] = edge_index;
1073 }
1074 }
1075
1076
1080 const_iterator
1081 begin() const
1082 {
1083 return edge_indices;
1084 }
1085
1086
1090 const_iterator
1091 end() const
1092 {
1093 // check whether the current object stores zero, one, or two
1094 // indices, and use this to point to the element past the
1095 // last valid one
1097 return edge_indices;
1099 return edge_indices + 1;
1100 else
1101 return edge_indices + 2;
1102 }
1103
1104 private:
1108 unsigned int edge_indices[2];
1109 };
1110
1111
1112
1124 template <>
1125 class EdgeDeltaSet<3> : public std::set<unsigned int>
1126 {};
1127
1128
1129
1134 template <int dim>
1135 std::vector<Edge<dim>>
1136 build_edges(const std::vector<CellData<dim>> &cells)
1137 {
1138 // build the edge list for all cells. because each cell has
1139 // GeometryInfo<dim>::lines_per_cell edges, the total number
1140 // of edges is this many times the number of cells. of course
1141 // some of them will be duplicates, and we throw them out below
1142 std::vector<Edge<dim>> edge_list;
1143 edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
1144 for (unsigned int i = 0; i < cells.size(); ++i)
1145 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1146 edge_list.emplace_back(cells[i], l);
1147
1148 // next sort the edge list and then remove duplicates
1149 std::sort(edge_list.begin(), edge_list.end());
1150 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1151 edge_list.end());
1152
1153 return edge_list;
1154 }
1155
1156
1157
1162 template <int dim>
1163 std::vector<Cell<dim>>
1164 build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
1165 std::vector<Edge<dim>> &edges)
1166 {
1167 std::vector<Cell<dim>> cell_list;
1168 cell_list.reserve(cells.size());
1169 for (unsigned int i = 0; i < cells.size(); ++i)
1170 {
1171 // create our own data structure for the cells and let it
1172 // connect to the edges array
1173 cell_list.emplace_back(cells[i], edges);
1174
1175 // then also inform the edges that they are adjacent
1176 // to the current cell, and where within this cell
1177 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1178 edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back(
1179 AdjacentCell(i, l));
1180 }
1181 Assert(cell_list.size() == cells.size(), ExcInternalError());
1182
1183 return cell_list;
1184 }
1185
1186
1187
1192 template <int dim>
1193 unsigned int
1194 get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
1195 const std::vector<Edge<dim>> &edges,
1196 const unsigned int current_cell)
1197 {
1198 for (unsigned int c = current_cell; c < cells.size(); ++c)
1199 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1200 if (edges[cells[c].edge_indices[l]].orientation_status ==
1201 Edge<dim>::not_oriented)
1202 return c;
1203
1205 }
1206
1207
1208
1214 template <int dim>
1215 void
1216 orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> &cells,
1217 std::vector<Edge<dim>> &edges,
1218 const unsigned int cell,
1219 const unsigned int local_edge)
1220 {
1221 // choose the direction of the first edge. we have free choice
1222 // here and could simply choose "forward" if that's what pleases
1223 // us. however, for backward compatibility with the previous
1224 // implementation used till 2016, let us just choose the
1225 // direction so that it matches what we have in the given cell.
1226 //
1227 // in fact, in what can only be assumed to be a bug in the
1228 // original implementation, after orienting all edges, the code
1229 // that rotates the cells so that they match edge orientations
1230 // (see the rotate_cell() function below) rotated the cell two
1231 // more times by 90 degrees. this is ok -- it simply flips all
1232 // edge orientations, which leaves them valid. rather than do
1233 // the same in the current implementation, we can achieve the
1234 // same effect by modifying the rule above to choose the
1235 // direction of the starting edge of this parallel set
1236 // *opposite* to what it looks like in the current cell
1237 //
1238 // this bug only existed in the 2d implementation since there
1239 // were different implementations for 2d and 3d. consequently,
1240 // only replicate it for the 2d case and be "intuitive" in 3d.
1241 if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1243 local_edge, 0)])
1244 // orient initial edge *opposite* to the way it is in the cell
1245 // (see above for the reason)
1246 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1247 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1248 else
1249 {
1250 Assert(
1251 edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1252 cells[cell].vertex_indices
1255 Assert(
1256 edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
1257 cells[cell].vertex_indices
1260
1261 // orient initial edge *opposite* to the way it is in the cell
1262 // (see above for the reason)
1263 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1264 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1265 }
1266
1267 // walk outward from the given edge as described in
1268 // the algorithm in the paper that documents all of
1269 // this
1270 //
1271 // note that in 2d, each of the Deltas can at most
1272 // contain two elements, whereas in 3d it can be arbitrarily many
1273 EdgeDeltaSet<dim> Delta_k;
1274 EdgeDeltaSet<dim> Delta_k_minus_1;
1275 Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]);
1276
1277 while (Delta_k_minus_1.begin() !=
1278 Delta_k_minus_1.end()) // while set is not empty
1279 {
1280 Delta_k.clear();
1281
1282 for (typename EdgeDeltaSet<dim>::const_iterator delta =
1283 Delta_k_minus_1.begin();
1284 delta != Delta_k_minus_1.end();
1285 ++delta)
1286 {
1287 Assert(edges[*delta].orientation_status !=
1288 Edge<dim>::not_oriented,
1290
1291 // now go through the cells adjacent to this edge
1292 for (typename AdjacentCells<dim>::const_iterator adjacent_cell =
1293 edges[*delta].adjacent_cells.begin();
1294 adjacent_cell != edges[*delta].adjacent_cells.end();
1295 ++adjacent_cell)
1296 {
1297 const unsigned int K = adjacent_cell->cell_index;
1298 const unsigned int delta_is_edge_in_K =
1299 adjacent_cell->edge_within_cell;
1300
1301 // figure out the direction of delta with respect to the cell
1302 // K (in the orientation in which the user has given it to us)
1303 const unsigned int first_edge_vertex =
1304 (edges[*delta].orientation_status == Edge<dim>::forward ?
1305 edges[*delta].vertex_indices[0] :
1306 edges[*delta].vertex_indices[1]);
1307 const unsigned int first_edge_vertex_in_K =
1308 cells[K]
1310 delta_is_edge_in_K, 0)];
1311 Assert(
1312 first_edge_vertex == first_edge_vertex_in_K ||
1313 first_edge_vertex ==
1315 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1317
1318 // now figure out which direction the each of the "opposite"
1319 // edges needs to be oriented into.
1320 for (unsigned int o_e = 0;
1321 o_e < ParallelEdges<dim>::n_other_parallel_edges;
1322 ++o_e)
1323 {
1324 // get the index of the opposite edge and select which its
1325 // first vertex needs to be based on how the current edge
1326 // is oriented in the current cell
1327 const unsigned int opposite_edge =
1328 cells[K].edge_indices[ParallelEdges<
1329 dim>::parallel_edges[delta_is_edge_in_K][o_e]];
1330 const unsigned int first_opposite_edge_vertex =
1331 cells[K].vertex_indices
1333 ParallelEdges<
1334 dim>::parallel_edges[delta_is_edge_in_K][o_e],
1335 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1336 1))];
1337
1338 // then determine the orientation of the edge based on
1339 // whether the vertex we want to be the edge's first
1340 // vertex is already the first vertex of the edge, or
1341 // whether it points in the opposite direction
1342 const typename Edge<dim>::OrientationStatus
1343 opposite_edge_orientation =
1344 (edges[opposite_edge].vertex_indices[0] ==
1345 first_opposite_edge_vertex ?
1346 Edge<dim>::forward :
1347 Edge<dim>::backward);
1348
1349 // see if the opposite edge (there is only one in 2d) has
1350 // already been oriented.
1351 if (edges[opposite_edge].orientation_status ==
1352 Edge<dim>::not_oriented)
1353 {
1354 // the opposite edge is not yet oriented. do orient it
1355 // and add it to Delta_k
1356 edges[opposite_edge].orientation_status =
1357 opposite_edge_orientation;
1358 Delta_k.insert(opposite_edge);
1359 }
1360 else
1361 {
1362 // this opposite edge has already been oriented. it
1363 // should be consistent with the current one in 2d,
1364 // while in 3d it may in fact be mis-oriented, and in
1365 // that case the mesh will not be orientable. indicate
1366 // this by throwing an exception that we can catch
1367 // further up; this has the advantage that we can
1368 // propagate through a couple of functions without
1369 // having to do error checking and without modifying
1370 // the 'cells' array that the user gave us
1371 if (dim == 2)
1372 {
1373 Assert(edges[opposite_edge].orientation_status ==
1374 opposite_edge_orientation,
1376 }
1377 else if (dim == 3)
1378 {
1379 if (edges[opposite_edge].orientation_status !=
1380 opposite_edge_orientation)
1381 throw ExcMeshNotOrientable();
1382 }
1383 else
1385 }
1386 }
1387 }
1388 }
1389
1390 // finally copy the new set to the previous one
1391 // (corresponding to increasing 'k' by one in the
1392 // algorithm)
1393 Delta_k_minus_1 = Delta_k;
1394 }
1395 }
1396
1397
1405 template <int dim>
1406 void
1407 rotate_cell(const std::vector<Cell<dim>> &cell_list,
1408 const std::vector<Edge<dim>> &edge_list,
1409 const unsigned int cell_index,
1410 std::vector<CellData<dim>> &raw_cells)
1411 {
1412 // find the first vertex of the cell. this is the vertex where dim edges
1413 // originate, so for each of the edges record which the starting vertex is
1414 unsigned int starting_vertex_of_edge[GeometryInfo<dim>::lines_per_cell];
1415 for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
1416 {
1417 Assert(edge_list[cell_list[cell_index].edge_indices[e]]
1418 .orientation_status != Edge<dim>::not_oriented,
1420 if (edge_list[cell_list[cell_index].edge_indices[e]]
1421 .orientation_status == Edge<dim>::forward)
1422 starting_vertex_of_edge[e] =
1423 edge_list[cell_list[cell_index].edge_indices[e]]
1424 .vertex_indices[0];
1425 else
1426 starting_vertex_of_edge[e] =
1427 edge_list[cell_list[cell_index].edge_indices[e]]
1428 .vertex_indices[1];
1429 }
1430
1431 // find the vertex number that appears dim times. this will then be
1432 // the vertex at which we want to locate the origin of the cell's
1433 // coordinate system (i.e., vertex 0)
1434 unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int;
1435 switch (dim)
1436 {
1437 case 2:
1438 {
1439 // in 2d, we can simply enumerate the possibilities where the
1440 // origin may be located because edges zero and one don't share
1441 // any vertices, and the same for edges two and three
1442 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1443 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1444 origin_vertex_of_cell = starting_vertex_of_edge[0];
1445 else if ((starting_vertex_of_edge[1] ==
1446 starting_vertex_of_edge[2]) ||
1447 (starting_vertex_of_edge[1] ==
1448 starting_vertex_of_edge[3]))
1449 origin_vertex_of_cell = starting_vertex_of_edge[1];
1450 else
1452
1453 break;
1454 }
1455
1456 case 3:
1457 {
1458 // one could probably do something similar in 3d, but that seems
1459 // more complicated than one wants to write down. just go
1460 // through the list of possible starting vertices and check
1461 for (origin_vertex_of_cell = 0;
1462 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1463 ++origin_vertex_of_cell)
1464 if (std::count(starting_vertex_of_edge,
1465 starting_vertex_of_edge +
1467 cell_list[cell_index]
1468 .vertex_indices[origin_vertex_of_cell]) == dim)
1469 break;
1470 Assert(origin_vertex_of_cell <
1473
1474 break;
1475 }
1476
1477 default:
1479 }
1480
1481 // now rotate raw_cells[cell_index] in such a way that its orientation
1482 // matches that of cell_list[cell_index]
1483 switch (dim)
1484 {
1485 case 2:
1486 {
1487 // in 2d, we can literally rotate the cell until its origin
1488 // matches the one that we have determined above should be
1489 // the origin vertex
1490 //
1491 // when doing a rotation, take into account the ordering of
1492 // vertices (not in clockwise or counter-clockwise sense)
1493 while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell)
1494 {
1495 const unsigned int tmp = raw_cells[cell_index].vertices[0];
1496 raw_cells[cell_index].vertices[0] =
1497 raw_cells[cell_index].vertices[1];
1498 raw_cells[cell_index].vertices[1] =
1499 raw_cells[cell_index].vertices[3];
1500 raw_cells[cell_index].vertices[3] =
1501 raw_cells[cell_index].vertices[2];
1502 raw_cells[cell_index].vertices[2] = tmp;
1503 }
1504 break;
1505 }
1506
1507 case 3:
1508 {
1509 // in 3d, the situation is a bit more complicated. from above, we
1510 // now know which vertex is at the origin (because 3 edges
1511 // originate from it), but that still leaves 3 possible rotations
1512 // of the cube. the important realization is that we can choose
1513 // any of them: in all 3 rotations, all edges originate from the
1514 // one vertex, and that fixes the directions of all 12 edges in
1515 // the cube because these 3 cover all 3 equivalence classes!
1516 // consequently, we can select an arbitrary one among the
1517 // permutations -- for example the following ones:
1518 static const unsigned int cube_permutations[8][8] = {
1519 {0, 1, 2, 3, 4, 5, 6, 7},
1520 {1, 5, 3, 7, 0, 4, 2, 6},
1521 {2, 6, 0, 4, 3, 7, 1, 5},
1522 {3, 2, 1, 0, 7, 6, 5, 4},
1523 {4, 0, 6, 2, 5, 1, 7, 3},
1524 {5, 4, 7, 6, 1, 0, 3, 2},
1525 {6, 7, 4, 5, 2, 3, 0, 1},
1526 {7, 3, 5, 1, 6, 2, 4, 0}};
1527
1528 unsigned int
1529 temp_vertex_indices[GeometryInfo<dim>::vertices_per_cell];
1530 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1531 temp_vertex_indices[v] =
1532 raw_cells[cell_index]
1533 .vertices[cube_permutations[origin_vertex_of_cell][v]];
1534 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1535 raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
1536
1537 break;
1538 }
1539
1540 default:
1541 {
1543 }
1544 }
1545 }
1546
1547
1553 template <int dim>
1554 void
1555 reorient(std::vector<CellData<dim>> &cells)
1556 {
1557 // first build the arrays that connect cells to edges and the other
1558 // way around
1559 std::vector<Edge<dim>> edge_list = build_edges(cells);
1560 std::vector<Cell<dim>> cell_list =
1561 build_cells_and_connect_edges(cells, edge_list);
1562
1563 // then loop over all cells and start orienting parallel edge sets
1564 // of cells that still have non-oriented edges
1565 unsigned int next_cell_with_unoriented_edge = 0;
1566 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1567 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1569 {
1570 // see which edge sets are still not oriented
1571 //
1572 // we do not need to look at each edge because if we orient edge
1573 // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
1574 // will be 3 other edges that are also oriented). there are only
1575 // dim independent sets of edges, so loop over these.
1576 //
1577 // we need to check whether each one of these starter edges may
1578 // already be oriented because the line (sheet) that connects
1579 // globally parallel edges may be self-intersecting in the
1580 // current cell
1581 for (unsigned int l = 0; l < dim; ++l)
1582 if (edge_list[cell_list[next_cell_with_unoriented_edge]
1583 .edge_indices[ParallelEdges<dim>::starter_edges[l]]]
1584 .orientation_status == Edge<dim>::not_oriented)
1585 orient_one_set_of_parallel_edges(
1586 cell_list,
1587 edge_list,
1588 next_cell_with_unoriented_edge,
1589 ParallelEdges<dim>::starter_edges[l]);
1590
1591 // ensure that we have really oriented all edges now, not just
1592 // the starter edges
1593 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1594 Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
1595 .edge_indices[l]]
1596 .orientation_status != Edge<dim>::not_oriented,
1598 }
1599
1600 // now that we have oriented all edges, we need to rotate cells
1601 // so that the edges point in the right direction with the now
1602 // rotated coordinate system
1603 for (unsigned int c = 0; c < cells.size(); ++c)
1604 rotate_cell(cell_list, edge_list, c, cells);
1605 }
1606
1607
1608 // overload of the function above for 1d -- there is nothing
1609 // to orient in that case
1610 void
1611 reorient(std::vector<CellData<1>> &)
1612 {}
1613 } // namespace
1614
1615
1616
1617 template <int dim>
1618 void
1620 {
1621 Assert(cells.size() != 0,
1622 ExcMessage(
1623 "List of elements to orient must have at least one cell"));
1624
1625 // there is nothing for us to do in 1d
1626 if (dim == 1)
1627 return;
1628
1629 // check if grids are already consistent. if so, do
1630 // nothing. if not, then do the reordering
1631 if (!is_consistent(cells))
1632 try
1633 {
1634 reorient(cells);
1635 }
1636 catch (const ExcMeshNotOrientable &)
1637 {
1638 // the mesh is not orientable. this is acceptable if we are in 3d,
1639 // as class Triangulation knows how to handle this, but it is
1640 // not in 2d; in that case, re-throw the exception
1641 if (dim < 3)
1642 throw;
1643 }
1644 }
1645
1646
1647
1648 template <int dim, int spacedim>
1649 std::map<unsigned int, Point<spacedim>>
1651 {
1652 std::map<unsigned int, Point<spacedim>> vertex_map;
1654 cell = tria.begin_active(),
1655 endc = tria.end();
1656 for (; cell != endc; ++cell)
1657 {
1658 for (const unsigned int i : cell->face_indices())
1659 {
1660 const typename Triangulation<dim, spacedim>::face_iterator &face =
1661 cell->face(i);
1662 if (face->at_boundary())
1663 {
1664 for (unsigned j = 0; j < face->n_vertices(); ++j)
1665 {
1666 const Point<spacedim> &vertex = face->vertex(j);
1667 const unsigned int vertex_index = face->vertex_index(j);
1668 vertex_map[vertex_index] = vertex;
1669 }
1670 }
1671 }
1672 }
1673 return vertex_map;
1674 }
1675
1676
1677
1678 template <int dim, int spacedim>
1679 void
1681 const bool isotropic,
1682 const unsigned int max_iterations)
1683 {
1684 unsigned int iter = 0;
1685 bool continue_refinement = true;
1686
1687 while (continue_refinement && (iter < max_iterations))
1688 {
1689 if (max_iterations != numbers::invalid_unsigned_int)
1690 ++iter;
1691 continue_refinement = false;
1692
1693 for (const auto &cell : tria.active_cell_iterators())
1694 for (const unsigned int j : cell->face_indices())
1695 if (cell->at_boundary(j) == false &&
1696 cell->neighbor(j)->has_children())
1697 {
1698 if (isotropic)
1699 {
1700 cell->set_refine_flag();
1701 continue_refinement = true;
1702 }
1703 else
1704 continue_refinement |= cell->flag_for_face_refinement(j);
1705 }
1706
1708 }
1709 }
1710
1711
1712
1713 template <int dim, int spacedim>
1714 void
1716 const double max_ratio,
1717 const unsigned int max_iterations)
1718 {
1719 unsigned int iter = 0;
1720 bool continue_refinement = true;
1721
1722 while (continue_refinement && (iter < max_iterations))
1723 {
1724 ++iter;
1725 continue_refinement = false;
1726 for (const auto &cell : tria.active_cell_iterators())
1727 {
1728 std::pair<unsigned int, double> info =
1729 GridTools::get_longest_direction<dim, spacedim>(cell);
1730 if (info.second > max_ratio)
1731 {
1732 cell->set_refine_flag(
1733 RefinementCase<dim>::cut_axis(info.first));
1734 continue_refinement = true;
1735 }
1736 }
1738 }
1739 }
1740
1741
1742
1743 template <int dim, int spacedim>
1744 std::map<unsigned int, Point<spacedim>>
1746 const Mapping<dim, spacedim> &mapping)
1747 {
1748 std::map<unsigned int, Point<spacedim>> result;
1749 for (const auto &cell : container.active_cell_iterators())
1750 {
1751 if (!cell->is_artificial())
1752 {
1753 const auto vs = mapping.get_vertices(cell);
1754 for (unsigned int i = 0; i < vs.size(); ++i)
1755 result[cell->vertex_index(i)] = vs[i];
1756 }
1757 }
1758 return result;
1759 }
1760
1761
1762
1763 template <int dim, int spacedim>
1764 std::vector<
1765 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1767 {
1768 std::vector<
1769 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1770 vertex_to_cell_map(triangulation.n_vertices());
1772 cell = triangulation.begin_active(),
1773 endc = triangulation.end();
1774 for (; cell != endc; ++cell)
1775 for (const unsigned int i : cell->vertex_indices())
1776 vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
1777
1778 // Check if mesh has hanging nodes. Do this only locally to
1779 // prevent communication and possible deadlock.
1780 if (triangulation.Triangulation<dim, spacedim>::has_hanging_nodes())
1781 {
1782 Assert(triangulation.all_reference_cells_are_hyper_cube(),
1784
1785 // Take care of hanging nodes
1786 cell = triangulation.begin_active();
1787 for (; cell != endc; ++cell)
1788 {
1789 for (const unsigned int i : cell->face_indices())
1790 {
1791 if ((cell->at_boundary(i) == false) &&
1792 (cell->neighbor(i)->is_active()))
1793 {
1795 adjacent_cell = cell->neighbor(i);
1796 for (unsigned int j = 0; j < cell->face(i)->n_vertices();
1797 ++j)
1798 vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(
1799 adjacent_cell);
1800 }
1801 }
1802
1803 // in 3d also loop over the edges
1804 if (dim == 3)
1805 {
1806 for (unsigned int i = 0; i < cell->n_lines(); ++i)
1807 if (cell->line(i)->has_children())
1808 // the only place where this vertex could have been
1809 // hiding is on the mid-edge point of the edge we
1810 // are looking at
1811 vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)]
1812 .insert(cell);
1813 }
1814 }
1815 }
1816
1817 return vertex_to_cell_map;
1818 }
1819
1820
1821
1822 template <int dim, int spacedim>
1823 void
1826 DynamicSparsityPattern &cell_connectivity)
1827 {
1828 cell_connectivity.reinit(triangulation.n_active_cells(),
1829 triangulation.n_active_cells());
1830
1831 // loop over all cells and their neighbors to build the sparsity
1832 // pattern. note that it's a bit hard to enter all the connections when a
1833 // neighbor has children since we would need to find out which of its
1834 // children is adjacent to the current cell. this problem can be omitted
1835 // if we only do something if the neighbor has no children -- in that case
1836 // it is either on the same or a coarser level than we are. in return, we
1837 // have to add entries in both directions for both cells
1838 for (const auto &cell : triangulation.active_cell_iterators())
1839 {
1840 const unsigned int index = cell->active_cell_index();
1841 cell_connectivity.add(index, index);
1842 for (auto f : cell->face_indices())
1843 if ((cell->at_boundary(f) == false) &&
1844 (cell->neighbor(f)->has_children() == false))
1845 {
1846 const unsigned int other_index =
1847 cell->neighbor(f)->active_cell_index();
1848 cell_connectivity.add(index, other_index);
1849 cell_connectivity.add(other_index, index);
1850 }
1851 }
1852 }
1853
1854
1855
1856 template <int dim, int spacedim>
1857 void
1860 DynamicSparsityPattern &cell_connectivity)
1861 {
1862 // The choice of 16 or fewer neighbors here is based on empirical
1863 // measurements.
1864 //
1865 // Vertices in a structured hexahedral mesh have 8 adjacent cells. In a
1866 // structured tetrahedral mesh, about 98% of vertices have 16 neighbors or
1867 // fewer. Similarly, in an unstructured tetrahedral mesh, if we count the
1868 // number of neighbors each vertex has we obtain the following distribution:
1869 //
1870 // 3, 1
1871 // 4, 728
1872 // 5, 4084
1873 // 6, 7614
1874 // 7, 17329
1875 // 8, 31145
1876 // 9, 46698
1877 // 10, 64193
1878 // 11, 68269
1879 // 12, 63574
1880 // 13, 57016
1881 // 14, 50476
1882 // 15, 41886
1883 // 16, 31820
1884 // 17, 21269
1885 // 18, 12217
1886 // 19, 6072
1887 // 20, 2527
1888 // 21, 825
1889 // 22, 262
1890 // 23, 61
1891 // 24, 12
1892 // 26, 1
1893 //
1894 // so about 86% of vertices have 16 neighbors or fewer. Hence, we picked 16
1895 // neighbors here to cover most cases without allocation.
1896 std::vector<boost::container::small_vector<unsigned int, 16>>
1897 vertex_to_cell(triangulation.n_vertices());
1898 for (const auto &cell : triangulation.active_cell_iterators())
1899 {
1900 for (const unsigned int v : cell->vertex_indices())
1901 vertex_to_cell[cell->vertex_index(v)].push_back(
1902 cell->active_cell_index());
1903 }
1904
1905 cell_connectivity.reinit(triangulation.n_active_cells(),
1906 triangulation.n_active_cells());
1907 std::vector<types::global_dof_index> neighbors;
1908 for (const auto &cell : triangulation.active_cell_iterators())
1909 {
1910 neighbors.clear();
1911 for (const unsigned int v : cell->vertex_indices())
1912 neighbors.insert(neighbors.end(),
1913 vertex_to_cell[cell->vertex_index(v)].begin(),
1914 vertex_to_cell[cell->vertex_index(v)].end());
1915 std::sort(neighbors.begin(), neighbors.end());
1916 cell_connectivity.add_entries(cell->active_cell_index(),
1917 neighbors.begin(),
1918 std::unique(neighbors.begin(),
1919 neighbors.end()),
1920 true);
1921 }
1922 }
1923
1924
1925 template <int dim, int spacedim>
1926 void
1929 const unsigned int level,
1930 DynamicSparsityPattern &cell_connectivity)
1931 {
1932 std::vector<boost::container::small_vector<unsigned int, 16>>
1933 vertex_to_cell(triangulation.n_vertices());
1935 triangulation.begin(level);
1936 cell != triangulation.end(level);
1937 ++cell)
1938 {
1939 for (const unsigned int v : cell->vertex_indices())
1940 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
1941 }
1942
1943 cell_connectivity.reinit(triangulation.n_cells(level),
1944 triangulation.n_cells(level));
1945 std::vector<types::global_dof_index> neighbors;
1946 for (const auto &cell : triangulation.cell_iterators_on_level(level))
1947 {
1948 neighbors.clear();
1949 for (const unsigned int v : cell->vertex_indices())
1950 neighbors.insert(neighbors.end(),
1951 vertex_to_cell[cell->vertex_index(v)].begin(),
1952 vertex_to_cell[cell->vertex_index(v)].end());
1953 std::sort(neighbors.begin(), neighbors.end());
1954 cell_connectivity.add_entries(cell->index(),
1955 neighbors.begin(),
1956 std::unique(neighbors.begin(),
1957 neighbors.end()),
1958 true);
1959 }
1960 }
1961} /* namespace GridTools */
1962
1963// explicit instantiations
1964#include "grid/grid_tools_topology.inst"
1965
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
std::size_t size() const
Definition array_view.h:689
Number side_length(const unsigned int direction) const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
void add(const size_type i, const size_type j)
void insert_face_data(const FaceIteratorType &)
void insert_face_data(const FaceIteratorType &face)
std::set< CellData< dim - 1 >, internal::CellDataComparator< dim - 1 > > face_data
Abstract base class for mapping classes.
Definition mapping.h:320
virtual boost::container::small_vector< Point< spacedim >, ReferenceCells::max_n_vertices< dim >() > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition point.h:113
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_levels() const
cell_iterator end() const
virtual void execute_coarsening_and_refinement()
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int level
Definition grid_out.cc:4635
unsigned int edge_within_cell
unsigned int edge_indices[GeometryInfo< dim >::lines_per_cell]
AdjacentCell adjacent_cells[2]
static const unsigned int n_other_parallel_edges
static const unsigned int starter_edges[dim]
static const unsigned int parallel_edges[GeometryInfo< dim >::lines_per_cell][n_other_parallel_edges]
OrientationStatus orientation_status
unsigned int vertex_indices[2]
const unsigned int v0
const unsigned int v1
unsigned int cell_index
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcMeshNotOrientable()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcGridHasInvalidCell(int arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void consistently_order_cells(std::vector< CellData< dim > > &cells)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void append_face_data(const CellData< 1 > &face_data, SubCellData &subcell_data)
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
void delete_unused_vertices(std::vector< Point< spacedim > > &vertices, std::vector< CellData< dim > > &cells, SubCellData &subcelldata)
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
std::size_t invert_cells_with_negative_measure(const std::vector< Point< spacedim > > &all_vertices, std::vector< CellData< dim > > &cells)
void delete_duplicated_vertices(std::vector< Point< spacedim > > &all_vertices, std::vector< CellData< dim > > &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
void invert_all_negative_measure_cells(const std::vector< Point< spacedim > > &all_vertices, std::vector< CellData< dim > > &cells)
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
double cell_measure(const std::vector< Point< dim > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::boundary_id internal_face_boundary_id
Definition types.h:323
constexpr types::manifold_id flat_manifold_id
Definition types.h:336
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< unsigned int > vertices
Definition cell_data.h:85
types::manifold_id manifold_id
Definition cell_data.h:126
types::material_id material_id
Definition cell_data.h:104
types::boundary_id boundary_id
Definition cell_data.h:115
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
bool operator()(const CellData< structdim > &a, const CellData< structdim > &b) const
std::vector< CellData< 2 > > boundary_quads
Definition cell_data.h:248
bool check_consistency(const unsigned int dim) const
std::vector< CellData< 1 > > boundary_lines
Definition cell_data.h:232
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)