From: Wolfgang Bangerth Date: Fri, 12 Jun 1998 08:42:37 +0000 (+0000) Subject: Changes related to the implementation of coarsening. X-Git-Tag: v8.0.0~22858 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3fe25ad7e9abcb0aaa96dacb8a1b9a91be3207d7;p=dealii.git Changes related to the implementation of coarsening. git-svn-id: https://svn.dealii.org/trunk@396 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index c2bf205983..803ce6d472 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -263,7 +263,7 @@ class ConstraintMatrix { */ DeclException1 (ExcLineExists, unsigned int, - << "The lines " << arg1 + << "The line " << arg1 << " which is to be created already exists."); /** * Exception diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 4afb2ad8ac..137c8f56cd 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -614,6 +614,10 @@ class DoFHandler : public DoFDimensionInfo { * * To condense a given sparsity pattern, * use #ConstraintMatrix::condense#. + * + * This function uses the user flags for + * the faces. If you need the user flags, + * store them beforehand. */ void make_constraint_matrix (ConstraintMatrix &) const; diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 7d1a1c1e44..253344676f 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -860,7 +860,7 @@ enum MeshSmoothing { * interface between regions of different materials. * * - * \subsection{Refinement and of a triangulation} + * \subsection{Refinement and coarsening of a triangulation} * * Refinement of a triangulation may be done through several ways. The most * low-level way is directly through iterators: let #i# be an iterator to @@ -881,6 +881,17 @@ enum MeshSmoothing { * structures and algorithms much much easier. To be honest, this is mostly * an algorithmic step than one needed by the finite element method. * + * To coarsen a grid, the same way as above is possible by using + * #i->set_coarsen_flag# and calling #execute_coarsening#. You can use + * #execute_coarsening_and_refinement# to get both actions done, first + * coarsening and refinement. The reason for this order is that the + * refinement usually adds some additional cells to keep the triangulation + * regular and thus satifies all refinement requests, while the coarsening + * does not delete cells not requested for; therefore the refinement will + * often revert some effects of coarsening while the opposite is not true. + * The stated order of coarsening before refinement will thus normally + * lead to a result closer to the intended one. + * * Marking cells for refinement 'by hand' through iterators is one way to * produce a new grid, especially if you know what kind of grid you are * looking for, e.g. if you want to have a grid successively refined @@ -895,20 +906,22 @@ enum MeshSmoothing { * vector of values, one per active cell, which denote the criterion according * to which the triangulation is to be refined. It marks all cells for which * the criterion is greater than the threshold being given as the second - * argument. + * argument. Analogously, + * #coarsen (const dVector &criterion, const double threshold)# flags those + * cells for coarsening for which the criterion is less than the treshold. * - * There are two variations of this function, which rely on #refine# by - * computing the threshold from other information: + * There are two variations of these functions, which rely on #refine# and + * coarsen by computing the thresholds from other information: * \begin{itemize} - * \item #refine_fixed_number#: this function takes a vector as above and - * a value between zero and one denoting the fraction of cells to be - * refined. For this purpose, it sorts the criteria per cell and takes - * the threshold to be the one belonging to the cell with the + * \item #refine_and_coarsen_fixed_number#: this function takes a vector as + * above and two values between zero and one denoting the fractions of cells to + * be refined and coarsened. For this purpose, it sorts the criteria per cell + * and takes the threshold to be the one belonging to the cell with the * #fraction times n_active_cells# highest criterion. For example, if * the fraction is $0.3$, the threshold is computed to a value such that * 30 per cent of cells have a criterion higher than the threshold and are * thus flagged for refinement. The flagging for refinement is done through - * the central #refine# function. + * the central #refine# function. For coarsening, the same holds. * * The sorting of criteria is not done actually, since we only need one * value, in the example above the criterion of the cell which is at @@ -919,15 +932,15 @@ enum MeshSmoothing { * than #N log N# for sorting all values. * * A typical value for the fraction of cells to be refined is 0.3. - * However, for singular functions or error functionals, you may want to - * chose a smaller value to avoid overrefinement in regions which do not - * contribute much to the error. - * - * \item #refine_fixed_fraction#: this function computes the threshold such - * that the number of cells getting flagged for refinement makes up for a - * certain fraction of the total error. If this fraction is 50 per cent, - * for example, the threshold is computed such that the cells with a - * criterion greater than the threshold together account for half of the + * However, for singular functions or singular error functionals, you may + * want to chose a smaller value to avoid overrefinement in regions which + * do not contribute much to the error. + * + * \item #refine_and_coarsen_fixed_fraction#: this function computes the + * threshold such that the number of cells getting flagged for refinement + * makes up for a certain fraction of the total error. If this fraction is 50 + * per cent, for example, the threshold is computed such that the cells with + * a criterion greater than the threshold together account for half of the * total error. The definition of the fraction is a bit unintuitive, since * the total error is the sum over all cells of the local contribution * squared. We define that the fraction $\alpha$ be such that those @@ -938,37 +951,19 @@ enum MeshSmoothing { * indicator with $\eta^2 = \sum \eta_K^2$, with here the sum running over * all cells. * + * For the bottom fraction the same holds: the treshold for coarsening is + * computed such that the cells with criterion less than the threshold + * together make up for the fraction of the total error specified. + * * This strategy is more suited for singular functions and error * functionals, but may lead to very slow convergence of the grid * if only few cells are refined in each step. * * From the implementational point, this time we really need to - * sort the array of criteria. However, it is not necessary to sort - * the whole array, since for example if you chose the fraction at - * 50 per cent of the total error, it is only necessary to sort at - * most the 50 per cent of cells ranking topmost in the list of error - * per cell. It is thus reasonable to use an algorithm like - * #partial_sort# of the C++ standard library, which only sorts part - * of the array and lets the rest unsorted. However, in many cases - * much fewer than 50 per cent of the cells account for 50 per cent - * of the error, so it may be possible to get away with sorting less - * than 50 per cent of the cells. We therefore divide the whole lot - * of 50 per cent of cells into, say, 5 parts, first sort for the - * 10 per cent with highest error; look whether they together make up - * for 50 per cent and if so thats ok, we can leave the rest unsorted; - * if not, sort the next 10 per cent, and so on. The default is to - * devide the maximum number of cells which may get refined (which - * equals the fraction of the total error, as explained above) into - * five parts, but this value may be given as a parameter to the - * #refine_fixed_fraction# function. For highly singular error - * functionals, it may be more efficient to chose a greater number - * than five. Chosing a value which is too large should not lead to - * a large performance drawback; chosing too small a value however - * may lead to significantly higher computational costs for sorting - * than necessary. - * + * sort the array of criteria. * Just like the other strategy described above, this function only - * computes the threshold value and then passes over to #refine#. + * computes the threshold values and then passes over to #refine# and + * #coarsen#. * * A typical value for the fraction of the total error is 0.5. * \end{itemize} @@ -982,6 +977,9 @@ enum MeshSmoothing { * element, such that the square of the total error is the sum over the * squares of the criteria on the cells. The criteria shall be positive. * + * You can suppress coarsening or refining by giving zero as the fraction + * for one of the operations. + * * * \subsection{Smoothing of a triangulation} * @@ -1179,6 +1177,8 @@ enum MeshSmoothing { * }; * \end{verbatim} * + * The same scheme is employed for coarsening and the coarsening flags. + * * You may write other information to the output file between different sets * of refinement information, as long as you read it upon re-creation of the * grid. You should make sure that the other information in the new @@ -1585,10 +1585,21 @@ class Triangulation : public TriaDimensionInfo { void refine (const dVector &criteria, const double threshold); + /** + * Analogue to the #refine# function: + * flag all cells for coarsening for + * which the criterion is less than the + * given threshold. + */ + void coarsen (const dVector &criteria, + const double threshold); + /** * Refine the triangulation by refining - * a certain fraction #fraction_of_cells# - * with the highest error. To actually + * a certain fraction #top_fraction_of_cells# + * with the highest error. Likewise coarsen + * the fraction #bottom_fraction_of_cells# + * with the least error. To actually * perform the refinement, call * #execute_refinement#. * @@ -1598,27 +1609,28 @@ class Triangulation : public TriaDimensionInfo { * Refer to the general doc of this class * for more information. */ - void refine_fixed_number (const dVector &criteria, - const double fraction_of_cells); + void refine_and_coarsen_fixed_number (const dVector &criteria, + const double top_fraction_of_cells, + const double bottom_fraction_of_cells); /** * Refine the triangulation by flagging * those cells which make up a certain - * #fraction_of_error# of the total error. + * #top_fraction# of the total error. + * Likewise, coarsen all cells which + * make up only #bottom_fraction#. * To actually perform the refinement, call - * #execute_refinement#. + * #execute_coarsening_and_refinement#. * - * #fraction_of_error# shall be a value + * #*_fraction# shall be a values * between zero and one. - * #n_sorting_parts# shall be one or - * greater. * * Refer to the general doc of this class * for more information. */ - void refine_fixed_fraction (const dVector &criteria, - const double fraction_of_error, - const unsigned int n_sorting_parts = 5); + void refine_and_coarsen_fixed_fraction (const dVector &criteria, + const double top_fraction, + const double bottom_fraction); /** * Refine all cells on all levels which @@ -1656,7 +1668,7 @@ class Triangulation : public TriaDimensionInfo { * Execute both refinement and coarsening * of the triangulation. */ - void execute_refinement_and_coarsening (); + void execute_coarsening_and_refinement (); /*@}*/ /** @@ -1689,6 +1701,17 @@ class Triangulation : public TriaDimensionInfo { * #save_refine_flags#. */ void load_refine_flags (istream &in); + + /** + * Analogue to #save_refine_flags#. + */ + void save_coarsen_flags (ostream &out) const; + + /** + * Analogue to #load_refine_flags#. + */ + void load_coarsen_flags (istream &out); + /*@}*/ @@ -2432,6 +2455,14 @@ class Triangulation : public TriaDimensionInfo { */ void prepare_coarsening (); + /** + * Actually delete a cell, which is the + * main step for the coarsening process. + * This is the dimension dependent part + * of #execute_coarsening#. + */ + void delete_cell (cell_iterator &cell); + /** * Array of pointers pointing to the * #TriangulationLevel# objects diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index 7d60ac272c..bcd59f9096 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -1138,7 +1138,26 @@ void DoFHandler<1>::make_constraint_matrix (ConstraintMatrix &cm) const { template <> void DoFHandler<2>::make_constraint_matrix (ConstraintMatrix &constraints) const { + const unsigned int dim = 2; + constraints.clear (); + + // first mark all faces which are subject + // to constraints. We do so by looping + // over all active cells and checking + // whether any of the faces are refined + // which can only be from the neighboring + // cell because this one is active. In that + // case, the face is subject to constraints + tria->clear_user_flags (); + Triangulation::active_cell_iterator cell = tria->begin_active(), + endc = tria->end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children()) + cell->face(face)->set_user_flag(); + + line_iterator line = begin_line(), endl = end_line(); diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 672055579e..4e5eb697a4 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -729,6 +729,74 @@ void Triangulation::load_refine_flags (istream &in) { }; + +template +void Triangulation::save_coarsen_flags (ostream &out) const { + unsigned int N = n_active_cells(); + active_cell_iterator cell = begin_active(), + endc = end(); + + unsigned char *flags = new unsigned char[N/8+1]; + for (unsigned int i=0; icoarsen_flag_set() ? (1<<(position%8)) : 0); + + // format: + // 0. magic number + // 1. number of active cells + // 2. the flags + // 3. magic number 0xabcd + out << mn_tria_coarsen_flags_begin << " " << N << endl; + for (unsigned int i=0; i(flags[i]) << " "; + + out << endl; + out << mn_tria_coarsen_flags_end << endl; + + delete[] flags; +}; + + + +template +void Triangulation::load_coarsen_flags (istream &in) { + unsigned int magic_number; + in >> magic_number; + Assert (magic_number==mn_tria_coarsen_flags_begin, ExcGridReadError()); + + unsigned int N; + in >> N; + Assert (N==n_active_cells(), ExcGridsDoNotMatch(N, n_active_cells())); + + unsigned char *flags = new unsigned char[N/8+1]; + unsigned short int tmp; + for (unsigned int i=0; i> tmp; + flags[i] = tmp; + }; + + + active_cell_iterator cell = begin_active(), + endc = end(); + unsigned int position=0; + for (; cell!=endc; ++cell, ++position) + if (flags[position/8] & (1<<(position%8))) + cell->set_coarsen_flag(); + else + cell->clear_coarsen_flag(); + + Assert (position==N, ExcGridReadError()); + + in >> magic_number; + Assert (magic_number==mn_tria_coarsen_flags_end, ExcGridReadError()); + + delete[] flags; +}; + + + #if deal_II_dimension == 1 template <> @@ -1910,23 +1978,47 @@ void Triangulation::refine (const dVector &criteria, template -void Triangulation::refine_fixed_number (const dVector &criteria, - const double fraction) { +void Triangulation::coarsen (const dVector &criteria, + const double threshold) { + Assert (criteria.size() == n_active_cells(), + ExcInvalidVectorSize(criteria.size(), n_active_cells())); + + active_cell_iterator cell = begin_active(); + const unsigned int n_cells = criteria.size(); + + for (unsigned int index=0; indexset_coarsen_flag(); +}; + + + +template +void Triangulation::refine_and_coarsen_fixed_number (const dVector &criteria, + const double top_fraction, + const double bottom_fraction) { // correct number of cells is // checked in #refine# - Assert ((fraction>0) && (fraction<=1), ExcInvalidParameterValue()); - - // refine at least one cell - const int refine_cells = max(static_cast(fraction*criteria.size()), + Assert ((top_fraction>0) && (top_fraction<=1), ExcInvalidParameterValue()); + Assert ((bottom_fraction>0) && (bottom_fraction<=1), ExcInvalidParameterValue()); + Assert (top_fraction+bottom_fraction <= 1, ExcInvalidParameterValue()); + // refine at least one cell; algorithmic + // simplification + const int refine_cells = max(static_cast(top_fraction*criteria.size()), 1); + const int coarsen_cells = max(static_cast(bottom_fraction*criteria.size()), + 1); dVector tmp(criteria); - nth_element (tmp.begin(), - tmp.begin()+refine_cells, + nth_element (tmp.begin(), tmp.begin()+refine_cells, tmp.end(), greater()); - refine (criteria, *(tmp.begin() + refine_cells)); + + nth_element (tmp.begin(), tmp.begin()+tmp.size()-coarsen_cells, + tmp.end(), + greater()); + coarsen (criteria, *(tmp.begin() + tmp.size() - coarsen_cells)); }; @@ -1940,35 +2032,16 @@ double sqr(double a) { template -void Triangulation::refine_fixed_fraction (const dVector &criteria, - const double fraction_of_error, - const unsigned int n_sorting_steps) { +void +Triangulation::refine_and_coarsen_fixed_fraction (const dVector &criteria, + const double top_fraction, + const double bottom_fraction) { // correct number of cells is // checked in #refine# - Assert ((fraction_of_error>0) && (fraction_of_error<=1), - ExcInvalidParameterValue()); - - // rename variable since we have to change it - unsigned n_sorting_parts = n_sorting_steps; - + Assert ((top_fraction>0) && (top_fraction<=1), ExcInvalidParameterValue()); + Assert ((bottom_fraction>0) && (bottom_fraction<=1), ExcInvalidParameterValue()); + Assert (top_fraction+bottom_fraction <= 1, ExcInvalidParameterValue()); - // number of cells to be sorted per part - unsigned cells_per_part - = static_cast(rint(fraction_of_error * criteria.size() / n_sorting_parts)); - - // if number of elements is so small or the - // fraction so high that we will get into trouble - // with the maximum number of elements to be - // sorted, fall back to only one sorting step. - // Do so also if cells_per_part was rounded - // to zero - if ((cells_per_part*n_sorting_parts > criteria.size()) || - (cells_per_part == 0)) - { - cells_per_part = criteria.size(); - n_sorting_parts = 1; - }; - // let tmp be the cellwise square of the // error, which is what we have to sum // up and compare with @@ -1978,70 +2051,46 @@ void Triangulation::refine_fixed_fraction (const dVector &criteria, const double total_error = tmp.l1_norm(); dVector partial_sums(criteria.size()); - for (unsigned int part=0; part()); - // compute partial sum of the range - // as yet sorted. In principle it - // would be sufficient to only sum up - // the newly sorted part and give the - // partial sum an initial value equal - // to the previously last partial sum, - // but at present I do not know how - // to do so in an easy way. Think - // about it and fix it if you want! - // (This way doesn't eat up much - // computing time anyway, much less - // than the sorting, so I don't care - // about fixing this myself.) - partial_sum (tmp.begin(), - tmp.begin()+(part+1)*cells_per_part, - partial_sums.begin()); - - // check whether the sorted - // region already is enough - if (*(partial_sums.begin()+(part+1)*cells_per_part-1) >= - (fraction_of_error*total_error)) - { - // find first entry in the partial - // sum which is greater than the - // fraction of the error. We only - // need to search the newly created - //region - const dVector::const_iterator threshold_ptr - = lower_bound (partial_sums.begin()+part*cells_per_part, - partial_sums.begin()+(part+1)*cells_per_part, - fraction_of_error*total_error); - Assert (threshold_ptr()); + partial_sum (tmp.begin(), tmp.end(), partial_sums.begin()); + + // compute thresholds + dVector::const_iterator p; + double top_threshold, bottom_threshold; + p = lower_bound (partial_sums.begin(), partial_sums.end(), + top_fraction*total_error); + if (p==partial_sums.begin()) + top_threshold = sqrt(*p); + else + top_threshold = sqrt(*p - *(p-1)); + + p = upper_bound (partial_sums.begin(), partial_sums.end(), + total_error*(1-bottom_fraction)); + if (p==partial_sums.end()) + bottom_threshold = 0; + 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 + // error indicator). That would mean that + // all cells will be flagged for + // refinement or coarsening, but some will + // be flagged for both, namely those for + // which the indicator equals the + // thresholds. This is forbidden, however. // - // Only exception: there are so few cells - // that fraction*n_cells == n_cells - // (integer arithmetic!) - Assert (n_sorting_parts * cells_per_part == criteria.size(), - ExcInternalError()); + // In that case we arbitrarily reduce the + // bottom threshold by one permille. + if (bottom_threshold==top_threshold) + bottom_threshold *= 0.999; + + // actually flag cells + refine (criteria, top_threshold); + coarsen (criteria, bottom_threshold); }; @@ -2978,20 +3027,30 @@ void Triangulation::prepare_refinement () { void TriangulationLevel<0>::reserve_space (const unsigned int total_cells, const unsigned int dimension) { - refine_flags.reserve (total_cells); - refine_flags.insert (refine_flags.end(), - total_cells - refine_flags.size(), - false); - - coarsen_flags.reserve (total_cells); - coarsen_flags.insert (coarsen_flags.end(), - total_cells - coarsen_flags.size(), - false); - - neighbors.reserve (total_cells*(2*dimension)); - neighbors.insert (neighbors.end(), - total_cells*(2*dimension) - neighbors.size(), - make_pair(-1,-1)); + // we need space for total_cells + // cells. Maybe we have more already + // with those cells which are unused, + // so only allocate new space if needed. + // + // note that all arrays should have equal + // sizes (checked by #monitor_memory# + if (total_cells > refine_flags.size()) + { + refine_flags.reserve (total_cells); + refine_flags.insert (refine_flags.end(), + total_cells - refine_flags.size(), + false); + + coarsen_flags.reserve (total_cells); + coarsen_flags.insert (coarsen_flags.end(), + total_cells - coarsen_flags.size(), + false); + + neighbors.reserve (total_cells*(2*dimension)); + neighbors.insert (neighbors.end(), + total_cells*(2*dimension) - neighbors.size(), + make_pair(-1,-1)); + }; }; @@ -3026,35 +3085,39 @@ void TriangulationLevel<1>::reserve_space (const unsigned int new_lines) { for (; u!=e; ++u) ++used_lines; - unsigned int new_size = used_lines + new_lines; + const unsigned int new_size = used_lines + new_lines; + // same as in #reserve_space<0>#: only + // allocate space if necessary + if (new_size>lines.lines.size()) + { // cout << " lines: pre: siz=" << lines.lines.size() << ", cap=" << lines.lines.capacity(); - lines.lines.reserve (new_size); + lines.lines.reserve (new_size); // cout << " inter: siz=" << lines.lines.size() << ", cap=" << lines.lines.capacity() // << " (newsize=" << new_size << ")"; - lines.lines.insert (lines.lines.end(), new_size-lines.lines.size(), Line()); + lines.lines.insert (lines.lines.end(), new_size-lines.lines.size(), Line()); // cout << " post: siz=" << lines.lines.size() << ", cap=" << lines.lines.capacity() << endl; // cout << " used : pre: siz=" << lines.used.size() << ", cap=" << lines.used.capacity(); - lines.used.reserve (new_size); + lines.used.reserve (new_size); // cout << " inter: siz=" << lines.used.size() << ", cap=" << lines.used.capacity() // << " (newsize=" << new_size << ")"; - lines.used.insert (lines.used.end(), new_size-lines.used.size(), false); + lines.used.insert (lines.used.end(), new_size-lines.used.size(), false); // cout << " post: siz=" << lines.used.size() << ", cap=" << lines.used.capacity() << endl; - lines.user_flags.reserve (new_size); - lines.user_flags.insert (lines.user_flags.end(), - new_size-lines.user_flags.size(), false); - - lines.children.reserve (new_size); - lines.children.insert (lines.children.end(), new_size-lines.children.size(), - -1); - - lines.material_id.reserve (new_size); - lines.material_id.insert (lines.material_id.end(), - new_size-lines.material_id.size(), - 255); - + lines.user_flags.reserve (new_size); + lines.user_flags.insert (lines.user_flags.end(), + new_size-lines.user_flags.size(), false); + + lines.children.reserve (new_size); + lines.children.insert (lines.children.end(), new_size-lines.children.size(), + -1); + + lines.material_id.reserve (new_size); + lines.material_id.insert (lines.material_id.end(), + new_size-lines.material_id.size(), + 255); + }; }; @@ -3085,8 +3148,6 @@ void TriangulationLevel<1>::monitor_memory (const unsigned int true_dimension) c ExcMemoryInexact (lines.lines.size(), lines.children.size())); Assert (lines.lines.size() == lines.material_id.size(), ExcMemoryInexact (lines.lines.size(), lines.material_id.size())); - Assert (lines.used[lines.used.size()-1]==true , - ExcUnusedMemoryAtEnd()); TriangulationLevel<0>::monitor_memory (true_dimension); }; @@ -3103,26 +3164,30 @@ void TriangulationLevel<2>::reserve_space (const unsigned int new_quads) { for (; u!=e; ++u) ++used_quads; - unsigned int new_size = used_quads + new_quads; - - quads.quads.reserve (new_size); - quads.quads.insert (quads.quads.end(), new_size-quads.quads.size(), Quad()); - - quads.used.reserve (new_size); - quads.used.insert (quads.used.end(), new_size-quads.used.size(), false); + const unsigned int new_size = used_quads + new_quads; + + // see above... + if (new_size>quads.quads.size()) + { + quads.quads.reserve (new_size); + quads.quads.insert (quads.quads.end(), new_size-quads.quads.size(), Quad()); + + quads.used.reserve (new_size); + quads.used.insert (quads.used.end(), new_size-quads.used.size(), false); - quads.user_flags.reserve (new_size); - quads.user_flags.insert (quads.user_flags.end(), - new_size-quads.user_flags.size(), false); + quads.user_flags.reserve (new_size); + quads.user_flags.insert (quads.user_flags.end(), + new_size-quads.user_flags.size(), false); - quads.children.reserve (new_size); - quads.children.insert (quads.children.end(), new_size-quads.children.size(), - -1); - - quads.material_id.reserve (new_size); - quads.material_id.insert (quads.material_id.end(), - new_size-quads.material_id.size(), - 255); + quads.children.reserve (new_size); + quads.children.insert (quads.children.end(), new_size-quads.children.size(), + -1); + + quads.material_id.reserve (new_size); + quads.material_id.insert (quads.material_id.end(), + new_size-quads.material_id.size(), + 255); + }; }; @@ -3152,8 +3217,6 @@ void TriangulationLevel<2>::monitor_memory (const unsigned int true_dimension) c ExcMemoryInexact (quads.quads.size(), quads.children.size())); Assert (quads.quads.size() == quads.material_id.size(), ExcMemoryInexact (quads.quads.size(), quads.material_id.size())); - Assert (quads.used[quads.used.size()-1]==true , - ExcUnusedMemoryAtEnd()); TriangulationLevel<1>::monitor_memory (true_dimension); };