From 45b69febf0cb1972d0d84bfe6f7ddef6dbdce701 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 23 Nov 2009 04:06:17 +0000 Subject: [PATCH] Cache Triangulation::n_levels(). git-svn-id: https://svn.dealii.org/trunk@20149 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/tria.h | 37 + .../source/grid/tria.all_dimensions.cc | 17 +- deal.II/deal.II/source/grid/tria.cc | 1598 +++++++++-------- deal.II/doc/news/changes.h | 17 +- 4 files changed, 867 insertions(+), 802 deletions(-) diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 5ebfb137db..966560291a 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -212,6 +212,13 @@ namespace internal template <> struct NumberCache<1> { + /** + * The number of levels on + * which we have used + * objects. + */ + unsigned int n_levels; + /** * Number of used lines in the whole * triangulation. @@ -3344,6 +3351,36 @@ Triangulation::vertex_used(const unsigned int index) const } + +template +inline +unsigned int Triangulation::n_levels () const +{ + return number_cache.n_levels; +} + + + +template +inline +unsigned int +Triangulation::n_vertices () const +{ + return vertices.size(); +} + + + +template +inline +const std::vector > & +Triangulation::get_vertices () const +{ + return vertices; +} + + + /* -------------- declaration of explicit specializations ------------- */ template <> unsigned int Triangulation<1,1>::n_raw_lines (const unsigned int level) const; diff --git a/deal.II/deal.II/source/grid/tria.all_dimensions.cc b/deal.II/deal.II/source/grid/tria.all_dimensions.cc index 00cb017d0b..8a7cee1efe 100644 --- a/deal.II/deal.II/source/grid/tria.all_dimensions.cc +++ b/deal.II/deal.II/source/grid/tria.all_dimensions.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -37,7 +37,7 @@ DEAL_II_NAMESPACE_OPEN bool SubCellData::check_consistency (const unsigned int dim) const { - switch (dim) + switch (dim) { case 1: return ((boundary_lines.size() == 0) && @@ -54,9 +54,11 @@ namespace internal namespace Triangulation { - NumberCache<1>::NumberCache () : + NumberCache<1>::NumberCache () + : + n_levels (0), n_lines (0), - n_active_lines (0) + n_active_lines (0) // all other fields are // default constructed {} @@ -66,7 +68,8 @@ namespace internal unsigned int NumberCache<1>::memory_consumption () const { - return (MemoryConsumption::memory_consumption (n_lines) + + return (MemoryConsumption::memory_consumption (n_levels) + + MemoryConsumption::memory_consumption (n_lines) + MemoryConsumption::memory_consumption (n_lines_level) + MemoryConsumption::memory_consumption (n_active_lines) + MemoryConsumption::memory_consumption (n_active_lines_level)); @@ -75,7 +78,7 @@ namespace internal NumberCache<2>::NumberCache () : n_quads (0), - n_active_quads (0) + n_active_quads (0) // all other fields are // default constructed {} @@ -96,7 +99,7 @@ namespace internal NumberCache<3>::NumberCache () : n_hexes (0), - n_active_hexes (0) + n_active_hexes (0) // all other fields are // default constructed {} diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 1d86aaf1a1..5380ba9dc0 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -56,9 +56,9 @@ namespace return (n_active_children == 0) || (n_active_children == cell->n_children()); } - - - + + + // return, wheter a given @p cell will be // coarsened, which is the case if all // children are active and have their coarsen @@ -69,12 +69,12 @@ namespace { // only cells with children should be // considered for coarsening - + if (cell->has_children()) { unsigned int children_to_coarsen=0; const unsigned int n_children=cell->n_children(); - + for (unsigned int c=0; cchild(c)->active() && cell->child(c)->coarsen_flag_set()) @@ -95,7 +95,7 @@ namespace // make sure we do not have to do this at // all... Assert(cell->has_children(), ExcInternalError()); - // ... and then simply return false + // ... and then simply return false return false; } @@ -130,14 +130,14 @@ namespace // expected_face_ref_case, which is no // refinement at all expected_face_ref_case=RefinementCase::no_refinement; - + const typename Triangulation::cell_iterator neighbor=cell->neighbor(face_no); // If we are at the boundary, there is no // neighbor which could refine the face if (neighbor.state()!=IteratorState::valid) return false; - + if (neighbor->has_children()) { // if the neighbor is refined, he may be @@ -156,7 +156,7 @@ namespace return true; } } - + // now, the neighbor is not refined, but // perhaps he will be const RefinementCase nb_ref_flag=neighbor->refine_flag_set(); @@ -194,7 +194,7 @@ namespace } else { - + // the neighbor is coarser. // this is the most complicated // case. It might be, that the @@ -239,7 +239,7 @@ namespace return true; } } - + // so, obviously we were not // one of the children, but a // grandchild. This is only @@ -295,7 +295,7 @@ namespace triangulation.n_levels()); std::vector max_adjacent_cell_level (triangulation.n_vertices(), 0); - + for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) @@ -380,24 +380,59 @@ namespace * to deal with the fact that lines have * levels, whereas for higher dimensions * they do not. + * + * The second argument indicates + * for how many levels the + * Triangulation has objects, + * though the highest levels need + * not contain active cells if they + * have previously all been + * coarsened away. */ template void compute_number_cache (const Triangulation &triangulation, - internal::Triangulation::NumberCache<1> &number_cache) + const unsigned int level_objects, + internal::Triangulation::NumberCache<1> &number_cache) { typedef typename Triangulation::line_iterator line_iterator; typedef typename Triangulation::active_line_iterator active_line_iterator; - + + number_cache.n_levels = 0; + if (level_objects > 0) + { + // check whether there are + // cells on the highest + // levels (there need not be, + // since they might all have + // been coarsened away) + typename Triangulation::raw_cell_iterator + cell = triangulation.last_raw (level_objects-1), + endc = triangulation.end(); + for (; cell!=endc; --cell) + if (cell->used()) + { + // return level of most + // refined existing cell + // (+1 because of + // counting conventions) + number_cache.n_levels = cell->level()+1; + break; + } + + // no cells at all? + Assert (number_cache.n_levels > 0, ExcInternalError()); + } + /////////////////////////////////// // update the number of lines // on the different levels in // the cache - number_cache.n_lines_level.resize (triangulation.n_levels()); + number_cache.n_lines_level.resize (number_cache.n_levels); number_cache.n_lines = 0; - number_cache.n_active_lines_level.resize (triangulation.n_levels()); + number_cache.n_active_lines_level.resize (number_cache.n_levels); number_cache.n_active_lines = 0; // for 1d, lines have levels so take @@ -405,25 +440,25 @@ namespace // globally if (dim == 1) { - for (unsigned int level=0; level void compute_number_cache (const Triangulation &triangulation, - internal::Triangulation::NumberCache<2> &number_cache) + const unsigned int level_objects, + internal::Triangulation::NumberCache<2> &number_cache) { + // update lines and n_levels + compute_number_cache (triangulation, + level_objects, + static_cast&> + (number_cache)); + typedef typename Triangulation::quad_iterator quad_iterator; typedef typename Triangulation::active_quad_iterator active_quad_iterator; - + /////////////////////////////////// // update the number of quads // on the different levels in // the cache - number_cache.n_quads_level.resize (triangulation.n_levels()); + number_cache.n_quads_level.resize (number_cache.n_levels); number_cache.n_quads = 0; - number_cache.n_active_quads_level.resize (triangulation.n_levels()); + number_cache.n_active_quads_level.resize (number_cache.n_levels); number_cache.n_active_quads = 0; // for 2d, quads have levels so take @@ -492,25 +542,25 @@ namespace // globally if (dim == 2) { - for (unsigned int level=0; level&>(number_cache)); } /** @@ -555,6 +601,14 @@ namespace * levels, whereas for higher dimensions * they do not. * + * The second argument indicates + * for how many levels the + * Triangulation has objects, + * though the highest levels need + * not contain active cells if they + * have previously all been + * coarsened away. + * * At the end of the function, we call the * respective function to update the number * cache for quads, which will in turn call @@ -562,21 +616,28 @@ namespace */ template void compute_number_cache (const Triangulation &triangulation, - internal::Triangulation::NumberCache<3> &number_cache) + const unsigned int level_objects, + internal::Triangulation::NumberCache<3> &number_cache) { + // update quads, lines and n_levels + compute_number_cache (triangulation, + level_objects, + static_cast&> + (number_cache)); + typedef typename Triangulation::hex_iterator hex_iterator; typedef typename Triangulation::active_hex_iterator active_hex_iterator; - + /////////////////////////////////// // update the number of hexes // on the different levels in // the cache - number_cache.n_hexes_level.resize (triangulation.n_levels()); + number_cache.n_hexes_level.resize (number_cache.n_levels); number_cache.n_hexes = 0; - number_cache.n_active_hexes_level.resize (triangulation.n_levels()); + number_cache.n_active_hexes_level.resize (number_cache.n_levels); number_cache.n_active_hexes = 0; // for 3d, hexes have levels so take @@ -584,25 +645,25 @@ namespace // globally if (dim == 3) { - for (unsigned int level=0; level&>(number_cache)); } @@ -661,7 +717,7 @@ namespace // nothing to do here: the format // hasn't changed for 1d } - + void reorder_compatibility (std::vector > &cells, @@ -670,7 +726,7 @@ namespace for (unsigned int cell=0; cell > &cells, @@ -721,7 +777,7 @@ namespace return line->child(0)->vertex_index(1); return numbers::invalid_unsigned_int; } - + template unsigned int @@ -743,7 +799,7 @@ namespace } return numbers::invalid_unsigned_int; } - + template unsigned int @@ -886,8 +942,8 @@ namespace { return false; } - - + + /** * For a given triangulation: set up the @@ -931,11 +987,11 @@ namespace { // quadrilateral {{0,1}, // face 0, face_orientation = false and true - {1,0}, // face 1, face_orientation = false and true - {1,0}, // face 2, face_orientation = false and true - {0,1}, // face 3, face_orientation = false and true - {0,0}, // face 4, invalid face - {0,0}},// face 5, invalid face + {1,0}, // face 1, face_orientation = false and true + {1,0}, // face 2, face_orientation = false and true + {0,1}, // face 3, face_orientation = false and true + {0,0}, // face 4, invalid face + {0,0}},// face 5, invalid face // hexahedron {{0,1}, {1,0}, @@ -1005,7 +1061,7 @@ namespace (cell->active() || GeometryInfo::face_refinement_case(cell->refinement_case(),f) == RefinementCase::isotropic_refinement)) { - + for (unsigned int c=0; cn_children(); ++c) adjacent_cells[2*face->child(c)->index() + offset] = cell; if (face->child(0)->has_children()) @@ -1019,7 +1075,7 @@ namespace adjacent_cells[2*face->child(1)->child(1)->index() + offset] = cell; } } // if cell active and face refined - } // else -> dim==3 + } // else -> dim==3 } // for all faces of all cells // now loop again over all cells and set the @@ -1139,7 +1195,7 @@ namespace internal /** * Exception * @ingroup Exceptions - */ + */ DeclException0 (ExcInteriorLineCantBeBoundary); /** * Exception @@ -1154,7 +1210,7 @@ namespace internal int, int, << "In SubCellData the line info of the line with vertex indices " << arg1 << " and " << arg2 << " is multiply set."); - + /** * A class into which we put many of the functions that implement @@ -1267,13 +1323,13 @@ namespace internal // dimension, the @p{subcelldata} // field is ignored. (only used for // error checking, which is a good - // idea in any case) + // idea in any case) const unsigned int dim=1; - + // copy vertices triangulation.vertices = v; triangulation.vertices_used = std::vector (v.size(), true); - + // store the indices of the lines // which are adjacent to a given // vertex @@ -1283,15 +1339,15 @@ namespace internal triangulation.levels.push_back (new internal::Triangulation::TriaLevel); triangulation.levels[0]->reserve_space (cells.size(), dim); triangulation.levels[0]->cells.reserve_space (0,cells.size()); - + // make up cells typename Triangulation::raw_line_iterator next_free_line = triangulation.begin_raw_line (); - for (unsigned int cell=0; cellused()) ++next_free_line; - + next_free_line->set (internal::Triangulation ::TriaObject<1> (cells[cell].vertices[0], cells[cell].vertices[1])); @@ -1299,7 +1355,7 @@ namespace internal next_free_line->set_material_id (cells[cell].material_id); next_free_line->clear_user_data (); next_free_line->set_subdomain_id (0); - + // note that this cell is // adjacent to these vertices lines_at_vertex[cells[cell].vertices[0]].push_back (cell); @@ -1311,7 +1367,7 @@ namespace internal { unsigned int boundary_nodes = 0; for (unsigned int i=0; i 1) && (boundary_nodes <= 2)), ExcMessage("The Triangulation has too many end points")); } - + // update neighborship info @@ -1357,7 +1413,7 @@ namespace internal // the second adjacent cell and // vice versa if (lines_at_vertex[line->vertex_index(vertex)][0] == line->index()) - if (lines_at_vertex[line->vertex_index(vertex)].size() == 2) + if (lines_at_vertex[line->vertex_index(vertex)].size() == 2) { const typename Triangulation::cell_iterator neighbor (&triangulation, @@ -1470,7 +1526,7 @@ namespace internal == needed_lines.end(), ExcGridHasInvalidCell(cell)); - + // insert line, with // invalid iterator if line // already exists, then @@ -1486,7 +1542,7 @@ namespace internal std::vector vertex_touch_count (v.size(), 0); typename std::map, typename Triangulation::line_iterator>::iterator i; - for (i=needed_lines.begin(); i!=needed_lines.end(); i++) + for (i=needed_lines.begin(); i!=needed_lines.end(); i++) { // touch the vertices of // this line @@ -1502,7 +1558,7 @@ namespace internal vertex_touch_count.end())) >= 2, ExcGridHasInvalidVertices()); } - + // reserve enough space triangulation.levels.push_back (new internal::Triangulation::TriaLevel); triangulation.faces = new internal::Triangulation::TriaFaces; @@ -1517,7 +1573,7 @@ namespace internal typename std::map, typename Triangulation::line_iterator>::iterator i; for (i = needed_lines.begin(); - line!=triangulation.end_line(); ++line, ++i) + line!=triangulation.end_line(); ++line, ++i) { line->set (internal::Triangulation::TriaObject<1>(i->first.first, i->first.second)); @@ -1546,17 +1602,17 @@ namespace internal lines[line]=needed_lines[std::make_pair( cells[c].vertices[GeometryInfo::line_to_cell_vertices(line, 0)], cells[c].vertices[GeometryInfo::line_to_cell_vertices(line, 1)])]; - + cell->set (internal::Triangulation::TriaObject<2> (lines[0]->index(), lines[1]->index(), lines[2]->index(), lines[3]->index())); - + cell->set_used_flag (); cell->set_material_id (cells[c].material_id); cell->clear_user_data (); cell->set_subdomain_id (0); - + // note that this cell is // adjacent to the four // lines @@ -1568,7 +1624,7 @@ namespace internal for (typename Triangulation::line_iterator line=triangulation.begin_line(); - line!=triangulation.end_line(); ++line) + line!=triangulation.end_line(); ++line) { const unsigned int n_adj_cells = adjacent_cells[line->index()].size(); // assert that every line has @@ -1576,7 +1632,7 @@ namespace internal AssertThrow ((n_adj_cells >= 1) && (n_adj_cells <= 2), ExcInternalError()); - + // if only one cell: line is at // boundary -> give it the // boundary indicator zero by @@ -1594,7 +1650,7 @@ namespace internal = subcelldata.boundary_lines.begin(); std::vector >::const_iterator end_boundary_line = subcelldata.boundary_lines.end(); - for (; boundary_line!=end_boundary_line; ++boundary_line) + for (; boundary_line!=end_boundary_line; ++boundary_line) { typename Triangulation::line_iterator line; std::pair line_vertices(std::make_pair(boundary_line->vertices[0], @@ -1603,14 +1659,14 @@ namespace internal // line found in this // direction line = needed_lines[line_vertices]; - else + else { // look whether it exists // in reverse direction std::swap (line_vertices.first, line_vertices.second); if (needed_lines.find(line_vertices) != needed_lines.end()) line = needed_lines[line_vertices]; - else + else // line does not exist AssertThrow (false, ExcLineInexistant(line_vertices.first, line_vertices.second)); @@ -1622,7 +1678,7 @@ namespace internal line->boundary_indicator() != 255), ExcMultiplySetLineInfoOfLine(line_vertices.first, line_vertices.second)); - + // Assert that only exterior lines // are given a boundary indicator AssertThrow (! (line->boundary_indicator() == 255), @@ -1651,7 +1707,7 @@ namespace internal // we are looking for else cell->set_neighbor (side, - adjacent_cells[cell->line(side)->index()][0]); + adjacent_cells[cell->line(side)->index()][0]); } @@ -1748,7 +1804,7 @@ namespace internal AssertThrow (cells[cell].vertices[vertex] < triangulation.vertices.size(), ExcInvalidVertexIndex (cell, cells[cell].vertices[vertex], triangulation.vertices.size())); - + for (unsigned int line=0; line::lines_per_cell; ++line) { // given a line vertex number @@ -1760,7 +1816,7 @@ namespace internal std::pair line_vertices( cells[cell].vertices[GeometryInfo::line_to_cell_vertices(line, 0)], cells[cell].vertices[GeometryInfo::line_to_cell_vertices(line, 1)]); - + // if that line was already inserted // in reverse order do nothing, else // insert the line @@ -1788,7 +1844,7 @@ namespace internal std::vector vertex_touch_count (v.size(), 0); typename std::map, typename Triangulation::line_iterator>::iterator i; - for (i=needed_lines.begin(); i!=needed_lines.end(); i++) + for (i=needed_lines.begin(); i!=needed_lines.end(); i++) { // touch the vertices of // this line @@ -1821,7 +1877,7 @@ namespace internal line = triangulation.begin_raw_line(); typename std::map, typename Triangulation::line_iterator>::iterator i; - for (i = needed_lines.begin(); line!=triangulation.end_line(); ++line, ++i) + for (i = needed_lines.begin(); line!=triangulation.end_line(); ++line, ++i) { line->set (internal::Triangulation::TriaObject<1>(i->first.first, i->first.second)); @@ -1855,7 +1911,7 @@ namespace internal std_cxx1x::array::lines_per_face> >, QuadComparator> needed_quads; - for (unsigned int cell=0; cell::lines_per_cell]; unsigned int face_line_list[GeometryInfo::lines_per_face]; std_cxx1x::array::lines_per_face> orientation; - + for (unsigned int line=0; line::lines_per_cell; ++line) { line_list[line]=std::pair ( @@ -1894,7 +1950,7 @@ namespace internal cells[cell].vertices[GeometryInfo::line_to_cell_vertices(line, 1)], cells[cell].vertices[GeometryInfo::line_to_cell_vertices(line, 0)]); } - + for (unsigned int face=0; face::faces_per_cell; ++face) { // set up a list of the lines to be @@ -1919,14 +1975,14 @@ namespace internal face_to_cell_lines(face,l)]]->index(); orientation[l]=false; } - - - internal::Triangulation::TriaObject<2> + + + internal::Triangulation::TriaObject<2> quad(face_line_list[0], face_line_list[1], face_line_list[2], face_line_list[3]); - + // insert quad, with // invalid iterator // @@ -1998,7 +2054,7 @@ namespace internal // // first reserve enough space triangulation.faces->quads.reserve_space (0,needed_quads.size()); - + { typename Triangulation::raw_quad_iterator quad = triangulation.begin_raw_quad(); @@ -2018,7 +2074,7 @@ namespace internal quad->set_line_orientation(1,q->second.second[1]); quad->set_line_orientation(2,q->second.second[2]); quad->set_line_orientation(3,q->second.second[3]); - + // now set the iterator for // this quad @@ -2028,7 +2084,7 @@ namespace internal ///////////////////////////////// // finally create the cells - triangulation.levels[0]->cells.reserve_space (cells.size()); + triangulation.levels[0]->cells.reserve_space (cells.size()); // store for each quad index the // adjacent cells @@ -2064,7 +2120,7 @@ namespace internal cells[c].vertices[GeometryInfo::line_to_cell_vertices(line, 1)], cells[c].vertices[GeometryInfo::line_to_cell_vertices(line, 0)]); } - + // get the iterators // corresponding to the // faces. also store @@ -2085,13 +2141,13 @@ namespace internal else face_line_list[l]=needed_lines[inverse_line_list[GeometryInfo:: face_to_cell_lines(face,l)]]->index(); - - internal::Triangulation::TriaObject<2> + + internal::Triangulation::TriaObject<2> quad(face_line_list[0], face_line_list[1], face_line_list[2], face_line_list[3]); - + if (needed_quads.find (quad) != needed_quads.end()) { // face is in standard @@ -2186,10 +2242,10 @@ namespace internal // so something went // wrong above Assert(false,ExcInternalError()); - + } }// for all faces - + // make the cell out of // these iterators cell->set (internal::Triangulation @@ -2199,7 +2255,7 @@ namespace internal face_iterator[3]->index(), face_iterator[4]->index(), face_iterator[5]->index())); - + cell->set_used_flag (); cell->set_material_id (cells[c].material_id); cell->clear_user_flag (); @@ -2214,15 +2270,15 @@ namespace internal cell->set_face_flip (quad, face_flip[quad]); cell->set_face_rotation (quad, face_rotation[quad]); } - - + + // note that this cell is // adjacent to the six // quads for (unsigned int quad=0; quad::faces_per_cell; ++quad) adjacent_cells[face_iterator[quad]->index()].push_back (cell); -#ifdef DEBUG +#ifdef DEBUG // make some checks on the // lines and their // ordering @@ -2247,7 +2303,7 @@ namespace internal std::pair (face,line))); std::multimap >::const_iterator map_iter=cell_to_face_lines.begin(); - + for (; map_iter!=cell_to_face_lines.end(); ++map_iter) { const unsigned int cell_line=map_iter->first; @@ -2338,7 +2394,7 @@ namespace internal if (quad->at_boundary()) for (unsigned int l=0; l<4; ++l) quad->line(l)->set_boundary_indicator (0); - + /////////////////////////////////////// // now set boundary indicators // where given @@ -2357,7 +2413,7 @@ namespace internal // line found in this // direction line = needed_lines[line_vertices]; - + else { // look wether it exists in @@ -2385,11 +2441,11 @@ namespace internal AssertThrow (line->boundary_indicator() == boundary_line->material_id, ExcMessage ("Duplicate boundary lines are only allowed " "if they carry the same boundary indicator.")); - + line->set_boundary_indicator (boundary_line->material_id); } - + // now go on with boundary faces std::vector >::const_iterator boundary_quad = subcelldata.boundary_quads.begin(); @@ -2399,7 +2455,7 @@ namespace internal { typename Triangulation::quad_iterator quad; typename Triangulation::line_iterator line[4]; - + // first find the lines that // are made up of the given // vertices, then build up a @@ -2412,7 +2468,7 @@ namespace internal std::pair line_vertices( boundary_quad->vertices[GeometryInfo::line_to_cell_vertices(i,0)], boundary_quad->vertices[GeometryInfo::line_to_cell_vertices(i,1)]); - + // check whether line // already exists if (needed_lines.find(line_vertices) != needed_lines.end()) @@ -2431,8 +2487,8 @@ namespace internal line_vertices.second)); } } - - + + // Set up 2 quads that are // built up from the lines for // reasons of comparison to @@ -2445,13 +2501,13 @@ namespace internal // and because boundary quad // orientation does not carry // any information. - internal::Triangulation::TriaObject<2> + internal::Triangulation::TriaObject<2> quad_compare_1(line[0]->index(), line[1]->index(), line[2]->index(), line[3]->index()); - internal::Triangulation::TriaObject<2> + internal::Triangulation::TriaObject<2> quad_compare_2(line[2]->index(), line[3]->index(), line[0]->index(), line[1]->index()); - + // try to find the quad with // lines situated as // constructed above. if it @@ -2475,7 +2531,7 @@ namespace internal line_counterclock[lex2cclock[i]]=line[i]; unsigned int n_rotations=0; bool not_found_quad_1; - while ( (not_found_quad_1=(needed_quads.find(quad_compare_1) == needed_quads.end())) && + while ( (not_found_quad_1=(needed_quads.find(quad_compare_1) == needed_quads.end())) && ( needed_quads.find(quad_compare_2) == needed_quads.end()) && (n_rotations<4)) { @@ -2490,10 +2546,10 @@ namespace internal quad_compare_1.set_face(i, line_counterclock[lex2cclock[i]]->index()); quad_compare_2.set_face((i+2)%4, line_counterclock[lex2cclock[i]]->index()); } - + ++n_rotations; } - + AssertThrow (n_rotations!=4, ExcQuadInexistant(line[0]->index(), line[1]->index(), line[2]->index(), line[3]->index())); @@ -2517,8 +2573,8 @@ namespace internal AssertThrow (quad->boundary_indicator() == boundary_quad->material_id, ExcMessage ("Duplicate boundary quads are only allowed " "if they carry the same boundary indicator.")); - - quad->set_boundary_indicator (boundary_quad->material_id); + + quad->set_boundary_indicator (boundary_quad->material_id); } @@ -2558,7 +2614,7 @@ namespace internal Triangulation<1,spacedim> &triangulation) { const unsigned int dim = 1; - + // if spacedim>1 we need to // make sure that we perturb // points but keep them on @@ -2572,7 +2628,7 @@ namespace internal // only difference being the // correction for split faces which // is not necessary in 1D - + // find the smallest length of the // lines adjacent to the // vertex. take the initial value @@ -2585,13 +2641,13 @@ namespace internal for (typename Triangulation::cell_iterator cell=triangulation.begin(0); cell!=triangulation.end(0); ++cell) almost_infinite_length += cell->diameter(); - + std::vector minimal_length (triangulation.vertices.size(), almost_infinite_length); // also note if a vertex is at // the boundary std::vector at_boundary (triangulation.vertices.size(), false); - + for (typename Triangulation::active_line_iterator line=triangulation.begin_active_line(); line != triangulation.end_line(); ++line) @@ -2601,7 +2657,7 @@ namespace internal at_boundary[line->vertex_index(0)] = true; at_boundary[line->vertex_index(1)] = true; } - + minimal_length[line->vertex_index(0)] = std::min(line->diameter(), minimal_length[line->vertex_index(0)]); @@ -2613,8 +2669,8 @@ namespace internal const unsigned int n_vertices = triangulation.vertices.size(); Point shift_vector; - - for (unsigned int vertex=0; vertex::cell_iterator cell=triangulation.begin(0); cell!=triangulation.end(0); ++cell) almost_infinite_length += cell->diameter(); - + std::vector minimal_length (triangulation.vertices.size(), almost_infinite_length); // also note if a vertex is at the // boundary std::vector at_boundary (triangulation.vertices.size(), false); - + for (typename Triangulation::active_line_iterator line=triangulation.begin_active_line(); line != triangulation.end_line(); ++line) @@ -2697,7 +2753,7 @@ namespace internal at_boundary[line->vertex_index(0)] = true; at_boundary[line->vertex_index(1)] = true; } - + minimal_length[line->vertex_index(0)] = std::min(line->diameter(), minimal_length[line->vertex_index(0)]); @@ -2709,8 +2765,8 @@ namespace internal const unsigned int n_vertices = triangulation.vertices.size(); Point shift_vector; - - for (unsigned int vertex=0; vertex::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->has_children() && !cell->face(face)->at_boundary()) @@ -2786,7 +2842,7 @@ namespace internal std::vector &) { const unsigned int dim = 1; - + // first we need to reset the // neighbor pointers of the // neighbors of this cell's @@ -2799,7 +2855,7 @@ namespace internal Assert (!cell->child(0)->has_children() && !cell->child(1)->has_children(), ExcInternalError()); - + // first do it for the cells to the // left if (cell->neighbor(0).state() == IteratorState::valid) @@ -2808,7 +2864,7 @@ namespace internal typename Triangulation::cell_iterator neighbor = cell->neighbor(0); Assert (neighbor->level() == cell->level(), ExcInternalError()); - + // right child neighbor = neighbor->child(1); while (1) @@ -2816,7 +2872,7 @@ namespace internal Assert (neighbor->neighbor(1) == cell->child(0), ExcInternalError()); neighbor->set_neighbor (1, cell); - + // move on to further // children on the // boundary between this @@ -2836,7 +2892,7 @@ namespace internal typename Triangulation::cell_iterator neighbor = cell->neighbor(1); Assert (neighbor->level() == cell->level(), ExcInternalError()); - + // left child neighbor = neighbor->child(0); while (1) @@ -2844,7 +2900,7 @@ namespace internal Assert (neighbor->neighbor(0) == cell->child(1), ExcInternalError()); neighbor->set_neighbor (0, cell); - + // move on to further // children on the // boundary between this @@ -2880,7 +2936,7 @@ namespace internal } - + template static void @@ -2898,9 +2954,9 @@ namespace internal // may be deleted std::vector::line_iterator> lines_to_delete(0); - + lines_to_delete.reserve(4*2+4); - + // now we decrease the counters for // lines contained in the child // cells @@ -2933,12 +2989,12 @@ namespace internal else { unsigned int inner_face_no=ref_case==RefinementCase::cut_x ? 1 : 3; - + // the inner line will not be // used any more lines_to_delete.push_back(cell->child(0)->line(inner_face_no)); } - + // invalidate children for (unsigned int child=0; childn_children(); ++child) { @@ -2972,29 +3028,29 @@ namespace internal (line_cell_count[line->child_index(0)] > 0 && line_cell_count[line->child_index(1)] > 0), ExcInternalError()); - + if (line_cell_count[line->child_index(0)]==0) { for (unsigned int c=0; c<2; ++c) Assert (!line->child(c)->has_children(), ExcInternalError()); - + // we may delete the line's // children and the middle vertex // as no cell references them // anymore triangulation.vertices_used[line->child(0)->vertex_index(1)] = false; - + lines_to_delete.push_back(line->child(0)); lines_to_delete.push_back(line->child(1)); - + line->clear_children(); } } } // finally, delete unneeded lines - + // clear user pointers, to avoid that // they may appear at unwanted places // later on... @@ -3022,7 +3078,7 @@ namespace internal std::vector &quad_cell_count) { const unsigned int dim=3; - + Assert(line_cell_count.size()==triangulation.n_raw_lines(), ExcInternalError()); Assert(quad_cell_count.size()==triangulation.n_raw_quads(), ExcInternalError()); @@ -3035,10 +3091,10 @@ namespace internal lines_to_delete(0); std::vector::quad_iterator> quads_to_delete(0); - + lines_to_delete.reserve(12*2+6*4+6); quads_to_delete.reserve(6*4+12); - + // now we decrease the counters for lines and // quads contained in the child cells for (unsigned int c=0; cn_children(); ++c) @@ -3058,7 +3114,7 @@ namespace internal // // for append quads and lines: only append // them to the list of objects to be deleted - + switch (ref_case) { case RefinementCase::cut_x: @@ -3099,17 +3155,17 @@ namespace internal quads_to_delete.push_back(cell->child(2)->face(1)); quads_to_delete.push_back(cell->child(4)->face(1)); quads_to_delete.push_back(cell->child(6)->face(1)); - + quads_to_delete.push_back(cell->child(0)->face(3)); quads_to_delete.push_back(cell->child(1)->face(3)); quads_to_delete.push_back(cell->child(4)->face(3)); quads_to_delete.push_back(cell->child(5)->face(3)); - + quads_to_delete.push_back(cell->child(0)->face(5)); quads_to_delete.push_back(cell->child(1)->face(5)); quads_to_delete.push_back(cell->child(2)->face(5)); - quads_to_delete.push_back(cell->child(3)->face(5)); - + quads_to_delete.push_back(cell->child(3)->face(5)); + lines_to_delete.push_back(cell->child(0)->line(5)); lines_to_delete.push_back(cell->child(0)->line(7)); lines_to_delete.push_back(cell->child(0)->line(11)); @@ -3129,7 +3185,7 @@ namespace internal Assert(false, ExcInternalError()); break; } - + // invalidate children for (unsigned int child=0; childn_children(); ++child) @@ -3146,7 +3202,7 @@ namespace internal cell->child(child)->set_face_flip(f,false); cell->child(child)->set_face_rotation(f,false); } - + cell->child(child)->clear_used_flag(); } @@ -3163,7 +3219,7 @@ namespace internal // still needing these objects. oherwise we // can delete them. first for quads (and // their inner lines). - + for (unsigned int quad_no=0; quad_no::faces_per_cell; ++quad_no) { typename Triangulation::quad_iterator @@ -3196,7 +3252,7 @@ namespace internal // delete possible grand_children unsigned int deleted_grandchildren=0; unsigned int number_of_child_refinements=0; - + for (unsigned int c=0; c<2; ++c) if (quad->child(c)->has_children()) { @@ -3327,7 +3383,7 @@ namespace internal Assert(!switch_1->has_children(), ExcInternalError()); Assert(!switch_2->has_children(), ExcInternalError()); - + const int switch_1_index=switch_1->index(); const int switch_2_index=switch_2->index(); for (unsigned int l=0; lclear_refinement_case(); quad->set_refinement_case(RefinementCase<2>::cut_xy); quad->set_children(0,child_0); - quad->set_children(2,child_2); + quad->set_children(2,child_2); } else { quad->clear_children(); quad->clear_refinement_case(); } - - + + } break; } @@ -3442,7 +3498,7 @@ namespace internal quad_cell_count[quad->child_index(2)] > 0 && quad_cell_count[quad->child_index(3)] > 0), ExcInternalError()); - + if (quad_cell_count[quad->child_index(0)]==0) { // we may delete the quad's @@ -3456,9 +3512,9 @@ namespace internal for (unsigned int child=0; childn_children(); ++child) quads_to_delete.push_back(quad->child(child)); - + triangulation.vertices_used[quad->child(0)->vertex_index(3)] = false; - + quad->clear_children(); quad->clear_refinement_case(); } @@ -3469,12 +3525,12 @@ namespace internal Assert(false, ExcInternalError()); break; } - + } // now we repeat a similar procedure // for the outer lines of this cell. - + // if in debug mode: check that each // of the lines for which we consider // deleting the children in fact has @@ -3488,7 +3544,7 @@ namespace internal Assert((GeometryInfo::line_refinement_case(ref_case,line_no) && line->has_children()) || GeometryInfo::line_refinement_case(ref_case,line_no)==RefinementCase<1>::no_refinement, ExcInternalError()); - + if (line->has_children()) { // if one of the cell counters is @@ -3499,29 +3555,29 @@ namespace internal (line_cell_count[line->child_index(0)] > 0 && line_cell_count[line->child_index(1)] > 0), ExcInternalError()); - + if (line_cell_count[line->child_index(0)]==0) { for (unsigned int c=0; c<2; ++c) Assert (!line->child(c)->has_children(), ExcInternalError()); - + // we may delete the line's // children and the middle vertex // as no cell references them // anymore triangulation.vertices_used[line->child(0)->vertex_index(1)] = false; - + lines_to_delete.push_back(line->child(0)); lines_to_delete.push_back(line->child(1)); - + line->clear_children(); } } } // finally, delete unneeded quads and lines - + // clear user pointers, to avoid that // they may appear at unwanted places // later on... @@ -3547,10 +3603,10 @@ namespace internal (*quad)->clear_refinement_case(); (*quad)->clear_user_flag(); (*quad)->clear_used_flag(); - } + } } - - + + /** * Create the children of a 2d * cell. The arguments indicate @@ -3586,11 +3642,11 @@ namespace internal are (unlike above) only two possibilities: a neighbor cell is on the same level or one level up (in both cases, it may or may not be refined later on, but we don't care here). - + First: Set up an array of the 3x3 vertices, which are distributed on the cell (the array consists of indices into the @p{vertices} std::vector - + 2--7--3 | | | 4--8--5 @@ -3599,35 +3655,35 @@ namespace internal note: in case of cut_x or cut_y not all these vertices are needed for the new cells - - Second: + + Second: Set up an array of the new lines (the array consists of iterator pointers into the lines arrays) - + .-6-.-7-. The directions are: .->-.->-. 1 9 3 ^ ^ ^ .-10.11-. .->-.->-. 0 8 2 ^ ^ ^ .-4-.-5-. .->-.->-. - cut_x: - .-4-.-5-. - | | | - 0 6 1 - | | | + cut_x: + .-4-.-5-. + | | | + 0 6 1 + | | | .-2-.-3-. - + cut_y: - .---5---. - 1 3 - .---6---. - 0 2 - .---4---. + .---5---. + 1 3 + .---6---. + 0 2 + .---4---. + - Third: Set up an array of neighbors: - + 6 7 .--.--. 1| | |3 @@ -3680,7 +3736,7 @@ namespace internal for (unsigned int line_no=0; line_no<4; ++line_no) if (cell->line(line_no)->has_children()) new_vertices[4+line_no]=cell->line(line_no)->child(0)->vertex_index(1); - + if (ref_case==RefinementCase::cut_xy) { // find the next @@ -3710,13 +3766,13 @@ namespace internal // to the information // stored in the // boundary class - new_point = + new_point = triangulation.boundary[cell->material_id()] ->get_new_point_on_quad (cell); } - + triangulation.vertices[new_vertices[8]] = new_point; - + // if the user_flag is set, i.e. if the // cell is at the boundary, use a @@ -3754,7 +3810,7 @@ namespace internal // do nothing in the following bound_face=GeometryInfo::faces_per_cell+1; } - + if (bound_face::faces_per_cell) // reset the cell's middle vertex // to the middle of the straight @@ -3769,8 +3825,8 @@ namespace internal ->child(0)->vertex(1)); } } - - + + // Now the lines: typename Triangulation::raw_line_iterator new_lines[12]; unsigned int lmin=8; @@ -3780,7 +3836,7 @@ namespace internal lmin=6; lmax=7; } - + for (unsigned int l=lmin; lused() == true) @@ -3791,14 +3847,14 @@ namespace internal Assert (new_lines[l]->used() == false, ExcCellShouldBeUnused()); } - + if (ref_case==RefinementCase::cut_xy) { - // .-6-.-7-. - // 1 9 3 - // .-10.11-. - // 0 8 2 - // .-4-.-5-. + // .-6-.-7-. + // 1 9 3 + // .-10.11-. + // 0 8 2 + // .-4-.-5-. // lines 0-7 already // exist, create only @@ -3821,10 +3877,10 @@ namespace internal } else if (ref_case==RefinementCase::cut_x) { - // .-4-.-5-. - // | | | - // 0 6 1 - // | | | + // .-4-.-5-. + // | | | + // 0 6 1 + // | | | // .-2-.-3-. new_lines[0]=cell->line(0); new_lines[1]=cell->line(1); @@ -3838,10 +3894,10 @@ namespace internal else { Assert(ref_case==RefinementCase::cut_y, ExcInternalError()); - // .---5---. - // 1 3 - // .---6---. - // 0 2 + // .---5---. + // 1 3 + // .---6---. + // 0 2 // .---4---. new_lines[0]=cell->line(0)->child(0); new_lines[1]=cell->line(0)->child(1); @@ -3869,10 +3925,10 @@ namespace internal subcells[GeometryInfo::max_children_per_cell]; while (next_unused_cell->used() == true) ++next_unused_cell; - + const unsigned int n_children= GeometryInfo::n_children(ref_case); - for (unsigned int i=0; iused() == false, ExcCellShouldBeUnused()); @@ -3892,11 +3948,11 @@ namespace internal // |0 | 1| // .--.--. // lines: - // .-6-.-7-. - // 1 9 3 - // .-10.11-. - // 0 8 2 - // .-4-.-5-. + // .-6-.-7-. + // 1 9 3 + // .-10.11-. + // 0 8 2 + // .-4-.-5-. subcells[0]->set (internal::Triangulation:: TriaObject<2>(new_lines[0]->index(), new_lines[8]->index(), @@ -3927,10 +3983,10 @@ namespace internal // | | | // .--.--. // lines: - // .-4-.-5-. - // | | | - // 0 6 1 - // | | | + // .-4-.-5-. + // | | | + // 0 6 1 + // | | | // .-2-.-3-. subcells[0]->set (internal::Triangulation:: TriaObject<2>(new_lines[0]->index(), @@ -3953,10 +4009,10 @@ namespace internal // | 0 | // .-----. // lines: - // .---5---. - // 1 3 - // .---6---. - // 0 2 + // .---5---. + // 1 3 + // .---6---. + // 0 2 // .---4---. subcells[0]->set (internal::Triangulation:: TriaObject<2>(new_lines[0]->index(), @@ -3969,17 +4025,17 @@ namespace internal new_lines[6]->index(), new_lines[5]->index())); } - + for (unsigned int i=0; iset_used_flag(); - subcells[i]->clear_refine_flag(); + subcells[i]->clear_refine_flag(); subcells[i]->clear_user_flag(); subcells[i]->clear_user_data(); subcells[i]->clear_children(); // inherit material - // properties + // properties subcells[i]->set_material_id (cell->material_id()); subcells[i]->set_subdomain_id (cell->subdomain_id()); } @@ -3990,15 +4046,15 @@ namespace internal for (unsigned int i=0; iset_children (2*i, subcells[2*i]->index()); // set the refine case - cell->set_refinement_case(ref_case); + cell->set_refinement_case(ref_case); // note that the // refinement flag was // already cleared at the - // beginning of this function + // beginning of this function } - - + + /** * A function that performs the @@ -4010,7 +4066,7 @@ namespace internal execute_refinement (Triangulation<1,spacedim> &triangulation) { const unsigned int dim = 1; - + // check whether a new level is // needed we have to check for this // on the highest level only (on @@ -4023,7 +4079,7 @@ namespace internal endc = triangulation.end(); for (; cell != endc; ++cell) if (cell->used()) - if (cell->refine_flag_set()) + if (cell->refine_flag_set()) { triangulation.levels .push_back (new internal::Triangulation::TriaLevel); @@ -4049,7 +4105,7 @@ namespace internal typename Triangulation::active_cell_iterator acell = triangulation.begin_active(level), aendc = triangulation.begin_active(level+1); - for (; acell!=aendc; ++acell) + for (; acell!=aendc; ++acell) if (acell->refine_flag_set()) ++flagged_cells; @@ -4078,7 +4134,7 @@ namespace internal .reserve_space (GeometryInfo<1>::max_children_per_cell * flagged_cells, 0); - + needed_vertices += flagged_cells; } @@ -4108,18 +4164,18 @@ namespace internal // index of next unused vertex unsigned int next_unused_vertex = 0; - - for (int level=triangulation.levels.size()-2; level>=0; --level) + + for (int level=triangulation.levels.size()-2; level>=0; --level) { typename Triangulation::active_cell_iterator cell = triangulation.begin_active(level), endc = triangulation.begin_active(level+1); - + typename Triangulation::raw_cell_iterator next_unused_cell = triangulation.begin_raw (level+1); - for (; (cell!=endc) && (cell->level()==level); ++cell) - if (cell->refine_flag_set()) + for (; (cell!=endc) && (cell->level()==level); ++cell) + if (cell->refine_flag_set()) { // clear refinement flag cell->clear_refine_flag (); @@ -4130,7 +4186,7 @@ namespace internal ++next_unused_vertex; Assert (next_unused_vertex < triangulation.vertices.size(), ExcTooFewVerticesAllocated()); - + // first insert new // vertex. if dim==spacedim // then simply use the @@ -4141,7 +4197,7 @@ namespace internal triangulation.vertices[next_unused_vertex] = (cell->vertex(0) + cell->vertex(1)) / 2; else - triangulation.vertices[next_unused_vertex] = + triangulation.vertices[next_unused_vertex] = triangulation.boundary[cell->material_id()] ->get_new_point_on_line(cell); triangulation.vertices_used[next_unused_vertex] = true; @@ -4172,8 +4228,8 @@ namespace internal ::TriaObject<1> (cell->vertex_index(0), next_unused_vertex)); first_child->set_material_id (cell->material_id()); - first_child->set_subdomain_id (cell->subdomain_id()); - + first_child->set_subdomain_id (cell->subdomain_id()); + // reset neighborship info (refer // to // internal::Triangulation::TriaLevel<0> @@ -4221,7 +4277,7 @@ namespace internal left_neighbor->set_neighbor (1, first_child); } } - + // insert second child second_child->clear_children (); second_child->set (internal::Triangulation @@ -4229,7 +4285,7 @@ namespace internal cell->vertex_index(1))); second_child->set_neighbor (0, first_child); second_child->set_material_id (cell->material_id()); - second_child->set_subdomain_id (cell->subdomain_id()); + second_child->set_subdomain_id (cell->subdomain_id()); if (cell->neighbor(1).state() != IteratorState::valid) second_child->set_neighbor (1, cell->neighbor(1)); else @@ -4245,7 +4301,7 @@ namespace internal // above { second_child->set_neighbor (1, cell->neighbor(1)->child(0)); - + typename Triangulation::cell_iterator right_neighbor = cell->neighbor(1); while (right_neighbor->has_children()) @@ -4254,7 +4310,7 @@ namespace internal right_neighbor->set_neighbor (0, second_child); } } - } + } } // in 1d, we can not have @@ -4280,7 +4336,7 @@ namespace internal execute_refinement (Triangulation<2,spacedim> &triangulation) { const unsigned int dim = 2; - + // check whether a new level is // needed we have to check for this // on the highest level only (on @@ -4294,7 +4350,7 @@ namespace internal endc = triangulation.end(); for (; cell != endc; ++cell) if (cell->used()) - if (cell->refine_flag_set()) + if (cell->refine_flag_set()) { triangulation.levels.push_back (new internal::Triangulation::TriaLevel); break; @@ -4345,13 +4401,13 @@ namespace internal typename Triangulation::active_cell_iterator cell = triangulation.begin_active(level), endc = triangulation.begin_active(level+1); - for (; cell!=endc; ++cell) - if (cell->refine_flag_set()) + for (; cell!=endc; ++cell) + if (cell->refine_flag_set()) { if (cell->refine_flag_set()==RefinementCase::cut_xy) { needed_cells += 4; - + // new vertex at // center of cell is // needed in any case @@ -4367,17 +4423,17 @@ namespace internal // anisotropic refinement is // used for at least one cell triangulation.anisotropic_refinement = true; - + needed_cells += 2; // no vertex at center - + // the inner line can // be stored as // single n_single_lines += 1; } - + // mark all faces // (lines) for // refinement; @@ -4406,7 +4462,7 @@ namespace internal line->set_user_flag (); //TODO[WB]: we overwrite the user_index here because we later on need to find // out which boundary object we have to ask to refine this line. find a better -// way to do that +// way to do that if (spacedim > dim) line->set_user_index(cell->material_id()); } @@ -4414,14 +4470,14 @@ namespace internal } } - + // count number of used cells // on the next higher level const unsigned int used_cells = std::count_if (triangulation.levels[level+1]->cells.used.begin(), triangulation.levels[level+1]->cells.used.end(), std::bind2nd (std::equal_to(), true)); - + // reserve space for the // used_cells cells already @@ -4462,7 +4518,7 @@ namespace internal // of the current cells triangulation.faces->lines. reserve_space (n_lines_in_pairs, 0); - + // add to needed vertices how many // vertices are already in use needed_vertices += std::count_if (triangulation.vertices_used.begin(), triangulation.vertices_used.end(), @@ -4480,7 +4536,7 @@ namespace internal } - // Do REFINEMENT + // Do REFINEMENT // on every level; exclude highest // level as above @@ -4539,7 +4595,7 @@ namespace internal // answer triangulation.vertices[next_unused_vertex] = triangulation.boundary[line->user_index()]->get_new_point_on_line (line); - + // now that we created // the right point, make // up the two child @@ -4569,7 +4625,7 @@ namespace internal // pointer of the present // line line->set_children (0, next_unused_line->index()); - + // set the two new lines const typename Triangulation::raw_line_iterator children[2] = { next_unused_line, @@ -4581,14 +4637,14 @@ namespace internal // will fail Assert (children[0]->used() == false, ExcCellShouldBeUnused()); Assert (children[1]->used() == false, ExcCellShouldBeUnused()); - + children[0]->set (internal::Triangulation ::TriaObject<1>(line->vertex_index(0), next_unused_vertex)); children[1]->set (internal::Triangulation ::TriaObject<1>(next_unused_vertex, line->vertex_index(1))); - + children[0]->set_used_flag(); children[1]->set_used_flag(); children[0]->clear_children(); @@ -4600,7 +4656,7 @@ namespace internal children[0]->set_boundary_indicator (line->boundary_indicator()); children[1]->set_boundary_indicator (line->boundary_indicator()); - + // finally clear flag // indicating the need // for refinement @@ -4608,7 +4664,7 @@ namespace internal } } - + // Now set up the new cells // reserve space for inner @@ -4619,40 +4675,40 @@ namespace internal typename Triangulation<2,spacedim>::DistortedCellList cells_with_distorted_children; - + // reset next_unused_line, as // now also single empty places // in the vector can be used typename Triangulation::raw_line_iterator next_unused_line = triangulation.begin_raw_line (); - - for (int level=0; level(triangulation.levels.size())-1; ++level) + + for (int level=0; level(triangulation.levels.size())-1; ++level) { - + // Remember: as we don't operate // on the finest level, begin_*(level+1) // is allowed typename Triangulation::active_cell_iterator cell = triangulation.begin_active(level), endc = triangulation.begin_active(level+1); - + typename Triangulation::raw_cell_iterator next_unused_cell = triangulation.begin_raw (level+1); - for (; cell!=endc; ++cell) - if (cell->refine_flag_set()) + for (; cell!=endc; ++cell) + if (cell->refine_flag_set()) { // set the user flag to // indicate, that at least one // line is at the boundary - + // TODO[Tobias Leicht] find a // better place to set this flag, // so that we do not need so much // time to check each cell here if (cell->at_boundary()) cell->set_user_flag(); - + // actually set up the children and // update neighbor information create_children (triangulation, @@ -4706,7 +4762,7 @@ namespace internal endc = triangulation.end(); for (; cell != endc; ++cell) if (cell->used()) - if (cell->refine_flag_set()) + if (cell->refine_flag_set()) { triangulation.levels.push_back (new internal::Triangulation::TriaLevel); break; @@ -4719,14 +4775,14 @@ namespace internal // to flag which lines and quads // need refinement triangulation.faces->quads.clear_user_data(); - + for (typename Triangulation::line_iterator line=triangulation.begin_line(); line!=triangulation.end_line(); ++line) line->clear_user_flag(); for (typename Triangulation::quad_iterator quad=triangulation.begin_quad(); quad!=triangulation.end_quad(); ++quad) quad->clear_user_flag(); - + // create an array of face refine cases. User // indices of faces will be set to values // corresponding with indices in this array. @@ -4735,7 +4791,7 @@ namespace internal RefinementCase::cut_x, RefinementCase::cut_y, RefinementCase::cut_xy}; - + // check how much space is needed // on every level // we need not check the highest @@ -4767,15 +4823,15 @@ namespace internal // vertices and new lines will // be needed unsigned int new_cells = 0; - + typename Triangulation::active_cell_iterator acell = triangulation.begin_active(level), aendc = triangulation.begin_active(level+1); - for (; acell!=aendc; ++acell) - if (acell->refine_flag_set()) + for (; acell!=aendc; ++acell) + if (acell->refine_flag_set()) { RefinementCase ref_case=acell->refine_flag_set(); - + // now for interior vertices, lines // and quads, which are needed in // any case @@ -4808,7 +4864,7 @@ namespace internal // we should never get here Assert(false, ExcInternalError()); } - + // mark all faces for refinement; // checking locally // if and how the neighbor @@ -4871,16 +4927,16 @@ namespace internal } } }// for all faces - + // flag all lines, that have to be // refined for (unsigned int line=0; line::lines_per_cell; ++line) if (GeometryInfo::line_refinement_case(ref_case,line) && !acell->line(line)->has_children()) acell->line(line)->set_user_flag(); - + }// if refine_flag set and for all cells on this level - + // count number of used cells on // the next higher level @@ -4956,7 +5012,7 @@ namespace internal } } } - + for (typename Triangulation::line_iterator line=triangulation.begin_line(); line!=triangulation.end_line(); ++line) if (line->user_flag_set()) @@ -5066,7 +5122,7 @@ namespace internal else triangulation.vertices[next_unused_vertex] = (line->vertex(0) + line->vertex(1)) / 2; - + // now that we created // the right point, make // up the two child lines @@ -5075,7 +5131,7 @@ namespace internal next_unused_line=triangulation.faces->lines.next_free_pair_line(triangulation); Assert(next_unused_line.state() == IteratorState::valid, ExcInternalError()); - + // now we found // two consecutive unused // lines, such that the @@ -5085,12 +5141,12 @@ namespace internal // pointer of the present // line line->set_children (0, next_unused_line->index()); - + // set the two new lines const typename Triangulation::raw_line_iterator children[2] = { next_unused_line, ++next_unused_line }; - + // some tests; if any of // the iterators should // be invalid, then @@ -5098,14 +5154,14 @@ namespace internal // will fail Assert (children[0]->used() == false, ExcCellShouldBeUnused()); Assert (children[1]->used() == false, ExcCellShouldBeUnused()); - + children[0]->set (internal::Triangulation ::TriaObject<1>(line->vertex_index(0), next_unused_vertex)); children[1]->set (internal::Triangulation ::TriaObject<1>(next_unused_vertex, line->vertex_index(1))); - + children[0]->set_used_flag(); children[1]->set_used_flag(); children[0]->clear_children(); @@ -5117,7 +5173,7 @@ namespace internal children[0]->set_boundary_indicator (line->boundary_indicator()); children[1]->set_boundary_indicator (line->boundary_indicator()); - + // finally clear flag // indicating the need // for refinement @@ -5190,11 +5246,11 @@ namespace internal // processed and refined. if (aniso_quad_ref_case == quad->refinement_case()) continue; - + Assert(quad->refinement_case()==RefinementCase::cut_xy || quad->refinement_case()==RefinementCase::no_refinement, ExcInternalError()); - + // this quad needs to be refined // anisotropically Assert(quad->user_index() == RefinementCase::cut_x || @@ -5204,7 +5260,7 @@ namespace internal // make the new line interior to // the quad typename Triangulation::raw_line_iterator new_line; - + new_line=triangulation.faces->lines.next_free_single_line(triangulation); Assert (new_line->used() == false, ExcCellShouldBeUnused()); @@ -5254,16 +5310,16 @@ namespace internal const unsigned int index[2][2]= {{1,0}, // child 0, line_orientation=false and true {0,1}}; // child 1, line_orientation=false and true - + // find some space (consecutive) // for the two newly to be // created quads. typename Triangulation::raw_quad_iterator new_quads[2]; - + next_unused_quad=triangulation.faces->quads.next_free_pair_quad(triangulation); new_quads[0] = next_unused_quad; Assert (new_quads[0]->used() == false, ExcCellShouldBeUnused()); - + ++next_unused_quad; new_quads[1] = next_unused_quad; Assert (new_quads[1]->used() == false, ExcCellShouldBeUnused()); @@ -5295,7 +5351,7 @@ namespace internal new_line->index(), quad->line_index(3))); } - + for (unsigned int i=0; i<2; ++i) { new_quads[i]->set_used_flag(); @@ -5375,11 +5431,11 @@ namespace internal else { Assert(aniso_quad_ref_case==RefinementCase::cut_y, ExcInternalError()); - + old_child[0]=quad->child(0)->line(3); old_child[1]=quad->child(1)->line(3); } - + if (old_child[0]->index()+1 != old_child[1]->index()) { // this is exactly the @@ -5389,13 +5445,13 @@ namespace internal // two new lines and copy // all info typename Triangulation::raw_line_iterator new_child[2]; - + new_child[0]=new_child[1]=triangulation.faces->lines.next_free_pair_line(triangulation); ++new_child[1]; new_child[0]->set_used_flag(); new_child[1]->set_used_flag(); - + const int old_index_0=old_child[0]->index(), old_index_1=old_child[1]->index(), new_index_0=new_child[0]->index(), @@ -5419,7 +5475,7 @@ namespace internal for (unsigned int i=0; i<2; ++i) { Assert(!old_child[i]->has_children(), ExcInternalError()); - + new_child[i]->set(internal::Triangulation::TriaObject<1>(old_child[i]->vertex_index(0), old_child[i]->vertex_index(1))); new_child[i]->set_boundary_indicator(old_child[i]->boundary_indicator()); @@ -5428,9 +5484,9 @@ namespace internal new_child[i]->set_user_flag(); else new_child[i]->clear_user_flag(); - + new_child[i]->clear_children(); - + old_child[i]->clear_user_flag(); old_child[i]->clear_user_index(); old_child[i]->clear_used_flag(); @@ -5545,7 +5601,7 @@ namespace internal switch_1->set_children(0, switch_2->child_index(0)); if (switch_2->refinement_case()==RefinementCase::cut_xy) switch_1->set_children(2, switch_2->child_index(2)); - + switch_2->set(internal::Triangulation::TriaObject<2>(switch_1_lines[0], switch_1_lines[1], switch_1_lines[2], @@ -5583,13 +5639,13 @@ namespace internal } quad->clear_children(); } - + // note these quads as children // to the present one quad->set_children (0, new_quads[0]->index()); quad->set_refinement_case(aniso_quad_ref_case); - + // finally clear flag // indicating the need // for refinement @@ -5612,7 +5668,7 @@ namespace internal ++next_unused_vertex; Assert (next_unused_vertex < triangulation.vertices.size(), ExcTooFewVerticesAllocated()); - + // now: if the quad is refined // anisotropically already, set // the anisotropic refinement @@ -5622,7 +5678,7 @@ namespace internal // as it is an outer line of the // two (anisotropic) children const RefinementCase quad_ref_case=quad->refinement_case(); - + if (quad_ref_case==RefinementCase::cut_x || quad_ref_case==RefinementCase::cut_y) { @@ -5656,21 +5712,21 @@ namespace internal triangulation.vertices[next_unused_vertex] = (middle_line->vertex(0) + middle_line->vertex(1)) / 2; triangulation.vertices_used[next_unused_vertex] = true; - + // now search a slot for the two // child lines next_unused_line=triangulation.faces->lines.next_free_pair_line(triangulation); - + // set the child // pointer of the present // line middle_line->set_children (0, next_unused_line->index()); - + // set the two new lines const typename Triangulation::raw_line_iterator children[2] = { next_unused_line, ++next_unused_line }; - + // some tests; if any of // the iterators should // be invalid, then @@ -5678,14 +5734,14 @@ namespace internal // will fail Assert (children[0]->used() == false, ExcCellShouldBeUnused()); Assert (children[1]->used() == false, ExcCellShouldBeUnused()); - + children[0]->set (internal::Triangulation:: TriaObject<1>(middle_line->vertex_index(0), next_unused_vertex)); children[1]->set (internal::Triangulation:: TriaObject<1>(next_unused_vertex, middle_line->vertex_index(1))); - + children[0]->set_used_flag(); children[1]->set_used_flag(); children[0]->clear_children(); @@ -5694,7 +5750,7 @@ namespace internal children[1]->clear_user_data(); children[0]->clear_user_flag(); children[1]->clear_user_flag(); - + children[0]->set_boundary_indicator (middle_line->boundary_indicator()); children[1]->set_boundary_indicator (middle_line->boundary_indicator()); } @@ -5716,7 +5772,7 @@ namespace internal // set the middle vertex // appropriately - if (quad->at_boundary()) + if (quad->at_boundary()) triangulation.vertices[next_unused_vertex] = triangulation.boundary[quad->boundary_indicator()]->get_new_point_on_quad (quad); else @@ -5821,15 +5877,15 @@ namespace internal // *2-*-3* // | 0 | // *--*--* - + const unsigned int vertex_indices[5] = { quad->line(0)->child(0)->vertex_index(1), quad->line(1)->child(0)->vertex_index(1), quad->line(2)->child(0)->vertex_index(1), quad->line(3)->child(0)->vertex_index(1), - next_unused_vertex + next_unused_vertex }; - + new_lines[0]->set (internal::Triangulation:: TriaObject<1>(vertex_indices[2], vertex_indices[4])); new_lines[1]->set (internal::Triangulation:: @@ -5847,7 +5903,7 @@ namespace internal new_lines[i]->clear_children(); new_lines[i]->set_boundary_indicator(quad->boundary_indicator()); } - + // now for the // quads. again, first // collect some data @@ -5871,7 +5927,7 @@ namespace internal const unsigned int index[2][2]= {{1,0}, // child 0, line_orientation=false and true {0,1}}; // child 1, line_orientation=false and true - + const unsigned int line_indices[12] = { quad->line(0)->child(index[0][quad->line_orientation(0)])->index(), quad->line(0)->child(index[1][quad->line_orientation(0)])->index(), @@ -5884,16 +5940,16 @@ namespace internal new_lines[0]->index(), new_lines[1]->index(), new_lines[2]->index(), - new_lines[3]->index() + new_lines[3]->index() }; - + // find some space (consecutive) // for the first two newly to be // created quads. typename Triangulation::raw_quad_iterator new_quads[4]; next_unused_quad=triangulation.faces->quads.next_free_pair_quad(triangulation); - + new_quads[0] = next_unused_quad; Assert (new_quads[0]->used() == false, ExcCellShouldBeUnused()); @@ -5919,7 +5975,7 @@ namespace internal line_indices[8], line_indices[4], line_indices[10])); - + quad->set_refinement_case(RefinementCase<2>::cut_xy); new_quads[0]->set (internal::Triangulation @@ -5971,7 +6027,7 @@ namespace internal new_quads[2]->set_line_orientation(3,quad->line_orientation(3)); new_quads[3]->set_line_orientation(1,quad->line_orientation(1)); new_quads[3]->set_line_orientation(3,quad->line_orientation(3)); - + // finally clear flag // indicating the need // for refinement @@ -6010,7 +6066,7 @@ namespace internal { // this hex needs to be // refined - + // clear flag indicating // the need for // refinement. do it here @@ -6061,7 +6117,7 @@ namespace internal Assert(false, ExcInternalError()); break; } - + // find some space for the newly to // be created interior lines and // initialize them. @@ -6070,7 +6126,7 @@ namespace internal for (unsigned int i=0; ilines.next_free_single_line(triangulation); - + Assert (new_lines[i]->used() == false, ExcCellShouldBeUnused()); new_lines[i]->set_used_flag(); @@ -6089,7 +6145,7 @@ namespace internal for (unsigned int i=0; iquads.next_free_single_quad(triangulation); - + Assert (new_quads[i]->used() == false, ExcCellShouldBeUnused()); new_quads[i]->set_used_flag(); @@ -6117,9 +6173,9 @@ namespace internal next_unused_hex=triangulation.levels[level+1]->cells.next_free_hex(triangulation,level+1); else ++next_unused_hex; - + new_hexes[i]=next_unused_hex; - + Assert (new_hexes[i]->used() == false, ExcCellShouldBeUnused()); new_hexes[i]->set_used_flag(); @@ -6156,7 +6212,7 @@ namespace internal // note these hexes as // children to the // present cell - for (unsigned int i=0; iset_children (2*i, new_hexes[2*i]->index()); // we have to take into account @@ -6173,7 +6229,7 @@ namespace internal hex->face_orientation (3), hex->face_orientation (4), hex->face_orientation (5) }; - + // face_flip const bool f_fl[6] = { hex->face_flip (0), @@ -6182,7 +6238,7 @@ namespace internal hex->face_flip (3), hex->face_flip (4), hex->face_flip (5) }; - + // face_rotation const bool f_ro[6] = { hex->face_rotation (0), @@ -6217,7 +6273,7 @@ namespace internal { 1, 1 }}, // RefinementCase::cut_x, face_flip=true, face_rotation=false and true { { 0, 1 }, // RefinementCase::cut_y, face_flip=false, face_rotation=false and true { 1, 0 }}};// RefinementCase::cut_y, face_flip=true, face_rotation=false and true - + /////////////////////////////////////// // // in the following we will do the @@ -6261,13 +6317,13 @@ namespace internal // | | | / // | | |/ // *----*----* - // + // // again, first // collect some data // about the indices of // the lines, with the // following numbering: - + // face 2: front plane // (note: x,y exchanged) // *---*---* @@ -6280,9 +6336,9 @@ namespace internal // (note: x,y exchanged) // m1 // *---*---* - // | | | + // | | | // | 1 | - // | | | + // | | | // *---*---* // face 4: bottom plane // *---*---* @@ -6315,11 +6371,11 @@ namespace internal hex->face(5)->child(0) ->line((hex->face(5)->refinement_case() == RefinementCase<2>::cut_x) ? 1 : 3) //3 }; - + lines=&lines_x[0]; - + unsigned int line_indices_x[4]; - + for (unsigned int i=0; i<4; ++i) line_indices_x[i]=lines[i]->index(); line_indices=&line_indices_x[0]; @@ -6358,7 +6414,7 @@ namespace internal ExcInternalError()); line_orientation_x[i]=false; } - + line_orientation=&line_orientation_x[0]; // set up the new quad, line @@ -6374,7 +6430,7 @@ namespace internal new_quads[0]->set_line_orientation(1,line_orientation[1]); new_quads[0]->set_line_orientation(2,line_orientation[2]); new_quads[0]->set_line_orientation(3,line_orientation[3]); - + // the quads are // numbered as follows: // @@ -6382,7 +6438,7 @@ namespace internal // of the old hex: // * // /| - // / | x + // / | x // / | *-------* *---------* // * | | | / / // | 0 | | | / / @@ -6412,23 +6468,23 @@ namespace internal const unsigned int quad_indices_x[11] = { new_quads[0]->index(), //0 - + hex->face(0)->index(), //1 - + hex->face(1)->index(), //2 hex->face(2)->child_index( child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]), //3 hex->face(2)->child_index(1-child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]), - + hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), //5 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), - + hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), //7 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), - + hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]), //9 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]) - + }; quad_indices=&quad_indices_x[0]; @@ -6468,37 +6524,37 @@ namespace internal // | | * // | |/ // *---------* - // + // // again, first // collect some data // about the indices of // the lines, with the // following numbering: - + // face 0: left plane // * - // /| + // /| // * | // /| | // * | | // | 0 | // | | * - // | |/ + // | |/ // | *m0 - // |/ - // * + // |/ + // * // face 1: right plane // * - // /| + // /| // m1* | // /| | // * | | // | 1 | // | | * - // | |/ - // | * - // |/ - // * + // | |/ + // | * + // |/ + // * // face 4: bottom plane // *-------* // / / @@ -6528,11 +6584,11 @@ namespace internal hex->face(5)->child(0) ->line((hex->face(5)->refinement_case() == RefinementCase<2>::cut_x) ? 1 : 3) //3 }; - + lines=&lines_y[0]; - + unsigned int line_indices_y[4]; - + for (unsigned int i=0; i<4; ++i) line_indices_y[i]=lines[i]->index(); line_indices=&line_indices_y[0]; @@ -6570,7 +6626,7 @@ namespace internal ExcInternalError()); line_orientation_y[i]=false; } - + line_orientation=&line_orientation_y[0]; // set up the new quad, line @@ -6594,7 +6650,7 @@ namespace internal // of the old hex: // * // /| - // / | x + // / | x // / | *-------* *---------* // * | | | / / // | | | 0 | / / @@ -6624,23 +6680,23 @@ namespace internal const unsigned int quad_indices_y[11] = { new_quads[0]->index(), //0 - + hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), //1 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), - + hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), //3 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), - + hex->face(2)->index(), //5 - + hex->face(3)->index(), //6 hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), //7 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), - + hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]), //9 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]) - + }; quad_indices=&quad_indices_y[0]; @@ -6680,37 +6736,37 @@ namespace internal // | | / // | |/ // *---------* - // + // // again, first // collect some data // about the indices of // the lines, with the // following numbering: - + // face 0: left plane // * - // /| + // /| // / | // / * // * /| // | 0 | // |/ * - // m0* / - // | / - // |/ - // * + // m0* / + // | / + // |/ + // * // face 1: right plane // * - // /| + // /| // / | // / *m1 // * /| // | 1 | // |/ * - // * / - // | / - // |/ - // * + // * / + // | / + // |/ + // * // face 2: front plane // (note: x,y exchanged) // *-------* @@ -6721,9 +6777,9 @@ namespace internal // face 3: back plane // (note: x,y exchanged) // *-------* - // | | + // | | // *---3---*m1 - // | | + // | | // *-------* // set up a list of line iterators @@ -6742,11 +6798,11 @@ namespace internal hex->face(3)->child(0) ->line((hex->face(3)->refinement_case() == RefinementCase<2>::cut_x) ? 1 : 3) //3 }; - + lines=&lines_z[0]; - + unsigned int line_indices_z[4]; - + for (unsigned int i=0; i<4; ++i) line_indices_z[i]=lines[i]->index(); line_indices=&line_indices_z[0]; @@ -6784,7 +6840,7 @@ namespace internal ExcInternalError()); line_orientation_z[i]=false; } - + line_orientation=&line_orientation_z[0]; // set up the new quad, line @@ -6795,7 +6851,7 @@ namespace internal line_indices[1], line_indices[2], line_indices[3])); - + new_quads[0]->set_line_orientation(0,line_orientation[0]); new_quads[0]->set_line_orientation(1,line_orientation[1]); new_quads[0]->set_line_orientation(2,line_orientation[2]); @@ -6808,7 +6864,7 @@ namespace internal // of the old hex: // * // /| - // / | x + // / | x // / | *-------* *---------* // * | | | / / // | | | | / 0 / @@ -6839,21 +6895,21 @@ namespace internal const unsigned int quad_indices_z[11] = { new_quads[0]->index(), //0 - + hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), //1 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), - + hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), //3 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), - + hex->face(2)->child_index( child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]), //5 hex->face(2)->child_index(1-child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]), - + hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), //7 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), - + hex->face(4)->index(), //9 - + hex->face(5)->index() //10 }; quad_indices=&quad_indices_z[0]; @@ -6894,7 +6950,7 @@ namespace internal // | | | * // | | |/ // *----*----* - // + // // first, create the new // internal line @@ -6907,31 +6963,31 @@ namespace internal // about the indices of // the lines, with the // following numbering: - + // face 0: left plane // * - // /| + // /| // * | // /| | // * | | // | 0 | // | | * - // | |/ - // | * - // |/ - // * + // | |/ + // | * + // |/ + // * // face 1: right plane // * - // /| + // /| // * | // /| | // * | | // | 1 | // | | * - // | |/ - // | * - // |/ - // * + // | |/ + // | * + // |/ + // * // face 2: front plane // (note: x,y exchanged) // *---*---* @@ -6942,9 +6998,9 @@ namespace internal // face 3: back plane // (note: x,y exchanged) // *---*---* - // | | | + // | | | // | 3 | - // | | | + // | | | // *---*---* // face 4: bottom plane // *---*---* @@ -6989,7 +7045,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[4],f_fl[4],f_ro[4])), //6 hex->face(4)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[4],f_fl[4],f_ro[4])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[4],f_fl[4],f_ro[4])), //7 - + hex->face(5)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[5],f_fl[5],f_ro[5])) ->line(GeometryInfo::standard_to_real_face_line(1,f_or[5],f_fl[5],f_ro[5])), //8 hex->face(5)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[5],f_fl[5],f_ro[5])) @@ -7001,11 +7057,11 @@ namespace internal new_lines[0] //12 }; - + lines=&lines_xy[0]; - + unsigned int line_indices_xy[13]; - + for (unsigned int i=0; i<13; ++i) line_indices_xy[i]=lines[i]->index(); line_indices=&line_indices_xy[0]; @@ -7044,7 +7100,7 @@ namespace internal ExcInternalError()); line_orientation_xy[i]=false; } - + // note: for lines 4 to 11 // (inner lines of the outer quads) // the following holds: the second @@ -7084,11 +7140,11 @@ namespace internal // * * // /| 9| // * | * | - // y/| | 8| 3 - // * |1| * | | + // y/| | 8| 3 + // * |1| * | | // | | |x | 12| - // |0| * | | * - // | |/ 2 |5 + // |0| * | | * + // | |/ 2 |5 // | * | * // |/ |4 // * * @@ -7096,10 +7152,10 @@ namespace internal // x // *---*---* *10-*-11* // | | | | | | - // | 2 | 3 | 0 12 1 + // | 2 | 3 | 0 12 1 // | | | | | | // *---*---*y *-6-*-7-* - + new_quads[0]->set (internal::Triangulation ::TriaObject<2>(line_indices[2], line_indices[12], @@ -7124,11 +7180,11 @@ namespace internal new_quads[0]->set_line_orientation(0,line_orientation[2]); new_quads[0]->set_line_orientation(2,line_orientation[4]); new_quads[0]->set_line_orientation(3,line_orientation[8]); - + new_quads[1]->set_line_orientation(1,line_orientation[3]); new_quads[1]->set_line_orientation(2,line_orientation[5]); new_quads[1]->set_line_orientation(3,line_orientation[9]); - + new_quads[2]->set_line_orientation(0,line_orientation[6]); new_quads[2]->set_line_orientation(1,line_orientation[10]); new_quads[2]->set_line_orientation(2,line_orientation[0]); @@ -7144,7 +7200,7 @@ namespace internal // of the old hex: // * // /| - // * | x + // * | x // /| | *---*---* *---------* // * |1| | | | / / // | | | | 2 | 3 | / / @@ -7177,24 +7233,24 @@ namespace internal new_quads[1]->index(), new_quads[2]->index(), new_quads[3]->index(), - + hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), //4 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), //6 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), - + hex->face(2)->child_index( child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]), //8 hex->face(2)->child_index(1-child_at_origin[hex->face(2)->refinement_case()-1][f_fl[2]][f_ro[2]]), - + hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), //10 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), - + hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[4],f_fl[4],f_ro[4])), //12 hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[4],f_fl[4],f_ro[4])), hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[4],f_fl[4],f_ro[4])), hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[4],f_fl[4],f_ro[4])), - + hex->face(5)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[5],f_fl[5],f_ro[5])), //16 hex->face(5)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[5],f_fl[5],f_ro[5])), hex->face(5)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[5],f_fl[5],f_ro[5])), @@ -7252,7 +7308,7 @@ namespace internal // | | | / // | | |/ // *----*----* - // + // // first, create the new // internal line @@ -7265,31 +7321,31 @@ namespace internal // about the indices of // the lines, with the // following numbering: - + // face 0: left plane // * - // /| + // /| // / | // / * // * /| // | 0 | // |/ * - // * / - // | / - // |/ - // * + // * / + // | / + // |/ + // * // face 1: right plane // * - // /| + // /| // / | // / * // * /| // | 1 | // |/ * - // * / - // | / - // |/ - // * + // * / + // | / + // |/ + // * // face 2: front plane // (note: x,y exchanged) // *---*---* @@ -7300,9 +7356,9 @@ namespace internal // face 3: back plane // (note: x,y exchanged) // *---*---* - // | 9 | + // | 9 | // *10-*-11* - // | 8 | + // | 8 | // *---*---* // face 4: bottom plane // *---*---* @@ -7347,7 +7403,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(1,f_or[2],f_fl[2],f_ro[2])), //6 hex->face(2)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[2],f_fl[2],f_ro[2])) ->line(GeometryInfo::standard_to_real_face_line(0,f_or[2],f_fl[2],f_ro[2])), //7 - + hex->face(3)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[3],f_fl[3],f_ro[3])) ->line(GeometryInfo::standard_to_real_face_line(3,f_or[3],f_fl[3],f_ro[3])), //8 hex->face(3)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[3],f_fl[3],f_ro[3])) @@ -7359,11 +7415,11 @@ namespace internal new_lines[0] //12 }; - + lines=&lines_xz[0]; - + unsigned int line_indices_xz[13]; - + for (unsigned int i=0; i<13; ++i) line_indices_xz[i]=lines[i]->index(); line_indices=&line_indices_xz[0]; @@ -7402,7 +7458,7 @@ namespace internal ExcInternalError()); line_orientation_xz[i]=false; } - + // note: for lines 4 to 11 // (inner lines of the outer quads) // the following holds: the second @@ -7443,14 +7499,14 @@ namespace internal // // * * // /| /| - // / | 3 9 - // y/ * / * + // / | 3 9 + // y/ * / * // * 3/| * /| // | / |x 5 12|8 - // |/ * |/ * - // * 2/ * / + // |/ * |/ * + // * 2/ * / // | / 4 2 - // |/ |/ + // |/ |/ // * * // // y @@ -7459,7 +7515,7 @@ namespace internal // / 0 / 1 / 0 12 1 // / / / / / / // *----*----*x *--6-*--7-* - + new_quads[0]->set (internal::Triangulation ::TriaObject<2>(line_indices[0], line_indices[12], @@ -7484,11 +7540,11 @@ namespace internal new_quads[0]->set_line_orientation(0,line_orientation[0]); new_quads[0]->set_line_orientation(2,line_orientation[6]); new_quads[0]->set_line_orientation(3,line_orientation[10]); - + new_quads[1]->set_line_orientation(1,line_orientation[1]); new_quads[1]->set_line_orientation(2,line_orientation[7]); new_quads[1]->set_line_orientation(3,line_orientation[11]); - + new_quads[2]->set_line_orientation(0,line_orientation[4]); new_quads[2]->set_line_orientation(1,line_orientation[8]); new_quads[2]->set_line_orientation(2,line_orientation[2]); @@ -7504,7 +7560,7 @@ namespace internal // of the old hex: // * // /| - // / | x + // / | x // /3 * *-------* *----*----* // * /| | | / / / // | / | | | / 0 / 1 / @@ -7538,18 +7594,18 @@ namespace internal new_quads[1]->index(), new_quads[2]->index(), new_quads[3]->index(), - + hex->face(0)->child_index( child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), //4 hex->face(0)->child_index(1-child_at_origin[hex->face(0)->refinement_case()-1][f_fl[0]][f_ro[0]]), hex->face(1)->child_index( child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), //6 hex->face(1)->child_index(1-child_at_origin[hex->face(1)->refinement_case()-1][f_fl[1]][f_ro[1]]), - - hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[2],f_fl[2],f_ro[2])), //8 + + hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[2],f_fl[2],f_ro[2])), //8 hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[2],f_fl[2],f_ro[2])), hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[2],f_fl[2],f_ro[2])), hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[2],f_fl[2],f_ro[2])), - + hex->face(3)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[3],f_fl[3],f_ro[3])), //12 hex->face(3)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[3],f_fl[3],f_ro[3])), hex->face(3)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[3],f_fl[3],f_ro[3])), @@ -7557,7 +7613,7 @@ namespace internal hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), //16 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), - + hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]), //18 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]) }; @@ -7623,7 +7679,7 @@ namespace internal // | | * // | |/ // *---------* - // + // // first, create the new // internal line @@ -7640,10 +7696,10 @@ namespace internal // 1 each are shown twice // for better // readability) - + // face 0: left plane // * * - // /| /| + // /| /| // * | * | // /| * /| * // * 5/| * |7| @@ -7651,11 +7707,11 @@ namespace internal // |/| * |6| * // * 4/ * |/ // | * | * - // |/ |/ + // |/ |/ // * * // face 1: right plane // * * - // /| /| + // /| /| // * | * | // /| * /| * // * 9/| * |11 @@ -7675,9 +7731,9 @@ namespace internal // face 3: back plane // (note: x,y exchanged) // *-------* - // | | + // | | // *---1---* - // | | + // | | // *-------* // face 4: bottom plane // *-------* @@ -7722,7 +7778,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[0],f_fl[0],f_ro[0])), //6 hex->face(0)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[0],f_fl[0],f_ro[0])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[0],f_fl[0],f_ro[0])), //7 - + hex->face(1)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[1],f_fl[1],f_ro[1])) ->line(GeometryInfo::standard_to_real_face_line(1,f_or[1],f_fl[1],f_ro[1])), //8 hex->face(1)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[1],f_fl[1],f_ro[1])) @@ -7734,11 +7790,11 @@ namespace internal new_lines[0] //12 }; - + lines=&lines_yz[0]; - + unsigned int line_indices_yz[13]; - + for (unsigned int i=0; i<13; ++i) line_indices_yz[i]=lines[i]->index(); line_indices=&line_indices_yz[0]; @@ -7777,7 +7833,7 @@ namespace internal ExcInternalError()); line_orientation_yz[i]=false; } - + // note: for lines 4 to 11 // (inner lines of the outer quads) // the following holds: the second @@ -7827,7 +7883,7 @@ namespace internal // *---------* *----12---* // / 0 / 6 10 // *---------*x *----0----* - + new_quads[0]->set (internal::Triangulation ::TriaObject<2>(line_indices[6], line_indices[10], @@ -7852,11 +7908,11 @@ namespace internal new_quads[0]->set_line_orientation(0,line_orientation[6]); new_quads[0]->set_line_orientation(1,line_orientation[10]); new_quads[0]->set_line_orientation(2,line_orientation[0]); - + new_quads[1]->set_line_orientation(0,line_orientation[7]); new_quads[1]->set_line_orientation(1,line_orientation[11]); new_quads[1]->set_line_orientation(3,line_orientation[1]); - + new_quads[2]->set_line_orientation(0,line_orientation[2]); new_quads[2]->set_line_orientation(2,line_orientation[4]); new_quads[2]->set_line_orientation(3,line_orientation[8]); @@ -7872,7 +7928,7 @@ namespace internal // of the old hex: // * // /| - // / | x + // / | x // / | *-------* *---------* // * | | 3 | / 1 / // | | *-------* *---------* @@ -7906,13 +7962,13 @@ namespace internal new_quads[1]->index(), new_quads[2]->index(), new_quads[3]->index(), - - hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[0],f_fl[0],f_ro[0])), //4 + + hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[0],f_fl[0],f_ro[0])), //4 hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[0],f_fl[0],f_ro[0])), hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[0],f_fl[0],f_ro[0])), hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[0],f_fl[0],f_ro[0])), - - hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[1],f_fl[1],f_ro[1])), //8 + + hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[1],f_fl[1],f_ro[1])), //8 hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[1],f_fl[1],f_ro[1])), hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[1],f_fl[1],f_ro[1])), hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[1],f_fl[1],f_ro[1])), @@ -7922,10 +7978,10 @@ namespace internal hex->face(3)->child_index( child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), //14 hex->face(3)->child_index(1-child_at_origin[hex->face(3)->refinement_case()-1][f_fl[3]][f_ro[3]]), - + hex->face(4)->child_index( child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), //16 hex->face(4)->child_index(1-child_at_origin[hex->face(4)->refinement_case()-1][f_fl[4]][f_ro[4]]), - + hex->face(5)->child_index( child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]), //18 hex->face(5)->child_index(1-child_at_origin[hex->face(5)->refinement_case()-1][f_fl[5]][f_ro[5]]) }; @@ -7982,7 +8038,7 @@ namespace internal // | | | * // | | |/ // *----*----* - // + // // find the next unused // vertex and set it @@ -7992,7 +8048,7 @@ namespace internal Assert (next_unused_vertex < triangulation.vertices.size(), ExcTooFewVerticesAllocated()); triangulation.vertices_used[next_unused_vertex] = true; - + // the new vertex is // definitely in the // interior, so we need not @@ -8059,10 +8115,10 @@ namespace internal middle_vertex_index(hex->face(3)), middle_vertex_index(hex->face(4)), middle_vertex_index(hex->face(5)), - next_unused_vertex + next_unused_vertex }; vertex_indices=&vertex_indices_xyz[0]; - + new_lines[0]->set (internal::Triangulation:: TriaObject<1>(vertex_indices[2], vertex_indices[6])); new_lines[1]->set (internal::Triangulation:: @@ -8085,10 +8141,10 @@ namespace internal // 1 each are shown twice // for better // readability) - + // face 0: left plane // * * - // /| /| + // /| /| // * | * | // /| * /| * // * 1/| * |3| @@ -8096,11 +8152,11 @@ namespace internal // |/| * |2| * // * 0/ * |/ // | * | * - // |/ |/ + // |/ |/ // * * // face 1: right plane // * * - // /| /| + // /| /| // * | * | // /| * /| * // * 5/| * |7| @@ -8120,9 +8176,9 @@ namespace internal // face 3: back plane // (note: x,y exchanged) // *---*---* - // | 15 | + // | 15 | // *12-*-13* - // | 14 | + // | 14 | // *---*---* // face 4: bottom plane // *---*---* @@ -8142,7 +8198,7 @@ namespace internal // *26-*-27* *26-*-27* // / 24 / | 28 | // *---*---* *---*---* - + // set up a list of line iterators // first. from this, construct // lists of line_indices and @@ -8158,7 +8214,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[0],f_fl[0],f_ro[0])), //2 hex->face(0)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[0],f_fl[0],f_ro[0])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[0],f_fl[0],f_ro[0])), //3 - + hex->face(1)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[1],f_fl[1],f_ro[1])) ->line(GeometryInfo::standard_to_real_face_line(1,f_or[1],f_fl[1],f_ro[1])), //4 hex->face(1)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[1],f_fl[1],f_ro[1])) @@ -8167,7 +8223,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[1],f_fl[1],f_ro[1])), //6 hex->face(1)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[1],f_fl[1],f_ro[1])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[1],f_fl[1],f_ro[1])), //7 - + hex->face(2)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[2],f_fl[2],f_ro[2])) ->line(GeometryInfo::standard_to_real_face_line(1,f_or[2],f_fl[2],f_ro[2])), //8 hex->face(2)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[2],f_fl[2],f_ro[2])) @@ -8176,7 +8232,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[2],f_fl[2],f_ro[2])), //10 hex->face(2)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[2],f_fl[2],f_ro[2])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[2],f_fl[2],f_ro[2])), //11 - + hex->face(3)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[3],f_fl[3],f_ro[3])) ->line(GeometryInfo::standard_to_real_face_line(1,f_or[3],f_fl[3],f_ro[3])), //12 hex->face(3)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[3],f_fl[3],f_ro[3])) @@ -8185,7 +8241,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[3],f_fl[3],f_ro[3])), //14 hex->face(3)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[3],f_fl[3],f_ro[3])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[3],f_fl[3],f_ro[3])), //15 - + hex->face(4)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[4],f_fl[4],f_ro[4])) ->line(GeometryInfo::standard_to_real_face_line(1,f_or[4],f_fl[4],f_ro[4])), //16 hex->face(4)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[4],f_fl[4],f_ro[4])) @@ -8194,7 +8250,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[4],f_fl[4],f_ro[4])), //18 hex->face(4)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[4],f_fl[4],f_ro[4])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[4],f_fl[4],f_ro[4])), //19 - + hex->face(5)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(0,f_or[5],f_fl[5],f_ro[5])) ->line(GeometryInfo::standard_to_real_face_line(1,f_or[5],f_fl[5],f_ro[5])), //20 hex->face(5)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[5],f_fl[5],f_ro[5])) @@ -8203,7 +8259,7 @@ namespace internal ->line(GeometryInfo::standard_to_real_face_line(3,f_or[5],f_fl[5],f_ro[5])), //22 hex->face(5)->isotropic_child(GeometryInfo::standard_to_real_face_vertex(3,f_or[5],f_fl[5],f_ro[5])) ->line(GeometryInfo::standard_to_real_face_line(2,f_or[5],f_fl[5],f_ro[5])), //23 - + new_lines[0], //24 new_lines[1], //25 new_lines[2], //26 @@ -8211,14 +8267,14 @@ namespace internal new_lines[4], //28 new_lines[5] //29 }; - + lines=&lines_xyz[0]; - + unsigned int line_indices_xyz[30]; for (unsigned int i=0; i<30; ++i) line_indices_xyz[i]=lines[i]->index(); line_indices=&line_indices_xyz[0]; - + // the orientation of lines for the // inner quads is quite tricky. as // these lines are newly created @@ -8228,7 +8284,7 @@ namespace internal // fill it with the respective // values bool line_orientation_xyz[30]; - + // note: for the first 24 lines // (inner lines of the outer quads) // the following holds: the second @@ -8275,11 +8331,11 @@ namespace internal // * * // /| 21| // * | * 15 - // y/|3* 20| * - // * |/| * |/| + // y/|3* 20| * + // * |/| * |/| // |2* |x 11 * 14 - // |/|1* |/| * - // * |/ * |17 + // |/|1* |/| * + // * |/ * |17 // |0* 10 * // |/ |16 // * * @@ -8297,7 +8353,7 @@ namespace internal // *----*----* *-26-*-27-* // / 8 / 9 / 2 24 6 // *----*----*x *--8-*--9-* - + new_quads[0]->set (internal::Triangulation ::TriaObject<2>(line_indices[10], line_indices[28], @@ -8358,44 +8414,44 @@ namespace internal line_indices[7], line_indices[27], line_indices[13])); - + // now reset the line_orientation // flags of outer lines as they // cannot be set in a loop (at // least not easily) new_quads[0]->set_line_orientation(0,line_orientation[10]); new_quads[0]->set_line_orientation(2,line_orientation[16]); - + new_quads[1]->set_line_orientation(1,line_orientation[14]); new_quads[1]->set_line_orientation(2,line_orientation[17]); - + new_quads[2]->set_line_orientation(0,line_orientation[11]); new_quads[2]->set_line_orientation(3,line_orientation[20]); new_quads[3]->set_line_orientation(1,line_orientation[15]); new_quads[3]->set_line_orientation(3,line_orientation[21]); - + new_quads[4]->set_line_orientation(0,line_orientation[18]); new_quads[4]->set_line_orientation(2,line_orientation[0]); - + new_quads[5]->set_line_orientation(1,line_orientation[22]); new_quads[5]->set_line_orientation(2,line_orientation[1]); - + new_quads[6]->set_line_orientation(0,line_orientation[19]); new_quads[6]->set_line_orientation(3,line_orientation[4]); - + new_quads[7]->set_line_orientation(1,line_orientation[23]); new_quads[7]->set_line_orientation(3,line_orientation[5]); - + new_quads[8]->set_line_orientation(0,line_orientation[2]); new_quads[8]->set_line_orientation(2,line_orientation[8]); - + new_quads[9]->set_line_orientation(1,line_orientation[6]); new_quads[9]->set_line_orientation(2,line_orientation[9]); - + new_quads[10]->set_line_orientation(0,line_orientation[3]); new_quads[10]->set_line_orientation(3,line_orientation[12]); - + new_quads[11]->set_line_orientation(1,line_orientation[7]); new_quads[11]->set_line_orientation(3,line_orientation[13]); @@ -8407,7 +8463,7 @@ namespace internal // need the indices of a // whole lotta // quads. - + // the quads are // numbered as follows: // @@ -8456,7 +8512,7 @@ namespace internal new_quads[9]->index(), new_quads[10]->index(), new_quads[11]->index(), //11 - + hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[0],f_fl[0],f_ro[0])), //12 hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[0],f_fl[0],f_ro[0])), hex->face(0)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[0],f_fl[0],f_ro[0])), @@ -8466,22 +8522,22 @@ namespace internal hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[1],f_fl[1],f_ro[1])), hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[1],f_fl[1],f_ro[1])), hex->face(1)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[1],f_fl[1],f_ro[1])), - + hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[2],f_fl[2],f_ro[2])), //20 hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[2],f_fl[2],f_ro[2])), hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[2],f_fl[2],f_ro[2])), hex->face(2)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[2],f_fl[2],f_ro[2])), - + hex->face(3)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[3],f_fl[3],f_ro[3])), //24 hex->face(3)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[3],f_fl[3],f_ro[3])), hex->face(3)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[3],f_fl[3],f_ro[3])), hex->face(3)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[3],f_fl[3],f_ro[3])), - + hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[4],f_fl[4],f_ro[4])), //28 hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[4],f_fl[4],f_ro[4])), hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[4],f_fl[4],f_ro[4])), hex->face(4)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(3,f_or[4],f_fl[4],f_ro[4])), - + hex->face(5)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(0,f_or[5],f_fl[5],f_ro[5])), //32 hex->face(5)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(1,f_or[5],f_fl[5],f_ro[5])), hex->face(5)->isotropic_child_index(GeometryInfo::standard_to_real_face_vertex(2,f_or[5],f_fl[5],f_ro[5])), @@ -8518,7 +8574,7 @@ namespace internal quad_indices[26], quad_indices[31], quad_indices[11])); - + // top children new_hexes[4]->set (internal::Triangulation ::TriaObject<3>(quad_indices[14], @@ -8565,7 +8621,7 @@ namespace internal Assert(false, ExcInternalError()); break; }//switch (ref_case) - + // and set face orientation // flags. note that new faces in // the interior of the mother cell @@ -8584,7 +8640,7 @@ namespace internal // flags. however, the latter two // are set to false by default as // this is the standard value - + // loop over all faces and all // (relevant) subfaces of that in // order to set the correct values @@ -8628,7 +8684,7 @@ namespace internal internal::int2type(), internal::int2type())) cells_with_distorted_children.distorted_cells.push_back (hex); - + // note that the // refinement flag was // already cleared at the @@ -8719,13 +8775,13 @@ namespace internal = StaticMappingQ1::mapping. transform_real_to_unit_cell(cell, new_bound); - + // Now, we have to // calculate the // distance from the // face in the unit // cell. - + // take the correct // coordinate direction (0 // for faces 0 and 1, 1 for @@ -8743,14 +8799,14 @@ namespace internal // face for isotropic // refinement const double allowed=0.25; - + if (dist>allowed) cell->flag_for_face_refinement(face_no); }//if flagged for anistropic refinement }//if (cell->face(face)->at_boundary()) }//for all cells } - + /** * Some dimension dependent stuff for @@ -8786,7 +8842,7 @@ namespace internal prepare_refinement_dim_dependent (Triangulation<3,spacedim> &triangulation) { const unsigned int dim = 3; - + // first clear flags on lines, // since we need them to determine // which lines will be refined @@ -8813,7 +8869,7 @@ namespace internal // sets these flags: triangulation.fix_coarsen_flags (); - + // flag those lines that are refined and // will not be coarsened and those that // will be refined @@ -8872,12 +8928,12 @@ namespace internal // is not already // flagged) bool offending_line_found = false; - + for (unsigned int c=0; c<2; ++c) { Assert (cell->line(line)->child(c)->has_children() == false, ExcInternalError()); - + if (cell->line(line)->child(c)->user_flag_set () && (GeometryInfo::line_refinement_case(cell->refine_flag_set(), line) @@ -8899,7 +8955,7 @@ namespace internal cell->flag_for_line_refinement(line); else cell->set_refine_flag(); - + for (unsigned int l=0; l::lines_per_cell; ++l) if (GeometryInfo::line_refinement_case(cell->refine_flag_set(), line) ==RefinementCase<1>::cut_x) @@ -8912,7 +8968,7 @@ namespace internal // changed // the grid offending_line_found = true; - + // it may save us several // loop iterations if we // flag all lines of @@ -8927,7 +8983,7 @@ namespace internal line) !=RefinementCase<1>::no_refinement)) cell->line(line)->set_user_flag(); - + break; } } @@ -8939,7 +8995,7 @@ namespace internal } } } - + // there is another thing here: // if any of the lines will be @@ -8990,11 +9046,11 @@ namespace internal /** * Helper function for - * @p fix_coarsen_flags. Return wether + * @p fix_coarsen_flags. Return wether * coarsening of this cell is allowed. * Coarsening can be forbidden if the * neighboring cells are or will be - * refined twice along the common face. + * refined twice along the common face. */ template static @@ -9007,7 +9063,7 @@ namespace internal // constraint there if (dim == 1) return true; - + const RefinementCase ref_case = cell->refinement_case(); for (unsigned int n=0; n::faces_per_cell; ++n) { @@ -9015,13 +9071,13 @@ namespace internal // along that face, coarsening // will not change anything, so // do nothing. the same applies, - // if the face is at the boandary + // if the face is at the boandary const RefinementCase face_ref_case = GeometryInfo::face_refinement_case(cell->refinement_case(), n); - + const unsigned int n_subfaces = GeometryInfo::n_children(face_ref_case); - + if (n_subfaces == 0 || cell->at_boundary(n)) continue; for (unsigned int c=0; cchild(GeometryInfo:: child_cell_on_face(ref_case, n,c)); - + const typename Triangulation::cell_iterator child_neighbor = child->neighbor(n); if (!child->neighbor_is_coarser(n)) @@ -9128,7 +9184,7 @@ Triangulation::~Triangulation () template -void Triangulation::clear () +void Triangulation::clear () { clear_despite_subscriptions(); } @@ -9151,7 +9207,7 @@ Triangulation::set_boundary (const unsigned int number, const Boundary& boundary_object) { Assert(number<255, ExcIndexRange(number,0,255)); - + boundary[number] = &boundary_object; } @@ -9168,22 +9224,22 @@ Triangulation::set_boundary (const unsigned int number) template const Boundary & -Triangulation::get_boundary (const unsigned int number) const +Triangulation::get_boundary (const unsigned int number) const { Assert(number<255, ExcIndexRange(number,0,255)); - + return *(boundary[number]); } template std::vector -Triangulation::get_boundary_indicators () const +Triangulation::get_boundary_indicators () const { // in 1d, the situation is simpler: if (dim == 1) return std::vector (0);; - + std::vector bi_exists(255, false); active_cell_iterator cell=begin_active(); for (; cell!=end(); ++cell) @@ -9218,14 +9274,14 @@ void Triangulation::copy_triangulation (const Triangulation(*old_tria.faces); for (unsigned i=0;i<255;++i) @@ -9236,7 +9292,7 @@ void Triangulation::copy_triangulation (const Triangulation(*old_tria.levels[level])); number_cache = old_tria.number_cache; - + // inform RefinementListeners of old_tria of // the copy operation typename std::list::iterator ref_listener = @@ -9244,7 +9300,7 @@ void Triangulation::copy_triangulation (const Triangulationcopy_notification (old_tria, *this); - + // note that we need not copy the // subscriptor! } @@ -9283,7 +9339,7 @@ create_triangulation (const std::vector > &v, // check that no forbidden arrays // are used Assert (subcelldata.check_consistency(dim), ExcInternalError()); - + // try to create a triangulation; // if this fails, we still want to // throw an exception but if we @@ -9300,8 +9356,8 @@ create_triangulation (const std::vector > &v, clear_despite_subscriptions(); throw; } - - compute_number_cache (*this, number_cache); + + compute_number_cache (*this, levels.size(), number_cache); // now verify that there are indeed // no distorted cells. as per the @@ -9404,14 +9460,14 @@ template void Triangulation::load_refine_flags (const std::vector &v) { AssertThrow (v.size() == dim*n_active_cells(), ExcGridReadError()); - + active_cell_iterator cell = begin_active(), endc = end(); std::vector::const_iterator i = v.begin(); for (; cell!=endc; ++cell) { unsigned int ref_case=0; - + for(unsigned int j=0; j::load_refine_flags (const std::vector &v cell->clear_refine_flag(); } - Assert (i == v.end(), ExcInternalError()); + Assert (i == v.end(), ExcInternalError()); } @@ -9469,7 +9525,7 @@ template void Triangulation::load_coarsen_flags (const std::vector &v) { Assert (v.size() == n_active_cells(), ExcGridReadError()); - + active_cell_iterator cell = begin_active(), endc = end(); std::vector::const_iterator i = v.begin(); @@ -9494,7 +9550,7 @@ bool Triangulation::get_anisotropic_refinement_flag() const /*-------------------- user data/flags -------------------------*/ -namespace +namespace { // clear user data of cells template @@ -9509,20 +9565,20 @@ namespace void clear_user_data (internal::Triangulation::TriaFaces<1> *) { // nothing to do in 1d - } + } void clear_user_data (internal::Triangulation::TriaFaces<2> *faces) { faces->lines.clear_user_data(); } - + void clear_user_data (internal::Triangulation::TriaFaces<3> *faces) { faces->lines.clear_user_data(); faces->quads.clear_user_data(); - } + } } @@ -9542,7 +9598,7 @@ namespace internal::Triangulation::TriaFaces<1> *) { for (unsigned int level=0; levelcells.clear_user_flags(); + levels[level]->cells.clear_user_flags(); } template @@ -9574,9 +9630,9 @@ namespace internal::Triangulation::TriaFaces<2> *) { for (unsigned int level=0; levelcells.clear_user_flags(); + levels[level]->cells.clear_user_flags(); } - + template void clear_user_flags_quad (std::vector*> &, internal::Triangulation::TriaFaces *faces) @@ -9608,12 +9664,12 @@ namespace { // nothing to do in 2d } - + void clear_user_flags_hex (std::vector*> &levels, internal::Triangulation::TriaFaces<3> *) { for (unsigned int level=0; levelcells.clear_user_flags(); + levels[level]->cells.clear_user_flags(); } } @@ -9649,10 +9705,10 @@ template void Triangulation::save_user_flags (std::ostream &out) const { save_user_flags_line (out); - + if (dim>=2) save_user_flags_quad (out); - + if (dim>=3) save_user_flags_hex (out); @@ -9679,12 +9735,12 @@ void Triangulation::save_user_flags (std::vector &v) const save_user_flags_quad (tmp); v.insert (v.end(), tmp.begin(), tmp.end()); } - + if (dim >= 3) { save_user_flags_hex (tmp); v.insert (v.end(), tmp.begin(), tmp.end()); - } + } if (dim >= 4) Assert (false, ExcNotImplemented()); @@ -9696,10 +9752,10 @@ template void Triangulation::load_user_flags (std::istream &in) { load_user_flags_line (in); - + if (dim>=2) load_user_flags_quad (in); - + if (dim>=3) load_user_flags_hex (in); @@ -9729,14 +9785,14 @@ void Triangulation::load_user_flags (const std::vector &v) v.begin()+n_lines(), v.begin()+n_lines()+n_quads()); load_user_flags_quad (tmp); } - + if (dim >= 3) { tmp.clear(); tmp.insert (tmp.end(), v.begin()+n_lines()+n_quads(), v.begin()+n_lines()+n_quads()+n_hexs()); load_user_flags_hex (tmp); - } + } if (dim >= 4) Assert (false, ExcNotImplemented()); @@ -9785,7 +9841,7 @@ template void Triangulation::load_user_flags_line (const std::vector &v) { Assert (v.size() == n_lines(), ExcGridReadError()); - + line_iterator line = begin_line(), endl = end_line(); std::vector::const_iterator i = v.begin(); @@ -9897,7 +9953,7 @@ void Triangulation::load_user_flags_quad (const std::vector Assert (v.size() == n_quads(), ExcGridReadError()); if (dim >= 2) - { + { quad_iterator quad = begin_quad(), endq = end_quad(); std::vector::const_iterator i = v.begin(); @@ -9919,7 +9975,7 @@ void Triangulation::save_user_flags_hex (std::vector &v) co v.resize (n_hexs(), false); if (dim >= 3) - { + { std::vector::iterator i = v.begin(); hex_iterator hex = begin_hex(), endh = end_hex(); @@ -9958,9 +10014,9 @@ template void Triangulation::load_user_flags_hex (const std::vector &v) { Assert (v.size() == n_hexs(), ExcGridReadError()); - + if (dim >= 3) - { + { hex_iterator hex = begin_hex(), endh = end_hex(); std::vector::const_iterator i = v.begin(); @@ -9970,7 +10026,7 @@ void Triangulation::load_user_flags_hex (const std::vector else clear_user_flag(hex); - Assert (i == v.end(), ExcInternalError()); + Assert (i == v.end(), ExcInternalError()); } } @@ -9993,12 +10049,12 @@ void Triangulation::save_user_indices (std::vector save_user_indices_quad (tmp); v.insert (v.end(), tmp.begin(), tmp.end()); } - + if (dim >= 3) { save_user_indices_hex (tmp); v.insert (v.end(), tmp.begin(), tmp.end()); - } + } if (dim >= 4) Assert (false, ExcNotImplemented()); @@ -10026,14 +10082,14 @@ void Triangulation::load_user_indices (const std::vector= 3) { tmp.clear (); tmp.insert (tmp.end(), v.begin()+n_lines()+n_quads(), v.begin()+n_lines()+n_quads()+n_hexs()); load_user_indices_hex (tmp); - } + } if (dim >= 4) Assert (false, ExcNotImplemented()); @@ -10095,7 +10151,7 @@ template void Triangulation::load_user_indices_line (const std::vector &v) { Assert (v.size() == n_lines(), ExcGridReadError()); - + line_iterator line = begin_line(), endl = end_line(); std::vector::const_iterator i = v.begin(); @@ -10110,7 +10166,7 @@ void Triangulation::save_user_indices_quad (std::vector= 2) - { + { std::vector::iterator i = v.begin(); quad_iterator quad = begin_quad(), endq = end_quad(); @@ -10125,9 +10181,9 @@ template void Triangulation::load_user_indices_quad (const std::vector &v) { Assert (v.size() == n_quads(), ExcGridReadError()); - + if (dim >= 2) - { + { quad_iterator quad = begin_quad(), endq = end_quad(); std::vector::const_iterator i = v.begin(); @@ -10143,7 +10199,7 @@ void Triangulation::save_user_indices_hex (std::vector= 3) - { + { std::vector::iterator i = v.begin(); hex_iterator hex = begin_hex(), endh = end_hex(); @@ -10158,9 +10214,9 @@ template void Triangulation::load_user_indices_hex (const std::vector &v) { Assert (v.size() == n_hexs(), ExcGridReadError()); - + if (dim >= 3) - { + { hex_iterator hex = begin_hex(), endh = end_hex(); std::vector::const_iterator i = v.begin(); @@ -10228,12 +10284,12 @@ void Triangulation::save_user_pointers (std::vector &v) c save_user_pointers_quad (tmp); v.insert (v.end(), tmp.begin(), tmp.end()); } - + if (dim >= 3) { save_user_pointers_hex (tmp); v.insert (v.end(), tmp.begin(), tmp.end()); - } + } if (dim >= 4) Assert (false, ExcNotImplemented()); @@ -10261,14 +10317,14 @@ void Triangulation::load_user_pointers (const std::vector v.begin()+n_lines(), v.begin()+n_lines()+n_quads()); load_user_pointers_quad (tmp); } - + if (dim >= 3) { tmp.clear (); tmp.insert (tmp.end(), v.begin()+n_lines()+n_quads(), v.begin()+n_lines()+n_quads()+n_hexs()); load_user_pointers_hex (tmp); - } + } if (dim >= 4) Assert (false, ExcNotImplemented()); @@ -10293,7 +10349,7 @@ template void Triangulation::load_user_pointers_line (const std::vector &v) { Assert (v.size() == n_lines(), ExcGridReadError()); - + line_iterator line = begin_line(), endl = end_line(); std::vector::const_iterator i = v.begin(); @@ -10309,7 +10365,7 @@ void Triangulation::save_user_pointers_quad (std::vector v.resize (n_quads(), 0); if (dim >= 2) - { + { std::vector::iterator i = v.begin(); quad_iterator quad = begin_quad(), endq = end_quad(); @@ -10324,7 +10380,7 @@ template void Triangulation::load_user_pointers_quad (const std::vector &v) { Assert (v.size() == n_quads(), ExcGridReadError()); - + if (dim >= 2) { quad_iterator quad = begin_quad(), @@ -10357,7 +10413,7 @@ template void Triangulation::load_user_pointers_hex (const std::vector &v) { Assert (v.size() == n_hexs(), ExcGridReadError()); - + if (dim >= 3) { hex_iterator hex = begin_hex(), @@ -10992,7 +11048,7 @@ Triangulation::begin_raw_quad (const unsigned int level) const if (levels[level]->cells.cells.size() == 0) return end_quad(); - + return raw_quad_iterator (const_cast*>(this), level, 0); @@ -11006,11 +11062,11 @@ Triangulation::begin_raw_quad (const unsigned int level) const 0, 0); } - - + + default: Assert (false, ExcNotImplemented()); - return raw_hex_iterator(); + return raw_hex_iterator(); } } @@ -11073,7 +11129,7 @@ Triangulation::end_quad (const unsigned int level) const quad_iterator(end_quad()) : begin_quad (level+1)); else - return quad_iterator(end_quad()); + return quad_iterator(end_quad()); } @@ -11128,7 +11184,7 @@ Triangulation::last_raw_quad (const unsigned int level) const default: Assert (false, ExcNotImplemented()); return raw_quad_iterator(); - } + } } @@ -11219,15 +11275,15 @@ Triangulation::begin_raw_hex (const unsigned int level) const if (levels[level]->cells.cells.size() == 0) return end_hex(); - + return raw_hex_iterator (const_cast*>(this), level, 0); } - + default: Assert (false, ExcNotImplemented()); - return raw_hex_iterator(); + return raw_hex_iterator(); } } @@ -11368,7 +11424,7 @@ Triangulation::last_hex () const } - + template typename Triangulation::active_hex_iterator Triangulation::last_active_hex (const unsigned int level) const @@ -11607,7 +11663,7 @@ unsigned int Triangulation::n_active_lines (const unsigned int le Assert (level < number_cache.n_lines_level.size(), ExcIndexRange (level, 0, number_cache.n_lines_level.size())); Assert (dim == 1, ExcFacesHaveNoLevel()); - + return number_cache.n_active_lines_level[level]; } @@ -11713,7 +11769,7 @@ unsigned int Triangulation::n_quads (const unsigned int level) co { Assert (dim == 2, ExcFacesHaveNoLevel()); Assert (level < number_cache.n_quads_level.size(), - ExcIndexRange (level, 0, number_cache.n_quads_level.size())); + ExcIndexRange (level, 0, number_cache.n_quads_level.size())); return number_cache.n_quads_level[level]; } @@ -11783,7 +11839,7 @@ unsigned int Triangulation::n_active_quads (const unsigned int le Assert (level < number_cache.n_quads_level.size(), ExcIndexRange (level, 0, number_cache.n_quads_level.size())); Assert (dim == 2, ExcFacesHaveNoLevel()); - + return number_cache.n_active_quads_level[level]; } @@ -11844,7 +11900,7 @@ unsigned int Triangulation::n_hexs (const unsigned int level) con { Assert (level < number_cache.n_hexes_level.size(), ExcIndexRange (level, 0, number_cache.n_hexes_level.size())); - + return number_cache.n_hexes_level[level]; } @@ -11871,7 +11927,7 @@ unsigned int Triangulation::n_active_hexs (const unsigned int lev { Assert (level < number_cache.n_hexes_level.size(), ExcIndexRange (level, 0, number_cache.n_hexes_level.size())); - + return number_cache.n_active_hexes_level[level]; } @@ -11879,53 +11935,9 @@ unsigned int Triangulation::n_active_hexs (const unsigned int lev #endif -template -unsigned int Triangulation::n_levels () const -{ - if (levels.size() == 0) - return 0; - // check whether there are cells on - // the highest levels (there need - // not be, since they might all - // have been coarsened away) - raw_cell_iterator cell = last_raw (levels.size()-1), - endc = end(); - for (; cell!=endc; --cell) - if (cell->used()) - // return level of most refined - // existing cell (+1 because of - // counting conventions) - return cell->level()+1; - - // no cells at all? - Assert (false, ExcInternalError()); - // just to make the compiler happy: - return 0; -} - - - -template -unsigned int -Triangulation::n_vertices () const -{ - return vertices.size(); -} - - - -template -const std::vector > & -Triangulation::get_vertices () const -{ - return vertices; -} - - - template unsigned int -Triangulation::n_used_vertices () const +Triangulation::n_used_vertices () const { return std::count_if (vertices_used.begin(), vertices_used.end(), std::bind2nd (std::equal_to(), true)); @@ -11935,7 +11947,7 @@ Triangulation::n_used_vertices () const template const std::vector & -Triangulation::get_used_vertices () const +Triangulation::get_used_vertices () const { return vertices_used; } @@ -12016,7 +12028,7 @@ Triangulation::execute_coarsening_and_refinement () // information update_neighbors(*this); #endif - + // Inform RefinementListeners // about end of refinement. for (ref_listener = refinement_listeners.begin (); @@ -12038,16 +12050,16 @@ Triangulation::clear_despite_subscriptions() for (unsigned int i=0; i(); } @@ -12063,10 +12075,10 @@ Triangulation::execute_refinement () // re-compute number of lines - compute_number_cache (*this, number_cache); - + compute_number_cache (*this, levels.size(), number_cache); + #ifdef DEBUG - for (unsigned int level=0; levelcells.monitor_memory (dim); // check whether really all @@ -12081,7 +12093,7 @@ Triangulation::execute_refinement () cell_iterator cell = begin(), endc = end(); while (cell != endc) - Assert (!(cell++)->refine_flag_set(), ExcInternalError ()); + Assert (!(cell++)->refine_flag_set(), ExcInternalError ()); #endif return cells_with_distorted_children; @@ -12099,7 +12111,7 @@ void Triangulation::execute_coarsening () // other dimensions it will be ignored std::vector line_cell_count = count_cells_bounded_by_line (*this); std::vector quad_cell_count = count_cells_bounded_by_quad (*this); - + // loop over all cells. Flag all // cells of which all children are // flagged for @@ -12120,7 +12132,7 @@ void Triangulation::execute_coarsening () cell_iterator cell = begin(), endc = end(); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell) if (!cell->active()) if (cell->child(0)->coarsen_flag_set()) { @@ -12160,7 +12172,7 @@ void Triangulation::execute_coarsening () // re-compute number of lines and // quads - compute_number_cache (*this, number_cache); + compute_number_cache (*this, levels.size(), number_cache); // in principle no user flags // should be @@ -12203,7 +12215,7 @@ void Triangulation::fix_coarsen_flags () save_coarsen_flags (previous_coarsen_flags); bool continue_iterating = true; - + do { if (smooth_grid & limit_level_difference_at_vertices) @@ -12212,7 +12224,7 @@ void Triangulation::fix_coarsen_flags () ExcMessage("In case of anisotropic refinement the " "limit_level_difference_at_vertices flag for " "mesh smoothing must not be set!")); - + // store highest level one // of the cells adjacent to // a vertex belongs to @@ -12258,7 +12270,7 @@ void Triangulation::fix_coarsen_flags () cell->level()-1); } } - + // loop over all cells in reverse // order. do so because we can then @@ -12294,7 +12306,7 @@ void Triangulation::fix_coarsen_flags () cell->level()+1) { cell->set_refine_flag(); - + for (unsigned int v=0; v::vertices_per_cell; ++v) vertex_level[cell->vertex_index(v)] @@ -12310,7 +12322,7 @@ void Triangulation::fix_coarsen_flags () } } } - + // loop over all cells. Flag all // cells of which all children are // flagged for coarsening and @@ -12342,35 +12354,35 @@ void Triangulation::fix_coarsen_flags () end_ac = end_active(0); for (; acell!=end_ac; ++acell) acell->clear_coarsen_flag(); - + cell_iterator cell = begin(), endc = end(); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell) { // nothing to do if we are // already on the finest level - if (cell->active()) + if (cell->active()) continue; const unsigned int n_children=cell->n_children(); unsigned int flagged_children=0; for (unsigned int child=0; childchild(child)->active() && - cell->child(child)->coarsen_flag_set()) + cell->child(child)->coarsen_flag_set()) { ++flagged_children; // clear flag since we // don't need it anymore cell->child(child)->clear_coarsen_flag(); } - + // flag this cell for // coarsening if all children // were flagged if (flagged_children == n_children) cell->set_user_flag(); } - + // in principle no coarsen flags // should be set any more at this // point @@ -12412,15 +12424,15 @@ void Triangulation::fix_coarsen_flags () { Assert (cell->child(c)->refine_flag_set()==false, ExcInternalError()); - + cell->child(c)->set_coarsen_flag(); } - + // clear all user flags again, now that we // don't need them any more clear_user_flags (); - + // now see if anything has // changed in the last // iteration of this function @@ -12448,10 +12460,10 @@ bool Triangulation<1,1>::prepare_coarsening_and_refinement () // do nothing in 1d, except setting // the coarsening flags correctly fix_coarsen_flags (); - + std::vector flags_after; save_coarsen_flags (flags_after); - + return (flags_before != flags_after); } @@ -12468,10 +12480,10 @@ bool Triangulation<1,2>::prepare_coarsening_and_refinement () // do nothing in 1d, except setting // the coarsening flags correctly fix_coarsen_flags (); - + std::vector flags_after; save_coarsen_flags (flags_after); - + return (flags_before != flags_after); } @@ -12508,7 +12520,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // break this loop at an arbitrary // number of runs, but that would // not be a clean solution, since - // we would either have to + // we would either have to // 1/ break the loop too early, in which // case the promise that a second // call to this function immediately @@ -12647,12 +12659,12 @@ bool Triangulation::prepare_coarsening_and_refinement () active_cell_iterator cell=begin_active(1), endc=end_active(1); - + for (; cell!=endc; ++cell) cell->clear_coarsen_flag(); } - + bool mesh_changed_in_this_loop = false; do { @@ -12670,7 +12682,7 @@ bool Triangulation::prepare_coarsening_and_refinement () { cell_iterator cell; const cell_iterator endc = end(); - + for (cell=begin(); cell!=endc; ++cell) { if (!cell->active()) @@ -12715,7 +12727,7 @@ bool Triangulation::prepare_coarsening_and_refinement () } // for (all cells) } // if (smooth_grid & ...) - + ////////////////////////////////////// // STEP 2: // eliminate refined islands in the @@ -12728,11 +12740,11 @@ bool Triangulation::prepare_coarsening_and_refinement () // be automatically fulfilled. if (smooth_grid & (eliminate_refined_inner_islands | eliminate_refined_boundary_islands) && - !(smooth_grid & patch_level_1)) + !(smooth_grid & patch_level_1)) { cell_iterator cell; const cell_iterator endc = end(); - + for (cell=begin(); cell!=endc; ++cell) if (!cell->active() || (cell->active() && cell->refine_flag_set())) { @@ -12753,13 +12765,13 @@ bool Triangulation::prepare_coarsening_and_refinement () bool all_children_active = true; if (!cell->active()) for (unsigned int c=0; cn_children(); ++c) - if (!cell->child(c)->active()) + if (!cell->child(c)->active()) { all_children_active = false; break; } - if (all_children_active) + if (all_children_active) { // count number // of refined and @@ -12780,7 +12792,7 @@ bool Triangulation::prepare_coarsening_and_refinement () unsigned int unrefined_neighbors = 0, total_neighbors = 0; - for (unsigned int n=0; n::faces_per_cell; ++n) + for (unsigned int n=0; n::faces_per_cell; ++n) { const cell_iterator neighbor = cell->neighbor(n); if (neighbor.state() == IteratorState::valid) @@ -12790,9 +12802,9 @@ bool Triangulation::prepare_coarsening_and_refinement () if (!face_will_be_refined_by_neighbor(cell,n)) ++unrefined_neighbors; } - + } - + // if all // neighbors // unrefined: @@ -12840,7 +12852,7 @@ bool Triangulation::prepare_coarsening_and_refinement () cell->child(c)->clear_refine_flag (); cell->child(c)->set_coarsen_flag (); } - else + else cell->clear_refine_flag(); } } @@ -12864,7 +12876,7 @@ bool Triangulation::prepare_coarsening_and_refinement () ExcMessage("In case of anisotropic refinement the " "limit_level_difference_at_vertices flag for " "mesh smoothing must not be set!")); - + // store highest level one // of the cells adjacent to // a vertex belongs to @@ -12902,7 +12914,7 @@ bool Triangulation::prepare_coarsening_and_refinement () cell->level()-1); } } - + // loop over all cells in reverse // order. do so because we can then @@ -12938,7 +12950,7 @@ bool Triangulation::prepare_coarsening_and_refinement () cell->level()+1) { cell->set_refine_flag(); - + for (unsigned int v=0; v::vertices_per_cell; ++v) vertex_level[cell->vertex_index(v)] @@ -12956,7 +12968,7 @@ bool Triangulation::prepare_coarsening_and_refinement () } ///////////////////////////////////// - // STEP 4: + // STEP 4: // eliminate unrefined // islands. this has higher // priority since this @@ -12969,7 +12981,7 @@ bool Triangulation::prepare_coarsening_and_refinement () { active_cell_iterator cell = begin_active(), endc = end(); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell) { // if cell is already // flagged for (isotropic) @@ -12995,7 +13007,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // current cell is flagged to be // refined in an according // direction. - + if (!(smooth_grid & allow_anisotropic_smoothing)) { // use first algorithm @@ -13004,7 +13016,7 @@ bool Triangulation::prepare_coarsening_and_refinement () for (unsigned int face=0; face::faces_per_cell; ++face) if (!cell->at_boundary(face)) { - if (face_will_be_refined_by_neighbor(cell,face)) + if (face_will_be_refined_by_neighbor(cell,face)) ++refined_neighbors; else ++unrefined_neighbors; @@ -13018,13 +13030,13 @@ bool Triangulation::prepare_coarsening_and_refinement () } } else - { + { // variable to store the cell // refine case needed to // fulfill all smoothing // requirements RefinementCase smoothing_cell_refinement_case=RefinementCase::no_refinement; - + // use second algorithm, do // the check individually for // each direction @@ -13037,7 +13049,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // in the same way as the // neighbors do... RefinementCase directional_cell_refinement_case=RefinementCase::isotropic_refinement; - + for (unsigned int face_index=0; face_index<2; ++face_index) { unsigned int face=2*face_pair+face_index; @@ -13047,7 +13059,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // face under // consideration RefinementCase expected_face_ref_case=RefinementCase::no_refinement; - + if (cell->neighbor(face).state() == IteratorState::valid) face_will_be_refined_by_neighbor(cell,face,expected_face_ref_case); // now extract which @@ -13120,7 +13132,7 @@ bool Triangulation::prepare_coarsening_and_refinement () cell->clear_coarsen_flag(); cell->set_refine_flag(cell->refine_flag_set() | smoothing_cell_refinement_case); } - }//else -> allow_anisotropic_smoothing + }//else -> allow_anisotropic_smoothing }// for all cells } @@ -13155,7 +13167,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // sets refinement flags and // does not set coarsening // flags. - if (smooth_grid & patch_level_1) + if (smooth_grid & patch_level_1) { // An important assumption // (A) is that before @@ -13190,7 +13202,7 @@ bool Triangulation::prepare_coarsening_and_refinement () ExcInternalError()); if (cell->child(0)->has_children() == true) continue; - + // cell is found to // be a patch. // combine the refine @@ -13204,7 +13216,7 @@ bool Triangulation::prepare_coarsening_and_refinement () for (unsigned int i=0; in_children(); ++i) { cell_iterator child = cell->child(i); - + child->clear_coarsen_flag(); child->set_refine_flag(combined_ref_case); } @@ -13272,7 +13284,7 @@ bool Triangulation::prepare_coarsening_and_refinement () if (has_active_grandchildren == false) continue; - + // ok, there are active // grandchildren. see @@ -13280,7 +13292,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // none of them are // flagged for // coarsening - unsigned int n_grandchildren=0; + unsigned int n_grandchildren=0; // count all coarsen // flags of the // grandchildren. @@ -13405,9 +13417,9 @@ bool Triangulation::prepare_coarsening_and_refinement () changed=false; active_cell_iterator cell=last_active(), endc=end(); - + for (; cell != endc; --cell) - if (cell->refine_flag_set()) + if (cell->refine_flag_set()) { // loop over neighbors of cell for (unsigned int i=0; i::faces_per_cell; ++i) @@ -13497,7 +13509,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // coarser at the given face. Now the upper child element of the lower left // wants to refine according to cut_z, such that there is a 'horizontal' // refinement of the face marked with ##### -// +// // / / // / / // *---------------* @@ -13509,8 +13521,8 @@ bool Triangulation::prepare_coarsening_and_refinement () // | | / // | |/ // *---------------* -// -// +// +// // *---------------* // /| /| // / | ##### / | @@ -13526,7 +13538,7 @@ bool Triangulation::prepare_coarsening_and_refinement () // this introduces too many hanging nodes and the neighboring (coarser) cell // (upper right) has to be refined. If it is only refined according to cut_z, // then everything is ok: -// +// // / / // / / // *---------------* @@ -13538,8 +13550,8 @@ bool Triangulation::prepare_coarsening_and_refinement () // | | / // | |/ // *---------------* -// -// +// +// // *---------------* // /| /| // / *---------------* @@ -13568,8 +13580,8 @@ bool Triangulation::prepare_coarsening_and_refinement () // | | | / // | | |/ // *-------*-------* -// -// +// +// // *---------------* // /| /| // / *---------------* @@ -13596,8 +13608,8 @@ bool Triangulation::prepare_coarsening_and_refinement () // | | | / // | | |/ // *-------*-------* -// -// +// +// // *-------*-------* // /| /| /| // / *-------*-------* @@ -13618,7 +13630,7 @@ bool Triangulation::prepare_coarsening_and_refinement () to_be_refined_along_y=0; const int this_face_index=cell->face_index(i); - + // step 1: detect, along which axis the face is currently refined if ((this_face_index == cell->neighbor(i)->face(nb_indices.first)->child_index(0)) || @@ -13670,8 +13682,8 @@ bool Triangulation::prepare_coarsening_and_refinement () ++to_be_refined_along_x; if (flagged_frc & RefinementCase::cut_y) ++to_be_refined_along_y; - -// step 3: set the refine flag of the (coarser and active) neighbor. + +// step 3: set the refine flag of the (coarser and active) neighbor. if ((smooth_grid & allow_anisotropic_smoothing) || cell->neighbor(i)->refine_flag_set()) { @@ -13688,7 +13700,7 @@ bool Triangulation::prepare_coarsening_and_refinement () changed=true; cell->neighbor(i)->set_refine_flag(); } - + // step 4: if necessary (see above) add to the refine flag of the current cell cell_iterator nb=cell->neighbor(i); RefinementCase nb_frc @@ -13764,14 +13776,14 @@ bool Triangulation::prepare_coarsening_and_refinement () } } } - + ////////////////////////////////////// // STEP 7: // take care that no double refinement // is done at each line in 3d or higher // dimensions. internal::Triangulation::Implementation::prepare_refinement_dim_dependent (*this); - + ////////////////////////////////////// // STEP 8: // make sure that all children of each @@ -13797,7 +13809,7 @@ bool Triangulation::prepare_coarsening_and_refinement () } while (mesh_changed_in_this_loop); - + // find out whether something was really // changed in this function. Note that // @p{flags_before_loop} represents the @@ -13819,23 +13831,23 @@ void Triangulation::write_bool_vector (const unsigned int magic_ const unsigned int N = v.size(); unsigned char *flags = new unsigned char[N/8+1]; for (unsigned int i=0; i(flags[i]) << ' '; - + out << std::endl << magic_number2 << std::endl; - + delete[] flags; AssertThrow (out, ExcIO()); @@ -13860,7 +13872,7 @@ void Triangulation::read_bool_vector (const unsigned int magic_n unsigned char *flags = new unsigned char[N/8+1]; unsigned short int tmp; - for (unsigned int i=0; i> tmp; flags[i] = tmp; @@ -13881,7 +13893,7 @@ void Triangulation::read_bool_vector (const unsigned int magic_n template unsigned int -Triangulation::memory_consumption () const +Triangulation::memory_consumption () const { unsigned int mem = 0; mem += MemoryConsumption::memory_consumption(levels); @@ -13893,7 +13905,7 @@ Triangulation::memory_consumption () const mem += sizeof(smooth_grid); mem += MemoryConsumption::memory_consumption (number_cache); mem += sizeof (faces); - mem += MemoryConsumption::memory_consumption (*faces); + mem += MemoryConsumption::memory_consumption (*faces); return mem; } @@ -13980,7 +13992,7 @@ template const Boundary<1,3> & Triangulation<1,3>::get_boundary (const unsigned int number) const; -#else +#else # if deal_II_dimension == 2 template class Triangulation<2,3>; # endif diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c2df68f6ae..8cb12751a2 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -359,9 +359,9 @@ inconvenience this causes.
  • - Improved: The Vector class has been equipped with an improved way to + Improved: The Vector class has been equipped with an improved way to calculate sums in inner products and norms. This reduces the accumulation - of round-off errors. Especially the solution with float vectors should + of round-off errors. Especially the solution with float vectors should profit from the new implementation.


    @@ -511,6 +511,19 @@ inconvenience this causes.

    deal.II

      +
    1. +

      + Fixed: The function Triangulation::n_levels() accidentally turned out to be + quite expensive, in particular if the mesh has been coarsened significantly + in the past. Since Triangulation::n_active_cells() calls Triangulation::n_levels() + repeatedly, by consequence it also is surprisingly expensive, which is + particularly annoying since this is a frequently called function. This has + now been fixed by caching the result of this function. +
      + (Wolfgang Bangerth 2009/11/22) +

      +
    2. +
    3. Fixed: The function DataOut class got confused when its DataOut::first_cell() and -- 2.39.5