From 275843aea198cb2fd7cd7417a540659651aaba3c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 18 Sep 2019 08:02:23 -0600 Subject: [PATCH] Add a test. --- tests/feinterface/face_orientation_01.cc | 106 +++++++++++++++++++ tests/feinterface/face_orientation_01.output | 5 + 2 files changed, 111 insertions(+) create mode 100644 tests/feinterface/face_orientation_01.cc create mode 100644 tests/feinterface/face_orientation_01.output diff --git a/tests/feinterface/face_orientation_01.cc b/tests/feinterface/face_orientation_01.cc new file mode 100644 index 0000000000..019d9d455d --- /dev/null +++ b/tests/feinterface/face_orientation_01.cc @@ -0,0 +1,106 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 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. +// +// --------------------------------------------------------------------- + + +// Check how things go when there are faces with odd face +// orientations. If the quadrature points were ordered differently on +// the two sides of an interface, i.e., if each FEFaceValues used the +// coordinate system of the cell it's on rather than the common +// coordinate system of the face, then one would get different values +// for the same shape function on the q'th quadrature point from both +// sides. We can check that by using a continuous element and +// computing the jump -- it ought to be zero at every point. + +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + + +template +void +test() +{ + for (unsigned int twist = 0; twist < 4; ++twist) + { + Triangulation triangulation; + GridGenerator::moebius(triangulation, 7, twist, 1.0, 0.2); + + DoFHandler dofh(triangulation); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + + for (auto cell : dofh.active_cell_iterators()) + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + if (cell->at_boundary(f) == false) + { + fiv.reinit(cell, + f, + numbers::invalid_unsigned_int, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + numbers::invalid_unsigned_int); + + Assert(fiv.get_fe_face_values(0).get_cell() == cell, + ExcInternalError()); + Assert(fiv.get_fe_face_values(1).get_cell() == cell->neighbor(f), + ExcInternalError()); + + // We have a continuous element, so the set of DoFs that + // are active on the interface is twice the number of + // DoFs per cell (namely, the two adjacent cells) minus + // the common set of DoFs on the face itself. + AssertDimension(fiv.n_current_interface_dofs(), + 2 * fe.dofs_per_cell - fe.dofs_per_face); + Assert(!fiv.at_boundary(), ExcInternalError()); + + for (unsigned int q = 0; q < fiv.get_quadrature().size(); ++q) + for (unsigned int i = 0; i < fiv.n_current_interface_dofs(); + ++i) + Assert(std::abs(fiv.jump(i, q)) <= 1e-12, ExcInternalError()); + } + + deallog << "Twist " << twist << "*pi/2 seems to work just fine..." + << std::endl; + } +} + + + +int +main() +{ + initlog(); + test<3>(); +} diff --git a/tests/feinterface/face_orientation_01.output b/tests/feinterface/face_orientation_01.output new file mode 100644 index 0000000000..ee8064635d --- /dev/null +++ b/tests/feinterface/face_orientation_01.output @@ -0,0 +1,5 @@ + +DEAL::Twist 0*pi/2 seems to work just fine... +DEAL::Twist 1*pi/2 seems to work just fine... +DEAL::Twist 2*pi/2 seems to work just fine... +DEAL::Twist 3*pi/2 seems to work just fine... -- 2.39.5