deal.II version GIT relicensing-2233-g9bd70d8389 2024-12-08 16:10: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 <algorithm>
24#include <map>
25#include <set>
26#include <vector>
27
29
30namespace GridTools
31{
32 // Generic functions for appending face data in 2d or 3d. TODO: we can
33 // remove these once we have 'if constexpr'.
34 namespace internal
35 {
36 inline void
37 append_face_data(const CellData<1> &face_data, SubCellData &subcell_data)
38 {
39 subcell_data.boundary_lines.push_back(face_data);
40 }
41
42
43
44 inline void
45 append_face_data(const CellData<2> &face_data, SubCellData &subcell_data)
46 {
47 subcell_data.boundary_quads.push_back(face_data);
48 }
49
50
51
52 // Lexical comparison for sorting CellData objects.
53 template <int structdim>
55 {
56 bool
58 const CellData<structdim> &b) const
59 {
60 // Check vertices:
61 if (std::lexicographical_compare(std::begin(a.vertices),
62 std::end(a.vertices),
63 std::begin(b.vertices),
64 std::end(b.vertices)))
65 return true;
66 // it should never be necessary to check the material or manifold
67 // ids as a 'tiebreaker' (since they must be equal if the vertex
68 // indices are equal). Assert it anyway:
69#ifdef DEBUG
70 if (std::equal(std::begin(a.vertices),
71 std::end(a.vertices),
72 std::begin(b.vertices)))
73 {
74 Assert(a.material_id == b.material_id &&
75 a.manifold_id == b.manifold_id,
77 "Two CellData objects with equal vertices must "
78 "have the same material/boundary ids and manifold "
79 "ids."));
80 }
81#endif
82 return false;
83 }
84 };
85
86
96 template <int dim>
98 {
99 public:
103 template <typename FaceIteratorType>
104 void
105 insert_face_data(const FaceIteratorType &face)
106 {
107 CellData<dim - 1> face_cell_data(face->n_vertices());
108 for (unsigned int vertex_n = 0; vertex_n < face->n_vertices();
109 ++vertex_n)
110 face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
111 face_cell_data.boundary_id = face->boundary_id();
112 face_cell_data.manifold_id = face->manifold_id();
113
114 face_data.insert(std::move(face_cell_data));
115 }
116
122 {
123 SubCellData subcell_data;
124
125 for (const CellData<dim - 1> &face_cell_data : face_data)
126 internal::append_face_data(face_cell_data, subcell_data);
127 return subcell_data;
128 }
129
130
131 private:
132 std::set<CellData<dim - 1>, internal::CellDataComparator<dim - 1>>
134 };
135
136
137 // Do nothing for dim=1:
138 template <>
140 {
141 public:
142 template <typename FaceIteratorType>
143 void
144 insert_face_data(const FaceIteratorType &)
145 {}
146
149 {
150 return SubCellData();
151 }
152 };
153 } // namespace internal
154
155
156
157 template <int dim, int spacedim>
158 std::
159 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
161 {
162 Assert(tria.n_levels() >= 1,
163 ExcMessage("The input triangulation must be non-empty."));
164
165 std::vector<Point<spacedim>> vertices = tria.get_vertices();
166 std::vector<CellData<dim>> cells;
167
169 std::set<CellData<1>, internal::CellDataComparator<1>>
170 line_data; // only used in 3d
171
172 for (const auto &cell : tria.cell_iterators_on_level(0))
173 {
174 // Save cell data
175 CellData<dim> cell_data(cell->n_vertices());
176 for (const unsigned int cell_vertex_n : cell->vertex_indices())
177 {
178 Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
180 cell_data.vertices[cell_vertex_n] =
181 cell->vertex_index(cell_vertex_n);
182 }
183 cell_data.material_id = cell->material_id();
184 cell_data.manifold_id = cell->manifold_id();
185 cells.emplace_back(std::move(cell_data));
186
187 // Save face data
188 if (dim > 1)
189 {
190 for (const unsigned int face_n : cell->face_indices())
191 // We don't need to insert anything if we have default values
192 {
193 const auto face = cell->face(face_n);
194 if (face->boundary_id() != numbers::internal_face_boundary_id ||
195 face->manifold_id() != numbers::flat_manifold_id)
196 face_data.insert_face_data(face);
197 }
198 }
199 // Save line data
200 if (dim == 3)
201 {
202 for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
203 {
204 const auto line = cell->line(line_n);
205 // We don't need to insert anything if we have default values
206 if (line->boundary_id() != numbers::internal_face_boundary_id ||
207 line->manifold_id() != numbers::flat_manifold_id)
208 {
209 CellData<1> line_cell_data(line->n_vertices());
210 for (const unsigned int vertex_n : line->vertex_indices())
211 line_cell_data.vertices[vertex_n] =
212 line->vertex_index(vertex_n);
213 line_cell_data.boundary_id = line->boundary_id();
214 line_cell_data.manifold_id = line->manifold_id();
215 line_data.insert(std::move(line_cell_data));
216 }
217 }
218 }
219 }
220
221 SubCellData subcell_data = face_data.get();
222
223 if (dim == 3)
224 for (const CellData<1> &face_line_data : line_data)
225 subcell_data.boundary_lines.push_back(face_line_data);
226
227 // We end up with a 'vertices' array that uses some of the entries,
228 // but not all -- specifically, all vertices referenced by level-0
229 // cells. We can compress the array:
230 GridTools::delete_unused_vertices(vertices, cells, subcell_data);
231
232 return std::tuple<std::vector<Point<spacedim>>,
233 std::vector<CellData<dim>>,
234 SubCellData>(std::move(vertices),
235 std::move(cells),
236 std::move(subcell_data));
237 }
238
239
240
241 template <int dim, int spacedim>
242 void
244 std::vector<CellData<dim>> &cells,
245 SubCellData &subcelldata)
246 {
247 Assert(
248 subcelldata.check_consistency(dim),
250 "Invalid SubCellData supplied according to ::check_consistency(). "
251 "This is caused by data containing objects for the wrong dimension."));
252
253 // first check which vertices are actually used
254 std::vector<bool> vertex_used(vertices.size(), false);
255 for (unsigned int c = 0; c < cells.size(); ++c)
256 for (unsigned int v = 0; v < cells[c].vertices.size(); ++v)
257 {
258 Assert(cells[c].vertices[v] < vertices.size(),
259 ExcMessage("Invalid vertex index encountered! cells[" +
260 Utilities::int_to_string(c) + "].vertices[" +
261 Utilities::int_to_string(v) + "]=" +
262 Utilities::int_to_string(cells[c].vertices[v]) +
263 " is invalid, because only " +
264 Utilities::int_to_string(vertices.size()) +
265 " vertices were supplied."));
266 vertex_used[cells[c].vertices[v]] = true;
267 }
268
269
270 // then renumber the vertices that are actually used in the same order as
271 // they were beforehand
272 const unsigned int invalid_vertex = numbers::invalid_unsigned_int;
273 std::vector<unsigned int> new_vertex_numbers(vertices.size(),
274 invalid_vertex);
275 unsigned int next_free_number = 0;
276 for (unsigned int i = 0; i < vertices.size(); ++i)
277 if (vertex_used[i] == true)
278 {
279 new_vertex_numbers[i] = next_free_number;
280 ++next_free_number;
281 }
282
283 // next replace old vertex numbers by the new ones
284 for (unsigned int c = 0; c < cells.size(); ++c)
285 for (auto &v : cells[c].vertices)
286 v = new_vertex_numbers[v];
287
288 // same for boundary data
289 for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT
290 ++c)
291 for (unsigned int v = 0;
292 v < subcelldata.boundary_lines[c].vertices.size();
293 ++v)
294 {
295 Assert(subcelldata.boundary_lines[c].vertices[v] <
296 new_vertex_numbers.size(),
298 "Invalid vertex index in subcelldata.boundary_lines. "
299 "subcelldata.boundary_lines[" +
300 Utilities::int_to_string(c) + "].vertices[" +
301 Utilities::int_to_string(v) + "]=" +
303 subcelldata.boundary_lines[c].vertices[v]) +
304 " is invalid, because only " +
305 Utilities::int_to_string(vertices.size()) +
306 " vertices were supplied."));
307 subcelldata.boundary_lines[c].vertices[v] =
308 new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
309 }
310
311 for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT
312 ++c)
313 for (unsigned int v = 0;
314 v < subcelldata.boundary_quads[c].vertices.size();
315 ++v)
316 {
317 Assert(subcelldata.boundary_quads[c].vertices[v] <
318 new_vertex_numbers.size(),
320 "Invalid vertex index in subcelldata.boundary_quads. "
321 "subcelldata.boundary_quads[" +
322 Utilities::int_to_string(c) + "].vertices[" +
323 Utilities::int_to_string(v) + "]=" +
325 subcelldata.boundary_quads[c].vertices[v]) +
326 " is invalid, because only " +
327 Utilities::int_to_string(vertices.size()) +
328 " vertices were supplied."));
329
330 subcelldata.boundary_quads[c].vertices[v] =
331 new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
332 }
333
334 // finally copy over the vertices which we really need to a new array and
335 // replace the old one by the new one
336 std::vector<Point<spacedim>> tmp;
337 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(), true));
338 for (unsigned int v = 0; v < vertices.size(); ++v)
339 if (vertex_used[v] == true)
340 tmp.push_back(vertices[v]);
341 swap(vertices, tmp);
342 }
343
344
345
346 template <int dim, int spacedim>
347 void
349 std::vector<CellData<dim>> &cells,
350 SubCellData &subcelldata,
351 std::vector<unsigned int> &considered_vertices,
352 const double tol)
353 {
354 if (tol == 0.0)
355 return; // nothing to do per definition
356
357 AssertIndexRange(2, vertices.size());
358 std::vector<unsigned int> new_vertex_numbers(vertices.size());
359 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
360
361 // if the considered_vertices vector is empty, consider all vertices
362 if (considered_vertices.empty())
363 considered_vertices = new_vertex_numbers;
364 Assert(considered_vertices.size() <= vertices.size(), ExcInternalError());
365
366 // The algorithm below improves upon the naive O(n^2) algorithm by first
367 // sorting vertices by their value in one component and then only
368 // comparing vertices for equality which are nearly equal in that
369 // component. For example, if @p vertices form a cube, then we will only
370 // compare points that have the same x coordinate when we try to find
371 // duplicated vertices.
372
373 // Start by finding the longest coordinate direction. This minimizes the
374 // number of points that need to be compared against each-other in a
375 // single set for typical geometries.
376 const BoundingBox<spacedim> bbox(vertices);
377
378 unsigned int longest_coordinate_direction = 0;
379 double longest_coordinate_length = bbox.side_length(0);
380 for (unsigned int d = 1; d < spacedim; ++d)
381 {
382 const double coordinate_length = bbox.side_length(d);
383 if (longest_coordinate_length < coordinate_length)
384 {
385 longest_coordinate_length = coordinate_length;
386 longest_coordinate_direction = d;
387 }
388 }
389
390 // Sort vertices (while preserving their vertex numbers) along that
391 // coordinate direction:
392 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
393 sorted_vertices.reserve(vertices.size());
394 for (const unsigned int vertex_n : considered_vertices)
395 {
396 AssertIndexRange(vertex_n, vertices.size());
397 sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]);
398 }
399 std::sort(sorted_vertices.begin(),
400 sorted_vertices.end(),
401 [&](const std::pair<unsigned int, Point<spacedim>> &a,
402 const std::pair<unsigned int, Point<spacedim>> &b) {
403 return a.second[longest_coordinate_direction] <
404 b.second[longest_coordinate_direction];
405 });
406
407 auto within_tolerance = [=](const Point<spacedim> &a,
408 const Point<spacedim> &b) {
409 for (unsigned int d = 0; d < spacedim; ++d)
410 if (std::abs(a[d] - b[d]) > tol)
411 return false;
412 return true;
413 };
414
415 // Find a range of numbers that have the same component in the longest
416 // coordinate direction:
417 auto range_start = sorted_vertices.begin();
418 while (range_start != sorted_vertices.end())
419 {
420 auto range_end = range_start + 1;
421 while (range_end != sorted_vertices.end() &&
422 std::abs(range_end->second[longest_coordinate_direction] -
423 range_start->second[longest_coordinate_direction]) <
424 tol)
425 ++range_end;
426
427 // preserve behavior with older versions of this function by replacing
428 // higher vertex numbers by lower vertex numbers
429 std::sort(range_start,
430 range_end,
431 [](const std::pair<unsigned int, Point<spacedim>> &a,
432 const std::pair<unsigned int, Point<spacedim>> &b) {
433 return a.first < b.first;
434 });
435
436 // Now de-duplicate [range_start, range_end)
437 //
438 // We have identified all points that are within a strip of width 'tol'
439 // in one coordinate direction. Now we need to figure out which of these
440 // are also close in other coordinate directions. If two are close, we
441 // can mark the second one for deletion.
442 for (auto reference = range_start; reference != range_end; ++reference)
443 {
444 if (reference->first != numbers::invalid_unsigned_int)
445 for (auto it = reference + 1; it != range_end; ++it)
446 {
447 if (within_tolerance(reference->second, it->second))
448 {
449 new_vertex_numbers[it->first] = reference->first;
450 // skip the replaced vertex in the future
452 }
453 }
454 }
455 range_start = range_end;
456 }
457
458 // now we got a renumbering list. simply renumber all vertices
459 // (non-duplicate vertices get renumbered to themselves, so nothing bad
460 // happens). after that, the duplicate vertices will be unused, so call
461 // delete_unused_vertices() to do that part of the job.
462 for (auto &cell : cells)
463 for (auto &vertex_index : cell.vertices)
464 vertex_index = new_vertex_numbers[vertex_index];
465 for (auto &quad : subcelldata.boundary_quads)
466 for (auto &vertex_index : quad.vertices)
467 vertex_index = new_vertex_numbers[vertex_index];
468 for (auto &line : subcelldata.boundary_lines)
469 for (auto &vertex_index : line.vertices)
470 vertex_index = new_vertex_numbers[vertex_index];
471
472 delete_unused_vertices(vertices, cells, subcelldata);
473 }
474
475
476
477 template <int dim>
478 void
480 const double tol)
481 {
482 if (vertices.empty())
483 return;
484
485 // 1) map point to local vertex index
486 std::map<Point<dim>, unsigned int, FloatingPointComparator<double>>
487 map_point_to_local_vertex_index{FloatingPointComparator<double>(tol)};
488
489 // 2) initialize map with existing points uniquely
490 for (unsigned int i = 0; i < vertices.size(); ++i)
491 map_point_to_local_vertex_index[vertices[i]] = i;
492
493 // no duplicate points are found
494 if (map_point_to_local_vertex_index.size() == vertices.size())
495 return;
496
497 // 3) remove duplicate entries from vertices
498 vertices.resize(map_point_to_local_vertex_index.size());
499 {
500 unsigned int j = 0;
501 for (const auto &p : map_point_to_local_vertex_index)
502 vertices[j++] = p.first;
503 }
504 }
505
506
507
508 template <int dim, int spacedim>
509 std::size_t
511 const std::vector<Point<spacedim>> &all_vertices,
512 std::vector<CellData<dim>> &cells)
513 {
514 // This function is presently only implemented for volumetric (codimension
515 // 0) elements.
516
517 if (dim == 1)
518 return 0;
519 if (dim == 2 && spacedim == 3)
521
522 std::size_t n_negative_cells = 0;
523 std::size_t cell_no = 0;
524 for (auto &cell : cells)
525 {
526 const ArrayView<const unsigned int> vertices(cell.vertices);
527 // Some pathologically twisted cells can have exactly zero measure but
528 // we can still fix them
529 if (GridTools::cell_measure(all_vertices, vertices) <= 0)
530 {
531 ++n_negative_cells;
532 const auto reference_cell =
534
535 if (reference_cell.is_hyper_cube())
536 {
537 if (dim == 2)
538 {
539 // flip the cell across the y = x line in 2d
540 std::swap(cell.vertices[1], cell.vertices[2]);
541 }
542 else if (dim == 3)
543 {
544 // swap the front and back faces in 3d
545 std::swap(cell.vertices[0], cell.vertices[2]);
546 std::swap(cell.vertices[1], cell.vertices[3]);
547 std::swap(cell.vertices[4], cell.vertices[6]);
548 std::swap(cell.vertices[5], cell.vertices[7]);
549 }
550 }
551 else if (reference_cell.is_simplex())
552 {
553 // By basic rules for computing determinants we can just swap
554 // two vertices to fix a negative volume. Arbitrarily pick the
555 // last two.
556 std::swap(cell.vertices[cell.vertices.size() - 2],
557 cell.vertices[cell.vertices.size() - 1]);
558 }
559 else if (reference_cell == ReferenceCells::Wedge)
560 {
561 // swap the two triangular faces
562 std::swap(cell.vertices[0], cell.vertices[3]);
563 std::swap(cell.vertices[1], cell.vertices[4]);
564 std::swap(cell.vertices[2], cell.vertices[5]);
565 }
566 else if (reference_cell == ReferenceCells::Pyramid)
567 {
568 // Try swapping two vertices in the base - perhaps things were
569 // read in the UCD (counter-clockwise) order instead of lexical
570 std::swap(cell.vertices[2], cell.vertices[3]);
571 }
572 else
573 {
575 }
576 // Check whether the resulting cell is now ok.
577 // If not, then the grid is seriously broken and
578 // we just give up.
579 AssertThrow(GridTools::cell_measure(all_vertices, vertices) > 0,
580 ExcGridHasInvalidCell(cell_no));
581 }
582 ++cell_no;
583 }
584 return n_negative_cells;
585 }
586
587
588 template <int dim, int spacedim>
589 void
591 const std::vector<Point<spacedim>> &all_vertices,
592 std::vector<CellData<dim>> &cells)
593 {
594 const std::size_t n_negative_cells =
595 invert_cells_with_negative_measure(all_vertices, cells);
596
597 // We assume that all cells of a grid have
598 // either positive or negative volumes but
599 // not both mixed. Although above reordering
600 // might work also on single cells, grids
601 // with both kind of cells are very likely to
602 // be broken. Check for this here.
603 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
605 std::string(
606 "This function assumes that either all cells have positive "
607 "volume, or that all cells have been specified in an "
608 "inverted vertex order so that their volume is negative. "
609 "(In the latter case, this class automatically inverts "
610 "every cell.) However, the mesh you have specified "
611 "appears to have both cells with positive and cells with "
612 "negative volume. You need to check your mesh which "
613 "cells these are and how they got there.\n"
614 "As a hint, of the total ") +
615 std::to_string(cells.size()) + " cells in the mesh, " +
616 std::to_string(n_negative_cells) +
617 " appear to have a negative volume."));
618 }
619
620
621
622 // Functions and classes for consistently_order_cells
623 namespace
624 {
630 struct CheapEdge
631 {
635 CheapEdge(const unsigned int v0, const unsigned int v1)
636 : v0(v0)
637 , v1(v1)
638 {}
639
644 bool
645 operator<(const CheapEdge &e) const
646 {
647 return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
648 }
649
650 private:
654 const unsigned int v0, v1;
655 };
656
657
666 template <int dim>
667 bool
668 is_consistent(const std::vector<CellData<dim>> &cells)
669 {
670 std::set<CheapEdge> edges;
671
672 for (typename std::vector<CellData<dim>>::const_iterator c =
673 cells.begin();
674 c != cells.end();
675 ++c)
676 {
677 // construct the edges in reverse order. for each of them,
678 // ensure that the reverse edge is not yet in the list of
679 // edges (return false if the reverse edge already *is* in
680 // the list) and then add the actual edge to it; std::set
681 // eliminates duplicates automatically
682 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
683 {
684 const CheapEdge reverse_edge(
687 if (edges.find(reverse_edge) != edges.end())
688 return false;
689
690
691 // ok, not. insert edge in correct order
692 const CheapEdge correct_edge(
695 edges.insert(correct_edge);
696 }
697 }
698
699 // no conflicts found, so return true
700 return true;
701 }
702
703
710 template <int dim>
711 struct ParallelEdges
712 {
718 static const unsigned int starter_edges[dim];
719
724 static const unsigned int n_other_parallel_edges = (1 << (dim - 1)) - 1;
725 static const unsigned int
728 };
729
730 template <>
731 const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2};
732
733 template <>
734 const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1},
735 {0},
736 {3},
737 {2}};
738
739 template <>
740 const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8};
741
742 template <>
743 const unsigned int ParallelEdges<3>::parallel_edges[12][3] = {
744 {1, 4, 5}, // line 0
745 {0, 4, 5}, // line 1
746 {3, 6, 7}, // line 2
747 {2, 6, 7}, // line 3
748 {0, 1, 5}, // line 4
749 {0, 1, 4}, // line 5
750 {2, 3, 7}, // line 6
751 {2, 3, 6}, // line 7
752 {9, 10, 11}, // line 8
753 {8, 10, 11}, // line 9
754 {8, 9, 11}, // line 10
755 {8, 9, 10} // line 11
756 };
757
758
763 struct AdjacentCell
764 {
768 AdjacentCell()
769 : cell_index(numbers::invalid_unsigned_int)
770 , edge_within_cell(numbers::invalid_unsigned_int)
771 {}
772
776 AdjacentCell(const unsigned int cell_index,
777 const unsigned int edge_within_cell)
780 {}
781
782
783 unsigned int cell_index;
784 unsigned int edge_within_cell;
785 };
786
787
788
789 template <int dim>
790 class AdjacentCells;
791
797 template <>
798 class AdjacentCells<2>
799 {
800 public:
805 using const_iterator = const AdjacentCell *;
806
815 void
816 push_back(const AdjacentCell &adjacent_cell)
817 {
819 adjacent_cells[0] = adjacent_cell;
820 else
821 {
825 adjacent_cells[1] = adjacent_cell;
826 }
827 }
828
829
834 const_iterator
835 begin() const
836 {
837 return adjacent_cells;
838 }
839
840
846 const_iterator
847 end() const
848 {
849 // check whether the current object stores zero, one, or two
850 // adjacent cells, and use this to point to the element past the
851 // last valid one
853 return adjacent_cells;
855 return adjacent_cells + 1;
856 else
857 return adjacent_cells + 2;
858 }
859
860 private:
867 AdjacentCell adjacent_cells[2];
868 };
869
870
871
879 template <>
880 class AdjacentCells<3> : public std::vector<AdjacentCell>
881 {};
882
883
893 template <int dim>
894 class Edge
895 {
896 public:
902 Edge(const CellData<dim> &cell, const unsigned int edge_number)
903 : orientation_status(not_oriented)
904 {
907
908 // copy vertices for this particular line
909 vertex_indices[0] =
910 cell
912 vertex_indices[1] =
913 cell
915
916 // bring them into standard orientation
917 if (vertex_indices[0] > vertex_indices[1])
918 std::swap(vertex_indices[0], vertex_indices[1]);
919 }
920
925 bool
926 operator<(const Edge<dim> &e) const
927 {
928 return ((vertex_indices[0] < e.vertex_indices[0]) ||
929 ((vertex_indices[0] == e.vertex_indices[0]) &&
930 (vertex_indices[1] < e.vertex_indices[1])));
931 }
932
936 bool
937 operator==(const Edge<dim> &e) const
938 {
939 return ((vertex_indices[0] == e.vertex_indices[0]) &&
940 (vertex_indices[1] == e.vertex_indices[1]));
941 }
942
947 unsigned int vertex_indices[2];
948
953 enum OrientationStatus
954 {
955 not_oriented,
956 forward,
957 backward
958 };
959
960 OrientationStatus orientation_status;
961
966 AdjacentCells<dim> adjacent_cells;
967 };
968
969
970
975 template <int dim>
976 struct Cell
977 {
983 Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
984 {
985 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
986 vertex_indices[i] = c.vertices[i];
987
988 // now for each of the edges of this cell, find the location inside the
989 // given edge_list array and store than index
990 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
991 {
992 const Edge<dim> e(c, l);
993 edge_indices[l] =
994 (std::lower_bound(edge_list.begin(), edge_list.end(), e) -
995 edge_list.begin());
996 Assert(edge_indices[l] < edge_list.size(), ExcInternalError());
997 Assert(edge_list[edge_indices[l]] == e, ExcInternalError());
998 }
999 }
1000
1005
1011 };
1012
1013
1014
1015 template <int dim>
1016 class EdgeDeltaSet;
1017
1027 template <>
1028 class EdgeDeltaSet<2>
1029 {
1030 public:
1034 using const_iterator = const unsigned int *;
1035
1040 EdgeDeltaSet()
1041 {
1043 }
1044
1045
1049 void
1050 clear()
1051 {
1053 }
1054
1059 void
1060 insert(const unsigned int edge_index)
1061 {
1063 edge_indices[0] = edge_index;
1064 else
1065 {
1068 edge_indices[1] = edge_index;
1069 }
1070 }
1071
1072
1076 const_iterator
1077 begin() const
1078 {
1079 return edge_indices;
1080 }
1081
1082
1086 const_iterator
1087 end() const
1088 {
1089 // check whether the current object stores zero, one, or two
1090 // indices, and use this to point to the element past the
1091 // last valid one
1093 return edge_indices;
1095 return edge_indices + 1;
1096 else
1097 return edge_indices + 2;
1098 }
1099
1100 private:
1104 unsigned int edge_indices[2];
1105 };
1106
1107
1108
1120 template <>
1121 class EdgeDeltaSet<3> : public std::set<unsigned int>
1122 {};
1123
1124
1125
1130 template <int dim>
1131 std::vector<Edge<dim>>
1132 build_edges(const std::vector<CellData<dim>> &cells)
1133 {
1134 // build the edge list for all cells. because each cell has
1135 // GeometryInfo<dim>::lines_per_cell edges, the total number
1136 // of edges is this many times the number of cells. of course
1137 // some of them will be duplicates, and we throw them out below
1138 std::vector<Edge<dim>> edge_list;
1139 edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
1140 for (unsigned int i = 0; i < cells.size(); ++i)
1141 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1142 edge_list.emplace_back(cells[i], l);
1143
1144 // next sort the edge list and then remove duplicates
1145 std::sort(edge_list.begin(), edge_list.end());
1146 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1147 edge_list.end());
1148
1149 return edge_list;
1150 }
1151
1152
1153
1158 template <int dim>
1159 std::vector<Cell<dim>>
1160 build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
1161 std::vector<Edge<dim>> &edges)
1162 {
1163 std::vector<Cell<dim>> cell_list;
1164 cell_list.reserve(cells.size());
1165 for (unsigned int i = 0; i < cells.size(); ++i)
1166 {
1167 // create our own data structure for the cells and let it
1168 // connect to the edges array
1169 cell_list.emplace_back(cells[i], edges);
1170
1171 // then also inform the edges that they are adjacent
1172 // to the current cell, and where within this cell
1173 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1174 edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back(
1175 AdjacentCell(i, l));
1176 }
1177 Assert(cell_list.size() == cells.size(), ExcInternalError());
1178
1179 return cell_list;
1180 }
1181
1182
1183
1188 template <int dim>
1189 unsigned int
1190 get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
1191 const std::vector<Edge<dim>> &edges,
1192 const unsigned int current_cell)
1193 {
1194 for (unsigned int c = current_cell; c < cells.size(); ++c)
1195 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1196 if (edges[cells[c].edge_indices[l]].orientation_status ==
1197 Edge<dim>::not_oriented)
1198 return c;
1199
1201 }
1202
1203
1204
1210 template <int dim>
1211 void
1212 orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> &cells,
1213 std::vector<Edge<dim>> &edges,
1214 const unsigned int cell,
1215 const unsigned int local_edge)
1216 {
1217 // choose the direction of the first edge. we have free choice
1218 // here and could simply choose "forward" if that's what pleases
1219 // us. however, for backward compatibility with the previous
1220 // implementation used till 2016, let us just choose the
1221 // direction so that it matches what we have in the given cell.
1222 //
1223 // in fact, in what can only be assumed to be a bug in the
1224 // original implementation, after orienting all edges, the code
1225 // that rotates the cells so that they match edge orientations
1226 // (see the rotate_cell() function below) rotated the cell two
1227 // more times by 90 degrees. this is ok -- it simply flips all
1228 // edge orientations, which leaves them valid. rather than do
1229 // the same in the current implementation, we can achieve the
1230 // same effect by modifying the rule above to choose the
1231 // direction of the starting edge of this parallel set
1232 // *opposite* to what it looks like in the current cell
1233 //
1234 // this bug only existed in the 2d implementation since there
1235 // were different implementations for 2d and 3d. consequently,
1236 // only replicate it for the 2d case and be "intuitive" in 3d.
1237 if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1239 local_edge, 0)])
1240 // orient initial edge *opposite* to the way it is in the cell
1241 // (see above for the reason)
1242 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1243 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1244 else
1245 {
1246 Assert(
1247 edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
1248 cells[cell].vertex_indices
1251 Assert(
1252 edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
1253 cells[cell].vertex_indices
1256
1257 // orient initial edge *opposite* to the way it is in the cell
1258 // (see above for the reason)
1259 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1260 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1261 }
1262
1263 // walk outward from the given edge as described in
1264 // the algorithm in the paper that documents all of
1265 // this
1266 //
1267 // note that in 2d, each of the Deltas can at most
1268 // contain two elements, whereas in 3d it can be arbitrarily many
1269 EdgeDeltaSet<dim> Delta_k;
1270 EdgeDeltaSet<dim> Delta_k_minus_1;
1271 Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]);
1272
1273 while (Delta_k_minus_1.begin() !=
1274 Delta_k_minus_1.end()) // while set is not empty
1275 {
1276 Delta_k.clear();
1277
1278 for (typename EdgeDeltaSet<dim>::const_iterator delta =
1279 Delta_k_minus_1.begin();
1280 delta != Delta_k_minus_1.end();
1281 ++delta)
1282 {
1283 Assert(edges[*delta].orientation_status !=
1284 Edge<dim>::not_oriented,
1286
1287 // now go through the cells adjacent to this edge
1288 for (typename AdjacentCells<dim>::const_iterator adjacent_cell =
1289 edges[*delta].adjacent_cells.begin();
1290 adjacent_cell != edges[*delta].adjacent_cells.end();
1291 ++adjacent_cell)
1292 {
1293 const unsigned int K = adjacent_cell->cell_index;
1294 const unsigned int delta_is_edge_in_K =
1295 adjacent_cell->edge_within_cell;
1296
1297 // figure out the direction of delta with respect to the cell
1298 // K (in the orientation in which the user has given it to us)
1299 const unsigned int first_edge_vertex =
1300 (edges[*delta].orientation_status == Edge<dim>::forward ?
1301 edges[*delta].vertex_indices[0] :
1302 edges[*delta].vertex_indices[1]);
1303 const unsigned int first_edge_vertex_in_K =
1304 cells[K]
1306 delta_is_edge_in_K, 0)];
1307 Assert(
1308 first_edge_vertex == first_edge_vertex_in_K ||
1309 first_edge_vertex ==
1311 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1313
1314 // now figure out which direction the each of the "opposite"
1315 // edges needs to be oriented into.
1316 for (unsigned int o_e = 0;
1317 o_e < ParallelEdges<dim>::n_other_parallel_edges;
1318 ++o_e)
1319 {
1320 // get the index of the opposite edge and select which its
1321 // first vertex needs to be based on how the current edge
1322 // is oriented in the current cell
1323 const unsigned int opposite_edge =
1324 cells[K].edge_indices[ParallelEdges<
1325 dim>::parallel_edges[delta_is_edge_in_K][o_e]];
1326 const unsigned int first_opposite_edge_vertex =
1327 cells[K].vertex_indices
1329 ParallelEdges<
1330 dim>::parallel_edges[delta_is_edge_in_K][o_e],
1331 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1332 1))];
1333
1334 // then determine the orientation of the edge based on
1335 // whether the vertex we want to be the edge's first
1336 // vertex is already the first vertex of the edge, or
1337 // whether it points in the opposite direction
1338 const typename Edge<dim>::OrientationStatus
1339 opposite_edge_orientation =
1340 (edges[opposite_edge].vertex_indices[0] ==
1341 first_opposite_edge_vertex ?
1342 Edge<dim>::forward :
1343 Edge<dim>::backward);
1344
1345 // see if the opposite edge (there is only one in 2d) has
1346 // already been oriented.
1347 if (edges[opposite_edge].orientation_status ==
1348 Edge<dim>::not_oriented)
1349 {
1350 // the opposite edge is not yet oriented. do orient it
1351 // and add it to Delta_k
1352 edges[opposite_edge].orientation_status =
1353 opposite_edge_orientation;
1354 Delta_k.insert(opposite_edge);
1355 }
1356 else
1357 {
1358 // this opposite edge has already been oriented. it
1359 // should be consistent with the current one in 2d,
1360 // while in 3d it may in fact be mis-oriented, and in
1361 // that case the mesh will not be orientable. indicate
1362 // this by throwing an exception that we can catch
1363 // further up; this has the advantage that we can
1364 // propagate through a couple of functions without
1365 // having to do error checking and without modifying
1366 // the 'cells' array that the user gave us
1367 if (dim == 2)
1368 {
1369 Assert(edges[opposite_edge].orientation_status ==
1370 opposite_edge_orientation,
1372 }
1373 else if (dim == 3)
1374 {
1375 if (edges[opposite_edge].orientation_status !=
1376 opposite_edge_orientation)
1377 throw ExcMeshNotOrientable();
1378 }
1379 else
1381 }
1382 }
1383 }
1384 }
1385
1386 // finally copy the new set to the previous one
1387 // (corresponding to increasing 'k' by one in the
1388 // algorithm)
1389 Delta_k_minus_1 = Delta_k;
1390 }
1391 }
1392
1393
1401 template <int dim>
1402 void
1403 rotate_cell(const std::vector<Cell<dim>> &cell_list,
1404 const std::vector<Edge<dim>> &edge_list,
1405 const unsigned int cell_index,
1406 std::vector<CellData<dim>> &raw_cells)
1407 {
1408 // find the first vertex of the cell. this is the vertex where dim edges
1409 // originate, so for each of the edges record which the starting vertex is
1410 unsigned int starting_vertex_of_edge[GeometryInfo<dim>::lines_per_cell];
1411 for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
1412 {
1413 Assert(edge_list[cell_list[cell_index].edge_indices[e]]
1414 .orientation_status != Edge<dim>::not_oriented,
1416 if (edge_list[cell_list[cell_index].edge_indices[e]]
1417 .orientation_status == Edge<dim>::forward)
1418 starting_vertex_of_edge[e] =
1419 edge_list[cell_list[cell_index].edge_indices[e]]
1420 .vertex_indices[0];
1421 else
1422 starting_vertex_of_edge[e] =
1423 edge_list[cell_list[cell_index].edge_indices[e]]
1424 .vertex_indices[1];
1425 }
1426
1427 // find the vertex number that appears dim times. this will then be
1428 // the vertex at which we want to locate the origin of the cell's
1429 // coordinate system (i.e., vertex 0)
1430 unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int;
1431 switch (dim)
1432 {
1433 case 2:
1434 {
1435 // in 2d, we can simply enumerate the possibilities where the
1436 // origin may be located because edges zero and one don't share
1437 // any vertices, and the same for edges two and three
1438 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1439 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1440 origin_vertex_of_cell = starting_vertex_of_edge[0];
1441 else if ((starting_vertex_of_edge[1] ==
1442 starting_vertex_of_edge[2]) ||
1443 (starting_vertex_of_edge[1] ==
1444 starting_vertex_of_edge[3]))
1445 origin_vertex_of_cell = starting_vertex_of_edge[1];
1446 else
1448
1449 break;
1450 }
1451
1452 case 3:
1453 {
1454 // one could probably do something similar in 3d, but that seems
1455 // more complicated than one wants to write down. just go
1456 // through the list of possible starting vertices and check
1457 for (origin_vertex_of_cell = 0;
1458 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1459 ++origin_vertex_of_cell)
1460 if (std::count(starting_vertex_of_edge,
1461 starting_vertex_of_edge +
1463 cell_list[cell_index]
1464 .vertex_indices[origin_vertex_of_cell]) == dim)
1465 break;
1466 Assert(origin_vertex_of_cell <
1469
1470 break;
1471 }
1472
1473 default:
1475 }
1476
1477 // now rotate raw_cells[cell_index] in such a way that its orientation
1478 // matches that of cell_list[cell_index]
1479 switch (dim)
1480 {
1481 case 2:
1482 {
1483 // in 2d, we can literally rotate the cell until its origin
1484 // matches the one that we have determined above should be
1485 // the origin vertex
1486 //
1487 // when doing a rotation, take into account the ordering of
1488 // vertices (not in clockwise or counter-clockwise sense)
1489 while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell)
1490 {
1491 const unsigned int tmp = raw_cells[cell_index].vertices[0];
1492 raw_cells[cell_index].vertices[0] =
1493 raw_cells[cell_index].vertices[1];
1494 raw_cells[cell_index].vertices[1] =
1495 raw_cells[cell_index].vertices[3];
1496 raw_cells[cell_index].vertices[3] =
1497 raw_cells[cell_index].vertices[2];
1498 raw_cells[cell_index].vertices[2] = tmp;
1499 }
1500 break;
1501 }
1502
1503 case 3:
1504 {
1505 // in 3d, the situation is a bit more complicated. from above, we
1506 // now know which vertex is at the origin (because 3 edges
1507 // originate from it), but that still leaves 3 possible rotations
1508 // of the cube. the important realization is that we can choose
1509 // any of them: in all 3 rotations, all edges originate from the
1510 // one vertex, and that fixes the directions of all 12 edges in
1511 // the cube because these 3 cover all 3 equivalence classes!
1512 // consequently, we can select an arbitrary one among the
1513 // permutations -- for example the following ones:
1514 static const unsigned int cube_permutations[8][8] = {
1515 {0, 1, 2, 3, 4, 5, 6, 7},
1516 {1, 5, 3, 7, 0, 4, 2, 6},
1517 {2, 6, 0, 4, 3, 7, 1, 5},
1518 {3, 2, 1, 0, 7, 6, 5, 4},
1519 {4, 0, 6, 2, 5, 1, 7, 3},
1520 {5, 4, 7, 6, 1, 0, 3, 2},
1521 {6, 7, 4, 5, 2, 3, 0, 1},
1522 {7, 3, 5, 1, 6, 2, 4, 0}};
1523
1524 unsigned int
1525 temp_vertex_indices[GeometryInfo<dim>::vertices_per_cell];
1526 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1527 temp_vertex_indices[v] =
1528 raw_cells[cell_index]
1529 .vertices[cube_permutations[origin_vertex_of_cell][v]];
1530 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1531 raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
1532
1533 break;
1534 }
1535
1536 default:
1537 {
1539 }
1540 }
1541 }
1542
1543
1549 template <int dim>
1550 void
1551 reorient(std::vector<CellData<dim>> &cells)
1552 {
1553 // first build the arrays that connect cells to edges and the other
1554 // way around
1555 std::vector<Edge<dim>> edge_list = build_edges(cells);
1556 std::vector<Cell<dim>> cell_list =
1557 build_cells_and_connect_edges(cells, edge_list);
1558
1559 // then loop over all cells and start orienting parallel edge sets
1560 // of cells that still have non-oriented edges
1561 unsigned int next_cell_with_unoriented_edge = 0;
1562 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1563 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1565 {
1566 // see which edge sets are still not oriented
1567 //
1568 // we do not need to look at each edge because if we orient edge
1569 // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
1570 // will be 3 other edges that are also oriented). there are only
1571 // dim independent sets of edges, so loop over these.
1572 //
1573 // we need to check whether each one of these starter edges may
1574 // already be oriented because the line (sheet) that connects
1575 // globally parallel edges may be self-intersecting in the
1576 // current cell
1577 for (unsigned int l = 0; l < dim; ++l)
1578 if (edge_list[cell_list[next_cell_with_unoriented_edge]
1579 .edge_indices[ParallelEdges<dim>::starter_edges[l]]]
1580 .orientation_status == Edge<dim>::not_oriented)
1581 orient_one_set_of_parallel_edges(
1582 cell_list,
1583 edge_list,
1584 next_cell_with_unoriented_edge,
1585 ParallelEdges<dim>::starter_edges[l]);
1586
1587 // ensure that we have really oriented all edges now, not just
1588 // the starter edges
1589 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1590 Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
1591 .edge_indices[l]]
1592 .orientation_status != Edge<dim>::not_oriented,
1594 }
1595
1596 // now that we have oriented all edges, we need to rotate cells
1597 // so that the edges point in the right direction with the now
1598 // rotated coordinate system
1599 for (unsigned int c = 0; c < cells.size(); ++c)
1600 rotate_cell(cell_list, edge_list, c, cells);
1601 }
1602
1603
1604 // overload of the function above for 1d -- there is nothing
1605 // to orient in that case
1606 void
1607 reorient(std::vector<CellData<1>> &)
1608 {}
1609 } // namespace
1610
1611
1612
1613 template <int dim>
1614 void
1616 {
1617 Assert(cells.size() != 0,
1618 ExcMessage(
1619 "List of elements to orient must have at least one cell"));
1620
1621 // there is nothing for us to do in 1d
1622 if (dim == 1)
1623 return;
1624
1625 // check if grids are already consistent. if so, do
1626 // nothing. if not, then do the reordering
1627 if (!is_consistent(cells))
1628 try
1629 {
1630 reorient(cells);
1631 }
1632 catch (const ExcMeshNotOrientable &)
1633 {
1634 // the mesh is not orientable. this is acceptable if we are in 3d,
1635 // as class Triangulation knows how to handle this, but it is
1636 // not in 2d; in that case, re-throw the exception
1637 if (dim < 3)
1638 throw;
1639 }
1640 }
1641
1642
1643
1644 template <int dim, int spacedim>
1645 std::map<unsigned int, Point<spacedim>>
1647 {
1648 std::map<unsigned int, Point<spacedim>> vertex_map;
1650 cell = tria.begin_active(),
1651 endc = tria.end();
1652 for (; cell != endc; ++cell)
1653 {
1654 for (const unsigned int i : cell->face_indices())
1655 {
1656 const typename Triangulation<dim, spacedim>::face_iterator &face =
1657 cell->face(i);
1658 if (face->at_boundary())
1659 {
1660 for (unsigned j = 0; j < face->n_vertices(); ++j)
1661 {
1662 const Point<spacedim> &vertex = face->vertex(j);
1663 const unsigned int vertex_index = face->vertex_index(j);
1664 vertex_map[vertex_index] = vertex;
1665 }
1666 }
1667 }
1668 }
1669 return vertex_map;
1670 }
1671
1672
1673
1674 template <int dim, int spacedim>
1675 void
1677 const bool isotropic,
1678 const unsigned int max_iterations)
1679 {
1680 unsigned int iter = 0;
1681 bool continue_refinement = true;
1682
1683 while (continue_refinement && (iter < max_iterations))
1684 {
1685 if (max_iterations != numbers::invalid_unsigned_int)
1686 ++iter;
1687 continue_refinement = false;
1688
1689 for (const auto &cell : tria.active_cell_iterators())
1690 for (const unsigned int j : cell->face_indices())
1691 if (cell->at_boundary(j) == false &&
1692 cell->neighbor(j)->has_children())
1693 {
1694 if (isotropic)
1695 {
1696 cell->set_refine_flag();
1697 continue_refinement = true;
1698 }
1699 else
1700 continue_refinement |= cell->flag_for_face_refinement(j);
1701 }
1702
1704 }
1705 }
1706
1707
1708
1709 template <int dim, int spacedim>
1710 void
1712 const double max_ratio,
1713 const unsigned int max_iterations)
1714 {
1715 unsigned int iter = 0;
1716 bool continue_refinement = true;
1717
1718 while (continue_refinement && (iter < max_iterations))
1719 {
1720 ++iter;
1721 continue_refinement = false;
1722 for (const auto &cell : tria.active_cell_iterators())
1723 {
1724 std::pair<unsigned int, double> info =
1725 GridTools::get_longest_direction<dim, spacedim>(cell);
1726 if (info.second > max_ratio)
1727 {
1728 cell->set_refine_flag(
1729 RefinementCase<dim>::cut_axis(info.first));
1730 continue_refinement = true;
1731 }
1732 }
1734 }
1735 }
1736
1737
1738
1739 template <int dim, int spacedim>
1740 std::map<unsigned int, Point<spacedim>>
1742 const Mapping<dim, spacedim> &mapping)
1743 {
1744 std::map<unsigned int, Point<spacedim>> result;
1745 for (const auto &cell : container.active_cell_iterators())
1746 {
1747 if (!cell->is_artificial())
1748 {
1749 const auto vs = mapping.get_vertices(cell);
1750 for (unsigned int i = 0; i < vs.size(); ++i)
1751 result[cell->vertex_index(i)] = vs[i];
1752 }
1753 }
1754 return result;
1755 }
1756
1757
1758
1759 template <int dim, int spacedim>
1760 std::vector<
1761 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1763 {
1764 std::vector<
1765 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1766 vertex_to_cell_map(triangulation.n_vertices());
1768 cell = triangulation.begin_active(),
1769 endc = triangulation.end();
1770 for (; cell != endc; ++cell)
1771 for (const unsigned int i : cell->vertex_indices())
1772 vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
1773
1774 // Check if mesh has hanging nodes. Do this only locally to
1775 // prevent communication and possible deadlock.
1776 if (triangulation.Triangulation<dim, spacedim>::has_hanging_nodes())
1777 {
1778 Assert(triangulation.all_reference_cells_are_hyper_cube(),
1780
1781 // Take care of hanging nodes
1782 cell = triangulation.begin_active();
1783 for (; cell != endc; ++cell)
1784 {
1785 for (const unsigned int i : cell->face_indices())
1786 {
1787 if ((cell->at_boundary(i) == false) &&
1788 (cell->neighbor(i)->is_active()))
1789 {
1791 adjacent_cell = cell->neighbor(i);
1792 for (unsigned int j = 0; j < cell->face(i)->n_vertices();
1793 ++j)
1794 vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(
1795 adjacent_cell);
1796 }
1797 }
1798
1799 // in 3d also loop over the edges
1800 if (dim == 3)
1801 {
1802 for (unsigned int i = 0; i < cell->n_lines(); ++i)
1803 if (cell->line(i)->has_children())
1804 // the only place where this vertex could have been
1805 // hiding is on the mid-edge point of the edge we
1806 // are looking at
1807 vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)]
1808 .insert(cell);
1809 }
1810 }
1811 }
1812
1813 return vertex_to_cell_map;
1814 }
1815
1816
1817
1818 template <int dim, int spacedim>
1819 void
1822 DynamicSparsityPattern &cell_connectivity)
1823 {
1824 cell_connectivity.reinit(triangulation.n_active_cells(),
1825 triangulation.n_active_cells());
1826
1827 // loop over all cells and their neighbors to build the sparsity
1828 // pattern. note that it's a bit hard to enter all the connections when a
1829 // neighbor has children since we would need to find out which of its
1830 // children is adjacent to the current cell. this problem can be omitted
1831 // if we only do something if the neighbor has no children -- in that case
1832 // it is either on the same or a coarser level than we are. in return, we
1833 // have to add entries in both directions for both cells
1834 for (const auto &cell : triangulation.active_cell_iterators())
1835 {
1836 const unsigned int index = cell->active_cell_index();
1837 cell_connectivity.add(index, index);
1838 for (auto f : cell->face_indices())
1839 if ((cell->at_boundary(f) == false) &&
1840 (cell->neighbor(f)->has_children() == false))
1841 {
1842 const unsigned int other_index =
1843 cell->neighbor(f)->active_cell_index();
1844 cell_connectivity.add(index, other_index);
1845 cell_connectivity.add(other_index, index);
1846 }
1847 }
1848 }
1849
1850
1851
1852 template <int dim, int spacedim>
1853 void
1856 DynamicSparsityPattern &cell_connectivity)
1857 {
1858 // The choice of 16 or fewer neighbors here is based on empirical
1859 // measurements.
1860 //
1861 // Vertices in a structured hexahedral mesh have 8 adjacent cells. In a
1862 // structured tetrahedral mesh, about 98% of vertices have 16 neighbors or
1863 // fewer. Similarly, in an unstructured tetrahedral mesh, if we count the
1864 // number of neighbors each vertex has we obtain the following distribution:
1865 //
1866 // 3, 1
1867 // 4, 728
1868 // 5, 4084
1869 // 6, 7614
1870 // 7, 17329
1871 // 8, 31145
1872 // 9, 46698
1873 // 10, 64193
1874 // 11, 68269
1875 // 12, 63574
1876 // 13, 57016
1877 // 14, 50476
1878 // 15, 41886
1879 // 16, 31820
1880 // 17, 21269
1881 // 18, 12217
1882 // 19, 6072
1883 // 20, 2527
1884 // 21, 825
1885 // 22, 262
1886 // 23, 61
1887 // 24, 12
1888 // 26, 1
1889 //
1890 // so about 86% of vertices have 16 neighbors or fewer. Hence, we picked 16
1891 // neighbors here to cover most cases without allocation.
1892 std::vector<boost::container::small_vector<unsigned int, 16>>
1893 vertex_to_cell(triangulation.n_vertices());
1894 for (const auto &cell : triangulation.active_cell_iterators())
1895 {
1896 for (const unsigned int v : cell->vertex_indices())
1897 vertex_to_cell[cell->vertex_index(v)].push_back(
1898 cell->active_cell_index());
1899 }
1900
1901 cell_connectivity.reinit(triangulation.n_active_cells(),
1902 triangulation.n_active_cells());
1903 std::vector<types::global_dof_index> neighbors;
1904 for (const auto &cell : triangulation.active_cell_iterators())
1905 {
1906 neighbors.clear();
1907 for (const unsigned int v : cell->vertex_indices())
1908 neighbors.insert(neighbors.end(),
1909 vertex_to_cell[cell->vertex_index(v)].begin(),
1910 vertex_to_cell[cell->vertex_index(v)].end());
1911 std::sort(neighbors.begin(), neighbors.end());
1912 cell_connectivity.add_entries(cell->active_cell_index(),
1913 neighbors.begin(),
1914 std::unique(neighbors.begin(),
1915 neighbors.end()),
1916 true);
1917 }
1918 }
1919
1920
1921 template <int dim, int spacedim>
1922 void
1925 const unsigned int level,
1926 DynamicSparsityPattern &cell_connectivity)
1927 {
1928 std::vector<boost::container::small_vector<unsigned int, 16>>
1929 vertex_to_cell(triangulation.n_vertices());
1931 triangulation.begin(level);
1932 cell != triangulation.end(level);
1933 ++cell)
1934 {
1935 for (const unsigned int v : cell->vertex_indices())
1936 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
1937 }
1938
1939 cell_connectivity.reinit(triangulation.n_cells(level),
1940 triangulation.n_cells(level));
1941 std::vector<types::global_dof_index> neighbors;
1942 for (const auto &cell : triangulation.cell_iterators_on_level(level))
1943 {
1944 neighbors.clear();
1945 for (const unsigned int v : cell->vertex_indices())
1946 neighbors.insert(neighbors.end(),
1947 vertex_to_cell[cell->vertex_index(v)].begin(),
1948 vertex_to_cell[cell->vertex_index(v)].end());
1949 std::sort(neighbors.begin(), neighbors.end());
1950 cell_connectivity.add_entries(cell->index(),
1951 neighbors.begin(),
1952 std::unique(neighbors.begin(),
1953 neighbors.end()),
1954 true);
1955 }
1956 }
1957} /* namespace GridTools */
1958
1959// explicit instantiations
1960#include "grid_tools_topology.inst"
1961
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
std::size_t size() const
Definition array_view.h:684
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:318
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition point.h:111
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
unsigned int level
Definition grid_out.cc:4632
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)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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:470
const types::boundary_id internal_face_boundary_id
Definition types.h:312
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::manifold_id flat_manifold_id
Definition types.h:325
::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
types::manifold_id manifold_id
types::material_id material_id
types::boundary_id boundary_id
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
bool check_consistency(const unsigned int dim) const
std::vector< CellData< 1 > > boundary_lines
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)