From 275843aea198cb2fd7cd7417a540659651aaba3c Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+  for (unsigned int twist = 0; twist < 4; ++twist)
+    {
+      Triangulation<dim> triangulation;
+      GridGenerator::moebius(triangulation, 7, twist, 1.0, 0.2);
+
+      DoFHandler<dim> dofh(triangulation);
+      FE_Q<dim>       fe(1);
+      dofh.distribute_dofs(fe);
+
+      MappingQ<dim> mapping(1);
+      UpdateFlags   update_flags = update_values;
+
+      FEInterfaceValues<dim> fiv(mapping,
+                                 fe,
+                                 QGauss<dim - 1>(fe.degree + 1),
+                                 update_flags);
+
+
+      for (auto cell : dofh.active_cell_iterators())
+        for (unsigned int f = 0; f < GeometryInfo<dim>::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