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.

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 8396ea01f40e1b65cc71eb9e4d372ee1e94954af..d304ffbbc316d01ff0d4b5f8e37d54de4cb7f4a2 100644 (file)
--- 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 <tt>v_1</tt> and <tt>v_2</tt> are considered
-   * equal, if <code> (v_1 + offset) - v_2</code> is parallel to the unit
-   * vector in @p direction.
+   * relation.
+   *
+   * Hereby, two vertices <tt>v_1</tt> and <tt>v_2</tt> are considered
+   * equal, if $M\cdot v_1 + offset - v_2</code> 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<double> &matrix = FullMatrix<double>());
 
 
   /**
@@ -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<double> &matrix = FullMatrix<double>());
 
 
   /**
@@ -1240,7 +1248,7 @@ namespace GridTools
    const int                                                          direction,
    std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
    const Tensor<1,CONTAINER::space_dimension>                        &offset = dealii::Tensor<1,CONTAINER::space_dimension>(),
-   const FullMatrix<double>                                          &matrix = FullMatrix<double>(IdentityMatrix(CONTAINER::space_dimension)),
+   const FullMatrix<double>                                          &matrix = FullMatrix<double>(),
    const std::vector<unsigned int>                                   &first_vector_components = std::vector<unsigned int>());
 
 

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 8b5ad5d775dd7a12765821dd910edc610d1aec35..0a5a0c8322f86a2e7b4baffaa3dd47b9cc4f711e 100644 (file)
--- 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 <direction>
    */
   template<int spacedim>
-  inline bool orthogonal_equality (const dealii::Point<spacedim> &point1,
-                                   const dealii::Point<spacedim> &point2,
-                                   const int                     direction,
-                                   const dealii::Tensor<1,spacedim> &offset,
-                                   const FullMatrix<double>      &matrix)
+  inline bool orthogonal_equality (const Point<spacedim>    &point1,
+                                   const Point<spacedim>    &point2,
+                                   const int                 direction,
+                                   const Tensor<1,spacedim> &offset,
+                                   const FullMatrix<double> &matrix)
   {
     Assert (0<=direction && direction<spacedim,
             ExcIndexRange (direction, 0, spacedim));
+
+    Assert(matrix.m() == matrix.n(), ExcInternalError());
+
+    Point<spacedim> 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<double> &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<double> &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<CellIterator> matched_face
-                = {{cell1, cell2},{face_idx1, face_idx2}, orientation, matrix, first_vector_components};
+              const PeriodicFacePair<CellIterator> matched_face = {
+                  {cell1, cell2},
+                  {face_idx1, face_idx2},
+                  orientation,
+                  matrix,
+                  first_vector_components};
                 matched_pairs.push_back(matched_face);
                 pairs2.erase(it2);
                 ++n_matches;