From 035699818b98bd1adc13e52f358cb725715b0b2a Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 31 Oct 2019 22:55:33 +0100 Subject: [PATCH] Fix Matrixfree for periodic faces and hanging nodes --- doc/news/changes/minor/20191103PeterMunch | 5 + .../deal.II/matrix_free/face_setup_internal.h | 7 +- source/fe/fe_values.cc | 19 +++- tests/fe/fe_subface_values_01.cc | 96 +++++++++++++++++++ tests/fe/fe_subface_values_01.output | 9 ++ 5 files changed, 133 insertions(+), 3 deletions(-) create mode 100644 doc/news/changes/minor/20191103PeterMunch create mode 100644 tests/fe/fe_subface_values_01.cc create mode 100644 tests/fe/fe_subface_values_01.output diff --git a/doc/news/changes/minor/20191103PeterMunch b/doc/news/changes/minor/20191103PeterMunch new file mode 100644 index 0000000000..fbf68480c5 --- /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. +
+(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 8f72fc56be..74f7fe168e 100644 --- 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::cell_iterator diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 914733969f..7bebc73782 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -5029,8 +5029,23 @@ FESubfaceValues::reinit( { Assert(face_no < GeometryInfo::faces_per_cell, ExcIndexRange(face_no, 0, GeometryInfo::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 index 0000000000..b881f54f31 --- /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::reinit( +// const typename Triangulation::cell_iterator &, +// const unsigned int, const unsigned int) for periodic faces with hanging +// nodes. + +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +#define PRECISION 4 + +template +void +run(unsigned int degree, unsigned int n_q_points) +{ + Triangulation 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::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 fe(degree); // dummy + QGaussLobatto quad(n_q_points); // dummy + UpdateFlags flags = update_quadrature_points; // dummy + + FESubfaceValues 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 index 0000000000..e27f15dca6 --- /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:: -- 2.39.5