]> https://gitweb.dealii.org/ - dealii.git/commitdiff
introduce FEInterfaceValues 8660/head
authorTimo Heister <timo.heister@gmail.com>
Sun, 25 Aug 2019 02:01:14 +0000 (22:01 -0400)
committerTimo Heister <timo.heister@gmail.com>
Tue, 17 Sep 2019 00:39:34 +0000 (20:39 -0400)
19 files changed:
doc/news/changes/major/20190824Heister [new file with mode: 0644]
include/deal.II/fe/fe_interface_values.h [new file with mode: 0644]
tests/feinterface/CMakeLists.txt [new file with mode: 0644]
tests/feinterface/fe_interface_values_01.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_01.output [new file with mode: 0644]
tests/feinterface/fe_interface_values_02.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_02.output [new file with mode: 0644]
tests/feinterface/fe_interface_values_03.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_03.output [new file with mode: 0644]
tests/feinterface/fe_interface_values_04.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_04.output [new file with mode: 0644]
tests/feinterface/fe_interface_values_05.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_05.output [new file with mode: 0644]
tests/feinterface/fe_interface_values_06.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_06.output [new file with mode: 0644]
tests/feinterface/fe_interface_values_07.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_07.output [new file with mode: 0644]
tests/feinterface/step-12.cc [new file with mode: 0644]
tests/feinterface/step-12.output [new file with mode: 0644]

diff --git a/doc/news/changes/major/20190824Heister b/doc/news/changes/major/20190824Heister
new file mode 100644 (file)
index 0000000..dd3d0eb
--- /dev/null
@@ -0,0 +1,5 @@
+New: The FEInterfaceValues class provides a new abstraction to assemble
+interface terms between two neighboring cells. This is commonly used in
+Discontinous Galerkin methods.
+<br>
+(Timo Heister, 2019/08/24)
diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h
new file mode 100644 (file)
index 0000000..fef1461
--- /dev/null
@@ -0,0 +1,694 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_fe_interface_values_h
+#define dealii_fe_interface_values_h
+
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * FEInterfaceValues is a data structure to access and assemble finite element
+ * data on interfaces between two cells of a mesh.
+ *
+ * It provides a way to access averages, jump terms, and similar operations used
+ * in Discontinuous Galerkin methods on a face between two neighboring cells.
+ * This allows the computation of typical mesh-dependent linear or bilinear
+ * forms in a similar way as FEValues does for cells and FEFaceValues does for
+ * faces. In
+ * the literature, the faces between neighboring cells are called "inner
+ * interfaces" or "facets".
+ *
+ * Internally, this class provides an abstraction for two FEFaceValues
+ * objects (or FESubfaceValues when using adaptive refinement). The class
+ * introduces a new "interface dof index" that walks over
+ * the union of the dof indices of the two FEFaceValues objects. Helper
+ * functions allow translating between the new "interface dof index" and the
+ * corresponding "cell index" (0 for the first cell, 1 for the second cell)
+ * and "dof index" within that cell.
+ *
+ * The class is made to be used inside MeshWorker::mesh_loop(). It is intended
+ * to be a low level replacement for MeshWorker and LocalIntegrators and a
+ * higher level abstraction compared to assembling face terms manually.
+ *
+ * @author Timo Heister, 2019.
+ */
+template <int dim, int spacedim = dim>
+class FEInterfaceValues
+{
+public:
+  /**
+   * Construct the FEInterfaceValues with a single FiniteElement (same on both
+   * sides of the facet). The FEFaceValues objects will be initialized with
+   * the given @p mapping, @p quadrature, and @p update_flags.
+   */
+  FEInterfaceValues(const Mapping<dim, spacedim> &      mapping,
+                    const FiniteElement<dim, spacedim> &fe,
+                    const Quadrature<dim - 1> &         quadrature,
+                    const UpdateFlags                   update_flags);
+
+  /**
+   * Construct the FEInterfaceValues with a single FiniteElement and
+   * a Q1 Mapping.
+   *
+   * See the constructor above.
+   */
+  FEInterfaceValues(const FiniteElement<dim, spacedim> &fe,
+                    const Quadrature<dim - 1> &         quadrature,
+                    const UpdateFlags                   update_flags);
+
+  /**
+   * Re-initialize this object to be used on a new interface given by two faces
+   * of two neighboring cells. The `cell` and `cell_neighbor` cells will be
+   * referred to through `cell_index` zero and one after this call in all places
+   * where one needs to identify the two cells adjacent to the interface.
+   *
+   * Use numbers::invalid_unsigned_int for @p sub_face_no or @p
+   * sub_face_no_neighbor to indicate that you want to work on the entire face,
+   * not a sub-face.
+   *
+   * The arguments (including their order) are identical to the @p face_worker arguments
+   * in MeshWorker::mesh_loop().
+   */
+  template <class CellIteratorType>
+  void
+  reinit(const CellIteratorType &cell,
+         const unsigned int &    face_no,
+         const unsigned int &    sub_face_no,
+         const CellIteratorType &cell_neighbor,
+         const unsigned int &    face_no_neighbor,
+         const unsigned int &    sub_face_no_neighbor);
+
+  /**
+   * Re-initialize this object to be used on an interface given by a single face
+   * @p face_no of the cell @p cell. This is useful to use FEInterfaceValues
+   * on boundaries of the domain.
+   *
+   * As a consequence, members like jump() will assume a value of zero for the
+   * values on the "other" side. Note that no sub_face_number is needed as a
+   * boundary face can not neighbor a finer cell.
+   *
+   * After calling this function at_boundary() will return true.
+   */
+  template <class CellIteratorType>
+  void
+  reinit(const CellIteratorType &cell, const unsigned int &face_no);
+
+  /**
+   * Return a reference to the FEFaceValues or FESubfaceValues object
+   * of the specified cell of the interface.
+   *
+   * The @p cell_index is either 0 or 1 and corresponds to the cell index
+   * returned by interface_dof_to_cell_and_dof_index().
+   */
+  const FEFaceValuesBase<dim, spacedim> &
+  get_fe_face_values(const unsigned int cell_index) const;
+
+  /**
+   * Return a reference to the quadrature object in use.
+   */
+  const Quadrature<dim - 1> &
+  get_quadrature() const;
+
+  /**
+   * Return the update flags set.
+   */
+  const UpdateFlags
+  get_update_flags() const;
+
+  /**
+   * @name Functions to query information on a given interface
+   * @{
+   */
+
+  /**
+   * Return if the current interface is a boundary face or an internal
+   * face with two adjacent cells.
+   *
+   * See the corresponding reinit() functions for details.
+   */
+  bool
+  at_boundary() const;
+
+  /**
+   * Return the vector of JxW values for each quadrature point.
+   *
+   * @dealiiRequiresUpdateFlags{update_JxW_values}
+   */
+  const std::vector<double> &
+  get_JxW_values() const;
+
+  /**
+   * Return the normal vector of the interface in each quadrature point.
+   *
+   * The return value is identical to get_fe_face_values(0).get_normal_vectors()
+   * and therefore, are outside normal vectors from the perspective of the
+   * first cell of this interface.
+   *
+   * @dealiiRequiresUpdateFlags{update_normal_vectors}
+   */
+  const std::vector<Tensor<1, spacedim>> &
+  get_normal_vectors() const;
+
+  /**
+   * Return a reference to the quadrature points in real space.
+   *
+   * @dealiiRequiresUpdateFlags{update_quadrature_points}
+   */
+  const std::vector<Point<spacedim>> &
+  get_quadrature_points() const;
+
+
+  /**
+   * Return the number of DoFs (or shape functions) on the current interface.
+   *
+   * @note This number is only available after a call to reinit() and can change
+   * from one call to reinit() to the next. For example, on a boundary interface
+   * it is equal to the number of dofs of the single FEFaceValues object, while
+   * it is twice that for an interior interface for a DG element. For a
+   * continuous element, it is slightly smaller because the two cells on the
+   * interface share some of the dofs.
+   */
+  unsigned
+  n_current_interface_dofs() const;
+
+  /**
+   * Return the set of joint DoF indices. This includes indices from both cells.
+   *
+   * @note This function is only available after a call to reinit() and can change
+   * from one call to reinit() to the next.
+   */
+  std::vector<types::global_dof_index>
+  get_interface_dof_indices() const;
+
+  /**
+   * Convert an interface dof index into the corresponding local DoF indices of
+   * the two cells. If an interface DoF is only active on one of the
+   * cells, the other index will be numbers::invalid_unsigned_int.
+   *
+   * For discontinuous finite elements each interface dof will correspond to
+   * exactly one DoF index.
+   *
+   * @note This function is only available after a call to reinit() and can change
+   * from one call to reinit() to the next.
+   */
+  std::array<unsigned int, 2>
+  interface_dof_to_dof_indices(const unsigned int interface_dof_index) const;
+
+  /**
+   * Return the normal in a given quadrature point.
+   *
+   * The normal points in outwards direction as seen from the first cell of
+   * this interface.
+   *
+   * @dealiiRequiresUpdateFlags{update_normal_vectors}
+   */
+  Tensor<1, spacedim>
+  normal(const unsigned int q_point_index) const;
+
+  /**
+   * @}
+   */
+
+  /**
+   * @name Functions to evaluate data of the shape functions
+   * @{
+   */
+
+  /**
+   * Return component @p component of the value of the shape function
+   * with interface dof index @p interface_dof_index in
+   * quadrature point @p q_point.
+   *
+   * The argument @p here_or_there selects between the value on cell 0 (here, @p true)
+   * and cell 1 (there, @p false). You can also interpret it as "upstream" (@p true)
+   * and "downstream" (@p false) as defined by the direction of the normal
+   * vector
+   * in this quadrature point. If @p here_or_there is true, the shape
+   * functions from the first cell of the interface is used.
+   *
+   * In other words, this function returns the limit of the value of the shape
+   * function in the given quadrature point when approaching it from one of the
+   * two cells of the interface.
+   *
+   * @note This function is typically used to pick the upstream or downstream
+   * value based on a direction. This can be achieved by using
+   * <code>(direction * normal)>0</code> as the first argument of this
+   * function.
+   */
+  double
+  shape_value(const bool         here_or_there,
+              const unsigned int interface_dof_index,
+              const unsigned int q_point,
+              const unsigned int component = 0) const;
+
+  /**
+   * Return the jump $[u]=u_{\text{cell0}} - u_{\text{cell1}}$ on the
+   * interface
+   * for the shape function @p interface_dof_index at the quadrature point
+   * @p q_point of component @p component.
+   *
+   * Note that one can define the jump in
+   * different ways (the value "there" minus the value "here", or the other way
+   * way around; both are used in the finite element literature). The definition
+   * here uses "value here minus value there", as seen from the first cell.
+   *
+   * If this is a boundary face (at_boundary() returns true), then
+   * $[u]=u_{\text{cell0}}$.
+   */
+  double
+  jump(const unsigned int interface_dof_index,
+       const unsigned int q_point,
+       const unsigned int component = 0) const;
+
+  /**
+   * Return the average $\{u\}=\frac{1}{2}u_{\text{cell0}} +
+   * \frac{1}{2}u_{\text{cell1}}$ on the interface
+   * for the shape function @p interface_dof_index at the quadrature point
+   * @p q_point of component @p component.
+   *
+   * If this is a boundary face (at_boundary() returns true), then
+   * $\{u\}=u_{\text{cell0}}$.
+   */
+  double
+  average(const unsigned int interface_dof_index,
+          const unsigned int q_point,
+          const unsigned int component = 0) const;
+
+
+  /**
+   * @}
+   */
+
+private:
+  /**
+   * The list of DoF indices for the current interface, filled in reinit().
+   */
+  std::vector<types::global_dof_index> interface_dof_indices;
+
+  /**
+   * The mapping from interface dof to the two local dof indices of the
+   * FeFaceValues objects. If an interface DoF is only active on one of the
+   * cells, the other one will have numbers::invalid_unsigned_int.
+   */
+  std::vector<std::array<unsigned int, 2>> dofmap;
+
+  /**
+   * The FEFaceValues object for the current cell.
+   */
+  FEFaceValues<dim> internal_fe_face_values;
+
+  /**
+   * The FEFaceValues object for the current cell if the cell is refined.
+   */
+  FESubfaceValues<dim> internal_fe_subface_values;
+
+  /**
+   * The FEFaceValues object for the neighboring cell.
+   */
+  FEFaceValues<dim> internal_fe_face_values_neighbor;
+
+  /**
+   * The FEFaceValues object for the neighboring cell if the cell is refined.
+   */
+  FESubfaceValues<dim> internal_fe_subface_values_neighbor;
+
+  /**
+   * Pointer to internal_fe_face_values or internal_fe_subface_values,
+   * respectively as determined in reinit().
+   */
+  FEFaceValuesBase<dim> *fe_face_values;
+
+  /**
+   * Pointer to internal_fe_face_values_neighbor,
+   * internal_fe_subface_values_neighbor, or nullptr, respectively
+   * as determined in reinit().
+   */
+  FEFaceValuesBase<dim> *fe_face_values_neighbor;
+};
+
+
+
+#ifndef DOXYGEN
+
+/*---------------------- Inline functions ---------------------*/
+
+template <int dim, int spacedim>
+FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
+  const Mapping<dim, spacedim> &      mapping,
+  const FiniteElement<dim, spacedim> &fe,
+  const Quadrature<dim - 1> &         quadrature,
+  const UpdateFlags                   update_flags)
+  : internal_fe_face_values(mapping, fe, quadrature, update_flags)
+  , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
+  , internal_fe_face_values_neighbor(mapping, fe, quadrature, update_flags)
+  , internal_fe_subface_values_neighbor(mapping, fe, quadrature, update_flags)
+  , fe_face_values(nullptr)
+  , fe_face_values_neighbor(nullptr)
+{}
+
+
+
+template <int dim, int spacedim>
+FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
+  const FiniteElement<dim, spacedim> &fe,
+  const Quadrature<dim - 1> &         quadrature,
+  const UpdateFlags                   update_flags)
+  : internal_fe_face_values(StaticMappingQ1<dim, spacedim>::mapping,
+                            fe,
+                            quadrature,
+                            update_flags)
+  , internal_fe_subface_values(StaticMappingQ1<dim, spacedim>::mapping,
+                               fe,
+                               quadrature,
+                               update_flags)
+  , internal_fe_face_values_neighbor(StaticMappingQ1<dim, spacedim>::mapping,
+                                     fe,
+                                     quadrature,
+                                     update_flags)
+  , internal_fe_subface_values_neighbor(StaticMappingQ1<dim, spacedim>::mapping,
+                                        fe,
+                                        quadrature,
+                                        update_flags)
+  , fe_face_values(nullptr)
+  , fe_face_values_neighbor(nullptr)
+{}
+
+
+
+template <int dim, int spacedim>
+template <class CellIteratorType>
+void
+FEInterfaceValues<dim, spacedim>::reinit(
+  const CellIteratorType &cell,
+  const unsigned int &    face_no,
+  const unsigned int &    sub_face_no,
+  const CellIteratorType &cell_neighbor,
+  const unsigned int &    face_no_neighbor,
+  const unsigned int &    sub_face_no_neighbor)
+{
+  if (sub_face_no == numbers::invalid_unsigned_int)
+    {
+      internal_fe_face_values.reinit(cell, face_no);
+      fe_face_values = &internal_fe_face_values;
+    }
+  else
+    {
+      internal_fe_subface_values.reinit(cell, face_no, sub_face_no);
+      fe_face_values = &internal_fe_subface_values;
+    }
+  if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
+    {
+      internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor);
+      fe_face_values_neighbor = &internal_fe_face_values_neighbor;
+    }
+  else
+    {
+      internal_fe_subface_values_neighbor.reinit(cell_neighbor,
+                                                 face_no_neighbor,
+                                                 sub_face_no_neighbor);
+      fe_face_values_neighbor = &internal_fe_subface_values_neighbor;
+    }
+
+  // Set up dof mapping and remove duplicates (for continuous elements).
+  {
+    // Get dof indices first:
+    std::vector<types::global_dof_index> v(
+      fe_face_values->get_fe().n_dofs_per_cell());
+    cell->get_dof_indices(v);
+    std::vector<types::global_dof_index> v2(
+      fe_face_values_neighbor->get_fe().n_dofs_per_cell());
+    cell_neighbor->get_dof_indices(v2);
+
+    // Fill a map from the global dof index to the left and right
+    // local index.
+    std::map<types::global_dof_index, std::pair<unsigned int, unsigned int>>
+                                          tempmap;
+    std::pair<unsigned int, unsigned int> invalid_entry(
+      numbers::invalid_unsigned_int, numbers::invalid_unsigned_int);
+
+    for (unsigned int i = 0; i < v.size(); ++i)
+      {
+        // If not already existing, add an invalid entry:
+        auto result = tempmap.insert(std::make_pair(v[i], invalid_entry));
+        result.first->second.first = i;
+      }
+
+    for (unsigned int i = 0; i < v2.size(); ++i)
+      {
+        // If not already existing, add an invalid entry:
+        auto result = tempmap.insert(std::make_pair(v2[i], invalid_entry));
+        result.first->second.second = i;
+      }
+
+    // Transfer from the map to the sorted std::vectors.
+    dofmap.resize(tempmap.size());
+    interface_dof_indices.resize(tempmap.size());
+    unsigned int idx = 0;
+    for (auto &x : tempmap)
+      {
+        interface_dof_indices[idx] = x.first;
+        dofmap[idx]                = {x.second.first, x.second.second};
+        ++idx;
+      }
+  }
+}
+
+
+
+template <int dim, int spacedim>
+template <class CellIteratorType>
+void
+FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
+                                         const unsigned int &    face_no)
+{
+  internal_fe_face_values.reinit(cell, face_no);
+  fe_face_values          = &internal_fe_face_values;
+  fe_face_values_neighbor = nullptr;
+
+  interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
+  cell->get_dof_indices(interface_dof_indices);
+
+
+  dofmap.resize(interface_dof_indices.size());
+
+  for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
+    {
+      dofmap[i] = {i, numbers::invalid_unsigned_int};
+    }
+}
+
+
+
+template <int dim, int spacedim>
+const std::vector<double> &
+FEInterfaceValues<dim, spacedim>::get_JxW_values() const
+{
+  Assert(fe_face_values != nullptr,
+         ExcMessage("This call requires a call to reinit() first."));
+  return fe_face_values->get_JxW_values();
+}
+
+
+
+template <int dim, int spacedim>
+const std::vector<Tensor<1, spacedim>> &
+FEInterfaceValues<dim, spacedim>::get_normal_vectors() const
+{
+  Assert(fe_face_values != nullptr,
+         ExcMessage("This call requires a call to reinit() first."));
+  return fe_face_values->get_normal_vectors();
+}
+
+
+
+template <int dim, int spacedim>
+const Quadrature<dim - 1> &
+FEInterfaceValues<dim, spacedim>::get_quadrature() const
+{
+  return internal_fe_face_values.get_quadrature();
+}
+
+
+
+template <int dim, int spacedim>
+const std::vector<Point<spacedim>> &
+FEInterfaceValues<dim, spacedim>::get_quadrature_points() const
+{
+  Assert(fe_face_values != nullptr,
+         ExcMessage("This call requires a call to reinit() first."));
+  return fe_face_values->get_quadrature_points();
+}
+
+
+
+template <int dim, int spacedim>
+const UpdateFlags
+FEInterfaceValues<dim, spacedim>::get_update_flags() const
+{
+  return internal_fe_face_values.get_update_flags();
+}
+
+
+
+template <int dim, int spacedim>
+unsigned
+FEInterfaceValues<dim, spacedim>::n_current_interface_dofs() const
+{
+  Assert(
+    interface_dof_indices.size() > 0,
+    ExcMessage(
+      "n_current_interface_dofs() is only available after a call to reinit()."));
+  return interface_dof_indices.size();
+}
+
+
+
+template <int dim, int spacedim>
+bool
+FEInterfaceValues<dim, spacedim>::at_boundary() const
+{
+  return fe_face_values_neighbor == nullptr;
+}
+
+
+
+template <int dim, int spacedim>
+std::vector<types::global_dof_index>
+FEInterfaceValues<dim, spacedim>::get_interface_dof_indices() const
+{
+  return interface_dof_indices;
+}
+
+
+
+template <int dim, int spacedim>
+std::array<unsigned int, 2>
+FEInterfaceValues<dim, spacedim>::interface_dof_to_dof_indices(
+  const unsigned int interface_dof_index) const
+{
+  AssertIndexRange(interface_dof_index, n_current_interface_dofs());
+  return dofmap[interface_dof_index];
+}
+
+
+
+template <int dim, int spacedim>
+const FEFaceValuesBase<dim, spacedim> &
+FEInterfaceValues<dim, spacedim>::get_fe_face_values(
+  const unsigned int cell_index) const
+{
+  AssertIndexRange(cell_index, 2);
+  Assert(
+    cell_index == 0 || !at_boundary(),
+    ExcMessage(
+      "You are on a boundary, so you can only ask for the first FEFaceValues object."));
+
+  return (cell_index == 0) ? *fe_face_values : *fe_face_values_neighbor;
+}
+
+
+
+template <int dim, int spacedim>
+Tensor<1, spacedim>
+FEInterfaceValues<dim, spacedim>::normal(const unsigned int q_point_index) const
+{
+  return fe_face_values->normal_vector(q_point_index);
+}
+
+
+
+template <int dim, int spacedim>
+double
+FEInterfaceValues<dim, spacedim>::shape_value(
+  const bool         here_or_there,
+  const unsigned int interface_dof_index,
+  const unsigned int q_point,
+  const unsigned int component) const
+{
+  const auto dofpair = dofmap[interface_dof_index];
+
+  if (here_or_there && dofpair[0] != numbers::invalid_unsigned_int)
+    return get_fe_face_values(0).shape_value_component(dofpair[0],
+                                                       q_point,
+                                                       component);
+  if (!here_or_there && dofpair[1] != numbers::invalid_unsigned_int)
+    return get_fe_face_values(1).shape_value_component(dofpair[1],
+                                                       q_point,
+                                                       component);
+
+  return 0.0;
+}
+
+
+
+template <int dim, int spacedim>
+double
+FEInterfaceValues<dim, spacedim>::jump(const unsigned int interface_dof_index,
+                                       const unsigned int q_point,
+                                       const unsigned int component) const
+{
+  const auto dofpair = dofmap[interface_dof_index];
+
+  double value = 0.0;
+
+  if (dofpair[0] != numbers::invalid_unsigned_int)
+    value += get_fe_face_values(0).shape_value_component(dofpair[0],
+                                                         q_point,
+                                                         component);
+  if (dofpair[1] != numbers::invalid_unsigned_int)
+    value -= get_fe_face_values(1).shape_value_component(dofpair[1],
+                                                         q_point,
+                                                         component);
+  return value;
+}
+
+
+
+template <int dim, int spacedim>
+double
+FEInterfaceValues<dim, spacedim>::average(
+  const unsigned int interface_dof_index,
+  const unsigned int q_point,
+  const unsigned int component) const
+{
+  const auto dofpair = dofmap[interface_dof_index];
+
+  if (at_boundary())
+    return 1.0 * get_fe_face_values(0).shape_value_component(dofpair[0],
+                                                             q_point,
+                                                             component);
+
+  double value = 0.0;
+
+  if (dofpair[0] != numbers::invalid_unsigned_int)
+    value += 0.5 * get_fe_face_values(0).shape_value_component(dofpair[0],
+                                                               q_point,
+                                                               component);
+  if (dofpair[1] != numbers::invalid_unsigned_int)
+    value += 0.5 * get_fe_face_values(1).shape_value_component(dofpair[1],
+                                                               q_point,
+                                                               component);
+
+  return value;
+}
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/tests/feinterface/CMakeLists.txt b/tests/feinterface/CMakeLists.txt
new file mode 100644 (file)
index 0000000..d33eafa
--- /dev/null
@@ -0,0 +1,4 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+DEAL_II_PICKUP_TESTS()
diff --git a/tests/feinterface/fe_interface_values_01.cc b/tests/feinterface/fe_interface_values_01.cc
new file mode 100644 (file)
index 0000000..4565e20
--- /dev/null
@@ -0,0 +1,174 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test basic properties of FEInterfaceValues, global refinement, no-hp
+
+#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/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+inspect_fiv(FEInterfaceValues<dim> &fiv)
+{
+  deallog << "at_boundary(): " << fiv.at_boundary() << "\n"
+          << "n_current_interface_dofs(): " << fiv.n_current_interface_dofs()
+          << "\n";
+
+  std::vector<types::global_dof_index> indices =
+    fiv.get_interface_dof_indices();
+  Assert(indices.size() == fiv.n_current_interface_dofs(), ExcInternalError());
+
+  deallog << "interface_dof_indices: ";
+  for (auto i : indices)
+    deallog << i << " ";
+  deallog << "\n";
+
+
+  unsigned int idx = 0;
+  for (auto v : indices)
+    {
+      deallog << "  index " << idx << " global_dof_index:" << v << ":\n";
+
+      const auto pair = fiv.interface_dof_to_dof_indices(idx);
+      deallog << "    dof indices: " << pair[0] << " | " << pair[1] << "\n";
+
+      ++idx;
+    }
+
+  deallog << "\n";
+  deallog << std::endl;
+}
+
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test()
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+
+  DoFHandler<dim> dofh(tria);
+  FE_DGQ<dim>     fe(1);
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+
+  auto cell = dofh.begin();
+
+  deallog << "** interface between cell 0 and 1 **\n";
+
+  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    if (!cell->at_boundary(f))
+      {
+        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());
+        Assert(fiv.n_current_interface_dofs() == 2 * fe.n_dofs_per_cell(),
+               ExcInternalError());
+        Assert(!fiv.at_boundary(), ExcInternalError());
+
+        auto mycell = cell;
+        for (unsigned int c = 0; c < 2; ++c)
+          {
+            std::vector<types::global_dof_index> indices(fe.n_dofs_per_cell());
+            mycell->get_dof_indices(indices);
+            deallog << "cell " << c << ": ";
+            for (auto i : indices)
+              deallog << i << " ";
+            deallog << "\n";
+            ++mycell;
+          }
+
+        inspect_fiv(fiv);
+      }
+
+  deallog << "** boundary interface on cell 1 **\n";
+
+  {
+    ++cell;
+    fiv.reinit(cell, 1);
+    Assert(fiv.get_fe_face_values(0).get_cell() == cell, ExcInternalError());
+    Assert(fiv.n_current_interface_dofs() == fe.n_dofs_per_cell(),
+           ExcInternalError());
+    Assert(fiv.at_boundary(), ExcInternalError());
+    inspect_fiv(fiv);
+  }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>();
+}
diff --git a/tests/feinterface/fe_interface_values_01.output b/tests/feinterface/fe_interface_values_01.output
new file mode 100644 (file)
index 0000000..263e929
--- /dev/null
@@ -0,0 +1,39 @@
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3 
+cell 1: 4 5 6 7 
+at_boundary(): 0
+n_current_interface_dofs(): 8
+interface_dof_indices: 0 1 2 3 4 5 6 7 
+  index 0 global_dof_index:0:
+    dof indices: 0 | 4294967295
+  index 1 global_dof_index:1:
+    dof indices: 1 | 4294967295
+  index 2 global_dof_index:2:
+    dof indices: 2 | 4294967295
+  index 3 global_dof_index:3:
+    dof indices: 3 | 4294967295
+  index 4 global_dof_index:4:
+    dof indices: 4294967295 | 0
+  index 5 global_dof_index:5:
+    dof indices: 4294967295 | 1
+  index 6 global_dof_index:6:
+    dof indices: 4294967295 | 2
+  index 7 global_dof_index:7:
+    dof indices: 4294967295 | 3
+
+
+DEAL::** boundary interface on cell 1 **
+at_boundary(): 1
+n_current_interface_dofs(): 4
+interface_dof_indices: 4 5 6 7 
+  index 0 global_dof_index:4:
+    dof indices: 0 | 4294967295
+  index 1 global_dof_index:5:
+    dof indices: 1 | 4294967295
+  index 2 global_dof_index:6:
+    dof indices: 2 | 4294967295
+  index 3 global_dof_index:7:
+    dof indices: 3 | 4294967295
+
+
diff --git a/tests/feinterface/fe_interface_values_02.cc b/tests/feinterface/fe_interface_values_02.cc
new file mode 100644 (file)
index 0000000..b6c9a79
--- /dev/null
@@ -0,0 +1,140 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate jump(), average(), shape_value() of FEInterfaceValues
+
+#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/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(unsigned int fe_degree)
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+
+  DoFHandler<dim> dofh(tria);
+  FE_DGQ<dim>     fe(fe_degree);
+  deallog << fe.get_name() << std::endl;
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+
+  auto cell = dofh.begin();
+
+  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    if (!cell->at_boundary(f))
+      {
+        fiv.reinit(cell,
+                   f,
+                   numbers::invalid_unsigned_int,
+                   cell->neighbor(f),
+                   cell->neighbor_of_neighbor(f),
+                   numbers::invalid_unsigned_int);
+
+        const unsigned int n_dofs = fiv.n_current_interface_dofs();
+        Vector<double>     cell_vector(n_dofs);
+
+        const auto &q_points = fiv.get_quadrature_points();
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(true, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(true): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(false, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(false): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.jump(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "jump(): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.average(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "average(): " << cell_vector;
+      }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>(0);
+  test<2>(1);
+  test<3>(0);
+  test<3>(1);
+}
diff --git a/tests/feinterface/fe_interface_values_02.output b/tests/feinterface/fe_interface_values_02.output
new file mode 100644 (file)
index 0000000..0a6b89e
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::FE_DGQ<2>(0)
+DEAL::shape_value(true): 1.00000 0.00000 
+DEAL::shape_value(false): 0.00000 1.00000 
+DEAL::jump(): 1.00000 -1.00000 
+DEAL::average(): 0.500000 0.500000 
+DEAL::FE_DGQ<2>(1)
+DEAL::shape_value(true): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 0.00000 0.00000 
+DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.500000 0.00000 
+DEAL::jump(): 0.00000 0.500000 0.00000 0.500000 -0.500000 0.00000 -0.500000 0.00000 
+DEAL::average(): 0.00000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.00000 
+DEAL::FE_DGQ<3>(0)
+DEAL::shape_value(true): 1.00000 0.00000 
+DEAL::shape_value(false): 0.00000 1.00000 
+DEAL::jump(): 1.00000 -1.00000 
+DEAL::average(): 0.500000 0.500000 
+DEAL::FE_DGQ<3>(1)
+DEAL::shape_value(true): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 
+DEAL::jump(): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 
+DEAL::average(): 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 
diff --git a/tests/feinterface/fe_interface_values_03.cc b/tests/feinterface/fe_interface_values_03.cc
new file mode 100644 (file)
index 0000000..a649ada
--- /dev/null
@@ -0,0 +1,159 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate jump(), average(), shape_value() of FEInterfaceValues on an adaptive
+// mesh
+
+#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/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(unsigned int fe_degree)
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+  tria.begin()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  DoFHandler<dim> dofh(tria);
+  FE_DGQ<dim>     fe(fe_degree);
+  deallog << fe.get_name() << std::endl;
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+  auto cell = dofh.begin(1);
+  ++cell;
+
+  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    if (!cell->at_boundary(f))
+      {
+        if (!cell->neighbor_is_coarser(f))
+          continue;
+
+        auto nn = cell->neighbor_of_coarser_neighbor(f);
+        fiv.reinit(cell,
+                   f,
+                   numbers::invalid_unsigned_int,
+                   cell->neighbor(f),
+                   nn.first,
+                   nn.second);
+
+        const unsigned int n_dofs = fiv.n_current_interface_dofs();
+        Vector<double>     cell_vector(n_dofs);
+
+        const auto &q_points = fiv.get_quadrature_points();
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          deallog << "qpoint " << qpoint << ": " << q_points[qpoint]
+                  << std::endl;
+
+        for (unsigned int idx = 0; idx < n_dofs; ++idx)
+          {
+            const auto pair = fiv.interface_dof_to_dof_indices(idx);
+            deallog << "  idx: " << idx
+                    << " global: " << fiv.get_interface_dof_indices()[idx]
+                    << " dof indices: " << pair[0] << " | " << pair[1]
+                    << std::endl;
+          }
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(true, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(true): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(false, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(false): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.jump(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "jump(): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.average(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "average(): " << cell_vector;
+      }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>(0);
+  test<2>(1);
+  test<3>(0);
+  test<3>(1);
+}
diff --git a/tests/feinterface/fe_interface_values_03.output b/tests/feinterface/fe_interface_values_03.output
new file mode 100644 (file)
index 0000000..8be1c3f
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL::FE_DGQ<2>(0)
+DEAL::qpoint 0: 1.00000 0.250000
+DEAL::  idx: 0 global: 0 dof indices: 4294967295 | 0
+DEAL::  idx: 1 global: 2 dof indices: 0 | 4294967295
+DEAL::shape_value(true): 0.00000 0.500000 
+DEAL::shape_value(false): 0.500000 0.00000 
+DEAL::jump(): -0.500000 0.500000 
+DEAL::average(): 0.250000 0.250000 
+DEAL::FE_DGQ<2>(1)
+DEAL::qpoint 0: 1.00000 0.105662
+DEAL::qpoint 1: 1.00000 0.394338
+DEAL::  idx: 0 global: 0 dof indices: 4294967295 | 0
+DEAL::  idx: 1 global: 1 dof indices: 4294967295 | 1
+DEAL::  idx: 2 global: 2 dof indices: 4294967295 | 2
+DEAL::  idx: 3 global: 3 dof indices: 4294967295 | 3
+DEAL::  idx: 4 global: 8 dof indices: 0 | 4294967295
+DEAL::  idx: 5 global: 9 dof indices: 1 | 4294967295
+DEAL::  idx: 6 global: 10 dof indices: 2 | 4294967295
+DEAL::  idx: 7 global: 11 dof indices: 3 | 4294967295
+DEAL::shape_value(true): 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000 
+DEAL::shape_value(false): 0.375000 0.00000 0.125000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::jump(): -0.375000 0.00000 -0.125000 0.00000 0.00000 0.250000 0.00000 0.250000 
+DEAL::average(): 0.187500 0.00000 0.0625000 0.00000 0.00000 0.125000 0.00000 0.125000 
+DEAL::FE_DGQ<3>(0)
+DEAL::qpoint 0: 1.00000 0.250000 0.250000
+DEAL::  idx: 0 global: 0 dof indices: 4294967295 | 0
+DEAL::  idx: 1 global: 2 dof indices: 0 | 4294967295
+DEAL::shape_value(true): 0.00000 0.250000 
+DEAL::shape_value(false): 0.250000 0.00000 
+DEAL::jump(): -0.250000 0.250000 
+DEAL::average(): 0.125000 0.125000 
+DEAL::FE_DGQ<3>(1)
+DEAL::qpoint 0: 1.00000 0.105662 0.105662
+DEAL::qpoint 1: 1.00000 0.394338 0.105662
+DEAL::qpoint 2: 1.00000 0.105662 0.394338
+DEAL::qpoint 3: 1.00000 0.394338 0.394338
+DEAL::  idx: 0 global: 0 dof indices: 4294967295 | 0
+DEAL::  idx: 1 global: 1 dof indices: 4294967295 | 1
+DEAL::  idx: 2 global: 2 dof indices: 4294967295 | 2
+DEAL::  idx: 3 global: 3 dof indices: 4294967295 | 3
+DEAL::  idx: 4 global: 4 dof indices: 4294967295 | 4
+DEAL::  idx: 5 global: 5 dof indices: 4294967295 | 5
+DEAL::  idx: 6 global: 6 dof indices: 4294967295 | 6
+DEAL::  idx: 7 global: 7 dof indices: 4294967295 | 7
+DEAL::  idx: 8 global: 16 dof indices: 0 | 4294967295
+DEAL::  idx: 9 global: 17 dof indices: 1 | 4294967295
+DEAL::  idx: 10 global: 18 dof indices: 2 | 4294967295
+DEAL::  idx: 11 global: 19 dof indices: 3 | 4294967295
+DEAL::  idx: 12 global: 20 dof indices: 4 | 4294967295
+DEAL::  idx: 13 global: 21 dof indices: 5 | 4294967295
+DEAL::  idx: 14 global: 22 dof indices: 6 | 4294967295
+DEAL::  idx: 15 global: 23 dof indices: 7 | 4294967295
+DEAL::shape_value(true): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 
+DEAL::shape_value(false): 0.140625 0.00000 0.0468750 0.00000 0.0468750 0.00000 0.0156250 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::jump(): -0.140625 0.00000 -0.0468750 0.00000 -0.0468750 0.00000 -0.0156250 0.00000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 
+DEAL::average(): 0.0703125 0.00000 0.0234375 0.00000 0.0234375 0.00000 0.00781250 0.00000 0.00000 0.0312500 0.00000 0.0312500 0.00000 0.0312500 0.00000 0.0312500 
diff --git a/tests/feinterface/fe_interface_values_04.cc b/tests/feinterface/fe_interface_values_04.cc
new file mode 100644 (file)
index 0000000..f7c2acb
--- /dev/null
@@ -0,0 +1,104 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate jump(), average(), shape_value() of FEInterfaceValues for
+// at_boundary()
+
+#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/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(unsigned int fe_degree)
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+
+  DoFHandler<dim> dofh(tria);
+  FE_DGQ<dim>     fe(fe_degree);
+  deallog << fe.get_name() << std::endl;
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+
+  auto               cell = dofh.begin();
+  const unsigned int f    = 2;
+  fiv.reinit(cell, f);
+
+  const unsigned int n_dofs = fiv.n_current_interface_dofs();
+  Vector<double>     cell_vector(n_dofs);
+
+  const auto &q_points = fiv.get_quadrature_points();
+
+  cell_vector = 0.0;
+  for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+    for (unsigned int i = 0; i < n_dofs; ++i)
+      cell_vector(i) +=
+        fiv.shape_value(true, i, qpoint) * fiv.get_JxW_values()[qpoint];
+  deallog << "shape_value(true): " << cell_vector;
+
+  cell_vector = 0.0;
+  for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+    for (unsigned int i = 0; i < n_dofs; ++i)
+      cell_vector(i) +=
+        fiv.shape_value(false, i, qpoint) * fiv.get_JxW_values()[qpoint];
+  deallog << "shape_value(false): " << cell_vector;
+
+  cell_vector = 0.0;
+  for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+    for (unsigned int i = 0; i < n_dofs; ++i)
+      cell_vector(i) += fiv.jump(i, qpoint) * fiv.get_JxW_values()[qpoint];
+  deallog << "jump(): " << cell_vector;
+
+  cell_vector = 0.0;
+  for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+    for (unsigned int i = 0; i < n_dofs; ++i)
+      cell_vector(i) += fiv.average(i, qpoint) * fiv.get_JxW_values()[qpoint];
+  deallog << "average(): " << cell_vector;
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>(0);
+  test<2>(1);
+  test<3>(0);
+  test<3>(1);
+}
diff --git a/tests/feinterface/fe_interface_values_04.output b/tests/feinterface/fe_interface_values_04.output
new file mode 100644 (file)
index 0000000..d7c239c
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::FE_DGQ<2>(0)
+DEAL::shape_value(true): 1.00000 
+DEAL::shape_value(false): 0.00000 
+DEAL::jump(): 1.00000 
+DEAL::average(): 1.00000 
+DEAL::FE_DGQ<2>(1)
+DEAL::shape_value(true): 0.500000 0.500000 0.00000 0.00000 
+DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 
+DEAL::jump(): 0.500000 0.500000 0.00000 0.00000 
+DEAL::average(): 0.500000 0.500000 0.00000 0.00000 
+DEAL::FE_DGQ<3>(0)
+DEAL::shape_value(true): 1.00000 
+DEAL::shape_value(false): 0.00000 
+DEAL::jump(): 1.00000 
+DEAL::average(): 1.00000 
+DEAL::FE_DGQ<3>(1)
+DEAL::shape_value(true): 0.250000 0.250000 0.00000 0.00000 0.250000 0.250000 0.00000 0.00000 
+DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::jump(): 0.250000 0.250000 0.00000 0.00000 0.250000 0.250000 0.00000 0.00000 
+DEAL::average(): 0.250000 0.250000 0.00000 0.00000 0.250000 0.250000 0.00000 0.00000 
diff --git a/tests/feinterface/fe_interface_values_05.cc b/tests/feinterface/fe_interface_values_05.cc
new file mode 100644 (file)
index 0000000..5525c2e
--- /dev/null
@@ -0,0 +1,172 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// continuous element, otherwise like fe_interface_values_01
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+inspect_fiv(FEInterfaceValues<dim> &fiv)
+{
+  deallog << "at_boundary(): " << fiv.at_boundary() << "\n"
+          << "n_current_interface_dofs(): " << fiv.n_current_interface_dofs()
+          << "\n";
+
+  std::vector<types::global_dof_index> indices =
+    fiv.get_interface_dof_indices();
+  Assert(indices.size() == fiv.n_current_interface_dofs(), ExcInternalError());
+
+  deallog << "interface_dof_indices: ";
+  for (auto i : indices)
+    deallog << i << " ";
+  deallog << "\n";
+
+
+  unsigned int idx = 0;
+  for (auto v : indices)
+    {
+      deallog << "  index " << idx << " global_dof_index:" << v << ":\n";
+
+      const auto pair = fiv.interface_dof_to_dof_indices(idx);
+      deallog << "    dof indices: " << pair[0] << " | " << pair[1] << "\n";
+
+      ++idx;
+    }
+
+  deallog << "\n";
+  deallog << std::endl;
+}
+
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test()
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(1);
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+
+  auto cell = dofh.begin();
+
+  deallog << "** interface between cell 0 and 1 **\n";
+
+  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    if (!cell->at_boundary(f))
+      {
+        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());
+        Assert(!fiv.at_boundary(), ExcInternalError());
+
+        auto mycell = cell;
+        for (unsigned int c = 0; c < 2; ++c)
+          {
+            std::vector<types::global_dof_index> indices(fe.n_dofs_per_cell());
+            mycell->get_dof_indices(indices);
+            deallog << "cell " << c << ": ";
+            for (auto i : indices)
+              deallog << i << " ";
+            deallog << "\n";
+            ++mycell;
+          }
+
+        inspect_fiv(fiv);
+      }
+
+  deallog << "** boundary interface on cell 1 **\n";
+
+  {
+    ++cell;
+    fiv.reinit(cell, 1);
+    Assert(fiv.get_fe_face_values(0).get_cell() == cell, ExcInternalError());
+    Assert(fiv.n_current_interface_dofs() == fe.n_dofs_per_cell(),
+           ExcInternalError());
+    Assert(fiv.at_boundary(), ExcInternalError());
+    inspect_fiv(fiv);
+  }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>();
+}
diff --git a/tests/feinterface/fe_interface_values_05.output b/tests/feinterface/fe_interface_values_05.output
new file mode 100644 (file)
index 0000000..f7814e1
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3 
+cell 1: 1 4 3 5 
+at_boundary(): 0
+n_current_interface_dofs(): 6
+interface_dof_indices: 0 1 2 3 4 5 
+  index 0 global_dof_index:0:
+    dof indices: 0 | 4294967295
+  index 1 global_dof_index:1:
+    dof indices: 1 | 0
+  index 2 global_dof_index:2:
+    dof indices: 2 | 4294967295
+  index 3 global_dof_index:3:
+    dof indices: 3 | 2
+  index 4 global_dof_index:4:
+    dof indices: 4294967295 | 1
+  index 5 global_dof_index:5:
+    dof indices: 4294967295 | 3
+
+
+DEAL::** boundary interface on cell 1 **
+at_boundary(): 1
+n_current_interface_dofs(): 4
+interface_dof_indices: 1 4 3 5 
+  index 0 global_dof_index:1:
+    dof indices: 0 | 4294967295
+  index 1 global_dof_index:4:
+    dof indices: 1 | 4294967295
+  index 2 global_dof_index:3:
+    dof indices: 2 | 4294967295
+  index 3 global_dof_index:5:
+    dof indices: 3 | 4294967295
+
+
diff --git a/tests/feinterface/fe_interface_values_06.cc b/tests/feinterface/fe_interface_values_06.cc
new file mode 100644 (file)
index 0000000..8a0b49c
--- /dev/null
@@ -0,0 +1,157 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate jump(), average(), shape_value() of FEInterfaceValues on an adaptive
+// mesh for a continuous element
+
+#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/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(unsigned int fe_degree)
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+  tria.begin()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(fe_degree);
+  deallog << fe.get_name() << std::endl;
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+  auto cell = dofh.begin(1);
+  ++cell;
+
+  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    if (!cell->at_boundary(f))
+      {
+        if (!cell->neighbor_is_coarser(f))
+          continue;
+
+        auto nn = cell->neighbor_of_coarser_neighbor(f);
+        fiv.reinit(cell,
+                   f,
+                   numbers::invalid_unsigned_int,
+                   cell->neighbor(f),
+                   nn.first,
+                   nn.second);
+
+        const unsigned int n_dofs = fiv.n_current_interface_dofs();
+        Vector<double>     cell_vector(n_dofs);
+
+        const auto &q_points = fiv.get_quadrature_points();
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          deallog << "qpoint " << qpoint << ": " << q_points[qpoint]
+                  << std::endl;
+
+        for (unsigned int idx = 0; idx < n_dofs; ++idx)
+          {
+            const auto pair = fiv.interface_dof_to_dof_indices(idx);
+            deallog << "  idx: " << idx
+                    << " global: " << fiv.get_interface_dof_indices()[idx]
+                    << " dof indices: " << pair[0] << " | " << pair[1]
+                    << std::endl;
+          }
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(true, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(true): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(false, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(false): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.jump(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "jump(): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.average(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "average(): " << cell_vector;
+      }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>(1);
+  test<3>(1);
+}
diff --git a/tests/feinterface/fe_interface_values_06.output b/tests/feinterface/fe_interface_values_06.output
new file mode 100644 (file)
index 0000000..aba856b
--- /dev/null
@@ -0,0 +1,39 @@
+
+DEAL::FE_Q<2>(1)
+DEAL::qpoint 0: 1.00000 0.105662
+DEAL::qpoint 1: 1.00000 0.394338
+DEAL::  idx: 0 global: 0 dof indices: 1 | 0
+DEAL::  idx: 1 global: 1 dof indices: 4294967295 | 1
+DEAL::  idx: 2 global: 2 dof indices: 4294967295 | 2
+DEAL::  idx: 3 global: 3 dof indices: 4294967295 | 3
+DEAL::  idx: 4 global: 5 dof indices: 0 | 4294967295
+DEAL::  idx: 5 global: 7 dof indices: 2 | 4294967295
+DEAL::  idx: 6 global: 8 dof indices: 3 | 4294967295
+DEAL::shape_value(true): 0.250000 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 
+DEAL::shape_value(false): 0.375000 0.00000 0.125000 0.00000 0.00000 0.00000 0.00000 
+DEAL::jump(): -0.125000 0.00000 -0.125000 0.00000 0.00000 0.00000 0.250000 
+DEAL::average(): 0.312500 0.00000 0.0625000 0.00000 0.00000 0.00000 0.125000 
+DEAL::FE_Q<3>(1)
+DEAL::qpoint 0: 1.00000 0.105662 0.105662
+DEAL::qpoint 1: 1.00000 0.394338 0.105662
+DEAL::qpoint 2: 1.00000 0.105662 0.394338
+DEAL::qpoint 3: 1.00000 0.394338 0.394338
+DEAL::  idx: 0 global: 0 dof indices: 1 | 0
+DEAL::  idx: 1 global: 1 dof indices: 4294967295 | 1
+DEAL::  idx: 2 global: 2 dof indices: 4294967295 | 2
+DEAL::  idx: 3 global: 3 dof indices: 4294967295 | 3
+DEAL::  idx: 4 global: 4 dof indices: 4294967295 | 4
+DEAL::  idx: 5 global: 5 dof indices: 4294967295 | 5
+DEAL::  idx: 6 global: 6 dof indices: 4294967295 | 6
+DEAL::  idx: 7 global: 7 dof indices: 4294967295 | 7
+DEAL::  idx: 8 global: 9 dof indices: 0 | 4294967295
+DEAL::  idx: 9 global: 11 dof indices: 2 | 4294967295
+DEAL::  idx: 10 global: 13 dof indices: 4 | 4294967295
+DEAL::  idx: 11 global: 15 dof indices: 6 | 4294967295
+DEAL::  idx: 12 global: 16 dof indices: 3 | 4294967295
+DEAL::  idx: 13 global: 17 dof indices: 5 | 4294967295
+DEAL::  idx: 14 global: 18 dof indices: 7 | 4294967295
+DEAL::shape_value(true): 0.0625000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.0625000 0.0625000 0.0625000 
+DEAL::shape_value(false): 0.140625 0.00000 0.0468750 0.00000 0.0468750 0.00000 0.0156250 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::jump(): -0.0781250 0.00000 -0.0468750 0.00000 -0.0468750 0.00000 -0.0156250 0.00000 0.00000 0.00000 0.00000 0.00000 0.0625000 0.0625000 0.0625000 
+DEAL::average(): 0.101562 0.00000 0.0234375 0.00000 0.0234375 0.00000 0.00781250 0.00000 0.00000 0.00000 0.00000 0.00000 0.0312500 0.0312500 0.0312500 
diff --git a/tests/feinterface/fe_interface_values_07.cc b/tests/feinterface/fe_interface_values_07.cc
new file mode 100644 (file)
index 0000000..4e71510
--- /dev/null
@@ -0,0 +1,139 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate jump(), average(), shape_value() of FEInterfaceValues
+// like fe_interface_values_02, but for a continuous element
+
+#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/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(unsigned int fe_degree)
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(fe_degree);
+  deallog << fe.get_name() << std::endl;
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+
+  auto cell = dofh.begin();
+
+  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    if (!cell->at_boundary(f))
+      {
+        fiv.reinit(cell,
+                   f,
+                   numbers::invalid_unsigned_int,
+                   cell->neighbor(f),
+                   cell->neighbor_of_neighbor(f),
+                   numbers::invalid_unsigned_int);
+
+        const unsigned int n_dofs = fiv.n_current_interface_dofs();
+        Vector<double>     cell_vector(n_dofs);
+
+        const auto &q_points = fiv.get_quadrature_points();
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(true, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(true): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(false, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(false): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.jump(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "jump(): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.average(i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "average(): " << cell_vector;
+      }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>(1);
+  test<3>(1);
+}
diff --git a/tests/feinterface/fe_interface_values_07.output b/tests/feinterface/fe_interface_values_07.output
new file mode 100644 (file)
index 0000000..ba20331
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::FE_Q<2>(1)
+DEAL::shape_value(true): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 
+DEAL::shape_value(false): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 
+DEAL::jump(): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::average(): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 
+DEAL::FE_Q<3>(1)
+DEAL::shape_value(true): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 
+DEAL::shape_value(false): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 
+DEAL::jump(): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::average(): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 
diff --git a/tests/feinterface/step-12.cc b/tests/feinterface/step-12.cc
new file mode 100644 (file)
index 0000000..3ea3efd
--- /dev/null
@@ -0,0 +1,546 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 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.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Guido Kanschat, Texas A&M University, 2009
+ * Timo Heister, Clemson University, 2019
+ */
+
+// a version of step-12 using FEInterfaceValues
+
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/derivative_approximation.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step12
+{
+  using namespace dealii;
+
+
+
+  template <int dim>
+  struct ScratchData
+  {
+    ScratchData(const Mapping<dim> &      mapping,
+                const FiniteElement<dim> &fe,
+                const unsigned int        quadrature_degree,
+                const UpdateFlags         update_flags = update_values |
+                                                 update_gradients |
+                                                 update_quadrature_points |
+                                                 update_JxW_values,
+                const UpdateFlags interface_update_flags =
+                  update_values | update_gradients | update_quadrature_points |
+                  update_JxW_values | update_normal_vectors)
+      : fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags)
+      , fe_interface_values(mapping,
+                            fe,
+                            QGauss<dim - 1>(quadrature_degree),
+                            interface_update_flags)
+    {}
+
+
+    ScratchData(const ScratchData<dim> &scratch_data)
+      : fe_values(scratch_data.fe_values.get_mapping(),
+                  scratch_data.fe_values.get_fe(),
+                  scratch_data.fe_values.get_quadrature(),
+                  scratch_data.fe_values.get_update_flags())
+      , fe_interface_values(
+          scratch_data.fe_values
+            .get_mapping(), // TODO: implement for fe_interface_values
+          scratch_data.fe_values.get_fe(),
+          scratch_data.fe_interface_values.get_quadrature(),
+          scratch_data.fe_interface_values.get_update_flags())
+    {}
+
+    FEValues<dim>          fe_values;
+    FEInterfaceValues<dim> fe_interface_values;
+  };
+
+
+
+  struct CopyDataFace
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> joint_dof_indices;
+  };
+
+
+
+  struct CopyData
+  {
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+    std::vector<CopyDataFace>            face_data;
+
+    template <class Iterator>
+    void
+    reinit(const Iterator &cell, unsigned int dofs_per_cell)
+    {
+      cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+      cell_rhs.reinit(dofs_per_cell);
+
+      local_dof_indices.resize(dofs_per_cell);
+      cell->get_dof_indices(local_dof_indices);
+    }
+  };
+
+
+
+  template <class MatrixType, class VectorType>
+  inline void
+  copy(const CopyData &                 c,
+       const AffineConstraints<double> &constraints,
+       MatrixType &                     system_matrix,
+       VectorType &                     system_rhs)
+  {
+    constraints.distribute_local_to_global(c.cell_matrix,
+                                           c.cell_rhs,
+                                           c.local_dof_indices,
+                                           system_matrix,
+                                           system_rhs);
+    for (auto &cdf : c.face_data)
+      {
+        const unsigned int dofs_per_cell = cdf.joint_dof_indices.size();
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          for (unsigned int k = 0; k < dofs_per_cell; ++k)
+            system_matrix.add(cdf.joint_dof_indices[i],
+                              cdf.joint_dof_indices[k],
+                              cdf.cell_matrix(i, k));
+      }
+  }
+
+
+
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    BoundaryValues() = default;
+    virtual void
+    value_list(const std::vector<Point<dim>> &points,
+               std::vector<double> &          values,
+               const unsigned int             component = 0) const override;
+  };
+
+  template <int dim>
+  void
+  BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
+                                  std::vector<double> &          values,
+                                  const unsigned int component) const
+  {
+    (void)component;
+    AssertIndexRange(component, 1);
+    Assert(values.size() == points.size(),
+           ExcDimensionMismatch(values.size(), points.size()));
+
+    for (unsigned int i = 0; i < values.size(); ++i)
+      {
+        if (points[i](0) < 0.5)
+          values[i] = 1.;
+        else
+          values[i] = 0.;
+      }
+  }
+
+
+  template <int dim>
+  Tensor<1, dim>
+  beta(const Point<dim> &p)
+  {
+    Assert(dim >= 2, ExcNotImplemented());
+
+    Point<dim> wind_field;
+    wind_field(0) = -p(1);
+    wind_field(1) = p(0);
+    wind_field /= wind_field.norm();
+
+    return wind_field;
+  }
+
+
+
+  template <int dim>
+  class AdvectionProblem
+  {
+  public:
+    AdvectionProblem();
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+    void
+    output_results(const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    const MappingQ1<dim> mapping;
+
+    FE_DGQ<dim>     fe;
+    DoFHandler<dim> dof_handler;
+
+    AffineConstraints<> constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double> solution;
+    Vector<double> right_hand_side;
+  };
+
+
+  template <int dim>
+  AdvectionProblem<dim>::AdvectionProblem()
+    : mapping()
+    , fe(1)
+    , dof_handler(triangulation)
+  {}
+
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+    sparsity_pattern.copy_from(dsp);
+
+    system_matrix.reinit(sparsity_pattern);
+    solution.reinit(dof_handler.n_dofs());
+    right_hand_side.reinit(dof_handler.n_dofs());
+  }
+
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::assemble_system()
+  {
+    typedef decltype(dof_handler.begin_active()) Iterator;
+    BoundaryValues<dim>                          boundary_function;
+
+    auto cell_worker = [&](const Iterator &  cell,
+                           ScratchData<dim> &scratch_data,
+                           CopyData &        copy_data) {
+      const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell;
+      copy_data.reinit(cell, n_dofs);
+      scratch_data.fe_values.reinit(cell);
+
+      const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+      const FEValues<dim> &      fe_v = scratch_data.fe_values;
+      const std::vector<double> &JxW  = fe_v.get_JxW_values();
+
+      for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+        {
+          auto beta_q = beta(q_points[point]);
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                copy_data.cell_matrix(i, j) +=
+                  -beta_q                      // -\beta
+                  * fe_v.shape_grad(i, point)  // \nabla \phi_i
+                  * fe_v.shape_value(j, point) // \phi_j
+                  * JxW[point];                // dx
+              }
+        }
+    };
+
+    auto boundary_worker = [&](const Iterator &    cell,
+                               const unsigned int &face_no,
+                               ScratchData<dim> &  scratch_data,
+                               CopyData &          copy_data) {
+      scratch_data.fe_interface_values.reinit(cell, face_no);
+      const FEFaceValuesBase<dim> &fe_face =
+        scratch_data.fe_interface_values.get_fe_face_values(0);
+
+      const auto &q_points = fe_face.get_quadrature_points();
+
+      const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
+      const std::vector<double> &        JxW     = fe_face.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
+
+      std::vector<double> g(q_points.size());
+      boundary_function.value_list(q_points, g);
+
+      for (unsigned int point = 0; point < q_points.size(); ++point)
+        {
+          const double beta_n = beta(q_points[point]) * normals[point];
+
+          if (beta_n > 0)
+            {
+              for (unsigned int i = 0; i < n_facet_dofs; ++i)
+                for (unsigned int j = 0; j < n_facet_dofs; ++j)
+                  copy_data.cell_matrix(i, j) +=
+                    fe_face.shape_value(i, point)   // \phi_i
+                    * fe_face.shape_value(j, point) // \phi_j
+                    * beta_n                        // \beta . n
+                    * JxW[point];                   // dx
+            }
+          else
+            for (unsigned int i = 0; i < n_facet_dofs; ++i)
+              copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
+                                       * g[point]                     // g
+                                       * beta_n      // \beta . n
+                                       * JxW[point]; // dx
+        }
+    };
+
+    auto face_worker = [&](const Iterator &    cell,
+                           const unsigned int &f,
+                           const unsigned int &sf,
+                           const Iterator &    ncell,
+                           const unsigned int &nf,
+                           const unsigned int &nsf,
+                           ScratchData<dim> &  scratch_data,
+                           CopyData &          copy_data) {
+      FEInterfaceValues<dim> &fe_facet = scratch_data.fe_interface_values;
+      fe_facet.reinit(cell, f, sf, ncell, nf, nsf);
+      const auto &q_points = fe_facet.get_quadrature_points();
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+      const unsigned int n_dofs        = fe_facet.n_current_interface_dofs();
+      copy_data_face.joint_dof_indices = fe_facet.get_interface_dof_indices();
+
+      copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+      const std::vector<double> &        JxW = fe_facet.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals =
+        fe_facet.get_normal_vectors();
+
+      for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+        {
+          const double beta_n = beta(q_points[qpoint]) * normals[qpoint];
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              copy_data_face.cell_matrix(i, j) +=
+                fe_facet.jump(i, qpoint)                        // [\phi_i]
+                * fe_facet.shape_value((beta_n > 0), j, qpoint) // phi_j^{UP}
+                * beta_n                                        // (\beta . n)
+                * JxW[qpoint];                                  // dx
+        }
+    };
+
+    auto copier = [&](const CopyData &c) {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             system_matrix,
+                                             right_hand_side);
+
+      for (auto &cdf : c.face_data)
+        {
+          constraints.distribute_local_to_global(cdf.cell_matrix,
+                                                 cdf.joint_dof_indices,
+                                                 system_matrix);
+        }
+    };
+
+    const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
+
+    ScratchData<dim> scratch_data(mapping, fe, n_gauss_points);
+    CopyData         copy_data;
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          copy_data,
+                          MeshWorker::assemble_own_cells |
+                            MeshWorker::assemble_boundary_faces |
+                            MeshWorker::assemble_own_interior_faces_once,
+                          boundary_worker,
+                          face_worker);
+  }
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::solve()
+  {
+    SolverControl      solver_control(1000, 1e-12, false, false);
+    SolverRichardson<> solver(solver_control);
+
+    PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
+
+    preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, right_hand_side, preconditioner),
+      solver_control.last_step(),
+      4,
+      17);
+  }
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::refine_grid()
+  {
+    Vector<float> gradient_indicator(triangulation.n_active_cells());
+
+    DerivativeApproximation::approximate_gradient(mapping,
+                                                  dof_handler,
+                                                  solution,
+                                                  gradient_indicator);
+
+    unsigned int cell_no = 0;
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      gradient_indicator(cell_no++) *=
+        std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
+
+    deallog << "gradient_indicator l-infty: "
+            << gradient_indicator.linfty_norm() << std::endl;
+
+
+    GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+                                                    gradient_indicator,
+                                                    0.3,
+                                                    0.1);
+
+    triangulation.execute_coarsening_and_refinement();
+  }
+
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::output_results(const unsigned int cycle) const
+  {
+    const std::string filename = "solution-" + std::to_string(cycle) + ".vtk";
+    deallog << "Writing solution to <" << filename << ">" << std::endl;
+    std::ofstream output(filename);
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "u");
+
+    data_out.build_patches();
+
+    data_out.write_vtk(output);
+  }
+
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::run()
+  {
+    for (unsigned int cycle = 0; cycle < 3; ++cycle)
+      {
+        deallog << "Cycle " << cycle << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(3);
+          }
+        else
+          refine_grid();
+
+        deallog << "Number of active cells:       "
+                << triangulation.n_active_cells() << std::endl;
+
+        setup_system();
+
+        deallog << "Number of degrees of freedom: " << dof_handler.n_dofs()
+                << std::endl;
+
+        assemble_system();
+        solve();
+
+        output_results(cycle);
+      }
+  }
+} // namespace Step12
+
+
+int
+main()
+{
+  initlog();
+
+  try
+    {
+      Step12::AdvectionProblem<2> dgmethod;
+      dgmethod.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
diff --git a/tests/feinterface/step-12.output b/tests/feinterface/step-12.output
new file mode 100644 (file)
index 0000000..87d0505
--- /dev/null
@@ -0,0 +1,18 @@
+
+DEAL::Cycle 0
+DEAL::Number of active cells:       64
+DEAL::Number of degrees of freedom: 256
+DEAL::Solver stopped within 4 - 17 iterations
+DEAL::Writing solution to <solution-0.vtk>
+DEAL::Cycle 1
+DEAL::gradient_indicator l-infty: 0.160661
+DEAL::Number of active cells:       112
+DEAL::Number of degrees of freedom: 448
+DEAL::Solver stopped within 4 - 17 iterations
+DEAL::Writing solution to <solution-1.vtk>
+DEAL::Cycle 2
+DEAL::gradient_indicator l-infty: 0.0785475
+DEAL::Number of active cells:       214
+DEAL::Number of degrees of freedom: 856
+DEAL::Solver stopped within 4 - 17 iterations
+DEAL::Writing solution to <solution-2.vtk>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.