+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2003 - 2016 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii__grid_reordering_internal_h
-#define dealii__grid_reordering_internal_h
-
-
-#include <deal.II/base/config.h>
-#include <deal.II/grid/tria.h>
-
-#include <vector>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-namespace internal
-{
- /**
- * Implement the algorithm described in the documentation of the
- * GridReordering<3> class.
- *
- * @author Michael Anderson, 2003
- */
- namespace GridReordering3d
- {
- /**
- * A structure indicating the direction of an edge. In the implementation
- * file, we define three objects, <tt>unoriented_edge</tt>,
- * <tt>forward_edge</tt>, and <tt>backward_edge</tt>, that denote whether
- * an edge has already been oriented, whether it is in standard
- * orientation, or whether it has reverse direction. The state that each
- * of these objects encode is stored in the <tt>orientation</tt> member
- * variable -- we would really need only three such values, which we pick
- * in the implementation file, and make sure when we compare such objects
- * that only these three special values are actually used.
- *
- * The reason for this way of implementing things is as follows. Usually,
- * such a property would be implemented as an enum. However, in the
- * previous implementation, a signed integer was used with unoriented=0,
- * forward=+1, and backward=-1. A number of operations, such as equality
- * of ordered edges were mapped to checking whether the product of two
- * edge orientations equals +1. Such arithmetic isn't always portable and
- * sometimes flagged when using -ftrapv with gcc. Using this class instead
- * makes sure that there isn't going to be any arithmetic going on on edge
- * orientations, just comparisons for equality or inequality.
- *
- * @author Wolfgang Bangerth, 2005
- */
- struct EdgeOrientation
- {
- /**
- * A value indicating the orientation.
- */
- char orientation;
-
- /**
- * Comparison operator.
- */
- bool operator == (const EdgeOrientation &edge_orientation) const;
-
- /**
- * Comparison operator.
- */
- bool operator != (const EdgeOrientation &edge_orientation) const;
- };
-
- /**
- * During building the connectivity information we don't need all the
- * heavy duty information about edges that we will need later. So we can
- * save memory and time by using a light-weight class for edges. It stores
- * the two vertices, but no direction, so we make the optimization to
- * store the vertex number in sorted order to allow for easier comparison
- * of edge objects.
- */
- struct CheapEdge
- {
- /**
- * The first node
- */
- const unsigned int node0;
-
- /**
- * The second node
- */
- const unsigned int node1;
-
- /**
- * Constructor. Take the vertex numbers and store them sorted.
- */
- CheapEdge (const unsigned int n0,
- const unsigned int n1);
-
- /**
- * Need a partial ordering for sorting algorithms
- */
- bool operator< (const CheapEdge &e2) const;
- };
-
-
-
- /**
- * A connectivity and orientation aware edge class.
- */
- struct Edge
- {
- /**
- * Simple constructor
- */
- Edge (const unsigned int n0,
- const unsigned int n1);
-
- /**
- * The IDs for the end nodes
- */
- unsigned int nodes[2];
-
- /**
- * Whether the edge has not already been oriented, points from node 0 to
- * node 1, or the reverse. The initial state of this flag is unoriented.
- */
- EdgeOrientation orientation_flag;
-
- /**
- * Used to determine which "sheet" or equivalence class of parallel
- * edges the edge falls in when oriented. numbers::invalid_unsigned_int
- * means not yet decided. This is also the default value after
- * construction. Each edge will later be assigned an index greater than
- * zero.
- */
- unsigned int group;
-
- /**
- * Indices of neighboring cubes.
- */
- std::vector<unsigned int> neighboring_cubes;
- };
-
- /**
- * A connectivity and orientation aware cell.
- *
- * The connectivity of the cell is not contained within. (This was for
- * flexibility in using deal.II's ordering of edges or the XDA format
- * etc.) For this information we need the ElemInfo class.
- *
- * One thing we do know is that the first four edges in the edge class are
- * parallel, as are the second four, and the third four.
- *
- * TODO: Need to move connectivity information out of cell and into edge.
- */
- struct Cell
- {
- /**
- * Default Constructor
- */
- Cell ();
-
- /**
- * The IDs for each of the edges.
- */
- unsigned int edges[GeometryInfo<3>::lines_per_cell];
-
- /**
- * The IDs for each of the nodes.
- */
- unsigned int nodes[GeometryInfo<3>::vertices_per_cell];
-
- /**
- * Which way do the edges point. Whether node 0 of the edge is the base
- * of the edge in local element (1) or node 1 is the base (-1).
- */
- EdgeOrientation local_orientation_flags[GeometryInfo<3>::lines_per_cell];
-
- /**
- * An internal flag used to determine whether the cell is in the queue
- * of cells to be oriented in the current sheet.
- */
- bool waiting_to_be_processed;
- };
-
-
- /**
- * This holds all the pieces for orientation together.
- *
- * Contains lists of nodes, edges and cells. As well as the information
- * about how they all connect together.
- */
- class Mesh
- {
- public:
- /**
- * Default Constructor
- */
- Mesh (const std::vector<CellData<3> > &incubes);
-
- /**
- * Export the data of this object to the deal.II format that the
- * Triangulation class wants as input.
- */
- void
- export_to_deal_format (std::vector<CellData<3> > &outcubes) const;
-
- private:
- /**
- * The list of edges
- */
- std::vector<Edge> edge_list;
-
- /**
- * The list of cells
- */
- std::vector<Cell> cell_list;
-
- /**
- * Check whether every cell in the mesh is sensible.
- */
- void sanity_check() const;
-
- /**
- * Given the cell list, build the edge list and all the connectivity
- * information and other stuff that we will need later.
- */
- void build_connectivity ();
-
- /**
- * Unimplemented private copy constructor to disable it.
- */
- Mesh (const Mesh &);
-
- /**
- * Unimplemented private assignment operator to disable it.
- */
- Mesh &operator=(const Mesh &);
-
- /**
- * Checks that each edge going into a node is correctly set up.
- */
- void sanity_check_node (const Cell &cell,
- const unsigned int local_node_num) const;
-
- /**
- * Let the orienter access out private fields.
- */
- friend class Orienter;
- };
-
-
- /**
- * The class that orients the edges of a triangulation in 3d. The member
- * variables basically only store the present state of the algorithm.
- */
- class Orienter
- {
- public:
- /**
- * Orient the given mesh. Creates an object of the present type and lets
- * that toil away at the task.
- *
- * This function is the single entry point to the functionality of this
- * class.
- *
- * Return, whether a consistent orientation of lines was possible for
- * the given mesh.
- */
- static
- bool
- orient_mesh (std::vector<CellData<3> > &incubes);
-
- private:
- /**
- * Internal representation of the given list of cells, including
- * connectivity information and the like.
- */
- Mesh mesh;
-
- /**
- * The cube we're looking at presently.
- */
- unsigned int cur_posn;
-
- /**
- * We have fully oriented all cubes before this one.
- */
- unsigned int marker_cube;
-
- /**
- * The index of the sheet or equivalence class we are presently
- * processing.
- */
- unsigned int cur_edge_group;
-
- /**
- * Indices of the cells to be processed within the present sheet. If a
- * cell is being processed presently, it is taken from this list.
- */
- std::vector<int> sheet_to_process;
-
-
- /**
- * Which edges of the current cell have been oriented during the current
- * iteration. Is reset when moving on to the next cube.
- */
- bool edge_orient_array[12];
-
- /**
- * Constructor. Take a list of cells and set up the internal data
- * structures of the mesh member variable.
- *
- * Since it is private, the only entry point of this class is the static
- * function orient_mesh().
- */
- Orienter (const std::vector<CellData<3> > &incubes);
-
- /**
- * Orient all the edges of a mesh.
- *
- * Return, whether this action was carried out successfully.
- */
- bool orient_edges ();
-
- /**
- * Given oriented edges, rotate the cubes so that the edges are in
- * standard direction.
- */
- void orient_cubes ();
-
- bool get_next_unoriented_cube ();
-
- /**
- * Return whether the cell with cell number @p cell_num is fully
- * oriented.
- */
- bool is_oriented (const unsigned int cell_num) const;
-
- bool orient_edges_in_current_cube ();
- bool orient_edge_set_in_current_cube (const unsigned int edge_set);
- bool orient_next_unoriented_edge ();
-
- /**
- * Return whether the cell is consistently oriented at present (i.e.
- * only considering those edges that are already oriented. This is a
- * sanity check that should be called from inside an assert macro.
- */
- bool cell_is_consistent (const unsigned int cell_num) const;
-
-
- void get_adjacent_cubes ();
- bool get_next_active_cube ();
- };
- } // namespace GridReordering3d
-} // namespace internal
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
#include <deal.II/grid/grid_reordering.h>
-#include <deal.II/grid/grid_reordering_internal.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/std_cxx11/bind.h>
+#include <deal.II/base/timer.h>
#include <algorithm>
#include <set>
{
namespace GridReordering2d
{
+ DeclExceptionMsg (ExcMeshNotOrientable,
+ "The edges of the mesh are not consistently orientable.");
+
/**
* A simple data structure denoting an edge, i.e., the ordered pair
* of its vertex indices. This is only used in the is_consistent()
template <int dim>
struct ParallelEdges
{
+ /**
+ * An array that contains the indices of dim edges that can
+ * serve as (arbitrarily chosen) starting points for the
+ * dim sets of parallel edges within each cell.
+ */
static const unsigned int starter_edges[dim];
- static const unsigned int parallel_edges[GeometryInfo<dim>::lines_per_cell][(1<<(dim-1)) - 1];
+
+ /**
+ * Number and indices of all of those edges parallel to each of the
+ * edges in a cell.
+ */
+ static const unsigned int n_other_parallel_edges = (1<<(dim-1)) - 1;
+ static const unsigned int parallel_edges[GeometryInfo<dim>::lines_per_cell][n_other_parallel_edges];
};
template <>
class AdjacentCells<2>
{
public:
+ /**
+ * An iterator that allows iterating over all cells adjacent
+ * to the edge represented by the current object.
+ */
typedef const AdjacentCell *const_iterator;
/**
};
+
+ /**
+ * A class that represents all of the cells adjacent to a given edge.
+ * This class corresponds to the 3d case where each edge can have an
+ * arbitrary number of adjacent cells. We represent this as a
+ * std::vector<AdjacentCell>, from which class the current one is
+ * derived and from which it inherits all of its member functions.
+ */
+ template <>
+ class AdjacentCells<3> : public std::vector<AdjacentCell>
+ {};
+
+
/**
* A class that describes all of the relevant properties of an
* edge. For the purpose of what we do here, that includes the
* base class.
*/
template <>
- class EdgeDeltaSet<3> : std::set<unsigned int>
+ class EdgeDeltaSet<3> : public std::set<unsigned int>
{};
* (global) parallel to the one identified by the @p cell and
* within it the one with index @p local_edge.
*/
+ template <int dim>
void
- orient_one_set_of_parallel_edges (const std::vector<Cell<2> > &cells,
- std::vector<Edge<2> > &edges,
- const unsigned int cell,
- const unsigned int local_edge)
+ 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)
{
- const unsigned int dim = 2;
-
// choose the direction of the first edge. we have free choice
// here and could simply choose "forward" if that's what pleases
// us. however, for backward compatibility with the previous
// same effect by modifying the rule above to choose the
// direction of the starting edge of this parallel set
// *opposite* to what it looks like in the current cell
+ //
+ // this bug only existed in the 2d implementation since there
+ // were different implementations for 2d and 3d. consequently,
+ // only replicate it for the 2d case and be "intuitive" in 3d.
if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0]
==
cells[cell].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices (local_edge, 0)])
// orient initial edge *opposite* to the way it is in the cell
// (see above for the reason)
- edges[cells[cell].edge_indices[local_edge]].orientation_status = Edge<dim>::backward;
+ edges[cells[cell].edge_indices[local_edge]].orientation_status = (dim == 2 ?
+ Edge<dim>::backward :
+ Edge<dim>::forward);
else
{
Assert (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0]
// orient initial edge *opposite* to the way it is in the cell
// (see above for the reason)
- edges[cells[cell].edge_indices[local_edge]].orientation_status = Edge<dim>::forward;
+ edges[cells[cell].edge_indices[local_edge]].orientation_status = (dim == 2 ?
+ Edge<dim>::forward :
+ Edge<dim>::backward);
}
// walk outward from the given edge as described in
{
Delta_k.clear ();
- for (EdgeDeltaSet<dim>::const_iterator delta = Delta_k_minus_1.begin();
+ for (typename EdgeDeltaSet<dim>::const_iterator delta = Delta_k_minus_1.begin();
delta != Delta_k_minus_1.end(); ++delta)
{
Assert (edges[*delta].orientation_status != Edge<dim>::not_oriented,
ExcInternalError());
// now go through the cells adjacent to this edge
- for (AdjacentCells<dim>::const_iterator
+ for (typename AdjacentCells<dim>::const_iterator
adjacent_cell = edges[*delta].adjacent_cells.begin();
adjacent_cell != edges[*delta].adjacent_cells.end(); ++adjacent_cell)
{
edges[*delta].vertex_indices[0]
:
edges[*delta].vertex_indices[1]);
- const unsigned int first_edge_vertex_in_K = cells[K].vertex_indices[GeometryInfo<dim>::face_to_cell_vertices(delta_is_edge_in_K, 0)];
+ const unsigned int first_edge_vertex_in_K
+ = cells[K].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(delta_is_edge_in_K, 0)];
Assert (first_edge_vertex == first_edge_vertex_in_K
||
- first_edge_vertex == cells[K].vertex_indices[GeometryInfo<dim>::face_to_cell_vertices(delta_is_edge_in_K, 1)],
+ first_edge_vertex == cells[K].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
ExcInternalError());
- // now figure out which direction the opposite edge needs to be into.
- const unsigned int opposite_edge
- = cells[K].edge_indices[ParallelEdges<2>::parallel_edges[delta_is_edge_in_K][0]];
- const unsigned int first_opposite_edge_vertex
- = cells[K].vertex_indices[GeometryInfo<dim>::face_to_cell_vertices(
- ParallelEdges<dim>::parallel_edges[delta_is_edge_in_K][0],
- (first_edge_vertex == first_edge_vertex_in_K
- ?
- 0
- :
- 1))];
-
- const Edge<dim>::OrientationStatus opposite_edge_orientation
- = (edges[opposite_edge].vertex_indices[0]
- ==
- first_opposite_edge_vertex
- ?
- Edge<dim>::forward
- :
- Edge<dim>::backward);
-
- // see if the opposite edge (there is only one in 2d) has already been
- // oriented.
- if (edges[opposite_edge].orientation_status == Edge<dim>::not_oriented)
+ // now figure out which direction the each of the "opposite" edges
+ // needs to be oriented into.
+ for (unsigned int o_e=0; o_e<ParallelEdges<dim>::n_other_parallel_edges; ++o_e)
{
- // the opposite edge is not yet oriented. do orient it and add it to
- // Delta_k
- edges[opposite_edge].orientation_status = opposite_edge_orientation;
- Delta_k.insert (opposite_edge);
- }
- else
- {
- // the opposite edge has already been oriented. assert that it is
- // consistent with the current one
- Assert (edges[opposite_edge].orientation_status == opposite_edge_orientation,
- ExcInternalError());
+ // get the index of the opposite edge and select which its first
+ // vertex needs to be based on how the current edge is oriented
+ // in the current cell
+ const unsigned int opposite_edge
+ = cells[K].edge_indices[ParallelEdges<dim>::parallel_edges[delta_is_edge_in_K][o_e]];
+ const unsigned int first_opposite_edge_vertex
+ = cells[K].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
+ ParallelEdges<dim>::parallel_edges[delta_is_edge_in_K][o_e],
+ (first_edge_vertex == first_edge_vertex_in_K
+ ?
+ 0
+ :
+ 1))];
+
+ // then determine the orientation of the edge based on
+ // whether the vertex we want to be the edge's first
+ // vertex is already the first vertex of the edge, or
+ // whether it points in the opposite direction
+ const typename Edge<dim>::OrientationStatus opposite_edge_orientation
+ = (edges[opposite_edge].vertex_indices[0]
+ ==
+ first_opposite_edge_vertex
+ ?
+ Edge<dim>::forward
+ :
+ Edge<dim>::backward);
+
+ // see if the opposite edge (there is only one in 2d) has already been
+ // oriented.
+ if (edges[opposite_edge].orientation_status == Edge<dim>::not_oriented)
+ {
+ // the opposite edge is not yet oriented. do orient it and add it to
+ // Delta_k
+ edges[opposite_edge].orientation_status = opposite_edge_orientation;
+ Delta_k.insert (opposite_edge);
+ }
+ else
+ {
+ // this opposite edge has already been oriented. it should be
+ // consistent with the current one in 2d, while in 3d it may in fact
+ // be mis-oriented, and in that case the mesh will not be
+ // orientable. indicate this by throwing an exception that we can
+ // catch further up; this has the advantage that we can propagate
+ // through a couple of functions without having to do error
+ // checking and without modifying the 'cells' array that the
+ // user gave us
+ if (dim == 2)
+ {
+ Assert (edges[opposite_edge].orientation_status == opposite_edge_orientation,
+ ExcMeshNotOrientable());
+ }
+ else if (dim == 3)
+ {
+ if (edges[opposite_edge].orientation_status != opposite_edge_orientation)
+ throw ExcMeshNotOrientable ();
+ }
+ else
+ Assert (false, ExcNotImplemented());
+ }
}
}
}
* system matches the ones of the adjacent edges. Store the
* rotated order of vertices in <code>raw_cells[cell_index]</code>.
*/
+ template <int dim>
void
- rotate_cell (const std::vector<Cell<2> > &cell_list,
- const std::vector<Edge<2> > &edge_list,
- const unsigned int cell_index,
- std::vector<CellData<2> > &raw_cells)
+ 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)
{
- // find the first vertex of the cell. this is the
- // vertex where two edges originate, so for
- // each of the four edges record which the
- // starting vertex is
- unsigned int starting_vertex_of_edge[4];
- for (unsigned int e=0; e<4; ++e)
+ // find the first vertex of the cell. this is the vertex where dim edges
+ // originate, so for each of the edges record which the starting vertex is
+ unsigned int starting_vertex_of_edge[GeometryInfo<dim>::lines_per_cell];
+ for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_cell; ++e)
{
Assert (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status
- != Edge<2>::not_oriented,
+ != Edge<dim>::not_oriented,
ExcInternalError());
- if (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status == Edge<2>::forward)
+ if (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status == Edge<dim>::forward)
starting_vertex_of_edge[e] = edge_list[cell_list[cell_index].edge_indices[e]].vertex_indices[0];
else
starting_vertex_of_edge[e] = edge_list[cell_list[cell_index].edge_indices[e]].vertex_indices[1];
}
- // find the vertex number that appears twice. this must either be
- // the first, second, or third vertex in the list. because edges
- // zero and one don't share any vertices, and the same for edges
- // two and three, the possibilities can easily be enumerated
- unsigned int starting_vertex_of_cell = numbers::invalid_unsigned_int;
- if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2])
- ||
- (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
- starting_vertex_of_cell = starting_vertex_of_edge[0];
- else if ((starting_vertex_of_edge[1] == starting_vertex_of_edge[2])
- ||
- (starting_vertex_of_edge[1] == starting_vertex_of_edge[3]))
- starting_vertex_of_cell = starting_vertex_of_edge[1];
- else
- Assert (false, ExcInternalError());
+ // find the vertex number that appears dim times. this will then be
+ // the vertex at which we want to locate the origin of the cell's
+ // coordinate system (i.e., vertex 0)
+ unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int;
+ switch (dim)
+ {
+ case 2:
+ {
+ // in 2d, we can simply enumerate the possibilities where the
+ // origin may be located because edges zero and one don't share
+ // any vertices, and the same for edges two and three
+ if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2])
+ ||
+ (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
+ origin_vertex_of_cell = starting_vertex_of_edge[0];
+ else if ((starting_vertex_of_edge[1] == starting_vertex_of_edge[2])
+ ||
+ (starting_vertex_of_edge[1] == starting_vertex_of_edge[3]))
+ origin_vertex_of_cell = starting_vertex_of_edge[1];
+ else
+ Assert (false, ExcInternalError());
+
+ break;
+ }
+
+ case 3:
+ {
+ // one could probably do something similar in 3d, but that seems
+ // more complicated than one wants to write down. just go
+ // through the list of possible starting vertices and check
+ for (origin_vertex_of_cell=0;
+ origin_vertex_of_cell<GeometryInfo<dim>::vertices_per_cell;
+ ++origin_vertex_of_cell)
+ if (std::count (&starting_vertex_of_edge[0],
+ &starting_vertex_of_edge[0]+GeometryInfo<dim>::lines_per_cell,
+ cell_list[cell_index].vertex_indices[origin_vertex_of_cell])
+ == dim)
+ break;
+ Assert (origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell,
+ ExcInternalError());
+
+ break;
+ }
- // now rotate raw_cells[cell_index] until the starting indices match.
- // take into account the ordering of vertices (not in clockwise
- // or counter-clockwise sense)
- while (raw_cells[cell_index].vertices[0] != starting_vertex_of_cell)
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ // now rotate raw_cells[cell_index] in such a way that its orientation
+ // matches that of cell_list[cell_index]
+ switch (dim)
+ {
+ case 2:
+ {
+ // in 2d, we can literally rotate the cell until its origin
+ // matches the one that we have determined above should be
+ // the origin vertex
+ //
+ // when doing a rotation, take into account the ordering of
+ // vertices (not in clockwise or counter-clockwise sense)
+ while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell)
+ {
+ const unsigned int tmp = raw_cells[cell_index].vertices[0];
+ raw_cells[cell_index].vertices[0] = raw_cells[cell_index].vertices[1];
+ raw_cells[cell_index].vertices[1] = raw_cells[cell_index].vertices[3];
+ raw_cells[cell_index].vertices[3] = raw_cells[cell_index].vertices[2];
+ raw_cells[cell_index].vertices[2] = tmp;
+ }
+ break;
+ }
+
+ case 3:
{
- const unsigned int tmp = raw_cells[cell_index].vertices[0];
- raw_cells[cell_index].vertices[0] = raw_cells[cell_index].vertices[1];
- raw_cells[cell_index].vertices[1] = raw_cells[cell_index].vertices[3];
- raw_cells[cell_index].vertices[3] = raw_cells[cell_index].vertices[2];
- raw_cells[cell_index].vertices[2] = tmp;
+ // in 3d, the situation is a bit more complicated. from above, we
+ // now know which vertex is at the origin (because 3 edges originate
+ // from it), but that still leaves 3 possible rotations of the cube.
+ // the important realization is that we can choose any of them:
+ // in all 3 rotations, all edges originate from the one vertex,
+ // and that fixes the directions of all 12 edges in the cube because
+ // these 3 cover all 3 equivalence classes! consequently, we can
+ // select an arbitrary one among the permutations -- for
+ // example the following ones:
+ static const unsigned int cube_permutations[8][8] =
+ {
+ {0,1,2,3,4,5,6,7},
+ {1,5,3,7,0,4,2,6},
+ {2,6,0,4,3,7,1,5},
+ {3,2,1,0,7,6,5,4},
+ {4,0,6,2,5,1,7,3},
+ {5,4,7,6,1,0,3,2},
+ {6,7,4,5,2,3,0,1},
+ {7,3,5,1,6,2,4,0}
+ };
+
+ unsigned int temp_vertex_indices[GeometryInfo<dim>::vertices_per_cell];
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ temp_vertex_indices[v]
+ = raw_cells[cell_index].vertices[cube_permutations[origin_vertex_of_cell][v]];
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
+
+ break;
+ }
+
+ default:
+ {
+ Assert (false, ExcNotImplemented());
+ }
}
}
// see which edge sets are still not oriented
//
// we do not need to look at each edge because if we orient edge
- // 0, we will end up with edge 1 also oriented. there are only
- // dim independent sets of edges
+ // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
+ // will be 3 other edges that are also oriented). there are only
+ // dim independent sets of edges, so loop over these.
+ //
+ // we need to check whether each one of these starter edges may
+ // already be oriented because the line (sheet) that connects
+ // globally parallel edges may be self-intersecting in the
+ // current cell
for (unsigned int l=0; l<dim; ++l)
if (edge_list[cell_list[next_cell_with_unoriented_edge].edge_indices[ParallelEdges<dim>::starter_edges[l]]].orientation_status
== Edge<dim>::not_oriented)
GridReordering<2,3>::reorder_cells (std::vector<CellData<2> > &cells,
const bool use_new_style_ordering)
{
- // if necessary, convert to old-style format
- if (use_new_style_ordering)
- reorder_new_to_old_style(cells);
+ // if necessary, convert to old (compatibility) to new-style format
+ if (!use_new_style_ordering)
+ reorder_old_to_new_style(cells);
- GridReordering<2>::reorder_cells(cells);
+ // check if grids are already
+ // consistent. if so, do
+ // nothing. if not, then do the
+ // reordering
+ if (!internal::GridReordering2d::is_consistent (cells))
+ internal::GridReordering2d::reorient(cells);
// and convert back if necessary
- if (use_new_style_ordering)
- reorder_old_to_new_style(cells);
+ if (!use_new_style_ordering)
+ reorder_new_to_old_style(cells);
}
-namespace internal
-{
- namespace GridReordering3d
- {
- DeclException1 (ExcGridOrientError,
- char *,
- << "Grid Orientation Error: " << arg1);
-
- const EdgeOrientation unoriented_edge = {'u'};
- const EdgeOrientation forward_edge = {'f'};
- const EdgeOrientation backward_edge = {'b'};
-
-
- inline
- bool
- EdgeOrientation::
- operator == (const EdgeOrientation &edge_orientation) const
- {
- Assert ((orientation == 'u') || (orientation == 'f') || (orientation == 'b'),
- ExcInternalError());
- return orientation == edge_orientation.orientation;
- }
-
-
-
- inline
- bool
- EdgeOrientation::
- operator != (const EdgeOrientation &edge_orientation) const
- {
- return ! (*this == edge_orientation);
- }
-
-
-
- namespace ElementInfo
- {
- /**
- * The numbers of the edges
- * coming into node i are
- * given by
- * edge_to_node[i][k] where
- * k=0,1,2.
- */
- static const unsigned int edge_to_node[8][3] =
- {
- {0,4,8},
- {0,5,9},
- {3,5,10},
- {3,4,11},
- {1,7,8},
- {1,6,9},
- {2,6,10},
- {2,7,11}
- };
-
-
- /**
- * The orientation of edge
- * coming into node i is
- * given by
- * edge_to_node_orient[i][k]
- * where k=0,1,2. 1 means the
- * given node is the start of
- * the edge -1 means the end
- * of the edge.
- */
- static const EdgeOrientation edge_to_node_orient[8][3] =
- {
- {forward_edge, forward_edge, forward_edge},
- {backward_edge, forward_edge, forward_edge},
- {backward_edge, backward_edge, forward_edge},
- {forward_edge, backward_edge, forward_edge},
- {forward_edge, forward_edge, backward_edge},
- {backward_edge, forward_edge, backward_edge},
- {backward_edge, backward_edge, backward_edge},
- {forward_edge, backward_edge, backward_edge}
- };
-
- /**
- * nodesonedge[i][0] is the
- * start node for edge i.
- * nodesonedge[i][1] is the
- * end node for edge i.
- */
- static const unsigned int nodes_on_edge[12][2] =
- {
- {0,1},
- {4,5},
- {7,6},
- {3,2},
- {0,3},
- {1,2},
- {5,6},
- {4,7},
- {0,4},
- {1,5},
- {2,6},
- {3,7}
- };
- }
-
-
- CheapEdge::CheapEdge (const unsigned int n0,
- const unsigned int n1)
- :
- // sort the
- // entries so
- // that
- // node0<node1
- node0(std::min (n0, n1)),
- node1(std::max (n0, n1))
- {}
-
-
-
- bool CheapEdge::operator< (const CheapEdge &e2) const
- {
- if (node0 < e2.node0) return true;
- if (node0 > e2.node0) return false;
- if (node1 < e2.node1) return true;
- return false;
- }
-
-
- Edge::Edge (const unsigned int n0,
- const unsigned int n1)
- :
- orientation_flag (unoriented_edge),
- group (numbers::invalid_unsigned_int)
- {
- nodes[0] = n0;
- nodes[1] = n1;
- }
-
-
-
- Cell::Cell ()
- {
- for (unsigned int i=0; i<GeometryInfo<3>::lines_per_cell; ++i)
- {
- edges[i] = numbers::invalid_unsigned_int;
- local_orientation_flags[i] = forward_edge;
- }
-
- for (unsigned int i=0; i<GeometryInfo<3>::vertices_per_cell; ++i)
- nodes[i] = numbers::invalid_unsigned_int;
-
- waiting_to_be_processed = false;
- }
-
-
-
- Mesh::Mesh (const std::vector<CellData<3> > &incubes)
- {
- // copy the cells into our own
- // internal data format.
- const unsigned int numelems = incubes.size();
- for (unsigned int i=0; i<numelems; ++i)
- {
- Cell the_cell;
- std::copy (&incubes[i].vertices[0],
- &incubes[i].vertices[GeometryInfo<3>::vertices_per_cell],
- &the_cell.nodes[0]);
-
- cell_list.push_back(the_cell);
- }
-
- // then build edges and
- // connectivity
- build_connectivity ();
- }
-
-
-
- void
- Mesh::sanity_check () const
- {
- for (unsigned int i=0; i<cell_list.size(); ++i)
- for (unsigned int j=0; j<8; ++j)
- sanity_check_node (cell_list[i], j);
- }
-
-
-
- void
- Mesh::sanity_check_node (const Cell &c,
- const unsigned int local_node_num) const
- {
-#ifdef DEBUG
- // check that every edge
- // coming into a node has the
- // same node value
-
- // Get the Local Node Numbers
- // of the incoming edges
- const unsigned int e0 = ElementInfo::edge_to_node[local_node_num][0];
- const unsigned int e1 = ElementInfo::edge_to_node[local_node_num][1];
- const unsigned int e2 = ElementInfo::edge_to_node[local_node_num][2];
-
- // Global Edge Numbers
- const unsigned int ge0 = c.edges[e0];
- const unsigned int ge1 = c.edges[e1];
- const unsigned int ge2 = c.edges[e2];
-
- const EdgeOrientation or0 = ElementInfo::edge_to_node_orient[local_node_num][0] ==
- c.local_orientation_flags[e0] ?
- forward_edge : backward_edge;
- const EdgeOrientation or1 = ElementInfo::edge_to_node_orient[local_node_num][1] ==
- c.local_orientation_flags[e1] ?
- forward_edge : backward_edge;
- const EdgeOrientation or2 = ElementInfo::edge_to_node_orient[local_node_num][2] ==
- c.local_orientation_flags[e2] ?
- forward_edge : backward_edge;
-
- // Make sure that edges agree
- // what the current node should
- // be.
- Assert ((edge_list[ge0].nodes[or0 == forward_edge ? 0 : 1] ==
- edge_list[ge1].nodes[or1 == forward_edge ? 0 : 1])
- &&
- (edge_list[ge1].nodes[or1 == forward_edge ? 0 : 1] ==
- edge_list[ge2].nodes[or2 == forward_edge ? 0 : 1]),
- ExcMessage ("This message does not satisfy the internal "
- "consistency check"));
-#else
- (void)c;
- (void)local_node_num;
-#endif
- }
-
-
-
- // This is the guts of the matter...
- void Mesh::build_connectivity ()
- {
- const unsigned int n_cells = cell_list.size();
-
- unsigned int n_edges = 0;
- // Correctly build the edge
- // list
- {
- // edge_map stores the
- // edge_number associated
- // with a given CheapEdge
- std::map<CheapEdge,unsigned int> edge_map;
- unsigned int ctr = 0;
- for (unsigned int cur_cell_id = 0;
- cur_cell_id<n_cells;
- ++cur_cell_id)
- {
- // Get the local node
- // numbers on edge
- // edge_num
- const Cell &cur_cell = cell_list[cur_cell_id];
-
- for (unsigned short int edge_num = 0;
- edge_num<12;
- ++edge_num)
- {
- unsigned int gl_edge_num = 0;
- EdgeOrientation l_edge_orient = forward_edge;
-
- // Construct the
- // CheapEdge
- const unsigned int
- node0 = cur_cell.nodes[ElementInfo::nodes_on_edge[edge_num][0]],
- node1 = cur_cell.nodes[ElementInfo::nodes_on_edge[edge_num][1]];
- const CheapEdge cur_edge (node0, node1);
-
- if (edge_map.count(cur_edge) == 0)
- // Edge not in map
- {
- // put edge in
- // hash map with
- // ctr value;
- edge_map[cur_edge] = ctr;
- gl_edge_num = ctr;
-
- // put the edge
- // into the
- // global edge
- // list
- edge_list.push_back(Edge(node0,node1));
- ctr++;
- }
- else
- {
- // get edge_num
- // from hash_map
- gl_edge_num = edge_map[cur_edge];
- if (edge_list[gl_edge_num].nodes[0] != node0)
- l_edge_orient = backward_edge;
- }
- // set edge number to
- // edgenum
- cell_list[cur_cell_id].edges[edge_num] = gl_edge_num;
- cell_list[cur_cell_id].local_orientation_flags[edge_num]
- = l_edge_orient;
- }
- }
- n_edges = ctr;
- }
-
- // Count each of the edges.
- {
- std::vector<int> edge_count(n_edges,0);
-
-
- // Count every time an edge
- // occurs in a cube.
- for (unsigned int cur_cell_id=0; cur_cell_id<n_cells; ++cur_cell_id)
- for (unsigned short int edge_num = 0; edge_num<12; ++edge_num)
- ++edge_count[cell_list[cur_cell_id].edges[edge_num]];
-
- // So we now know how many
- // cubes contain a given
- // edge. Just need to store
- // the list of cubes in the
- // edge
-
- // Allocate the space for the
- // neighbor list
- for (unsigned int cur_edge_id=0; cur_edge_id<n_edges; ++cur_edge_id)
- edge_list[cur_edge_id].neighboring_cubes
- .resize (edge_count[cur_edge_id]);
-
- // Store the position of the
- // current neighbor in the
- // edge's neighbor list
- std::vector<int> cur_cell_edge_list_posn(n_edges,0);
- for (unsigned int cur_cell_id=0; cur_cell_id<n_cells; ++cur_cell_id)
- for (unsigned short int edge_num=0; edge_num<12; ++edge_num)
- {
- const unsigned int
- gl_edge_id = cell_list[cur_cell_id].edges[edge_num];
- Edge &cur_edge = edge_list[gl_edge_id];
- cur_edge.neighboring_cubes[cur_cell_edge_list_posn[gl_edge_id]]
- = cur_cell_id;
- cur_cell_edge_list_posn[gl_edge_id]++;
- }
- }
- }
-
-
-
- void
- Mesh::export_to_deal_format (std::vector<CellData<3> > &outcubes) const
- {
- Assert (outcubes.size() == cell_list.size(),
- ExcInternalError());
-
- // simply overwrite the output
- // array with the new
- // information
- for (unsigned int i=0; i<cell_list.size(); ++i)
- std::copy (&cell_list[i].nodes[0],
- &cell_list[i].nodes[GeometryInfo<3>::vertices_per_cell],
- &outcubes[i].vertices[0]);
- }
-
-
-
- Orienter::Orienter (const std::vector<CellData<3> > &incubes)
- :
- mesh (incubes),
- cur_posn (0),
- marker_cube (0),
- cur_edge_group (0)
- {
- for (unsigned int i = 0; i<12; ++i)
- edge_orient_array[i] = false;
- }
-
-
-
- bool Orienter::orient_mesh (std::vector<CellData<3> > &incubes)
- {
- Orienter orienter (incubes);
-
- // First check that the mesh is
- // sensible
- orienter.mesh.sanity_check ();
-
- // Orient the mesh
-
- // if not successful, break here, else go
- // on
- if (!orienter.orient_edges ())
- return false;
-
- // Now we have a bunch of oriented
- // edges int the structure we only
- // have to turn the cubes so they
- // match the edge orientation.
- orienter.orient_cubes ();
-
- // Copy the elements from our
- // internal structure back into
- // their original location.
- orienter.mesh.export_to_deal_format (incubes);
- // reordering was successful
- return true;
- }
-
- /**
- * This assigns an orientation
- * to each edge so that every
- * cube is a rotated Deal.II
- * cube.
- */
- bool Orienter::orient_edges ()
- {
- // While there are still cubes
- // to orient
- while (get_next_unoriented_cube())
- // And there are edges in
- // the cube to orient
- while (orient_next_unoriented_edge())
- {
- // Make all the sides
- // in the current set
- // match
- orient_edges_in_current_cube();
-
- // Add the adjacent
- // cubes to the list
- // for processing
- get_adjacent_cubes();
- // Start working on
- // this list of cubes
- while (get_next_active_cube())
- {
- // Make sure the
- // Cube doesn't
- // have a
- // contradiction
- if (!cell_is_consistent(cur_posn))
- return false;
-
- // If we needed to
- // orient any edges
- // in the current
- // cube then we may
- // have to process
- // the neighbor.
- if (orient_edges_in_current_cube())
- get_adjacent_cubes();
- }
-
- // start the next sheet
- // (equivalence class
- // of edges)
- ++cur_edge_group;
- }
- return true;
- }
-
-
-
- bool Orienter::get_next_unoriented_cube ()
- {
- // The last cube in the list
- const unsigned int n_cubes = mesh.cell_list.size();
- // Keep shifting along the list
- // until we find a cube which
- // is not fully oriented or the
- // end.
- while ( (marker_cube<n_cubes) &&
- (is_oriented(marker_cube)) )
- ++marker_cube;
- cur_posn = marker_cube;
- // Return true if we now point
- // at a valid cube.
- return (cur_posn < n_cubes);
- }
-
-
-
- bool Orienter::is_oriented (const unsigned int cell_num) const
- {
- for (unsigned int i=0; i<12; ++i)
- if (mesh.edge_list[mesh.cell_list[cell_num].edges[i]].orientation_flag
- == unoriented_edge)
- return false;
- return true;
- }
-
-
-
- bool
- Orienter::cell_is_consistent(const unsigned int cell_num) const
- {
-
- const Cell &c = mesh.cell_list[cell_num];
-
- // Checks that all oriented
- // edges in the group are
- // oriented consistently.
- for (unsigned int group=0; group<3; ++group)
- {
- // When a nonzero
- // orientation is first
- // encountered in the group
- // it is stored in this
- EdgeOrientation value = unoriented_edge;
- // Loop over all parallel
- // edges
- for (unsigned int i=4*group; i<4*(group+1); ++i)
- {
- // If the edge has
- // orientation
- if ((c.local_orientation_flags[i] !=
- unoriented_edge)
- &&
- (mesh.edge_list[c.edges[i]].orientation_flag !=
- unoriented_edge))
- {
- const EdgeOrientation this_edge_direction
- = (c.local_orientation_flags[i]
- == mesh.edge_list[c.edges[i]].orientation_flag ?
- forward_edge : backward_edge);
-
- // If we haven't
- // seen an oriented
- // edge before,
- // then store its
- // value:
- if (value == unoriented_edge)
- value = this_edge_direction;
- else
- // If we have
- // seen an
- // oriented edge
- // in this group
- // we'd better
- // have the same
- // orientation.
- if (value != this_edge_direction)
- return false;
- }
- }
- }
- return true;
- }
-
-
-
- bool Orienter::orient_next_unoriented_edge ()
- {
- cur_posn = marker_cube;
- const Cell &c = mesh.cell_list[cur_posn];
- unsigned int edge = 0;
-
- // search for the unoriented
- // side
- while ((edge<12) &&
- (mesh.edge_list[c.edges[edge]].orientation_flag !=
- unoriented_edge))
- ++edge;
-
- // if we found none then return
- // false
- if (edge == 12)
- return false;
-
- // Which edge group we're in.
- const unsigned int edge_group = edge/4;
-
- // A sanity check that none of
- // the other edges in the group
- // have been oriented yet Each
- // of the edges in the group
- // should be un-oriented
- for (unsigned int j = edge_group*4; j<edge_group*4+4; ++j)
- Assert (mesh.edge_list[c.edges[j]].orientation_flag ==
- unoriented_edge,
- ExcGridOrientError("Tried to orient edge when other edges "
- "in group are already oriented!"));
-
- // Make the edge alignment
- // match that of the local
- // cube.
- mesh.edge_list[c.edges[edge]].orientation_flag
- = c.local_orientation_flags[edge];
- mesh.edge_list[c.edges[edge]].group = cur_edge_group;
-
- // Remember that we have oriented
- // this edge in the current cell.
- edge_orient_array[edge] = true;
-
- return true;
- }
-
-
-
- bool Orienter::orient_edges_in_current_cube ()
- {
- for (unsigned int edge_group=0; edge_group<3; ++edge_group)
- if (orient_edge_set_in_current_cube(edge_group) == true)
- return true;
-
- return false;
- }
-
-
-
- bool
- Orienter::orient_edge_set_in_current_cube (const unsigned int n)
- {
- const Cell &c = mesh.cell_list[cur_posn];
-
- // Check if any edge is
- // oriented
- unsigned int n_oriented = 0;
- EdgeOrientation glorient = unoriented_edge;
- unsigned int edge_flags = 0;
- unsigned int cur_flag = 1;
- for (unsigned int i = 4*n; i<4*(n+1); ++i, cur_flag<<=1)
- {
- if ((mesh.edge_list[c.edges[i]].orientation_flag !=
- unoriented_edge)
- &&
- (c.local_orientation_flags[i] !=
- unoriented_edge))
- {
- ++n_oriented;
-
- const EdgeOrientation orient
- = (mesh.edge_list[c.edges[i]].orientation_flag ==
- c.local_orientation_flags[i] ?
- forward_edge : backward_edge);
-
- if (glorient == unoriented_edge)
- glorient = orient;
- else
- AssertThrow(orient == glorient,
- ExcGridOrientError("Attempted to Orient Misaligned cube"));
- }
- else
- edge_flags |= cur_flag;
- }
-
- // were any of the sides
- // oriented? were they all
- // already oriented?
- if ((glorient == unoriented_edge) || (n_oriented == 4))
- return false;
-
- // If so orient all edges
- // consistently.
- cur_flag = 1;
- for (unsigned int i=4*n; i<4*(n+1); ++i, cur_flag<<=1)
- if ((edge_flags & cur_flag) != 0)
- {
- mesh.edge_list[c.edges[i]].orientation_flag
- = (c.local_orientation_flags[i] == glorient ?
- forward_edge : backward_edge);
-
- mesh.edge_list[c.edges[i]].group = cur_edge_group;
- // Remember that we have oriented
- // this edge in the current cell.
- edge_orient_array[i] = true;
- }
-
- return true;
- }
-
-
-
- void Orienter::get_adjacent_cubes ()
- {
- const Cell &c = mesh.cell_list[cur_posn];
- for (unsigned int e=0; e<12; ++e)
- // Only need to add the adjacent
- // cubes for edges we recently
- // oriented
- if (edge_orient_array[e] == true)
- {
- const Edge &the_edge = mesh.edge_list[c.edges[e]];
- for (unsigned int local_cube_num = 0;
- local_cube_num < the_edge.neighboring_cubes.size();
- ++local_cube_num)
- {
- const unsigned int
- global_cell_num = the_edge.neighboring_cubes[local_cube_num];
- Cell &ncell = mesh.cell_list[global_cell_num];
-
- // If the cell is waiting to be
- // processed we dont want to add
- // it to the list a second time.
- if (!ncell.waiting_to_be_processed)
- {
- sheet_to_process.push_back(global_cell_num);
- ncell.waiting_to_be_processed = true;
- }
- }
- }
- // we're done with this cube so
- // clear its processing flags.
- for (unsigned int e=0; e<12; ++e)
- edge_orient_array[e] = false;
-
- }
-
-
-
- bool Orienter::get_next_active_cube ()
- {
- // Mark the curent Cube as
- // finished with.
- Cell &c = mesh.cell_list[cur_posn];
- c.waiting_to_be_processed = false;
- if (sheet_to_process.empty() == false)
- {
- cur_posn = sheet_to_process.back();
- sheet_to_process.pop_back();
- return true;
- }
- return false;
- }
-
-
- void Orienter::orient_cubes ()
- {
- // We assume that the mesh has
- // all edges oriented already.
-
- // This is a list of
- // permutations that take node
- // 0 to node i but only rotate
- // the cube. (This set is far
- // from unique (there are 3 for
- // each node - for our
- // algorithm it doesn't matter
- // which of the three we use)
- static const unsigned int CubePermutations[8][8] =
- {
- {0,1,2,3,4,5,6,7},
- {1,2,3,0,5,6,7,4},
- {2,3,0,1,6,7,4,5},
- {3,0,1,2,7,4,5,6},
- {4,7,6,5,0,3,2,1},
- {5,4,7,6,1,0,3,2},
- {6,5,4,7,2,1,0,3},
- {7,6,5,4,3,2,1,0}
- };
-
- // So now we need to work out
- // which node needs to be
- // mapped to the zero node.
- // The trick is that the node
- // that should be the local
- // zero node has three edges
- // coming into it.
- for (unsigned int i=0; i<mesh.cell_list.size(); ++i)
- {
- Cell &the_cell = mesh.cell_list[i];
-
- // This stores whether the
- // global oriented edge
- // points in the same
- // direction as it's local
- // edge on the current
- // cube. (for each edge on
- // the curent cube)
- EdgeOrientation local_edge_orientation[12];
- for (unsigned int j = 0; j<12; ++j)
- {
- // get the global edge
- const Edge &the_edge = mesh.edge_list[the_cell.edges[j]];
- // All edges should be
- // oriented at this
- // stage..
- Assert (the_edge.orientation_flag != unoriented_edge,
- ExcGridOrientError ("Unoriented edge encountered"));
- // calculate whether it
- // points the right way
- // or not
- local_edge_orientation[j] = (the_cell.local_orientation_flags[j] ==
- the_edge.orientation_flag ?
- forward_edge : backward_edge);
- }
-
- // Here the number of
- // incoming edges is
- // tallied for each node.
- unsigned int perm_num = numbers::invalid_unsigned_int;
- for (unsigned int node_num=0; node_num<8; ++node_num)
- {
- // The local edge
- // numbers coming into
- // the node
- const unsigned int e0 = ElementInfo::edge_to_node[node_num][0];
- const unsigned int e1 = ElementInfo::edge_to_node[node_num][1];
- const unsigned int e2 = ElementInfo::edge_to_node[node_num][2];
-
- // The local
- // orientation of the
- // edge coming into the
- // node.
- const EdgeOrientation sign0 = ElementInfo::edge_to_node_orient[node_num][0];
- const EdgeOrientation sign1 = ElementInfo::edge_to_node_orient[node_num][1];
- const EdgeOrientation sign2 = ElementInfo::edge_to_node_orient[node_num][2];
-
- // Add one to the total
- // for each edge
- // pointing in
- Assert (local_edge_orientation[e0] != unoriented_edge,
- ExcInternalError());
- Assert (local_edge_orientation[e1] != unoriented_edge,
- ExcInternalError());
- Assert (local_edge_orientation[e2] != unoriented_edge,
- ExcInternalError());
-
- const unsigned int
- total = (((local_edge_orientation[e0] == sign0) ? 1 : 0)
- +((local_edge_orientation[e1] == sign1) ? 1 : 0)
- +((local_edge_orientation[e2] == sign2) ? 1 : 0));
-
- if (total == 3)
- {
- Assert (perm_num == numbers::invalid_unsigned_int,
- ExcGridOrientError("More than one node with 3 incoming "
- "edges found in curent hex."));
- perm_num = node_num;
- }
- }
- // We should now have a
- // valid permutation number
- Assert (perm_num != numbers::invalid_unsigned_int,
- ExcGridOrientError("No node having 3 incoming edges found in curent hex."));
-
- // So use the appropriate
- // rotation to get the new
- // cube
- unsigned int temp[8];
- for (unsigned int v=0; v<8; ++v)
- temp[v] = the_cell.nodes[CubePermutations[perm_num][v]];
- for (unsigned int v=0; v<8; ++v)
- the_cell.nodes[v] = temp[v];
- }
- }
- } // namespace GridReordering3d
-} // namespace internal
-
-
-
template<>
void
GridReordering<3>::reorder_cells (std::vector<CellData<3> > &cells,
Assert (cells.size() != 0,
ExcMessage("List of elements to orient must have at least one cell"));
- // if necessary, convert to old-style format
- if (use_new_style_ordering)
- reorder_new_to_old_style(cells);
-
- // create a backup to use if GridReordering
- // was not successful
- std::vector<CellData<3> > backup=cells;
-
- // This does the real work
- const bool success=
- internal::GridReordering3d::Orienter::orient_mesh (cells);
+ // if necessary, convert to new-style format
+ if (use_new_style_ordering == false)
+ reorder_old_to_new_style(cells);
- // if reordering was not successful use
- // original connectivity, otherwise do
- // nothing (i.e. use the reordered
- // connectivity)
- if (!success)
- cells=backup;
+ // check if grids are already consistent. if so, do
+ // nothing. if not, then do the reordering
+ if (!internal::GridReordering2d::is_consistent (cells))
+ try
+ {
+ internal::GridReordering2d::reorient(cells);
+ }
+ catch (const internal::GridReordering2d::ExcMeshNotOrientable &)
+ {
+ // the mesh is not orientable
+ }
// and convert back if necessary
- if (use_new_style_ordering)
- reorder_old_to_new_style(cells);
+ if (use_new_style_ordering == false)
+ reorder_new_to_old_style(cells);
}
projects / dealii.git / commitdiff