* and any combination of that...
* @endcode
*
- * The optional parameter @p rotation is a rotation matrix that is
- * applied to all vector valued components of an FESystem on @p face_1
- * prior to periodically constraining them to the DoFs of @p face_2.
+ * Optionally a rotation matrix @p matrix along with a flag @p
+ * rotation_matrix can be specified that describes how DoFs on @p face_1
+ * should be modified prior to constraining to the DoFs of @p face_2. If
+ * rotation_matrix is true the matrix is enterpreted as a rotation matrix
+ * applied to all vector valued blocks of an FESystem. For more details
+ * see make_periodicity_constraints and the glossary @ref
+ * GlossPeriodicConstraints "glossary entry on periodic boundary
+ * conditions".
*
* Detailed information can be found in the @see @ref
* GlossPeriodicConstraints "Glossary entry on periodic boundary
const bool face_orientation = true,
const bool face_flip = false,
const bool face_rotation = false,
- const FullMatrix<double> &rotation = IdentityMatrix(FaceIterator::space_dimension));
+ const FullMatrix<double> &matrix = IdentityMatrix(FaceIterator::space_dimension),
+ const bool rotation_matrix = true);
struct PeriodicFacePair
{
/**
- * the cells associated with the two 'periodic' faces
+ * The cells associated with the two 'periodic' faces.
*/
CellIterator cell[2];
/**
- * the local face indices (with respect to the specified cells) of the
- * two 'periodic' faces
+ * The local face indices (with respect to the specified cells) of the
+ * two 'periodic' faces.
*/
unsigned int face_idx[2];
/**
- * the relative orientation of the first face with respect to the
+ * The relative orientation of the first face with respect to the
* second face as described in orthogonal_equality and
* make_periodicity_constraints (and stored as a bitset).
*/
std::bitset<3> orientation;
/**
- * a rotation matrix that describes how vector valued DoFs of the first face
- * should be rotated prior to constraining to the DoFs of the second
- * face, see make_periodicity_constraints.
+ * A matrix that describes how vector valued DoFs of the first face
+ * should be modified prior to constraining to the DoFs of the second
+ * face. If rotation_matrix is true the matrix is enterpreted as a
+ * rotation matrix applied to all vector valued blocks of an FESystem.
+ * For more details see make_periodicity_constraints and the glossary
+ * @ref GlossPeriodicConstraints "glossary entry on periodic
+ * boundary conditions".
*/
- FullMatrix<double> rotation;
+ FullMatrix<double> matrix;
+
+ /**
+ * A flag indicating whether matrix should be interpreted as rotation
+ * matrix or as an interpolation matrix.
+ */
+ bool rotation_matrix;
};
* 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)$.
*
- * Optionally a rotation matrix @p rotation can be specified that
- * describes how vector valued DoFs of the first boundary @p b_id1 should
- * be rotated prior to constraining to the DoFs of the second boundary @p
- * b_id2, see make_periodicity_constraints.
+ * Optionally a rotation matrix @p matrix along with a flag @p
+ * rotation_matrix can be specified that describes how vector valued DoFs
+ * of the first face should be modified prior to constraining to the DoFs
+ * of the second face. If rotation_matrix is true the matrix is
+ * enterpreted as a rotation matrix applied to all vector valued blocks
+ * of an FESystem. For more details see make_periodicity_constraints and
+ * the glossary @ref GlossPeriodicConstraints "glossary entry on periodic
+ * boundary conditions".
*
* @note The created std::vector can be used in
* DoFTools::make_periodicity_constraints and in
const int direction,
std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
const Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(),
- const FullMatrix<double> &rotation = IdentityMatrix(CONTAINER::space_dimension));
+ const FullMatrix<double> &matrix = IdentityMatrix(CONTAINER::space_dimension),
+ const bool rotation_matrix = true);
/**
const int direction,
std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
const dealii::Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(),
- const FullMatrix<double> &rotation = IdentityMatrix(CONTAINER::space_dimension));
+ const FullMatrix<double> &matrix = IdentityMatrix(CONTAINER::space_dimension),
+ const bool rotation_matrix = true);
const unsigned int j =
cell_to_rotated_face_index[fe.face_to_cell_index(identity_constraint_target,
0, /* It doesn't really matter, just assume
- * we're on the first face...
- */
+ * we're on the first face...
+ */
face_orientation, face_flip, face_rotation)];
// if the two aren't already identity constrained (whichever way
const bool face_orientation,
const bool face_flip,
const bool face_rotation,
- const FullMatrix<double> &rotation)
+ const FullMatrix<double> &matrix,
+ const bool rotation_matrix)
{
static const int dim = FaceIterator::AccessorType::dimension;
face_orientation,
face_flip,
face_rotation,
- rotation);
+ matrix,
+ rotation_matrix);
}
}
else
// otherwise at least one of the two faces is active and
// we need to enter the constraints
{
+ // Build up the transformation matrix:
+
+ if (rotation_matrix)
+ AssertThrow(false, ExcNotImplemented());
+// FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face));
+
if (face_2->has_children() == false)
- set_periodicity_constraints(face_2, face_1,
- FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), // TODO
- constraint_matrix,
- component_mask,
- face_orientation, face_flip, face_rotation);
+ {
+ // TODO: This is horrible...
+ FullMatrix<double> inverse;
+ inverse.invert(matrix);
+ set_periodicity_constraints(face_2, face_1,
+ inverse,
+ constraint_matrix,
+ component_mask,
+ face_orientation, face_flip, face_rotation);
+ }
else
+ {
set_periodicity_constraints(face_1, face_2,
- FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), // TODO
+ matrix,
constraint_matrix,
component_mask,
face_orientation, face_flip, face_rotation);
+ }
}
}
void
make_periodicity_constraints
(const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
- &periodic_faces,
+ &periodic_faces,
dealii::ConstraintMatrix &constraint_matrix,
const ComponentMask &component_mask)
{
it->orientation[0],
it->orientation[1],
it->orientation[2],
- it->rotation);
+ it->matrix,
+ it->rotation_matrix);
}
}
dealii::ConstraintMatrix &,
const ComponentMask &,
bool, bool, bool,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const bool);
template
void
match_periodic_face_pairs
(std::set<std::pair<CellIterator, unsigned int> > &pairs1,
std::set<std::pair<typename identity<CellIterator>::type, unsigned int> > &pairs2,
- const int direction,
- std::vector<PeriodicFacePair<CellIterator> > &matched_pairs,
+ const int direction,
+ std::vector<PeriodicFacePair<CellIterator> > &matched_pairs,
const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset,
- const FullMatrix<double> &rotation)
+ const FullMatrix<double> &matrix,
+ const bool rotation_matrix)
{
static const int space_dim = CellIterator::AccessorType::space_dimension;
Assert (0<=direction && direction<space_dim,
// remove the matched cell in pairs2 to speed up the
// matching:
const PeriodicFacePair<CellIterator> matched_face
- = {{cell1, cell2},{face_idx1, face_idx2}, orientation, rotation};
+ = {{cell1, cell2},{face_idx1, face_idx2}, orientation, matrix, rotation_matrix};
matched_pairs.push_back(matched_face);
pairs2.erase(it2);
++n_matches;
template<typename CONTAINER>
void
collect_periodic_faces
- (const CONTAINER &container,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
+ (const CONTAINER &container,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
const Tensor<1,CONTAINER::space_dimension> &offset,
- const FullMatrix<double> &rotation)
+ const FullMatrix<double> &matrix,
+ const bool rotation_matrix)
{
static const int dim = CONTAINER::dimension;
static const int space_dim = CONTAINER::space_dimension;
ExcMessage ("Unmatched faces on periodic boundaries"));
// and call match_periodic_face_pairs that does the actual matching:
- match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset, rotation);
+ match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs,
+ offset, matrix, rotation_matrix);
}
template<typename CONTAINER>
void
collect_periodic_faces
- (const CONTAINER &container,
- const types::boundary_id b_id,
- const int direction,
+ (const CONTAINER &container,
+ const types::boundary_id b_id,
+ const int direction,
std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
const Tensor<1,CONTAINER::space_dimension> &offset,
- const FullMatrix<double> &rotation)
+ const FullMatrix<double> &matrix,
+ const bool rotation_matrix)
{
static const int dim = CONTAINER::dimension;
static const int space_dim = CONTAINER::space_dimension;
#endif
// and call match_periodic_face_pairs that does the actual matching:
- match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset, rotation);
+ match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs,
+ offset, matrix, rotation_matrix);
#ifdef DEBUG
//check for standard orientation
const int,
std::vector<PeriodicFacePair<X::cell_iterator> > &,
const Tensor<1,X::space_dimension> &,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const bool);
template
void collect_periodic_faces<X> (const X &,
const int,
std::vector<PeriodicFacePair<X::cell_iterator> > &,
const Tensor<1,X::space_dimension> &,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const bool);
#endif
#if deal_II_dimension >= 2
namespace GridTools \{
-
template
void
collect_periodic_faces<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> >
const int,
std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const bool);
template
void
const int,
std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
- const FullMatrix<double> &);
-
+ const FullMatrix<double> &,
+ const bool);
\}
#endif
#endif