From e65e8e981d7fb0c3877cf8c53c9b13eb47254b7b Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 6 Sep 2013 09:02:24 +0000 Subject: [PATCH] Added support for periodic boundary conditions on distributed meshes git-svn-id: https://svn.dealii.org/trunk@30624 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/distributed/tria.h | 52 +- deal.II/include/deal.II/dofs/dof_tools.h | 31 +- deal.II/source/distributed/tria.cc | 650 ++++++++++++++++--- deal.II/source/dofs/dof_handler_policy.cc | 29 +- deal.II/source/dofs/dof_tools_constraints.cc | 47 +- 5 files changed, 671 insertions(+), 138 deletions(-) diff --git a/deal.II/include/deal.II/distributed/tria.h b/deal.II/include/deal.II/distributed/tria.h index 6275739a61..a1088fba3b 100644 --- a/deal.II/include/deal.II/distributed/tria.h +++ b/deal.II/include/deal.II/distributed/tria.h @@ -25,7 +25,9 @@ #include #include +#include +#include #include #include #include @@ -694,9 +696,30 @@ namespace parallel * in hierarchical ordering is the ith deal cell starting * from begin(0). */ - const std::vector & + const std::vector & get_p4est_tree_to_coarse_cell_permutation() const; + + /** + * Join faces in the p4est forest due to periodic boundary conditions. + * + * The entries in the std::vector should have the form + * std_cxx1x::tuple. + * + * The vector can be filled by the function + * DoFTools::identify_periodic_face_pairs. + * + * @note Before this function can be used the triangulation has to be + * initialized and must not be refined. + * Calling this function more than once is possible, but not recommended: + * The function destroys and rebuilds the p4est forest each time it is called. + */ + void + add_periodicity + (const std::vector>&); + + private: /** * MPI communicator to be @@ -759,6 +782,11 @@ namespace parallel * triangulation. */ typename dealii::internal::p4est::types::forest *parallel_forest; + /** + * A data structure that holds some + * information about the ghost cells of the triangulation. + */ + typename dealii::internal::p4est::types::ghost *parallel_ghost; /** * A flag that indicates @@ -839,8 +867,8 @@ namespace parallel * by p4est is located on geometrically * close coarse grid cells. */ - std::vector coarse_cell_to_p4est_tree_permutation; - std::vector p4est_tree_to_coarse_cell_permutation; + std::vector coarse_cell_to_p4est_tree_permutation; + std::vector p4est_tree_to_coarse_cell_permutation; /** * Return a pointer to the p4est @@ -892,6 +920,15 @@ namespace parallel */ void attach_mesh_data(); + /** + * fills a map that, for each vertex, lists all the processors whose + * subdomains are adjacent to that vertex. Used by + * DoFHandler::Policy::ParallelDistributed. + */ + void + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors); template friend class dealii::internal::DoFHandler::Policy::ParallelDistributed; }; @@ -990,6 +1027,15 @@ namespace parallel */ Settings settings; + /** + * Like above, this method, which is only implemented for dim = 2 or 3, + * needs a stub because it is used in dof_handler_policy.cc + */ + void + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors); + }; } } diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index b7e1d21a7c..def3858884 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -1078,9 +1078,10 @@ namespace DoFTools * More information on the topic can be found in the * @ref GlossFaceOrientation "glossary" article. * - * @note This function will not work for DoFHandler objects that are - * built on a parallel::distributed::Triangulation object unless both - * faces (or their children) are owned by the current processor. + * @note For DoFHandler objects that are built on a + * parallel::distributed::Triangulation object + * parallel::distributed::Triangulation::add_periodicity has to be called + * before. * * @author Matthias Maier, 2012 */ @@ -1131,8 +1132,10 @@ namespace DoFTools * function will be used for which the respective flag was set in the * component mask. * - * @note This function will not work for DoFHandler objects that are - * built on a parallel::distributed::Triangulation object. + @ note For DoFHandler objects that are built on a + * parallel::distributed::Triangulation object + * parallel::distributed::Triangulation::add_periodicity has to be called + * before. * * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" * @@ -1157,8 +1160,10 @@ namespace DoFTools * @p orthogonal_equality. This can be used to implement conditions such * as $u(0,y)=u(1,y+1)$. * - * @note This function will not work for DoFHandler objects that are - * built on a parallel::distributed::Triangulation object. + * @note For DoFHandler objects that are built on a + * parallel::distributed::Triangulation object + * parallel::distributed::Triangulation::add_periodicity has to be called + * before. * * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" * @@ -1190,8 +1195,10 @@ namespace DoFTools * meshes with cells not in @ref GlossFaceOrientation * "standard orientation". * - * @note This function will not work for DoFHandler objects that are - * built on a parallel::distributed::Triangulation object. + * @note For DoFHandler objects that are built on a + * parallel::distributed::Triangulation object + * parallel::distributed::Triangulation::add_periodicity has to be called + * before. * * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ @@ -1217,8 +1224,10 @@ namespace DoFTools * meshes with cells not in @ref GlossFaceOrientation * "standard orientation". * - * @note This function will not work for DoFHandler objects that are - * built on a parallel::distributed::Triangulation object. + * @note For DoFHandler objects that are built on a + * parallel::distributed::Triangulation object + * parallel::distributed::Triangulation::add_periodicity has to be called + * before. * * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc index 4269403891..e4d7bfd274 100644 --- a/deal.II/source/distributed/tria.cc +++ b/deal.II/source/distributed/tria.cc @@ -32,12 +32,14 @@ # include # include # include +# include # include # include # include # include # include +# include #endif #include @@ -104,6 +106,10 @@ namespace internal int (&quadrant_is_ancestor) (const types<2>::quadrant *q1, const types<2>::quadrant *q2); + static + int (&quadrant_ancestor_id) (const types<2>::quadrant *q, + int level); + static int (&comm_find_owner) (types<2>::forest *p4est, const types<2>::locidx which_tree, @@ -116,6 +122,16 @@ namespace internal types<2>::topidx num_corners, types<2>::topidx num_vtt); + static + void (&connectivity_join_faces) (types<2>::connectivity *conn, + types<2>::topidx tree_left, + types<2>::topidx tree_right, + int face_left, + int face_right, + int orientation); + + + static void (&connectivity_destroy) (p4est_connectivity_t *connectivity); @@ -171,6 +187,9 @@ namespace internal void (&connectivity_save) (const char *filename, types<2>::connectivity *connectivity); + static + int (&connectivity_is_valid) (types<2>::connectivity *connectivity); + static types<2>::connectivity *(&connectivity_load) (const char *filename, long *length); @@ -233,6 +252,10 @@ namespace internal const types<2>::quadrant *q2) = p4est_quadrant_is_ancestor; + int (&functions<2>::quadrant_ancestor_id) (const types<2>::quadrant *q, + int level) + = p4est_quadrant_ancestor_id; + int (&functions<2>::comm_find_owner) (types<2>::forest *p4est, const types<2>::locidx which_tree, const types<2>::quadrant *q, @@ -245,6 +268,16 @@ namespace internal types<2>::topidx num_vtt) = p4est_connectivity_new; +#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1) + void (&functions<2>::connectivity_join_faces) (types<2>::connectivity *conn, + types<2>::topidx tree_left, + types<2>::topidx tree_right, + int face_left, + int face_right, + int orientation) + = p4est_connectivity_join_faces; +#endif + void (&functions<2>::connectivity_destroy) (p4est_connectivity_t *connectivity) = p4est_connectivity_destroy; @@ -301,6 +334,10 @@ namespace internal types<2>::connectivity *connectivity) = p4est_connectivity_save; + int (&functions<2>::connectivity_is_valid) (types<2>::connectivity + *connectivity) + = p4est_connectivity_is_valid; + types<2>::connectivity * (&functions<2>::connectivity_load) (const char *filename, long *length) @@ -365,6 +402,10 @@ namespace internal int (&quadrant_is_ancestor) (const types<3>::quadrant *q1, const types<3>::quadrant *q2); + static + int (&quadrant_ancestor_id) (const types<3>::quadrant *q, + int level); + static int (&comm_find_owner) (types<3>::forest *p4est, const types<3>::locidx which_tree, @@ -379,6 +420,14 @@ namespace internal types<3>::topidx num_corners, types<3>::topidx num_ctt); + static + void (&connectivity_join_faces) (types<3>::connectivity *conn, + types<3>::topidx tree_left, + types<3>::topidx tree_right, + int face_left, + int face_right, + int orientation); + static void (&connectivity_destroy) (p8est_connectivity_t *connectivity); @@ -434,6 +483,9 @@ namespace internal void (&connectivity_save) (const char *filename, types<3>::connectivity *connectivity); + static + int (&connectivity_is_valid) (types<3>::connectivity *connectivity); + static types<3>::connectivity *(&connectivity_load) (const char *filename, long *length); @@ -497,6 +549,10 @@ namespace internal const types<3>::quadrant *q2) = p8est_quadrant_is_ancestor; + int (&functions<3>::quadrant_ancestor_id) (const types<3>::quadrant *q, + int level) + = p8est_quadrant_ancestor_id; + int (&functions<3>::comm_find_owner) (types<3>::forest *p4est, const types<3>::locidx which_tree, const types<3>::quadrant *q, @@ -514,6 +570,16 @@ namespace internal void (&functions<3>::connectivity_destroy) (p8est_connectivity_t *connectivity) = p8est_connectivity_destroy; +#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1) + void (&functions<3>::connectivity_join_faces) (types<3>::connectivity *conn, + types<3>::topidx tree_left, + types<3>::topidx tree_right, + int face_left, + int face_right, + int orientation) + = p8est_connectivity_join_faces; +#endif + types<3>::forest *(&functions<3>::new_forest) (MPI_Comm mpicomm, types<3>::connectivity *connectivity, types<3>::locidx min_quadrants, @@ -567,6 +633,10 @@ namespace internal types<3>::connectivity *connectivity) = p8est_connectivity_save; + int (&functions<3>::connectivity_is_valid) (types<3>::connectivity + *connectivity) + = p8est_connectivity_is_valid; + types<3>::connectivity * (&functions<3>::connectivity_load) (const char *filename, long *length) @@ -668,6 +738,37 @@ namespace internal { return functions::quadrant_is_ancestor(&q1, &q2); } + + /** + * This struct templatizes the p4est iterate structs and function + * prototypes, which are used to execute callback functions for faces, + * edges, and corners that require local neighborhood information, i.e. + * the neighboring cells */ + template struct iter; + + template <> struct iter<2> + { + typedef p4est_iter_corner_info_t corner_info; + typedef p4est_iter_corner_side_t corner_side; + typedef p4est_iter_corner_t corner_iter; + typedef p4est_iter_face_info_t face_info; + typedef p4est_iter_face_side_t face_side; + typedef p4est_iter_face_t face_iter; + }; + + template <> struct iter<3> + { + typedef p8est_iter_corner_info_t corner_info; + typedef p8est_iter_corner_side_t corner_side; + typedef p8est_iter_corner_t corner_iter; + typedef p8est_iter_edge_info_t edge_info; + typedef p8est_iter_edge_side_t edge_side; + typedef p8est_iter_edge_t edge_iter; + typedef p8est_iter_face_info_t face_info; + typedef p8est_iter_face_side_t face_side; + typedef p8est_iter_face_t face_iter; + }; + } } @@ -736,7 +837,7 @@ namespace const std::vector::active_cell_iterator,unsigned int> > > & vertex_to_cell, - const std::vector &coarse_cell_to_p4est_tree_permutation, + const std::vector &coarse_cell_to_p4est_tree_permutation, const bool set_vertex_info, typename internal::p4est::types::connectivity *connectivity) { @@ -1072,17 +1173,24 @@ namespace const typename internal::p4est::types::forest &forest, const types::subdomain_id my_subdomain) { + ssize_t sidx; // check if this cell exists in // the local p4est cell - if (sc_array_bsearch(const_cast(&tree.quadrants), + if ((sidx = sc_array_bsearch(const_cast(&tree.quadrants), &p4est_cell, - internal::p4est::functions::quadrant_compare) + internal::p4est::functions::quadrant_compare)) != -1) { // yes, cell found in local part of p4est delete_all_children (dealii_cell); if (!dealii_cell->has_children()) dealii_cell->set_subdomain_id(my_subdomain); + + typename internal::p4est::types::quadrant *match = + static_cast::quadrant *> + (sc_array_index_ssize_t(const_cast(&tree.quadrants), + sidx)); + match->p.user_int = dealii_cell->index(); } else { @@ -1167,66 +1275,36 @@ namespace template void - match_quadrant_recursively (const typename internal::p4est::types::quadrant &this_quadrant, - const typename Triangulation::cell_iterator &dealii_cell, - const typename internal::p4est::types::quadrant &ghost_quadrant, - const typename internal::p4est::types::forest &forest, - unsigned int ghost_owner) + match_quadrant (const dealii::Triangulation *tria, + unsigned int dealii_index, + typename internal::p4est::types::quadrant &ghost_quadrant, + unsigned int ghost_owner) { - if (internal::p4est::functions::quadrant_is_equal(&this_quadrant, &ghost_quadrant)) - { - // this is the ghostcell + int i, child_id; + int l = ghost_quadrant.level; - if (dealii_cell->has_children()) - delete_all_children (dealii_cell); - else + for (i = 0; i < l; i++) + { + typename Triangulation::cell_iterator cell (tria, i, dealii_index); + if (cell->has_children () == false) { - dealii_cell->clear_coarsen_flag(); - dealii_cell->set_subdomain_id(ghost_owner); + cell->clear_coarsen_flag(); + cell->set_refine_flag (); + return; } - return; - } - if (! internal::p4est::functions::quadrant_is_ancestor ( &this_quadrant, &ghost_quadrant)) - { - return; - } - - if (dealii_cell->has_children () == false) - { - dealii_cell->clear_coarsen_flag(); - dealii_cell->set_refine_flag (); - return; + child_id = internal::p4est::functions::quadrant_ancestor_id (&ghost_quadrant, i + 1); + dealii_index = cell->child_index(child_id); } - typename internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - for (unsigned int c=0; - c::max_children_per_cell; ++c) - switch (dim) - { - case 2: - P4EST_QUADRANT_INIT(&p4est_child[c]); - break; - case 3: - P8EST_QUADRANT_INIT(&p4est_child[c]); - break; - default: - Assert (false, ExcNotImplemented()); - } - - - internal::p4est::functions:: - quadrant_childrenv (&this_quadrant, p4est_child); - - for (unsigned int c=0; - c::max_children_per_cell; ++c) - + typename Triangulation::cell_iterator cell (tria, l, dealii_index); + if (cell->has_children()) + delete_all_children (cell); + else { - match_quadrant_recursively (p4est_child[c], dealii_cell->child(c), ghost_quadrant, forest, ghost_owner); + cell->clear_coarsen_flag(); + cell->set_subdomain_id(ghost_owner); } - - } @@ -1504,7 +1582,7 @@ namespace { public: RefineAndCoarsenList (const Triangulation &triangulation, - const std::vector &p4est_tree_to_coarse_cell_permutation, + const std::vector &p4est_tree_to_coarse_cell_permutation, const types::subdomain_id my_subdomain, typename internal::p4est::types::forest &forest); @@ -1578,8 +1656,8 @@ namespace template RefineAndCoarsenList:: RefineAndCoarsenList (const Triangulation &triangulation, - const std::vector &p4est_tree_to_coarse_cell_permutation, - const types::subdomain_id my_subdomain, + const std::vector &p4est_tree_to_coarse_cell_permutation, + const types::subdomain_id my_subdomain, typename internal::p4est::types::forest &forest) : forest(forest) @@ -1943,6 +2021,8 @@ namespace parallel // for several different space dimensions dealii::internal::p4est::InitFinalize::do_initialize (); + parallel_ghost = 0; + number_cache.n_locally_owned_active_cells .resize (Utilities::MPI::n_mpi_processes (mpi_communicator)); } @@ -2024,6 +2104,12 @@ namespace parallel { triangulation_has_content = false; + if (parallel_ghost != 0) + { + dealii::internal::p4est::functions::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } + if (parallel_forest != 0) { dealii::internal::p4est::functions::destroy (parallel_forest); @@ -2121,6 +2207,10 @@ namespace parallel Triangulation:: load(const char *filename) { + if (parallel_ghost != 0) { + dealii::internal::p4est::functions::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } dealii::internal::p4est::functions::destroy (parallel_forest); parallel_forest = 0; dealii::internal::p4est::functions::connectivity_destroy (connectivity); @@ -2557,8 +2647,12 @@ namespace parallel // query p4est for the ghost cells - typename dealii::internal::p4est::types::ghost *ghostlayer; - ghostlayer = dealii::internal::p4est::functions::ghost_new (parallel_forest, + if (parallel_ghost != 0) + { + dealii::internal::p4est::functions::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } + parallel_ghost = dealii::internal::p4est::functions::ghost_new (parallel_forest, (dim == 2 ? typename dealii::internal::p4est::types:: @@ -2567,7 +2661,7 @@ namespace parallel typename dealii::internal::p4est::types:: balance_type(P8EST_BALANCE_CORNER))); - Assert (ghostlayer, ExcInternalError()); + Assert (parallel_ghost, ExcInternalError()); // set all cells to artificial. we will later set it to the correct @@ -2632,26 +2726,22 @@ namespace parallel // and recurse. typename dealii::internal::p4est::types::quadrant *quadr; unsigned int ghost_owner=0; + typename dealii::internal::p4est::types::topidx ghost_tree=0; - for (unsigned int g_idx=0; g_idxghosts.elem_count; ++g_idx) + for (unsigned int g_idx=0; g_idxghosts.elem_count; ++g_idx) { - while (g_idx >= (unsigned int)ghostlayer->proc_offsets[ghost_owner+1]) + while (g_idx >= (unsigned int)parallel_ghost->proc_offsets[ghost_owner+1]) ++ghost_owner; + while (g_idx >= (unsigned int)parallel_ghost->tree_offsets[ghost_tree+1]) + ++ghost_tree; quadr = static_cast::quadrant *> - ( sc_array_index(&ghostlayer->ghosts, g_idx) ); + ( sc_array_index(¶llel_ghost->ghosts, g_idx) ); unsigned int coarse_cell_index = - p4est_tree_to_coarse_cell_permutation[quadr->p.piggy3.which_tree]; - - const typename Triangulation::cell_iterator - cell (this, 0U, coarse_cell_index); - - typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; - dealii::internal::p4est::init_coarse_quadrant (p4est_coarse_cell); + p4est_tree_to_coarse_cell_permutation[ghost_tree]; - match_quadrant_recursively (p4est_coarse_cell, cell, *quadr, - *parallel_forest, ghost_owner); + match_quadrant (this, coarse_cell_index, *quadr, ghost_owner); } // fix all the flags to make @@ -2710,7 +2800,7 @@ namespace parallel ++num_ghosts; } - Assert( num_ghosts == ghostlayer->ghosts.elem_count, ExcInternalError()); + Assert( num_ghosts == parallel_ghost->ghosts.elem_count, ExcInternalError()); #endif @@ -2866,8 +2956,6 @@ namespace parallel } - dealii::internal::p4est::functions::ghost_destroy (ghostlayer); - this->smooth_grid = save_smooth; } @@ -2950,6 +3038,11 @@ namespace parallel ExcInternalError()); parallel_forest->user_pointer = &refine_and_coarsen_list; + if (parallel_ghost != 0) + { + dealii::internal::p4est::functions::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } dealii::internal::p4est::functions:: refine (parallel_forest, /* refine_recursive */ false, &RefineAndCoarsenList::refine_callback, @@ -3216,7 +3309,307 @@ namespace parallel return p4est_tree_to_coarse_cell_permutation; } + namespace + { + /** + * This is the callback data structure used to fill + * vertices_with_ghost_neighbors via the p4est_iterate tool + */ + template + struct find_ghosts + { + 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); + struct find_ghosts *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); + struct find_ghosts *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); + struct find_ghosts *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; + } + } + + template <> + void + Triangulation<1,1>:: + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors) + { + Assert (false, ExcNotImplemented()); + } + + template <> + void + Triangulation<1,2>:: + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors) + { + Assert (false, ExcNotImplemented()); + } + + template <> + void + Triangulation<1,3>:: + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors) + { + Assert (false, ExcNotImplemented()); + } + + /** + * Determine the neighboring subdomains that are adjacent to each vertex. + * This is achieved via the p4est_iterate tool + */ + template <> + void + Triangulation<2,2>:: + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors) + { + struct find_ghosts<2,2> fg; + + fg.subids = sc_array_new (sizeof (dealii::types::subdomain_id)); + fg.triangulation = this; + fg.vertices_with_ghost_neighbors = &(vertices_with_ghost_neighbors); + + p4est_iterate (this->parallel_forest, this->parallel_ghost, static_cast(&fg), + NULL, find_ghosts_face<2,2>, find_ghosts_corner<2,2>); + + sc_array_destroy (fg.subids); + } + + /** + * Determine the neighboring subdomains that are adjacent to each vertex. + * This is achieved via the p4est_iterate tool + */ + template <> + void + Triangulation<2,3>:: + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors) + { + struct find_ghosts<2,3> fg; + + fg.subids = sc_array_new (sizeof (dealii::types::subdomain_id)); + fg.triangulation = this; + fg.vertices_with_ghost_neighbors = &(vertices_with_ghost_neighbors); + + p4est_iterate (this->parallel_forest, this->parallel_ghost, static_cast(&fg), + NULL, find_ghosts_face<2,3>, find_ghosts_corner<2,3>); + + sc_array_destroy (fg.subids); + } + + /** + * Determine the neighboring subdomains that are adjacent to each vertex. + * This is achieved via the p8est_iterate tool + */ + template <> + void + Triangulation<3,3>:: + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors) + { + struct find_ghosts<3,3> fg; + + fg.subids = sc_array_new (sizeof (dealii::types::subdomain_id)); + fg.triangulation = this; + fg.vertices_with_ghost_neighbors = &(vertices_with_ghost_neighbors); + p8est_iterate (this->parallel_forest, this->parallel_ghost, static_cast(&fg), + NULL, find_ghosts_face<3,3>, find_ghosts_edge<3,3>, find_ghosts_corner<3,3>); + + sc_array_destroy (fg.subids); + } template MPI_Comm @@ -3226,6 +3619,89 @@ namespace parallel } + template + void + Triangulation::add_periodicity + (const std::vector>& + periodicity_vector) + { +#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1) + Assert (triangulation_has_content == true, + ExcMessage ("The triangulation is empty!")); + Assert (this->n_levels() == 1, + ExcMessage ("The triangulation is refined!")); + + typename std::vector + >::const_iterator periodic_tuple; + periodic_tuple = periodicity_vector.begin(); + + typename std::vector + >::const_iterator periodic_end; + periodic_end = periodicity_vector.end(); + + for (; periodic_tuple(*periodic_tuple); + const cell_iterator second_cell=std::get<2>(*periodic_tuple); + const unsigned int face_right=std::get<3>(*periodic_tuple); + const unsigned int face_left=std::get<1>(*periodic_tuple); + + //respective cells of the matching faces in p4est + const unsigned int tree_left + = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(), + first_cell)]; + const unsigned int tree_right + = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(), + second_cell)]; + + dealii::internal::p4est::functions:: + connectivity_join_faces (connectivity, + tree_left, + tree_right, + face_left, + face_right, + /* orientation */ 0); + /* The orientation parameter above describes the difference between + * the cell on the left and the cell on the right would number of the + * corners of the face. In the periodic domains most users will want, + * this orientation will be 0, i.e. the two cells would number the + * corners the same way. More exotic periodicity, like moebius strips + * or converting an unstructured quad/hex mesh into a periodic domain, + * are not supported right now, and undefined behavior will occur if + * users try to make them periodic. This may be addressed at a later + * date. + */ + } + + + Assert(dealii::internal::p4est::functions::connectivity_is_valid + (connectivity) == 1, + ExcInternalError()); + + // now create a forest out of the + // connectivity data structure + dealii::internal::p4est::functions::destroy (parallel_forest); + parallel_forest + = dealii::internal::p4est::functions:: + new_forest (mpi_communicator, + connectivity, + /* minimum initial number of quadrants per tree */ 0, + /* minimum level of upfront refinement */ 0, + /* use uniform upfront refinement */ 1, + /* user_data_size = */ 0, + /* user_data_constructor = */ NULL, + /* user_pointer */ this); + +#else + Assert(false, ExcMessage ("Need p4est version >=0.3.4.1!")); +#endif + } + template std::size_t @@ -3385,6 +3861,30 @@ namespace parallel } } + 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; + const std::vector perm = triangulation->get_p4est_tree_to_coarse_cell_permutation (); + unsigned int dealii_index = perm[treeidx]; + + for (i = 0; i < l; i++) { + typename dealii::Triangulation::cell_iterator cell (triangulation, i, dealii_index); + child_id = 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; + } + // TODO: again problems with specialization in diff --git a/deal.II/source/dofs/dof_handler_policy.cc b/deal.II/source/dofs/dof_handler_policy.cc index 496a9e16f4..4f6f4d0927 100644 --- a/deal.II/source/dofs/dof_handler_policy.cc +++ b/deal.II/source/dofs/dof_handler_policy.cc @@ -2092,25 +2092,14 @@ namespace internal if (!cell->is_artificial()) cell->set_user_flag(); - //mark the vertices we are interested - //in, i.e. belonging to own and marked cells - const std::vector locally_active_vertices - = mark_locally_active_vertices (*tr); - // add each ghostcells' // subdomain to the vertex and // keep track of interesting // neighbors std::map > vertices_with_ghost_neighbors; - for (typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (cell->is_ghost ()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (locally_active_vertices[cell->vertex_index(v)]) - vertices_with_ghost_neighbors[cell->vertex_index(v)] - .insert (cell->subdomain_id()); + + tr->fill_vertices_with_ghost_neighbors (vertices_with_ghost_neighbors); /* Send and receive cells. After this, @@ -2581,10 +2570,6 @@ namespace internal for (cell = dof_handler.begin_active(); cell != endc; ++cell) if (!cell->is_artificial()) cell->set_user_flag(); - //mark the vertices we are interested - //in, i.e. belonging to own and marked cells - const std::vector locally_active_vertices - = mark_locally_active_vertices (*tr); // add each ghostcells' // subdomain to the vertex and @@ -2592,14 +2577,8 @@ namespace internal // neighbors std::map > vertices_with_ghost_neighbors; - for (typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (cell->is_ghost ()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (locally_active_vertices[cell->vertex_index(v)]) - vertices_with_ghost_neighbors[cell->vertex_index(v)] - .insert (cell->subdomain_id()); + + tr->fill_vertices_with_ghost_neighbors (vertices_with_ghost_neighbors); // Send and receive cells. After this, only // the local cells are marked, that received diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc index 6e21770c30..402557bee3 100644 --- a/deal.II/source/dofs/dof_tools_constraints.cc +++ b/deal.II/source/dofs/dof_tools_constraints.cc @@ -1710,6 +1710,29 @@ namespace DoFTools face_1->get_dof_indices(dofs_1, face_1_index); face_2->get_dof_indices(dofs_2, face_2_index); + for (unsigned int i=0; i < dofs_per_face; i++) { + if (dofs_1[i] == numbers::invalid_dof_index || + dofs_2[i] == numbers::invalid_dof_index) { + /* If either of these faces have no indices, stop. This is so + * that there is no attempt to match artificial cells of + * parallel distributed triangulations. + * + * While it seems like we ought to be able to avoid even calling + * set_periodicity_constraints for artificial faces, this + * situation can arise when a face that is being made periodic + * is only partially touched by the local subdomain. + * make_periodicity_constraints will be called recursively even + * for the section of the face that is not touched by the local + * subdomain. + * + * Until there is a better way to determine if the cells that + * neighbor a face are artificial, we simply test to see if the + * face does not have a valid dof initialization. + */ + return; + } + } + // Well, this is a hack: // // There is no @@ -1983,18 +2006,6 @@ namespace DoFTools Assert (0<=direction && direction PTRIA; - const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); - Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), - ExcMessage ("This function can not be used with distributed triangulations." - "See the documentation for more information.")); - } -#endif - Assert (b_id1 != b_id2, ExcMessage ("The boundary indicators b_id1 and b_id2 must be" "different to denote different boundaries.")); @@ -2077,18 +2088,6 @@ namespace DoFTools Assert(dim == space_dim, ExcNotImplemented()); -#if defined(DEBUG) && defined(DEAL_II_WITH_P4EST) - // Check whether we run on a non parallel mesh or on a - // parallel::distributed::Triangulation in serial - { - typedef typename parallel::distributed::Triangulation PTRIA; - const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); - Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), - ExcMessage ("This function can not be used with distributed triangulations." - "See the documentation for more information.")); - } -#endif - typedef typename DH::face_iterator FaceIterator; typedef std::map FaceMap; -- 2.39.5