}
}
-
+//forward declaration of the data type for periodic face pairs
+namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
namespace parallel
{
get_p4est_tree_to_coarse_cell_permutation() const;
+
/**
* Join faces in the p4est forest due to periodic boundary conditions.
*
+ * The vector can be filled by the function
+ * GridTools::collect_periodic_faces.
+ *
+ * @todo At the moment just default orientation is implemented.
+ *
+ * @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<GridTools::PeriodicFacePair<cell_iterator> >&);
+
+ /**
+ * Same as the function above, but takes a different argument.
+ *
* The entries in the std::vector should have the form
* std_cxx1x::tuple<cell1, face_no1, cell2, face_no2>.
*
* The vector can be filled by the function
- * DoFTools::identify_periodic_face_pairs.
- *
+ * GridTools::identify_periodic_face_pairs.
+ *
+ * @note This function can only be used if the faces are in
+ * default orientation.
+ *
* @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:
cell_iterator, unsigned int> >&);
+
private:
/**
* MPI communicator to be
template <class GridClass> class InterGridMap;
template <int dim, int spacedim> class Mapping;
+namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
//TODO: map_support_points_to_dofs should generate a multimap, rather than just a map, since several dofs may be located at the same support point
*/
template<typename FaceIterator>
void
- make_periodicity_constraints (const FaceIterator &face_1,
- const typename identity<FaceIterator>::type &face_2,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask(),
- const bool face_orientation = true,
- const bool face_flip = false,
- const bool face_rotation = false);
+ make_periodicity_constraints
+ (const FaceIterator &face_1,
+ const typename identity<FaceIterator>::type &face_2,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask(),
+ const bool face_orientation = true,
+ const bool face_flip = false,
+ const bool face_rotation = false);
+
/**
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
/**
*
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*
- * @author Matthias Maier, 2012
+ * @author Daniel Arndt, 2013, Matthias Maier, 2012
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- dealii::Tensor<1,DH::space_dimension> &offset,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
+ dealii::Tensor<1,DH::space_dimension> &offset,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
/**
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id,
- const int direction,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const int direction,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
/**
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id,
- const int direction,
- dealii::Tensor<1,DH::space_dimension> &offset,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const int direction,
+ dealii::Tensor<1,DH::space_dimension> &offset,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
+
+ /**
+ * Same as above but the periodicity information is given by
+ * @p periodic_faces. This std::vector can be created by
+ * GridTools::collect_periodic_faces.
+ *
+ * @note For DoFHandler objects that are built on a
+ * parallel::distributed::Triangulation object
+ * parallel::distributed::Triangulation::add_periodicity has to be called
+ * before.
+ */
+ template<typename DH>
+ void
+ make_periodicity_constraints
+ (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+ &periodic_faces,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
//@}
const std::set<types::boundary_id> &boundary_ids
= std::set<types::boundary_id>());
-
+ /**
+ * Data type that provides all the information that is needed
+ * to create periodicity constraints and a periodic p4est forest
+ * with respect to two periodic cell faces
+ */
+ template<typename CellIterator>
+ struct PeriodicFacePair
+ {
+ CellIterator cell[2];
+ unsigned int face_idx[2];
+ std::bitset<3> orientation;
+ };
/**
* An orthogonal equality test for faces.
/**
- * This function loops over all faces from @p begin to @p end and
- * collects a set of periodic face pairs for @p direction:
+ * This function will collect periodic face pairs on the highest (i.e.
+ * coarsest) mesh level.
*
* Define a 'first' boundary as all boundary faces having boundary_id
* @p b_id1 and a 'second' boundary consisting of all faces belonging
* boundary with faces belonging to the second boundary with the help
* of orthogonal_equality().
*
- * The bitset that is returned together with the second face encodes the
+ * The bitset that is returned inside of PeriodicFacePair encodes the
* _relative_ orientation of the first face with respect to the second
* face, see the documentation of orthogonal_equality for further details.
*
+ * The @p direction refers to the space direction in which periodicity
+ * is enforced.
+ *
* The @p offset is a vector tangential to the faces that is added to the
* location of vertices of the 'first' boundary when attempting to match
* them to the corresponding vertices of the 'second' boundary. This can
* be used to implement conditions such as $u(0,y)=u(1,y+1)$.
*
- * @author Matthias Maier, 2012
- */
- template<typename FaceIterator>
- std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
- collect_periodic_face_pairs (const FaceIterator &begin,
- const typename identity<FaceIterator>::type &end,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset);
-
-
- /**
- * Same function as above, but accepts a Triangulation or DoFHandler
- * object @p container (a container is a collection of objects, here a
- * collection of cells) instead of an explicit face iterator range.
- *
- * This function will collect periodic face pairs on the highest (i.e.
- * coarsest) mesh level.
+ * @note The created std::vector can be used in
+ * DoFTools::make_periodicity_constraints and in
+ * parallel::distributed::Triangulation::add_periodicity to enforce
+ * periodicity algebraically.
*
- * @author Matthias Maier, 2012
+ * @author Daniel Arndt, 2013
*/
template<typename DH>
- std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
- collect_periodic_face_pairs (const DH &container,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- const dealii::Tensor<1,DH::space_dimension> &offset);
-
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset);
/**
* This compatibility version of collect_periodic_face_pairs only works
* meshes with cells not in @ref GlossFaceOrientation
* "standard orientation".
*
+ * @author Daniel Arndt, 2013
+ */
+ template<typename DH>
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces
+ (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset);
+
+ /**
+ * This function does the same as collect_periodic_faces but returns a
+ * different data type.
+ *
+ * @author Matthias Maier, 2012
+ *
+ * @note The returned data type is not compatible with
+ * DoFTools::make_periodicity_constraints
+ *
+ * @deprecated
+ */
+ template<typename DH>
+ std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
+ collect_periodic_face_pairs
+ (const DH &container,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED;
+
+ /**
+ * This compatibility version of collect_periodic_face_pairs only works
+ * on grids with cells in @ref GlossFaceOrientation "standard orientation".
+ *
* @author Matthias Maier, 2012
+ *
+ * @note The returned data type is not compatible with
+ * DoFTools::make_periodicity_constraints
+ *
+ * @deprecated
*/
template<typename DH>
std::map<typename DH::face_iterator, typename DH::face_iterator>
- collect_periodic_face_pairs (const DH &dof_handler, /*TODO: Name*/
- const types::boundary_id b_id,
- const int direction,
- const dealii::Tensor<1,DH::space_dimension> &offset);
+ collect_periodic_face_pairs
+ (const DH &dof_handler, /*TODO: Name*/
+ const types::boundary_id b_id,
+ int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED;
/**
- * Add periodicity information to the @p periodicity_vector for the
- * boundaries with boundary id @p b_id1 and @p b_id2 in cartesian
- * direction @p direction.
+ * This version loops over all faces from @p begin to @p end
+ * instead of accepting a DoFHandler or a Triangulation.
*
- * This function tries to match all faces belonging to the first
- * boundary with faces belonging to the second boundary
- * by comparing the center of the cell faces. To find the correct
- * corresponding faces, the direction argument indicates in which
- * cartesian direction periodicity should be set.
- * ((0,1,2) -> (x,y,z)-direction)
+ * @author Matthias Maier, 2012
+ *
+ * @note This function cannot produce the return as the other
+ * collect_periodic_faces functions.
+ *
+ * @deprecated
+ */
+ template<typename FaceIterator>
+ std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
+ collect_periodic_face_pairs
+ (const FaceIterator &begin,
+ const typename identity<FaceIterator>::type &end,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ DEAL_II_DEPRECATED;
+
+ /**
+ * This function does the same as collect_periodic_faces but returns a
+ * different data type.
+ *
+ * @author Daniel Arndt, 2013
*
- * The output of this function can be used in
+ * @note The returned data type is not compatible with
+ * DoFTools::make_periodicity_constraints, but with a version of
* parallel::distributed::Triangulation::add_periodicity
+ *
+ * @note Use collect_periodic_faces instead.
+ *
+ * @deprecated
*/
template<typename DH>
void
identify_periodic_face_pairs
- (const DH &dof_handler,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const unsigned int direction,
- std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
- typename DH::cell_iterator, unsigned int> >
- &periodicity_vector);
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+ typename DH::cell_iterator, unsigned int> >
+ &periodicity_vector) DEAL_II_DEPRECATED;
/**
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/distributed/tria.h>
#include <deal.II/grid/grid_tools.h>
+#include <deal.II/distributed/tria.h>
#ifdef DEAL_II_WITH_P4EST
# include <p4est_bits.h>
return mpi_communicator;
}
-
- template <int dim, int spacedim>
+ template<int dim, int spacedim>
void
Triangulation<dim,spacedim>::add_periodicity
- (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
- cell_iterator, unsigned int> >&
- periodicity_vector)
+ (const std::vector<GridTools::PeriodicFacePair<cell_iterator> >&
+ periodicity_vector)
{
-#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+ #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
- <typename std_cxx1x::tuple
- <cell_iterator, unsigned int,
- cell_iterator, unsigned int> >::const_iterator periodic_tuple;
- periodic_tuple = periodicity_vector.begin();
-
- typename std::vector
- <typename std_cxx1x::tuple
- <cell_iterator, unsigned int,
- cell_iterator, unsigned int> >::const_iterator periodic_end;
- periodic_end = periodicity_vector.end();
-
- for (; periodic_tuple<periodic_end; ++periodic_tuple)
+
+ typedef std::vector<GridTools::PeriodicFacePair<cell_iterator> >
+ FaceVector;
+ typename FaceVector::const_iterator it, periodic_end;
+ it = periodicity_vector.begin();
+
+ for (; it<periodic_end; ++it)
{
- const cell_iterator first_cell=std_cxx1x::get<0>(*periodic_tuple);
- const cell_iterator second_cell=std_cxx1x::get<2>(*periodic_tuple);
- const unsigned int face_right=std_cxx1x::get<3>(*periodic_tuple);
- const unsigned int face_left=std_cxx1x::get<1>(*periodic_tuple);
-
+ const cell_iterator first_cell = it->cell[0];
+ const cell_iterator second_cell = it->cell[1];
+ const unsigned int face_left = it->face_idx[0];
+ const unsigned int face_right = it->face_idx[1];
+
//respective cells of the matching faces in p4est
const unsigned int tree_left
= coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
= coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
second_cell)];
+ //TODO Add support for non default orientation.
+ Assert(it->orientation == 1,
+ ExcMessage("Found a face match with non standard orientation. "
+ "This function is only suitable for meshes with "
+ "cells in default orientation"));
+
dealii::internal::p4est::functions<dim>::
connectivity_join_faces (connectivity,
tree_left,
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,
* date.
*/
}
-
-
+
+
Assert(dealii::internal::p4est::functions<dim>::connectivity_is_valid
- (connectivity) == 1,
- ExcInternalError());
-
- // now create a forest out of the connectivity data structure
+ (connectivity) == 1, ExcInternalError());
+
+ // now create a forest out of the connectivity data structure
dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
parallel_forest
= dealii::internal::p4est::functions<dim>::
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);
+ /* 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 <int dim, int spacedim>
+ void
+ Triangulation<dim,spacedim>::add_periodicity
+ (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >&
+ periodicity_vector)
+ {
+#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+ typedef std::vector<std_cxx1x::tuple <cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >
+ FaceVector;
+
+ typename FaceVector::const_iterator it, end_periodic;
+ it = periodicity_vector.begin();
+ end_periodic = periodicity_vector.end();
+
+ std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_faces;
+
+ for(; it!=end_periodic; ++it)
+ {
+ const cell_iterator cell1 = std::get<0> (*it);
+ const cell_iterator cell2 = std::get<2> (*it);
+ const unsigned int face_idx1 = std::get<1> (*it);
+ const unsigned int face_idx2 = std::get<3> (*it);
+ const GridTools::PeriodicFacePair<cell_iterator> matched_face
+ = {{cell1, cell2},{face_idx1, face_idx2}, 1};
+ periodic_faces.push_back(matched_face);
+ }
+ add_periodicity(periodic_faces);
#else
Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
#endif
}
-
+
template<typename DH>
void
make_periodicity_constraints (const DH &dof_handler,
ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
"different to denote different boundaries."));
- typedef typename DH::face_iterator FaceIterator;
- typedef std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > > FaceMap;
+ typedef std::vector<GridTools::PeriodicFacePair
+ <typename DH::cell_iterator> > FaceVector;
// Collect matching periodic cells on the coarsest level:
- FaceMap matched_cells =
- GridTools::collect_periodic_face_pairs(dof_handler,
- b_id1, b_id2,
- direction, offset);
+ FaceVector matched_cells =
+ GridTools::collect_periodic_faces(dof_handler, b_id1, b_id2,
+ direction, offset);
// And apply the low level make_periodicity_constraints function to
// every matching pair:
- for (typename FaceMap::iterator it = matched_cells.begin();
+ for (typename FaceVector::iterator it = matched_cells.begin();
it != matched_cells.end(); ++it)
{
typedef typename DH::face_iterator FaceIterator;
- const FaceIterator &face_1 = it->first;
- const FaceIterator &face_2 = it->second.first;
- const std::bitset<3> &orientation = it->second.second;
+ const FaceIterator &face_1 = it->cell[0]->face(it->face_idx[0]);
+ const FaceIterator &face_2 = it->cell[1]->face(it->face_idx[1]);
+ const std::bitset<3> &orientation = it->orientation;
Assert(face_1->at_boundary() && face_2->at_boundary(),
ExcInternalError());
Assert(dim == space_dim,
ExcNotImplemented());
- typedef typename DH::face_iterator FaceIterator;
- typedef std::map<FaceIterator, FaceIterator> FaceMap;
+ typedef std::vector<GridTools::PeriodicFacePair
+ <typename DH::cell_iterator> > FaceVector;
// Collect matching periodic cells on the coarsest level:
- FaceMap matched_cells =
- GridTools::collect_periodic_face_pairs(dof_handler,
- b_id,
- direction, offset);
+ const FaceVector matched_cells =
+ GridTools::collect_periodic_faces(dof_handler, b_id,
+ direction, offset);
// And apply the low level make_periodicity_constraints function to
// every matching pair:
- for (typename FaceMap::iterator it = matched_cells.begin();
+ for (typename FaceVector::const_iterator it = matched_cells.begin();
it != matched_cells.end(); ++it)
{
typedef typename DH::face_iterator FaceIterator;
- const FaceIterator &face_1 = it->first;
- const FaceIterator &face_2 = it->second;
+ const FaceIterator &face_1 = it->cell[0]->face(it->face_idx[0]);
+ const FaceIterator &face_2 = it->cell[1]->face(it->face_idx[1]);
Assert(face_1->at_boundary() && face_2->at_boundary(),
ExcInternalError());
+ template<typename DH>
+ void
+ make_periodicity_constraints
+ (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+ &periodic_faces,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask)
+ {
+ typedef std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+ FaceVector;
+ typename FaceVector::const_iterator it, end_periodic;
+ it = periodic_faces.begin();
+ end_periodic = periodic_faces.end();
+
+ for(; it!=end_periodic; ++it)
+ {
+ typedef typename DH::face_iterator FaceIterator;
+ const FaceIterator face_1 = it->cell[0]->face(it->face_idx[0]);
+ const FaceIterator face_2 = it->cell[1]->face(it->face_idx[1]);
+
+ make_periodicity_constraints(face_1,
+ face_2,
+ constraint_matrix,
+ component_mask,
+ it->orientation[0],
+ it->orientation[1],
+ it->orientation[2]);
+ }
+ }
+
+
+
namespace internal
{
namespace
dealii::Tensor<1,DH::space_dimension> &,
dealii::ConstraintMatrix &,
const ComponentMask &);
+
+ template
+ void
+ DoFTools::make_periodicity_constraints<DH>
+ (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> > &,
+ dealii::ConstraintMatrix &,
+ const ComponentMask &);
+
#endif
}
}
-
/*
* Internally used in orthogonal_equality
*
*
* See the comment on the next function as well as the detailed
* documentation of make_periodicity_constraints and
- * collect_periodic_face_pairs for details
+ * collect_periodic_faces for details
*/
template<int dim> struct OrientationLookupTable {};
/*
- * Internally used in collect_periodic_face_pairs
+ * Internally used in collect_periodic_faces
+ */
+ template<typename CellIterator>
+ std::vector<PeriodicFacePair<CellIterator> >
+ match_periodic_face_pairs
+ (std::set<std::pair<CellIterator, unsigned int> >
+ &pairs1,
+ std::set<std::pair<typename identity<CellIterator>::type, unsigned int> >
+ &pairs2,
+ int direction,
+ const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
+ &offset)
+ {
+ static const int space_dim = CellIterator::AccessorType::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ Assert (pairs1.size() == pairs2.size(),
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ typename std::vector<PeriodicFacePair<CellIterator> > matched_faces;
+
+ // Match with a complexity of O(n^2). This could be improved...
+ std::bitset<3> orientation;
+ typedef typename std::set
+ <std::pair<CellIterator, unsigned int> >::const_iterator PairIterator;
+ for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
+ {
+ for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
+ {
+ const CellIterator cell1 = it1->first;
+ const CellIterator cell2 = it2->first;
+ const unsigned int face_idx1 = it1->second;
+ const unsigned int face_idx2 = it2->second;
+ if (GridTools::orthogonal_equality(orientation,
+ cell1->face(face_idx1),
+ cell2->face(face_idx2),
+ direction, offset))
+ {
+ // We have a match, so insert the matching pairs and
+ // remove the matched cell in pairs2 to speed up the
+ // matching:
+ const PeriodicFacePair<CellIterator> matched_face
+ = {{cell1, cell2},{face_idx1, face_idx2}, orientation};
+ matched_faces.push_back(matched_face);
+ pairs2.erase(it2);
+ break;
+ }
+ }
+ }
+
+ AssertThrow (matched_faces.size() == pairs1.size() && pairs2.size() == 0,
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ return matched_faces;
+ }
+
+ /* Deprecated version of the function above with different return value.
+ * It is used the deprecated collect_periodic_face_pairs.
*/
template<typename FaceIterator>
std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
- match_periodic_face_pairs (std::set<FaceIterator> &faces1, /* not const! */
- std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
- int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ match_periodic_face_pairs
+ (std::set<FaceIterator> &faces1, /* not const! */
+ std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
+ int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ DEAL_II_DEPRECATED;
+
+ template<typename FaceIterator>
+ std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
+ match_periodic_face_pairs
+ (std::set<FaceIterator> &faces1, /* not const! */
+ std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
+ int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
{
static const int space_dim = FaceIterator::AccessorType::space_dimension;
Assert (0<=direction && direction<space_dim,
}
+ template<typename DH>
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset)
+ {
+ static const unsigned int dim = DH::dimension;
+ static const unsigned int space_dim = DH::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ // Loop over all cells on the highest level and collect all boundary
+ // faces belonging to b_id1 and b_id2:
+
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs1;
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs2;
+
+ for (typename DH::cell_iterator cell = dof_handler.begin(0);
+ cell != dof_handler.end(0); ++cell)
+ {
+ for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+ {
+ const typename DH::face_iterator face = cell->face(i);
+ if (face->at_boundary() && face->boundary_indicator() == b_id1)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair1
+ = std::make_pair(cell, i);
+ pairs1.insert(pair1);
+ }
+
+ if (face->at_boundary() && face->boundary_indicator() == b_id2)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair2
+ = std::make_pair(cell, i);
+ pairs2.insert(pair2);
+ }
+ }
+ }
+
+ Assert (pairs1.size() == pairs2.size(),
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ // and call match_periodic_face_pairs that does the actual matching:
+ return match_periodic_face_pairs(pairs1, pairs2, direction, offset);
+ }
+
+ template<typename DH>
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset)
+ {
+ static const unsigned int dim = DH::dimension;
+ static const unsigned int space_dim = DH::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ Assert(dim == space_dim,
+ ExcNotImplemented());
+
+ // Loop over all cells on the highest level and collect all boundary
+ // faces 2*direction and 2*direction*1:
+
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs1;
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs2;
+
+ for (typename DH::cell_iterator cell = dof_handler.begin(0);
+ cell != dof_handler.end(0); ++cell)
+ {
+ const typename DH::face_iterator face_1 = cell->face(2*direction);
+ const typename DH::face_iterator face_2 = cell->face(2*direction+1);
+
+ if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair1
+ = std::make_pair(cell, 2*direction);
+ pairs1.insert(pair1);
+ }
+
+ if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair1
+ = std::make_pair(cell, 2*direction+1);
+ pairs1.insert(pair1);
+ }
+ }
+
+ Assert (pairs1.size() == pairs2.size(),
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ // and call match_periodic_face_pairs that does the actual matching:
+
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ FaceVector;
+
+ FaceVector matching = match_periodic_face_pairs(pairs1, pairs2,
+ direction, offset);
+
+ for (typename FaceVector::iterator pairing = matching.begin();
+ pairing != matching.end(); ++pairing)
+ {
+ Assert(pairing->orientation == 1,
+ ExcMessage("Found a face match with non standard orientation. "
+ "This function is only suitable for meshes with cells "
+ "in default orientation"));
+ }
+
+ return matching;
+ }
template<typename FaceIterator>
std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
return match_periodic_face_pairs(faces1, faces2, direction, offset);
}
-
-
template<typename DH>
std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
collect_periodic_face_pairs (const DH &dof_handler,
int direction,
const dealii::Tensor<1,DH::space_dimension> &offset)
{
- static const int dim = DH::dimension;
- static const int space_dim = DH::space_dimension;
- Assert (0<=direction && direction<space_dim,
- ExcIndexRange (direction, 0, space_dim));
-
- // Loop over all cells on the highest level and collect all boundary
- // faces belonging to b_id1 and b_id2:
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> > FaceVector;
- std::set<typename DH::face_iterator> faces1;
- std::set<typename DH::face_iterator> faces2;
+ const FaceVector face_vector
+ = collect_periodic_faces (dof_handler, b_id1, b_id2, direction, offset);
- for (typename DH::cell_iterator cell = dof_handler.begin(0);
- cell != dof_handler.end(0); ++cell)
- {
- for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
- {
- const typename DH::face_iterator face = cell->face(i);
-
- if (face->at_boundary() && face->boundary_indicator() == b_id1)
- faces1.insert(face);
-
- if (face->at_boundary() && face->boundary_indicator() == b_id2)
- faces2.insert(face);
- }
- }
-
- Assert (faces1.size() == faces2.size(),
- ExcMessage ("Unmatched faces on periodic boundaries"));
+ std::map<typename DH::face_iterator,
+ std::pair<typename DH::face_iterator, std::bitset<3> > >
+ return_value;
+ for(typename FaceVector::const_iterator it = face_vector.begin();
+ it != face_vector.end(); ++it)
+ {
+ const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]);
+ const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]);
+ return_value[face1] = std::make_pair(face2, it->orientation);
+ }
- // and call match_periodic_face_pairs that does the actual matching:
- return match_periodic_face_pairs(faces1, faces2, direction, offset);
+ return return_value;
}
-
-
+
template<typename DH>
std::map<typename DH::face_iterator, typename DH::face_iterator>
collect_periodic_face_pairs (const DH &dof_handler,
int direction,
const dealii::Tensor<1,DH::space_dimension> &offset)
{
- static const int dim = DH::dimension;
- static const int space_dim = DH::space_dimension;
- Assert (0<=direction && direction<space_dim,
- ExcIndexRange (direction, 0, space_dim));
-
- Assert(dim == space_dim,
- ExcNotImplemented());
-
- // Loop over all cells on the highest level and collect all boundary
- // faces 2*direction and 2*direction*1:
-
- std::set<typename DH::face_iterator> faces1;
- std::set<typename DH::face_iterator> faces2;
-
- for (typename DH::cell_iterator cell = dof_handler.begin(0);
- cell != dof_handler.end(0); ++cell)
- {
- const typename DH::face_iterator face_1 = cell->face(2*direction);
- const typename DH::face_iterator face_2 = cell->face(2*direction+1);
-
- if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
- faces1.insert(face_1);
-
- if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
- faces2.insert(face_2);
- }
-
- Assert (faces1.size() == faces2.size(),
- ExcMessage ("Unmatched faces on periodic boundaries"));
-
- // and call match_periodic_face_pairs that does the actual matching:
-
- typedef std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > > FaceMap;
- FaceMap matching = match_periodic_face_pairs(faces1, faces2, direction, offset);
-
- std::map<typename DH::face_iterator, typename DH::face_iterator>
- return_value;
-
- for (typename FaceMap::iterator pairing = matching.begin();
- pairing != matching.end(); ++pairing)
- {
- Assert(pairing->second.second == 1,
- ExcMessage("Found a face match with non standard orientation. "
- "This function is only suitable for meshes with cells "
- "in default orientation"));
-
- return_value[pairing->first] = pairing->second.first;
- }
-
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> > FaceVector;
+
+ const FaceVector face_vector
+ = collect_periodic_faces (dof_handler, b_id, direction, offset);
+
+ std::map<typename DH::face_iterator, typename DH::face_iterator> return_value;
+ for(typename FaceVector::const_iterator it = face_vector.begin();
+ it != face_vector.end(); ++it)
+ {
+ const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]);
+ const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]);
+ return_value[face1] = face2;
+ }
+
return return_value;
}
-
+
template<typename DH>
void
identify_periodic_face_pairs
typename DH::cell_iterator, unsigned int> >
&periodicity_vector)
{
- static const unsigned int dim = DH::dimension;
- static const unsigned int spacedim = DH::space_dimension;
-
- Assert (0<=direction && direction<spacedim,
- ExcIndexRange(direction, 0, spacedim));
-
- Assert (b_id1 != b_id2,
- ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
- "different to denote different boundaries."));
-
- std::map<std::pair<typename DH::cell_iterator,
- unsigned int>,
- Point<spacedim> > face_locations;
-
- // Collect faces with boundary_indicator b_id1
- typename DH::cell_iterator cell = dof_handler.begin();
- typename DH::cell_iterator endc = dof_handler.end();
- for (; cell != endc; ++cell)
- for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
- if(cell->face(f)->boundary_indicator() == b_id1)
- {
- Point<spacedim> face_center (cell->face(f)->center());
- face_center(direction)=0.;
- const std::pair<typename DH::cell_iterator, unsigned int>
- cell_face_pair = std::make_pair(cell, f);
- face_locations[cell_face_pair]=face_center;
- }
-
- // Match faces with boundary_indicator b_id2 to the ones in
- //face_locations
- cell = dof_handler.begin();
- for (; cell != endc; ++cell)
- for(unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
- if(cell->face(f)->boundary_indicator() == b_id2)
- {
- typename std::map<std::pair<typename DH::cell_iterator,
- unsigned int>,
- Point<spacedim> >::iterator p
- = face_locations.begin();
-
- Point<spacedim> center2 (cell->face(f)->center());
- center2(direction)=0.;
-
- for (; p != face_locations.end(); ++p)
- {
- if (center2.distance(p->second) < 1e-4*cell->face(f)->diameter())
- {
- const std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
- typename DH::cell_iterator, unsigned int>
- periodic_tuple (p->first.first,
- p->first.second,
- cell,
- f);
-
- periodicity_vector.push_back(periodic_tuple);
-
- face_locations.erase(p);
- break;
- }
- Assert (p != face_locations.end(),
- ExcMessage ("No corresponding face was found!"));
- }
- }
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ FaceVector;
+ const FaceVector periodic_faces
+ = collect_periodic_faces(dof_handler, b_id1, b_id2, direction,
+ dealii::Tensor<1,DH::space_dimension> ());
+
+ typename FaceVector::const_iterator it, end_faces;
+ it = periodic_faces.begin();
+ end_faces = periodic_faces.end();
+
+ for(; it!=end_faces; ++it)
+ {
- Assert (face_locations.size() == 0,
- ExcMessage ("There are unmatched faces!"));
+ const std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+ typename DH::cell_iterator, unsigned int>
+ periodic_tuple (it->cell[0], it->face_idx[0],
+ it->cell[1], it->face_idx[1]);
+
+ periodicity_vector.push_back(periodic_tuple);
+ }
}
#if deal_II_dimension >= 2
+ template
+ std::vector<PeriodicFacePair<X::cell_iterator> >
+ collect_periodic_faces
+ (const X &,
+ const types::boundary_id,
+ const types::boundary_id,
+ unsigned int,
+ const Tensor<1,deal_II_space_dimension> &);
+
+ template<typename DH>
+ std::vector<PeriodicFacePair<X::cell_iterator> >
+ collect_periodic_faces
+ (const X &,
+ const types::boundary_id,
+ unsigned int,
+ const Tensor<1,deal_II_space_dimension> &);
+
template
std::map<X::face_iterator, std::pair<X::face_iterator, std::bitset<3> > >
collect_periodic_face_pairs
const types::boundary_id,
int,
const Tensor<1,deal_II_space_dimension> &);
-
template
std::map<X::active_face_iterator, std::pair<X::active_face_iterator, std::bitset<3> > >
collect_periodic_face_pairs
const types::boundary_id,
int,
const Tensor<1,deal_II_space_dimension> &);
-
+
template
void
identify_periodic_face_pairs
projects / dealii-svn.git / commitdiff