--- /dev/null
+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)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+DEAL_II_PICKUP_TESTS()
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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
+
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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
+
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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 ©_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;
+}
--- /dev/null
+
+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>
projects / dealii.git / commitdiff