Fix Matrixfree for periodic faces and hanging nodes

diff --git a/doc/news/changes/minor/20191103PeterMunch b/doc/news/changes/minor/20191103PeterMunch
new file mode 100644 (file)
index 0000000..fbf6848
--- /dev/null
+++ b/doc/news/changes/minor/20191103PeterMunch
@@ -0,0 +1,5 @@
+Fixed: FESubfaceValues::reinit() would previously throw an exception when 
+called on a face behind a refined periodic neighbor. This case is now 
+handled properly.
+<br>
+(Peter Munch, 2019/11/03)

diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index 8f72fc56bedd0b7e4474c0487e8768df962f2497..74f7fe168e5abca1f6d5e9972663698052d78225 100644 (file)
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -222,7 +222,12 @@ namespace internal
                   // faces properly
                   const CellId id_mine = dcell->id();
                   if (use_active_cells && neighbor->has_children())
-                    for (unsigned int c = 0; c < dcell->face(f)->n_children();
+                    for (unsigned int c = 0;
+                         c < (dcell->has_periodic_neighbor(f) ?
+                                dcell->periodic_neighbor(f)
+                                  ->face(dcell->periodic_neighbor_face_no(f))
+                                  ->n_children() :
+                                dcell->face(f)->n_children());
                          ++c)
                       {
                         typename dealii::Triangulation<dim>::cell_iterator

diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index 914733969f809cd3d9a51916ee1efea55972a98b..7bebc737820b16a737e300a9a8ec7aa6fab63517 100644 (file)
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -5029,8 +5029,23 @@ FESubfaceValues<dim, spacedim>::reinit(
 {
   Assert(face_no < GeometryInfo<dim>::faces_per_cell,
          ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
-  Assert(subface_no < cell->face(face_no)->n_children(),
-         ExcIndexRange(subface_no, 0, cell->face(face_no)->n_children()));
+  // We would like to check for subface_no < cell->face(face_no)->n_children(),
+  // but unfortunately the current function is also called for
+  // faces without children for periodic faces, which have hanging nodes on
+  // the other side (see
+  // include/deal.II/matrix_free/mapping_info.templates.h).
+  Assert(subface_no < cell->face(face_no)->n_children() ||
+           (cell->has_periodic_neighbor(face_no) &&
+            subface_no < cell->periodic_neighbor(face_no)
+                           ->face(cell->periodic_neighbor_face_no(face_no))
+                           ->n_children()),
+         ExcIndexRange(subface_no,
+                       0,
+                       (cell->has_periodic_neighbor(face_no) ?
+                          cell->periodic_neighbor(face_no)
+                            ->face(cell->periodic_neighbor_face_no(face_no))
+                            ->n_children() :
+                          cell->face(face_no)->n_children())));
 
   this->maybe_invalidate_previous_present_cell(cell);
   reset_pointer_in_place_if_possible<

diff --git a/tests/fe/fe_subface_values_01.cc b/tests/fe/fe_subface_values_01.cc
new file mode 100644 (file)
index 0000000..b881f54
--- /dev/null
+++ b/tests/fe/fe_subface_values_01.cc
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test FESubfaceValues<dim, spacedim>::reinit(
+// const typename Triangulation<dim, spacedim>::cell_iterator &,
+// const unsigned int, const unsigned int) for periodic faces with hanging
+// nodes.
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include "../tests.h"
+
+#define PRECISION 4
+
+template <int dim>
+void
+run(unsigned int degree, unsigned int n_q_points)
+{
+  Triangulation<dim> tria;
+
+  // Create grid
+  GridGenerator::hyper_cube(tria);
+
+  tria.begin_active()->face(0)->set_boundary_id(10);
+  tria.begin_active()->face(1)->set_boundary_id(11);
+  tria.begin_active()->face(2)->set_boundary_id(12);
+  tria.begin_active()->face(3)->set_boundary_id(13);
+
+  // Add periodic boundary conds
+  std::vector<
+    GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
+    periodic_faces;
+  GridTools::collect_periodic_faces(tria, 10, 11, 0, periodic_faces);
+  GridTools::collect_periodic_faces(tria, 12, 13, 1, periodic_faces);
+  tria.add_periodicity(periodic_faces);
+
+  // refine grid
+  tria.refine_global(1);
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  // Output grid file for each processor and global
+  GridOut grid_out;
+  grid_out.write_mesh_per_processor_as_vtu(tria, "mesh");
+
+  FE_Q<dim>              fe(degree);                       // dummy
+  QGaussLobatto<dim - 1> quad(n_q_points);                 // dummy
+  UpdateFlags            flags = update_quadrature_points; // dummy
+
+  FESubfaceValues<dim> sub_face(fe, quad, flags);
+  sub_face.reinit(CellId(0, {1}).to_cell(tria), 0, 0);
+  for (auto p : sub_face.get_quadrature_points())
+    deallog << p << std::endl;
+  deallog << std::endl;
+
+  sub_face.reinit(CellId(0, {1}).to_cell(tria), 1, 0);
+  for (auto p : sub_face.get_quadrature_points())
+    deallog << p << std::endl;
+  deallog << std::endl;
+}
+
+int
+main()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(PRECISION);
+  deallog.attach(logfile);
+
+  run<2>(2 /*=degree*/, 3 /*qpoints*/);
+}

diff --git a/tests/fe/fe_subface_values_01.output b/tests/fe/fe_subface_values_01.output
new file mode 100644 (file)
index 0000000..e27f15d
--- /dev/null
+++ b/tests/fe/fe_subface_values_01.output
@@ -0,0 +1,9 @@
+
+DEAL::0.5000 0.000
+DEAL::0.5000 0.1250
+DEAL::0.5000 0.2500
+DEAL::
+DEAL::1.000 0.000
+DEAL::1.000 0.1250
+DEAL::1.000 0.2500
+DEAL::