From: Wolfgang Bangerth Date: Thu, 4 Jun 1998 15:20:14 +0000 (+0000) Subject: Add more smoothing functionality to the refinement of triangulations X-Git-Tag: v8.0.0~22875 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=97ebc9198237ec5a5a0b937bfc985cfe920016c1;p=dealii.git Add more smoothing functionality to the refinement of triangulations git-svn-id: https://svn.dealii.org/trunk@379 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 781f87189a..53474f58ea 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -543,6 +543,20 @@ class TriaDimensionInfo<2> { +/** + * Declare some symbolic names for mesh smoothing algorithms. The meaning of + * these flags is documented in the #Triangulation# class. + */ +enum MeshSmoothing { + none = 0x0, + limit_level_difference_at_vertices = 0x1, + eliminate_unrefined_islands = 0x2, + maximum_smoothing = 0xffffffff +}; + + + + /*------------------------------------------------------------------------*/ @@ -840,7 +854,7 @@ class TriaDimensionInfo<2> { * interface between regions of different materials. * * - * \subsection{Refinement of a triangulation} + * \subsection{Refinement and 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 @@ -861,63 +875,6 @@ class TriaDimensionInfo<2> { * structures and algorithms much much easier. To be honest, this is mostly * an algorithmic step than one needed by the finite element method. * - * However, some degradation of approximation properties has been observed - * for grids which were refined more than once across a face, as descibed - * above, so there is also a practical justification for the above. - * It can also be shown, that such degradation occurs if the - * triangulation contains vertices which are member of cells with levels - * differing by more than one. One such example is the following: - * \begin{verbatim} - * | | | | - * x-----x-----x--x--x-- - * | | | | | - * | | x--x--x - * | | | | | - * x-----x-----x--x--x-- - * | | | - * | | | - * | | | - * | x-----x-- - * | | | - * | | | - * | | | - * x-----------x-----x-- - * \end{verbatim} - * It seems that in two space dimensions, the maximum jump in levels between - * cells sharing a common vertex is two (as in the example above). This is - * not true if more than four cells meet at a vertex. It is not uncommon - * that a coarse (initial) mesh contains vertices at which six or even eight - * cells meet, when small features of the domain have to be resolved even on - * the coarsest mesh. In that case, the maximum difference in levels is - * three or four, respectively. The problem gets even worse in three space - * dimensions. - * - * Looking at an interpolation of the second derivative of the finite - * element solution (asuming bilinear finite elements), one sees that the - * numerical solution is almost totally wrong, compared with the true second - * derivative. Indeed, on regular meshes, there exist sharp estimations that - * the $H^2$-error is only $O(1)$, so we should not be suprised; however, the - * numerical solution may show a value for the second derivative which may - * be a factor of ten away from the true value. These problems are located - * on the small cell adjacent to the center vertex, where cells of - * non-subsequent levels meet, as well as on the upper and right neighbor - * of this cell (but with a less degree of deviation from the true value). - * - * Due to the approximational problems described above, the - * #Triangulation# constructor takes an argument specifying whether a - * smoothing step shall be performed on the grid each time #execute_refinement# - * is called. The default is that such a step not be done, since this results - * in additional cells being produced, which may not be necessary in all - * cases. If switched on, calling #execute_refinement# results in - * flagging additional cells for refinement to avoid - * vertices as the ones mentioned. The algorithms for both regularisation - * and smoothing of triangulations are described below in the section on - * technical issues. The reason why this parameter must be given to the - * constructor rather than to #execute_refinement# is that it would result - * in algorithmic problems if you called #execute_refinement# once without - * and once with smoothing, since then in some refinement steps would need - * to be refined twice. - * * 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 @@ -1012,6 +969,102 @@ class TriaDimensionInfo<2> { * squares of the criteria on the cells. The criteria shall be positive. * * + * \subsection{Smoothing of a triangulation} + * + * Some degradation of approximation properties has been observed + * for grids which are too unstructured. Therefore, the #prepare_refinement# + * function which is automatically called by the #execute_refinement# function + * can do some smoothing of the triangulation. For this purpose the + * #Triangulation# constructor takes an argument specifying whether a + * smoothing step shall be performed on the grid each time #execute_refinement# + * is called. The default is that such a step not be done, since this results + * in additional cells being produced, which may not be necessary in all + * cases. If switched on, calling #execute_refinement# results in + * flagging additional cells for refinement to avoid + * vertices as the ones mentioned. The algorithms for both regularisation + * and smoothing of triangulations are described below in the section on + * technical issues. The reason why this parameter must be given to the + * constructor rather than to #execute_refinement# is that it would result + * in algorithmic problems if you called #execute_refinement# once without + * and once with smoothing, since then in some refinement steps would need + * to be refined twice. + * + * The parameter taken by the constructor is an integer which may be composed + * bitwise by the constants defined in the #enum MeshSmoothing#. The meaning + * of these constants is explained in the following: + * \begin{itemize} + * \item #limit_level_difference_at_vertices#: + * It can be shown, that degradation of approximation occurs if the + * triangulation contains vertices which are member of cells with levels + * differing by more than one. One such example is the following: + * \begin{verbatim} + * | | | | + * x-----x-----x--x--x-- + * | | | | | + * | | x--x--x + * | | | | | + * x-----x-----x--x--x-- + * | | | + * | | | + * | | | + * | x-----x-- + * | | | + * | | | + * | | | + * x-----------x-----x-- + * \end{verbatim} + * It seems that in two space dimensions, the maximum jump in levels between + * cells sharing a common vertex is two (as in the example above). This is + * not true if more than four cells meet at a vertex. It is not uncommon + * that a coarse (initial) mesh contains vertices at which six or even eight + * cells meet, when small features of the domain have to be resolved even on + * the coarsest mesh. In that case, the maximum difference in levels is + * three or four, respectively. The problem gets even worse in three space + * dimensions. + * + * Looking at an interpolation of the second derivative of the finite + * element solution (asuming bilinear finite elements), one sees that the + * numerical solution is almost totally wrong, compared with the true second + * derivative. Indeed, on regular meshes, there exist sharp estimations that + * the $H^2$-error is only $O(1)$, so we should not be suprised; however, the + * numerical solution may show a value for the second derivative which may + * be a factor of ten away from the true value. These problems are located + * on the small cell adjacent to the center vertex, where cells of + * non-subsequent levels meet, as well as on the upper and right neighbor + * of this cell (but with a less degree of deviation from the true value). + * + * If the smoothing indicator given to the constructor contains the bit for + * #limit_level_difference_at_vertices#, situations as the above one are + * eliminated by also marking the lower left cell for refinement. + * + * \item #eliminate_unrefined_islands#: + * Single cells which are not refined and are surrounded by cells which are + * refined usually also lead to a sharp decline in approximation properties + * locally. The reason is that the nodes on the faces between unrefined and + * refined cells are not real degrees of freedom but carry constraints. The + * patch without additional degrees of freedom is thus significantly larger + * then the unrefined cell itself. If in the parameter passed to the + * constructor the bit for #eliminate_unrefined_islands# is set, all cells + * which are not flagged for refinement but which are surrounded by more + * refined cells than unrefined cells are flagged for refinement. Cells + * which are not yet refined but flagged for that are accounted for the + * number of refined neighbors. Cells on the boundary are not accounted for + * at all. An unrefined island is, by this definition + * also a cell which (in 2D) is surrounded by three refined cells and one + * unrefined one, or one surrounded by two refined cells, one unrefined one + * and is at the boundary on one side. It is thus not a true island, as the + * name of the flag may indicate. However, no better name came to mind to + * the author by now. + * + * \item #maximum_smoothing#: + * This flag includes all the above ones and therefore combines all + * smoothing algorithms implemented. + * + * \item #none#: + * Select no smoothing at all. + * \end{itemize} + * + * * \subsection{Material and boundary information} * * Each line, quad, etc stores one byte of information denoting the material @@ -1183,31 +1236,44 @@ class TriaDimensionInfo<2> { * * \begin{itemize} * \item {\it Regularisation:} The algorithm walks over all cells checking - * whether the present cell is flagged for refinement and a neighbor of the - * present cell is refined once less than the present one. If so, flag the - * neighbor for refinement. Because of the induction above, there may be no - * neighbor with level two less than the present one. - * - * The neighbor thus flagged for refinement may induce more cells which need - * to be refined. However, such cells which need additional refinement always - * are on one level lower than the present one, so we can get away with only - * one sweep over all cells if we do the loop in the reverse way, starting - * with those on the highest level. This way, we may flag additional cells - * on lower levels, but if these induce more refinement needed, this is - * performed later on when we visit them in out backward running loop. - * - * \item {\it Smoothing:} First a list is set up which stores for each vertex - * the highest level one of the adjacent cells belongs to. Now, since we did - * smoothing in the previous refinement steps also, each cell may only have - * vertices with levels at most one greater than the level of the present - * cell. - * - * However, if we store the level plus one for cells marked for refinement, - * we may end up with cells which have vertices of level two greater than - * the cells level. We need to refine this cell also, and need thus also - * update the levels of its vertices. This itself may lead to cells needing - * refinement, but these are on lower levels, as above, which is why we - * may do all kinds of additional flagging in one loop only. + * whether the present cell is flagged for refinement and a neighbor of the + * present cell is refined once less than the present one. If so, flag the + * neighbor for refinement. Because of the induction above, there may be no + * neighbor with level two less than the present one. + * + * The neighbor thus flagged for refinement may induce more cells which need + * to be refined. However, such cells which need additional refinement always + * are on one level lower than the present one, so we can get away with only + * one sweep over all cells if we do the loop in the reverse way, starting + * with those on the highest level. This way, we may flag additional cells + * on lower levels, but if these induce more refinement needed, this is + * performed later on when we visit them in out backward running loop. + * + * \item {\it Smoothing:} + * \begin{itemize} + * \item #limit_level_difference_at_vertices#: + * First a list is set up which stores for each vertex + * the highest level one of the adjacent cells belongs to. Now, since we did + * smoothing in the previous refinement steps also, each cell may only have + * vertices with levels at most one greater than the level of the present + * cell. + * + * However, if we store the level plus one for cells marked for refinement, + * we may end up with cells which have vertices of level two greater than + * the cells level. We need to refine this cell also, and need thus also + * update the levels of its vertices. This itself may lead to cells needing + * refinement, but these are on lower levels, as above, which is why we + * may do all kinds of additional flagging in one loop only. + * + * \item #eliminate_unrefined_islands#: + * For each cell we count the number of neighbors which are refined or + * flagged for refinement. If this exceeds the total number of neighbors + * (which is the number of faces minus the number of faces of this cell + * which are located on the boundary), then this cell is flagged for + * refinement. Since this may lead to cells on the same level which also + * will need refinement, we will need additional loops of regularisation + * and smoothing over all cells until nothing changes any more. + * \end{itemize} * \end{itemize} * * Regularisation and smoothing are a bit complementary in that we check @@ -1287,7 +1353,7 @@ class Triangulation : public TriaDimensionInfo { * the first level of the hierarchy. * Do not create any cells. */ - Triangulation (const bool smooth_grid = false); + Triangulation (const MeshSmoothing smooth_grid = none); /** * Copy constructor. You should really @@ -2328,7 +2394,7 @@ class Triangulation : public TriaDimensionInfo { * the general doc of this class for * more information about this. */ - bool smooth_grid; + MeshSmoothing smooth_grid; // Friendship includes local classes. friend class TriaAccessor; diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index c91b3bdf10..41963fa5c9 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -15,7 +15,7 @@ template -Triangulation::Triangulation (const bool smooth_grid) : +Triangulation::Triangulation (const MeshSmoothing smooth_grid) : smooth_grid(smooth_grid) { static StraightBoundary default_boundary; @@ -2839,74 +2839,116 @@ void Triangulation::prepare_refinement () { // there are no such cells. if (dim>=2) { - // store highest level one of the cells - // adjacent to a vertex belongs to; do - // so only if mesh smoothing is - // required - vector vertex_level; - if (smooth_grid) + // store whether some cells were flagged + // or deflagged for refinement in this + // loop and loop until everything is + // settled. + bool mesh_changed_in_this_loop; + do { - vertex_level.resize (vertices.size(), 0); - active_cell_iterator cell = begin_active(), - endc = end(); - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; vertex::vertices_per_cell; - ++vertex) - if (cell->refine_flag_set()) - vertex_level[cell->vertex_index(vertex)] - = max (vertex_level[cell->vertex_index(vertex)], - cell->level()+1); - else - vertex_level[cell->vertex_index(vertex)] - = max (vertex_level[cell->vertex_index(vertex)], - cell->level()); - }; + mesh_changed_in_this_loop = false; + + // store highest level one of the cells + // adjacent to a vertex belongs to; do + // so only if mesh smoothing is + // required + vector vertex_level; + if (smooth_grid & limit_level_difference_at_vertices) + { + vertex_level.resize (vertices.size(), 0); + active_cell_iterator cell = begin_active(), + endc = end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + if (cell->refine_flag_set()) + vertex_level[cell->vertex_index(vertex)] + = max (vertex_level[cell->vertex_index(vertex)], + cell->level()+1); + else + vertex_level[cell->vertex_index(vertex)] + = max (vertex_level[cell->vertex_index(vertex)], + cell->level()); + }; - active_cell_iterator cell = last_active(), - endc = end(); + active_cell_iterator cell = last_active(), + endc = end(); - // loop over active cells - for (; cell != endc; --cell) - if (cell->refine_flag_set() == true) - { - // loop over neighbors of cell - for (unsigned int i=0; i::faces_per_cell; ++i) - if (cell->neighbor(i).state() == valid) - { - // regularisation? - if ((cell->neighbor_level(i) == cell->level()-1) - && - (cell->neighbor(i)->refine_flag_set() == false)) - cell->neighbor(i)->set_refine_flag(); - }; - } - else - // smoothing? - if (smooth_grid) - for (unsigned int vertex=0; - vertex::vertices_per_cell; ++vertex) - if (vertex_level[cell->vertex_index(vertex)] > - cell->level()+1) + // loop over active cells + for (; cell != endc; --cell) + if (cell->refine_flag_set() == true) + { + // loop over neighbors of cell + for (unsigned int i=0; i::faces_per_cell; ++i) + if (cell->neighbor(i).state() == valid) + { + // regularisation? + if ((cell->neighbor_level(i) == cell->level()-1) + && + (cell->neighbor(i)->refine_flag_set() == false)) + { + cell->neighbor(i)->set_refine_flag(); + mesh_changed_in_this_loop = true; + }; + }; + } + else + // smoothing? + if (smooth_grid & limit_level_difference_at_vertices) + for (unsigned int vertex=0; + vertex::vertices_per_cell; ++vertex) + if (vertex_level[cell->vertex_index(vertex)] > + cell->level()+1) + { + // if we did not make an + // error, the level diff + // should not be more than + // two + Assert (vertex_level[cell->vertex_index(vertex)] == + cell->level()+2, + ExcInternalError()); + + // refine cell and + // update vertex levels + cell->set_refine_flag(); + mesh_changed_in_this_loop = true; + for (unsigned int v=0; v::vertices_per_cell; + ++v) + vertex_level[cell->vertex_index(v)] + = max (vertex_level[cell->vertex_index(v)], + cell->level()+1); + }; + + // additional smoothing + if (smooth_grid & eliminate_unrefined_islands) + { + active_cell_iterator cell = begin_active(), + endc = end(); + for (; cell!=endc; ++cell) { - // if we did not make an - // error, the level diff - // should not be more than - // two - Assert (vertex_level[cell->vertex_index(vertex)] == - cell->level()+2, - ExcInternalError()); - - // refine cell and - // update vertex levels - cell->set_refine_flag(); - for (unsigned int v=0; v::vertices_per_cell; - ++v) - vertex_level[cell->vertex_index(v)] - = max (vertex_level[cell->vertex_index(v)], - cell->level()+1); + unsigned int refined_neighbors = 0, + unrefined_neighbors = 0; + for (unsigned int face=0; face::faces_per_cell; ++face) + if (!cell->at_boundary(face)) + // neighbor may only be on + // the same level or one + // level below because of + // the regularisation above + if ((cell->neighbor_level(face) == cell->level()) && + (cell->neighbor(face)->refine_flag_set())) + ++refined_neighbors; + else + ++unrefined_neighbors; + + if (unrefined_neighbors < refined_neighbors) + cell->set_refine_flag (); }; + }; + } + while (mesh_changed_in_this_loop == true); + }; - + // check whether a new level is needed raw_cell_iterator cell = begin_active (levels.size()-1),