From: maier Date: Mon, 21 May 2012 22:21:25 +0000 (+0000) Subject: Add a helper function GridTools::collect_periodic_cell_pairs for periodic X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f53709e9ff6be40678b002a0c1d6fba3d3fa19d4;p=dealii-svn.git Add a helper function GridTools::collect_periodic_cell_pairs for periodic boundary conditions git-svn-id: https://svn.dealii.org/trunk@25529 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index 9358faf344..2bd4b431ef 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -901,6 +901,74 @@ namespace GridTools const std::set &boundary_ids = std::set()); + + /** + * This function loops over all cells + * from @p begin to @p end and collects a + * set of periodic cell pairs for + * @p direction: + * + * Given a @p direction, + * define a 'left' boundary as all + * boundary faces belonging to + * @p boundary_component with corresponding + * unit normal (of the @ref + * GlossReferenceCell "reference cell") in + * negative @p direction, i.e. all boundary + * faces with + * face(2*direction)->at_boundary()==true. + * Analogously, a 'right' boundary + * consisting of all faces of @p + * boundary_component with unit normal + * in positive @p direction, + * face(2*direction+1)->at_boundary()==true. + * + * This function tries to match + * all cells with faces belonging to the + * left' boundary with the cells with faces + * belonging to the 'right' boundary + * with the help of an + * orthogonal equality relation: + * Two cells do match if the vertices of + * the respective boundary faces + * can be transformed into each other by + * parallel translation in @p direction. + * + * The @p offset is a vector tangential to + * the faces that is added to the location + * of vertices of the 'left' boundary when + * attempting to match them to the + * corresponding vertices of the 'right' + * boundary. This can be used to implement + * conditions such as $u(0,y)=u(1,y+1)$. + */ + template + std::map + collect_periodic_cell_pairs (const CellIterator &begin, + const typename identity::type &end, + const types::boundary_id_t boundary_component, + int direction, + const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> + &offset = dealii::Tensor<1,CellIterator::AccessorType::space_dimension>()); + + + /** + * Same as above but matches all active + * cells, i.e. this function calls the + * above function with + * @p dof_handler.begin_active() and + * @p dof_handler.end() as the first two + * arguments. + */ + template + std::map + collect_periodic_cell_pairs (const DH &dof_handler, + const types::boundary_id_t boundary_component, + int direction, + const dealii::Tensor<1,DH::space_dimension> + &offset = Tensor<1,DH::space_dimension>()); + + /** * Exception */ diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index 8c6cea3de2..25ffb101e4 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -2009,7 +2009,139 @@ namespace GridTools return surface_to_volume_mapping; } -} + // Internally used in + // collect_periodic_cell_pairs. + template + inline bool orthogonal_equality (const CellIterator &cell1, + const CellIterator &cell2, + const int direction, + const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> + &offset) + { + // An orthogonal equality test for + // cells: + // Two cells are equal if the + // corresponding vertices of the + // boundary faces can be transformed + // into each other parallel to the + // given direction. + // It is assumed that + // cell1->face(2*direction) and + // cell2->face(2*direction+1) are the + // corresponding boundary faces. + // The optional argument offset will + // be added to each vertex of cell1 + + static const int dim = CellIterator::AccessorType::dimension; + static const int space_dim = CellIterator::AccessorType::space_dimension; + Assert(dim == space_dim, + ExcNotImplemented()); + // ... otherwise we would need two + // directions: One for selecting the + // faces, thus living in dim, + // the other direction for selecting the + // direction for the orthogonal + // projection, thus living in + // space_dim. + + for (int i = 0; i < space_dim; ++i) { + // Only compare + // coordinate-components != direction: + if (i == direction) + continue; + + for (unsigned int v = 0; v < GeometryInfo::vertices_per_face; ++v) + if (fabs ( cell1->face(2*direction)->vertex(v)(i) + + offset[i] + - cell2->face(2*direction+1)->vertex(v)(i)) + > 1.e-10) + return false; + } + return true; + } + + + template + std::map + collect_periodic_cell_pairs (const CellIterator &begin, + const typename identity::type &end, + const types::boundary_id_t boundary_component, + int direction, + const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> + &offset = dealii::Tensor<1,CellIterator::AccessorType::space_dimension>()) + { + static const int space_dim = CellIterator::AccessorType::space_dimension; + Assert (0<=direction && direction cells1; + std::set cells2; + + // collect boundary cells, i.e. cells + // with face(2*direction), resp. + // face(2*direction+1) being on the + // boundary and having the correct + // color: + CellIterator cell = begin; + for (; cell!= end; ++cell) { + if (cell->face(2*direction)->at_boundary() && + cell->face(2*direction)->boundary_indicator() == boundary_component) { + cells1.insert(cell); + } + if (cell->face(2*direction+1)->at_boundary() && + cell->face(2*direction+1)->boundary_indicator() == boundary_component) { + cells2.insert(cell); + } + } + + Assert (cells1.size() == cells2.size(), + ExcMessage ("Unmatched cells on periodic boundaries")); + + std::map matched_cells; + + // Match with a complexity of O(n^2). + // This could be improved... + typedef typename std::set::const_iterator SetIterator; + for (SetIterator it1 = cells1.begin(); it1 != cells1.end(); ++it1) { + for (SetIterator it2 = cells2.begin(); it2 != cells2.end(); ++it2) { + if (orthogonal_equality(*it1, *it2, direction, offset)) { + // We have a match, so insert the + // matching pairs and remove the matched + // cell in cells2 to speed up the + // matching: + matched_cells[*it1] = *it2; + cells2.erase(it2); + break; + } + } + } + + Assert (matched_cells.size() == cells1.size() && + cells2.size() == 0, + ExcMessage ("Unmatched cells on periodic boundaries")); + + return matched_cells; + } + + + template + std::map + collect_periodic_cell_pairs (const DH &dof_handler, + const types::boundary_id_t boundary_component, + int direction, + const dealii::Tensor<1,DH::space_dimension> + &offset = Tensor<1,DH::space_dimension>()) + { + return collect_periodic_cell_pairs + (dof_handler.begin_active(), + dof_handler.end(), + boundary_component, + direction, + offset); + } + + +} /* namespace GridTools */ // explicit instantiations diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index e9da14bd40..63a370af00 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -176,3 +176,30 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS #endif } + + +for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) +{ +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools \{ + + template + std::map + collect_periodic_cell_pairs(const X::cell_iterator &, + const X::cell_iterator &, + const types::boundary_id_t, + int, + const Tensor<1,deal_II_space_dimension> &); + + template + std::map + collect_periodic_cell_pairs (const X&, + const types::boundary_id_t, + int, + const Tensor<1,deal_II_space_dimension> &); + + + \} + #endif +} +