* 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
* 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
* 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
* 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
* 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}
* 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}
*
* };
* \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
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#.
*
* 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
* Execute both refinement and coarsening
* of the triangulation.
*/
- void execute_refinement_and_coarsening ();
+ void execute_coarsening_and_refinement ();
/*@}*/
/**
* #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);
+
/*@}*/
*/
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
};
+
+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 <>
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));
};
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
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);
};
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));
+ };
};
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);
+ };
};
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);
};
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);
+ };
};
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);
};