From 2e0523933fdf8bcaa27a8b81270b58b0657036cb Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sat, 24 Aug 2019 22:01:14 -0400 Subject: [PATCH] introduce FEInterfaceValues --- doc/news/changes/major/20190824Heister | 5 + include/deal.II/fe/fe_interface_values.h | 694 ++++++++++++++++++ tests/feinterface/CMakeLists.txt | 4 + tests/feinterface/fe_interface_values_01.cc | 174 +++++ .../feinterface/fe_interface_values_01.output | 39 + tests/feinterface/fe_interface_values_02.cc | 140 ++++ .../feinterface/fe_interface_values_02.output | 21 + tests/feinterface/fe_interface_values_03.cc | 159 ++++ .../feinterface/fe_interface_values_03.output | 57 ++ tests/feinterface/fe_interface_values_04.cc | 104 +++ .../feinterface/fe_interface_values_04.output | 21 + tests/feinterface/fe_interface_values_05.cc | 172 +++++ .../feinterface/fe_interface_values_05.output | 35 + tests/feinterface/fe_interface_values_06.cc | 157 ++++ .../feinterface/fe_interface_values_06.output | 39 + tests/feinterface/fe_interface_values_07.cc | 139 ++++ .../feinterface/fe_interface_values_07.output | 11 + tests/feinterface/step-12.cc | 546 ++++++++++++++ tests/feinterface/step-12.output | 18 + 19 files changed, 2535 insertions(+) create mode 100644 doc/news/changes/major/20190824Heister create mode 100644 include/deal.II/fe/fe_interface_values.h create mode 100644 tests/feinterface/CMakeLists.txt create mode 100644 tests/feinterface/fe_interface_values_01.cc create mode 100644 tests/feinterface/fe_interface_values_01.output create mode 100644 tests/feinterface/fe_interface_values_02.cc create mode 100644 tests/feinterface/fe_interface_values_02.output create mode 100644 tests/feinterface/fe_interface_values_03.cc create mode 100644 tests/feinterface/fe_interface_values_03.output create mode 100644 tests/feinterface/fe_interface_values_04.cc create mode 100644 tests/feinterface/fe_interface_values_04.output create mode 100644 tests/feinterface/fe_interface_values_05.cc create mode 100644 tests/feinterface/fe_interface_values_05.output create mode 100644 tests/feinterface/fe_interface_values_06.cc create mode 100644 tests/feinterface/fe_interface_values_06.output create mode 100644 tests/feinterface/fe_interface_values_07.cc create mode 100644 tests/feinterface/fe_interface_values_07.output create mode 100644 tests/feinterface/step-12.cc create mode 100644 tests/feinterface/step-12.output diff --git a/doc/news/changes/major/20190824Heister b/doc/news/changes/major/20190824Heister new file mode 100644 index 0000000000..dd3d0ebaa7 --- /dev/null +++ b/doc/news/changes/major/20190824Heister @@ -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. +
+(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 index 0000000000..fef1461668 --- /dev/null +++ b/include/deal.II/fe/fe_interface_values.h @@ -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 +#include + +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 +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 & mapping, + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags update_flags); + + /** + * Construct the FEInterfaceValues with a single FiniteElement and + * a Q1 Mapping. + * + * See the constructor above. + */ + FEInterfaceValues(const FiniteElement &fe, + const Quadrature & 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 + 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 + 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 & + get_fe_face_values(const unsigned int cell_index) const; + + /** + * Return a reference to the quadrature object in use. + */ + const Quadrature & + 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 & + 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> & + get_normal_vectors() const; + + /** + * Return a reference to the quadrature points in real space. + * + * @dealiiRequiresUpdateFlags{update_quadrature_points} + */ + const std::vector> & + 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 + 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 + 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 + * (direction * normal)>0 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 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> dofmap; + + /** + * The FEFaceValues object for the current cell. + */ + FEFaceValues internal_fe_face_values; + + /** + * The FEFaceValues object for the current cell if the cell is refined. + */ + FESubfaceValues internal_fe_subface_values; + + /** + * The FEFaceValues object for the neighboring cell. + */ + FEFaceValues internal_fe_face_values_neighbor; + + /** + * The FEFaceValues object for the neighboring cell if the cell is refined. + */ + FESubfaceValues internal_fe_subface_values_neighbor; + + /** + * Pointer to internal_fe_face_values or internal_fe_subface_values, + * respectively as determined in reinit(). + */ + FEFaceValuesBase *fe_face_values; + + /** + * Pointer to internal_fe_face_values_neighbor, + * internal_fe_subface_values_neighbor, or nullptr, respectively + * as determined in reinit(). + */ + FEFaceValuesBase *fe_face_values_neighbor; +}; + + + +#ifndef DOXYGEN + +/*---------------------- Inline functions ---------------------*/ + +template +FEInterfaceValues::FEInterfaceValues( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature & 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 +FEInterfaceValues::FEInterfaceValues( + const FiniteElement &fe, + const Quadrature & quadrature, + const UpdateFlags update_flags) + : internal_fe_face_values(StaticMappingQ1::mapping, + fe, + quadrature, + update_flags) + , internal_fe_subface_values(StaticMappingQ1::mapping, + fe, + quadrature, + update_flags) + , internal_fe_face_values_neighbor(StaticMappingQ1::mapping, + fe, + quadrature, + update_flags) + , internal_fe_subface_values_neighbor(StaticMappingQ1::mapping, + fe, + quadrature, + update_flags) + , fe_face_values(nullptr) + , fe_face_values_neighbor(nullptr) +{} + + + +template +template +void +FEInterfaceValues::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 v( + fe_face_values->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(v); + std::vector 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> + tempmap; + std::pair 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 +template +void +FEInterfaceValues::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 +const std::vector & +FEInterfaceValues::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 +const std::vector> & +FEInterfaceValues::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 +const Quadrature & +FEInterfaceValues::get_quadrature() const +{ + return internal_fe_face_values.get_quadrature(); +} + + + +template +const std::vector> & +FEInterfaceValues::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 +const UpdateFlags +FEInterfaceValues::get_update_flags() const +{ + return internal_fe_face_values.get_update_flags(); +} + + + +template +unsigned +FEInterfaceValues::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 +bool +FEInterfaceValues::at_boundary() const +{ + return fe_face_values_neighbor == nullptr; +} + + + +template +std::vector +FEInterfaceValues::get_interface_dof_indices() const +{ + return interface_dof_indices; +} + + + +template +std::array +FEInterfaceValues::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 +const FEFaceValuesBase & +FEInterfaceValues::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 +Tensor<1, spacedim> +FEInterfaceValues::normal(const unsigned int q_point_index) const +{ + return fe_face_values->normal_vector(q_point_index); +} + + + +template +double +FEInterfaceValues::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 +double +FEInterfaceValues::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 +double +FEInterfaceValues::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 index 0000000000..d33eafa20b --- /dev/null +++ b/tests/feinterface/CMakeLists.txt @@ -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 index 0000000000..4565e209cc --- /dev/null +++ b/tests/feinterface/fe_interface_values_01.cc @@ -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 + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + + +template +void +inspect_fiv(FEInterfaceValues &fiv) +{ + deallog << "at_boundary(): " << fiv.at_boundary() << "\n" + << "n_current_interface_dofs(): " << fiv.n_current_interface_dofs() + << "\n"; + + std::vector 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 +void +make_2_cells(Triangulation &tria); + +template <> +void make_2_cells<2>(Triangulation<2> &tria) +{ + const unsigned int dim = 2; + std::vector repetitions = {2, 1}; + Point p1; + Point 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 repetitions = {2, 1, 1}; + Point p1; + Point p2(2.0, 1.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + + +template +void +test() +{ + Triangulation tria; + make_2_cells(tria); + + DoFHandler dofh(tria); + FE_DGQ fe(1); + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + + auto cell = dofh.begin(); + + deallog << "** interface between cell 0 and 1 **\n"; + + for (unsigned int f = 0; f < GeometryInfo::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 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 index 0000000000..263e929668 --- /dev/null +++ b/tests/feinterface/fe_interface_values_01.output @@ -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 index 0000000000..b6c9a794d6 --- /dev/null +++ b/tests/feinterface/fe_interface_values_02.cc @@ -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 + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +make_2_cells(Triangulation &tria); + +template <> +void make_2_cells<2>(Triangulation<2> &tria) +{ + const unsigned int dim = 2; + std::vector repetitions = {2, 1}; + Point p1; + Point 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 repetitions = {2, 1, 1}; + Point p1; + Point p2(2.0, 1.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + + +template +void +test(unsigned int fe_degree) +{ + Triangulation tria; + make_2_cells(tria); + + DoFHandler dofh(tria); + FE_DGQ fe(fe_degree); + deallog << fe.get_name() << std::endl; + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + + auto cell = dofh.begin(); + + for (unsigned int f = 0; f < GeometryInfo::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 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 index 0000000000..0a6b89e53b --- /dev/null +++ b/tests/feinterface/fe_interface_values_02.output @@ -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 index 0000000000..a649ada606 --- /dev/null +++ b/tests/feinterface/fe_interface_values_03.cc @@ -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 + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +make_2_cells(Triangulation &tria); + +template <> +void make_2_cells<2>(Triangulation<2> &tria) +{ + const unsigned int dim = 2; + std::vector repetitions = {2, 1}; + Point p1; + Point 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 repetitions = {2, 1, 1}; + Point p1; + Point p2(2.0, 1.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + + +template +void +test(unsigned int fe_degree) +{ + Triangulation tria; + make_2_cells(tria); + tria.begin()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + DoFHandler dofh(tria); + FE_DGQ fe(fe_degree); + deallog << fe.get_name() << std::endl; + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + auto cell = dofh.begin(1); + ++cell; + + for (unsigned int f = 0; f < GeometryInfo::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 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 index 0000000000..8be1c3fd7a --- /dev/null +++ b/tests/feinterface/fe_interface_values_03.output @@ -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 index 0000000000..f7c2acb7d0 --- /dev/null +++ b/tests/feinterface/fe_interface_values_04.cc @@ -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 + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + + +template +void +test(unsigned int fe_degree) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dofh(tria); + FE_DGQ fe(fe_degree); + deallog << fe.get_name() << std::endl; + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(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 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 index 0000000000..d7c239c56f --- /dev/null +++ b/tests/feinterface/fe_interface_values_04.output @@ -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 index 0000000000..5525c2e8ba --- /dev/null +++ b/tests/feinterface/fe_interface_values_05.cc @@ -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 + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + + +template +void +inspect_fiv(FEInterfaceValues &fiv) +{ + deallog << "at_boundary(): " << fiv.at_boundary() << "\n" + << "n_current_interface_dofs(): " << fiv.n_current_interface_dofs() + << "\n"; + + std::vector 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 +void +make_2_cells(Triangulation &tria); + +template <> +void make_2_cells<2>(Triangulation<2> &tria) +{ + const unsigned int dim = 2; + std::vector repetitions = {2, 1}; + Point p1; + Point 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 repetitions = {2, 1, 1}; + Point p1; + Point p2(2.0, 1.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + + +template +void +test() +{ + Triangulation tria; + make_2_cells(tria); + + DoFHandler dofh(tria); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + + auto cell = dofh.begin(); + + deallog << "** interface between cell 0 and 1 **\n"; + + for (unsigned int f = 0; f < GeometryInfo::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 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 index 0000000000..f7814e1ad5 --- /dev/null +++ b/tests/feinterface/fe_interface_values_05.output @@ -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 index 0000000000..8a0b49c48f --- /dev/null +++ b/tests/feinterface/fe_interface_values_06.cc @@ -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 + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +make_2_cells(Triangulation &tria); + +template <> +void make_2_cells<2>(Triangulation<2> &tria) +{ + const unsigned int dim = 2; + std::vector repetitions = {2, 1}; + Point p1; + Point 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 repetitions = {2, 1, 1}; + Point p1; + Point p2(2.0, 1.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + + +template +void +test(unsigned int fe_degree) +{ + Triangulation tria; + make_2_cells(tria); + tria.begin()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + DoFHandler dofh(tria); + FE_Q fe(fe_degree); + deallog << fe.get_name() << std::endl; + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + auto cell = dofh.begin(1); + ++cell; + + for (unsigned int f = 0; f < GeometryInfo::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 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 index 0000000000..aba856b592 --- /dev/null +++ b/tests/feinterface/fe_interface_values_06.output @@ -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 index 0000000000..4e7151039e --- /dev/null +++ b/tests/feinterface/fe_interface_values_07.cc @@ -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 + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +make_2_cells(Triangulation &tria); + +template <> +void make_2_cells<2>(Triangulation<2> &tria) +{ + const unsigned int dim = 2; + std::vector repetitions = {2, 1}; + Point p1; + Point 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 repetitions = {2, 1, 1}; + Point p1; + Point p2(2.0, 1.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + + +template +void +test(unsigned int fe_degree) +{ + Triangulation tria; + make_2_cells(tria); + + DoFHandler dofh(tria); + FE_Q fe(fe_degree); + deallog << fe.get_name() << std::endl; + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + + auto cell = dofh.begin(); + + for (unsigned int f = 0; f < GeometryInfo::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 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 index 0000000000..ba20331f05 --- /dev/null +++ b/tests/feinterface/fe_interface_values_07.output @@ -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 index 0000000000..3ea3efd69d --- /dev/null +++ b/tests/feinterface/step-12.cc @@ -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 +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +namespace Step12 +{ + using namespace dealii; + + + + template + struct ScratchData + { + ScratchData(const Mapping & mapping, + const FiniteElement &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(quadrature_degree), update_flags) + , fe_interface_values(mapping, + fe, + QGauss(quadrature_degree), + interface_update_flags) + {} + + + ScratchData(const ScratchData &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 fe_values; + FEInterfaceValues fe_interface_values; + }; + + + + struct CopyDataFace + { + FullMatrix cell_matrix; + std::vector joint_dof_indices; + }; + + + + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + + template + 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 + inline void + copy(const CopyData & c, + const AffineConstraints &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 + class BoundaryValues : public Function + { + public: + BoundaryValues() = default; + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const override; + }; + + template + void + BoundaryValues::value_list(const std::vector> &points, + std::vector & 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 + Tensor<1, dim> + beta(const Point &p) + { + Assert(dim >= 2, ExcNotImplemented()); + + Point wind_field; + wind_field(0) = -p(1); + wind_field(1) = p(0); + wind_field /= wind_field.norm(); + + return wind_field; + } + + + + template + 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 triangulation; + const MappingQ1 mapping; + + FE_DGQ fe; + DoFHandler dof_handler; + + AffineConstraints<> constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector right_hand_side; + }; + + + template + AdvectionProblem::AdvectionProblem() + : mapping() + , fe(1) + , dof_handler(triangulation) + {} + + + template + void + AdvectionProblem::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 + void + AdvectionProblem::assemble_system() + { + typedef decltype(dof_handler.begin_active()) Iterator; + BoundaryValues boundary_function; + + auto cell_worker = [&](const Iterator & cell, + ScratchData &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 & fe_v = scratch_data.fe_values; + const std::vector &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 & scratch_data, + CopyData & copy_data) { + scratch_data.fe_interface_values.reinit(cell, face_no); + const FEFaceValuesBase &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 & JxW = fe_face.get_JxW_values(); + const std::vector> &normals = fe_face.get_normal_vectors(); + + std::vector 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 & scratch_data, + CopyData & copy_data) { + FEInterfaceValues &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 & JxW = fe_facet.get_JxW_values(); + const std::vector> &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 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 + void + AdvectionProblem::solve() + { + SolverControl solver_control(1000, 1e-12, false, false); + SolverRichardson<> solver(solver_control); + + PreconditionBlockSSOR> 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 + void + AdvectionProblem::refine_grid() + { + Vector 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 + void + AdvectionProblem::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 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 + void + AdvectionProblem::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 index 0000000000..87d050540c --- /dev/null +++ b/tests/feinterface/step-12.output @@ -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 +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 +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 -- 2.39.5