From 48d8293bf02c7de5a925826580a08c73c0d371b5 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 26 Nov 2014 22:53:04 +0100 Subject: [PATCH] Bugfix: Conditionally use matrix in orthogonal_equality Only apply the parameter matrix in orthogonal_equality if it is a spacedim x spacedim matrix and can be interpreted as a rotation. Further bugfixes and documentation updates. --- include/deal.II/grid/grid_tools.h | 20 ++++++++++---- source/grid/grid_tools.cc | 46 +++++++++++++++++++++---------- 2 files changed, 45 insertions(+), 21 deletions(-) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 8396ea01f4..d304ffbbc3 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1102,9 +1102,13 @@ namespace GridTools * * @p face1 and @p face2 are considered equal, if a one to one matching * between its vertices can be achieved via an orthogonal equality - * relation: Two vertices v_1 and v_2 are considered - * equal, if (v_1 + offset) - v_2 is parallel to the unit - * vector in @p direction. + * relation. + * + * Hereby, two vertices v_1 and v_2 are considered + * equal, if $M\cdot v_1 + offset - v_2 is parallel to the unit + * vector in unit direction @p direction. If the parameter @p matrix + * is a reference to a spacedim x spacedim matrix, $M$ is set to @p + * matrix, otherwise $M$ is the identity matrix. * * If the matching was successful, the _relative_ orientation of @p face1 * with respect to @p face2 is returned in the bitset @p orientation, @@ -1166,7 +1170,9 @@ namespace GridTools const FaceIterator &face1, const FaceIterator &face2, const int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset); + const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset + = Tensor<1,FaceIterator::AccessorType::space_dimension>(), + const FullMatrix &matrix = FullMatrix()); /** @@ -1177,7 +1183,9 @@ namespace GridTools orthogonal_equality (const FaceIterator &face1, const FaceIterator &face2, const int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset); + const Tensor<2,FaceIterator::AccessorType::space_dimension> &offset + = Tensor<1,FaceIterator::AccessorType::space_dimension>(), + const FullMatrix &matrix = FullMatrix()); /** @@ -1240,7 +1248,7 @@ namespace GridTools const int direction, 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 FullMatrix &matrix = FullMatrix(), const std::vector &first_vector_components = std::vector()); diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 8b5ad5d775..0a5a0c8322 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2372,26 +2372,35 @@ next_cell: * 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 FullMatrix &matrix) + inline bool orthogonal_equality (const Point &point1, + const Point &point2, + const int direction, + const Tensor<1,spacedim> &offset, + const FullMatrix &matrix) { Assert (0<=direction && direction distance; + + if (matrix.m() == spacedim) + for (int i = 0; i < spacedim; ++i) + for (int j = 0; j < spacedim; ++j) + distance(i) = matrix(i,j) * point1(j); + else + distance = point1; + + distance += offset - point2; + for (int i = 0; i < spacedim; ++i) { // Only compare coordinate-components != direction: if (i == direction) continue; - 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) + if (fabs(distance(i)) > 1.e-10) return false; } @@ -2483,9 +2492,12 @@ next_cell: const FaceIterator &face1, const FaceIterator &face2, const int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, + const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, const FullMatrix &matrix) { + Assert(matrix.m() == matrix.n(), + ExcMessage("The supplied matrix must be a square matrix")); + static const int dim = FaceIterator::AccessorType::dimension; // Do a full matching of the face vertices: @@ -2527,7 +2539,7 @@ next_cell: orthogonal_equality (const FaceIterator &face1, const FaceIterator &face2, const int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, + const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, const FullMatrix &matrix) { // Call the function above with a dummy orientation array @@ -2581,8 +2593,12 @@ next_cell: // 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, matrix, first_vector_components}; + const PeriodicFacePair matched_face = { + {cell1, cell2}, + {face_idx1, face_idx2}, + orientation, + matrix, + first_vector_components}; matched_pairs.push_back(matched_face); pairs2.erase(it2); ++n_matches; -- 2.39.5