From: bangerth Date: Thu, 16 May 2013 22:10:32 +0000 (+0000) Subject: Refactor a bit, in preparation to dealing with periodic faces not on X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6a2a1c654c61ae286cb3212797fd08333d061d5d;p=dealii-svn.git Refactor a bit, in preparation to dealing with periodic faces not on the same refinement level. git-svn-id: https://svn.dealii.org/trunk@29521 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index 45a2e0406f..675ba22fd0 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -2786,7 +2786,120 @@ namespace DoFTools - template + namespace + { + // 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 + // + // @param transformation A matrix that maps degrees of freedom from one face + // to another. If the DoFs on the two faces are supposed to match exactly, then + // the matrix so provided will be the identity matrix. if face 2 is once refined + // from face 1, then the matrix needs to be the interpolation matrix from a face + // to this particular child + template + void + set_periodicity_constraints (const FaceIterator &face_1, + const typename identity::type &face_2, + const FullMatrix &transformation, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask, + const bool face_orientation, + const bool face_flip, + const bool face_rotation) + { + static const int dim = FaceIterator::AccessorType::dimension; + static const int spacedim = FaceIterator::AccessorType::space_dimension; + + // 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 ... + 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")); + + const FiniteElement &fe = face_1->get_fe(face_1_index); + + Assert (component_mask.represents_n_components(fe.n_components()), + ExcMessage ("The number of components in the mask has to be either " + "zero or equal to the number of components in the finite " + "element.")); + + const unsigned int dofs_per_face = fe.dofs_per_face; + + std::vector dofs_1(dofs_per_face); + std::vector dofs_2(dofs_per_face); + + face_1->get_dof_indices(dofs_1, face_1_index); + face_2->get_dof_indices(dofs_2, face_2_index); + + + // 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 directly. But luckily we can match them back :-] + + std::map 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]) + { + if (!constraint_matrix.is_constrained(dofs_1[i])) + { + constraint_matrix.add_line(dofs_1[i]); + constraint_matrix.add_entry(dofs_1[i], dofs_2[j], 1.0); + } + } + } + } + } + + + template void make_periodicity_constraints (const FaceIterator &face_1, const typename identity::type &face_2, @@ -2882,93 +2995,15 @@ namespace DoFTools face_flip, face_rotation); } - return; - } - - // .. 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 .. - 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")); - - const dealii::FiniteElement &fe = face_1->get_fe(face_1_index); - - Assert (component_mask.represents_n_components(fe.n_components()), - ExcMessage ("The number of components in the mask has to be either " - "zero or equal to the number of components in the finite " - "element.")); - - const unsigned int dofs_per_face = fe.dofs_per_face; - - std::vector dofs_1(dofs_per_face); - std::vector dofs_2(dofs_per_face); - - face_1->get_dof_indices(dofs_1, face_1_index); - face_2->get_dof_indices(dofs_2, face_2_index); - - - // 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 directly. But luckily we can match them back :-] - - std::map 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]) - { - if (!constraint_matrix.is_constrained(dofs_1[i])) - { - constraint_matrix.add_line(dofs_1[i]); - constraint_matrix.add_entry(dofs_1[i], dofs_2[j], 1.0); - } - } } + else + // otherwise at least one of the two faces is active and + // we need to enter the constraints + 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); }