Reference documentation for deal.II version Git a285ab78cc 2020-08-04 19:09:24 +0200
\(\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\}}\)
grid_reordering.cc
Go to the documentation of this file.
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 
22 
23 #include <algorithm>
24 #include <fstream>
25 #include <functional>
26 #include <iostream>
27 #include <set>
28 
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(
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(
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])
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 (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 
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] ==
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] ==
653  local_edge, 1)],
654  ExcInternalError());
655  Assert(
656  edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
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 (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 
948  template <int dim>
949  void
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)) !=
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
1005  void reorient(std::vector<CellData<1>> &)
1006  {}
1007 } // namespace
1008 
1009 
1010 // anonymous namespace for internal helper functions
1011 namespace
1012 {
1022  void reorder_new_to_old_style(std::vector<CellData<1>> &)
1023  {}
1024 
1025 
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 
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 
1049  void reorder_old_to_new_style(std::vector<CellData<1>> &)
1050  {}
1051 
1052 
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 
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
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
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
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
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
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
1205  const std::vector<Point<3>> &,
1206  std::vector<CellData<2>> &)
1207 {
1208  Assert(false, ExcNotImplemented());
1209 }
1210 
1211 
1212 
1213 template <>
1214 void
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 
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1075
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:196
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
std::vector< unsigned int > vertices
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
Definition: point.h:110
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcMessage(std::string arg1)
double cell_measure< 3 >(const std::vector< Point< 3 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
#define Assert(cond, exc)
Definition: exceptions.h:1411
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
VectorType::value_type * end(VectorType &V)
std::string to_string(const T &t)
Definition: patterns.h:2341
Point< 3 > vertices[4]
__global__ void set(Number *val, const Number s, const size_type N)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
Definition: memory_space.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
VectorType::value_type * begin(VectorType &V)
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()
double cell_measure< 2 >(const std::vector< Point< 2 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()