1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16
17 #include <deal.II/base/timer.h>
18 #include <deal.II/base/utilities.h>
19
20 #include <deal.II/grid/grid_reordering.h>
21 #include <deal.II/grid/grid_tools.h>
22
23 #include <algorithm>
24 #include <fstream>
25 #include <functional>
26 #include <iostream>
27 #include <set>
28
29 DEAL_II_NAMESPACE_OPEN
30
31
32 namespace
33 {
34 /**
35 * A simple data structure denoting an edge, i.e., the ordered pair
36 * of its vertex indices. This is only used in the is_consistent()
37 * function.
38 */
39 struct CheapEdge
40 {
41 /**
42 * Construct an edge from the global indices of its two vertices.
43 */
CheapEdge__anon3d046c610111::CheapEdge44 CheapEdge(const unsigned int v0, const unsigned int v1)
45 : v0(v0)
46 , v1(v1)
47 {}
48
49 /**
50 * Comparison operator for edges. It compares based on the
51 * lexicographic ordering of the two vertex indices.
52 */
53 bool
operator <__anon3d046c610111::CheapEdge54 operator<(const CheapEdge &e) const
55 {
56 return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
57 }
58
59 private:
60 /**
61 * The global indices of the vertices that define the edge.
62 */
63 const unsigned int v0, v1;
64 };
65
66
67 /**
68 * A function that determines whether the edges in a mesh are
69 * already consistently oriented. It does so by adding all edges
70 * of all cells into a set (which automatically eliminates
71 * duplicates) but before that checks whether the reverse edge is
72 * already in the set -- which would imply that a neighboring cell
73 * is inconsistently oriented.
74 */
75 template <int dim>
76 bool
is_consistent(const std::vector<CellData<dim>> & cells)77 is_consistent(const std::vector<CellData<dim>> &cells)
78 {
79 std::set<CheapEdge> edges;
80
81 for (typename std::vector<CellData<dim>>::const_iterator c = cells.begin();
82 c != cells.end();
83 ++c)
84 {
85 // construct the edges in reverse order. for each of them,
86 // ensure that the reverse edge is not yet in the list of
87 // edges (return false if the reverse edge already *is* in
88 // the list) and then add the actual edge to it; std::set
89 // eliminates duplicates automatically
90 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
91 {
92 const CheapEdge reverse_edge(
93 c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)],
94 c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 0)]);
95 if (edges.find(reverse_edge) != edges.end())
96 return false;
97
98
99 // ok, not. insert edge in correct order
100 const CheapEdge correct_edge(
101 c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 0)],
102 c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)]);
103 edges.insert(correct_edge);
104 }
105 }
106
107 // no conflicts found, so return true
108 return true;
109 }
110
111
112 /**
113 * A structure that describes some properties of parallel edges
114 * such as what starter edges are (i.e., representative elements
115 * of the sets of parallel edges within a cell) and what the set
116 * of parallel edges to each edge is.
117 */
118 template <int dim>
119 struct ParallelEdges
120 {
121 /**
122 * An array that contains the indices of dim edges that can
123 * serve as (arbitrarily chosen) starting points for the
124 * dim sets of parallel edges within each cell.
125 */
126 static const unsigned int starter_edges[dim];
127
128 /**
129 * Number and indices of all of those edges parallel to each of the
130 * edges in a cell.
131 */
132 static const unsigned int n_other_parallel_edges = (1 << (dim - 1)) - 1;
133 static const unsigned int parallel_edges[GeometryInfo<dim>::lines_per_cell]
134 [n_other_parallel_edges];
135 };
136
137 template <>
138 const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2};
139
140 template <>
141 const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1},
142 {0},
143 {3},
144 {2}};
145
146 template <>
147 const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8};
148
149 template <>
150 const unsigned int ParallelEdges<3>::parallel_edges[12][3] = {
151 {1, 4, 5}, // line 0
152 {0, 4, 5}, // line 1
153 {3, 6, 7}, // line 2
154 {2, 6, 7}, // line 3
155 {0, 1, 5}, // line 4
156 {0, 1, 4}, // line 5
157 {2, 3, 7}, // line 6
158 {2, 3, 6}, // line 7
159 {9, 10, 11}, // line 8
160 {8, 10, 11}, // line 9
161 {8, 9, 11}, // line 10
162 {8, 9, 10} // line 11
163 };
164
165
166 /**
167 * A structure that store the index of a cell and, crucially, how a
168 * given edge relates to this cell.
169 */
170 struct AdjacentCell
171 {
172 /**
173 * Default constructor. Initialize the fields with invalid values.
174 */
AdjacentCell__anon3d046c610111::AdjacentCell175 AdjacentCell()
176 : cell_index(numbers::invalid_unsigned_int)
177 , edge_within_cell(numbers::invalid_unsigned_int)
178 {}
179
180 /**
181 * Constructor. Initialize the fields with the given values.
182 */
AdjacentCell__anon3d046c610111::AdjacentCell183 AdjacentCell(const unsigned int cell_index,
184 const unsigned int edge_within_cell)
185 : cell_index(cell_index)
186 , edge_within_cell(edge_within_cell)
187 {}
188
189
190 unsigned int cell_index;
191 unsigned int edge_within_cell;
192 };
193
194
195
196 template <int dim>
197 class AdjacentCells;
198
199 /**
200 * A class that represents all of the cells adjacent to a given edge.
201 * This class corresponds to the 2d case where each edge has at most
202 * two adjacent cells.
203 */
204 template <>
205 class AdjacentCells<2>
206 {
207 public:
208 /**
209 * An iterator that allows iterating over all cells adjacent
210 * to the edge represented by the current object.
211 */
212 using const_iterator = const AdjacentCell *;
213
214 /**
215 * Add the given cell to the collection of cells adjacent to
216 * the edge this object corresponds to. Since we are covering
217 * the 2d case, the set of adjacent cells currently
218 * represented by this object must have either zero or
219 * one element already, since we can not add more than two
220 * adjacent cells for each edge.
221 */
222 void
push_back(const AdjacentCell & adjacent_cell)223 push_back(const AdjacentCell &adjacent_cell)
224 {
225 if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int)
226 adjacent_cells[0] = adjacent_cell;
227 else
228 {
229 Assert(adjacent_cells[1].cell_index == numbers::invalid_unsigned_int,
230 ExcInternalError());
231 adjacent_cells[1] = adjacent_cell;
232 }
233 }
234
235
236 /**
237 * Return an iterator to the first valid cell stored as adjacent to the
238 * edge represented by the current object.
239 */
240 const_iterator
begin() const241 begin() const
242 {
243 return adjacent_cells;
244 }
245
246
247 /**
248 * Return an iterator to the element past the last valid cell stored
249 * as adjacent to the edge represented by the current object.
250 * @return
251 */
252 const_iterator
end() const253 end() const
254 {
255 // check whether the current object stores zero, one, or two
256 // adjacent cells, and use this to point to the element past the
257 // last valid one
258 if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int)
259 return adjacent_cells;
260 else if (adjacent_cells[1].cell_index == numbers::invalid_unsigned_int)
261 return adjacent_cells + 1;
262 else
263 return adjacent_cells + 2;
264 }
265
266 private:
267 /**
268 * References to the (at most) two cells that are adjacent to
269 * the edge this object corresponds to. Unused elements are
270 * default-initialized and have invalid values; in particular,
271 * their cell_index field equals numbers::invalid_unsigned_int.
272 */
273 AdjacentCell adjacent_cells[2];
274 };
275
276
277
278 /**
279 * A class that represents all of the cells adjacent to a given edge.
280 * This class corresponds to the 3d case where each edge can have an
281 * arbitrary number of adjacent cells. We represent this as a
282 * std::vector<AdjacentCell>, from which class the current one is
283 * derived and from which it inherits all of its member functions.
284 */
285 template <>
286 class AdjacentCells<3> : public std::vector<AdjacentCell>
287 {};
288
289
290 /**
291 * A class that describes all of the relevant properties of an
292 * edge. For the purpose of what we do here, that includes the
293 * indices of the two vertices, and the indices of the adjacent
294 * cells (together with a description *where* in each of the
295 * adjacent cells the edge is located). It also includes the
296 * (global) direction of the edge: either from the first vertex to
297 * the second, the other way around, or so far undetermined.
298 */
299 template <int dim>
300 class Edge
301 {
302 public:
303 /**
304 * Constructor. Create the edge based on the information given
305 * in @p cell, and selecting the edge with number @p edge_number
306 * within this cell. Initialize the edge as unoriented.
307 */
Edge(const CellData<dim> & cell,const unsigned int edge_number)308 Edge(const CellData<dim> &cell, const unsigned int edge_number)
309 : orientation_status(not_oriented)
310 {
311 Assert(edge_number < GeometryInfo<dim>::lines_per_cell,
312 ExcInternalError());
313
314 // copy vertices for this particular line
315 vertex_indices[0] =
316 cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 0)];
317 vertex_indices[1] =
318 cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 1)];
319
320 // bring them into standard orientation
321 if (vertex_indices[0] > vertex_indices[1])
322 std::swap(vertex_indices[0], vertex_indices[1]);
323 }
324
325 /**
326 * Comparison operator for edges. It compares based on the
327 * lexicographic ordering of the two vertex indices.
328 */
329 bool
operator <(const Edge<dim> & e) const330 operator<(const Edge<dim> &e) const
331 {
332 return ((vertex_indices[0] < e.vertex_indices[0]) ||
333 ((vertex_indices[0] == e.vertex_indices[0]) &&
334 (vertex_indices[1] < e.vertex_indices[1])));
335 }
336
337 /**
338 * Compare two edges for equality based on their vertex indices.
339 */
340 bool
operator ==(const Edge<dim> & e) const341 operator==(const Edge<dim> &e) const
342 {
343 return ((vertex_indices[0] == e.vertex_indices[0]) &&
344 (vertex_indices[1] == e.vertex_indices[1]));
345 }
346
347 /**
348 * The global indices of the two vertices that bound this edge. These
349 * will be ordered so that the first index is less than the second.
350 */
351 unsigned int vertex_indices[2];
352
353 /**
354 * An enum that indicates the direction of this edge with
355 * regard to the two vertices that bound it.
356 */
357 enum OrientationStatus
358 {
359 not_oriented,
360 forward,
361 backward
362 };
363
364 OrientationStatus orientation_status;
365
366 /**
367 * Store the set of cells adjacent to this edge (these cells then
368 * also store *where* in the cell the edge is located).
369 */
370 AdjacentCells<dim> adjacent_cells;
371 };
372
373
374
375 /**
376 * A data structure that represents a cell with all of its vertices
377 * and edges.
378 */
379 template <int dim>
380 struct Cell
381 {
382 /**
383 * Construct a Cell object from a CellData object. Also take a
384 * (sorted) list of edges and to point the edges of the current
385 * object into this list of edges.
386 */
Cell__anon3d046c610111::Cell387 Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
388 {
389 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
390 vertex_indices[i] = c.vertices[i];
391
392 // now for each of the edges of this cell, find the location inside the
393 // given edge_list array and store than index
394 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
395 {
396 const Edge<dim> e(c, l);
397 edge_indices[l] =
398 (std::lower_bound(edge_list.begin(), edge_list.end(), e) -
399 edge_list.begin());
400 Assert(edge_indices[l] < edge_list.size(), ExcInternalError());
401 Assert(edge_list[edge_indices[l]] == e, ExcInternalError())
402 }
403 }
404
405 /**
406 * A list of global indices for the vertices that bound this cell.
407 */
408 unsigned int vertex_indices[GeometryInfo<dim>::vertices_per_cell];
409
410 /**
411 * A list of indices into the 'edge_list' array passed to the constructor
412 * for the edges of the current cell.
413 */
414 unsigned int edge_indices[GeometryInfo<dim>::lines_per_cell];
415 };
416
417
418
419 template <int dim>
420 class EdgeDeltaSet;
421
422 /**
423 * A class that represents by how much the set of parallel edges
424 * grows in each step. In the graph orientation paper, this set is
425 * called $\Delta_k$, thus the name.
426 *
427 * In 2d, this set can only include zero, one, or two elements.
428 * Consequently, the appropriate data structure is one in which
429 * we store at most 2 elements in a fixed sized data structure.
430 */
431 template <>
432 class EdgeDeltaSet<2>
433 {
434 public:
435 /**
436 * Iterator type for the elements of the set.
437 */
438 using const_iterator = const unsigned int *;
439
440 /**
441 * Default constructor. Initialize both slots as unused, corresponding
442 * to an empty set.
443 */
EdgeDeltaSet()444 EdgeDeltaSet()
445 {
446 edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
447 }
448
449
450 /**
451 * Delete the elements of the set by marking both slots as unused.
452 */
453 void
clear()454 clear()
455 {
456 edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
457 }
458
459 /**
460 * Insert one element into the set. This will fail if the set already
461 * has two elements.
462 */
463 void
insert(const unsigned int edge_index)464 insert(const unsigned int edge_index)
465 {
466 if (edge_indices[0] == numbers::invalid_unsigned_int)
467 edge_indices[0] = edge_index;
468 else
469 {
470 Assert(edge_indices[1] == numbers::invalid_unsigned_int,
471 ExcInternalError());
472 edge_indices[1] = edge_index;
473 }
474 }
475
476
477 /**
478 * Return an iterator pointing to the first element of the set.
479 */
480 const_iterator
begin() const481 begin() const
482 {
483 return edge_indices;
484 }
485
486
487 /**
488 * Return an iterator pointing to the element past the last used one.
489 */
490 const_iterator
end() const491 end() const
492 {
493 // check whether the current object stores zero, one, or two
494 // indices, and use this to point to the element past the
495 // last valid one
496 if (edge_indices[0] == numbers::invalid_unsigned_int)
497 return edge_indices;
498 else if (edge_indices[1] == numbers::invalid_unsigned_int)
499 return edge_indices + 1;
500 else
501 return edge_indices + 2;
502 }
503
504 private:
505 /**
506 * Storage space to store the indices of at most two edges.
507 */
508 unsigned int edge_indices[2];
509 };
510
511
512
513 /**
514 * A class that represents by how much the set of parallel edges
515 * grows in each step. In the graph orientation paper, this set is
516 * called $\Delta_k$, thus the name.
517 *
518 * In 3d, this set can have arbitrarily many elements, unlike the
519 * 2d case specialized above. Consequently, we simply represent
520 * the data structure with a std::set. Class derivation ensures
521 * that we simply inherit all of the member functions of the
522 * base class.
523 */
524 template <>
525 class EdgeDeltaSet<3> : public std::set<unsigned int>
526 {};
527
528
529
530 /**
531 * From a list of cells, build a sorted vector that contains all of the edges
532 * that exist in the mesh.
533 */
534 template <int dim>
535 std::vector<Edge<dim>>
build_edges(const std::vector<CellData<dim>> & cells)536 build_edges(const std::vector<CellData<dim>> &cells)
537 {
538 // build the edge list for all cells. because each cell has
539 // GeometryInfo<dim>::lines_per_cell edges, the total number
540 // of edges is this many times the number of cells. of course
541 // some of them will be duplicates, and we throw them out below
542 std::vector<Edge<dim>> edge_list;
543 edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
544 for (unsigned int i = 0; i < cells.size(); ++i)
545 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
546 edge_list.emplace_back(cells[i], l);
547
548 // next sort the edge list and then remove duplicates
549 std::sort(edge_list.begin(), edge_list.end());
550 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
551 edge_list.end());
552
553 return edge_list;
554 }
555
556
557
558 /**
559 * Build the cell list. Update the edge array to let edges know
560 * which cells are adjacent to them.
561 */
562 template <int dim>
563 std::vector<Cell<dim>>
build_cells_and_connect_edges(const std::vector<CellData<dim>> & cells,std::vector<Edge<dim>> & edges)564 build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
565 std::vector<Edge<dim>> & edges)
566 {
567 std::vector<Cell<dim>> cell_list;
568 cell_list.reserve(cells.size());
569 for (unsigned int i = 0; i < cells.size(); ++i)
570 {
571 // create our own data structure for the cells and let it
572 // connect to the edges array
573 cell_list.emplace_back(cells[i], edges);
574
575 // then also inform the edges that they are adjacent
576 // to the current cell, and where within this cell
577 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
578 edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back(
579 AdjacentCell(i, l));
580 }
581 Assert(cell_list.size() == cells.size(), ExcInternalError());
582
583 return cell_list;
584 }
585
586
587
588 /**
589 * Return the index within 'cells' of the first cell that has at least one
590 * edge that is not yet oriented.
591 */
592 template <int dim>
593 unsigned int
get_next_unoriented_cell(const std::vector<Cell<dim>> & cells,const std::vector<Edge<dim>> & edges,const unsigned int current_cell)594 get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
595 const std::vector<Edge<dim>> &edges,
596 const unsigned int current_cell)
597 {
598 for (unsigned int c = current_cell; c < cells.size(); ++c)
599 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
600 if (edges[cells[c].edge_indices[l]].orientation_status ==
601 Edge<dim>::not_oriented)
602 return c;
603
604 return numbers::invalid_unsigned_int;
605 }
606
607
608
609 /**
610 * Given a set of cells and edges, orient all edges that are
611 * (global) parallel to the one identified by the @p cell and
612 * within it the one with index @p local_edge.
613 */
614 template <int dim>
615 void
orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> & cells,std::vector<Edge<dim>> & edges,const unsigned int cell,const unsigned int local_edge)616 orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> &cells,
617 std::vector<Edge<dim>> & edges,
618 const unsigned int cell,
619 const unsigned int local_edge)
620 {
621 // choose the direction of the first edge. we have free choice
622 // here and could simply choose "forward" if that's what pleases
623 // us. however, for backward compatibility with the previous
624 // implementation used till 2016, let us just choose the
625 // direction so that it matches what we have in the given cell.
626 //
627 // in fact, in what can only be assumed to be a bug in the
628 // original implementation, after orienting all edges, the code
629 // that rotates the cells so that they match edge orientations
630 // (see the rotate_cell() function below) rotated the cell two
631 // more times by 90 degrees. this is ok -- it simply flips all
632 // edge orientations, which leaves them valid. rather than do
633 // the same in the current implementation, we can achieve the
634 // same effect by modifying the rule above to choose the
635 // direction of the starting edge of this parallel set
636 // *opposite* to what it looks like in the current cell
637 //
638 // this bug only existed in the 2d implementation since there
639 // were different implementations for 2d and 3d. consequently,
640 // only replicate it for the 2d case and be "intuitive" in 3d.
641 if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
642 cells[cell].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
643 local_edge, 0)])
644 // orient initial edge *opposite* to the way it is in the cell
645 // (see above for the reason)
646 edges[cells[cell].edge_indices[local_edge]].orientation_status =
647 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
648 else
649 {
650 Assert(
651 edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
652 cells[cell].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
653 local_edge, 1)],
654 ExcInternalError());
655 Assert(
656 edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
657 cells[cell].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
658 local_edge, 0)],
659 ExcInternalError());
660
661 // orient initial edge *opposite* to the way it is in the cell
662 // (see above for the reason)
663 edges[cells[cell].edge_indices[local_edge]].orientation_status =
664 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
665 }
666
667 // walk outward from the given edge as described in
668 // the algorithm in the paper that documents all of
669 // this
670 //
671 // note that in 2d, each of the Deltas can at most
672 // contain two elements, whereas in 3d it can be arbitrarily many
673 EdgeDeltaSet<dim> Delta_k;
674 EdgeDeltaSet<dim> Delta_k_minus_1;
675 Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]);
676
677 while (Delta_k_minus_1.begin() !=
678 Delta_k_minus_1.end()) // while set is not empty
679 {
680 Delta_k.clear();
681
682 for (typename EdgeDeltaSet<dim>::const_iterator delta =
683 Delta_k_minus_1.begin();
684 delta != Delta_k_minus_1.end();
685 ++delta)
686 {
687 Assert(edges[*delta].orientation_status != Edge<dim>::not_oriented,
688 ExcInternalError());
689
690 // now go through the cells adjacent to this edge
691 for (typename AdjacentCells<dim>::const_iterator adjacent_cell =
692 edges[*delta].adjacent_cells.begin();
693 adjacent_cell != edges[*delta].adjacent_cells.end();
694 ++adjacent_cell)
695 {
696 const unsigned int K = adjacent_cell->cell_index;
697 const unsigned int delta_is_edge_in_K =
698 adjacent_cell->edge_within_cell;
699
700 // figure out the direction of delta with respect to the cell K
701 // (in the orientation in which the user has given it to us)
702 const unsigned int first_edge_vertex =
703 (edges[*delta].orientation_status == Edge<dim>::forward ?
704 edges[*delta].vertex_indices[0] :
705 edges[*delta].vertex_indices[1]);
706 const unsigned int first_edge_vertex_in_K =
707 cells[K]
708 .vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
709 delta_is_edge_in_K, 0)];
710 Assert(
711 first_edge_vertex == first_edge_vertex_in_K ||
712 first_edge_vertex ==
713 cells[K].vertex_indices[GeometryInfo<
714 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
715 ExcInternalError());
716
717 // now figure out which direction the each of the "opposite"
718 // edges needs to be oriented into.
719 for (unsigned int o_e = 0;
720 o_e < ParallelEdges<dim>::n_other_parallel_edges;
721 ++o_e)
722 {
723 // get the index of the opposite edge and select which its
724 // first vertex needs to be based on how the current edge is
725 // oriented in the current cell
726 const unsigned int opposite_edge =
727 cells[K].edge_indices[ParallelEdges<
728 dim>::parallel_edges[delta_is_edge_in_K][o_e]];
729 const unsigned int first_opposite_edge_vertex =
730 cells[K].vertex_indices
731 [GeometryInfo<dim>::line_to_cell_vertices(
732 ParallelEdges<dim>::parallel_edges[delta_is_edge_in_K]
733 [o_e],
734 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
735 1))];
736
737 // then determine the orientation of the edge based on
738 // whether the vertex we want to be the edge's first
739 // vertex is already the first vertex of the edge, or
740 // whether it points in the opposite direction
741 const typename Edge<dim>::OrientationStatus
742 opposite_edge_orientation =
743 (edges[opposite_edge].vertex_indices[0] ==
744 first_opposite_edge_vertex ?
745 Edge<dim>::forward :
746 Edge<dim>::backward);
747
748 // see if the opposite edge (there is only one in 2d) has
749 // already been oriented.
750 if (edges[opposite_edge].orientation_status ==
751 Edge<dim>::not_oriented)
752 {
753 // the opposite edge is not yet oriented. do orient it
754 // and add it to Delta_k
755 edges[opposite_edge].orientation_status =
756 opposite_edge_orientation;
757 Delta_k.insert(opposite_edge);
758 }
759 else
760 {
761 // this opposite edge has already been oriented. it
762 // should be consistent with the current one in 2d,
763 // while in 3d it may in fact be mis-oriented, and in
764 // that case the mesh will not be orientable. indicate
765 // this by throwing an exception that we can catch
766 // further up; this has the advantage that we can
767 // propagate through a couple of functions without
768 // having to do error checking and without modifying the
769 // 'cells' array that the user gave us
770 if (dim == 2)
771 {
772 Assert(edges[opposite_edge].orientation_status ==
773 opposite_edge_orientation,
774 ExcMeshNotOrientable());
775 }
776 else if (dim == 3)
777 {
778 if (edges[opposite_edge].orientation_status !=
779 opposite_edge_orientation)
780 throw ExcMeshNotOrientable();
781 }
782 else
783 Assert(false, ExcNotImplemented());
784 }
785 }
786 }
787 }
788
789 // finally copy the new set to the previous one
790 // (corresponding to increasing 'k' by one in the
791 // algorithm)
792 Delta_k_minus_1 = Delta_k;
793 }
794 }
795
796
797 /**
798 * Given data structures @p cell_list and @p edge_list, where
799 * all edges are already oriented, rotate the cell with
800 * index @p cell_index in such a way that its local coordinate
801 * system matches the ones of the adjacent edges. Store the
802 * rotated order of vertices in <code>raw_cells[cell_index]</code>.
803 */
804 template <int dim>
805 void
rotate_cell(const std::vector<Cell<dim>> & cell_list,const std::vector<Edge<dim>> & edge_list,const unsigned int cell_index,std::vector<CellData<dim>> & raw_cells)806 rotate_cell(const std::vector<Cell<dim>> &cell_list,
807 const std::vector<Edge<dim>> &edge_list,
808 const unsigned int cell_index,
809 std::vector<CellData<dim>> & raw_cells)
810 {
811 // find the first vertex of the cell. this is the vertex where dim edges
812 // originate, so for each of the edges record which the starting vertex is
813 unsigned int starting_vertex_of_edge[GeometryInfo<dim>::lines_per_cell];
814 for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
815 {
816 Assert(edge_list[cell_list[cell_index].edge_indices[e]]
817 .orientation_status != Edge<dim>::not_oriented,
818 ExcInternalError());
819 if (edge_list[cell_list[cell_index].edge_indices[e]]
820 .orientation_status == Edge<dim>::forward)
821 starting_vertex_of_edge[e] =
822 edge_list[cell_list[cell_index].edge_indices[e]].vertex_indices[0];
823 else
824 starting_vertex_of_edge[e] =
825 edge_list[cell_list[cell_index].edge_indices[e]].vertex_indices[1];
826 }
827
828 // find the vertex number that appears dim times. this will then be
829 // the vertex at which we want to locate the origin of the cell's
830 // coordinate system (i.e., vertex 0)
831 unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int;
832 switch (dim)
833 {
834 case 2:
835 {
836 // in 2d, we can simply enumerate the possibilities where the
837 // origin may be located because edges zero and one don't share
838 // any vertices, and the same for edges two and three
839 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
840 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
841 origin_vertex_of_cell = starting_vertex_of_edge[0];
842 else if ((starting_vertex_of_edge[1] ==
843 starting_vertex_of_edge[2]) ||
844 (starting_vertex_of_edge[1] == starting_vertex_of_edge[3]))
845 origin_vertex_of_cell = starting_vertex_of_edge[1];
846 else
847 Assert(false, ExcInternalError());
848
849 break;
850 }
851
852 case 3:
853 {
854 // one could probably do something similar in 3d, but that seems
855 // more complicated than one wants to write down. just go
856 // through the list of possible starting vertices and check
857 for (origin_vertex_of_cell = 0;
858 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
859 ++origin_vertex_of_cell)
860 if (std::count(starting_vertex_of_edge,
861 starting_vertex_of_edge +
862 GeometryInfo<dim>::lines_per_cell,
863 cell_list[cell_index]
864 .vertex_indices[origin_vertex_of_cell]) == dim)
865 break;
866 Assert(origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell,
867 ExcInternalError());
868
869 break;
870 }
871
872 default:
873 Assert(false, ExcNotImplemented());
874 }
875
876 // now rotate raw_cells[cell_index] in such a way that its orientation
877 // matches that of cell_list[cell_index]
878 switch (dim)
879 {
880 case 2:
881 {
882 // in 2d, we can literally rotate the cell until its origin
883 // matches the one that we have determined above should be
884 // the origin vertex
885 //
886 // when doing a rotation, take into account the ordering of
887 // vertices (not in clockwise or counter-clockwise sense)
888 while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell)
889 {
890 const unsigned int tmp = raw_cells[cell_index].vertices[0];
891 raw_cells[cell_index].vertices[0] =
892 raw_cells[cell_index].vertices[1];
893 raw_cells[cell_index].vertices[1] =
894 raw_cells[cell_index].vertices[3];
895 raw_cells[cell_index].vertices[3] =
896 raw_cells[cell_index].vertices[2];
897 raw_cells[cell_index].vertices[2] = tmp;
898 }
899 break;
900 }
901
902 case 3:
903 {
904 // in 3d, the situation is a bit more complicated. from above, we
905 // now know which vertex is at the origin (because 3 edges originate
906 // from it), but that still leaves 3 possible rotations of the cube.
907 // the important realization is that we can choose any of them:
908 // in all 3 rotations, all edges originate from the one vertex,
909 // and that fixes the directions of all 12 edges in the cube because
910 // these 3 cover all 3 equivalence classes! consequently, we can
911 // select an arbitrary one among the permutations -- for
912 // example the following ones:
913 static const unsigned int cube_permutations[8][8] = {
914 {0, 1, 2, 3, 4, 5, 6, 7},
915 {1, 5, 3, 7, 0, 4, 2, 6},
916 {2, 6, 0, 4, 3, 7, 1, 5},
917 {3, 2, 1, 0, 7, 6, 5, 4},
918 {4, 0, 6, 2, 5, 1, 7, 3},
919 {5, 4, 7, 6, 1, 0, 3, 2},
920 {6, 7, 4, 5, 2, 3, 0, 1},
921 {7, 3, 5, 1, 6, 2, 4, 0}};
922
923 unsigned int
924 temp_vertex_indices[GeometryInfo<dim>::vertices_per_cell];
925 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
926 temp_vertex_indices[v] =
927 raw_cells[cell_index]
928 .vertices[cube_permutations[origin_vertex_of_cell][v]];
929 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
930 raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
931
932 break;
933 }
934
935 default:
936 {
937 Assert(false, ExcNotImplemented());
938 }
939 }
940 }
941
942
943 /**
944 * Given a set of cells, find globally unique edge orientations
945 * and then rotate cells so that the coordinate system of the cell
946 * coincides with the coordinate systems of the adjacent edges.
947 */
948 template <int dim>
949 void
reorient(std::vector<CellData<dim>> & cells)950 reorient(std::vector<CellData<dim>> &cells)
951 {
952 // first build the arrays that connect cells to edges and the other
953 // way around
954 std::vector<Edge<dim>> edge_list = build_edges(cells);
955 std::vector<Cell<dim>> cell_list =
956 build_cells_and_connect_edges(cells, edge_list);
957
958 // then loop over all cells and start orienting parallel edge sets
959 // of cells that still have non-oriented edges
960 unsigned int next_cell_with_unoriented_edge = 0;
961 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
962 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
963 numbers::invalid_unsigned_int)
964 {
965 // see which edge sets are still not oriented
966 //
967 // we do not need to look at each edge because if we orient edge
968 // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
969 // will be 3 other edges that are also oriented). there are only
970 // dim independent sets of edges, so loop over these.
971 //
972 // we need to check whether each one of these starter edges may
973 // already be oriented because the line (sheet) that connects
974 // globally parallel edges may be self-intersecting in the
975 // current cell
976 for (unsigned int l = 0; l < dim; ++l)
977 if (edge_list[cell_list[next_cell_with_unoriented_edge]
978 .edge_indices[ParallelEdges<dim>::starter_edges[l]]]
979 .orientation_status == Edge<dim>::not_oriented)
980 orient_one_set_of_parallel_edges(
981 cell_list,
982 edge_list,
983 next_cell_with_unoriented_edge,
984 ParallelEdges<dim>::starter_edges[l]);
985
986 // ensure that we have really oriented all edges now, not just
987 // the starter edges
988 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
989 Assert(
990 edge_list[cell_list[next_cell_with_unoriented_edge].edge_indices[l]]
991 .orientation_status != Edge<dim>::not_oriented,
992 ExcInternalError());
993 }
994
995 // now that we have oriented all edges, we need to rotate cells
996 // so that the edges point in the right direction with the now
997 // rotated coordinate system
998 for (unsigned int c = 0; c < cells.size(); ++c)
999 rotate_cell(cell_list, edge_list, c, cells);
1000 }
1001
1002
1003 // overload of the function above for 1d -- there is nothing
1004 // to orient in that case
reorient(std::vector<CellData<1>> &)1005 void reorient(std::vector<CellData<1>> &)
1006 {}
1007 } // namespace
1008
1009
1010 // anonymous namespace for internal helper functions
1011 namespace
1012 {
1013 /**
1014 * A set of functions that
1015 * reorder the data from the
1016 * "current" to the "classic"
1017 * format of vertex numbering of
1018 * cells and faces. These functions
1019 * do the reordering of their
1020 * arguments in-place.
1021 */
reorder_new_to_old_style(std::vector<CellData<1>> &)1022 void reorder_new_to_old_style(std::vector<CellData<1>> &)
1023 {}
1024
1025
reorder_new_to_old_style(std::vector<CellData<2>> & cells)1026 void reorder_new_to_old_style(std::vector<CellData<2>> &cells)
1027 {
1028 for (auto &cell : cells)
1029 std::swap(cell.vertices[2], cell.vertices[3]);
1030 }
1031
1032
reorder_new_to_old_style(std::vector<CellData<3>> & cells)1033 void reorder_new_to_old_style(std::vector<CellData<3>> &cells)
1034 {
1035 unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
1036 for (auto &cell : cells)
1037 {
1038 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
1039 tmp[i] = cell.vertices[i];
1040 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
1041 cell.vertices[i] = tmp[GeometryInfo<3>::ucd_to_deal[i]];
1042 }
1043 }
1044
1045
1046 /**
1047 * And now also in the opposite direction.
1048 */
reorder_old_to_new_style(std::vector<CellData<1>> &)1049 void reorder_old_to_new_style(std::vector<CellData<1>> &)
1050 {}
1051
1052
reorder_old_to_new_style(std::vector<CellData<2>> & cells)1053 void reorder_old_to_new_style(std::vector<CellData<2>> &cells)
1054 {
1055 // just invert the permutation:
1056 reorder_new_to_old_style(cells);
1057 }
1058
1059
reorder_old_to_new_style(std::vector<CellData<3>> & cells)1060 void reorder_old_to_new_style(std::vector<CellData<3>> &cells)
1061 {
1062 // undo the ordering above
1063 unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
1064 for (auto &cell : cells)
1065 {
1066 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
1067 tmp[i] = cell.vertices[i];
1068 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
1069 cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
1070 }
1071 }
1072 } // namespace
1073
1074
1075
1076 template <int dim, int spacedim>
1077 void
reorder_cells(std::vector<CellData<dim>> & cells,const bool use_new_style_ordering)1078 GridReordering<dim, spacedim>::reorder_cells(std::vector<CellData<dim>> &cells,
1079 const bool use_new_style_ordering)
1080 {
1081 Assert(cells.size() != 0,
1082 ExcMessage("List of elements to orient must have at least one cell"));
1083
1084 // there is nothing for us to do in 1d
1085 if (dim == 1)
1086 return;
1087
1088 // if necessary, convert to new-style format
1089 if (use_new_style_ordering == false)
1090 reorder_old_to_new_style(cells);
1091
1092 // check if grids are already consistent. if so, do
1093 // nothing. if not, then do the reordering
1094 if (!is_consistent(cells))
1095 try
1096 {
1097 reorient(cells);
1098 }
1099 catch (const ExcMeshNotOrientable &)
1100 {
1101 // the mesh is not orientable. this is acceptable if we are in 3d,
1102 // as class Triangulation knows how to handle this, but it is
1103 // not in 2d; in that case, re-throw the exception
1104 if (dim < 3)
1105 throw;
1106 }
1107
1108 // and convert back if necessary
1109 if (use_new_style_ordering == false)
1110 reorder_new_to_old_style(cells);
1111 }
1112
1113
1114
1115 template <>
1116 void
invert_all_cells_of_negative_grid(const std::vector<Point<1>> &,std::vector<CellData<1>> &)1117 GridReordering<1>::invert_all_cells_of_negative_grid(
1118 const std::vector<Point<1>> &,
1119 std::vector<CellData<1>> &)
1120 {
1121 // nothing to be done in 1d
1122 }
1123
1124
1125
1126 template <>
1127 void
invert_all_cells_of_negative_grid(const std::vector<Point<2>> &,std::vector<CellData<1>> &)1128 GridReordering<1, 2>::invert_all_cells_of_negative_grid(
1129 const std::vector<Point<2>> &,
1130 std::vector<CellData<1>> &)
1131 {
1132 // nothing to be done in 1d
1133 }
1134
1135
1136
1137 template <>
1138 void
invert_all_cells_of_negative_grid(const std::vector<Point<3>> &,std::vector<CellData<1>> &)1139 GridReordering<1, 3>::invert_all_cells_of_negative_grid(
1140 const std::vector<Point<3>> &,
1141 std::vector<CellData<1>> &)
1142 {
1143 // nothing to be done in 1d
1144 }
1145
1146
1147 template <>
1148 void
invert_all_cells_of_negative_grid(const std::vector<Point<2>> & all_vertices,std::vector<CellData<2>> & cells)1149 GridReordering<2>::invert_all_cells_of_negative_grid(
1150 const std::vector<Point<2>> &all_vertices,
1151 std::vector<CellData<2>> & cells)
1152 {
1153 unsigned int vertices_lex[GeometryInfo<2>::vertices_per_cell];
1154 unsigned int n_negative_cells = 0;
1155 for (auto &cell : cells)
1156 {
1157 // GridTools::cell_measure
1158 // requires the vertices to be
1159 // in lexicographic ordering
1160 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1161 vertices_lex[GeometryInfo<2>::ucd_to_deal[i]] = cell.vertices[i];
1162 if (GridTools::cell_measure<2>(all_vertices, vertices_lex) < 0)
1163 {
1164 ++n_negative_cells;
1165 std::swap(cell.vertices[1], cell.vertices[3]);
1166
1167 // Check whether the resulting cell is now ok.
1168 // If not, then the grid is seriously broken and
1169 // we just give up.
1170 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1171 vertices_lex[GeometryInfo<2>::ucd_to_deal[i]] = cell.vertices[i];
1172 AssertThrow(GridTools::cell_measure<2>(all_vertices, vertices_lex) >
1173 0,
1174 ExcInternalError());
1175 }
1176 }
1177
1178 // We assume that all cells of a grid have
1179 // either positive or negative volumes but
1180 // not both mixed. Although above reordering
1181 // might work also on single cells, grids
1182 // with both kind of cells are very likely to
1183 // be broken. Check for this here.
1184 AssertThrow(
1185 n_negative_cells == 0 || n_negative_cells == cells.size(),
1186 ExcMessage(
1187 std::string("This class assumes that either all cells have positive "
1188 "volume, or that all cells have been specified in an "
1189 "inverted vertex order so that their volume is negative. "
1190 "(In the latter case, this class automatically inverts "
1191 "every cell.) However, the mesh you have specified "
1192 "appears to have both cells with positive and cells with "
1193 "negative volume. You need to check your mesh which "
1194 "cells these are and how they got there.\n"
1195 "As a hint, of the total ") +
1196 std::to_string(cells.size()) + " cells in the mesh, " +
1197 std::to_string(n_negative_cells) + " appear to have a negative volume."));
1198 }
1199
1200
1201
1202 template <>
1203 void
invert_all_cells_of_negative_grid(const std::vector<Point<3>> &,std::vector<CellData<2>> &)1204 GridReordering<2, 3>::invert_all_cells_of_negative_grid(
1205 const std::vector<Point<3>> &,
1206 std::vector<CellData<2>> &)
1207 {
1208 Assert(false, ExcNotImplemented());
1209 }
1210
1211
1212
1213 template <>
1214 void
invert_all_cells_of_negative_grid(const std::vector<Point<3>> & all_vertices,std::vector<CellData<3>> & cells)1215 GridReordering<3>::invert_all_cells_of_negative_grid(
1216 const std::vector<Point<3>> &all_vertices,
1217 std::vector<CellData<3>> & cells)
1218 {
1219 unsigned int vertices_lex[GeometryInfo<3>::vertices_per_cell];
1220 unsigned int n_negative_cells = 0;
1221 for (auto &cell : cells)
1222 {
1223 // GridTools::cell_measure
1224 // requires the vertices to be
1225 // in lexicographic ordering
1226 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
1227 vertices_lex[GeometryInfo<3>::ucd_to_deal[i]] = cell.vertices[i];
1228 if (GridTools::cell_measure<3>(all_vertices, vertices_lex) < 0)
1229 {
1230 ++n_negative_cells;
1231 // reorder vertices: swap front and back face
1232 for (unsigned int i = 0; i < 4; ++i)
1233 std::swap(cell.vertices[i], cell.vertices[i + 4]);
1234
1235 // Check whether the resulting cell is now ok.
1236 // If not, then the grid is seriously broken and
1237 // we just give up.
1238 for (const unsigned int i : GeometryInfo<3>::vertex_indices())
1239 vertices_lex[GeometryInfo<3>::ucd_to_deal[i]] = cell.vertices[i];
1240 AssertThrow(GridTools::cell_measure<3>(all_vertices, vertices_lex) >
1241 0,
1242 ExcInternalError());
1243 }
1244 }
1245
1246 // We assume that all cells of a
1247 // grid have either positive or
1248 // negative volumes but not both
1249 // mixed. Although above reordering
1250 // might work also on single cells,
1251 // grids with both kind of cells
1252 // are very likely to be
1253 // broken. Check for this here.
1254 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
1255 ExcMessage("While sorting the cells that will be passed for "
1256 "creating a Triangulation object, deal.II found that "
1257 "some but not all cells have a negative volume. (If "
1258 "all cells had a negative volume, they would simply "
1259 "all have been inverted.) This usually happens in "
1260 "hand-generated meshes if one accidentally uses an "
1261 "incorrect convention for ordering the vertices in "
1262 "one or more cells; in that case, you may want to "
1263 "double check that you specified the vertex indices "
1264 "in their correct order. If you are reading a mesh "
1265 "that was created by a mesh generator, then this "
1266 "exception indicates that some of the cells created "
1267 "are so badly distorted that their volume becomes "
1268 "negative; this commonly occurs at complex geometric "
1269 "features, and you may see if the problem can be "
1270 "fixed by playing with the parameters that control "
1271 "mesh properties in your mesh generator, such as "
1272 "the number of cells, the mesh density, etc."));
1273 }
1274
1275
1276
1277 /* ------------------------ explicit instantiations ------------------- */
1278 template class GridReordering<1, 1>;
1279 template class GridReordering<1, 2>;
1280 template class GridReordering<1, 3>;
1281 template class GridReordering<2, 2>;
1282 template class GridReordering<2, 3>;
1283 template class GridReordering<3, 3>;
1284
1285 DEAL_II_NAMESPACE_CLOSE
1286