Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
grid_reordering.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2019 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 {
39  struct CheapEdge
40  {
44  CheapEdge(const unsigned int v0, const unsigned int v1)
45  : v0(v0)
46  , v1(v1)
47  {}
48 
53  bool
54  operator<(const CheapEdge &e) const
55  {
56  return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
57  }
58 
59  private:
63  const unsigned int v0, v1;
64  };
65 
66 
75  template <int dim>
76  bool
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(
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 
118  template <int dim>
119  struct ParallelEdges
120  {
126  static const unsigned int starter_edges[dim];
127 
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 
170  struct AdjacentCell
171  {
175  AdjacentCell()
176  : cell_index(numbers::invalid_unsigned_int)
177  , edge_within_cell(numbers::invalid_unsigned_int)
178  {}
179 
183  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 
204  template <>
205  class AdjacentCells<2>
206  {
207  public:
212  using const_iterator = const AdjacentCell *;
213 
222  void
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 
240  const_iterator
241  begin() const
242  {
243  return adjacent_cells;
244  }
245 
246 
252  const_iterator
253  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:
273  AdjacentCell adjacent_cells[2];
274  };
275 
276 
277 
285  template <>
286  class AdjacentCells<3> : public std::vector<AdjacentCell>
287  {};
288 
289 
299  template <int dim>
300  class Edge
301  {
302  public:
308  Edge(const CellData<dim> &cell, const unsigned int edge_number)
309  : orientation_status(not_oriented)
310  {
312  ExcInternalError());
313 
314  // copy vertices for this particular line
315  vertex_indices[0] =
317  vertex_indices[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 
329  bool
330  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 
340  bool
341  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 
351  unsigned int vertex_indices[2];
352 
357  enum OrientationStatus
358  {
359  not_oriented,
360  forward,
361  backward
362  };
363 
364  OrientationStatus orientation_status;
365 
370  AdjacentCells<dim> adjacent_cells;
371  };
372 
373 
374 
379  template <int dim>
380  struct Cell
381  {
387  Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
388  {
389  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
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 
408  unsigned int vertex_indices[GeometryInfo<dim>::vertices_per_cell];
409 
414  unsigned int edge_indices[GeometryInfo<dim>::lines_per_cell];
415  };
416 
417 
418 
419  template <int dim>
420  class EdgeDeltaSet;
421 
431  template <>
432  class EdgeDeltaSet<2>
433  {
434  public:
438  using const_iterator = const unsigned int *;
439 
444  EdgeDeltaSet()
445  {
446  edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
447  }
448 
449 
453  void
454  clear()
455  {
456  edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
457  }
458 
463  void
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 
480  const_iterator
481  begin() const
482  {
483  return edge_indices;
484  }
485 
486 
490  const_iterator
491  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:
508  unsigned int edge_indices[2];
509  };
510 
511 
512 
524  template <>
525  class EdgeDeltaSet<3> : public std::set<unsigned int>
526  {};
527 
528 
529 
534  template <int dim>
535  std::vector<Edge<dim>>
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 
562  template <int dim>
563  std::vector<Cell<dim>>
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 
592  template <int dim>
593  unsigned int
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 
605  }
606 
607 
608 
614  template <int dim>
615  void
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]
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
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,
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 
804  template <int dim>
805  void
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 +
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 (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
926  ++v)
927  temp_vertex_indices[v] =
928  raw_cells[cell_index]
929  .vertices[cube_permutations[origin_vertex_of_cell][v]];
930  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
931  ++v)
932  raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
933 
934  break;
935  }
936 
937  default:
938  {
939  Assert(false, ExcNotImplemented());
940  }
941  }
942  }
943 
944 
950  template <int dim>
951  void
952  reorient(std::vector<CellData<dim>> &cells)
953  {
954  // first build the arrays that connect cells to edges and the other
955  // way around
956  std::vector<Edge<dim>> edge_list = build_edges(cells);
957  std::vector<Cell<dim>> cell_list =
958  build_cells_and_connect_edges(cells, edge_list);
959 
960  // then loop over all cells and start orienting parallel edge sets
961  // of cells that still have non-oriented edges
962  unsigned int next_cell_with_unoriented_edge = 0;
963  while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
964  cell_list, edge_list, next_cell_with_unoriented_edge)) !=
966  {
967  // see which edge sets are still not oriented
968  //
969  // we do not need to look at each edge because if we orient edge
970  // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
971  // will be 3 other edges that are also oriented). there are only
972  // dim independent sets of edges, so loop over these.
973  //
974  // we need to check whether each one of these starter edges may
975  // already be oriented because the line (sheet) that connects
976  // globally parallel edges may be self-intersecting in the
977  // current cell
978  for (unsigned int l = 0; l < dim; ++l)
979  if (edge_list[cell_list[next_cell_with_unoriented_edge]
980  .edge_indices[ParallelEdges<dim>::starter_edges[l]]]
981  .orientation_status == Edge<dim>::not_oriented)
982  orient_one_set_of_parallel_edges(
983  cell_list,
984  edge_list,
985  next_cell_with_unoriented_edge,
986  ParallelEdges<dim>::starter_edges[l]);
987 
988  // ensure that we have really oriented all edges now, not just
989  // the starter edges
990  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
991  Assert(
992  edge_list[cell_list[next_cell_with_unoriented_edge].edge_indices[l]]
993  .orientation_status != Edge<dim>::not_oriented,
994  ExcInternalError());
995  }
996 
997  // now that we have oriented all edges, we need to rotate cells
998  // so that the edges point in the right direction with the now
999  // rotated coordinate system
1000  for (unsigned int c = 0; c < cells.size(); ++c)
1001  rotate_cell(cell_list, edge_list, c, cells);
1002  }
1003 
1004 
1005  // overload of the function above for 1d -- there is nothing
1006  // to orient in that case
1007  void reorient(std::vector<CellData<1>> &)
1008  {}
1009 } // namespace
1010 
1011 
1012 // anonymous namespace for internal helper functions
1013 namespace
1014 {
1024  void reorder_new_to_old_style(std::vector<CellData<1>> &)
1025  {}
1026 
1027 
1028  void reorder_new_to_old_style(std::vector<CellData<2>> &cells)
1029  {
1030  for (auto &cell : cells)
1031  std::swap(cell.vertices[2], cell.vertices[3]);
1032  }
1033 
1034 
1035  void reorder_new_to_old_style(std::vector<CellData<3>> &cells)
1036  {
1037  unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
1038  for (auto &cell : cells)
1039  {
1040  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell; ++i)
1041  tmp[i] = cell.vertices[i];
1042  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell; ++i)
1043  cell.vertices[i] = tmp[GeometryInfo<3>::ucd_to_deal[i]];
1044  }
1045  }
1046 
1047 
1051  void reorder_old_to_new_style(std::vector<CellData<1>> &)
1052  {}
1053 
1054 
1055  void reorder_old_to_new_style(std::vector<CellData<2>> &cells)
1056  {
1057  // just invert the permutation:
1058  reorder_new_to_old_style(cells);
1059  }
1060 
1061 
1062  void reorder_old_to_new_style(std::vector<CellData<3>> &cells)
1063  {
1064  // undo the ordering above
1065  unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
1066  for (auto &cell : cells)
1067  {
1068  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell; ++i)
1069  tmp[i] = cell.vertices[i];
1070  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell; ++i)
1071  cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
1072  }
1073  }
1074 } // namespace
1075 
1076 
1077 
1078 template <int dim, int spacedim>
1079 void
1081  const bool use_new_style_ordering)
1082 {
1083  Assert(cells.size() != 0,
1084  ExcMessage("List of elements to orient must have at least one cell"));
1085 
1086  // there is nothing for us to do in 1d
1087  if (dim == 1)
1088  return;
1089 
1090  // if necessary, convert to new-style format
1091  if (use_new_style_ordering == false)
1092  reorder_old_to_new_style(cells);
1093 
1094  // check if grids are already consistent. if so, do
1095  // nothing. if not, then do the reordering
1096  if (!is_consistent(cells))
1097  try
1098  {
1099  reorient(cells);
1100  }
1101  catch (const ExcMeshNotOrientable &)
1102  {
1103  // the mesh is not orientable. this is acceptable if we are in 3d,
1104  // as class Triangulation knows how to handle this, but it is
1105  // not in 2d; in that case, re-throw the exception
1106  if (dim < 3)
1107  throw;
1108  }
1109 
1110  // and convert back if necessary
1111  if (use_new_style_ordering == false)
1112  reorder_new_to_old_style(cells);
1113 }
1114 
1115 
1116 
1117 template <>
1118 void
1120  const std::vector<Point<1>> &,
1121  std::vector<CellData<1>> &)
1122 {
1123  // nothing to be done in 1d
1124 }
1125 
1126 
1127 
1128 template <>
1129 void
1131  const std::vector<Point<2>> &,
1132  std::vector<CellData<1>> &)
1133 {
1134  // nothing to be done in 1d
1135 }
1136 
1137 
1138 
1139 template <>
1140 void
1142  const std::vector<Point<3>> &,
1143  std::vector<CellData<1>> &)
1144 {
1145  // nothing to be done in 1d
1146 }
1147 
1148 
1149 template <>
1150 void
1152  const std::vector<Point<2>> &all_vertices,
1153  std::vector<CellData<2>> & cells)
1154 {
1155  unsigned int vertices_lex[GeometryInfo<2>::vertices_per_cell];
1156  unsigned int n_negative_cells = 0;
1157  for (auto &cell : cells)
1158  {
1159  // GridTools::cell_measure
1160  // requires the vertices to be
1161  // in lexicographic ordering
1162  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
1163  vertices_lex[GeometryInfo<2>::ucd_to_deal[i]] = cell.vertices[i];
1164  if (GridTools::cell_measure<2>(all_vertices, vertices_lex) < 0)
1165  {
1166  ++n_negative_cells;
1167  std::swap(cell.vertices[1], cell.vertices[3]);
1168 
1169  // check whether the
1170  // resulting cell is now ok.
1171  // if not, then the grid is
1172  // seriously broken and
1173  // should be sticked into the
1174  // bin
1175  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
1176  vertices_lex[GeometryInfo<2>::ucd_to_deal[i]] = cell.vertices[i];
1177  AssertThrow(GridTools::cell_measure<2>(all_vertices, vertices_lex) >
1178  0,
1179  ExcInternalError());
1180  }
1181  }
1182 
1183  // We assume that all cells of a grid have
1184  // either positive or negative volumes but
1185  // not both mixed. Although above reordering
1186  // might work also on single cells, grids
1187  // with both kind of cells are very likely to
1188  // be broken. Check for this here.
1189  AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
1190  ExcMessage(
1191  std::string(
1192  "This class assumes that either all cells have positive "
1193  "volume, or that all cells have been specified in an "
1194  "inverted vertex order so that their volume is negative. "
1195  "(In the latter case, this class automatically inverts "
1196  "every cell.) However, the mesh you have specified "
1197  "appears to have both cells with positive and cells with "
1198  "negative volume. You need to check your mesh which "
1199  "cells these are and how they got there.\n"
1200  "As a hint, of the total ") +
1201  Utilities::to_string(cells.size()) + " cells in the mesh, " +
1202  Utilities::to_string(n_negative_cells) +
1203  " appear to have a negative volume."));
1204 }
1205 
1206 
1207 
1208 template <>
1209 void
1211  const std::vector<Point<3>> &,
1212  std::vector<CellData<2>> &)
1213 {
1214  Assert(false, ExcNotImplemented());
1215 }
1216 
1217 
1218 
1219 template <>
1220 void
1222  const std::vector<Point<3>> &all_vertices,
1223  std::vector<CellData<3>> & cells)
1224 {
1225  unsigned int vertices_lex[GeometryInfo<3>::vertices_per_cell];
1226  unsigned int n_negative_cells = 0;
1227  for (auto &cell : cells)
1228  {
1229  // GridTools::cell_measure
1230  // requires the vertices to be
1231  // in lexicographic ordering
1232  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell; ++i)
1233  vertices_lex[GeometryInfo<3>::ucd_to_deal[i]] = cell.vertices[i];
1234  if (GridTools::cell_measure<3>(all_vertices, vertices_lex) < 0)
1235  {
1236  ++n_negative_cells;
1237  // reorder vertices: swap front and back face
1238  for (unsigned int i = 0; i < 4; ++i)
1239  std::swap(cell.vertices[i], cell.vertices[i + 4]);
1240 
1241  // check whether the
1242  // resulting cell is now ok.
1243  // if not, then the grid is
1244  // seriously broken and
1245  // should be sticked into the
1246  // bin
1247  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell; ++i)
1248  vertices_lex[GeometryInfo<3>::ucd_to_deal[i]] = cell.vertices[i];
1249  AssertThrow(GridTools::cell_measure<3>(all_vertices, vertices_lex) >
1250  0,
1251  ExcInternalError());
1252  }
1253  }
1254 
1255  // We assume that all cells of a
1256  // grid have either positive or
1257  // negative volumes but not both
1258  // mixed. Although above reordering
1259  // might work also on single cells,
1260  // grids with both kind of cells
1261  // are very likely to be
1262  // broken. Check for this here.
1263  AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
1264  ExcMessage("While sorting the cells that will be passed for "
1265  "creating a Triangulation object, deal.II found that "
1266  "some but not all cells have a negative volume. (If "
1267  "all cells had a negative volume, they would simply "
1268  "all have been inverted.) This usually happens in "
1269  "hand-generated meshes if one accidentally uses an "
1270  "incorrect convention for ordering the vertices in "
1271  "one or more cells; in that case, you may want to "
1272  "double check that you specified the vertex indices "
1273  "in their correct order. If you are reading a mesh "
1274  "that was created by a mesh generator, then this "
1275  "exception indicates that some of the cells created "
1276  "are so badly distorted that their volume becomes "
1277  "negative; this commonly occurs at complex geometric "
1278  "features, and you may see if the problem can be "
1279  "fixed by playing with the parameters that control "
1280  "mesh properties in your mesh generator, such as "
1281  "the number of cells, the mesh density, etc."));
1282 }
1283 
1284 
1285 
1286 /* ------------------------ explicit instantiations ------------------- */
1287 template class GridReordering<1, 1>;
1288 template class GridReordering<1, 2>;
1289 template class GridReordering<1, 3>;
1290 template class GridReordering<2, 2>;
1291 template class GridReordering<2, 3>;
1292 template class GridReordering<3, 3>;
1293 
1294 DEAL_II_NAMESPACE_CLOSE
static void reorder_cells(std::vector< CellData< dim >> &original_cells, const bool use_new_style_ordering=false)
static const unsigned int invalid_unsigned_int
Definition: types.h:173
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
Definition: point.h:110
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:392
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcNotImplemented()
static void invert_all_cells_of_negative_grid(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &original_cells)
static ::ExceptionBase & ExcMeshNotOrientable()
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int vertices[GeometryInfo< structdim >::vertices_per_cell]
Definition: tria.h:141
static ::ExceptionBase & ExcInternalError()