From 7fb474eaa881abafb89a25bf3582ae625d027c37 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Thu, 8 Dec 2022 12:03:27 +0100 Subject: [PATCH] Add hp to FEInterfaceValues --- include/deal.II/fe/fe_interface_values.h | 357 +++++++++++++++--- tests/feinterface/fe_interface_values_11.cc | 196 ++++++++++ .../feinterface/fe_interface_values_11.output | 228 +++++++++++ tests/feinterface/fe_interface_values_12.cc | 173 +++++++++ .../feinterface/fe_interface_values_12.output | 57 +++ 5 files changed, 949 insertions(+), 62 deletions(-) create mode 100644 tests/feinterface/fe_interface_values_11.cc create mode 100644 tests/feinterface/fe_interface_values_11.output create mode 100644 tests/feinterface/fe_interface_values_12.cc create mode 100644 tests/feinterface/fe_interface_values_12.output diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index 52f430d44e..70287b4661 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -23,6 +23,8 @@ #include #include +#include +#include #include DEAL_II_NAMESPACE_OPEN @@ -1430,6 +1432,23 @@ public: const Quadrature & quadrature, const UpdateFlags update_flags); + /** + * Construct the FEInterfaceValues object with different FiniteElements + * assigned to different faces. + */ + FEInterfaceValues( + const hp::MappingCollection &mapping_collection, + const hp::FECollection & fe_collection, + const hp::QCollection & quadrature_collection, + const UpdateFlags update_flags); + + /** + * Same as above, but using a Q1 mapping. + */ + FEInterfaceValues(const hp::FECollection &fe_collection, + const hp::QCollection &quadrature_collection, + 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 @@ -2117,22 +2136,23 @@ private: /** * The FEFaceValues object for the current cell. */ - FEFaceValues internal_fe_face_values; + std::unique_ptr> internal_fe_face_values; /** * The FEFaceValues object for the current cell if the cell is refined. */ - FESubfaceValues internal_fe_subface_values; + std::unique_ptr> internal_fe_subface_values; /** * The FEFaceValues object for the neighboring cell. */ - FEFaceValues internal_fe_face_values_neighbor; + std::unique_ptr> internal_fe_face_values_neighbor; /** * The FEFaceValues object for the neighboring cell if the cell is refined. */ - FESubfaceValues internal_fe_subface_values_neighbor; + std::unique_ptr> + internal_fe_subface_values_neighbor; /** * Pointer to internal_fe_face_values or internal_fe_subface_values, @@ -2147,8 +2167,37 @@ private: */ FEFaceValuesBase *fe_face_values_neighbor; - /* Make the view classes friends of this class, since they access internal - * data. + /** + * An hp::FEValues object for the FEFaceValues on the + * present cell. + */ + std::unique_ptr> internal_hp_fe_face_values; + + /** + * An hp::FEValues object for the FESubfaceValues on the + * present cell. + */ + std::unique_ptr> + internal_hp_fe_subface_values; + + /** + * An hp::FEValues object for the FEFaceValues on the + * neighbor of the present cell. + */ + std::unique_ptr> + internal_hp_fe_face_values_neighbor; + + /** + * An hp::FEValues object for the FESubfaceValues on the + * neighboring cell. + */ + std::unique_ptr> + internal_hp_fe_subface_values_neighbor; + + + /* + * Make the view classes friends of this class, since they + * access internal data. */ template friend class FEInterfaceViews::Scalar; @@ -2169,10 +2218,26 @@ FEInterfaceValues::FEInterfaceValues( const Quadrature & quadrature, const UpdateFlags update_flags) : n_quadrature_points(quadrature.size()) - , 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) + , internal_fe_face_values( + std::make_unique>(mapping, + fe, + quadrature, + update_flags)) + , internal_fe_subface_values( + std::make_unique>(mapping, + fe, + quadrature, + update_flags)) + , internal_fe_face_values_neighbor( + std::make_unique>(mapping, + fe, + quadrature, + update_flags)) + , internal_fe_subface_values_neighbor( + std::make_unique>(mapping, + fe, + quadrature, + update_flags)) , fe_face_values(nullptr) , fe_face_values_neighbor(nullptr) {} @@ -2186,45 +2251,138 @@ FEInterfaceValues::FEInterfaceValues( const hp::QCollection & quadrature, const UpdateFlags update_flags) : n_quadrature_points(quadrature.max_n_quadrature_points()) - , 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[0], update_flags) - , internal_fe_subface_values_neighbor(mapping, - fe, - quadrature[0], - update_flags) + + , internal_fe_face_values( + std::make_unique>(mapping, + fe, + quadrature, + update_flags)) + , internal_fe_subface_values( + std::make_unique>(mapping, + fe, + quadrature, + update_flags)) + , internal_fe_face_values_neighbor( + std::make_unique>(mapping, + fe, + quadrature[0], + update_flags)) + , internal_fe_subface_values_neighbor( + std::make_unique>(mapping, + fe, + quadrature[0], + update_flags)) , fe_face_values(nullptr) , fe_face_values_neighbor(nullptr) {} +template +FEInterfaceValues::FEInterfaceValues( + const hp::MappingCollection &mapping_collection, + const hp::FECollection & fe_collection, + const hp::QCollection & quadrature_collection, + const UpdateFlags update_flags) + : n_quadrature_points(quadrature_collection.max_n_quadrature_points()) + , fe_face_values(nullptr) + , fe_face_values_neighbor(nullptr) + , internal_hp_fe_face_values( + std::make_unique>(mapping_collection, + fe_collection, + quadrature_collection, + update_flags)) + , internal_hp_fe_subface_values( + std::make_unique>( + mapping_collection, + fe_collection, + quadrature_collection, + update_flags)) + , internal_hp_fe_face_values_neighbor( + std::make_unique>(mapping_collection, + fe_collection, + quadrature_collection, + update_flags)) + , internal_hp_fe_subface_values_neighbor( + std::make_unique>( + mapping_collection, + fe_collection, + quadrature_collection, + update_flags)) +{ + AssertDimension(dim, spacedim); +} + + + +template +FEInterfaceValues::FEInterfaceValues( + const hp::FECollection &fe_collection, + const hp::QCollection & quadrature_collection, + const UpdateFlags update_flags) + : n_quadrature_points(quadrature_collection.max_n_quadrature_points()) + , fe_face_values(nullptr) + , fe_face_values_neighbor(nullptr) + , internal_hp_fe_face_values( + std::make_unique>( + fe_collection.get_reference_cell_default_linear_mapping(), + fe_collection, + quadrature_collection, + update_flags)) + , internal_hp_fe_subface_values( + std::make_unique>( + fe_collection.get_reference_cell_default_linear_mapping(), + fe_collection, + quadrature_collection, + update_flags)) + , internal_hp_fe_face_values_neighbor( + std::make_unique>( + fe_collection.get_reference_cell_default_linear_mapping(), + fe_collection, + quadrature_collection, + update_flags)) + , internal_hp_fe_subface_values_neighbor( + std::make_unique>( + fe_collection.get_reference_cell_default_linear_mapping(), + fe_collection, + quadrature_collection, + update_flags)) +{ + AssertDimension(dim, spacedim); +} + + + template FEInterfaceValues::FEInterfaceValues( const FiniteElement &fe, const Quadrature & quadrature, const UpdateFlags update_flags) : n_quadrature_points(quadrature.size()) - , internal_fe_face_values( + , internal_fe_face_values(std::make_unique>( fe.reference_cell().template get_default_linear_mapping(), fe, quadrature, - update_flags) - , internal_fe_subface_values( + update_flags)) + , internal_fe_subface_values(std::make_unique>( fe.reference_cell().template get_default_linear_mapping(), fe, quadrature, - update_flags) + update_flags)) , internal_fe_face_values_neighbor( - fe.reference_cell().template get_default_linear_mapping(), - fe, - quadrature, - update_flags) + std::make_unique>( + fe.reference_cell() + .template get_default_linear_mapping(), + fe, + quadrature, + update_flags)) , internal_fe_subface_values_neighbor( - fe.reference_cell().template get_default_linear_mapping(), - fe, - quadrature, - update_flags) + std::make_unique>( + fe.reference_cell() + .template get_default_linear_mapping(), + fe, + quadrature, + update_flags)) , fe_face_values(nullptr) , fe_face_values_neighbor(nullptr) {} @@ -2242,34 +2400,94 @@ FEInterfaceValues::reinit( 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) + Assert(internal_fe_face_values || internal_hp_fe_face_values, + ExcNotInitialized()); + + if (internal_fe_face_values) { - internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor); - fe_face_values_neighbor = &internal_fe_face_values_neighbor; + if (sub_face_no == numbers::invalid_unsigned_int) + { + internal_fe_face_values->reinit(cell, face_no); + fe_face_values = internal_fe_face_values.get(); + } + else + { + internal_fe_subface_values->reinit(cell, face_no, sub_face_no); + fe_face_values = internal_fe_subface_values.get(); + } + 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.get(); + } + 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.get(); + } + + AssertDimension(fe_face_values->n_quadrature_points, + fe_face_values_neighbor->n_quadrature_points); + + const_cast(this->n_quadrature_points) = + fe_face_values->n_quadrature_points; } - else + else if (internal_hp_fe_face_values) { - internal_fe_subface_values_neighbor.reinit(cell_neighbor, - face_no_neighbor, - sub_face_no_neighbor); - fe_face_values_neighbor = &internal_fe_subface_values_neighbor; - } + const unsigned int dominated_fe_index = + internal_hp_fe_face_values->get_fe_collection().find_dominated_fe( + {cell->active_fe_index(), cell_neighbor->active_fe_index()}); + + // Same as if above, but when hp is enabled. + if (sub_face_no == numbers::invalid_unsigned_int) + { + internal_hp_fe_face_values->reinit(cell, + face_no, + dominated_fe_index, + dominated_fe_index); + fe_face_values = &const_cast &>( + internal_hp_fe_face_values->get_present_fe_values()); + } + else + { + internal_hp_fe_subface_values->reinit( + cell, face_no, sub_face_no, dominated_fe_index, dominated_fe_index); + + fe_face_values = &const_cast &>( + internal_hp_fe_subface_values->get_present_fe_values()); + } + if (sub_face_no_neighbor == numbers::invalid_unsigned_int) + { + internal_hp_fe_face_values_neighbor->reinit(cell_neighbor, + face_no_neighbor, + dominated_fe_index, + dominated_fe_index); + + fe_face_values_neighbor = &const_cast &>( + internal_hp_fe_face_values_neighbor->get_present_fe_values()); + } + else + { + internal_hp_fe_subface_values_neighbor->reinit(cell_neighbor, + face_no_neighbor, + sub_face_no_neighbor, + dominated_fe_index, + dominated_fe_index); + + fe_face_values_neighbor = + &const_cast &>( + internal_hp_fe_subface_values_neighbor->get_present_fe_values()); + } - AssertDimension(fe_face_values->n_quadrature_points, - fe_face_values_neighbor->n_quadrature_points); + AssertDimension(fe_face_values->n_quadrature_points, + fe_face_values_neighbor->n_quadrature_points); - const_cast(this->n_quadrature_points) = - fe_face_values->n_quadrature_points; + const_cast(this->n_quadrature_points) = + fe_face_values->n_quadrature_points; + } // Set up dof mapping and remove duplicates (for continuous elements). { @@ -2323,15 +2541,30 @@ 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; + Assert(internal_fe_face_values || internal_hp_fe_face_values, + ExcNotInitialized()); - interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell()); - cell->get_active_or_mg_dof_indices(interface_dof_indices); + if (internal_fe_face_values) + { + internal_fe_face_values->reinit(cell, face_no); + fe_face_values = internal_fe_face_values.get(); + fe_face_values_neighbor = nullptr; - dofmap.resize(interface_dof_indices.size()); + interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell()); + cell->get_active_or_mg_dof_indices(interface_dof_indices); + } + else if (internal_hp_fe_face_values) + { + internal_hp_fe_face_values->reinit(cell, face_no); + fe_face_values = &const_cast &>( + internal_hp_fe_face_values->get_present_fe_values()); + fe_face_values_neighbor = nullptr; + interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell()); + cell->get_active_or_mg_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}}; @@ -2377,7 +2610,7 @@ template const Mapping & FEInterfaceValues::get_mapping() const { - return internal_fe_face_values.get_mapping(); + return internal_fe_face_values->get_mapping(); } @@ -2386,7 +2619,7 @@ template const FiniteElement & FEInterfaceValues::get_fe() const { - return internal_fe_face_values.get_fe(); + return internal_fe_face_values->get_fe(); } @@ -2395,7 +2628,7 @@ template const Quadrature & FEInterfaceValues::get_quadrature() const { - return internal_fe_face_values.get_quadrature(); + return internal_fe_face_values->get_quadrature(); } @@ -2424,7 +2657,7 @@ template UpdateFlags FEInterfaceValues::get_update_flags() const { - return internal_fe_face_values.get_update_flags(); + return internal_fe_face_values->get_update_flags(); } diff --git a/tests/feinterface/fe_interface_values_11.cc b/tests/feinterface/fe_interface_values_11.cc new file mode 100644 index 0000000000..52b09de88c --- /dev/null +++ b/tests/feinterface/fe_interface_values_11.cc @@ -0,0 +1,196 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2022 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 in the hp case. + +#include + +#include +#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: " << static_cast(pair[0]) << " | " + << static_cast(pair[1]) << "\n"; + + ++idx; + } + + 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(const unsigned int p) +{ + Triangulation tria; + make_2_cells(tria); + + DoFHandler dofh(tria); + hp::FECollection fe_collection; + hp::QCollection q_collection; + + fe_collection.push_back(FE_DGQ(p)); + fe_collection.push_back(FE_DGQ(p + 1)); + + q_collection.push_back(QGauss(p)); + q_collection.push_back(QGauss(p + 1)); + + // Set different finite elements spaces on the two cells. + unsigned int fe_index = 0; + for (const auto &cell : dofh.active_cell_iterators()) + { + cell->set_active_fe_index(fe_index); + ++fe_index; + } + + dofh.distribute_dofs(fe_collection); + + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(fe_collection, q_collection, update_flags); + + + auto cell = dofh.begin(); + + deallog << "** interface between cell 0 and 1 **\n"; + + for (const unsigned int f : GeometryInfo::face_indices()) + 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() == + fe_collection[0].n_dofs_per_cell() + + fe_collection[1].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_collection[c].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 0 **\n" << std::endl; + + { + fiv.reinit(cell, 0); + Assert(fiv.get_fe_face_values(0).get_cell() == cell, ExcInternalError()); + Assert(fiv.n_current_interface_dofs() == fe_collection[0].n_dofs_per_cell(), + ExcInternalError()); + Assert(fiv.at_boundary(), ExcInternalError()); + inspect_fiv(fiv); + } +} + + + +int +main() +{ + initlog(); + for (unsigned int p : {1, 2}) + test<2>(p); + for (unsigned int p : {1}) + test<3>(p); +} diff --git a/tests/feinterface/fe_interface_values_11.output b/tests/feinterface/fe_interface_values_11.output new file mode 100644 index 0000000000..081798974b --- /dev/null +++ b/tests/feinterface/fe_interface_values_11.output @@ -0,0 +1,228 @@ + +DEAL::** interface between cell 0 and 1 ** +cell 0: 0 1 2 3 +cell 1: 4 5 6 7 8 9 10 11 12 +at_boundary(): 0 +n_current_interface_dofs(): 13 +interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 + index 0 global_dof_index:0: + dof indices: 0 | -1 + index 1 global_dof_index:1: + dof indices: 1 | -1 + index 2 global_dof_index:2: + dof indices: 2 | -1 + index 3 global_dof_index:3: + dof indices: 3 | -1 + index 4 global_dof_index:4: + dof indices: -1 | 0 + index 5 global_dof_index:5: + dof indices: -1 | 1 + index 6 global_dof_index:6: + dof indices: -1 | 2 + index 7 global_dof_index:7: + dof indices: -1 | 3 + index 8 global_dof_index:8: + dof indices: -1 | 4 + index 9 global_dof_index:9: + dof indices: -1 | 5 + index 10 global_dof_index:10: + dof indices: -1 | 6 + index 11 global_dof_index:11: + dof indices: -1 | 7 + index 12 global_dof_index:12: + dof indices: -1 | 8 + +DEAL::** boundary interface on cell 0 ** + +DEAL::at_boundary(): 1 +n_current_interface_dofs(): 4 +interface_dof_indices: 0 1 2 3 + index 0 global_dof_index:0: + dof indices: 0 | -1 + index 1 global_dof_index:1: + dof indices: 1 | -1 + index 2 global_dof_index:2: + dof indices: 2 | -1 + index 3 global_dof_index:3: + dof indices: 3 | -1 + +DEAL::** interface between cell 0 and 1 ** +cell 0: 0 1 2 3 4 5 6 7 8 +cell 1: 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 +at_boundary(): 0 +n_current_interface_dofs(): 25 +interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 + index 0 global_dof_index:0: + dof indices: 0 | -1 + index 1 global_dof_index:1: + dof indices: 1 | -1 + index 2 global_dof_index:2: + dof indices: 2 | -1 + index 3 global_dof_index:3: + dof indices: 3 | -1 + index 4 global_dof_index:4: + dof indices: 4 | -1 + index 5 global_dof_index:5: + dof indices: 5 | -1 + index 6 global_dof_index:6: + dof indices: 6 | -1 + index 7 global_dof_index:7: + dof indices: 7 | -1 + index 8 global_dof_index:8: + dof indices: 8 | -1 + index 9 global_dof_index:9: + dof indices: -1 | 0 + index 10 global_dof_index:10: + dof indices: -1 | 1 + index 11 global_dof_index:11: + dof indices: -1 | 2 + index 12 global_dof_index:12: + dof indices: -1 | 3 + index 13 global_dof_index:13: + dof indices: -1 | 4 + index 14 global_dof_index:14: + dof indices: -1 | 5 + index 15 global_dof_index:15: + dof indices: -1 | 6 + index 16 global_dof_index:16: + dof indices: -1 | 7 + index 17 global_dof_index:17: + dof indices: -1 | 8 + index 18 global_dof_index:18: + dof indices: -1 | 9 + index 19 global_dof_index:19: + dof indices: -1 | 10 + index 20 global_dof_index:20: + dof indices: -1 | 11 + index 21 global_dof_index:21: + dof indices: -1 | 12 + index 22 global_dof_index:22: + dof indices: -1 | 13 + index 23 global_dof_index:23: + dof indices: -1 | 14 + index 24 global_dof_index:24: + dof indices: -1 | 15 + +DEAL::** boundary interface on cell 0 ** + +DEAL::at_boundary(): 1 +n_current_interface_dofs(): 9 +interface_dof_indices: 0 1 2 3 4 5 6 7 8 + index 0 global_dof_index:0: + dof indices: 0 | -1 + index 1 global_dof_index:1: + dof indices: 1 | -1 + index 2 global_dof_index:2: + dof indices: 2 | -1 + index 3 global_dof_index:3: + dof indices: 3 | -1 + index 4 global_dof_index:4: + dof indices: 4 | -1 + index 5 global_dof_index:5: + dof indices: 5 | -1 + index 6 global_dof_index:6: + dof indices: 6 | -1 + index 7 global_dof_index:7: + dof indices: 7 | -1 + index 8 global_dof_index:8: + dof indices: 8 | -1 + +DEAL::** interface between cell 0 and 1 ** +cell 0: 0 1 2 3 4 5 6 7 +cell 1: 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 +at_boundary(): 0 +n_current_interface_dofs(): 35 +interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 + index 0 global_dof_index:0: + dof indices: 0 | -1 + index 1 global_dof_index:1: + dof indices: 1 | -1 + index 2 global_dof_index:2: + dof indices: 2 | -1 + index 3 global_dof_index:3: + dof indices: 3 | -1 + index 4 global_dof_index:4: + dof indices: 4 | -1 + index 5 global_dof_index:5: + dof indices: 5 | -1 + index 6 global_dof_index:6: + dof indices: 6 | -1 + index 7 global_dof_index:7: + dof indices: 7 | -1 + index 8 global_dof_index:8: + dof indices: -1 | 0 + index 9 global_dof_index:9: + dof indices: -1 | 1 + index 10 global_dof_index:10: + dof indices: -1 | 2 + index 11 global_dof_index:11: + dof indices: -1 | 3 + index 12 global_dof_index:12: + dof indices: -1 | 4 + index 13 global_dof_index:13: + dof indices: -1 | 5 + index 14 global_dof_index:14: + dof indices: -1 | 6 + index 15 global_dof_index:15: + dof indices: -1 | 7 + index 16 global_dof_index:16: + dof indices: -1 | 8 + index 17 global_dof_index:17: + dof indices: -1 | 9 + index 18 global_dof_index:18: + dof indices: -1 | 10 + index 19 global_dof_index:19: + dof indices: -1 | 11 + index 20 global_dof_index:20: + dof indices: -1 | 12 + index 21 global_dof_index:21: + dof indices: -1 | 13 + index 22 global_dof_index:22: + dof indices: -1 | 14 + index 23 global_dof_index:23: + dof indices: -1 | 15 + index 24 global_dof_index:24: + dof indices: -1 | 16 + index 25 global_dof_index:25: + dof indices: -1 | 17 + index 26 global_dof_index:26: + dof indices: -1 | 18 + index 27 global_dof_index:27: + dof indices: -1 | 19 + index 28 global_dof_index:28: + dof indices: -1 | 20 + index 29 global_dof_index:29: + dof indices: -1 | 21 + index 30 global_dof_index:30: + dof indices: -1 | 22 + index 31 global_dof_index:31: + dof indices: -1 | 23 + index 32 global_dof_index:32: + dof indices: -1 | 24 + index 33 global_dof_index:33: + dof indices: -1 | 25 + index 34 global_dof_index:34: + dof indices: -1 | 26 + +DEAL::** boundary interface on cell 0 ** + +DEAL::at_boundary(): 1 +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 | -1 + index 1 global_dof_index:1: + dof indices: 1 | -1 + index 2 global_dof_index:2: + dof indices: 2 | -1 + index 3 global_dof_index:3: + dof indices: 3 | -1 + index 4 global_dof_index:4: + dof indices: 4 | -1 + index 5 global_dof_index:5: + dof indices: 5 | -1 + index 6 global_dof_index:6: + dof indices: 6 | -1 + index 7 global_dof_index:7: + dof indices: 7 | -1 + diff --git a/tests/feinterface/fe_interface_values_12.cc b/tests/feinterface/fe_interface_values_12.cc new file mode 100644 index 0000000000..e32d1f0a22 --- /dev/null +++ b/tests/feinterface/fe_interface_values_12.cc @@ -0,0 +1,173 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2021 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_in_shape_values(), average_of_shape_values(), shape_value() of +// FEInterfaceValues on an adaptive mesh in the hp scenario. + +#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(const unsigned int fe_degree0, const unsigned int fe_degree1 = 0) +{ + Triangulation tria; + make_2_cells(tria); + + tria.begin()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + + + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + DoFHandler dofh(tria); + hp::FECollection fe_collection; + hp::QCollection q_collection; + fe_collection.push_back(FE_DGQ(fe_degree0)); + fe_collection.push_back(FE_DGQ(fe_degree1)); + + q_collection.push_back(QGauss(fe_degree0 + 1)); + // Set different finite elements spaces on the two cells. + unsigned int fe_index = 0; + for (const auto &cell : dofh.active_cell_iterators_on_level(0)) + { + cell->set_active_fe_index(fe_index); + ++fe_index; + } + + deallog << fe_collection[0].get_name() << "-" << fe_collection[1].get_name() + << std::endl; + dofh.distribute_dofs(fe_collection); + + FEInterfaceValues fiv(fe_collection, q_collection, update_flags); + + auto cell = dofh.begin(1); + ++cell; + + for (const unsigned int f : GeometryInfo::face_indices()) + 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: " << static_cast(pair[0]) << " | " + << static_cast(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 << 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(false, i, qpoint) * fiv.get_JxW_values()[qpoint]; + deallog << "shape_value(false): " << cell_vector << 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.jump_in_shape_values(i, qpoint) * + fiv.get_JxW_values()[qpoint]; + deallog << "jump_in_shape_values(): " << cell_vector << 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.average_of_shape_values(i, qpoint) * + fiv.get_JxW_values()[qpoint]; + deallog << "average_of_shape_values(): " << cell_vector << std::endl; + } +} + + + +int +main() +{ + initlog(); + test<2>(0); + test<2>(1, 2); + test<3>(0); + test<3>(1, 2); +} diff --git a/tests/feinterface/fe_interface_values_12.output b/tests/feinterface/fe_interface_values_12.output new file mode 100644 index 0000000000..93d406c647 --- /dev/null +++ b/tests/feinterface/fe_interface_values_12.output @@ -0,0 +1,57 @@ + +DEAL::FE_DGQ<2>(0)-FE_DGQ<2>(0) +DEAL::qpoint 0: 1.00000 0.250000 +DEAL:: idx: 0 global: 0 dof indices: -1 | 0 +DEAL:: idx: 1 global: 2 dof indices: 0 | -1 +DEAL::shape_value(true): 0.00000 0.500000 +DEAL::shape_value(false): 0.500000 0.00000 +DEAL::jump_in_shape_values(): -0.500000 0.500000 +DEAL::average_of_shape_values(): 0.250000 0.250000 +DEAL::FE_DGQ<2>(1)-FE_DGQ<2>(2) +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: -1 | 1 +DEAL:: idx: 2 global: 2 dof indices: -1 | 2 +DEAL:: idx: 3 global: 3 dof indices: -1 | 3 +DEAL:: idx: 4 global: 8 dof indices: 0 | -1 +DEAL:: idx: 5 global: 9 dof indices: 1 | -1 +DEAL:: idx: 6 global: 10 dof indices: 2 | -1 +DEAL:: idx: 7 global: 11 dof indices: 3 | -1 +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_in_shape_values(): -0.375000 0.00000 -0.125000 0.00000 0.00000 0.250000 0.00000 0.250000 +DEAL::average_of_shape_values(): 0.187500 0.00000 0.0625000 0.00000 0.00000 0.125000 0.00000 0.125000 +DEAL::FE_DGQ<3>(0)-FE_DGQ<3>(0) +DEAL::qpoint 0: 1.00000 0.250000 0.250000 +DEAL:: idx: 0 global: 0 dof indices: -1 | 0 +DEAL:: idx: 1 global: 2 dof indices: 0 | -1 +DEAL::shape_value(true): 0.00000 0.250000 +DEAL::shape_value(false): 0.250000 0.00000 +DEAL::jump_in_shape_values(): -0.250000 0.250000 +DEAL::average_of_shape_values(): 0.125000 0.125000 +DEAL::FE_DGQ<3>(1)-FE_DGQ<3>(2) +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: -1 | 1 +DEAL:: idx: 2 global: 2 dof indices: -1 | 2 +DEAL:: idx: 3 global: 3 dof indices: -1 | 3 +DEAL:: idx: 4 global: 4 dof indices: -1 | 4 +DEAL:: idx: 5 global: 5 dof indices: -1 | 5 +DEAL:: idx: 6 global: 6 dof indices: -1 | 6 +DEAL:: idx: 7 global: 7 dof indices: -1 | 7 +DEAL:: idx: 8 global: 16 dof indices: 0 | -1 +DEAL:: idx: 9 global: 17 dof indices: 1 | -1 +DEAL:: idx: 10 global: 18 dof indices: 2 | -1 +DEAL:: idx: 11 global: 19 dof indices: 3 | -1 +DEAL:: idx: 12 global: 20 dof indices: 4 | -1 +DEAL:: idx: 13 global: 21 dof indices: 5 | -1 +DEAL:: idx: 14 global: 22 dof indices: 6 | -1 +DEAL:: idx: 15 global: 23 dof indices: 7 | -1 +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_in_shape_values(): -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_of_shape_values(): 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 -- 2.39.5