From 6c812b019ec470d2cdc307c99e542208fc0a08e5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 24 Nov 2016 09:55:41 -0700 Subject: [PATCH] Clean up some p4est-related code. Specifically, commit daf3146 (via #3625) introduced a dimension dependent dispatch mechanism that allowed writing some code in a more generic way. The idea was right, but the approach duplicated the dispatching because we already have a dispatching mechanism via the internal::p4est::functions classes. This patch simply merges the two approaches. It also allows to move some code out of distributed/tria.cc into distributed/p4est_wrappers.cc, at the cost of some code churn. The patch does not actually change any kind of functionality -- it just moves things. --- include/deal.II/distributed/p4est_wrappers.h | 38 +++ source/distributed/p4est_wrappers.cc | 275 +++++++++++++++++ source/distributed/p4est_wrappers.inst.in | 25 ++ source/distributed/tria.cc | 307 +------------------ 4 files changed, 342 insertions(+), 303 deletions(-) diff --git a/include/deal.II/distributed/p4est_wrappers.h b/include/deal.II/distributed/p4est_wrappers.h index 33deee0b0a..795a561947 100644 --- a/include/deal.II/distributed/p4est_wrappers.h +++ b/include/deal.II/distributed/p4est_wrappers.h @@ -33,8 +33,21 @@ # include # include +#include +#include + + DEAL_II_NAMESPACE_OPEN +namespace parallel +{ + namespace distributed + { + template class Triangulation; + } +} + + namespace internal { namespace p4est @@ -272,6 +285,11 @@ namespace internal static size_t (&connectivity_memory_used) (types<2>::connectivity *p4est); + template + static void iterate(dealii::internal::p4est::types<2>::forest *parallel_forest, + dealii::internal::p4est::types<2>::ghost *parallel_ghost, + void *user_data); + static const unsigned int max_level = P4EST_MAXLEVEL; }; @@ -454,6 +472,11 @@ namespace internal static size_t (&connectivity_memory_used) (types<3>::connectivity *p4est); + template + static void iterate(dealii::internal::p4est::types<3>::forest *parallel_forest, + dealii::internal::p4est::types<3>::ghost *parallel_ghost, + void *user_data); + static const unsigned int max_level = P8EST_MAXLEVEL; }; @@ -539,6 +562,21 @@ namespace internal bool tree_exists_locally (const typename types::forest *parallel_forest, const typename types::topidx coarse_grid_cell); + + + /** + * This is the callback data structure used to fill + * vertices_with_ghost_neighbors via the p4est_iterate tool + */ + template + struct FindGhosts + { + typename dealii::parallel::distributed::Triangulation *triangulation; + sc_array_t *subids; + std::map > *vertices_with_ghost_neighbors; + }; + + } } diff --git a/source/distributed/p4est_wrappers.cc b/source/distributed/p4est_wrappers.cc index d2e6d4c614..c88da22d20 100644 --- a/source/distributed/p4est_wrappers.cc +++ b/source/distributed/p4est_wrappers.cc @@ -1,4 +1,5 @@ #include +#include DEAL_II_NAMESPACE_OPEN @@ -8,6 +9,251 @@ namespace internal { namespace p4est { + namespace + { + template + typename dealii::Triangulation::cell_iterator + cell_from_quad + (typename dealii::parallel::distributed::Triangulation *triangulation, + typename dealii::internal::p4est::types::topidx treeidx, + typename dealii::internal::p4est::types::quadrant &quad) + { + int i, l = quad.level; + int child_id; + dealii::types::global_dof_index dealii_index = + triangulation->get_p4est_tree_to_coarse_cell_permutation()[treeidx]; + + for (i = 0; i < l; i++) + { + typename dealii::Triangulation::cell_iterator cell (triangulation, i, dealii_index); + child_id = dealii::internal::p4est::functions::quadrant_ancestor_id (&quad, i + 1); + Assert (cell->has_children (), ExcMessage ("p4est quadrant does not correspond to a cell!")); + dealii_index = cell->child_index(child_id); + } + + typename dealii::Triangulation::cell_iterator out_cell (triangulation, l, dealii_index); + + return out_cell; + } + + + /** At a corner (vertex), determine if any of the neighboring cells are + * ghosts. If there are, find out their subdomain ids, and if this is a + * local vertex, then add these subdomain ids to the map + * vertices_with_ghost_neighbors of that index + */ + template + void + find_ghosts_corner + (typename dealii::internal::p4est::iter::corner_info *info, + void *user_data) + { + int i, j; + int nsides = info->sides.elem_count; + typename dealii::internal::p4est::iter::corner_side *sides = + (typename dealii::internal::p4est::iter::corner_side *) + (info->sides.array); + FindGhosts *fg = static_cast *>(user_data); + sc_array_t *subids = fg->subids; + typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; + int nsubs; + dealii::types::subdomain_id *subdomain_ids; + std::map > + *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; + + subids->elem_count = 0; + for (i = 0; i < nsides; i++) + { + if (sides[i].is_ghost) + { + typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad)); + Assert (cell->is_ghost(), ExcMessage ("ghost quad did not find ghost cell")); + dealii::types::subdomain_id *subid = + static_cast(sc_array_push (subids)); + *subid = cell->subdomain_id(); + } + } + + if (!subids->elem_count) + { + return; + } + + nsubs = (int) subids->elem_count; + subdomain_ids = (dealii::types::subdomain_id *) (subids->array); + + for (i = 0; i < nsides; i++) + { + if (!sides[i].is_ghost) + { + typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad)); + + Assert (!cell->is_ghost(), ExcMessage ("local quad found ghost cell")); + + for (j = 0; j < nsubs; j++) + { + (*vertices_with_ghost_neighbors)[cell->vertex_index(sides[i].corner)] + .insert (subdomain_ids[j]); + } + } + } + + subids->elem_count = 0; + } + + /** Similar to find_ghosts_corner, but for the hanging vertex in the + * middle of an edge + */ + template + void + find_ghosts_edge + (typename dealii::internal::p4est::iter::edge_info *info, + void *user_data) + { + int i, j, k; + int nsides = info->sides.elem_count; + typename dealii::internal::p4est::iter::edge_side *sides = + (typename dealii::internal::p4est::iter::edge_side *) + (info->sides.array); + FindGhosts *fg = static_cast *>(user_data); + sc_array_t *subids = fg->subids; + typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; + int nsubs; + dealii::types::subdomain_id *subdomain_ids; + std::map > + *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; + + subids->elem_count = 0; + for (i = 0; i < nsides; i++) + { + if (sides[i].is_hanging) + { + for (j = 0; j < 2; j++) + { + if (sides[i].is.hanging.is_ghost[j]) + { + typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); + dealii::types::subdomain_id *subid = + static_cast(sc_array_push (subids)); + *subid = cell->subdomain_id(); + } + } + } + } + + if (!subids->elem_count) + { + return; + } + + nsubs = (int) subids->elem_count; + subdomain_ids = (dealii::types::subdomain_id *) (subids->array); + + for (i = 0; i < nsides; i++) + { + if (sides[i].is_hanging) + { + for (j = 0; j < 2; j++) + { + if (!sides[i].is.hanging.is_ghost[j]) + { + typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); + + for (k = 0; k < nsubs; k++) + { + (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_edge_corners[sides[i].edge][1^j])] + .insert (subdomain_ids[k]); + } + } + } + } + } + + subids->elem_count = 0; + } + + /** Similar to find_ghosts_corner, but for the hanging vertex in the + * middle of a face + */ + template + void + find_ghosts_face + (typename dealii::internal::p4est::iter::face_info *info, + void *user_data) + { + int i, j, k; + int nsides = info->sides.elem_count; + typename dealii::internal::p4est::iter::face_side *sides = + (typename dealii::internal::p4est::iter::face_side *) + (info->sides.array); + FindGhosts *fg = static_cast *>(user_data); + sc_array_t *subids = fg->subids; + typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; + int nsubs; + dealii::types::subdomain_id *subdomain_ids; + std::map > + *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; + int limit = (dim == 2) ? 2 : 4; + + subids->elem_count = 0; + for (i = 0; i < nsides; i++) + { + if (sides[i].is_hanging) + { + for (j = 0; j < limit; j++) + { + if (sides[i].is.hanging.is_ghost[j]) + { + typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); + dealii::types::subdomain_id *subid = + static_cast(sc_array_push (subids)); + *subid = cell->subdomain_id(); + } + } + } + } + + if (!subids->elem_count) + { + return; + } + + nsubs = (int) subids->elem_count; + subdomain_ids = (dealii::types::subdomain_id *) (subids->array); + + for (i = 0; i < nsides; i++) + { + if (sides[i].is_hanging) + { + for (j = 0; j < limit; j++) + { + if (!sides[i].is.hanging.is_ghost[j]) + { + typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); + + for (k = 0; k < nsubs; k++) + { + if (dim == 2) + { + (*vertices_with_ghost_neighbors)[cell->vertex_index(p4est_face_corners[sides[i].face][(limit - 1)^j])] + .insert (subdomain_ids[k]); + } + else + { + (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_face_corners[sides[i].face][(limit - 1)^j])] + .insert (subdomain_ids[k]); + } + } + } + } + } + } + + subids->elem_count = 0; + } + } + + int (&functions<2>::quadrant_compare) (const void *v1, const void *v2) = p4est_quadrant_compare; @@ -191,8 +437,22 @@ namespace internal size_t (&functions<2>::connectivity_memory_used) (types<2>::connectivity *p4est) = p4est_connectivity_memory_used; + template + void functions<2>::iterate(dealii::internal::p4est::types<2>::forest *parallel_forest, + dealii::internal::p4est::types<2>::ghost *parallel_ghost, + void *user_data) + { + p4est_iterate (parallel_forest, parallel_ghost, user_data, + NULL, + find_ghosts_face<2,spacedim>, + find_ghosts_corner<2,spacedim>); + } + const unsigned int functions<2>::max_level; + + + int (&functions<3>::quadrant_compare) (const void *v1, const void *v2) = p8est_quadrant_compare; @@ -379,8 +639,23 @@ namespace internal size_t (&functions<3>::connectivity_memory_used) (types<3>::connectivity *p4est) = p8est_connectivity_memory_used; + template + void functions<3>::iterate(dealii::internal::p4est::types<3>::forest *parallel_forest, + dealii::internal::p4est::types<3>::ghost *parallel_ghost, + void *user_data) + { + p8est_iterate (parallel_forest, parallel_ghost, user_data, + NULL, + find_ghosts_face<3,spacedim>, + find_ghosts_edge<3,spacedim>, + find_ghosts_corner<3,spacedim>); + } + const unsigned int functions<3>::max_level; + + + template void init_quadrant_children diff --git a/source/distributed/p4est_wrappers.inst.in b/source/distributed/p4est_wrappers.inst.in index 1806dc8201..464352a06c 100644 --- a/source/distributed/p4est_wrappers.inst.in +++ b/source/distributed/p4est_wrappers.inst.in @@ -58,3 +58,28 @@ for (deal_II_dimension : DIMENSIONS) \} # endif // DEAL_II_WITH_P4EST } + + +for (deal_II_dimension, deal_II_space_dimension : DIMENSIONS) +{ +# ifdef DEAL_II_WITH_P4EST + + namespace internal + \{ + namespace p4est + \{ +#if deal_II_dimension > 1 +#if deal_II_space_dimension >= deal_II_dimension + template + void + functions:: + iterate(dealii::internal::p4est::types::forest *parallel_forest, + dealii::internal::p4est::types::ghost *parallel_ghost, + void *user_data); +#endif +#endif + \} + \} +# endif // DEAL_II_WITH_P4EST +} + diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 646375d75c..3e175ffdb4 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -3380,280 +3380,6 @@ namespace parallel } - - namespace - { - /** - * This is the callback data structure used to fill - * vertices_with_ghost_neighbors via the p4est_iterate tool - */ - template - struct FindGhosts - { - typename dealii::parallel::distributed::Triangulation *triangulation; - sc_array_t *subids; - std::map > - *vertices_with_ghost_neighbors; - }; - - /** At a corner (vertex), determine if any of the neighboring cells are - * ghosts. If there are, find out their subdomain ids, and if this is a - * local vertex, then add these subdomain ids to the map - * vertices_with_ghost_neighbors of that index - */ - template - void - find_ghosts_corner - (typename dealii::internal::p4est::iter::corner_info *info, - void *user_data) - { - int i, j; - int nsides = info->sides.elem_count; - typename dealii::internal::p4est::iter::corner_side *sides = - (typename dealii::internal::p4est::iter::corner_side *) - (info->sides.array); - FindGhosts *fg = static_cast *>(user_data); - sc_array_t *subids = fg->subids; - typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; - int nsubs; - dealii::types::subdomain_id *subdomain_ids; - std::map > - *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; - - subids->elem_count = 0; - for (i = 0; i < nsides; i++) - { - if (sides[i].is_ghost) - { - typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad)); - Assert (cell->is_ghost(), ExcMessage ("ghost quad did not find ghost cell")); - dealii::types::subdomain_id *subid = - static_cast(sc_array_push (subids)); - *subid = cell->subdomain_id(); - } - } - - if (!subids->elem_count) - { - return; - } - - nsubs = (int) subids->elem_count; - subdomain_ids = (dealii::types::subdomain_id *) (subids->array); - - for (i = 0; i < nsides; i++) - { - if (!sides[i].is_ghost) - { - typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad)); - - Assert (!cell->is_ghost(), ExcMessage ("local quad found ghost cell")); - - for (j = 0; j < nsubs; j++) - { - (*vertices_with_ghost_neighbors)[cell->vertex_index(sides[i].corner)] - .insert (subdomain_ids[j]); - } - } - } - - subids->elem_count = 0; - } - - /** Similar to find_ghosts_corner, but for the hanging vertex in the - * middle of an edge - */ - template - void - find_ghosts_edge - (typename dealii::internal::p4est::iter::edge_info *info, - void *user_data) - { - int i, j, k; - int nsides = info->sides.elem_count; - typename dealii::internal::p4est::iter::edge_side *sides = - (typename dealii::internal::p4est::iter::edge_side *) - (info->sides.array); - FindGhosts *fg = static_cast *>(user_data); - sc_array_t *subids = fg->subids; - typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; - int nsubs; - dealii::types::subdomain_id *subdomain_ids; - std::map > - *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; - - subids->elem_count = 0; - for (i = 0; i < nsides; i++) - { - if (sides[i].is_hanging) - { - for (j = 0; j < 2; j++) - { - if (sides[i].is.hanging.is_ghost[j]) - { - typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); - dealii::types::subdomain_id *subid = - static_cast(sc_array_push (subids)); - *subid = cell->subdomain_id(); - } - } - } - } - - if (!subids->elem_count) - { - return; - } - - nsubs = (int) subids->elem_count; - subdomain_ids = (dealii::types::subdomain_id *) (subids->array); - - for (i = 0; i < nsides; i++) - { - if (sides[i].is_hanging) - { - for (j = 0; j < 2; j++) - { - if (!sides[i].is.hanging.is_ghost[j]) - { - typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); - - for (k = 0; k < nsubs; k++) - { - (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_edge_corners[sides[i].edge][1^j])] - .insert (subdomain_ids[k]); - } - } - } - } - } - - subids->elem_count = 0; - } - - /** Similar to find_ghosts_corner, but for the hanging vertex in the - * middle of a face - */ - template - void - find_ghosts_face - (typename dealii::internal::p4est::iter::face_info *info, - void *user_data) - { - int i, j, k; - int nsides = info->sides.elem_count; - typename dealii::internal::p4est::iter::face_side *sides = - (typename dealii::internal::p4est::iter::face_side *) - (info->sides.array); - FindGhosts *fg = static_cast *>(user_data); - sc_array_t *subids = fg->subids; - typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; - int nsubs; - dealii::types::subdomain_id *subdomain_ids; - std::map > - *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; - int limit = (dim == 2) ? 2 : 4; - - subids->elem_count = 0; - for (i = 0; i < nsides; i++) - { - if (sides[i].is_hanging) - { - for (j = 0; j < limit; j++) - { - if (sides[i].is.hanging.is_ghost[j]) - { - typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); - dealii::types::subdomain_id *subid = - static_cast(sc_array_push (subids)); - *subid = cell->subdomain_id(); - } - } - } - } - - if (!subids->elem_count) - { - return; - } - - nsubs = (int) subids->elem_count; - subdomain_ids = (dealii::types::subdomain_id *) (subids->array); - - for (i = 0; i < nsides; i++) - { - if (sides[i].is_hanging) - { - for (j = 0; j < limit; j++) - { - if (!sides[i].is.hanging.is_ghost[j]) - { - typename dealii::parallel::distributed::Triangulation::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); - - for (k = 0; k < nsubs; k++) - { - if (dim == 2) - { - (*vertices_with_ghost_neighbors)[cell->vertex_index(p4est_face_corners[sides[i].face][(limit - 1)^j])] - .insert (subdomain_ids[k]); - } - else - { - (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_face_corners[sides[i].face][(limit - 1)^j])] - .insert (subdomain_ids[k]); - } - } - } - } - } - } - - subids->elem_count = 0; - } - - - /** - * This struct wraps the p4est_iterate / p8est_iterate in a templated - * fashion so that we can avoid the compiler complaining about invalid - * types being referred, which would be the case inside switch or if - * statements. - */ - template - struct pXest_iterate_wrapper - { - static void iterate(typename dealii::internal::p4est::types::forest *parallel_forest, - typename dealii::internal::p4est::types::ghost *parallel_ghost, - void *user_data); - }; - - template - struct pXest_iterate_wrapper<2,spacedim> - { - static void iterate(dealii::internal::p4est::types<2>::forest *parallel_forest, - dealii::internal::p4est::types<2>::ghost *parallel_ghost, - void *user_data) - { - p4est_iterate (parallel_forest, parallel_ghost, user_data, - NULL, find_ghosts_face<2,spacedim>, - find_ghosts_corner<2,spacedim>); - } - }; - - template - struct pXest_iterate_wrapper<3,spacedim> - { - static void iterate(dealii::internal::p4est::types<3>::forest *parallel_forest, - dealii::internal::p4est::types<3>::ghost *parallel_ghost, - void *user_data) - { - p8est_iterate (parallel_forest, parallel_ghost, user_data, - NULL, find_ghosts_face<3,spacedim>, find_ghosts_edge<3,spacedim>, - find_ghosts_corner<3,spacedim>); - } - }; - } - - /** * Determine the neighboring subdomains that are adjacent to each vertex. * This is achieved via the p4est_iterate/p8est_iterate tool @@ -3667,14 +3393,14 @@ namespace parallel { Assert (dim>1, ExcNotImplemented()); - FindGhosts fg; + dealii::internal::p4est::FindGhosts fg; fg.subids = sc_array_new (sizeof (dealii::types::subdomain_id)); fg.triangulation = this; fg.vertices_with_ghost_neighbors = &vertices_with_ghost_neighbors; - pXest_iterate_wrapper::iterate (this->parallel_forest, - this->parallel_ghost, - static_cast(&fg)); + dealii::internal::p4est::functions::template iterate (this->parallel_forest, + this->parallel_ghost, + static_cast(&fg)); sc_array_destroy (fg.subids); } @@ -4138,31 +3864,6 @@ namespace parallel return weights; } - template - typename dealii::Triangulation::cell_iterator - cell_from_quad - (typename dealii::parallel::distributed::Triangulation *triangulation, - typename dealii::internal::p4est::types::topidx treeidx, - typename dealii::internal::p4est::types::quadrant &quad) - { - int i, l = quad.level; - int child_id; - types::global_dof_index dealii_index = - triangulation->get_p4est_tree_to_coarse_cell_permutation()[treeidx]; - - for (i = 0; i < l; i++) - { - typename dealii::Triangulation::cell_iterator cell (triangulation, i, dealii_index); - child_id = dealii::internal::p4est::functions::quadrant_ancestor_id (&quad, i + 1); - Assert (cell->has_children (), ExcMessage ("p4est quadrant does not correspond to a cell!")); - dealii_index = cell->child_index(child_id); - } - - typename dealii::Triangulation::cell_iterator out_cell (triangulation, l, dealii_index); - - return out_cell; - } - template -- 2.39.5