From a210889d74ef6352c7d732296361b81c9778e5da Mon Sep 17 00:00:00 2001 From: Spencer Patty Date: Thu, 11 Feb 2016 15:49:25 -0600 Subject: [PATCH] add code + documentation for support of vertices being perturbed to GridTools::build_triangulation_from_patches --- include/deal.II/grid/grid_tools.h | 7 ++++ source/grid/grid_tools.cc | 58 ++++++++++++++++++++++++++----- 2 files changed, 56 insertions(+), 9 deletions(-) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 9d64822122..ca673f41fc 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1019,6 +1019,13 @@ namespace GridTools * identify the cells in the output Triangulation and corresponding cells in * the input list. * + * The function copies the location of vertices of cells from the cells of the + * source triangulation to the triangulation that is built from the list of + * patch cells. This adds support for triangulations which have been + * perturbed or smoothed in some manner which makes the triangulation + * deviate from the standard deal.ii refinement strategy of placing new + * vertices at midpoints of faces or edges. + * * The operation implemented by this function is frequently used in the * definition of error estimators that need to solve "local" problems on * each cell and its neighbors. A similar construction is necessary in the diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index eb61d2164c..123a2e28f6 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2932,9 +2932,9 @@ next_cell: for (uniform_cell=uniform_cells.begin(); uniform_cell!=uniform_cells.end(); ++uniform_cell) { bool repeat_vertex; - for (unsigned int j=0; j< GeometryInfo::vertices_per_cell; ++j) + for (unsigned int v=0; v::vertices_per_cell; ++v) { - Point position=(*uniform_cell)->vertex (j); + Point position=(*uniform_cell)->vertex (v); repeat_vertex=false; for (unsigned int m=0; m::vertices_per_cell; ++v) + active_tria_cell->vertex(v) = patch[i]->vertex(v); + + Assert(active_tria_cell->center().distance(patch_to_global_tria_map_tmp[active_tria_cell]->center()) + <=1e-15*active_tria_cell->diameter(), ExcInternalError()); + active_tria_cell->set_user_flag(); break; } @@ -3014,17 +3025,30 @@ next_cell: // children may be added into the map, instead // these children may not yet be in the map - for (unsigned int c=0; c< cell ->n_children(); ++c) + for (unsigned int c=0; cn_children(); ++c) { if (patch_to_global_tria_map_tmp.find(cell->child(c)) == patch_to_global_tria_map_tmp.end()) { - patch_to_global_tria_map_tmp.insert (std::make_pair(cell ->child(c), + patch_to_global_tria_map_tmp.insert (std::make_pair(cell->child(c), patch_to_global_tria_map_tmp[cell]->child(c))); - Assert(cell->child(c)->center().distance( patch_to_global_tria_map_tmp[cell]->child(c)->center()) - <=1e-15*cell->child(c)->diameter(), - ExcInternalError()); + // One might be tempted to assert that the cell + // being added here has the same center as the + // equivalent cell in the global triangulation, + // but it may not be the case. For triangulations + // that have been perturbed or smoothed, the cell + // indices and levels may be the same, but the + // vertex locations may not. We adjust + // the vertices of the cells that have no + // children (ie the active cells) to be + // consistent with the global triangulation + // later on and add assertions at that time + // to guarantee the cells in the + // local_triangulation are physically at the same + // locations of the cells in the patch of the + // global triangulation. + } } // The parent cell whose children were added @@ -3037,6 +3061,22 @@ next_cell: } while (refinement_necessary); + + + // Last assertion check to make sure we have the right cells and centers + // in the map, and hence the correct vertices of the triangulation + for (typename Triangulation::cell_iterator + cell = local_triangulation.begin(); + cell != local_triangulation.end(); ++cell) + { + Assert(patch_to_global_tria_map_tmp.find(cell) != patch_to_global_tria_map_tmp.end(), + ExcInternalError() ); + + Assert(cell->center().distance( patch_to_global_tria_map_tmp[cell]->center())<=1e-15*cell->diameter(), + ExcInternalError()); + } + + typename std::map::cell_iterator, typename Container::cell_iterator>::iterator map_tmp_it = patch_to_global_tria_map_tmp.begin(),map_tmp_end = patch_to_global_tria_map_tmp.end(); -- 2.39.5