From ccfff7c44df8e73404bfcade8ec3886db6bedd15 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 27 Oct 2003 14:40:00 +0000 Subject: [PATCH] Final adjustments for conversion to deal.II-style coding. git-svn-id: https://svn.dealii.org/trunk@8166 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/grid/grid_reordering_internal.h | 225 +++++-- .../deal.II/source/grid/grid_reordering.cc | 603 +++++++++--------- 2 files changed, 458 insertions(+), 370 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_reordering_internal.h b/deal.II/deal.II/include/grid/grid_reordering_internal.h index fc0931020c..fffff97f5e 100644 --- a/deal.II/deal.II/include/grid/grid_reordering_internal.h +++ b/deal.II/deal.II/include/grid/grid_reordering_internal.h @@ -398,32 +398,41 @@ namespace internal /** * During building the - * conectivity information we + * 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. + * 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 */ - unsigned int node0; + const unsigned int node0; /** * The second node */ - unsigned int node1; + const unsigned int node1; /** - * A simple constructor + * Constructor. Take the + * vertex numbers and store + * them sorted. */ CheapEdge (const unsigned int n0, const unsigned int n1); - /** + /** * Need a partial ordering * for the STL */ @@ -433,7 +442,7 @@ namespace internal /** - * A conectivity and orientation + * A connectivity and orientation * aware edge class. */ struct Edge @@ -442,8 +451,7 @@ namespace internal * Simple constructor */ Edge (const unsigned int n0, - const unsigned int n1, - const int orient = 0); + const unsigned int n1); /** * The IDs for the end nodes @@ -451,21 +459,30 @@ namespace internal unsigned int nodes[2]; /** - * Whether the edge has been - * oriented (0), points from - * node 0 to node 1 (1), or - * the reverse (-1). + * Whether the edge has not + * already been oriented (0), + * points from node 0 to node + * 1 (1), or the reverse + * (-1). The initial state of + * this flag is zero. */ - int orientation_flag; + signed short int orientation_flag; /** * Used to determine which - * "sheet" of parallel edges + * "sheet" or equivalence + * class of parallel edges * the edge falls in when - * oriented. 0 means not yet - * decided. + * oriented. + * deal_II_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. */ - int group; + unsigned int group; /** * Indices of neighboring cubes. @@ -491,7 +508,7 @@ namespace internal * are the second four, and the * third four. * - * TODO: Need to move conectivity information out + * TODO: Need to move connectivity information out * of cell and into edge. */ struct Cell @@ -501,11 +518,6 @@ namespace internal */ Cell (); - /** - * Copy Constructor - */ - Cell (const Cell &c); - /** * The IDs for each of the edges. */ @@ -524,7 +536,7 @@ namespace internal * (1) or node 1 is the base * (-1). */ - int local_orientation_flags[GeometryInfo<3>::lines_per_cell]; + signed int local_orientation_flags[GeometryInfo<3>::lines_per_cell]; /** * An internal flag used to @@ -552,8 +564,19 @@ namespace internal /** * Default Constructor */ - Mesh (); + Mesh (const std::vector > &incubes); + /** + * Export the data of this + * object to the deal.II + * format that the + * @ref{Triangulation} class + * wants as input. + */ + void + export_to_deal_format (std::vector > &outcubes) const; + + private: /** * The list of edges */ @@ -570,7 +593,15 @@ namespace internal */ void sanity_check() const; - private: + /** + * 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. @@ -579,7 +610,7 @@ namespace internal /** * Unimplemented private - * assignemnet operator to + * assignment operator to * disable it. */ Mesh& operator=(const Mesh&); @@ -590,21 +621,55 @@ namespace internal * correctly set up. */ void sanity_check_node (const Cell &cell, - const unsigned int local_node_num) const; + 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: /** - * Constructor. + * 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. */ - Orienter(); + static + void + orient_mesh (std::vector > &incubes); + private: + /** + * Internal representation of + * the given list of cells, + * including connectivity + * information and the like. + */ + Mesh mesh; + /** - * The cube we're looking at now. + * The cube we're looking at + * presently. */ unsigned int cur_posn; @@ -614,41 +679,85 @@ namespace internal */ unsigned int marker_cube; - std::vector sheet_to_process; + /** + * The index of the sheet or + * equivalence class we are + * presently processing. + */ + unsigned int cur_edge_group; - int cur_edge_group; + /** + * Indices of the cells to be + * processed withing the + * present sheet. If a cell + * is being processed + * presently, it is taken + * from this list. + */ + std::vector sheet_to_process; bool edge_orient_array[12]; - - bool orient_edges (Mesh &m); - void orient_cubes (Mesh &m); - - bool get_next_unoriented_cube (Mesh &m); - bool is_oriented (const Mesh &m, - int cell_num); - bool orient_edges_in_current_cube (Mesh &m); - bool orient_edge_set_in_current_cube (Mesh &m, - const unsigned int edge_set); - bool orient_next_unoriented_edge (Mesh &m); - bool consistent (Mesh &m, - int cell_num); + /** + * 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 + * @ref{orient_mesh}. + */ + Orienter (const std::vector > &incubes); + /** + * Orient all the edges of a + * mesh. + */ + void orient_edges (); - void get_adjacent_cubes (Mesh &m); - bool get_next_active_cube (Mesh &m); - }; + /** + * 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; - /** - * Creates the connectivity - * information for the mesh m. - */ - void build_mesh (Mesh &m); + 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 + * consistenty 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 diff --git a/deal.II/deal.II/source/grid/grid_reordering.cc b/deal.II/deal.II/source/grid/grid_reordering.cc index 50629e6a87..7f13ca5622 100644 --- a/deal.II/deal.II/source/grid/grid_reordering.cc +++ b/deal.II/deal.II/source/grid/grid_reordering.cc @@ -59,7 +59,7 @@ namespace internal * used in the is_consistent * function. */ - struct Edge + struct Edge { Edge (const unsigned int v0, const unsigned int v1) @@ -85,10 +85,14 @@ namespace internal { // construct the four edges // in reverse order - const Edge reverse_edges[4] = { Edge (c->vertices[1], c->vertices[0]), - Edge (c->vertices[2], c->vertices[1]), - Edge (c->vertices[2], c->vertices[3]), - Edge (c->vertices[3], c->vertices[0]) }; + const Edge reverse_edges[4] = { Edge (c->vertices[1], + c->vertices[0]), + Edge (c->vertices[2], + c->vertices[1]), + Edge (c->vertices[2], + c->vertices[3]), + Edge (c->vertices[3], + c->vertices[0]) }; // for each of them, check // whether they are already // in the set @@ -583,15 +587,6 @@ namespace internal << "Grid Orientation Error: " << arg1); - // sort two integers - static inline void sort2 (unsigned int &v1, - unsigned int &v2) - { - if (v1>v2) - std::swap (v1, v2); - } - - namespace ElementInfo { /** @@ -678,12 +673,13 @@ namespace internal CheapEdge::CheapEdge (const unsigned int n0, const unsigned int n1) : - node0(n0), node1(n1) - { - // sort the entries so that - // node0::lines_per_cell; ++i) + { + edges[i] = static_cast(-1); + local_orientation_flags[i] = 1; + } + + for (unsigned int i=0; i::vertices_per_cell; ++i) + nodes[i] = static_cast(-1); + + waiting_to_be_processed = false; + } + + + + Mesh::Mesh (const std::vector > &incubes) + { + // copy the cells into our own + // internal data format. + const unsigned int numelems = incubes.size(); + for (unsigned int i=0; i::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 = m.cell_list; - std::vector & edge_list = m.edge_list; + // check that every edge + // coming into a node has the + // same node value - const unsigned int cell_list_length = cell_list.size(); + // 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 int or0 = ElementInfo::edge_to_node_orient[local_node_num][0] * + c.local_orientation_flags[e0]; + const int or1 = ElementInfo::edge_to_node_orient[local_node_num][1] * + c.local_orientation_flags[e1]; + const int or2 = ElementInfo::edge_to_node_orient[local_node_num][2] * + c.local_orientation_flags[e2]; + // Make sure that edges agree + // what the current node should + // be. + Assert ((edge_list[ge0].nodes[or0 == 1 ? 0 : 1] == + edge_list[ge1].nodes[or1 == 1 ? 0 : 1]) + && + (edge_list[ge1].nodes[or1 == 1 ? 0 : 1] == + edge_list[ge2].nodes[or2 == 1 ? 0 : 1]), + ExcMessage ("This message does not satisfy the internal " + "consistency check")); + } + + + + // 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 @@ -715,25 +807,26 @@ namespace internal std::map edge_map; unsigned int ctr = 0; for (unsigned int cur_cell_id = 0; - cur_cell_id cur_cell_edge_list_posn(n_edges,0); - for (unsigned int cur_cell_id = 0; - cur_cell_id > &outcubes) const { - for (unsigned int i=0; i::lines_per_cell; ++i) - { - edges[i] = c.edges[i]; - local_orientation_flags[i] = c.local_orientation_flags[i]; - } - - for (unsigned int i=0; i::vertices_per_cell; ++i) - nodes[i] = c.nodes[i]; - - waiting_to_be_processed = c.waiting_to_be_processed; - } - + Assert (outcubes.size() == cell_list.size(), + ExcInternalError()); - - Cell::Cell () - { - for (unsigned int i=0; i::lines_per_cell; ++i) - { - edges[i] = static_cast(-1); - local_orientation_flags[i] = 1; - } - - for (unsigned int i=0; i::vertices_per_cell; ++i) - nodes[i] = static_cast(-1); - - waiting_to_be_processed = false; + // simply overwrite the output + // array with the new + // information + for (unsigned int i=0; i::vertices_per_cell], + &outcubes[i].vertices[0]); } - - - - Mesh::Mesh () - {} - void - Mesh::sanity_check() const - { - for (unsigned int i=0; i > &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; } - - void - Mesh::sanity_check_node(const Cell &c, - const unsigned int local_node_num) const + + void Orienter::orient_mesh (std::vector > &incubes) { - // 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 int or0 = ElementInfo::edge_to_node_orient[local_node_num][0] * - c.local_orientation_flags[e0]; - const int or1 = ElementInfo::edge_to_node_orient[local_node_num][1] * - c.local_orientation_flags[e1]; - const int or2 = ElementInfo::edge_to_node_orient[local_node_num][2] * - c.local_orientation_flags[e2]; - - // Make sure that edges agree - // what the current node should - // be. - Assert ((edge_list[ge0].nodes[or0 == 1 ? 0 : 1] == - edge_list[ge1].nodes[or1 == 1 ? 0 : 1]) - && - (edge_list[ge1].nodes[or1 == 1 ? 0 : 1] == - edge_list[ge2].nodes[or2 == 1 ? 0 : 1]), - ExcMessage ("This message does not satisfy the internal " - "consistency check")); - } - + Orienter orienter (incubes); + + // First check that the mesh is + // sensible + orienter.mesh.sanity_check (); + // Orient the mesh + orienter.orient_edges (); - Orienter::Orienter () - { - for (unsigned int i = 0; i<12; ++i) - edge_orient_array[i] = 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); } /** - * This assignes an orientation + * This assigns an orientation * to each edge so that every * cube is a rotated Deal.II * cube. */ - bool Orienter::orient_edges (Mesh &m) + void Orienter::orient_edges () { - - // First check that the mesh is - // sensible - m.sanity_check (); - - // We start by looking only at - // the first cube. - cur_posn = 0; - marker_cube = 0; - - // We mark each edge with a - // group number (mostly for - // mesh debugging purposes) - cur_edge_group = 1; // While there are still cubes // to orient - while (get_next_unoriented_cube(m)) - { - // And there are edges in - // the cube to orient - while (orient_next_unoriented_edge(m)) - { - // Make all the sides - // in the current set - // match - orient_edges_in_current_cube(m); - // Add the adjacent - // cubes to the list - // for processing - get_adjacent_cubes(m); - // Start working on - // this list of cubes - while (get_next_active_cube(m)) - { - // Make sure the - // Cube doesn't - // have a - // contradiction - AssertThrow(consistent(m,cur_posn), - ExcGridOrientError("Mesh is Unorientable")); - // 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(m)) - get_adjacent_cubes(m); - } - cur_edge_group++; - } - } - return true; + 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 + AssertThrow(cell_is_consistent(cur_posn), + ExcGridOrientError("Mesh is Unorientable")); + + // 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; + } } - bool Orienter::get_next_unoriented_cube(Mesh &m) + bool Orienter::get_next_unoriented_cube () { // The last cube in the list - unsigned int end_cube_num = m.cell_list.size(); + 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(-1); for (unsigned int node_num=0; node_num<8; ++node_num) { // The local edge @@ -1320,7 +1338,7 @@ namespace internal if (total == 3) { - Assert (perm_num == -1, + Assert (perm_num == static_cast(-1), ExcGridOrientError("More than one node with 3 incoming " "edges found in curent hex.")); perm_num = node_num; @@ -1328,7 +1346,7 @@ namespace internal } // We should now have a // valid permutation number - Assert (perm_num != -1, + Assert (perm_num != static_cast(-1), ExcGridOrientError("No node having 3 incoming edges found in curent hex.")); // So use the apropriate @@ -1346,54 +1364,15 @@ namespace internal -void GridReordering<3>::reorder_cells (std::vector > &incubes) +void +GridReordering<3>::reorder_cells (std::vector > &incubes) { Assert (incubes.size() != 0, ExcMessage("List of elements to orient was of zero length")); // This does the real work - internal::GridReordering3d::Orienter orienter; - - // This is the internal store for - // all global connectivity - // information it starts prety much - // empty. - internal::GridReordering3d::Mesh the_mesh; - - // Copy the cells into our own - // internal data format. - const unsigned int numelems = incubes.size(); - for (unsigned int i=0; i