From 8e6ae92e121d84c10782a34d5411789120da2380 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 26 Nov 2014 23:38:39 +0100 Subject: [PATCH] Code cleanup and several bugfixes Refactor the creation of the final (possibly rotated) interpolation matrix in make_periodicity_constraints into its own function. Add some additional asserts to catch corner cases. Return early if nothing to do (fixes dof_tools_21* tests) --- source/dofs/dof_tools_constraints.cc | 230 ++++++++++++++++----------- 1 file changed, 133 insertions(+), 97 deletions(-) diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 1925471b5c..489c4f150d 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -1643,6 +1643,8 @@ namespace DoFTools namespace { + // Internally used in make_periodicity_constraints. + // // enter constraints for periodicity into the given ConstraintMatrix object. // this function is called when at least one of the two face iterators corresponds // to an active object without further children @@ -1869,10 +1871,101 @@ namespace DoFTools } } } - } - // Implementation of the low level interface: + // Internally used in make_periodicity_constraints. + // + // Build up a (possibly rotated) interpolation matrix that is used in + // set_periodicity_constraints with the help of user supplied matrix + // and first_vector_components. + template + FullMatrix compute_transformation( + const FiniteElement &fe, + const FullMatrix &matrix, + const std::vector &first_vector_components) + { + Assert(matrix.m() == matrix.n(), ExcInternalError()); + + const unsigned int n_dofs = fe.dofs_per_face; + + if (matrix.m() == n_dofs) + { + // In case of m == n == n_dofs the supplied matrix is already + // an interpolation matrix, so we use it directly: + return matrix; + } + + if (first_vector_components.empty() && matrix.m() == 0) + { + // Just the identity matrix in case no rotation is specified: + return IdentityMatrix(n_dofs); + } + + // 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 + + Assert(matrix.m() == (int)dim, ExcInternalError()) + + Quadrature quadrature (fe.get_unit_face_support_points()); + FEFaceValues fe_face_values (fe, 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; + + // start with a pristine interpolation matrix... + FullMatrix transformation = IdentityMatrix(n_dofs); + + for (unsigned int i=0; i < n_dofs; ++i) + { + std::vector::const_iterator comp_it + = std::find (first_vector_components.begin(), + first_vector_components.end(), + fe.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 <= fe.n_components(), + ExcMessage("Error: the finite element does not have enough components " + "to define rotated periodic boundaries.")); + + for (unsigned int k=0; k < n_dofs; ++k) + if ((k != i) + && + (quadrature.point(k) == quadrature.point(i)) + && + (fe.face_system_to_component_index(k).first >= + first_vector_component) + && + (fe.face_system_to_component_index(k).first < + first_vector_component + dim)) + vector_dofs[fe.face_system_to_component_index(k).first - + first_vector_component] + = k; + + // ... and rotate all dofs belonging to vector valued + // components that are selected by first_vector_components: + for (int i=0; i @@ -1886,6 +1979,7 @@ namespace DoFTools const bool face_rotation, const FullMatrix &matrix, const std::vector &first_vector_components) + { static const int dim = FaceIterator::AccessorType::dimension; @@ -1975,16 +2069,18 @@ namespace DoFTools }, }; - // 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::max_children_per_face && - face_2->n_children() == GeometryInfo::max_children_per_face, - ExcNotImplemented()); + // In the case that both faces have children, we loop over all + // children and apply make_periodicty_constrains recursively: - for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - { + Assert(face_1->n_children() == GeometryInfo::max_children_per_face && + face_2->n_children() == + GeometryInfo::max_children_per_face, + ExcNotImplemented()); + + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; + ++i) { // Lookup the index for the second face unsigned int j; switch (dim) @@ -2012,111 +2108,51 @@ namespace DoFTools } 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.m() == n_dofs) - { - // In case of m == n == n_dofs the supplied matrix is already - // an interpolation matrix, so we use it directly: - Assert(matrix.n() == n_dofs, ExcInternalError()); - transformation = matrix; - } - else if (!first_vector_components.empty()) - { - // 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 + // Otherwise at least one of the two faces is active and + // we need to do some work and enter the constraints! - const FiniteElement &fe1 - = face_1->get_fe(face_1->nth_active_fe_index(0)); + // The finite element that matters is the one on the active face: + const FiniteElement &fe = + face_1->has_children() + ? face_2->get_fe(face_2->nth_active_fe_index(0)) + : face_1->get_fe(face_1->nth_active_fe_index(0)); - Quadrature quadrature (fe1.get_unit_face_support_points()); + const unsigned int n_dofs = fe.dofs_per_face; - FEFaceValues fe_face_values - (fe1, quadrature, update_q_points); + // Sometimes we just have nothing to do (for all finite elements, + // or systems which accidentally don't have any dofs on the + // boundary). + if (n_dofs == 0) + return; - // have an array that stores the location of each vector-dof tuple - // we want to rotate. - typedef std_cxx1x::array DoFTuple; + const FullMatrix transformation = + compute_transformation(fe, matrix, first_vector_components); - // start with a pristine interpolation matrix... - transformation = IdentityMatrix(n_dofs); - - 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; - - // ... and rotate all dofs belonging to vector valued - // components that are selected by first_vector_components: - for (int i=0; ihas_children() == false) + if (!face_2->has_children()) { FullMatrix inverse(transformation.m()); inverse.invert(transformation); - set_periodicity_constraints(face_2, face_1, + set_periodicity_constraints(face_2, + face_1, inverse, constraint_matrix, component_mask, - face_orientation, face_flip, face_rotation); + face_orientation, + face_flip, + face_rotation); } else { - set_periodicity_constraints(face_1, face_2, + Assert(!face_1->has_children(), ExcInternalError()); + + set_periodicity_constraints(face_1, + face_2, transformation, constraint_matrix, component_mask, - face_orientation, face_flip, face_rotation); + face_orientation, + face_flip, + face_rotation); } } } -- 2.39.5