From 47ec3bffe224f679edf39805f40c395afe729caa Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 17 Jun 1998 09:47:50 +0000 Subject: [PATCH] Fix some really nasty bugs which came to the surface with the implementation of coarsening. Some of them should have wasted large amounts of memory :-( git-svn-id: https://svn.dealii.org/trunk@406 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/dofs/dof_handler.cc | 58 +++++---- deal.II/deal.II/source/grid/tria.cc | 133 ++++++++++++--------- 2 files changed, 106 insertions(+), 85 deletions(-) diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index bcd59f9096..bc9b73ee0e 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -779,6 +779,7 @@ void DoFHandler::distribute_dofs (const FiniteElementBase &fe) { unsigned int next_free_dof = 0; active_cell_iterator cell = begin_active(), endc = end(); + for (; cell != endc; ++cell) next_free_dof = distribute_dofs_on_cell (cell, next_free_dof); @@ -857,23 +858,19 @@ unsigned int DoFHandler<2>::distribute_dofs_on_cell (active_cell_iterator &cell, cell->set_vertex_dof_index (vertex, d, next_free_dof++); // for the four sides - for (unsigned int side=0; side::faces_per_cell; ++side) - { - line_iterator line = cell->line(side); - - // distribute dofs if necessary - if ((selected_fe->dofs_per_line > 0) && + if (selected_fe->dofs_per_line > 0) + for (unsigned int side=0; side::faces_per_cell; ++side) + { + line_iterator line = cell->line(side); + + // distribute dofs if necessary: // check whether line dof is already // numbered (check only first dof) - (line->dof_index(0) == -1)) - // if not: distribute dofs - for (unsigned int d=0; ddofs_per_line; ++d) - line->set_dof_index (d, next_free_dof++); - - // note if line is subject to constraints - if (line->child_index(0) != -1) - line->set_user_flag (); - }; + if (line->dof_index(0) == -1) + // if not: distribute dofs + for (unsigned int d=0; ddofs_per_line; ++d) + line->set_dof_index (d, next_free_dof++); + }; // dofs of quad @@ -1156,6 +1153,7 @@ void DoFHandler<2>::make_constraint_matrix (ConstraintMatrix &constraints) const for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->has_children()) cell->face(face)->set_user_flag(); + @@ -1834,16 +1832,17 @@ void DoFHandler<1>::reserve_space () { delete levels[i]; levels.resize (0); - vertex_dofs.erase (vertex_dofs.begin(), vertex_dofs.end()); - vertex_dofs.resize (tria->vertices.size(), -1); + vertex_dofs = vector(tria->vertices.size()* + selected_fe->dofs_per_vertex, + -1); for (unsigned int i=0; in_levels(); ++i) { levels.push_back (new DoFLevel<1>); - levels.back()->line_dofs.resize (tria->levels[i]->lines.lines.size() * - selected_fe->dofs_per_line, - -1); + levels.back()->line_dofs = vector(tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); }; }; @@ -1865,20 +1864,19 @@ void DoFHandler<2>::reserve_space () { delete levels[i]; levels.resize (0); - vertex_dofs.erase (vertex_dofs.begin(), vertex_dofs.end()); - vertex_dofs.resize (tria->vertices.size(), - -1); - + vertex_dofs = vector(tria->vertices.size()* + selected_fe->dofs_per_vertex, + -1); for (unsigned int i=0; in_levels(); ++i) { levels.push_back (new DoFLevel<2>); - levels.back()->line_dofs.resize (tria->levels[i]->lines.lines.size() * - selected_fe->dofs_per_line, - -1); - levels.back()->quad_dofs.resize (tria->levels[i]->quads.quads.size() * - selected_fe->dofs_per_quad, - -1); + levels.back()->line_dofs = vector (tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); + levels.back()->quad_dofs = vector (tria->levels[i]->quads.quads.size() * + selected_fe->dofs_per_quad, + -1); }; }; diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index e3a6d5f203..7ab1b8987f 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -605,7 +605,7 @@ void Triangulation<2>::create_hyper_L (const double a, const double b) { template <> void Triangulation<2>::create_hyper_ball (const Point<2> &p, const double radius) { const unsigned int dim=2; - const double a = 0.453; // equilibrate cell sizes + const double a = 1./(1+sqrt(2)/2); // equilibrate cell sizes const Point vertices[8] = { p+Point(-1,-1)*(radius/sqrt(2)), p+Point(+1,-1)*(radius/sqrt(2)), p+Point(-1,-1)*(radius/sqrt(2)*a), @@ -1815,7 +1815,26 @@ unsigned int Triangulation::n_active_quads (const unsigned int level) const template unsigned int Triangulation::n_levels () const { - return levels.size(); + 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; }; @@ -1922,6 +1941,14 @@ void Triangulation::refine (const dVector &criteria, Assert (criteria.size() == n_active_cells(), ExcInvalidVectorSize(criteria.size(), n_active_cells())); + // nothing to do; especially we + // do not want to flag with zero + // error since then we may get + // into conflict with coarsening + // in some cases + if (threshold==0) + return; + active_cell_iterator cell = begin_active(); const unsigned int n_cells = criteria.size(); @@ -2026,8 +2053,6 @@ Triangulation::refine_and_coarsen_fixed_fraction (const dVector &criteria, else bottom_threshold = sqrt(*p - *(p-1)); - Assert (bottom_threshold<=top_threshold, ExcInternalError()); - // in some rare cases it may happen that // both thresholds are the same (e.g. if // there are many cells with the same @@ -2038,10 +2063,16 @@ Triangulation::refine_and_coarsen_fixed_fraction (const dVector &criteria, // which the indicator equals the // thresholds. This is forbidden, however. // - // In that case we arbitrarily reduce the - // bottom threshold by one permille. - if (bottom_threshold==top_threshold) - bottom_threshold *= 0.999; + // In some rare cases with very few cells + // we also could get integer round off + // errors and get problems with + // the top and bottom fractions. + // + // In these case we arbitrarily reduce the + // bottom threshold by one permille below + // the top threshold + if (bottom_threshold>=top_threshold) + bottom_threshold = 0.999*top_threshold; // actually flag cells refine (criteria, top_threshold); @@ -2111,11 +2142,10 @@ void Triangulation<1>::execute_refinement () { // count number of used cells on // the next higher level - unsigned int used_cells = 0; - count_if (levels[level+1]->lines.used.begin(), - levels[level+1]->lines.used.end(), - bind2nd (equal_to(), true), - used_cells); + const unsigned int used_cells + = count_if (levels[level+1]->lines.used.begin(), + levels[level+1]->lines.used.end(), + bind2nd (equal_to(), true)); // reserve space for the used_cells // cells already existing on the next @@ -2137,12 +2167,18 @@ void Triangulation<1>::execute_refinement () { // add to needed vertices how // many vertices are already in use - count_if (vertices_used.begin(), vertices_used.end(), - bind2nd (equal_to(), true), - needed_vertices); - - vertices.resize (needed_vertices, Point<1>()); - vertices_used.resize (needed_vertices, false); + needed_vertices += count_if (vertices_used.begin(), vertices_used.end(), + bind2nd (equal_to(), true)); + // if we need more vertices: create them, + // if not: leave the array as is, since + // shrinking is not really possible because + // some of the vertices at the end may be + // in use + if (needed_vertices > vertices.size()) + { + vertices.resize (needed_vertices, Point<1>()); + vertices_used.resize (needed_vertices, false); + }; // Do REFINEMENT @@ -2401,21 +2437,12 @@ void Triangulation<2>::execute_refinement () { }; }; - // count number of used lines on - // the next higher level - unsigned int used_lines = 0; - count_if (levels[level+1]->lines.used.begin(), - levels[level+1]->lines.used.end(), - bind2nd (equal_to(), true), - used_lines); - // count number of used cells on // the next higher level - unsigned int used_cells = 0; - count_if (levels[level+1]->quads.used.begin(), + const unsigned int used_cells + = count_if (levels[level+1]->quads.used.begin(), levels[level+1]->quads.used.end(), - bind2nd (equal_to(), true), - used_cells); + bind2nd (equal_to(), true)); // reserve space for the used_cells @@ -2438,14 +2465,18 @@ void Triangulation<2>::execute_refinement () { // add to needed vertices how // many vertices are already in use - count_if (vertices_used.begin(), vertices_used.end(), - bind2nd (equal_to(), true), - needed_vertices); - - // reserve enough space for all vertices - vertices.resize (needed_vertices, Point<2>()); - vertices_used.resize (needed_vertices, false); - + needed_vertices += count_if (vertices_used.begin(), vertices_used.end(), + bind2nd (equal_to(), true)); + // if we need more vertices: create them, + // if not: leave the array as is, since + // shrinking is not really possible because + // some of the vertices at the end may be + // in use + if (needed_vertices > vertices.size()) + { + vertices.resize (needed_vertices, Point()); + vertices_used.resize (needed_vertices, false); + }; // Do REFINEMENT @@ -3537,14 +3568,10 @@ void TriangulationLevel<0>::monitor_memory (const unsigned int true_dimension) c void TriangulationLevel<1>::reserve_space (const unsigned int new_lines) { - vector::iterator u = lines.used.begin(), - e = lines.used.end(); - - unsigned int used_lines = 0; - for (; u!=e; ++u) - ++used_lines; - - const unsigned int new_size = used_lines + new_lines; + const unsigned int new_size = new_lines + + count_if (lines.used.begin(), + lines.used.end(), + bind2nd (equal_to(), true)); // same as in #reserve_space<0>#: only // allocate space if necessary @@ -3616,14 +3643,10 @@ void TriangulationLevel<1>::monitor_memory (const unsigned int true_dimension) c #if deal_II_dimension >= 2 void TriangulationLevel<2>::reserve_space (const unsigned int new_quads) { - vector::iterator u = quads.used.begin(), - e = quads.used.end(); - - unsigned int used_quads = 0; - for (; u!=e; ++u) - ++used_quads; - - const unsigned int new_size = used_quads + new_quads; + const unsigned int new_size = new_quads + + count_if (quads.used.begin(), + quads.used.end(), + bind2nd (equal_to(), true)); // see above... if (new_size>quads.quads.size()) -- 2.39.5