ConstraintMatrix &constraints);
//@}
+
+ /**
+ * @name Periodic Boundary Conditions
+ * @{
+ */
+
+ /**
+ * Insert the (algebraic) constraints due
+ * to periodic boundary conditions into
+ * a ConstraintMatrix @p
+ * constraint_matrix.
+ *
+ * Given a pair of not necessarily
+ * active faces @p face_1 and @p face_2,
+ * this functions constrains all DoFs
+ * associated with the boundary
+ * described by @p face_1 to the
+ * respective DoFs of the boundary
+ * described by @p face_2. More
+ * precisely:
+ *
+ * If @p face_1 and @p face_2 are both
+ * active faces it adds the DoFs of
+ * @p face_1 to the list of constrained
+ * DoFs in @p constraint_matrix and adds
+ * lines to constrain them to the
+ * corresponding values of the DoFs on
+ * @p face_2.
+ * Otherwise it loops recursively over
+ * the children of @p face_1 and @p face_2.
+ *
+ * For this to work @p face_1 and @p face_2
+ * must have the same refinement
+ * history, i.e. either @p face_1 and
+ * @p face_2 must be active faces or
+ * must be isotropically refined and
+ * have the same number of child faces
+ * that recursively obey this rule.
+ * (The anisotropic case is not yet
+ * implemented.)
+ *
+ * This routine only constrains DoFs that
+ * are not already constrained.
+ * If this routine encounters a DoF that
+ * already is constrained (for instance
+ * by Dirichlet boundary conditions),
+ * the old setting of the constraint
+ * (dofs the entry is constrained to,
+ * inhomogeneities) is kept and nothing
+ * happens.
+ *
+ * Furthermore, no DoFs belonging to (or
+ * belonging to any descendant of) @p
+ * face_2 get constrained or get marked
+ * as being constrained.
+ *
+ * The flags in the last parameter,
+ * @p component_mask denote which
+ * components of the finite element space
+ * shall be constrained with periodic
+ * boundary conditions. If it is left as
+ * specified by the default value (i.e. an
+ * empty array), all components are
+ * constrainted. If it is different from
+ * the default value, it is assumed that
+ * the number of entries equals the number
+ * of components in the boundary functions
+ * and the finite element, and those
+ * components in the given boundary
+ * function will be used for which the
+ * respective flag was set in the component
+ * mask.
+ */
+ template<typename FaceIterator>
+ void
+ make_periodicity_constraints (const FaceIterator &face_1,
+ const typename identity<FaceIterator>::type &face_2,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const std::vector<bool> &component_mask = std::vector<bool>());
+
+
+ /**
+ * Insert the (algebraic) constraints due
+ * to periodic boundary conditions into
+ * a ConstraintMatrix @p
+ * constraint_matrix.
+ *
+ * This function serves as a high level
+ * interface for the
+ * make_periodicity_constraints function
+ * that takes face_iterator arguments.
+ *
+ * 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>.
+ *
+ * This function tries to match
+ * all faces belonging to the 'left' and
+ * 'right' boundary with the help of an
+ * orthogonal equality relation:
+ * Two faces do match if their vertices
+ * can be transformed into each other by
+ * parallel translation in @p direction.
+ *
+ * If this matching is successfull it
+ * constrains all DoFs associated with the
+ * 'left' boundary to the respective DoFs
+ * of the 'right' boundary.
+ *
+ * This routine only constrains DoFs that
+ * are not already constrained.
+ * If this routine encounters a DoF that
+ * already is constrained (for instance
+ * by Dirichlet boundary conditions),
+ * the old setting of the constraint
+ * (dofs the entry is constrained to,
+ * inhomogeneities) is kept and nothing
+ * happens.
+ *
+ * Furthermore, no DoFs belonging to the
+ * 'right' boundary get constrained or get
+ * marked as being constrained.
+ *
+ * The flags in the last parameter,
+ * @p component_mask denote which
+ * components of the finite element space
+ * shall be constrained with periodic
+ * boundary conditions. If it is left as
+ * specified by the default value (i.e. an
+ * empty array), all components are
+ * constrainted. If it is different from
+ * the default value, it is assumed that
+ * the number of entries equals the number
+ * of components in the boundary functions
+ * and the finite element, and those
+ * components in the given boundary
+ * function will be used for which the
+ * respective flag was set in the component
+ * mask.
+ */
+ template<typename DH>
+ void
+ make_periodicity_constraints (const DH &dof_handler,
+ const types::boundary_id_t boundary_component,
+ const int direction,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const std::vector<bool> &component_mask = std::vector<bool>());
+
+
+ /**
+ * Same as above but with an optional
+ * argument @p offset.
+ * 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<typename DH>
+ void
+ make_periodicity_constraints (const DH &dof_handler,
+ const types::boundary_id_t boundary_component,
+ const int direction,
+ dealii::Tensor<1,DH::space_dimension>
+ &offset,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const std::vector<bool> &component_mask = std::vector<bool>());
+ //@}
+
+
/**
* Take a vector of values which live on
* cells (e.g. an error per cell) and
}
+ template<typename FaceIterator>
+ void
+ make_periodicity_constraints (const FaceIterator &face_1,
+ const typename identity<FaceIterator>::type &face_2,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const std::vector<bool> &component_mask = std::vector<bool>())
+ {
+ static const int dim = FaceIterator::AccessorType::dimension;
+
+ Assert(face_1->at_boundary() && face_2->at_boundary(),
+ ExcMessage ("Faces for periodicity constraints must be on the boundary"));
+
+
+ // In the case that both faces have
+ // children, we loop over all children
+ // and applu make_periodicty_constrains
+ // recursively:
+ if (face_1->has_children() && face_2->has_children()) {
+ Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
+ face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
+ ExcNotImplemented());
+
+ for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i) {
+ make_periodicity_constraints (face_1->child(i),
+ face_2->child(i),
+ constraint_matrix,
+ component_mask);
+ }
+ return;
+ }
+ // .. otherwise we should be in the case
+ // were both faces are active and have
+ // no children ..
+ Assert (!face_1->has_children() && !face_2->has_children(),
+ ExcNotImplemented());
+
+ Assert (face_1->n_active_fe_indices() == 1 && face_2->n_active_fe_indices() == 1,
+ ExcInternalError());
+
+ // .. then we match the
+ // corresponding DoFs of both faces ..
+ const unsigned int face_1_index = face_1->nth_active_fe_index(0);
+ const unsigned int face_2_index = face_2->nth_active_fe_index(0);
+ Assert ( face_1->get_fe(face_1_index)
+ == face_2->get_fe(face_1_index),
+ ExcMessage ("Matching periodic cells need to use the same finite element"));
+
+ const dealii::FiniteElement<dim> &fe = face_1->get_fe(face_1_index);
+
+ Assert (component_mask.size() == 0 ||
+ component_mask.size() == fe.n_components(),
+ ExcMessage ("The number of components in the mask has to be either "
+ "zero or equal to the number of components in the finite "
+ "element."));
+
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+
+ std::vector<unsigned int> dofs_1(dofs_per_face);
+ std::vector<unsigned int> dofs_2(dofs_per_face);
+
+ face_1->get_dof_indices(dofs_1, face_1_index);
+ face_2->get_dof_indices(dofs_2, face_2_index);
+
+ // .. and constrain them (respecting
+ // component_mask):
+ for (unsigned int i = 0; i < dofs_per_face; ++i) {
+ if (component_mask.size() == 0 ||
+ component_mask[fe.face_system_to_component_index(i).first]) {
+ if(!constraint_matrix.is_constrained(dofs_1[i])) {
+ constraint_matrix.add_line(dofs_1[i]);
+ constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0);
+ }
+ }
+ }
+ }
+
+
+ template<typename DH>
+ void
+ make_periodicity_constraints (const DH &dof_handler,
+ const types::boundary_id_t boundary_component,
+ const int direction,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const std::vector<bool> &component_mask = std::vector<bool>())
+ {
+ Tensor<1,DH::space_dimension> dummy;
+ make_periodicity_constraints (dof_handler,
+ boundary_component,
+ direction,
+ dummy,
+ constraint_matrix,
+ component_mask);
+ }
+
+
+ template<typename DH>
+ void
+ make_periodicity_constraints (const DH &dof_handler,
+ const types::boundary_id_t boundary_component,
+ const int direction,
+ dealii::Tensor<1,DH::space_dimension>
+ &offset,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const std::vector<bool> &component_mask = std::vector<bool>())
+ {
+ static const int space_dim = DH::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ typedef typename DH::cell_iterator CellIterator;
+
+ // We collect matching periodic cells on
+ // the coarsest level:
+ std::map<CellIterator, CellIterator>
+ matched_cells =
+ GridTools::collect_periodic_cell_pairs(dof_handler.begin(0),
+ dof_handler.end(0),
+ boundary_component,
+ direction,
+ offset);
+
+ // And apply the low level
+ // make_periodicity_constraints function
+ // to every matching pair:
+ for (auto it = matched_cells.begin(); it != matched_cells.end(); ++it) {
+ typedef typename DH::face_iterator FaceIterator;
+ FaceIterator face_1 = it->first->face(2*direction);
+ FaceIterator face_2 = it->second->face(2*direction+1);
+
+ Assert(face_1->at_boundary() && face_2->at_boundary(),
+ ExcInternalError());
+
+ Assert (face_1->boundary_indicator() == boundary_component &&
+ face_2->boundary_indicator() == boundary_component,
+ ExcInternalError());
+
+ make_periodicity_constraints(face_1,
+ face_2,
+ constraint_matrix,
+ component_mask);
+ }
+ }
+
+
namespace internal
{
projects / dealii-svn.git / commitdiff