]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Changes related to the implementation of coarsening.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 12 Jun 1998 08:42:37 +0000 (08:42 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 12 Jun 1998 08:42:37 +0000 (08:42 +0000)
git-svn-id: https://svn.dealii.org/trunk@396 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/grid/tria.cc

index c2bf2059831d251ee0be1af145d90cbecac3e9c9..803ce6d472c4833d62dcadf5b4241a0c83b18f7e 100644 (file)
@@ -263,7 +263,7 @@ class ConstraintMatrix {
                                      */
     DeclException1 (ExcLineExists,
                    unsigned int,
-                   << "The lines " << arg1
+                   << "The line " << arg1
                    << " which is to be created already exists.");
                                     /**
                                      * Exception
index 4afb2ad8ac3dc88851e2fe55e43a89bc33ddc571..137c8f56cdbdc48f3bdef5770e86fd8709c89a07 100644 (file)
@@ -614,6 +614,10 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      *
                                      * 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;
 
index 7d1a1c1e446d6a43b90ae0a2e609fbe5dd8cda3e..253344676fc27602fafef4e1f6451bc2fd0ce1b6 100644 (file)
@@ -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
+ *     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<dim> {
     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<dim> {
                                      * 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<dim> {
                                      * 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<dim> {
                                      *  #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<dim> {
                                       */
     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<dim># objects
index 7d60ac272ccd5fb8972cbf6dcaeed68500d27cb2..bcd59f9096cfa0f4d57ce9978c15719031fb27f7 100644 (file)
@@ -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<dim>::active_cell_iterator cell = tria->begin_active(),
+                                          endc = tria->end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      if (cell->face(face)->has_children()) 
+       cell->face(face)->set_user_flag();
+  
+
   
   line_iterator line = begin_line(),
                endl = end_line();
index 672055579ea9fc0a33afce52666b9fd1f731875f..4e5eb697a46b40da3be2cd803e75ee66567ec594 100644 (file)
@@ -729,6 +729,74 @@ void Triangulation<dim>::load_refine_flags (istream &in) {
 };
 
 
+
+template <int dim>
+void Triangulation<dim>::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; i<N/8+1; ++i) flags[i]=0;
+  
+  for (unsigned int position=0; cell!=endc; ++cell, ++position)
+    flags[position/8] |= (cell->coarsen_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<N/8+1; ++i) 
+    out << static_cast<unsigned int>(flags[i]) << " ";
+  
+  out << endl;
+  out << mn_tria_coarsen_flags_end << endl;
+
+  delete[] flags;
+};
+
+
+
+template <int dim>
+void Triangulation<dim>::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<N/8+1; ++i) 
+    {
+      in >> 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<dim>::refine (const dVector &criteria,
 
 
 template <int dim>
-void Triangulation<dim>::refine_fixed_number (const dVector &criteria,
-                                             const double   fraction) {
+void Triangulation<dim>::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; index<n_cells; ++cell, ++index)
+    if (criteria(index) <= threshold)
+      cell->set_coarsen_flag();
+};
+
+
+
+template <int dim>
+void Triangulation<dim>::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<int>(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<int>(top_fraction*criteria.size()),
                               1);
+  const int coarsen_cells = max(static_cast<int>(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<double>());
-
   refine (criteria, *(tmp.begin() + refine_cells));
+
+  nth_element (tmp.begin(), tmp.begin()+tmp.size()-coarsen_cells,
+              tmp.end(),
+              greater<double>());
+  coarsen (criteria, *(tmp.begin() + tmp.size() - coarsen_cells));
 };
 
 
@@ -1940,35 +2032,16 @@ double sqr(double a) {
 
 
 template <int dim>
-void Triangulation<dim>::refine_fixed_fraction (const dVector     &criteria,
-                                               const double       fraction_of_error,
-                                               const unsigned int n_sorting_steps) {
+void
+Triangulation<dim>::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<int>(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<dim>::refine_fixed_fraction (const dVector     &criteria,
   const double total_error = tmp.l1_norm();
   
   dVector partial_sums(criteria.size());
-  for (unsigned int part=0; part<n_sorting_parts; ++part)
-    {
-                                      // partially sort next part of range
-      partial_sort (tmp.begin()+part*cells_per_part,
-                   tmp.begin()+(part+1)*cells_per_part,
-                   tmp.end(),
-                   greater<double>());
-                                      // 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_sums.begin()+(part+1)*cells_per_part,
-                 ExcInternalError());
-
-                                          // now find the corresponding
-                                          // criterion and call refine
-         if (threshold_ptr==partial_sums.begin())
-           refine (criteria, sqrt(*threshold_ptr));
-         else
-           refine (criteria,
-                                                    // revert partial sum into
-                                                    // a single value
-                   sqrt(*threshold_ptr - *(threshold_ptr-1)));
-         return;
-       };
-    };
-                                  // this should not have happened: when
-                                  // we come to this point, we have either
-                                  // used more cells than the given fraction
-                                  // to reach the fraction_of_error, or
-                                  // something has gone terribly wrong.
+  sort (tmp.begin(), tmp.end(), greater<double>());
+  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<dim>::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);
 };

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.