*/
/**
- * 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 (see @ref GlossComponentMask)
- * 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 all components are
- * constrained. 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.
+ * Insert the (algebraic) constraints due to periodic boundary
+ * conditions into a ConstraintMatrix @p constraint_matrix.
*
- * @note This function will not work
- * for DoFHandler objects that are
- * built on a
- * parallel::distributed::Triangulation
- * object.
+ * 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. This happens on a purely algebraic level, meaning,
+ * the global DoF with (local face) index <tt>i</tt> on @p face_1 gets
+ * constraint to the DoF with (local face) index <tt>i</tt> on @p face_2
+ * (possibly corrected for orientation, see below).
+ *
+ * Otherwise, if @p face_1 and @p face_2 are not active faces, this
+ * function 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 @p component_mask (see @ref GlossComponentMask)
+ * 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 all components are constrained. 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.
+ *
+ * @p face_orientation, @p face_flip and @p face_rotation describe an
+ * orientation that should be applied to @p face_1 prior to matching and
+ * constraining DoFs. More precisely, this matches local face DoF indices
+ * in the following manner:
+ *
+ * In 2d: <tt>face_orientation</tt> must always be <tt>true</tt>,
+ * <tt>face_rotation</tt> is always <tt>false</tt>, and face_flip has the
+ * meaning of <tt>line_flip</tt>; this implies e.g. for <tt>Q1</tt>:
+ *
+ * @code
+ *
+ * face_orientation = true, face_flip = false, face_rotation = false:
+ *
+ * face1: face2:
+ *
+ * 1 1
+ * | <--> |
+ * 0 0
+ *
+ * Resulting constraints: 0 <-> 0, 1 <-> 1
+ *
+ * (Numbers denote local face DoF indices.)
+ *
+ *
+ * face_orientation = true, face_flip = true, face_rotation = false:
+ *
+ * face1: face2:
+ *
+ * 0 1
+ * | <--> |
+ * 1 0
+ *
+ * Resulting constraints: 1 <-> 0, 0 <-> 1
+ * @endcode
+ *
+ * And simliarly for the case of Q1 in 3d:
+ *
+ * @code
+ *
+ * face_orientation = true, face_flip = false, face_rotation = false:
+ *
+ * face1: face2:
+ *
+ * 2 - 3 2 - 3
+ * | | <--> | |
+ * 0 - 1 0 - 1
+ *
+ * Resulting constraints: 0 <-> 0, 1 <-> 1, 2 <-> 2, 3 <-> 3
+ *
+ * (Numbers denote local face DoF indices.)
+ *
+ *
+ * face_orientation = false, face_flip = false, face_rotation = false:
+ *
+ * face1: face2:
+ *
+ * 1 - 3 2 - 3
+ * | | <--> | |
+ * 0 - 2 0 - 1
+ *
+ * Resulting constraints: 0 <-> 0, 2 <-> 1, 1 <-> 2, 3 <-> 3
+ *
+ *
+ * face_orientation = true, face_flip = true, face_rotation = false:
+ *
+ * face1: face2:
+ *
+ * 1 - 0 2 - 3
+ * | | <--> | |
+ * 3 - 2 0 - 1
+ *
+ * Resulting constraints: 3 <-> 0, 2 <-> 1, 1 <-> 2, 0 <-> 3
+ *
+ *
+ * face_orientation = true, face_flip = false, face_rotation = true
+ *
+ * face1: face2:
+ *
+ * 0 - 2 2 - 3
+ * | | <--> | |
+ * 1 - 3 0 - 1
+ *
+ * Resulting constraints: 1 <-> 0, 3 <-> 1, 0 <-> 2, 2 <-> 3
+ *
+ * and any combination of that...
+ * @endcode
+ *
+ * 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.
+ *
+ * @author Matthias Maier, 2012
*/
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());
-
-
- /**
- * 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 (see @ref GlossComponentMask)
- * 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 all components are
- * constrained. 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.
+ const ComponentMask &component_mask = ComponentMask(),
+ bool face_orientation = true,
+ bool face_flip = false,
+ bool face_rotation = false);
+
+
+ /**
+ * Insert the (algebraic) constraints due to periodic boundary
+ * conditions into a ConstraintMatrix @p constraint_matrix.
*
- * @note This function will not work
- * for DoFHandler objects that are
- * built on a
- * parallel::distributed::Triangulation
- * object.
+ * This function serves as a high level interface for the
+ * make_periodicity_constraints function that takes face_iterator
+ * arguments.
+ *
+ * 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 @p orthogonal_equality.
+ *
+ * If this matching is successfull it constrains all DoFs associated
+ * with the 'first' boundary to the respective DoFs of the 'second'
+ * boundary respecting the relative orientation of the two faces.
+ *
+ * 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 'second' boundary get
+ * constrained or get marked as being constrained.
+ *
+ * The flags in the last parameter, @p component_mask (see @ref
+ * GlossComponentMask) 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 all components are constrained. 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.
+ *
+ * @note This function will not work for DoFHandler objects that are
+ * built on a parallel::distributed::Triangulation object.
+ *
+ * @author Matthias Maier, 2012
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id boundary_component,
- 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());
/**
- * 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)$.
+ * Same as above but with an optional argument @p offset.
*
- * @note This function will not work
- * for DoFHandler objects that are
- * built on a
- * parallel::distributed::Triangulation
- * object.
+ * 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 via
+ * @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.
+ *
+ * @author Matthias Maier, 2012
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id boundary_component,
- 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 FaceIterator>
void
make_periodicity_constraints (const FaceIterator &face_1,
const typename identity<FaceIterator>::type &face_2,
dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask)
+ const ComponentMask &component_mask,
+ bool face_orientation,
+ bool face_flip,
+ bool face_rotation)
{
static const int dim = FaceIterator::AccessorType::dimension;
+ Assert( (dim != 1)
+ ||
+ (face_orientation == true && face_flip == false && face_rotation == false),
+ ExcMessage ("The supplied orientation (face_orientation, face_flip, face_rotation) is invalid for 1D"));
+
+ Assert( (dim != 2)
+ ||
+ (face_orientation == true && face_rotation == false),
+ ExcMessage ("The supplied orientation (face_orientation, face_flip, face_rotation) is invalid for 2D"));
+
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 apply make_periodicty_constrains
- // recursively:
+ // A lookup table on how to go through the child faces depending on the
+ // orientation:
+
+ static const int lookup_table_2d[2][2]
+ // flip
+ = { {0, 1}, // false
+ {1, 0}, // true
+ };
+
+ static const int lookup_table_3d[2][2][2][4]
+ // orientation flip rotation
+ = { { { {0, 2, 1, 3}, // false false false
+ {2, 3, 0, 1}, }, // false false true
+ { {3, 1, 2, 0}, // false true false
+ {1, 0, 3, 2}, }, }, // false true true
+ { { {0, 1, 2, 3}, // true false false
+ {1, 3, 0, 2}, }, // true false true
+ { {3, 2, 1, 0}, // true true false
+ {2, 0, 3, 1}, }, }, // true true true
+ };
+
+ // In the case that both faces have children, we loop over all
+ // children and apply make_periodicty_constrains recursively:
if (face_1->has_children() && face_2->has_children())
{
Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
{
+ // Lookup the index for the second face
+ unsigned int j;
+ switch(dim)
+ {
+ case 2:
+ j = lookup_table_2d[face_flip][i];
+ break;
+ case 3:
+ j = lookup_table_3d[face_orientation][face_flip][face_rotation][i];
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+
make_periodicity_constraints (face_1->child(i),
- face_2->child(i),
+ face_2->child(j),
constraint_matrix,
- component_mask);
+ component_mask,
+ face_orientation,
+ face_flip,
+ face_rotation);
}
return;
}
- // .. otherwise we should be in the case
- // were both faces are active and have
- // no children ..
+
+ // .. 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 ..
+ // .. 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"));
+ 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);
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):
+
+ // Well, this is a hack:
+ //
+ // There is no
+ // face_to_face_index(face_index,
+ // face_orientation,
+ // face_flip,
+ // face_rotation)
+ // function in FiniteElementData, so we have to use
+ // face_to_cell_index(face_index, face
+ // face_orientation,
+ // face_flip,
+ // face_rotation)
+ // But this will give us an index on a cell - something we cannot work
+ // with. But luckily we can match them back :-]
+
+ std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
+
+ // Build up a cell to face index for face_2:
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int cell_index = fe.
+ face_to_cell_index(i,
+ 0, // It doesn't really matter, just assume
+ // we're on the first face...
+ true, false, false // default orientation
+ );
+ cell_to_rotated_face_index[cell_index] = i;
+ }
+
for (unsigned int i = 0; i < dofs_per_face; ++i)
{
+ // Query the correct face_index on face_2 respecting the given
+ // orientation:
+ const unsigned int j =
+ cell_to_rotated_face_index[fe.
+ face_to_cell_index(i,
+ 0, // It doesn't really matter, just assume
+ // we're on the first face...
+ face_orientation,
+ face_flip,
+ face_rotation)];
+
+ // And finally constrain the two DoFs respecting component_mask:
if ((component_mask.n_selected_components(fe.n_components()) == fe.n_components())
||
(component_mask[fe.face_system_to_component_index(i).first] == true))
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);
+ constraint_matrix.add_entry(dofs_1[i], dofs_2[j], 1.0);
}
}
}
}
+
template<typename DH>
void
make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id boundary_component,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
const int direction,
dealii::ConstraintMatrix &constraint_matrix,
const ComponentMask &component_mask)
{
Tensor<1,DH::space_dimension> dummy;
make_periodicity_constraints (dof_handler,
- boundary_component,
+ b_id1,
+ b_id2,
direction,
dummy,
constraint_matrix,
}
+
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id boundary_component,
- const int direction,
- dealii::Tensor<1,DH::space_dimension>
- &offset,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask)
+ 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)
{
static const int space_dim = DH::space_dimension;
Assert (0<=direction && direction<space_dim,
ExcMessage ("This function can not be used with distributed triangulations."
"See the documentation for more information."));
- 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 (typename std::map<CellIterator, CellIterator>::iterator it = matched_cells.begin();
+ typedef typename DH::face_iterator FaceIterator;
+ typedef std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > > FaceMap;
+
+ // Collect matching periodic cells on the coarsest level:
+ FaceMap matched_cells = GridTools::collect_periodic_face_pairs(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();
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);
+ const FaceIterator &face_1 = it->first;
+ const FaceIterator &face_2 = it->second.first;
+ const std::bitset<3> &orientation = it->second.second;
Assert(face_1->at_boundary() && face_2->at_boundary(),
ExcInternalError());
- Assert (face_1->boundary_indicator() == boundary_component &&
- face_2->boundary_indicator() == boundary_component,
+ Assert (face_1->boundary_indicator() == b_id1 &&
+ face_2->boundary_indicator() == b_id2,
ExcInternalError());
make_periodicity_constraints(face_1,
face_2,
constraint_matrix,
- component_mask);
+ component_mask,
+ orientation[0],
+ orientation[1],
+ orientation[2]);
}
}
projects / dealii-svn.git / commitdiff