* 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 the finite element. This
- * can be used to enforce periodicity in only one variable in a system
- * of equations.
+ * is different from the default value, it is assumed that the number of
+ * entries equals the number of components of the finite element. This
+ * can be used to enforce periodicity in only one variable in a system of
+ * equations.
*
* @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
* and any combination of that...
* @endcode
*
- * More information on the topic can be found in the
- * @ref GlossFaceOrientation "glossary" article.
+ * Optionally a matrix @p matrix along with an std::vector @p
+ * first_vector_components 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 the std::vector first_vector_components is non empty the
+ * matrix is interpreted as a rotation matrix that is applied to all
+ * vector valued blocks listed in first_vector_components of the
+ * FESystem.
*
- * @note For DoFHandler objects that are built on a
- * parallel::distributed::Triangulation object
- * parallel::distributed::Triangulation::add_periodicity has to be called
- * before.
+ * Detailed information can be found in the @see @ref
+ * GlossPeriodicConstraints "Glossary entry on periodic boundary
+ * conditions".
*
- * @author Matthias Maier, 2012
+ * @todo: Reference to soon be written example step and glossary article.
+ *
+ * @author Matthias Maier, 2012, 2014
*/
template<typename FaceIterator>
void
const ComponentMask &component_mask = ComponentMask(),
const bool face_orientation = true,
const bool face_flip = false,
- const bool face_rotation = false);
+ const bool face_rotation = false,
+ const FullMatrix<double> &matrix = FullMatrix<double>(),
+ const std::vector<unsigned int> &first_vector_components = std::vector<unsigned int>());
+
+
+
+ /**
+ * Insert the (algebraic) constraints due to periodic boundary
+ * conditions into a ConstraintMatrix @p constraint_matrix.
+ *
+ * This is the main high level interface for above low level variant of
+ * make_periodicity_constraints. It takes an std::vector @p
+ * periodic_faces as argument and applies above
+ * make_periodicity_constraints on each entry. The std::vector @p
+ * periodic_faces 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.
+ *
+ * @see @ref GlossPeriodicConstraints "Glossary entry on periodic
+ * boundary conditions" for further information.
+ *
+ * @author Daniel Arndt, Matthias Maier, 2013, 2014
+ */
+ 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());
* 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.
+ * make_periodicity_constraints function.
*
* Define a 'first' boundary as all boundary faces having boundary_id
* @p b_id1 and a 'second' boundary consisting of all faces belonging
* function will be used for which the respective flag was set in the
* component mask.
*
- * @note For DoFHandler objects that are built on a
- * parallel::distributed::Triangulation object
- * parallel::distributed::Triangulation::add_periodicity has to be called
- * before.
+ * @note: This function is a convenience wrapper. It internally calls
+ * GridTools::collect_periodic_faces with the supplied paramaters and
+ * feeds the output to above make_periodicity_constraints variant. If you
+ * need more functionality use GridTools::collect_periodic_faces
+ * directly.
*
- * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+ * @see @ref GlossPeriodicConstraints "Glossary entry on periodic
+ * boundary conditions" for further information.
*
* @author Matthias Maier, 2012
*/
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 '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 For DoFHandler objects that are built on a
- * parallel::distributed::Triangulation object
- * parallel::distributed::Triangulation::add_periodicity has to be called
- * before.
- *
- * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
- *
- * @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());
-
/**
* This compatibility version of make_periodicity_constraints only works
* meshes with cells not in @ref GlossFaceOrientation
* "standard orientation".
*
- * @note For DoFHandler objects that are built on a
- * parallel::distributed::Triangulation object
- * parallel::distributed::Triangulation::add_periodicity has to be called
- * before.
+ * @note: This function is a convenience wrapper. It internally calls
+ * GridTools::collect_periodic_faces with the supplied paramaters and
+ * feeds the output to above make_periodicity_constraints variant. If you
+ * need more functionality use GridTools::collect_periodic_faces
+ * directly.
*
- * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+ * @see @ref GlossPeriodicConstraints "Glossary entry on periodic
+ * boundary conditions" for further information.
*/
template<typename DH>
void
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 '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 version of make_periodicity_constraints will not work on
- * meshes with cells not in @ref GlossFaceOrientation
- * "standard orientation".
- *
- * @note For DoFHandler objects that are built on a
- * parallel::distributed::Triangulation object
- * parallel::distributed::Triangulation::add_periodicity has to be called
- * before.
- *
- * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+ * @deprecated This function is deprecated. Use
+ * GridTools::collect_periodic_faces in conjunction with
+ * make_periodicity_constraints instead.
*/
template<typename DH>
void
make_periodicity_constraints
(const DH &dof_handler,
- const types::boundary_id b_id,
+ 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());
+ const ComponentMask &component_mask = ComponentMask()) DEAL_II_DEPRECATED;
+
/**
- * Same as above but the periodicity information is given by
- * @p periodic_faces. This std::vector can be created by
- * GridTools::collect_periodic_faces.
+ * 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 For DoFHandler objects that are built on a
- * parallel::distributed::Triangulation object
- * parallel::distributed::Triangulation::add_periodicity has to be called
- * before.
+ * @note This version of make_periodicity_constraints will not work on
+ * meshes with cells not in @ref GlossFaceOrientation
+ * "standard orientation".
+ *
+ * @deprecated This function is deprecated. Use
+ * GridTools::collect_periodic_faces in conjunction with
+ * make_periodicity_constraints instead.
*/
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 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()) DEAL_II_DEPRECATED;
//@}
/*@{*/
/**
- * 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
+ * Data type that provides all information necessary to create
+ * periodicity constraints and a periodic p4est forest with respect to
+ * two 'periodic' cell faces.
*/
template<typename CellIterator>
struct PeriodicFacePair
{
+ /**
+ * 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.
+ */
unsigned int face_idx[2];
+
+ /**
+ * 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 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 the std::vector first_vector_components is non empty the
+ * matrix is interpreted as a rotation matrix that is applied to all
+ * vector valued blocks listed in first_vector_components of the
+ * FESystem. For more details see make_periodicity_constraints and the
+ * glossary @ref GlossPeriodicConstraints "glossary entry on periodic
+ * boundary conditions".
+ */
+ FullMatrix<double> matrix;
+
+ /**
+ * A vector of unsigned ints pointing to the first components of all
+ * vector valued blocks that should be rotated by matrix.
+ */
+ std::vector<unsigned int> first_vector_components;
};
+
/**
* An orthogonal equality test for faces.
*
* 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 matrix along with a vector @p
+ * first_vector_components 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 @p first_vector_components is non
+ * empty the matrix is interpreted as a rotation matrix that is applied
+ * to all vector valued blocks listet in @p first_vector_components of
+ * the FESystem. For more details see make_periodicity_constraints and
+ * the glossary @ref GlossPeriodicConstraints "glossary entry on periodic
+ * boundary conditions".
+ *
* @tparam Container A type that satisfies the
* requirements of a mesh container (see @ref GlossMeshAsAContainer).
*
* times with different boundary ids to generate a vector with all periodic
* pairs.
*
- * @author Daniel Arndt, Matthias Maier, 2013
+ * @author Daniel Arndt, Matthias Maier, 2013, 2014
*/
- template <typename Container>
+ template <typename CONTAINER>
void
collect_periodic_faces
- (const Container &container,
+ (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 dealii::Tensor<1,Container::space_dimension> &offset = dealii::Tensor<1,Container::space_dimension>());
+ 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> &matrix = FullMatrix<double>(IdentityMatrix(CONTAINER::space_dimension)),
+ const std::vector<unsigned int> &first_vector_components = std::vector<unsigned int>());
/**
* This function will collect periodic face pairs on the coarsest mesh level
* and add them to @p matched_pairs leaving the original contents intact.
*
+ * Optionally a rotation matrix @p matrix along with a vector @p
+ * first_vector_components 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 @p first_vector_components is non
+ * empty the matrix is interpreted as a rotation matrix that is applied
+ * to all vector valued blocks listet in @p first_vector_components of
+ * the FESystem. For more details see make_periodicity_constraints and
+ * the glossary @ref GlossPeriodicConstraints "glossary entry on periodic
+ * boundary conditions".
+ *
* @tparam Container A type that satisfies the
* requirements of a mesh container (see @ref GlossMeshAsAContainer).
*
* meshes with cells not in @ref GlossFaceOrientation
* "standard orientation".
*
- * @author Daniel Arndt, Matthias Maier, 2013
+ * @author Daniel Arndt, Matthias Maier, 2013, 2014
*/
- template <typename Container>
+ template <typename CONTAINER>
void
collect_periodic_faces
- (const Container &container,
+ (const CONTAINER &container,
const types::boundary_id b_id,
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>());
+ 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> &matrix = FullMatrix<double>(),
+ const std::vector<unsigned int> &first_vector_components = std::vector<unsigned int>());
+
/*@}*/
/**
//
// ---------------------------------------------------------------------
+#include <deal.II/base/std_cxx1x/array.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/table.h>
#include <deal.II/base/template_constraints.h>
{
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),
+ Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
ExcMessage ("Matching periodic cells need to use the same finite element"));
const FiniteElement<dim, spacedim> &fe = face_1->get_fe(face_1_index);
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
// 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(jj, 0, /* It doesn't really matter, just assume
- * we're on the first face...
- */
- face_orientation, face_flip, face_rotation)];
+ cell_to_rotated_face_index[fe.face_to_cell_index
+ (jj, 0, face_orientation, face_flip, face_rotation)];
// And finally constrain the two DoFs respecting component_mask:
if (transformation(i,jj) != 0)
}
+ // Implementation of the low level interface:
+
+
template <typename FaceIterator>
void
make_periodicity_constraints (const FaceIterator &face_1,
const ComponentMask &component_mask,
const bool face_orientation,
const bool face_flip,
- const bool face_rotation)
+ const bool face_rotation,
+ const FullMatrix<double> &matrix,
+ const std::vector<unsigned int> &first_vector_components)
{
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"));
+ Assert(first_vector_components.empty()
+ || (matrix.size(0) == (int)dim && matrix.size(1) == (int)dim),
+ ExcMessage ("first_vector_components is nonempty, so matrix must be a rotation "
+ "matrix exactly of size dim"));
+
+#ifdef DEBUG
+ if (!face_1->has_children())
+ {
+ Assert(face_1->n_active_fe_indices()==1, ExcInternalError());
+ const unsigned int n_dofs = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
+ Assert(first_vector_components.size() != 0 ||
+ ( (matrix.size(0) == 0 && matrix.size(1) == 0) ||
+ (matrix.size(0) == n_dofs && matrix.size(1) == n_dofs) ||
+ (matrix.size(0) == (int)dim && matrix.size(1) == (int)dim) ),
+ ExcMessage ("first_vector_components is empty, so matrix must have either "
+ "size 0 or dim or the size must be equal to the number of DoFs on the face"));
+ }
+#endif
// A lookup table on how to go through the child faces depending on the
// orientation:
component_mask,
face_orientation,
face_flip,
- face_rotation);
+ face_rotation,
+ matrix,
+ first_vector_components);
}
}
else
// otherwise at least one of the two faces is active and
// we need to enter the constraints
{
+ // Build up the transformation matrix:
+ FullMatrix<double> transformation;
+ const unsigned int n_dofs = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
+
+ if (matrix.size(0) == n_dofs)
+ transformation = matrix;
+ else if (!first_vector_components.empty() && matrix.size(0) == (int)dim)
+ {
+ // The matrix describes a rotation and we have to build a
+ // transformation matrix, we assume that for a 0° rotation
+ // we would have to build the identity matrix
+
+ const FiniteElement<dim> &fe1
+ = face_1->get_fe(face_1->nth_active_fe_index(0));
+
+ Quadrature<dim-1> quadrature (fe1.get_unit_face_support_points());
+
+ FEFaceValues<dim> fe_face_values
+ (fe1, quadrature, update_q_points);
+
+ // have an array that stores the location of each vector-dof tuple
+ // we want to rotate.
+ typedef std_cxx1x::array<unsigned int, dim> DoFTuple;
+
+ FullMatrix<double> rot_transformation
+ = IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face);
+
+ for (unsigned int i=0; i<fe1.dofs_per_face; ++i)
+ {
+ std::vector<unsigned int>::const_iterator comp_it
+ = std::find (first_vector_components.begin(),
+ first_vector_components.end(),
+ fe1.face_system_to_component_index(i).first);
+ if (comp_it != first_vector_components.end())
+ {
+ const unsigned int first_vector_component = *comp_it;
+
+ // find corresponding other components of vector
+ DoFTuple vector_dofs;
+ vector_dofs[0] = i;
+
+ Assert(*comp_it+dim<=fe1.n_components(),
+ ExcMessage("Error: the finite element does not have enough components "
+ "to define rotated periodic boundaries."));
+
+ for (unsigned int k=0; k<fe1.dofs_per_face; ++k)
+ if ((k != i)
+ &&
+ (quadrature.point(k) == quadrature.point(i))
+ &&
+ (fe1.face_system_to_component_index(k).first >=
+ first_vector_component)
+ &&
+ (fe1.face_system_to_component_index(k).first <
+ first_vector_component + dim))
+ vector_dofs[fe1.face_system_to_component_index(k).first -
+ first_vector_component]
+ = k;
+
+ for (int i=0; i<dim; ++i)
+ {
+ rot_transformation[vector_dofs[i]][vector_dofs[i]]=0.;
+ for (int j=0; j<dim; ++j)
+ rot_transformation[vector_dofs[i]][vector_dofs[j]]=matrix[i][j];
+ }
+ }
+ }
+
+ transformation = rot_transformation;
+ }
+ else
+ transformation = IdentityMatrix(n_dofs);
+
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)),
- constraint_matrix,
- component_mask,
- face_orientation, face_flip, face_rotation);
+ {
+ FullMatrix<double> inverse(transformation.size(0));
+ inverse.invert(transformation);
+ 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)),
- constraint_matrix,
- component_mask,
- face_orientation, face_flip, face_rotation);
+ {
+ set_periodicity_constraints(face_1, face_2,
+ transformation,
+ constraint_matrix,
+ component_mask,
+ face_orientation, face_flip, face_rotation);
+ }
}
}
+ 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();
+
+ // Loop over all periodic faces...
+ 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]);
+
+ Assert(face_1->at_boundary() && face_2->at_boundary(),
+ ExcInternalError());
+
+ Assert (face_1 != face_2,
+ ExcInternalError());
+
+ // ... and apply the low level make_periodicity_constraints function to
+ // every matching pair:
+ make_periodicity_constraints(face_1,
+ face_2,
+ constraint_matrix,
+ component_mask,
+ it->orientation[0],
+ it->orientation[1],
+ it->orientation[2],
+ it->matrix,
+ it->first_vector_components);
+ }
+ }
+
+
+ // High level interface variants:
+
+
template<typename DH>
void
make_periodicity_constraints (const DH &dof_handler,
- 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();
-
-
- // And apply the low level make_periodicity_constraints function to
- // every matching pair:
- 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]);
-
- Assert(face_1->at_boundary() && face_2->at_boundary(),
- ExcInternalError());
-
- Assert (face_1 != face_2,
- ExcInternalError());
-
- make_periodicity_constraints(face_1,
- face_2,
- constraint_matrix,
- component_mask,
- it->orientation[0],
- it->orientation[1],
- it->orientation[2]);
- }
- }
-
-
-
namespace internal
{
namespace Assembler
const DH::face_iterator &,
dealii::ConstraintMatrix &,
const ComponentMask &,
- bool, bool, bool);
+ bool, bool, bool,
+ const FullMatrix<double> &,
+ const std::vector<unsigned int> &);
+
+ template
+ void
+ DoFTools::make_periodicity_constraints<DH>
+ (const std::vector<GridTools::PeriodicFacePair<DH::cell_iterator> > &,
+ dealii::ConstraintMatrix &,
+ const ComponentMask &);
+
template
void
dealii::Tensor<1,DH::space_dimension> &,
dealii::ConstraintMatrix &,
const ComponentMask &);
-
- template
- void
- DoFTools::make_periodicity_constraints<DH>
- (const std::vector<GridTools::PeriodicFacePair<DH::cell_iterator> > &,
- dealii::ConstraintMatrix &,
- const ComponentMask &);
-
#endif
}
* An orthogonal equality test for points:
*
* point1 and point2 are considered equal, if
- * (point1 + offset) - point2
+ * matrix.(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)
+ const dealii::Tensor<1,spacedim> &offset,
+ const FullMatrix<double> &matrix)
{
Assert (0<=direction && direction<spacedim,
ExcIndexRange (direction, 0, spacedim));
if (i == direction)
continue;
- if (fabs(point1(i) + offset[i] - point2(i)) > 1.e-10)
+ double transformed_p1_comp=0.;
+
+ for (int j = 0; j < spacedim; ++j)
+ transformed_p1_comp += matrix(i,j)*point1(j) + offset[j];
+
+ if (fabs(transformed_p1_comp - point2(i)) > 1.e-10)
return false;
}
+
return true;
}
const FaceIterator &face1,
const FaceIterator &face2,
const int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
+ const FullMatrix<double> &matrix)
{
static const int dim = FaceIterator::AccessorType::dimension;
it++)
{
if (orthogonal_equality(face1->vertex(i),face2->vertex(*it),
- direction, offset))
+ direction, offset, matrix))
{
matching[i] = *it;
face2_vertices.erase(it);
orthogonal_equality (const FaceIterator &face1,
const FaceIterator &face2,
const int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
+ const FullMatrix<double> &matrix)
{
// Call the function above with a dummy orientation array
std::bitset<3> dummy;
- return orthogonal_equality (dummy, face1, face2, direction, offset);
+ return orthogonal_equality (dummy, face1, face2, direction, offset, matrix);
}
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 dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset)
+ const int direction,
+ std::vector<PeriodicFacePair<CellIterator> > &matched_pairs,
+ const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset,
+ const FullMatrix<double> &matrix,
+ const std::vector<unsigned int> &first_vector_components)
{
static const int space_dim = CellIterator::AccessorType::space_dimension;
Assert (0<=direction && direction<space_dim,
if (GridTools::orthogonal_equality(orientation,
cell1->face(face_idx1),
cell2->face(face_idx2),
- direction, offset))
+ direction, offset,
+ matrix))
{
// 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};
+ = {{cell1, cell2},{face_idx1, face_idx2}, orientation, matrix, first_vector_components};
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 dealii::Tensor<1,CONTAINER::space_dimension> &offset)
+ const Tensor<1,CONTAINER::space_dimension> &offset,
+ const FullMatrix<double> &matrix,
+ const std::vector<unsigned int> &first_vector_components)
{
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);
+ match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs,
+ offset, matrix, first_vector_components);
}
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 dealii::Tensor<1,CONTAINER::space_dimension> &offset)
+ const Tensor<1,CONTAINER::space_dimension> &offset,
+ const FullMatrix<double> &matrix,
+ const std::vector<unsigned int> &first_vector_components)
{
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);
+ match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs,
+ offset, matrix, first_vector_components);
#ifdef DEBUG
//check for standard orientation
const X::active_face_iterator&,
const X::active_face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
template
bool orthogonal_equality<X::face_iterator> (std::bitset<3> &,
const X::face_iterator&,
const X::face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
template
bool orthogonal_equality<X::active_face_iterator> (const X::active_face_iterator&,
const X::active_face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
template
bool orthogonal_equality<X::face_iterator> (const X::face_iterator&,
const X::face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
#if deal_II_dimension >= 2
const types::boundary_id,
const int,
std::vector<PeriodicFacePair<X::cell_iterator> > &,
- const Tensor<1,X::space_dimension> &);
+ const Tensor<1,X::space_dimension> &,
+ const FullMatrix<double> &,
+ const std::vector<unsigned int> &);
template
void collect_periodic_faces<X> (const X &,
const types::boundary_id,
const int,
std::vector<PeriodicFacePair<X::cell_iterator> > &,
- const Tensor<1,X::space_dimension> &);
+ const Tensor<1,X::space_dimension> &,
+ const FullMatrix<double> &,
+ const std::vector<unsigned int> &);
#endif
#if deal_II_dimension >= 2
namespace GridTools \{
-
template
void
collect_periodic_faces<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> >
const types::boundary_id,
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 Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
+ const FullMatrix<double> &,
+ const std::vector<unsigned int> &);
template
void
const types::boundary_id,
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 Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
+ const FullMatrix<double> &,
+ const std::vector<unsigned int> &);
\}
#endif
#endif
projects / dealii.git / commitdiff