#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/mapping_q1.h>
+#include <bitset>
#include <list>
DEAL_II_NAMESPACE_OPEN
* function which computes all
* active vertex patches by looping
* over cells.
- *
+ *
* Find and return a vector of
* iterators to active cells that
* surround a given vertex @p vertex.
= std::set<types::boundary_id>());
+
/**
- * This function loops over all cells
- * from @p begin to @p end and collects a
- * set of periodic cell pairs for
- * @p direction:
+ * An orthogonal equality test for faces.
+ *
+ * @p face1 and @p face2 are considered equal, if a one to one matching
+ * between its vertices can be achieved via an orthogonal equality
+ * relation: Two vertices <tt>v_1</tt> and <tt>v_2</tt> are considered
+ * equal, if <code> (v_1 + offset) - v_2</code> is parallel to the unit
+ * vector in @p direction.
+ *
+ * If the matching was successful, the _relative_ orientation of @p face1
+ * with respect to @p face2 is returned in the bitset @p orientation,
+ * where
+ * @code
+ * orientation[0] -> face_orientation
+ * orientation[1] -> face_flip
+ * orientation[2] -> face_rotation
+ * @endcode
+ *
+ * In 2D <tt>face_orientation</tt> is always <tt>true</tt>,
+ * <tt>face_rotation</tt> is always <tt>false</tt>, and face_flip has the
+ * meaning of <tt>line_flip</tt>. More precisely in 3d:
+ *
+ * <tt>face_orientation</tt>:
+ * <tt>true</tt> if @p face1 and @p face2 have the same orientation.
+ * Otherwise, the vertex indices of @p face1 match the vertex indices of
+ * @p face2 in the following manner:
+ *
+ * @code
+ * face1: face2:
+ *
+ * 1 - 3 2 - 3
+ * | | <--> | |
+ * 0 - 2 0 - 1
+ * @endcode
*
- * 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
- * <tt>face(2*direction)->at_boundary()==true</tt>.
- * Analogously, a 'right' boundary
- * consisting of all faces of @p
- * boundary_component with unit normal
- * in positive @p direction,
- * <tt>face(2*direction+1)->at_boundary()==true</tt>.
+ * <tt>face_flip</tt>:
+ * <tt>true</tt> if the matched vertices are rotated by 180 degrees:
*
- * 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.
+ * @code
+ * face1: face2:
*
- * 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)$.
+ * 1 - 0 2 - 3
+ * | | <--> | |
+ * 3 - 2 0 - 1
+ * @endcode
+ *
+ * <tt>face_rotation</tt>: <tt>true</tt> if the matched vertices are
+ * rotated by 90 degrees counterclockwise:
+ *
+ * @code
+ * face1: face2:
+ *
+ * 0 - 2 2 - 3
+ * | | <--> | |
+ * 1 - 3 0 - 1
+ * @endcode
+ *
+ * and any combination of that...
+ * More information on the topic can be found in the
+ * @ref GlossFaceOrientation "glossary" article.
+ *
+ * @author Matthias Maier, 2012
*/
- template<typename CellIterator>
- std::map<CellIterator, CellIterator>
- collect_periodic_cell_pairs (const CellIterator &begin,
- const typename identity<CellIterator>::type &end,
- const types::boundary_id boundary_component,
+ template<typename FaceIterator>
+ bool
+ orthogonal_equality (std::bitset<3> &orientation,
+ const FaceIterator &face1,
+ const FaceIterator &face2,
+ const int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset);
+
+
+ /**
+ * Same function as above, but doesn't return the actual orientation
+ */
+ template<typename FaceIterator>
+ bool
+ orthogonal_equality (const FaceIterator &face1,
+ const FaceIterator &face2,
+ const int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset);
+
+
+ /**
+ * This function loops over all faces from @p begin to @p end and
+ * collects a set of periodic face pairs for @p direction:
+ *
+ * Define a 'first' boundary as all boundary faces having boundary_id
+ * @p b_id1 and a 'second' boundary consisting of all faces belonging
+ * to @p b_id2.
+ *
+ * This function tries to match all faces belonging to the first
+ * boundary with faces belonging to the second boundary with the help
+ * of orthogonal_equality.
+ *
+ * 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,
int direction,
- const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
- &offset = dealii::Tensor<1,CellIterator::AccessorType::space_dimension>());
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset);
/**
- * 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.
+ * Same function as above, but accepts a Triangulation or DoFHandler
+ * object @p dof_handler instead of an explicit face iterator range.
+ *
+ * This function will collect periodic face pairs on the highest (i.e.
+ * coarsest) mesh level.
+ *
+ * @author Matthias Maier, 2012
*/
template<typename DH>
- std::map<typename DH::cell_iterator, typename DH::cell_iterator>
- collect_periodic_cell_pairs (const DH &dof_handler,
- const types::boundary_id boundary_component,
- int direction,
- const dealii::Tensor<1,DH::space_dimension>
- &offset = Tensor<1,DH::space_dimension>());
+ std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
+ collect_periodic_face_pairs (const DH &dof_handler, /*TODO: Name*/
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset);
+
/**
<< " is not used in the given triangulation");
-}
+} /*namespace GridTools*/
//---------------------------------------------------------------------------
+#include <deal.II/base/std_cxx1x/array.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
return surface_to_volume_mapping;
}
- // Internally used in
- // collect_periodic_cell_pairs.
- template<typename CellIterator>
- inline bool orthogonal_equality (const CellIterator &cell1,
- const CellIterator &cell2,
- const int direction,
- const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
- &offset)
+
+
+ /*
+ * Internally used in orthogonal_equality
+ *
+ * An orthogonal equality test for points:
+ *
+ * point1 and point2 are considered equal, if
+ * (point1 + offset) - point2
+ * is parallel to the unit vector in <direction>
+ */
+ template<int spacedim>
+ inline bool orthogonal_equality (const dealii::Point<spacedim> &point1,
+ const dealii::Point<spacedim> &point2,
+ const int direction,
+ const dealii::Tensor<1,spacedim> &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)
+ Assert (0<=direction && direction<spacedim,
+ ExcIndexRange (direction, 0, spacedim));
+ for (int i = 0; i < spacedim; ++i)
{
- // Only compare
- // coordinate-components != direction:
+ // Only compare coordinate-components != direction:
if (i == direction)
continue;
- for (unsigned int v = 0; v < GeometryInfo<dim>::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;
+ if (fabs(point1(i) + offset[i] - point2(i)) > 1.e-10)
+ return false;
}
return true;
}
- template<typename CellIterator>
- std::map<CellIterator, CellIterator>
- collect_periodic_cell_pairs (const CellIterator &begin,
- const typename identity<CellIterator>::type &end,
- const types::boundary_id boundary_component,
- int direction,
- const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
- &offset)
+
+ /*
+ * Internally used in orthogonal_equality
+ *
+ * A lookup table to transform vertex matchings to orientation flags of
+ * the form (face_orientation, face_flip, face_rotation)
+ *
+ * See the comment on the next function as well as the detailed
+ * documentation of make_periodicity_constraints and
+ * collect_periodic_face_pairs for details
+ */
+ template<int dim> struct OrientationLookupTable {};
+
+ template<> struct OrientationLookupTable<1>
{
- static const int space_dim = CellIterator::AccessorType::space_dimension;
- Assert (0<=direction && direction<space_dim,
- ExcIndexRange (direction, 0, space_dim));
+ typedef std_cxx1x::array<unsigned int, GeometryInfo<1>::vertices_per_face> MATCH_T;
+ static inline std::bitset<3> lookup (const MATCH_T &)
+ {
+ // The 1D case is trivial
+ return 4; // [true ,false,false]
+ }
+ };
+
+ template<> struct OrientationLookupTable<2>
+ {
+ typedef std_cxx1x::array<unsigned int, GeometryInfo<2>::vertices_per_face> MATCH_T;
+ static inline std::bitset<3> lookup (const MATCH_T &matching)
+ {
+ // In 2D matching faces (=lines) results in two cases: Either
+ // they are aligned or flipped. We store this "line_flip"
+ // property somewhat sloppy as "face_flip"
+ // (always: face_orientation = true, face_rotation = false)
+
+ static const MATCH_T m_tff = {{ 0 , 1 }};
+ if (matching == m_tff) return 1; // [true ,false,false]
+ static const MATCH_T m_ttf = {{ 1 , 0 }};
+ if (matching == m_ttf) return 3; // [true ,true ,false]
+ AssertThrow(false, ExcInternalError());
+ }
+ };
+
+ template<> struct OrientationLookupTable<3>
+ {
+ typedef std_cxx1x::array<unsigned int, GeometryInfo<3>::vertices_per_face> MATCH_T;
+ static inline std::bitset<3> lookup (const MATCH_T &matching)
+ {
+ // The full fledged 3D case. *Yay*
+ // See the documentation in include/deal.II/base/geometry_info.h
+ // as well as the actual implementation in source/grid/tria.cc
+ // for more details...
+
+ static const MATCH_T m_tff = {{ 0 , 1 , 2 , 3 }};
+ if (matching == m_tff) return 1; // [true ,false,false]
+ static const MATCH_T m_tft = {{ 1 , 3 , 0 , 2 }};
+ if (matching == m_tft) return 5; // [true ,false,true ]
+ static const MATCH_T m_ttf = {{ 3 , 2 , 1 , 0 }};
+ if (matching == m_ttf) return 3; // [true ,true ,false]
+ static const MATCH_T m_ttt = {{ 2 , 0 , 3 , 1 }};
+ if (matching == m_ttt) return 7; // [true ,true ,true ]
+ static const MATCH_T m_fff = {{ 0 , 2 , 1 , 3 }};
+ if (matching == m_fff) return 0; // [false,false,false]
+ static const MATCH_T m_fft = {{ 2 , 3 , 0 , 1 }};
+ if (matching == m_fft) return 4; // [false,false,true ]
+ static const MATCH_T m_ftf = {{ 3 , 1 , 2 , 0 }};
+ if (matching == m_ftf) return 2; // [false,true ,false]
+ static const MATCH_T m_ftt = {{ 1 , 0 , 3 , 2 }};
+ if (matching == m_ftt) return 6; // [false,true ,true ]
+ AssertThrow(false, ExcInternalError());
+ }
+ };
+
- std::set<CellIterator> cells1;
- std::set<CellIterator> 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)
+ template<typename FaceIterator>
+ inline bool
+ orthogonal_equality (std::bitset<3> &orientation,
+ const FaceIterator &face1,
+ const FaceIterator &face2,
+ const int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ {
+ static const int dim = FaceIterator::AccessorType::dimension;
+
+ // Do a full matching of the face vertices:
+
+ std_cxx1x::
+ array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
+
+ std::set<unsigned int> face2_vertices;
+ for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
+ face2_vertices.insert(i);
+
+ for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
{
- if (cell->face(2*direction)->at_boundary() &&
- cell->face(2*direction)->boundary_indicator() == boundary_component)
+ for (std::set<unsigned int>::iterator it = face2_vertices.begin();
+ it != face2_vertices.end();
+ it++)
{
- cells1.insert(cell);
- }
- if (cell->face(2*direction+1)->at_boundary() &&
- cell->face(2*direction+1)->boundary_indicator() == boundary_component)
- {
- cells2.insert(cell);
+ if (orthogonal_equality(face1->vertex(i),face2->vertex(*it),
+ direction, offset))
+ {
+ matching[i] = *it;
+ face2_vertices.erase(it);
+ break; // jump out of the innermost loop
+ }
}
}
- Assert (cells1.size() == cells2.size(),
- ExcMessage ("Unmatched cells on periodic boundaries"));
+ // And finally, a lookup to determine the ordering bitmask:
+ if (face2_vertices.empty())
+ orientation = OrientationLookupTable<dim>::lookup(matching);
+
+ return face2_vertices.empty();
+ }
+
+
+
+ template<typename FaceIterator>
+ inline bool
+ orthogonal_equality (const FaceIterator &face1,
+ const FaceIterator &face2,
+ const int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ {
+ // Call the function above with a dummy orientation array
+ std::bitset<3> dummy;
+ return orthogonal_equality (dummy, face1, face2, direction, offset);
+ }
+
+
- std::map<CellIterator, CellIterator> matched_cells;
+ /*
+ * Internally used in 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)
+ {
+ static const int space_dim = FaceIterator::AccessorType::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ Assert (faces1.size() == faces2.size(),
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ typedef std::pair<FaceIterator, std::bitset<3> > ResultPair;
+ std::map<FaceIterator, ResultPair> matched_faces;
- // Match with a complexity of O(n^2).
- // This could be improved...
- typedef typename std::set<CellIterator>::const_iterator SetIterator;
- for (SetIterator it1 = cells1.begin(); it1 != cells1.end(); ++it1)
+ // Match with a complexity of O(n^2). This could be improved...
+ std::bitset<3> orientation;
+ typedef typename std::set<FaceIterator>::const_iterator SetIterator;
+ for (SetIterator it1 = faces1.begin(); it1 != faces1.end(); ++it1)
{
- for (SetIterator it2 = cells2.begin(); it2 != cells2.end(); ++it2)
+ for (SetIterator it2 = faces2.begin(); it2 != faces2.end(); ++it2)
{
- if (orthogonal_equality(*it1, *it2, direction, offset))
+ if (GridTools::orthogonal_equality(orientation, *it1, *it2,
+ direction, offset))
{
- // We have a match, so insert the
- // matching pairs and remove the matched
- // cell in cells2 to speed up the
+ // We have a match, so insert the matching pairs and
+ // remove the matched cell in faces2 to speed up the
// matching:
- matched_cells[*it1] = *it2;
- cells2.erase(it2);
+ matched_faces[*it1] = std::make_pair(*it2, orientation);
+ faces2.erase(it2);
break;
}
}
}
- Assert (matched_cells.size() == cells1.size() &&
- cells2.size() == 0,
- ExcMessage ("Unmatched cells on periodic boundaries"));
+ AssertThrow (matched_faces.size() == faces1.size() && faces2.size() == 0,
+ ExcMessage ("Unmatched faces on periodic boundaries"));
- return matched_cells;
+ return matched_faces;
}
+
+ 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,
+ 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,
+ ExcIndexRange (direction, 0, space_dim));
+
+ // Collect all faces belonging to b_id1 and b_id2:
+
+ std::set<FaceIterator> faces1;
+ std::set<FaceIterator> faces2;
+
+ for (FaceIterator face = begin; face!= end; ++face)
+ {
+ 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"));
+
+ // and call match_periodic_face_pairs that does the actual matching:
+
+ return match_periodic_face_pairs(faces1, faces2, direction, offset);
+ }
+
+
+
template<typename DH>
- std::map<typename DH::cell_iterator, typename DH::cell_iterator>
- collect_periodic_cell_pairs (const DH &dof_handler,
- const types::boundary_id boundary_component,
- int direction,
- const dealii::Tensor<1,DH::space_dimension>
- &offset)
+ std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
+ collect_periodic_face_pairs (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset)
{
- return collect_periodic_cell_pairs<typename DH::cell_iterator>
- (dof_handler.begin_active(),
- dof_handler.end(),
- boundary_component,
- direction,
- 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:
+
+ 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)
+ {
+ 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"));
+
+ // and call match_periodic_face_pairs that does the actual matching:
+
+ return match_periodic_face_pairs(faces1, faces2, direction, offset);
}
projects / dealii-svn.git / commitdiff