From baa7a73698266bf7cc9a63778464c2b753637bf2 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 26 Nov 2014 23:38:17 +0100 Subject: [PATCH] restructure and clean up some code --- source/dofs/dof_tools_constraints.cc | 83 +++++++++++++++++++--------- 1 file changed, 56 insertions(+), 27 deletions(-) diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index cf52dcac5e..1925471b5c 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -1909,24 +1909,40 @@ namespace DoFTools "on the very same face")); Assert(face_1->at_boundary() && face_2->at_boundary(), - ExcMessage ("Faces for periodicity constraints must be on the 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")); + Assert(matrix.m() == matrix.n(), + ExcMessage("The supplied (rotation or interpolation) matrix must " + "be a square matrix")); + + Assert(first_vector_components.empty() || matrix.m() == (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")); + 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(matrix.m() == 0 || matrix.m() == n_dofs || + matrix.m() == (int)dim, + ExcMessage ("matrix must have either size 0 or dim or the size " + "must be equal to the number of DoFs on the face")); + } + + if (!face_2->has_children()) + { + Assert(face_2->n_active_fe_indices() == 1, ExcInternalError()); + const unsigned int n_dofs = + face_2->get_fe(face_2->nth_active_fe_index(0)).dofs_per_face; + + Assert(matrix.m() == 0 || matrix.m() == n_dofs || + matrix.m() == (int)dim, + ExcMessage ("matrix must have either size 0 or dim or the size " + "must be equal to the number of DoFs on the face")); } #endif @@ -1995,16 +2011,25 @@ namespace DoFTools } } else - // otherwise at least one of the two faces is active and - // we need to enter the constraints { + // 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) + 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 @@ -2022,8 +2047,8 @@ namespace DoFTools // 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); + // start with a pristine interpolation matrix... + transformation = IdentityMatrix(n_dofs); for (unsigned int i=0; ihas_children() == false) { - FullMatrix inverse(transformation.size(0)); + FullMatrix inverse(transformation.m()); inverse.invert(transformation); set_periodicity_constraints(face_2, face_1, inverse, -- 2.39.5