From faeded4b15f3bca73802c019d6cf1aa5ddb33e94 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Tue, 23 Jan 2018 01:27:01 +0100
Subject: [PATCH] Fix periodicity constraints for multiply constrained DoFs

---
 source/dofs/dof_tools_constraints.cc | 343 ++++++++++++++++-----------
 tests/bits/periodicity_06.cc         | 218 +++++++++++++++++
 tests/bits/periodicity_06.output     |   9 +
 tests/bits/periodicity_07.cc         | 249 +++++++++++++++++++
 tests/bits/periodicity_07.output     |   5 +
 5 files changed, 682 insertions(+), 142 deletions(-)
 create mode 100644 tests/bits/periodicity_06.cc
 create mode 100644 tests/bits/periodicity_06.output
 create mode 100644 tests/bits/periodicity_07.cc
 create mode 100644 tests/bits/periodicity_07.output

diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc
index 00718d0672..4b2fe3cce1 100644
--- a/source/dofs/dof_tools_constraints.cc
+++ b/source/dofs/dof_tools_constraints.cc
@@ -1873,167 +1873,226 @@ namespace DoFTools
 
           // loop over all dofs on face 2 and constrain them against the ones on face 1
           for (unsigned int i=0; i<dofs_per_face; ++i)
-            if (!constraint_matrix.is_constrained(dofs_2[i]))
-              if ((component_mask.n_selected_components(fe.n_components())
-                   == fe.n_components())
-                  ||
-                  component_mask[fe.face_system_to_component_index(i).first])
-                {
-                  // as mentioned in the comment above this function, we need
-                  // to be careful about treating identity constraints differently.
-                  // consequently, find out whether this dof 'i' will be
-                  // identity constrained
-                  //
-                  // to check whether this is the case, first see whether there are
-                  // any weights other than 0 and 1, then in a first stage make sure
-                  // that if so there is only one weight equal to 1
-                  //
-                  // afterwards do the same for constraints of type dof1=-dof2
-                  bool is_identity_constrained = true;
-                  const double eps = 1.e-13;
+            if ((component_mask.n_selected_components(fe.n_components())
+                 == fe.n_components())
+                ||
+                component_mask[fe.face_system_to_component_index(i).first])
+              {
+                // as mentioned in the comment above this function, we need
+                // to be careful about treating identity constraints differently.
+                // consequently, find out whether this dof 'i' will be
+                // identity constrained
+                //
+                // to check whether this is the case, first see whether there are
+                // any weights other than 0 and 1, then in a first stage make sure
+                // that if so there is only one weight equal to 1
+                //
+                // afterwards do the same for constraints of type dof1=-dof2
+                bool is_identity_constrained = true;
+                const double eps = 1.e-13;
+                for (unsigned int jj=0; jj<dofs_per_face; ++jj)
+                  if (((std::abs(transformation(i,jj)) < eps) ||
+                       (std::abs(transformation(i,jj)-1) < eps)) == false)
+                    {
+                      is_identity_constrained = false;
+                      break;
+                    }
+                unsigned int identity_constraint_target = numbers::invalid_unsigned_int;
+                if (is_identity_constrained == true)
+                  {
+                    bool one_identity_found = false;
+                    for (unsigned int jj=0; jj<dofs_per_face; ++jj)
+                      if (std::abs(transformation(i,jj)-1.) < eps)
+                        {
+                          if (one_identity_found == false)
+                            {
+                              one_identity_found = true;
+                              identity_constraint_target = jj;
+                            }
+                          else
+                            {
+                              is_identity_constrained = false;
+                              identity_constraint_target = numbers::invalid_unsigned_int;
+                              break;
+                            }
+                        }
+                  }
+
+                bool is_inverse_constrained = !is_identity_constrained;
+                unsigned int inverse_constraint_target = numbers::invalid_unsigned_int;
+                if (is_inverse_constrained)
                   for (unsigned int jj=0; jj<dofs_per_face; ++jj)
                     if (((std::abs(transformation(i,jj)) < eps) ||
-                         (std::abs(transformation(i,jj)-1) < eps)) == false)
+                         (std::abs(transformation(i,jj)+1) < eps)) == false)
                       {
-                        is_identity_constrained = false;
+                        is_inverse_constrained = false;
                         break;
                       }
-                  unsigned int identity_constraint_target = numbers::invalid_unsigned_int;
-                  if (is_identity_constrained == true)
-                    {
-                      bool one_identity_found = false;
-                      for (unsigned int jj=0; jj<dofs_per_face; ++jj)
-                        if (std::abs(transformation(i,jj)-1.) < eps)
+                if (is_inverse_constrained)
+                  {
+                    bool one_identity_found = false;
+                    for (unsigned int jj=0; jj<dofs_per_face; ++jj)
+                      if (std::abs(transformation(i,jj)+1) < eps)
+                        {
+                          if (one_identity_found == false)
+                            {
+                              one_identity_found = true;
+                              inverse_constraint_target = jj;
+                            }
+                          else
+                            {
+                              is_inverse_constrained = false;
+                              inverse_constraint_target = numbers::invalid_unsigned_int;
+                              break;
+                            }
+                        }
+                  }
+
+                const unsigned int target = is_identity_constrained
+                                            ? identity_constraint_target
+                                            : inverse_constraint_target;
+
+                // find out whether this dof also exists on face 1
+                // if this is true and the constraint is no identity
+                // constraint to itself, set it to zero
+                bool constraint_set = false;
+                for (unsigned int j=0; j<dofs_per_face; ++j)
+                  {
+                    if (dofs_2[i] == dofs_1[j])
+                      if (!(is_identity_constrained && target==i))
+                        {
+                          constraint_matrix.add_line(dofs_2[i]);
+                          constraint_set = true;
+                        }
+                  }
+
+                if (!constraint_set)
+                  {
+                    // now treat constraints, either as an equality constraint or
+                    // as a sequence of constraints
+                    if (is_identity_constrained == true || is_inverse_constrained == true)
+                      {
+                        // Query the correct face_index on face_1 respecting the given
+                        // orientation:
+                        const unsigned int j
+                          = cell_to_rotated_face_index[fe.face_to_cell_index(target,
+                                                                             0, /* It doesn't really matter, just assume
+                                                                                   * we're on the first face...
+                                                                                   */
+                                                                             face_orientation, face_flip, face_rotation)];
+
+                        if (constraint_matrix.is_constrained(dofs_2[i]))
                           {
-                            if (one_identity_found == false)
+                            // if the two aren't already identity constrained (whichever way
+                            // around) or already identical (in case of rotated periodicity constraints),
+                            // then enter the constraint. otherwise there is nothing for us still to do
+                            bool enter_constraint = false;
+                            // see if this would add an identity constraint cycle
+                            if (!constraint_matrix.is_constrained(dofs_1[j]))
                               {
-                                one_identity_found = true;
-                                identity_constraint_target = jj;
+                                types::global_dof_index new_dof = dofs_2[i];
+                                while (new_dof != dofs_1[j])
+                                  if (constraint_matrix.is_constrained(new_dof))
+                                    {
+                                      const std::vector<std::pair<types::global_dof_index, double > > *constraint_entries
+                                        = constraint_matrix.get_constraint_entries(new_dof);
+                                      if (constraint_entries->size()==1)
+                                        new_dof = (*constraint_entries)[0].first;
+                                      else
+                                        {
+                                          enter_constraint = true;
+                                          break;
+                                        }
+                                    }
+                                  else
+                                    {
+                                      enter_constraint = true;
+                                      break;
+                                    }
                               }
-                            else
+
+                            if (enter_constraint)
                               {
-                                is_identity_constrained = false;
-                                identity_constraint_target = numbers::invalid_unsigned_int;
-                                break;
+                                constraint_matrix.add_line(dofs_1[j]);
+                                constraint_matrix.add_entry(dofs_1[j], dofs_2[i], is_identity_constrained?1.0:-1.0);
                               }
                           }
-                    }
-
-                  bool is_inverse_constrained = !is_identity_constrained;
-                  unsigned int inverse_constraint_target = numbers::invalid_unsigned_int;
-                  if (is_inverse_constrained)
-                    for (unsigned int jj=0; jj<dofs_per_face; ++jj)
-                      if (((std::abs(transformation(i,jj)) < eps) ||
-                           (std::abs(transformation(i,jj)+1) < eps)) == false)
-                        {
-                          is_inverse_constrained = false;
-                          break;
-                        }
-                  if (is_inverse_constrained)
-                    {
-                      bool one_identity_found = false;
-                      for (unsigned int jj=0; jj<dofs_per_face; ++jj)
-                        if (std::abs(transformation(i,jj)+1) < eps)
+                        else
                           {
-                            if (one_identity_found == false)
+                            // if the two aren't already identity constrained (whichever way
+                            // around) or already identical (in case of rotated periodicity constraints),
+                            // then enter the constraint. otherwise there is nothing for us still to do
+                            bool enter_constraint = false;
+                            if (!constraint_matrix.is_constrained(dofs_1[j]))
                               {
-                                one_identity_found = true;
-                                inverse_constraint_target = jj;
+                                if (dofs_2[i] != dofs_1[j])
+                                  enter_constraint = true;
                               }
-                            else
+                            else //dofs_1[j] is constrained, is it identity or inverse constrained?
                               {
-                                is_inverse_constrained = false;
-                                inverse_constraint_target = numbers::invalid_unsigned_int;
-                                break;
+                                const std::vector<std::pair<types::global_dof_index, double > > *constraint_entries
+                                  = constraint_matrix.get_constraint_entries(dofs_1[j]);
+                                if (constraint_entries->size()==1 && (*constraint_entries)[0].first == dofs_2[i])
+                                  {
+                                    if ((is_identity_constrained && std::abs((*constraint_entries)[0].second-1) > eps) ||
+                                        (is_inverse_constrained && std::abs((*constraint_entries)[0].second+1) > eps))
+                                      {
+                                        //this pair of constraints means that both dofs have to be constrained to 0.
+                                        constraint_matrix.add_line(dofs_2[i]);
+                                      }
+                                  }
+                                else
+                                  {
+                                    // see if this would add an identity constraint cycle
+                                    types::global_dof_index new_dof = dofs_1[j];
+                                    while (new_dof != dofs_2[i])
+                                      if (constraint_matrix.is_constrained(new_dof))
+                                        {
+                                          const std::vector<std::pair<types::global_dof_index, double > > *constraint_entries
+                                            = constraint_matrix.get_constraint_entries(new_dof);
+                                          if (constraint_entries->size()==1)
+                                            new_dof = (*constraint_entries)[0].first;
+                                          else
+                                            {
+                                              enter_constraint = true;
+                                              break;
+                                            }
+                                        }
+                                      else
+                                        {
+                                          enter_constraint = true;
+                                          break;
+                                        }
+                                  }
                               }
-                          }
-                    }
 
-                  const unsigned int target = is_identity_constrained
-                                              ? identity_constraint_target
-                                              : inverse_constraint_target;
-
-                  // find out whether this dof also exists on face 1
-                  // if this is true and the constraint is no identity
-                  // constraint to itself, set it to zero
-                  bool constrained_set = false;
-                  for (unsigned int j=0; j<dofs_per_face; ++j)
-                    {
-                      if (dofs_2[i] == dofs_1[j])
-                        if (!(is_identity_constrained && target==i))
+                            if (enter_constraint)
+                              {
+                                constraint_matrix.add_line(dofs_2[i]);
+                                constraint_matrix.add_entry(dofs_2[i], dofs_1[j], is_identity_constrained?1.0:-1.0);
+                              }
+                          }
+                      }
+                    else if (!constraint_matrix.is_constrained(dofs_2[i]))
+                      {
+                        // this is just a regular constraint. enter it piece by piece
+                        constraint_matrix.add_line(dofs_2[i]);
+                        for (unsigned int jj=0; jj<dofs_per_face; ++jj)
                           {
-                            constraint_matrix.add_line(dofs_2[i]);
-                            constrained_set = true;
+                            // Query the correct face_index on face_1 respecting the given
+                            // orientation:
+                            const unsigned int j =
+                              cell_to_rotated_face_index[fe.face_to_cell_index
+                                                         (jj, 0, face_orientation, face_flip, face_rotation)];
+
+                            // And finally constrain the two DoFs respecting component_mask:
+                            if (transformation(i,jj) != 0)
+                              constraint_matrix.add_entry(dofs_2[i], dofs_1[j],
+                                                          transformation(i,jj));
                           }
-                    }
-
-                  if (!constrained_set)
-                    {
-                      // now treat constraints, either as an equality constraint or
-                      // as a sequence of constraints
-                      if (is_identity_constrained == true || is_inverse_constrained == true)
-                        {
-                          // Query the correct face_index on face_1 respecting the given
-                          // orientation:
-                          const unsigned int j
-                            = cell_to_rotated_face_index[fe.face_to_cell_index(target,
-                                                                               0, /* It doesn't really matter, just assume
-                                                                                 * we're on the first face...
-                                                                                 */
-                                                                               face_orientation, face_flip, face_rotation)];
-
-                          // if the two aren't already identity constrained (whichever way
-                          // around) or already identical (in case of rotated periodicity constraints),
-                          // then enter the constraint. otherwise there is nothing for us still to do
-                          bool enter_constraint = false;
-                          if (!constraint_matrix.is_constrained(dofs_1[j]))
-                            {
-                              if (dofs_2[i] != dofs_1[j])
-                                enter_constraint = true;
-                            }
-                          else //dofs_1[j] is constrained, is it identity or inverse constrained?
-                            {
-                              const std::vector<std::pair<types::global_dof_index, double > > *constraint_entries
-                                = constraint_matrix.get_constraint_entries(dofs_1[j]);
-                              if (constraint_entries->size()==1 && (*constraint_entries)[0].first == dofs_2[i])
-                                {
-                                  if ((is_identity_constrained && std::abs((*constraint_entries)[0].second-1) > eps) ||
-                                      (is_inverse_constrained && std::abs((*constraint_entries)[0].second+1) > eps))
-                                    {
-                                      //this pair of constraints means that both dofs have to be constrained to 0.
-                                      constraint_matrix.add_line(dofs_2[i]);
-                                    }
-                                }
-                              else
-                                enter_constraint = true;
-                            }
-
-                          if (enter_constraint)
-                            {
-                              constraint_matrix.add_line(dofs_2[i]);
-                              constraint_matrix.add_entry(dofs_2[i], dofs_1[j], is_identity_constrained?1.0:-1.0);
-                            }
-                        }
-                      else
-                        {
-                          // this is just a regular constraint. enter it piece by piece
-                          constraint_matrix.add_line(dofs_2[i]);
-                          for (unsigned int jj=0; jj<dofs_per_face; ++jj)
-                            {
-                              // Query the correct face_index on face_1 respecting the given
-                              // orientation:
-                              const unsigned int j =
-                                cell_to_rotated_face_index[fe.face_to_cell_index
-                                                           (jj, 0, face_orientation, face_flip, face_rotation)];
-
-                              // And finally constrain the two DoFs respecting component_mask:
-                              if (transformation(i,jj) != 0)
-                                constraint_matrix.add_entry(dofs_2[i], dofs_1[j],
-                                                            transformation(i,jj));
-                            }
-                        }
-                    }
-                }
+                      }
+                  }
+              }
         }
     }
 
diff --git a/tests/bits/periodicity_06.cc b/tests/bits/periodicity_06.cc
new file mode 100644
index 0000000000..71ec4602e5
--- /dev/null
+++ b/tests/bits/periodicity_06.cc
@@ -0,0 +1,218 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Make sure that periodic boundary conditions also work correctly if we have
+// multiple periodic boundary pairs that meet at an edge.
+// This tests the 2D case.
+
+
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+ConstraintMatrix make_constraint_matrix(const DoFHandler<2> &dof_handler, int version)
+{
+  constexpr int dim = 2;
+  ConstraintMatrix constraints;
+  constraints.clear();
+  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+  std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> > periodicity_vectorDof;
+  switch (version)
+    {
+    case 0:
+      GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof);
+    case 1:
+      GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof);
+    case 2:
+      GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof);
+    case 3:
+      GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof);
+    case 4:
+      GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof);
+    case 5:
+      GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof);
+    case 6:
+      GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof);
+    case 7:
+      GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof);
+    }
+
+  DoFTools::make_periodicity_constraints<DoFHandler<dim> >(periodicity_vectorDof, constraints);
+
+  constraints.close();
+  /*  std::map<types::global_dof_index, Point<dim> > support_points;
+    DoFTools::map_dofs_to_support_points (MappingQ<dim,dim>(1), dof_handler, support_points);
+    for (const auto &line: constraints.get_lines())
+      for (const auto &entry: line.entries)
+        std::cout << "DoF " << line.index << " at " << support_points[line.index]
+                  << " is constrained to " << " DoF " << entry.first << " at "
+                  << support_points[entry.first]
+                  << " with value " << entry.second << std::endl;*/
+  return constraints;
+}
+
+template <int dim>
+class PeriodicReference : public Function<dim>
+{
+public:
+  PeriodicReference () : Function<dim>() {}
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const override
+  {
+    if (dim==3)
+      return std::sin(p(0)+1.)*std::sin(p(1)+2.)*std::sin(p(2)+3.);
+    return std::sin(p(0)+1.)*std::sin(p(1)+2.);
+  }
+};
+
+
+template <int dim>
+void get_point_value
+(const DoFHandler<dim> &dof_handler,
+ const Point<dim> &point,
+ const Vector<double> &solution,
+ Vector<double> &value)
+{
+  VectorTools::point_value (dof_handler, solution,
+                            point, value);
+}
+
+
+void check_periodicity(const DoFHandler<2> &dof_handler, Vector<double> &solution, const unsigned int cycle)
+{
+  unsigned int n_points = 2;
+  for (unsigned int i = 0; i<cycle; i++)
+    n_points*=2;
+
+  //don't test exactly at the support points, since point_value is not stable there
+  const double eps = 1./(16.*n_points);
+
+  for (unsigned int i=1; i< n_points; i++)
+    {
+      Vector<double> value1(1);
+      Vector<double> value2(1);
+
+      Point <2> point1;
+      point1(0)=-numbers::PI+2.*i/n_points+eps;
+      point1(1)=-numbers::PI;
+      Point <2> point2;
+      point2(0)=-numbers::PI+2.*i/n_points+eps;
+      point2(1)=numbers::PI;
+
+      VectorTools::point_value (dof_handler, solution, point1, value1);
+      VectorTools::point_value (dof_handler, solution, point2, value2);
+
+      if (std::abs(value2[0]-value1[0])>1e-8)
+        {
+          std::cout << point1 << "\t" << "fail" << std::endl;
+          std::cout<<point1<< "\t" << value1[0] << std::endl;
+          std::cout<<point2<< "\t" << value2[0] << std::endl;
+          AssertThrow(false, ExcInternalError());
+        }
+      else
+        {
+          std::cout << point1 << "\t" << "pass" << std::endl;
+        }
+    }
+  for (unsigned int i=1; i< n_points; i++)
+    {
+      Vector<double> value1(1);
+      Vector<double> value2(1);
+
+      Point <2> point1;
+      point1(1)=-numbers::PI+2.*i/n_points+eps;
+      point1(0)=-numbers::PI;
+      Point <2> point2;
+      point2(1)=-numbers::PI+2.*i/n_points+eps;
+      point2(0)=numbers::PI;
+
+      VectorTools::point_value (dof_handler, solution, point1, value1);
+      VectorTools::point_value (dof_handler, solution, point2, value2);
+
+      if (std::abs(value2[0]-value1[0])>1e-8)
+        {
+          std::cout << point1 << "\t" << "fail" << std::endl;
+          std::cout<<point1<< "\t" << value1[0] << std::endl;
+          std::cout<<point2<< "\t" << value2[0] << std::endl;
+          Assert(false, ExcInternalError());
+        }
+      else
+        {
+          std::cout << point1 << "\t" << "pass" << std::endl;
+        }
+    }
+}
+
+
+int main (int argc, char *argv[])
+{
+  initlog();
+
+  constexpr int dim = 2;
+  const double L=numbers::PI;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube (triangulation, -L, L, true);
+
+  triangulation.refine_global(1);
+  typename Triangulation<dim>::active_cell_iterator cellBegin = triangulation.begin_active();
+  cellBegin->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement();
+
+  FE_Q<dim> fe(1);
+  DoFHandler<dim> dof_handler (triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  std::vector<ConstraintMatrix> constraints(8);
+
+  PeriodicReference<dim> periodic_function;
+
+  std::vector<Vector<double> > projection(8, Vector<double> (dof_handler.n_dofs()));
+
+  for (unsigned int i=0; i<8; ++i)
+    {
+      deallog << "Testing version " << i << std::endl;
+      constraints[i] = make_constraint_matrix (dof_handler, i);
+      VectorTools::project (dof_handler, constraints[i], QGauss<dim>(3), periodic_function, projection[i]);
+      check_periodicity(dof_handler, projection[i], i);
+    }
+}
+
diff --git a/tests/bits/periodicity_06.output b/tests/bits/periodicity_06.output
new file mode 100644
index 0000000000..d6ba424002
--- /dev/null
+++ b/tests/bits/periodicity_06.output
@@ -0,0 +1,9 @@
+
+DEAL::Testing version 0
+DEAL::Testing version 1
+DEAL::Testing version 2
+DEAL::Testing version 3
+DEAL::Testing version 4
+DEAL::Testing version 5
+DEAL::Testing version 6
+DEAL::Testing version 7
diff --git a/tests/bits/periodicity_07.cc b/tests/bits/periodicity_07.cc
new file mode 100644
index 0000000000..0f93f9eff1
--- /dev/null
+++ b/tests/bits/periodicity_07.cc
@@ -0,0 +1,249 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Make sure that periodic boundary conditions also work correctly if we have
+// multiple periodic boundary pairs that meet at an edge.
+// This tests the 3D case.
+
+
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+ConstraintMatrix make_constraint_matrix(const DoFHandler<3> &dof_handler, int version)
+{
+  constexpr int dim = 3;
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+  std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> > periodicity_vectorDof;
+  switch (version)
+    {
+    case 0:
+      GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 4, 5, 2, periodicity_vectorDof);
+    case 1:
+      GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 5, 4, 2, periodicity_vectorDof);
+    case 2:
+      GridTools::collect_periodic_faces(dof_handler, 4, 5, 2, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof);
+    case 3:
+      GridTools::collect_periodic_faces(dof_handler, 5, 4, 2, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof);
+      GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof);
+    }
+
+  DoFTools::make_periodicity_constraints<DoFHandler<dim> >(periodicity_vectorDof, constraints);
+
+  constraints.close();
+  /*std::map<types::global_dof_index, Point<dim> > support_points;
+  DoFTools::map_dofs_to_support_points (MappingQ<dim,dim>(1), dof_handler, support_points);
+  for (const auto &line: constraints.get_lines())
+    for (const auto &entry: line.entries)
+      std::cout << "DoF " << line.index << " at " << support_points[line.index]
+                << " is constrained to " << " DoF " << entry.first << " at "
+                << support_points[entry.first]
+                << " with value " << entry.second << std::endl;*/
+  return constraints;
+}
+
+template <int dim>
+class PeriodicReference : public Function<dim>
+{
+public:
+  PeriodicReference () : Function<dim>() {}
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const override
+  {
+    if (dim==3)
+      return std::sin(p(0)+1.)*std::sin(p(1)+2.)*std::sin(p(2)+3.);
+    return std::sin(p(0)+1.)*std::sin(p(1)+2.);
+  }
+};
+
+
+template <int dim>
+void get_point_value
+(const DoFHandler<dim> &dof_handler,
+ const Point<dim> &point,
+ const Vector<double> &solution,
+ Vector<double> &value)
+{
+  VectorTools::point_value (dof_handler, solution,
+                            point, value);
+}
+
+
+void check_periodicity(const DoFHandler<3> &dof_handler, Vector<double> &solution, const unsigned int cycle)
+{
+  unsigned int n_points = 2;
+  for (unsigned int i = 0; i<cycle; ++i)
+    n_points*=2;
+
+  //don't test exactly at the support points, since point_value is not stable there
+  const double eps = 1./(16.*n_points);
+
+  for (unsigned int i=1; i<n_points; ++i)
+    for (unsigned int j=1; j<n_points; ++j)
+      {
+        Vector<double> value1(1);
+        Vector<double> value2(1);
+
+        Point<3> point1;
+        point1(0)=-numbers::PI+2.*i/n_points+eps;
+        point1(1)=-numbers::PI;
+        point1(2)=-numbers::PI+2.*j/n_points+eps;
+        Point<3> point2;
+        point2(0)=-numbers::PI+2.*i/n_points+eps;
+        point2(1)=numbers::PI;
+        point2(2)=-numbers::PI+2.*j/n_points+eps;
+
+        VectorTools::point_value (dof_handler, solution, point1, value1);
+        VectorTools::point_value (dof_handler, solution, point2, value2);
+
+        if (std::abs(value2[0]-value1[0])>1e-8)
+          {
+            std::cout << point1 << "\t" << "fail" << std::endl;
+            std::cout<<point1<< "\t" << value1[0] << std::endl;
+            std::cout<<point2<< "\t" << value2[0] << std::endl;
+            Assert(false, ExcInternalError());
+          }
+        else
+          {
+            std::cout << point1 << "\t" << "pass" << std::endl;
+          }
+      }
+
+  for (unsigned int i=1; i<n_points; ++i)
+    for (unsigned int j=1; j<n_points; ++j)
+      {
+        Vector<double> value1(1);
+        Vector<double> value2(1);
+
+        Point <3> point1;
+        point1(2)=-numbers::PI+2.*j/n_points+eps;
+        point1(1)=-numbers::PI+2.*i/n_points+eps;
+        point1(0)=-numbers::PI;
+        Point <3> point2;
+        point2(2)=-numbers::PI+2.*j/n_points+eps;
+        point2(1)=-numbers::PI+2.*i/n_points+eps;
+        point2(0)=numbers::PI;
+
+        VectorTools::point_value (dof_handler, solution, point1, value1);
+        VectorTools::point_value (dof_handler, solution, point2, value2);
+
+        if (std::abs(value2[0]-value1[0])>1e-8)
+          {
+            std::cout << point1 << "\t" << "fail" << std::endl;
+            std::cout<<point1<< "\t" << value1[0] << std::endl;
+            std::cout<<point2<< "\t" << value2[0] << std::endl;
+            Assert(false, ExcInternalError());
+          }
+        else
+          {
+            std::cout << point1 << "\t" << "pass" << std::endl;
+          }
+      }
+
+  for (unsigned int i=1; i<n_points; ++i)
+    for (unsigned int j=1; j<n_points; ++j)
+      {
+        Vector<double> value1(1);
+        Vector<double> value2(1);
+
+        Point <3> point1;
+        point1(0)=-numbers::PI+2.*j/n_points+eps;
+        point1(1)=-numbers::PI+2.*i/n_points+eps;
+        point1(2)=-numbers::PI;
+        Point <3> point2;
+        point2(0)=-numbers::PI+2.*j/n_points+eps;
+        point2(1)=-numbers::PI+2.*i/n_points+eps;
+        point2(2)=numbers::PI;
+
+        VectorTools::point_value (dof_handler, solution, point1, value1);
+        VectorTools::point_value (dof_handler, solution, point2, value2);
+
+        if (std::abs(value2[0]-value1[0])>1e-8)
+          {
+            std::cout << point1 << "\t" << "fail" << std::endl;
+            std::cout<<point1<< "\t" << value1[0] << std::endl;
+            std::cout<<point2<< "\t" << value2[0] << std::endl;
+            Assert(false, ExcInternalError());
+          }
+        else
+          {
+            std::cout << point1 << "\t" << "pass" << std::endl;
+          }
+      }
+}
+
+
+int main (int argc, char *argv[])
+{
+  initlog();
+
+  constexpr int dim = 3;
+  const double L=numbers::PI;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube (triangulation, -L, L, true);
+
+  triangulation.refine_global(1);
+  typename Triangulation<dim>::active_cell_iterator cellBegin = triangulation.begin_active();
+  cellBegin->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement();
+
+  FE_Q<dim> fe(1);
+  DoFHandler<dim> dof_handler (triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  std::vector<ConstraintMatrix> constraints(4);
+
+  PeriodicReference<dim> periodic_function;
+
+  std::vector<Vector<double> > projection(4, Vector<double> (dof_handler.n_dofs()));
+
+  for (unsigned int i=0; i<4; ++i)
+    {
+      deallog << "Testing version " << i << std::endl;
+      constraints[i] = make_constraint_matrix (dof_handler, i);
+      VectorTools::project (dof_handler, constraints[i], QGauss<dim>(3), periodic_function, projection[i]);
+      check_periodicity(dof_handler, projection[i], i);
+    }
+}
+
diff --git a/tests/bits/periodicity_07.output b/tests/bits/periodicity_07.output
new file mode 100644
index 0000000000..0132ef2ce9
--- /dev/null
+++ b/tests/bits/periodicity_07.output
@@ -0,0 +1,5 @@
+
+DEAL::Testing version 0
+DEAL::Testing version 1
+DEAL::Testing version 2
+DEAL::Testing version 3
-- 
2.39.5