From 2d029a9872898a39b0df66d628ec85d21dc0ea83 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 19 Nov 2014 21:14:12 +0100 Subject: [PATCH] Implement rotation for vector valued dofs on periodic bc This commit completes the interface for periodic boundary conditions. It exposes the interpolation matrix (that is internally used to satisfy hanging node constraints on periodic boundaries) and provides a possibility to generate it from a rotation matrix and an index for the vector component in question (in a similar vein as we already do this for different vector valued finite elements). Signed-off-by: Matthias Maier --- include/deal.II/dofs/dof_tools.h | 169 ++++++++-------- include/deal.II/grid/grid_tools.h | 84 ++++++-- source/dofs/dof_tools_constraints.cc | 227 ++++++++++++++++------ source/dofs/dof_tools_constraints.inst.in | 20 +- source/grid/grid_tools.cc | 64 +++--- source/grid/grid_tools.inst.in | 30 ++- 6 files changed, 401 insertions(+), 193 deletions(-) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index 85bdd1a1ff..e7a6961463 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -848,10 +848,10 @@ namespace DoFTools * 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 @@ -945,15 +945,21 @@ namespace DoFTools * 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 void @@ -964,7 +970,39 @@ namespace DoFTools 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 &matrix = FullMatrix(), + const std::vector &first_vector_components = std::vector()); + + + + /** + * 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 + void + make_periodicity_constraints + (const std::vector > + &periodic_faces, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask = ComponentMask()); @@ -973,8 +1011,7 @@ namespace DoFTools * 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 @@ -1004,12 +1041,14 @@ namespace DoFTools * 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 */ @@ -1024,35 +1063,6 @@ namespace DoFTools 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 - 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 @@ -1069,12 +1079,14 @@ namespace DoFTools * 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 void @@ -1086,53 +1098,54 @@ namespace DoFTools 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 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 void make_periodicity_constraints - (const std::vector > - &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; //@} diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 460c988618..fc1b9507de 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1052,18 +1052,51 @@ namespace GridTools /*@{*/ /** - * 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 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 matrix; + + /** + * A vector of unsigned ints pointing to the first components of all + * vector valued blocks that should be rotated by matrix. + */ + std::vector first_vector_components; }; + /** * An orthogonal equality test for faces. * @@ -1173,6 +1206,16 @@ namespace GridTools * 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). * @@ -1186,17 +1229,19 @@ namespace GridTools * 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 + template 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 > &matched_pairs, - const dealii::Tensor<1,Container::space_dimension> &offset = dealii::Tensor<1,Container::space_dimension>()); + std::vector > &matched_pairs, + const Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(), + const FullMatrix &matrix = FullMatrix(IdentityMatrix(CONTAINER::space_dimension)), + const std::vector &first_vector_components = std::vector()); /** @@ -1213,6 +1258,16 @@ namespace GridTools * 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). * @@ -1220,16 +1275,19 @@ namespace GridTools * meshes with cells not in @ref GlossFaceOrientation * "standard orientation". * - * @author Daniel Arndt, Matthias Maier, 2013 + * @author Daniel Arndt, Matthias Maier, 2013, 2014 */ - template + template void collect_periodic_faces - (const Container &container, + (const CONTAINER &container, const types::boundary_id b_id, const int direction, - std::vector > &matched_pairs, - const dealii::Tensor<1,Container::space_dimension> &offset = dealii::Tensor<1,Container::space_dimension>()); + std::vector > &matched_pairs, + const dealii::Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(), + const FullMatrix &matrix = FullMatrix(), + const std::vector &first_vector_components = std::vector()); + /*@}*/ /** diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 56aaa3616e..cf52dcac5e 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -13,6 +13,7 @@ // // --------------------------------------------------------------------- +#include #include #include #include @@ -1711,7 +1712,7 @@ namespace DoFTools { 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 &fe = face_1->get_fe(face_1_index); @@ -1834,8 +1835,8 @@ namespace DoFTools 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 @@ -1856,10 +1857,8 @@ namespace DoFTools // 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) @@ -1873,6 +1872,9 @@ namespace DoFTools } + // Implementation of the low level interface: + + template void make_periodicity_constraints (const FaceIterator &face_1, @@ -1881,7 +1883,9 @@ namespace DoFTools const ComponentMask &component_mask, const bool face_orientation, const bool face_flip, - const bool face_rotation) + const bool face_rotation, + const FullMatrix &matrix, + const std::vector &first_vector_components) { static const int dim = FaceIterator::AccessorType::dimension; @@ -1907,6 +1911,24 @@ namespace DoFTools 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: @@ -1967,30 +1989,156 @@ namespace DoFTools 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 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 &fe1 + = face_1->get_fe(face_1->nth_active_fe_index(0)); + + Quadrature quadrature (fe1.get_unit_face_support_points()); + + FEFaceValues 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 DoFTuple; + + FullMatrix rot_transformation + = IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face); + + for (unsigned int i=0; i::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= + 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; ihas_children() == false) - set_periodicity_constraints(face_2, face_1, - FullMatrix(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 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(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 + void + make_periodicity_constraints + (const std::vector > + &periodic_faces, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask) + { + typedef std::vector > + 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 void make_periodicity_constraints (const DH &dof_handler, @@ -2093,47 +2241,6 @@ namespace DoFTools - template - void - make_periodicity_constraints - (const std::vector > - &periodic_faces, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask) - { - typedef std::vector > - 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 diff --git a/source/dofs/dof_tools_constraints.inst.in b/source/dofs/dof_tools_constraints.inst.in index 4b2e4e8e21..3017a6d0b0 100644 --- a/source/dofs/dof_tools_constraints.inst.in +++ b/source/dofs/dof_tools_constraints.inst.in @@ -30,7 +30,17 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) const DH::face_iterator &, dealii::ConstraintMatrix &, const ComponentMask &, - bool, bool, bool); + bool, bool, bool, + const FullMatrix &, + const std::vector &); + + template + void + DoFTools::make_periodicity_constraints + (const std::vector > &, + dealii::ConstraintMatrix &, + const ComponentMask &); + template void @@ -67,14 +77,6 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) dealii::Tensor<1,DH::space_dimension> &, dealii::ConstraintMatrix &, const ComponentMask &); - - template - void - DoFTools::make_periodicity_constraints - (const std::vector > &, - dealii::ConstraintMatrix &, - const ComponentMask &); - #endif } diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 47570805b7..72a48e641d 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2367,14 +2367,15 @@ next_cell: * 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 */ template inline bool orthogonal_equality (const dealii::Point &point1, const dealii::Point &point2, const int direction, - const dealii::Tensor<1,spacedim> &offset) + const dealii::Tensor<1,spacedim> &offset, + const FullMatrix &matrix) { Assert (0<=direction && direction 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; } @@ -2475,7 +2482,8 @@ next_cell: 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 &matrix) { static const int dim = FaceIterator::AccessorType::dimension; @@ -2495,7 +2503,7 @@ next_cell: it++) { if (orthogonal_equality(face1->vertex(i),face2->vertex(*it), - direction, offset)) + direction, offset, matrix)) { matching[i] = *it; face2_vertices.erase(it); @@ -2518,11 +2526,12 @@ next_cell: 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 &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); } @@ -2535,9 +2544,11 @@ next_cell: match_periodic_face_pairs (std::set > &pairs1, std::set::type, unsigned int> > &pairs2, - const int direction, - std::vector > &matched_pairs, - const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset) + const int direction, + std::vector > &matched_pairs, + const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset, + const FullMatrix &matrix, + const std::vector &first_vector_components) { static const int space_dim = CellIterator::AccessorType::space_dimension; Assert (0<=direction && directionface(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 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; @@ -2588,12 +2600,14 @@ next_cell: template 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 > &matched_pairs, - const dealii::Tensor<1,CONTAINER::space_dimension> &offset) + const Tensor<1,CONTAINER::space_dimension> &offset, + const FullMatrix &matrix, + const std::vector &first_vector_components) { static const int dim = CONTAINER::dimension; static const int space_dim = CONTAINER::space_dimension; @@ -2632,7 +2646,8 @@ next_cell: 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); } @@ -2640,11 +2655,13 @@ next_cell: template 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 > &matched_pairs, - const dealii::Tensor<1,CONTAINER::space_dimension> &offset) + const Tensor<1,CONTAINER::space_dimension> &offset, + const FullMatrix &matrix, + const std::vector &first_vector_components) { static const int dim = CONTAINER::dimension; static const int space_dim = CONTAINER::space_dimension; @@ -2690,7 +2707,8 @@ next_cell: #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 diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index b949d29f2d..1c2ee8d9ec 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -302,26 +302,30 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II 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 &); template bool orthogonal_equality (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 &); template bool orthogonal_equality (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 &); template bool orthogonal_equality (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 &); #if deal_II_dimension >= 2 @@ -331,14 +335,18 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const types::boundary_id, const int, std::vector > &, - const Tensor<1,X::space_dimension> &); + const Tensor<1,X::space_dimension> &, + const FullMatrix &, + const std::vector &); template void collect_periodic_faces (const X &, const types::boundary_id, const int, std::vector > &, - const Tensor<1,X::space_dimension> &); + const Tensor<1,X::space_dimension> &, + const FullMatrix &, + const std::vector &); #endif @@ -352,7 +360,6 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS #if deal_II_dimension >= 2 namespace GridTools \{ - template void collect_periodic_faces > @@ -361,7 +368,9 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS const types::boundary_id, const int, std::vector::cell_iterator> > &, - const Tensor<1,parallel::distributed::Triangulation::space_dimension> &); + const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, + const FullMatrix &, + const std::vector &); template void @@ -370,8 +379,9 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS const types::boundary_id, const int, std::vector::cell_iterator> > &, - const Tensor<1,parallel::distributed::Triangulation::space_dimension> &); - + const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, + const FullMatrix &, + const std::vector &); \} #endif #endif -- 2.39.5