diff --git a/CMakeLists.txt b/CMakeLists.txt index e50a4baef68aa81b4a0dc5e214faa2f3971feaac..347565bc405dbb35c21eebf7df653201647ec406 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,7 @@ # or GPL-2.0+ (http://opensource.org/licenses/gpl-license) # with "runtime exception" (http://www.dune-project.org/license.html) # Authors: -# Felix Schindler (2016) +# Felix Schindler (2016 - 2017) # Rene Milk (2016) # Tobias Leibner (2016) @@ -43,9 +43,20 @@ make_dependent_modules_sys_included() add_subdirectory(dune) add_subdirectory(doc) +if(dune-pybindxi_FOUND) + foreach(_file + dune/__init__.py + dune/xt/__init__.py + dune/xt/grid/__init__.py + ) + execute_process(COMMAND ln -s ${CMAKE_CURRENT_SOURCE_DIR}/python/${_file} ${CMAKE_CURRENT_BINARY_DIR}/${_file}) + endforeach() +endif() + # enable headercheck add_definitions("-DENABLE_HEADERCHECK=1") add_format(${CMAKE_CURRENT_SOURCE_DIR}) +add_tidy(${CMAKE_CURRENT_SOURCE_DIR}) # finalize the dune project, e.g. generating config.h etc. finalize_dune_project(GENERATE_CONFIG_H_CMAKE) diff --git a/dune.module b/dune.module index 59357cc70a7ddead29a0a0e54dc193b4da855fc9..adb482e9c53b051e880d17f9efcf5c4b4dcda7c5 100644 --- a/dune.module +++ b/dune.module @@ -11,4 +11,4 @@ Module: dune-xt-grid Version: 2.5.0-dev Maintainer: dune-xt-dev@listserv.uni-muenster.de Depends: dune-common (>= 2.5) dune-testtools (>= 2.5) dune-geometry (>= 2.5) dune-grid (>= 2.5) dune-xt-common (>=2.5) -Suggests: dune-pdelab dune-fem +Suggests: dune-alugrid (>= 2.5) dune-grid-glue dune-fem dune-pdelab dune-ug diff --git a/dune/xt/CMakeLists.txt b/dune/xt/CMakeLists.txt index 87395d7cc1063c82301c9c8693eed1cd1390ce27..624d57e776fa36c836433c19f1834573bfeb757c 100644 --- a/dune/xt/CMakeLists.txt +++ b/dune/xt/CMakeLists.txt @@ -5,13 +5,12 @@ # or GPL-2.0+ (http://opensource.org/licenses/gpl-license) # with "runtime exception" (http://www.dune-project.org/license.html) # Authors: -# Felix Schindler (2016) +# Felix Schindler (2016 - 2017) # Rene Milk (2016) # Tobias Leibner (2016) set(lib_dune_xt_grid_sources grid/boundaryinfo/interfaces.cc) dune_library_add_sources(dunextgrid SOURCES ${lib_dune_xt_grid_sources}) -include_directories(SYSTEM ${DUNE_XT_COMMON_TEST_DIR}/gtest) -add_subdirectory(grid/test EXCLUDE_FROM_ALL ) +add_subdirectory(grid) diff --git a/dune/xt/grid/CMakeLists.txt b/dune/xt/grid/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..5fa078538284e5c7d8bed7483636aac20aa5bfcf --- /dev/null +++ b/dune/xt/grid/CMakeLists.txt @@ -0,0 +1,23 @@ +# This file is part of the dune-xt-grid project: +# https://github.com/dune-community/dune-xt-grid +# The copyright lies with the authors of this file (see below). +# License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +# or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +# with "runtime exception" (http://www.dune-project.org/license.html) +# Authors: +# Felix Schindler (2016 - 2017) +# Rene Milk (2016) +# Tobias Leibner (2016) + +include_directories(SYSTEM ${DUNE_XT_COMMON_TEST_DIR}/gtest) +add_subdirectory(test EXCLUDE_FROM_ALL ) + +if(dune-pybindxi_FOUND) + dune_pybindxi_add_module(_grid bindings.cc) + target_link_dune_default_libraries(_grid) + if(DUNE_XT_WITH_PYTHON_BINDINGS) + add_custom_target(bindings ALL DEPENDS _grid) + else() + add_custom_target(bindings DEPENDS _grid) + endif() +endif() diff --git a/dune/xt/grid/bindings.cc b/dune/xt/grid/bindings.cc new file mode 100644 index 0000000000000000000000000000000000000000..209f279642cc86fff95eadc1f8e7ddf591151448 --- /dev/null +++ b/dune/xt/grid/bindings.cc @@ -0,0 +1,172 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2016 - 2017) + +#include "config.h" + +#if HAVE_DUNE_PYBINDXI + +#include <string> +#include <vector> + +#include <dune/pybindxi/pybind11.h> +#include <dune/pybindxi/stl.h> + +#include <dune/xt/common/timedlogging.hh> +#include <dune/xt/common/configuration.pbh> +#include <dune/xt/common/fvector.pbh> + +#include "grids.hh" +#include "layers.hh" +#include "intersection.hh" +#include "type_traits.hh" + +#include "boundaryinfo.pbh" +#include "gridprovider.pbh" +#include "walker.pbh" +#include "walker/apply-on.pbh" + +namespace py = pybind11; +using namespace pybind11::literals; + + +template <class I, bool required = true> +struct for_Grid_and_Intersection +{ + template <class G> + static void addbind(py::module& m, const std::string& grid_id, const std::string& id) + { + using namespace Dune::XT::Grid; + pybind11::class_<BoundaryInfo<I>>( + m, std::string("BoundaryInfo__" + grid_id + id).c_str(), std::string("BoundaryInfo__" + grid_id + id).c_str()); + bind_BoundaryInfo<AllDirichletBoundaryInfo<I>>(m, "AllDirichletBoundaryInfo__" + grid_id + id); + bind_BoundaryInfo<AllNeumannBoundaryInfo<I>>(m, "AllNeumannBoundaryInfo__" + grid_id + id); + bind_BoundaryInfo<NormalBasedBoundaryInfo<I>>(m, "NormalBasedBoundaryInfo__" + grid_id + id); + } +}; + +template <class I> +struct for_Grid_and_Intersection<I, false> +{ + template <class G> + static void addbind(py::module& /*m*/, const std::string& /*grid_id*/, const std::string& /*id*/) + { + } +}; + + +template <class G> +void addbind_for_Grid(py::module& m, const std::string& grid_id) +{ + using namespace Dune::XT; + using namespace Dune::XT::Grid; + + typedef typename Layer<G, Layers::level, Backends::view>::type LevelView; + typedef typename Layer<G, Layers::leaf, Backends::view>::type LeafView; + typedef typename Intersection<LeafView>::type FVI; + typedef typename Intersection<LevelView>::type LVI; +#if HAVE_DUNE_FEM + typedef typename Layer<G, Layers::level, Backends::part>::type LevelPart; + typedef typename Layer<G, Layers::leaf, Backends::part>::type LeafPart; + typedef typename Intersection<LeafPart>::type FPI; + typedef typename Intersection<LevelPart>::type LPI; +#endif + + bind_GridProvider<G>(m, grid_id); + bind_make_cube_grid<G>(m, grid_id); + + for_Grid_and_Intersection<FVI, true>::template addbind<G>(m, + grid_id, + (std::is_same<FVI, LVI>::value +#if HAVE_DUNE_FEM + && std::is_same<FVI, FPI>::value + && std::is_same<FVI, LPI>::value +#endif + ) + ? "" + : "_leaf_view"); + for_Grid_and_Intersection<LVI, !(std::is_same<LVI, FVI>::value)>::template addbind<G>(m, grid_id, "_level_view"); +#if HAVE_DUNE_FEM + for_Grid_and_Intersection<FPI, + !(std::is_same<FPI, FVI>::value + || std::is_same<FPI, LVI>::value)>::template addbind<G>(m, grid_id, "_leaf_part"); + for_Grid_and_Intersection<LPI, + !(std::is_same<LPI, FVI>::value || std::is_same<LPI, LVI>::value + || std::is_same<LPI, FPI>::value)>::template addbind<G>(m, grid_id, "_level_part"); +#endif // HAVE_DUNE_FEM + + +#define BIND(V, id) \ + bind_Walker<V>(m, grid_id + "_" + id); \ + \ + m.def(std::string(std::string("make_all_dirichlet_boundary_info__") + id).c_str(), \ + [](const GridProvider<typename extract_grid<V>::type>& /*grid_provider*/) { \ + return AllDirichletBoundaryInfo<typename Intersection<V>::Type>(); \ + }, \ + "grid_provider"_a); \ + m.def(std::string(std::string("make_all_neumann_boundary_info__") + id).c_str(), \ + [](const GridProvider<typename extract_grid<V>::type>& /*grid_provider*/) { \ + return AllNeumannBoundaryInfo<typename Intersection<V>::Type>(); \ + }, \ + "grid_provider"_a); \ + m.def(std::string(std::string("make_normalbased_boundary_info__") + id).c_str(), \ + [](const GridProvider<typename extract_grid<V>::type>& /*grid_provider*/, const Common::Configuration& cfg) { \ + return *NormalBasedBoundaryInfo<typename Intersection<V>::Type>::create(cfg); \ + }, \ + "grid_provider"_a, \ + "cfg"_a = normalbased_boundaryinfo_default_config()); \ + addbind_WhichIntersection<V>(m, grid_id, id); /* end BIND \ + */ + + BIND(LevelView, "level_view"); + BIND(LeafView, "leaf_view"); +#if HAVE_DUNE_FEM + BIND(LevelPart, "level_part"); + BIND(LeafPart, "leaf_part"); +#endif +#undef BIND +} // ... addbind_for_Grid(...) + + +PYBIND11_PLUGIN(_grid) +{ + py::module m("_grid", "dune-xt-grid"); + + py::module::import("dune.xt.common"); + + addbind_for_Grid<Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>>(m, "2d_cube_yaspgrid"); +#if HAVE_ALUGRID || HAVE_DUNE_ALUGRID + addbind_for_Grid<Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>>(m, "2d_simplex_aluconform"); +#endif +#if HAVE_UG + addbind_for_Grid<Dune::UGGrid<2>>(m, "2d_simplex_uggrid"); +#endif + + m.def("init_logger", + [](const ssize_t max_info_level, + const ssize_t max_debug_level, + const bool enable_warnings, + const bool enable_colors, + const std::string& info_color, + const std::string& debug_color, + const std::string& warning_color) { + Dune::XT::Common::TimedLogger().create( + max_info_level, max_debug_level, enable_warnings, enable_colors, info_color, debug_color, warning_color); + }, + "max_info_level"_a = -1, + "max_debug_level"_a = -1, + "enable_warnings"_a = true, + "enable_colors"_a = true, + "info_color"_a = "blue", + "debug_color"_a = "darkgray", + "warning_color"_a = "red"); + + return m.ptr(); +} + +#endif // HAVE_DUNE_PYBINDXI diff --git a/dune/xt/grid/boundaryinfo.pbh b/dune/xt/grid/boundaryinfo.pbh new file mode 100644 index 0000000000000000000000000000000000000000..857e07ed1f7f525bc8921b2dc252aeaf3ecc60ba --- /dev/null +++ b/dune/xt/grid/boundaryinfo.pbh @@ -0,0 +1,50 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2016 - 2017) + +#ifndef DUNE_XT_GRID_BOUNDARYINFO_PBH +#define DUNE_XT_GRID_BOUNDARYINFO_PBH +#if HAVE_DUNE_PYBINDXI + +#include <sstream> +#include <type_traits> + +#include <dune/pybindxi/pybind11.h> +#include <dune/pybindxi/operators.h> + +#include "gridprovider/provider.hh" +#include "boundaryinfo.hh" + +namespace Dune { +namespace XT { +namespace Grid { + + +template <class C> +pybind11::class_<C> bind_BoundaryInfo(pybind11::module& m, const std::string& id /*, const std::string& method_id*/) +{ + namespace py = pybind11; + using namespace pybind11::literals; + + py::class_<C, BoundaryInfo<typename C::IntersectionType>> c(m, id.c_str(), id.c_str(), py::metaclass()); + + c.def(py::init<>()); + + c.def_property_readonly_static("static_id", [](const C& /*self*/) { return C::static_id(); }); + c.def("__repr__", [id](const C& /*self*/) { return id + "(" + C::static_id() + ")"; }); + + return c; +} // ... bind_BoundaryInfo(...) + + +} // namespace Grid +} // namespace XT +} // namespace Dune + +#endif // HAVE_DUNE_PYBINDXI +#endif // DUNE_XT_GRID_BOUNDARYINFO_PBH diff --git a/dune/xt/grid/dd/glued.hh b/dune/xt/grid/dd/glued.hh new file mode 100644 index 0000000000000000000000000000000000000000..d2fff63922e5d188cffd9b39df5b8a6a20b7da53 --- /dev/null +++ b/dune/xt/grid/dd/glued.hh @@ -0,0 +1,1098 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2017) + +#ifndef DUNE_XT_GRID_DD_GLUED_HH +#define DUNE_XT_GRID_DD_GLUED_HH + +#include <array> +#include <algorithm> +#include <memory> +#include <map> + +#include <boost/numeric/conversion/cast.hpp> + +#include <dune/grid/common/rangegenerators.hh> + +#include <dune/grid-glue/extractors/codim1extractor.hh> +#include <dune/grid-glue/extractors/extractorpredicate.hh> +#include <dune/grid-glue/gridglue.hh> +#include <dune/grid-glue/merging/contactmerge.hh> + +#include <dune/xt/common/exceptions.hh> +#include <dune/xt/common/float_cmp.hh> +#include <dune/xt/common/ranges.hh> +#include <dune/xt/common/timedlogging.hh> +#include <dune/xt/grid/intersection.hh> +#include <dune/xt/grid/gridprovider/provider.hh> +#include <dune/xt/grid/gridprovider/cube.hh> +#include <dune/xt/grid/search.hh> + +namespace Dune { +namespace XT { +namespace Grid { +namespace DD { +namespace Exceptions { + + +class intersection_orientation_is_broken : public Dune::InvalidStateException +{ +}; + + +} // namespace Exceptions + + +#if HAVE_DUNE_GRID_GLUE + + +template <typename P0, typename P1> +size_t check_for_broken_coupling_intersections( + const GridGlue::GridGlue<P0, P1>& glue, + const typename P0::ctype& tolerance = 10 * XT::Common::FloatCmp::DefaultEpsilon<typename P0::ctype>::value()) +{ + const auto inside_grid_view = glue.template gridView<0>(); + size_t failures = 0; + // walk the coupling + const auto coupling_intersection_it_end = glue.template iend<0>(); + for (auto coupling_intersection_it = glue.template ibegin<0>(); + coupling_intersection_it != coupling_intersection_it_end; + ++coupling_intersection_it) { + const auto& coupling_intersection = *coupling_intersection_it; + const auto coupling_intersection_normal = coupling_intersection.centerUnitOuterNormal(); + const auto local_entity = coupling_intersection.inside(); + typename std::remove_const<decltype(coupling_intersection_normal)>::type local_intersection_normal(0.); + // find the intersection of the local inside entity that corresponds to the coupling intersection + size_t found = 0; + for (auto&& local_intersection : intersections(inside_grid_view, local_entity)) { + // the coupling intersection may be smaller than the local intersection, so check if all of the corners of the + // coupling intersection lie within this local intersection + int corners_inside = 0; + for (auto ii : XT::Common::value_range(coupling_intersection.geometry().corners())) + if (XT::Grid::contains(local_intersection, coupling_intersection.geometry().corner(ii))) + ++corners_inside; + if (corners_inside == coupling_intersection.geometry().corners()) { + // this is the one + ++found; + local_intersection_normal = local_intersection.centerUnitOuterNormal(); + } + } + if (found != 1) + DUNE_THROW(InvalidStateException, + "This should not happen!\n" + << "There were " + << found + << " local intersections which contain the coupling intersection, " + << "and there must not be more than one!"); + // now the expected normal is local_intersection_normal + // and we would like coupling_intersection_normal to point in the same direction + // since they have unit length, they should be identical + if ((local_intersection_normal - coupling_intersection_normal).infinity_norm() > tolerance) + ++failures; + } + return failures; +} // ... check_for_broken_coupling_intersections(...) + + +// forward +template <class MacroGridType, class LocalGridType> +class GluedVTKWriter; + + +template <class MacroGridImp, class LocalGridImp> +class Glued +{ + template <class G, bool anything = true> + struct allowed_macro_grid + { + static const bool value = true; + }; + + template <class G, bool anything = true> + struct allowed_local_grid + { + static const bool value = true; + }; + +#if HAVE_ALUGRID + template <class Comm, bool anything> + struct allowed_local_grid<ALUGrid<3, 3, simplex, conforming, Comm>, anything> + { + static const bool value = false; + }; + + template <class Comm, bool anything> + struct allowed_local_grid<ALUGrid<3, 3, cube, conforming, Comm>, anything> + { + static const bool value = false; + }; +#endif + + static_assert(allowed_macro_grid<MacroGridImp>::value, + "This macro grid is known to fail, enable on your onw risk by disabling this check!"); + static_assert(allowed_local_grid<LocalGridImp>::value, + "This local grid is known to fail, enable on your onw risk by disabling this check!"); + +public: + typedef MacroGridImp MacroGridType; + typedef XT::Grid::GridProvider<MacroGridType> MacroGridProviderType; + typedef typename MacroGridProviderType::LeafGridViewType MacroGridViewType; + typedef typename MacroGridType::template Codim<0>::Entity MacroEntityType; + + typedef LocalGridImp LocalGridType; + typedef XT::Grid::GridProvider<LocalGridType> LocalGridProviderType; + +#if HAVE_DUNE_FEM + typedef typename LocalGridProviderType::LevelGridPartType MicroGridPartType; +#endif + typedef typename LocalGridType::LevelGridView MicroGridViewType; + typedef typename MicroGridViewType::template Codim<0>::Entity MicroEntityType; + typedef typename MicroGridViewType::template Codim<0>::EntityPointer MicroEntityPointerType; + + typedef typename MacroGridType::ctype ctype; + static const size_t dimDomain = MacroGridType::dimension; + static const size_t dimWorld = MacroGridType::dimensionworld; + +private: + typedef typename LocalGridProviderType::LevelGridViewType LocalViewType; + typedef GridGlue::Codim1Extractor<LocalViewType> LocalExtractorType; + typedef GridGlue::GridGlue<LocalExtractorType, LocalExtractorType> GlueType; + + template <class GridView, class MacroIntersectionType> + class CouplingFaceDescriptor : public GridGlue::ExtractorPredicate<GridView, 1> + { + typedef typename GridView::Traits::template Codim<0>::Entity LocalEntityType; + typedef typename GridView::ctype ctype; + + public: + CouplingFaceDescriptor(const MacroIntersectionType& macro_intersection) + : macro_intersection_(macro_intersection) + { + } + + virtual bool contains(const LocalEntityType& element, unsigned int face) const override final + { + const auto local_intersection_ptr = element.template subEntity<1>(face); +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) + const auto& local_intersection = local_intersection_ptr; +#else + const auto& local_intersection = *local_intersection_ptr; +#endif + const auto& local_intersection_geometry = local_intersection.geometry(); + // Check if all corners of the local intersection lie within the macro intersection. + for (auto ii : XT::Common::value_range(local_intersection_geometry.corners())) + if (!XT::Grid::contains(macro_intersection_, local_intersection_geometry.corner(ii))) + return false; + return true; + } // ... contains(...) + + private: + const MacroIntersectionType& macro_intersection_; + }; // CouplingFaceDescriptor + + template <class GV, class MI> + static CouplingFaceDescriptor<GV, MI> create_descriptor(const GV& /*gv*/, const MI& mi) + { + return CouplingFaceDescriptor<GV, MI>(mi); + } + +public: + Glued(MacroGridProviderType& macro_grid_provider, + const size_t num_local_refinements = 0, + const bool prepare_glues = false, + const bool allow_for_broken_orientation_of_coupling_intersections = false, + const ctype& allowed_overlap = 10 * XT::Common::FloatCmp::DefaultEpsilon<ctype>::value()) + : macro_grid_(macro_grid_provider) + , allowed_overlap_(allowed_overlap) + , macro_leaf_view_(macro_grid_.leaf_view()) + , local_grids_(macro_leaf_view_.indexSet().size(0), nullptr) + , glues_(macro_leaf_view_.indexSet().size(0)) + { + setup_local_grids(); + if (num_local_refinements > 0) + for (auto& local_grid_provider : local_grids_) { + assert(local_grid_provider); + local_grid_provider->grid().globalRefine(boost::numeric_cast<int>(num_local_refinements)); + } + if (prepare_glues) + setup_glues(allow_for_broken_orientation_of_coupling_intersections); + } // Glued(...) + + const MacroGridViewType& macro_grid_view() const + { + return macro_leaf_view_; + } + + size_t num_subdomains() const + { + return macro_grid_view().indexSet().size(0); + } + + size_t subdomain(const MacroEntityType& macro_entity) const + { + assert(macro_leaf_view_.indexSet().contains(macro_entity)); + return macro_leaf_view_.indexSet().index(macro_entity); + } + + bool boundary(const MacroEntityType& macro_entity) const + { + assert(macro_leaf_view_.indexSet().contains(macro_entity)); + return macro_entity.hasBoundaryIntersections(); + } + + MicroGridViewType global_grid_view() + { + auto logger = XT::Common::TimedLogger().get("grid-multiscale.glued.global_grid_view"); + logger.warn() << "Requiring access to global micro grid!" << std::endl; + prepare_global_grid(); + return global_grid_->level_view(global_grid_->grid().maxLevel()); + } + +#if HAVE_DUNE_FEM + MicroGridPartType global_grid_part() + { + auto logger = XT::Common::TimedLogger().get("grid-multiscale.glued.global_grid_part"); + logger.warn() << "Requiring access to global micro grid!" << std::endl; + prepare_global_grid(); + return MicroGridPartType(global_grid_->grid(), global_grid_->grid().maxLevel()); + } +#endif // HAVE_DUNE_FEM + + const std::vector<std::vector<size_t>>& local_to_global_indices() + { + auto logger = XT::Common::TimedLogger().get("grid-multiscale.glued.local_to_global_indices"); + logger.warn() << "Requiring access to global micro grid!" << std::endl; + prepare_global_grid(); + return *local_to_global_indices_; + } + + MicroEntityPointerType local_to_global_entity(const size_t subd, const MicroEntityType& local_entity) + { + return local_to_global_entity(subd, + local_grids_[subd]->level_view(max_local_level(subd)).indexSet().index(local_entity)); + } + + MicroEntityPointerType local_to_global_entity(const size_t subd, const size_t local_entity_index) + { + auto logger = XT::Common::TimedLogger().get("grid-multiscale.glued.local_to_global_entity"); + logger.warn() << "Requiring access to global micro grid!" << std::endl; + prepare_global_grid(); + const size_t global_index_of_local_entity = local_to_global_indices_->operator[](subd)[local_entity_index]; + const auto global_micro_grid_view = global_grid_view(); + const auto entity_it_end = global_micro_grid_view.template end<0>(); + for (auto entity_it = global_micro_grid_view.template begin<0>(); entity_it != entity_it_end; ++entity_it) { + const auto& entity = *entity_it; + if (global_micro_grid_view.indexSet().index(entity) == global_index_of_local_entity) + return entity_it; + } + DUNE_THROW(XT::Common::Exceptions::wrong_input_given, + "subdomain: " << subd << "\n" + << "local_entity_index: " + << local_entity_index + << "\n" + << "global_index_of_local_entity: " + << global_index_of_local_entity); + return global_micro_grid_view.template begin<0>(); + } // ... local_to_global_entity(...) + + MicroEntityPointerType global_to_local_entity(const MicroEntityType& micro_entity) + { + prepare_global_grid(); + return global_to_local_entity(global_grid_view().indexSet().index(micro_entity)); + } + + MicroEntityPointerType global_to_local_entity(const size_t micro_entity_index) + { + auto logger = XT::Common::TimedLogger().get("grid-multiscale.glued.global_to_local_entity"); + logger.warn() << "Requiring access to global micro grid!" << std::endl; + prepare_global_grid(); + auto subdomain_and_local_entity_index = global_to_local_indices_->operator[](micro_entity_index); + const auto subd = subdomain_and_local_entity_index.first; + const auto local_index_of_global_entity = subdomain_and_local_entity_index.second; + const auto local_grid_view = local_grids_[subd]->level_view(max_local_level(subd)); + const auto entity_it_end = local_grid_view.template end<0>(); + for (auto entity_it = local_grid_view.template begin<0>(); entity_it != entity_it_end; ++entity_it) { + const auto& entity = *entity_it; + if (local_grid_view.indexSet().index(entity) == local_index_of_global_entity) + return MicroEntityPointerType(entity_it); + } + DUNE_THROW(XT::Common::Exceptions::wrong_input_given, + "micro_entity_index: " << micro_entity_index << "\n" + "subdomain: " + << subd + << "\n" + << "local_index_of_global_entity: " + << local_index_of_global_entity); + return MicroEntityPointerType(local_grid_view.template begin<0>()); + } // ... global_to_local_entity(...) + + const std::vector<std::pair<size_t, size_t>>& global_to_local_indices() + { + auto logger = XT::Common::TimedLogger().get("grid-multiscale.glued.global_to_local_indices"); + logger.warn() << "Requiring access to global micro grid!" << std::endl; + prepare_global_grid(); + assert(global_to_local_indices_); + return *global_to_local_indices_; + } + + const LocalGridProviderType& local_grid(const MacroEntityType& macro_entity) const + { + assert(macro_leaf_view_.indexSet().contains(macro_entity)); + return local_grid(macro_leaf_view_.indexSet().index(macro_entity)); + } + + LocalGridProviderType& local_grid(const MacroEntityType& macro_entity) + { + assert(macro_leaf_view_.indexSet().contains(macro_entity)); + return local_grid(macro_leaf_view_.indexSet().index(macro_entity)); + } + + LocalGridProviderType& local_grid(const size_t macro_entity_index) + { + return *(local_grids_.at(macro_entity_index)); + } + + const LocalGridProviderType& local_grid(const size_t macro_entity_index) const + { + return *(local_grids_.at(macro_entity_index)); + } + + /** + * \brief Returns (and creates, if it does not exist) the coupling glue between the local grid view of level + * local_level_macro_entity on macro_entity and the local grid view of level local_level_macro_neighbor on + * macro_neighbor. + * \note Access is not implemented efficiently. This could be improved by using std::map::find instead of + * std::map::operator[]. + */ + const GlueType& coupling(const MacroEntityType& macro_entity, + const int local_level_macro_entity, + const MacroEntityType& macro_neighbor, + const int local_level_macro_neighbor, + const bool allow_for_broken_orientation_of_coupling_intersections = false) + { + const auto& macro_index_set = macro_leaf_view_.indexSet(); + assert(macro_index_set.contains(macro_entity)); + assert(macro_index_set.contains(macro_neighbor)); + if (local_level_macro_entity > max_local_level(macro_entity)) + DUNE_THROW(XT::Common::Exceptions::you_are_using_this_wrong, + "max_local_level(macro_entity): " << max_local_level(macro_entity) << "\n" + << " local_level_macro_entity: " + << local_level_macro_entity); + if (local_level_macro_neighbor > max_local_level(macro_neighbor)) + DUNE_THROW(XT::Common::Exceptions::you_are_using_this_wrong, + "max_local_level(macro_neighbor): " << max_local_level(macro_neighbor) << "\n" + << " local_level_macro_neighbor: " + << local_level_macro_neighbor); + const auto entity_index = macro_index_set.index(macro_entity); + const auto neighbor_index = macro_index_set.index(macro_neighbor); + if (glues_[entity_index][neighbor_index][local_level_macro_entity][local_level_macro_neighbor] == nullptr) { + // find the corresponding macro intersection ... + const auto macro_intersection_it_end = macro_leaf_view_.iend(macro_entity); + for (auto macro_intersection_it = macro_leaf_view_.ibegin(macro_entity); + macro_intersection_it != macro_intersection_it_end; + ++macro_intersection_it) { + const auto& macro_intersection = *macro_intersection_it; + if (macro_intersection.neighbor() && !macro_intersection.boundary()) { + const auto real_neighbor_ptr = macro_intersection.outside(); +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) + const auto& real_neighbor = real_neighbor_ptr; +#else + const auto& real_neighbor = *real_neighbor_ptr; +#endif + if (macro_index_set.index(real_neighbor) == neighbor_index) + glues_[entity_index][neighbor_index][local_level_macro_entity][local_level_macro_neighbor] = create_glue( + macro_entity, macro_neighbor, macro_intersection, local_level_macro_entity, local_level_macro_neighbor); + } + } // ... find the corresponding macro intersection + } + + const auto& glue = *(glues_[entity_index][neighbor_index][local_level_macro_entity][local_level_macro_neighbor]); + if (!allow_for_broken_orientation_of_coupling_intersections) { + const size_t brocken_intersections = check_for_broken_coupling_intersections(glue); + if (brocken_intersections > 0) + DUNE_THROW(Exceptions::intersection_orientation_is_broken, + "The coupling glue between the grid views of\n" + << " level " + << local_level_macro_entity + << " on macro entity " + << macro_leaf_view_.indexSet().index(macro_entity) + << " and\n" + << " level " + << local_level_macro_neighbor + << " on macro neighbor " + << macro_leaf_view_.indexSet().index(macro_neighbor) + << "\n" + << " contains\n" + << " " + << brocken_intersections + << "/" + << glue.size() + << " intersections with wrong orientation!"); + } + return glue; + } // ... coupling(...) + + const GlueType& coupling(const MacroEntityType& macro_entity, + const MacroEntityType& macro_neighbor, + const bool allow_for_broken_orientation_of_coupling_intersections = false) + { + return coupling(macro_entity, + max_local_level(macro_entity), + macro_neighbor, + max_local_level(macro_neighbor), + allow_for_broken_orientation_of_coupling_intersections); + } + + int max_local_level(const MacroEntityType& macro_entity) const + { + return local_grid(macro_entity).grid().maxLevel(); + } + + int max_local_level(const size_t macro_entity_index) const + { + return local_grid(macro_entity_index).grid().maxLevel(); + } + + const std::vector<std::pair<MicroEntityPointerType, std::vector<int>>>& + local_boundary_entities(const MacroEntityType& macro_entity, const int local_level) + { + // auto logger = Stuff::Common::TimedLogger().get("grid-multiscale.glued.local_boundary_entities"); + assert(macro_leaf_view_.indexSet().contains(macro_entity)); + const size_t macro_entity_index = macro_leaf_view_.indexSet().index(macro_entity); + if (local_level > max_local_level(macro_entity)) + DUNE_THROW(XT::Common::Exceptions::you_are_using_this_wrong, + "macro_entity_index: " << macro_entity_index << "\n" + << "local_level: " + << local_level + << "\n" + << "max_local_level(macro_entity): " + << max_local_level(macro_entity)); + auto& local_level_to_boundary_entity_ptrs_with_local_intersections = + macro_entity_to_local_level_to_boundary_entity_ptrs_with_local_intersections_[macro_entity_index]; + auto& boundary_entity_ptrs_with_local_intersections = + local_level_to_boundary_entity_ptrs_with_local_intersections[local_level]; + // logger.debug() << "macro_entity: " << macro_entity_index << std::endl; + if (boundary(macro_entity) && boundary_entity_ptrs_with_local_intersections.empty()) { + // create the container, therefore + const auto local_leaf_view = local_grids_[macro_entity_index]->leaf_view(); + // * walk the local grid (manually, to have access to the entity pointer) + const auto local_entity_it_end = local_leaf_view.template end<0>(); + for (auto local_entity_it = local_leaf_view.template begin<0>(); local_entity_it != local_entity_it_end; + ++local_entity_it) { + const auto& local_entity = *local_entity_it; + // logger.debug() << "local_entity: " << local_leaf_view.indexSet().index(local_entity) << " "; + if (local_entity.hasBoundaryIntersections()) { + // logger.debug() << "(boundary entity)" << std::endl; + std::vector<int> local_boundary_intersections; + // This entity has intersections on the local grid boundary, those could either be the domain boundary (which + // we are looking for) or a boundary to another local grid (which we are not looking for). To find out + // * walk the intersections + for (auto&& local_intersection : intersections(local_leaf_view, local_entity)) { + if (local_intersection.boundary() && !local_intersection.neighbor()) { + // logger.debug() << "local_intersection: " << local_intersection.indexInInside() << ":" << + // std::endl; + const auto local_intersection_geometry = local_intersection.geometry(); + const size_t num_corners = boost::numeric_cast<size_t>(local_intersection_geometry.corners()); + // ** Check if all corners of the intersection lie on the domain boundary (aka the boundary intersection + // of + // the macro entity this local grid belongs to. Therefore + // *** walk the intersections of the macro entity + for (auto&& macro_intersection : intersections(macro_leaf_view_, macro_entity)) { + if (macro_intersection.boundary() && !macro_intersection.neighbor()) { + // This macro intersection lies on the domain boundary, check if the local intersection is contained. + size_t corners_lie_on_boundary = 0; + for (size_t ii = 0; ii < num_corners; ++ii) { + if (XT::Grid::contains(macro_intersection, + local_intersection_geometry.corner(boost::numeric_cast<int>(ii)))) + ++corners_lie_on_boundary; + } + if (corners_lie_on_boundary == num_corners) { + // add the information to the container + local_boundary_intersections.push_back(local_intersection.indexInInside()); + } // add the information to the container + } + } // *** walk the intersections of the macro entity + } + } // * walk the intersections + if (!local_boundary_intersections.empty()) { + // add this local entity and its local intersections to the container + boundary_entity_ptrs_with_local_intersections.emplace_back(MicroEntityPointerType(*local_entity_it), + local_boundary_intersections); + } + } /*else + logger.debug() << "(inner entity)" << std::endl;*/ + } // * walk the local grid + } // create the container + return boundary_entity_ptrs_with_local_intersections; + } // ... local_boundary_entities(...) + + template <class... Args> + void visualize(const std::string& filename = "grid.multiscale.glued", Args&&... args) + { + auto logger = XT::Common::TimedLogger().get("grid-multiscale.glued.visualize"); + macro_grid_.visualize(filename + ".macro", std::forward<Args>(args)...); + GluedVTKWriter<MacroGridType, LocalGridType> vtk_writer(*this); + const auto& macro_index_set = macro_leaf_view_.indexSet(); + std::vector<std::vector<double>> subdomain_visualization(macro_index_set.size(0)); + std::vector<std::vector<double>> boundary_visualization(macro_index_set.size(0)); + std::vector<std::vector<double>> inside_outside_coupling_visualization(macro_index_set.size(0)); + std::vector<std::vector<double>> outside_inside_coupling_visualization(macro_index_set.size(0)); + // walk the macro grid + for (auto&& macro_entity : elements(macro_leaf_view_)) { + const size_t macro_entity_index = macro_index_set.index(macro_entity); + logger.debug() << "macro_entity: " << macro_entity_index << " "; + const auto local_level = max_local_level(macro_entity); + const auto local_grid_view = local_grids_[macro_entity_index]->level_view(local_level); + subdomain_visualization[macro_entity_index] = + std::vector<double>(local_grid_view.indexSet().size(0), macro_entity_index); + boundary_visualization[macro_entity_index] = std::vector<double>(local_grid_view.indexSet().size(0), -1); + if (inside_outside_coupling_visualization[macro_entity_index].empty()) + inside_outside_coupling_visualization[macro_entity_index] = + std::vector<double>(local_grid_view.indexSet().size(0), -1); + if (outside_inside_coupling_visualization[macro_entity_index].empty()) + outside_inside_coupling_visualization[macro_entity_index] = + std::vector<double>(local_grid_view.indexSet().size(0), -1); + // local boundary entities + if (boundary(macro_entity)) { + logger.debug() << "(boundary entity)" << std::endl; + const auto& boundary_entity_ptrs_with_local_intersections = local_boundary_entities(macro_entity, local_level); + logger.debug() << " " << boundary_entity_ptrs_with_local_intersections.size() << "/" + << local_grid_view.indexSet().size(0) << " boundary entities" << std::endl; + for (const auto& element : boundary_entity_ptrs_with_local_intersections) { + const auto& local_entity_ptr = element.first; + const auto& local_intersections = element.second; + if (!local_intersections.empty()) { + const auto& local_entity = *local_entity_ptr; + const size_t local_entity_index = local_grid_view.indexSet().index(local_entity); + boundary_visualization[macro_entity_index][local_entity_index] = macro_entity_index; + } + } + } else + logger.debug() << "(inner entity)" << std::endl; + for (auto&& macro_intersection : intersections(macro_leaf_view_, macro_entity)) { + if (!macro_intersection.boundary() && macro_intersection.neighbor()) { + const auto macro_neighbor_ptr = macro_intersection.outside(); +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) + const auto& macro_neighbor = macro_neighbor_ptr; +#else + const auto& macro_neighbor = *macro_neighbor_ptr; +#endif + const size_t macro_neighbor_index = macro_leaf_view_.indexSet().index(macro_neighbor); + const auto local_neighbor_level = max_local_level(macro_neighbor); + const auto local_neighbor_grid_view = local_grids_[macro_neighbor_index]->level_view(local_neighbor_level); + if (inside_outside_coupling_visualization[macro_neighbor_index].empty()) + inside_outside_coupling_visualization[macro_neighbor_index] = + std::vector<double>(local_neighbor_grid_view.indexSet().size(0), -1); + if (outside_inside_coupling_visualization[macro_neighbor_index].empty()) + outside_inside_coupling_visualization[macro_neighbor_index] = + std::vector<double>(local_neighbor_grid_view.indexSet().size(0), -1); + // walk the coupling, where this is the inside + size_t num_coupling_intersections = 0; + const auto& in_out_coupling_glue = coupling(macro_entity, + local_level, + macro_neighbor, + local_neighbor_level, + /*allow_for_broken_orientation_of_coupling_intersections=*/true); + const auto in_out_coupling_intersection_it_end = in_out_coupling_glue.template iend<0>(); + for (auto in_out_coupling_intersection_it = in_out_coupling_glue.template ibegin<0>(); + in_out_coupling_intersection_it != in_out_coupling_intersection_it_end; + ++in_out_coupling_intersection_it) { + ++num_coupling_intersections; + const auto& coupling_intersection = *in_out_coupling_intersection_it; + const auto local_entity = coupling_intersection.inside(); + const size_t local_entity_index = local_grid_view.indexSet().index(local_entity); + inside_outside_coupling_visualization[macro_entity_index][local_entity_index] = macro_entity_index; + const auto local_neighbor = coupling_intersection.outside(); + const size_t local_neighbor_index = local_neighbor_grid_view.indexSet().index(local_neighbor); + inside_outside_coupling_visualization[macro_neighbor_index][local_neighbor_index] = macro_neighbor_index; + } + // walk the coupling, where this is the outside + size_t out_in_num_coupling_intersections = 0; + const auto& out_in_coupling_glue = coupling(macro_neighbor, + local_neighbor_level, + macro_entity, + local_level, + /*allow_for_broken_orientation_of_coupling_intersections=*/true); + const auto out_in_coupling_intersection_it_end = out_in_coupling_glue.template iend<0>(); + for (auto out_in_coupling_intersection_it = out_in_coupling_glue.template ibegin<0>(); + out_in_coupling_intersection_it != out_in_coupling_intersection_it_end; + ++out_in_coupling_intersection_it) { + ++out_in_num_coupling_intersections; + const auto& coupling_intersection = *out_in_coupling_intersection_it; + const auto local_entity = coupling_intersection.inside(); + const size_t local_entity_index = local_grid_view.indexSet().index(local_entity); + outside_inside_coupling_visualization[macro_neighbor_index][local_entity_index] = macro_neighbor_index; + const auto local_neighbor = coupling_intersection.outside(); + const size_t local_neighbor_index = local_neighbor_grid_view.indexSet().index(local_neighbor); + outside_inside_coupling_visualization[macro_entity_index][local_neighbor_index] = macro_entity_index; + } + if (num_coupling_intersections != out_in_num_coupling_intersections) + DUNE_THROW(XT::Common::Exceptions::internal_error, + "The coupling glue is broken!\n" + << "macro entity (local level): " + << macro_entity_index + << " (" + << local_level + << ")\n" + << "macro neighbor (local level): " + << macro_neighbor_index + << " (" + << local_neighbor_level + << ")"); + logger.debug() << " " << num_coupling_intersections << " coupling intersections with neighbor " + << macro_neighbor_index << std::endl; + } + } + } // walk the macro grid + vtk_writer.addCellData(subdomain_visualization, "subdomains"); + vtk_writer.addCellData(boundary_visualization, "local boundary entities"); + vtk_writer.addCellData(inside_outside_coupling_visualization, "local coupling entities (inside/outside)"); + vtk_writer.addCellData(outside_inside_coupling_visualization, "local coupling entities (outside/inside)"); + vtk_writer.write(filename, VTK::appendedraw); + } // ... visualize(...) + +private: + template <class MacroEntityType> + static std::shared_ptr<LocalGridProviderType> create_grid_of_simplex(const MacroEntityType& macro_entity) + { + try { + GridFactory<LocalGridType> subdomain_factory; + const auto num_vertices = macro_entity.subEntities(dimDomain); + std::vector<unsigned int> vertex_ids(num_vertices, 0); + for (auto&& local_vertex_id : XT::Common::value_range(num_vertices)) { + const auto vertex = macro_entity.template subEntity<dimDomain>(local_vertex_id).geometry().center(); + subdomain_factory.insertVertex(vertex); + vertex_ids[local_vertex_id] = local_vertex_id; + } + subdomain_factory.insertElement(macro_entity.geometry().type(), vertex_ids); + return std::make_shared<XT::Grid::GridProvider<LocalGridType>>(subdomain_factory.createGrid()); + } catch (GridError& ee) { + DUNE_THROW(GridError, + "It was not possible to create a grid for this simplex with the given GridType!\n\n" + << "GridType: " + << XT::Common::Typename<LocalGridType>::value() + << "\n\n" + << "This was the original error: " + << ee.what()); + } + } // ... create_grid_of_simplex(...) + + template <class MacroEntityType> + static std::shared_ptr<LocalGridProviderType> create_grid_of_cube(const MacroEntityType& macro_entity) + { + const auto num_vertices = macro_entity.subEntities(dimDomain); + FieldVector<ctype, dimDomain> lower_left(std::numeric_limits<ctype>::max()); + FieldVector<ctype, dimDomain> upper_right(std::numeric_limits<ctype>::min()); + for (auto&& local_vertex_id : XT::Common::value_range(num_vertices)) { + const auto vertex = macro_entity.template subEntity<dimDomain>(local_vertex_id).geometry().center(); + for (size_t dd = 0; dd < dimDomain; ++dd) { + lower_left[dd] = std::min(lower_left[dd], vertex[dd]); + upper_right[dd] = std::max(upper_right[dd], vertex[dd]); + } + } + std::array<unsigned int, dimDomain> num_refinements; + std::fill(num_refinements.begin(), num_refinements.end(), 1); + return std::make_shared<XT::Grid::GridProvider<LocalGridType>>( + XT::Grid::make_cube_grid<LocalGridType>(lower_left, upper_right, num_refinements)); + } // ... create_grid_of_cube(...) + + void setup_local_grids() + { + const auto& macro_index_set = macro_leaf_view_.indexSet(); + for (auto&& macro_entity : elements(macro_leaf_view_)) { + auto macro_entity_index = macro_index_set.index(macro_entity); + if (macro_entity.type().isSimplex()) + local_grids_[macro_entity_index] = create_grid_of_simplex(macro_entity); + else if (macro_entity.type().isCube()) + local_grids_[macro_entity_index] = create_grid_of_cube(macro_entity); + else + DUNE_THROW(GridError, "Unknown entity.type() encountered: " << macro_entity.type()); + } + } // ... setup_local_grids() + + template <class MacroIntersectionType> + std::shared_ptr<GlueType> create_glue(const MacroEntityType& macro_entity, + const MacroEntityType& macro_neighbor, + const MacroIntersectionType& macro_intersection, + const int local_entity_level, + const int local_neighbor_level) const + { + assert(local_entity_level >= 0); + assert(local_neighbor_level >= 0); + const auto& local_entity_grid = local_grid(macro_entity); + const auto& local_neighbor_grid = local_grid(macro_neighbor); + assert(local_entity_level <= local_entity_grid.grid().maxLevel()); + assert(local_neighbor_level <= local_neighbor_grid.grid().maxLevel()); + auto local_entity_view = local_entity_grid.level_view(local_entity_level); + auto local_neighbor_view = local_neighbor_grid.level_view(local_neighbor_level); + // create descriptors, these can be discarded after creating the extractors + auto entity_descriptor = create_descriptor(local_entity_view, macro_intersection); + auto neighbor_descriptor = create_descriptor(local_neighbor_view, macro_intersection); + // create extractors and merger as shared_ptr, so glue will handle memory + auto entity_extractor = std::make_shared<LocalExtractorType>(local_entity_view, entity_descriptor); + auto neighbor_extractor = std::make_shared<LocalExtractorType>(local_neighbor_view, neighbor_descriptor); + auto contact_merger = std::make_shared<GridGlue::ContactMerge<dimWorld, ctype>>(allowed_overlap_); + // create glue + auto glue = std::make_shared<GlueType>(entity_extractor, neighbor_extractor, contact_merger); + glue->build(); + if (glue->size() == 0) + DUNE_THROW(GridError, + "Something went wrong, the coupling glue is empty!\n" + << " macro_entity " + << macro_leaf_view_.indexSet().index(macro_entity) + << "\n" + << " local_entity_level " + << local_entity_level + << "\n" + << " macro_neighbor " + << macro_leaf_view_.indexSet().index(macro_neighbor) + << "\n" + << " local_neighbor_level " + << local_neighbor_level + << "\n"); + return glue; + } // ... create_glue(...) + + void setup_glues(const bool allow_for_broken_orientation_of_coupling_intersections = false) + { + const auto& macro_index_set = macro_leaf_view_.indexSet(); + for (auto&& macro_entity : elements(macro_leaf_view_)) { + const auto macro_entity_index = macro_index_set.index(macro_entity); + auto& entity_glues = glues_[macro_entity_index]; + // walk the neighbors ... + const auto macro_intersection_it_end = macro_leaf_view_.iend(macro_entity); + for (auto macro_intersection_it = macro_leaf_view_.ibegin(macro_entity); + macro_intersection_it != macro_intersection_it_end; + ++macro_intersection_it) { + const auto& macro_intersection = *macro_intersection_it; + if (macro_intersection.neighbor() && !macro_intersection.boundary()) { + const auto macro_neighbor_ptr = macro_intersection.outside(); +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) + const auto& macro_neighbor = macro_neighbor_ptr; +#else + const auto& macro_neighbor = *macro_neighbor_ptr; +#endif + const auto macro_neighbor_index = macro_index_set.index(macro_neighbor); + for (auto local_entity_level : + XT::Common::value_range(local_grids_[macro_entity_index]->grid().maxLevel() + 1)) + for (auto local_neighbor_level : + XT::Common::value_range(local_grids_[macro_neighbor_index]->grid().maxLevel() + 1)) { + auto glue = create_glue( + macro_entity, macro_neighbor, macro_intersection, local_entity_level, local_neighbor_level); + if (!allow_for_broken_orientation_of_coupling_intersections) { + const size_t brocken_intersections = check_for_broken_coupling_intersections(*glue); + if (brocken_intersections > 0) + DUNE_THROW(Exceptions::intersection_orientation_is_broken, + "The coupling glue between the grid views of\n" + << " level " + << local_entity_level + << " on macro entity " + << macro_leaf_view_.indexSet().index(macro_entity) + << " and\n" + << " level " + << local_neighbor_level + << " on macro neighbor " + << macro_leaf_view_.indexSet().index(macro_neighbor) + << "\n" + << "contains\n" + << " " + << brocken_intersections + << "/" + << glue->size() + << " intersections with wrong" + << "orientation!"); + } + entity_glues[macro_neighbor_index][local_entity_level][local_neighbor_level] = glue; + } + } + } // ... walk the neighbors + } + } // ... setup_glues(...) + + void prepare_global_grid() + { + if (global_grid_) + return; + const auto& macro_index_set = macro_leaf_view_.indexSet(); + std::vector<FieldVector<ctype, dimDomain>> vertices; + std::vector<std::vector<std::vector<unsigned int>>> entity_to_vertex_ids(local_grids_.size()); + std::vector<std::vector<GeometryType>> geometry_types(local_grids_.size()); + // walk the grid for the first time + for (auto&& macro_entity : elements(macro_leaf_view_)) { + const auto local_leaf_view = local_grid(macro_entity).leaf_view(); + const auto& local_index_set = local_leaf_view.indexSet(); + const size_t macro_index = macro_index_set.index(macro_entity); + entity_to_vertex_ids[macro_index] = std::vector<std::vector<unsigned int>>(local_index_set.size(0)); + geometry_types[macro_index] = std::vector<GeometryType>(local_index_set.size(0)); + for (auto&& micro_entity : elements(local_leaf_view)) { + const size_t micro_index = local_index_set.index(micro_entity); + const auto num_vertices = micro_entity. +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) + + subEntities(dimDomain); +#else + template count<dimDomain>(); +#endif + entity_to_vertex_ids[macro_index][micro_index] = std::vector<unsigned int>(num_vertices); + geometry_types[macro_index][micro_index] = micro_entity.geometry().type(); + for (unsigned int local_vertex_id = 0; local_vertex_id < num_vertices; ++local_vertex_id) { + const unsigned int global_vertex_id = find_insert_vertex(vertices, + micro_entity + .template subEntity<dimDomain>(local_vertex_id) +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) + . +#else + -> +#endif + geometry() + .center()); + entity_to_vertex_ids[macro_index][micro_index][local_vertex_id] = global_vertex_id; + } + } // walk the local grid + } // walk the macro grid + GridFactory<LocalGridType> global_factory; + for (const auto& vertex : vertices) + global_factory.insertVertex(vertex); + size_t II = 0; + size_t JJ = 0; + try { + for (size_t ii = 0; ii < local_grids_.size(); ++ii, II = ii) + for (size_t jj = 0; jj < entity_to_vertex_ids[ii].size(); ++jj, JJ = jj) + global_factory.insertElement(geometry_types[ii][jj], entity_to_vertex_ids[ii][jj]); + } catch (GridError& ee) { + DUNE_THROW(GridError, + "It was not possible to insert an element into the grid factory!\n\n" + << "GridType: " + << XT::Common::Typename<LocalGridType>::value() + << "\n" + << "GeometryType: " + << geometry_types[II][JJ] + << "\n" + << "\n" + << "This was the original error: " + << ee.what()); + } // try + global_grid_ = std::make_unique<XT::Grid::GridProvider<LocalGridType>>(global_factory.createGrid()); + + // build maps of indices, relating the local to the global grid + const auto global_view = global_grid_->leaf_view(); + local_to_global_indices_ = std::make_unique<std::vector<std::vector<size_t>>>(macro_leaf_view_.indexSet().size(0)); + auto& local_to_global_inds = *local_to_global_indices_; + global_to_local_indices_ = std::make_unique<std::vector<std::pair<size_t, size_t>>>(global_view.indexSet().size(0)); + auto& global_to_local_inds = *global_to_local_indices_; + // therefore + // * create a search on the global view: if all is fine, we only need to walk that one once + std::vector<FieldVector<ctype, dimDomain>> local_entity_center{FieldVector<ctype, dimDomain>(0.0)}; + auto global_search = XT::Grid::make_entity_in_level_search(global_view); + // * walk the macro grid + for (auto&& macro_entity : elements(macro_leaf_view_)) { + const size_t subd = macro_leaf_view_.indexSet().index(macro_entity); + const auto local_leaf_view = local_grid(macro_entity).leaf_view(); + const auto& local_index_set = local_leaf_view.indexSet(); + local_to_global_inds[subd] = std::vector<size_t>(local_index_set.size(0)); + // * walk the local grid + for (auto&& local_entity : elements(local_leaf_view)) { + const size_t local_entity_index = local_index_set.index(local_entity); + local_entity_center[0] = local_entity.geometry().center(); + const auto global_entity_ptr_unique_ptrs = global_search(local_entity_center); + // the search has to be successfull, since the global grid has been constructed to exactly contain each local + // entity + assert(global_entity_ptr_unique_ptrs.size() == 1); + const auto& global_entity_ptr_unique_ptr = global_entity_ptr_unique_ptrs.at(0); + assert(global_entity_ptr_unique_ptr); + const auto& global_entity_ptr = *global_entity_ptr_unique_ptr; + const auto& global_entity = *global_entity_ptr; + const size_t global_entity_index = global_view.indexSet().index(global_entity); + // store information + local_to_global_inds[subd][local_entity_index] = global_entity_index; + global_to_local_inds[global_entity_index] = {subd, local_entity_index}; + } // * walk the local grid + } // * walk the macro grid + } // ... prepare_global_grid(...) + + size_t find_insert_vertex(std::vector<FieldVector<ctype, dimDomain>>& vertices, + FieldVector<ctype, dimDomain>&& vertex) const + { + // check if vertex is already contained + for (size_t ii = 0; ii < vertices.size(); ++ii) + if (XT::Common::FloatCmp::eq(vertex, vertices[ii])) + return ii; + // if not, add it + vertices.emplace_back(std::move(vertex)); + return vertices.size() - 1; + } // ... find_insert_vertex(...) + + MacroGridProviderType& macro_grid_; + const ctype allowed_overlap_; + MacroGridViewType macro_leaf_view_; + std::vector<std::shared_ptr<LocalGridProviderType>> local_grids_; + std::vector<std::map<size_t, std::map<int, std::map<int, std::shared_ptr<GlueType>>>>> glues_; + std::map<size_t, std::map<int, std::vector<std::pair<MicroEntityPointerType, std::vector<int>>>>> + macro_entity_to_local_level_to_boundary_entity_ptrs_with_local_intersections_; + std::unique_ptr<LocalGridProviderType> global_grid_; + std::unique_ptr<std::vector<std::vector<size_t>>> local_to_global_indices_; + std::unique_ptr<std::vector<std::pair<size_t, size_t>>> global_to_local_indices_; +}; // class Glued + + +template <class MacroGridType, class LocalGridType> +class GluedVTKWriter +{ + typedef typename LocalGridType::LevelGridView LocalGridViewType; + + // we only need this class to access a protected pwrite method of VTKWriter + class LocalVTKWriter : public VTKWriter<LocalGridViewType> + { + typedef VTKWriter<LocalGridViewType> BaseType; + + public: + LocalVTKWriter(const LocalGridViewType& local_grid_view, const size_t subdomain, const size_t num_subdomains) + : BaseType(local_grid_view) + , commRank_(boost::numeric_cast<int>(subdomain)) + , commSize_(boost::numeric_cast<int>(num_subdomains)) + { + } + + void write_locally(const std::string& name, VTK::OutputType ot) + { + BaseType::pwrite(name, "", "", ot, commRank_, commSize_); + } + + private: + const int commRank_; + const int commSize_; + }; // class LocalVTKWriter + +public: + typedef Glued<MacroGridType, LocalGridType> GluedGridType; + + GluedVTKWriter(const GluedGridType& glued_grid, const int local_level = -1) + : glued_grid_(glued_grid) + , local_levels_(glued_grid_.num_subdomains(), -1) + { + // set each local level to its respective max + if (local_level < 0) + for (auto&& macro_entity : elements(glued_grid_.macro_grid_view())) + local_levels_[glued_grid_.subdomain(macro_entity)] = glued_grid_.max_local_level(macro_entity); + prepare_local_vtk_writers(); + } // GluedVTKWriter(...) + + GluedVTKWriter(const GluedGridType& glued_grid, const std::vector<int>& local_levels) + : glued_grid_(glued_grid) + , local_levels_(local_levels) + { + for (size_t ss = 0; ss < glued_grid_.num_subdomains(); ++ss) + if (local_levels_[ss] < 0) + local_levels_[ss] = glued_grid_.max_local_level(ss); + prepare_local_vtk_writers(); + } + + template <class V> + void addCellData(const std::vector<std::vector<V>>& vectors, const std::string& name, const int ncomps = 1) + { + if (vectors.size() != glued_grid_.num_subdomains()) + DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, + "vectors.size(): " << vectors.size() << "\n" + << "glued_grid_.num_subdomains(): " + << glued_grid_.num_subdomains()); + for (size_t ss = 0; ss < glued_grid_.num_subdomains(); ++ss) + local_vtk_writers_[ss]->addCellData(vectors[ss], name, ncomps); + } // ... addCellData(...) + + template <class VTKFunctionType> + void addCellData(const std::vector<std::shared_ptr<VTKFunctionType>>& functions) + { + if (functions.size() != glued_grid_.num_subdomains()) + DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, + "funcitons.size(): " << functions.size() << "\n" + << "glued_grid_.num_subdomains(): " + << glued_grid_.num_subdomains()); + for (size_t ss = 0; ss < glued_grid_.num_subdomains(); ++ss) + local_vtk_writers_[ss]->addCellData(functions[ss]); + } // ... addCellData(...) + + template <class V> + void addVertexData(const std::vector<std::vector<V>>& vectors, const std::string& name, const int ncomps = 1) + { + if (vectors.size() != glued_grid_.num_subdomains()) + DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, + "vectors.size(): " << vectors.size() << "\n" + << "glued_grid_.num_subdomains(): " + << glued_grid_.num_subdomains()); + for (size_t ss = 0; ss < glued_grid_.num_subdomains(); ++ss) + local_vtk_writers_[ss]->addVertexData(vectors[ss], name, ncomps); + } // ... addVertexData(...) + + template <class VTKFunctionType> + void addVertexData(const std::vector<std::shared_ptr<VTKFunctionType>>& functions) + { + if (functions.size() != glued_grid_.num_subdomains()) + DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, + "funcitons.size(): " << functions.size() << "\n" + << "glued_grid_.num_subdomains(): " + << glued_grid_.num_subdomains()); + for (size_t ss = 0; ss < glued_grid_.num_subdomains(); ++ss) + local_vtk_writers_[ss]->addVertexData(functions[ss]); + } // ... addVertexData(...) + + void clear() + { + for (auto& local_vtk_writer : local_vtk_writers_) + local_vtk_writer.clear(); + } + + void write(const std::string& name, VTK::OutputType type = VTK::ascii) + { + for (size_t ss = 0; ss < glued_grid_.num_subdomains(); ++ss) + local_vtk_writers_[ss]->write_locally(name, type); + } + +private: + void prepare_local_vtk_writers() + { + for (auto&& macro_entity : elements(glued_grid_.macro_grid_view())) { + const size_t subdomain = glued_grid_.subdomain(macro_entity); + local_vtk_writers_.emplace_back( + new LocalVTKWriter(glued_grid_.local_grid(macro_entity).level_view(local_levels_[subdomain]), + subdomain, + glued_grid_.num_subdomains())); + } + } // ... prepare_local_vtk_writers(...) + + const GluedGridType& glued_grid_; + std::vector<int> local_levels_; + std::vector<std::unique_ptr<LocalVTKWriter>> local_vtk_writers_; +}; // class GluedVTKWriter + + +#else // HAVE_DUNE_GRID_GLUE + + +template <class MacroGridType, class LocalGridType> +class GluedVTKWriter +{ + static_assert(AlwaysFalse<MacroGridImp>::value, "You are missing dune-grid-glue!"); +}; + + +template <class MacroGridImp, class LocalGridImp> +class Glued +{ + static_assert(AlwaysFalse<MacroGridImp>::value, "You are missing dune-grid-glue!"); +}; + + +#endif // HAVE_DUNE_GRID_GLUE + + +} // namespace DD +} // namespace Grid +} // namespace XT +} // namespace Dune + +#endif // DUNE_XT_GRID_DD_GLUED_HH diff --git a/dune/xt/grid/gridprovider.pbh b/dune/xt/grid/gridprovider.pbh new file mode 100644 index 0000000000000000000000000000000000000000..bec78aa49c73a312d748b1a75cf663b62a0eef6b --- /dev/null +++ b/dune/xt/grid/gridprovider.pbh @@ -0,0 +1,86 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2016 - 2017) + +#ifndef DUNE_XT_GRID_GRIDPROVIDER_PBH +#define DUNE_XT_GRID_GRIDPROVIDER_PBH +#if HAVE_DUNE_PYBINDXI + +#include <dune/pybindxi/pybind11.h> +#include <dune/pybindxi/stl.h> + +#include <dune/xt/common/configuration.pbh> +#include <dune/xt/common/fvector.pbh> + +#include "gridprovider.hh" + +namespace Dune { +namespace XT { +namespace Grid { + + +template <class G> +pybind11::class_<GridProvider<G>> bind_GridProvider(pybind11::module& m, const std::string& grid_id) +{ + namespace py = pybind11; + using namespace pybind11::literals; + + typedef GridProvider<G> C; + + py::class_<C> c(m, std::string("GridProvider__" + grid_id).c_str(), std::string("GridProvider__" + grid_id).c_str()); + + c.def("max_level", &C::max_level); + c.def("global_refine", &C::global_refine, "count"_a); + c.def("visualize", + [](const C& self, const std::string& filename, const Common::Configuration& boundary_info_cfg) { + self.visualize(filename, boundary_info_cfg); + }, + "filename"_a = C::static_id(), + "boundary_info_cfg"_a = Common::Configuration()); + + c.def_property_readonly("grid_type", [grid_id](const C& /*self*/) { return grid_id; }); + c.def_property_readonly("dim", [](const C& /*self*/) { return C::dimDomain; }); + + return c; +} // ... bind_GridProvider(...) + + +template <class G> +void bind_make_cube_grid(pybind11::module& m, const std::string& grid_id) +{ + namespace py = pybind11; + using namespace pybind11::literals; + + m.def(std::string("make_cube_grid__" + grid_id).c_str(), + [](const Common::Configuration& cfg) { return make_cube_grid<G>(cfg); }, + "cfg"_a = cube_gridprovider_default_config()); + + m.def(std::string("make_cube_grid__" + grid_id).c_str(), + [](const FieldVector<typename G::ctype, G::dimension>& lower_left, + const FieldVector<typename G::ctype, G::dimension>& upper_right, + const std::array<unsigned int, G::dimension>& num_elements, + const unsigned int num_refinements, + const std::array<unsigned int, G::dimension>& overlap_size) { + return make_cube_grid<G>(lower_left, upper_right, num_elements, num_refinements, overlap_size); + }, + "lower_left"_a, + "upper_right"_a, + "num_elements"_a = + cube_gridprovider_default_config().template get<std::vector<unsigned int>>("num_elements").at(0), + "num_refinements"_a = cube_gridprovider_default_config().template get<unsigned int>("num_refinements"), + "overlap_size"_a = + cube_gridprovider_default_config().template get<std::vector<unsigned int>>("overlap_size").at(0)); +} // ... bind_make_cube_grid(...) + + +} // namespace Grid +} // namespace XT +} // namespace Dune + +#endif // HAVE_DUNE_PYBINDXI +#endif // DUNE_XT_GRID_GRIDPROVIDER_PBH diff --git a/dune/xt/grid/grids.bindings.hh b/dune/xt/grid/grids.bindings.hh new file mode 100644 index 0000000000000000000000000000000000000000..102ed1a65a9b5c92865b666c2ff871e804297b9a --- /dev/null +++ b/dune/xt/grid/grids.bindings.hh @@ -0,0 +1,18 @@ +#ifndef DUNE_XT_GRID_GRIDS_BINDINGS_HH +#define DUNE_XT_GRID_GRIDS_BINDINGS_HH + +#include "grids.hh" + +namespace Dune { + + +// this is used by other headers +typedef YaspGrid<2, EquidistantOffsetCoordinates<double, 2>> YASP_2D_EQUIDISTANT_OFFSET; +#if HAVE_ALUGRID || HAVE_DUNE_ALUGRID +typedef ALUGrid<2, 2, simplex, conforming> ALU_2D_SIMPLEX_CONFORMING; +#endif + + +} // namespace Dune + +#endif // DUNE_XT_GRID_GRIDS_BINDINGS_HH diff --git a/dune/xt/grid/grids.hh b/dune/xt/grid/grids.hh index 551a82487dfcd66710a71c713749110a24048183..101db0a3d82520fb62045b80f903fb425625b736 100644 --- a/dune/xt/grid/grids.hh +++ b/dune/xt/grid/grids.hh @@ -10,12 +10,11 @@ #ifndef DUNE_XT_GRID_GRIDS_HH #define DUNE_XT_GRID_GRIDS_HH -#if HAVE_ALBERTA // clang-format off +#if HAVE_ALBERTA # include <dune/xt/common/disable_warnings.hh> # include <dune/grid/albertagrid.hh> # include <dune/xt/common/reenable_warnings.hh> #endif -#include <dune/grid/yaspgrid.hh> #if HAVE_DUNE_ALUGRID # include <dune/alugrid/grid.hh> @@ -25,8 +24,10 @@ # include <dune/grid/spgrid.hh> #endif -#if HAVE_DUNE_UGGRID +#if HAVE_DUNE_UGGRID || HAVE_UG # include <dune/grid/uggrid.hh> -#endif // clang-format on +#endif + +#include <dune/grid/yaspgrid.hh> #endif // DUNE_XT_GRID_GRIDS_HH diff --git a/dune/xt/grid/intersection.hh b/dune/xt/grid/intersection.hh index 68cad31993fee0ff483965c4bf474249c248af25..a0d8b775a0d239c9342fdc22da0742fe5672fc0c 100644 --- a/dune/xt/grid/intersection.hh +++ b/dune/xt/grid/intersection.hh @@ -81,7 +81,7 @@ void print_intersection(const IntersectionType& intersection, /** - * \brief Checks if intersection contains the given global_point. + * \brief Checks if intersection contains the given global_point (2d variant). * * Returns true, if global_point lies on the line between the corners of intersection. */ @@ -116,6 +116,44 @@ contains(const Dune::Intersection<G, I>& intersection, } // ... contains(...) +/** + * \brief Checks if intersection contains the given global_point (3d variant). + * + * Checks if global_point lies within the plane spanned by the first three corners of intersection. + * + * \sa + * http://math.stackexchange.com/questions/684141/check-if-a-point-is-on-a-plane-minimize-the-use-of-multiplications-and-divisio + */ +template <class G, class I, class D> +typename std::enable_if<Dune::Intersection<G, I>::dimension == 3, bool>::type +contains(const Dune::Intersection<G, I>& intersection, + const Dune::FieldVector<D, 3>& global_point, + const D& tolerance = Common::FloatCmp::DefaultEpsilon<D>::value()) +{ + const auto& geometry = intersection.geometry(); + // get the global coordinates of the intersections corners, there should be at least 3 (ignore the fourth if there is + // one, 3 points is enough in 3d) + assert(geometry.corners() >= 3); + std::vector<Dune::FieldVector<D, 3>> points(4, Dune::FieldVector<D, 3>(0.)); + for (size_t ii = 0; ii < 3; ++ii) + points[ii] = geometry.corner(ii); + points[3] = global_point; + // form a matrix of these points, where the top three entries of each column are given by the three entries of each + // point and the bottom row is given by one, i.e. + // a_0 b_0 c_0 d_0 + // a_1 b_1 c_1 d_1 + // a_2 b_2 c_2 d_2 + // 1 1 1 1 + FieldMatrix<D, 4, 4> matrix(1.); // ensures the 1 on the last row + for (size_t ii = 0; ii < 3; ++ii) // only set the first three rows + for (size_t jj = 0; jj < 4; ++jj) + matrix[ii][jj] = points[jj][ii]; + // the point lies on the plane given by the corners if the determinant of this matrix is zero + const D det = matrix.determinant(); + return std::abs(det) < tolerance; +} // ... contains(...) + + } // namespace Grid } // namespace XT } // namespace Dune diff --git a/dune/xt/grid/test/dd_glued.hh b/dune/xt/grid/test/dd_glued.hh new file mode 100644 index 0000000000000000000000000000000000000000..821eebd82c8981b5cc0c826d9a6abeb761c45f1a --- /dev/null +++ b/dune/xt/grid/test/dd_glued.hh @@ -0,0 +1,302 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2017) + +#ifndef DUNE_XT_GRID_TEST_DD_GLUED_HH +#define DUNE_XT_GRID_TEST_DD_GLUED_HH +#if HAVE_DUNE_GRID_GLUE + +#include <map> +#include <set> +#include <sstream> +#include <vector> + +#include <dune/grid/yaspgrid.hh> +#include <dune/grid/common/rangegenerators.hh> + +#if HAVE_DUNE_ALUGRID +#include <dune/alugrid/grid.hh> +#elif HAVE_ALUGRID +#include <dune/grid/alugrid.hh> +#endif + +#include <dune/xt/common/memory.hh> +#include <dune/xt/common/test/gtest/gtest.h> +#include <dune/xt/grid/gridprovider.hh> +#include <dune/xt/grid/dd/glued.hh> + + +template <class T> +std::string convert_to_initializer_list_str(const std::set<T>& results) +{ + std::stringstream out; + if (results.size() == 0) + out << "{}"; + else if (results.size() == 1) + out << "{" << *results.begin() << "}"; + else { + auto iterator = results.begin(); + out << "{" << *iterator; + ++iterator; + for (; iterator != results.end(); ++iterator) { + out << ", " << *iterator; + } + out << "}"; + } + return out.str(); +} + +template <class L, class R> +std::string convert_to_initializer_list_str(const std::pair<L, R>& results) +{ + std::stringstream out; + out << "{" << results.first << ", " << results.second << "}"; + return out.str(); +} + +template <class F, class S> +std::string convert_to_initializer_list_str(const std::map<F, S>& results) +{ + std::stringstream out; + if (results.size() == 0) + out << "{}" << std::endl; + else if (results.size() == 1) + out << "{{" << convert_to_initializer_list_str(results.begin()->first) << ", " << results.begin()->second << "}}"; + else { + auto iterator = results.begin(); + out << "{{" << convert_to_initializer_list_str(iterator->first) << ", " << iterator->second << "}"; + ++iterator; + for (; iterator != results.end(); ++iterator) { + out << ",\n {" << convert_to_initializer_list_str(iterator->first) << ", " << iterator->second << "}"; + } + out << "}"; + } + return out.str(); +} + + +namespace Dune { +namespace XT { +namespace Grid { + + +template <class M, class L, bool aything = true> +struct ExpectedResults +{ + static_assert(AlwaysFalse<M>::value, "Please add me for this grid!"); +}; + + +/// \note assumes that all macro entities contain local grids of same refiment levels +template <class GridTuple> +struct GluedMultiscaleGridTest : public ::testing::Test +{ + typedef typename std::tuple_element<0, GridTuple>::type MacroGridType; + typedef typename std::tuple_element<1, GridTuple>::type LocalGridType; + typedef ExpectedResults<MacroGridType, LocalGridType> Expectations; + + void setup() + { + if (!macro_grid_) + macro_grid_ = std::make_unique<GridProvider<MacroGridType>>( + make_cube_grid<MacroGridType>(0., 1., 3, Expectations::num_coarse_refinements())); + ASSERT_NE(macro_grid_, nullptr) << "This should not happen!"; + if (!multiscale_grid_) + multiscale_grid_ = std::make_unique<DD::Glued<MacroGridType, LocalGridType>>( + *macro_grid_, + Expectations::num_local_refinements(), + /*prepare_glues=*/false, + /*allow_for_broken_orientation_of_coupling_intersections=*/true); + ASSERT_NE(multiscale_grid_, nullptr) << "This should not happen!"; + for (auto&& macro_entity : Dune::elements(multiscale_grid_->macro_grid_view())) { + EXPECT_EQ(multiscale_grid_->max_local_level(macro_entity), Expectations::num_local_refinements()); + } + } // ... setup() + + void couplings_are_of_correct_size() + { + setup(); + ASSERT_NE(macro_grid_, nullptr) << "This should not happen!"; + ASSERT_NE(multiscale_grid_, nullptr) << "This should not happen!"; + + const auto& macro_grid_view = multiscale_grid_->macro_grid_view(); + for (auto&& macro_entity : Dune::elements(macro_grid_view)) { + const auto entity_index = macro_grid_view.indexSet().index(macro_entity); + for (auto&& macro_intersection : Dune::intersections(macro_grid_view, macro_entity)) { + if (macro_intersection.neighbor() && !macro_intersection.boundary()) { + const auto macro_neighbor = macro_intersection.outside(); + const auto neighbor_index = macro_grid_view.indexSet().index(macro_neighbor); + const auto& coupling = + multiscale_grid_->coupling(macro_entity, + 1, + macro_neighbor, + Expectations::num_local_refinements(), + /*allow_for_broken_orientation_of_coupling_intersections=*/true); + EXPECT_EQ(Expectations::num_local_couplings_intersections().count(coupling.size()), 1) + << "entity: " << entity_index << "\n" + << "neighbor: " << neighbor_index << "\n" + << "expected num_local_couplings_intersections: " + << convert_to_initializer_list_str(Expectations::num_local_couplings_intersections()) + << "\nactual num_local_couplings_intersections: " << coupling.size(); + } + } + } + } // ... couplings_are_of_correct_size(...) + + void visualize_is_callable() + { + setup(); + ASSERT_NE(macro_grid_, nullptr) << "This should not happen!"; + ASSERT_NE(multiscale_grid_, nullptr) << "This should not happen!"; + + multiscale_grid_->visualize(Expectations::id()); + } // ... visualize_is_callable(...) + + void check_intersection_orientation_for_equal_levels(const bool expect_failure = Expectations::failure_for_equal()) + { + setup(); + ASSERT_NE(macro_grid_, nullptr) << "This should not happen!"; + ASSERT_NE(multiscale_grid_, nullptr) << "This should not happen!"; + size_t failure = 0; + + for (ssize_t level = 0; level <= multiscale_grid_->max_local_level(0); ++level) + failure += check_intersections_for_levels(level, level, expect_failure); + + if (failure) + std::cout << "The actual numbers of broken intersections are\n" + << convert_to_initializer_list_str(count_wrong_intersections_on_all_levels()) << std::endl; + } // ... check_intersection_orientation_for_equal_levels(...) + + void check_intersection_orientation_for_higher_neighbor_levels( + const bool expect_failure = Expectations::failure_for_higher()) + { + setup(); + ASSERT_NE(macro_grid_, nullptr) << "This should not happen!"; + ASSERT_NE(multiscale_grid_, nullptr) << "This should not happen!"; + size_t failure = 0; + + for (ssize_t entity_level = 0; entity_level <= multiscale_grid_->max_local_level(0); ++entity_level) + for (ssize_t neighbor_level = entity_level; neighbor_level <= multiscale_grid_->max_local_level(0); + ++neighbor_level) + failure += check_intersections_for_levels(entity_level, neighbor_level, expect_failure); + + if (failure) + std::cout << "The actual numbers of broken intersections are\n" + << convert_to_initializer_list_str(count_wrong_intersections_on_all_levels()) << std::endl; + } // ... check_intersection_orientation_for_higher_neighbor_levels(...) + + void check_intersection_orientation_for_lower_or_equal_neighbor_levels( + const bool expect_failure = Expectations::failure_for_lower_or_equal()) + { + setup(); + ASSERT_NE(macro_grid_, nullptr) << "This should not happen!"; + ASSERT_NE(multiscale_grid_, nullptr) << "This should not happen!"; + size_t failure = 0; + + for (ssize_t entity_level = 0; entity_level <= multiscale_grid_->max_local_level(0); ++entity_level) + for (ssize_t neighbor_level = 0; neighbor_level <= entity_level; ++neighbor_level) + failure += check_intersections_for_levels(entity_level, neighbor_level, expect_failure); + + if (failure) + std::cout << "The actual numbers of broken intersections are\n" + << convert_to_initializer_list_str(count_wrong_intersections_on_all_levels()) << std::endl; + } // ... check_intersection_orientation_for_lower_or_equal_neighbor_levels(...) + + size_t + check_intersections_for_levels(const ssize_t entity_level, const ssize_t neighbor_level, const bool expect_failure) + { + size_t failure = 0; + const auto actual_num_wrongly_oriented_intersections = count_wrong_intersections(entity_level, neighbor_level); + if (expect_failure) { + const auto expected_results = Expectations::results(); + const auto search_for_levels_in_expected_results = + expected_results.find(std::make_pair(entity_level, neighbor_level)); + EXPECT_NE(search_for_levels_in_expected_results, expected_results.end()) + << "missing expected results for entity and neighbor level " + << convert_to_initializer_list_str(std::make_pair(entity_level, neighbor_level)) << ".\n" + << "-> PLEASE ADD A RECORD TO\n" + << " ExpectedResults<" << Common::Typename<MacroGridType>::value() << ", " + << Common::Typename<LocalGridType>::value() << ">!"; + if (search_for_levels_in_expected_results == expected_results.end()) { + DUNE_THROW(InvalidStateException, + "Cannot use ASSERT_EQ above, so we need to exit this way.\n\n" + << "The actual numbers of broken intersections are\n" + << convert_to_initializer_list_str(count_wrong_intersections_on_all_levels())); + } + const auto expected_num_wrongly_oriented_intersections = search_for_levels_in_expected_results->second; + if (expected_num_wrongly_oriented_intersections != actual_num_wrongly_oriented_intersections) + ++failure; + EXPECT_EQ(expected_num_wrongly_oriented_intersections, actual_num_wrongly_oriented_intersections) + << "\nTHIS IS A GOOD THING, AN ACTUAL NUMBER OF FAILURES WHICH IS LOWER THAN THE EXPECTED NUMBER OF " + << "FAILURES IS AN IMPROVEMENT!\n" + << "-> BE HAPPY AND UPDATE THE RECORD IN!\n" + << " ExpectedResults<" << Common::Typename<MacroGridType>::value() << ", " + << Common::Typename<LocalGridType>::value() << ">!\n" + << "IF THE UPDATED RECORDS DO NOT INDICATE FAILURES ANYMORE, UPDATE THE TESTS!\n"; + } else { + if (actual_num_wrongly_oriented_intersections != 0) + ++failure; + EXPECT_EQ(actual_num_wrongly_oriented_intersections, 0); + } + return failure; + } // ... check_intersections(...) + + std::map<std::pair<size_t, size_t>, size_t> count_wrong_intersections_on_all_levels() + { + setup(); + if (!macro_grid_) + DUNE_THROW(InvalidStateException, "This should not happen!"); // <- cannot use ASSERT_NE here, non void return + if (!multiscale_grid_) + DUNE_THROW(InvalidStateException, "This should not happen!"); // <- s.a. + + std::map<std::pair<size_t, size_t>, size_t> results; + for (ssize_t entity_level = 0; entity_level <= multiscale_grid_->max_local_level(0); ++entity_level) + for (ssize_t neighbor_level = 0; neighbor_level <= multiscale_grid_->max_local_level(0); ++neighbor_level) + results[std::make_pair(entity_level, neighbor_level)] = count_wrong_intersections(entity_level, neighbor_level); + return results; + } + + size_t count_wrong_intersections(const size_t entity_level, const size_t neighbor_level) + { + setup(); + if (!macro_grid_) + DUNE_THROW(InvalidStateException, "This should not happen!"); // <- cannot use ASSERT_NE here, non void return + if (!multiscale_grid_) + DUNE_THROW(InvalidStateException, "This should not happen!"); // <- s.a. + + const auto& macro_grid_view = multiscale_grid_->macro_grid_view(); + size_t failures = 0; + for (auto&& macro_entity : Dune::elements(macro_grid_view)) { + for (auto&& macro_intersection : Dune::intersections(macro_grid_view, macro_entity)) { + if (macro_intersection.neighbor() && !macro_intersection.boundary()) { + const auto macro_neighbor = macro_intersection.outside(); + // const auto local_grid_view = multiscale_grid_->local_grid(macro_entity).level_view(entity_level); + const auto& coupling_glue = + multiscale_grid_->coupling(macro_entity, + entity_level, + macro_neighbor, + neighbor_level, + /*allow_for_broken_orientation_of_coupling_intersections=*/true); + failures += DD::check_for_broken_coupling_intersections(coupling_glue); + } + } + } + return failures; + } // ... count_wrong_intersections(...) + + std::unique_ptr<GridProvider<MacroGridType>> macro_grid_; + std::unique_ptr<DD::Glued<MacroGridType, LocalGridType>> multiscale_grid_; +}; // struct GluedMultiscaleGridTest + + +} // namespace Grid +} // namespace XT +} // namespace Dune + +#endif // HAVE_DUNE_GRID_GLUE +#endif // DUNE_XT_GRID_TEST_DD_GLUED_HH diff --git a/dune/xt/grid/test/dd_glued_2d.cc b/dune/xt/grid/test/dd_glued_2d.cc new file mode 100644 index 0000000000000000000000000000000000000000..168d93b0eb054fdf303de620a75ef2ce561d2e59 --- /dev/null +++ b/dune/xt/grid/test/dd_glued_2d.cc @@ -0,0 +1,350 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2017) + +#include <dune/xt/common/test/main.hxx> // <- this one has to come first (includes the config.h)! + +#include "dd_glued.hh" + +#if HAVE_DUNE_GRID_GLUE + +namespace Dune { +namespace XT { +namespace Grid { + + +template <bool anything> +struct ExpectedResults<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, + YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, + anything> +{ + static int num_coarse_refinements() + { + return 0; + } + + static int num_local_refinements() + { + return 2; + } + + static std::string id() + { + return "2d_yaspgrid_yaspgrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { + return {4}; + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return false; + } + + static bool failure_for_higher() + { + return false; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { + return {{{0, 0}, 0}, + {{0, 1}, 0}, + {{0, 2}, 0}, + {{1, 0}, 24}, + {{1, 1}, 0}, + {{1, 2}, 0}, + {{2, 0}, 72}, + {{2, 1}, 48}, + {{2, 2}, 0}}; + } +}; // struct ExpectedResults<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, YaspGrid<2, +// EquidistantOffsetCoordinates<double, 2>>, anything> + +#if HAVE_DUNE_ALUGRID || HAVE_ALUGRID + +template <class Comm, bool anything> +struct ExpectedResults<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, + ALUGrid<2, 2, simplex, conforming, Comm>, + anything> +{ + static int num_coarse_refinements() + { +#if HAVE_ALUGRID + return 0; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return -1; +#endif + } + + static int num_local_refinements() + { +#if HAVE_ALUGRID + return 3; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return -1; +#endif + } + + static std::string id() + { + return "2d_yaspgrid_alugridconforming"; + } + + static std::set<int> num_local_couplings_intersections() + { +#if HAVE_ALUGRID + return {2}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return {0}; +#endif + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return false; + } + + static bool failure_for_higher() + { + return false; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { +#if HAVE_ALUGRID + return {{{0, 0}, 0}, + {{0, 1}, 0}, + {{0, 2}, 0}, + {{0, 3}, 0}, + {{1, 0}, 0}, + {{1, 1}, 0}, + {{1, 2}, 0}, + {{1, 3}, 0}, + {{2, 0}, 24}, + {{2, 1}, 24}, + {{2, 2}, 0}, + {{2, 3}, 0}, + {{3, 0}, 24}, + {{3, 1}, 24}, + {{3, 2}, 0}, + { {3, 3}, + 0 }}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return {}; +#endif + } +}; // struct ExpectedResults<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, ALUGrid<2, 2, simplex, conforming, +// Comm>, anything> + +template <class Comm, bool anything> +struct ExpectedResults<ALUGrid<2, 2, simplex, conforming, Comm>, ALUGrid<2, 2, simplex, conforming, Comm>, anything> +{ + static int num_coarse_refinements() + { +#if HAVE_ALUGRID + return 1; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return -1; +#endif + } + + static int num_local_refinements() + { +#if HAVE_ALUGRID + return 3; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return -1; +#endif + } + + static std::string id() + { + return "2d_alugridcoforming_alugridconforming"; + } + + static std::set<size_t> num_local_couplings_intersections() + { +#if HAVE_ALUGRID + return {2, 4}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return false; + } + + static bool failure_for_higher() + { + return false; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { +#if HAVE_ALUGRID + return {{{0, 0}, 0}, + {{0, 1}, 0}, + {{0, 2}, 0}, + {{0, 3}, 0}, + {{1, 0}, 24}, + {{1, 1}, 0}, + {{1, 2}, 0}, + {{1, 3}, 0}, + {{2, 0}, 96}, + {{2, 1}, 72}, + {{2, 2}, 0}, + {{2, 3}, 0}, + {{3, 0}, 144}, + {{3, 1}, 120}, + {{3, 2}, 48}, + { {3, 3}, + 0 }}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return {}; +#endif + } +}; // struct ExpectedResults<ALUGrid<2, 2, simplex, conforming, Comm>, ALUGrid<2, 2, simplex, conforming, Comm>, +// anything> + +#endif // HAVE_DUNE_ALUGRID || HAVE_ALUGRID +#if HAVE_UG + +template <bool anything> +struct ExpectedResults<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, UGGrid<2>, anything> +{ + static int num_coarse_refinements() + { + return 0; + } + + static int num_local_refinements() + { + return 2; + } + + static std::string id() + { + return "2d_yaspgrid_uggrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { + return {4}; + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return false; + } + + static bool failure_for_higher() + { + return false; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { + return {{{0, 0}, 0}, + {{0, 1}, 0}, + {{0, 2}, 0}, + {{1, 0}, 24}, + {{1, 1}, 0}, + {{1, 2}, 0}, + {{2, 0}, 72}, + {{2, 1}, 48}, + {{2, 2}, 0}}; + } +}; // struct ExpectedResults<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, UGGrid<2>, anything> + +#endif // HAVE_UG + + +} // namespace Grid +} // namespace XT +} // namespace Dune + + +using namespace Dune; +using namespace Dune::XT::Grid; + + +// clang-format off +typedef ::testing::Types< std::tuple<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, + YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>> +#if HAVE_ALUGRID + , std::tuple<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, + ALUGrid<2, 2, simplex, conforming>> + , std::tuple<ALUGrid<2, 2, simplex, conforming>, ALUGrid<2, 2, simplex, conforming>> +#endif // HAVE_ALUGRID +#if HAVE_UG + , std::tuple<YaspGrid<2, EquidistantOffsetCoordinates<double, 2>>, UGGrid<2>> +#endif + > GridTypes; // clang-format on + +TYPED_TEST_CASE(GluedMultiscaleGridTest, GridTypes); +TYPED_TEST(GluedMultiscaleGridTest, setup_works) +{ + this->setup(); +} +TYPED_TEST(GluedMultiscaleGridTest, visualize_is_callable) +{ + this->visualize_is_callable(); +} +TYPED_TEST(GluedMultiscaleGridTest, couplings_are_of_correct_size) +{ + this->couplings_are_of_correct_size(); +} +TYPED_TEST(GluedMultiscaleGridTest, intersections_are_correctly_oriented_for_equal_levels) +{ + this->check_intersection_orientation_for_equal_levels(); +} +TYPED_TEST(GluedMultiscaleGridTest, intersections_are_correctly_oriented_for_higher_neighbor_levels) +{ + this->check_intersection_orientation_for_higher_neighbor_levels(); +} +TYPED_TEST(GluedMultiscaleGridTest, + __STILL_BROKEN__intersection_orientation_is_wrong_for_lower_or_equal_neighbor_levels) +{ + this->check_intersection_orientation_for_lower_or_equal_neighbor_levels(); +} + + +#endif // HAVE_DUNE_GRID_GLUE diff --git a/dune/xt/grid/test/dd_glued_3d.cc b/dune/xt/grid/test/dd_glued_3d.cc new file mode 100644 index 0000000000000000000000000000000000000000..3f7be847c595e3a6d4d51152b51f0b61fad3a201 --- /dev/null +++ b/dune/xt/grid/test/dd_glued_3d.cc @@ -0,0 +1,502 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2017) + +#include <dune/xt/common/test/main.hxx> // <- this one has to come first (includes the config.h)! + +#include "dd_glued.hh" + +#if HAVE_DUNE_GRID_GLUE + +namespace Dune { +namespace XT { +namespace Grid { + + +template <bool anything> +struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + anything> +{ + static int num_coarse_refinements() + { + return 0; + } + + static int num_local_refinements() + { + return 2; + } + + static std::string id() + { + return "3d_yaspgrid_yaspgrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { + // we expect 16 rectangles, each containing two triangles + return {32}; + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return true; + } + + static bool failure_for_higher() + { + return true; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { + return {{{0, 0}, 216}, + {{0, 1}, 864}, + {{0, 2}, 3456}, + {{1, 0}, 108}, + {{1, 1}, 864}, + {{1, 2}, 3456}, + {{2, 0}, 108}, + {{2, 1}, 108}, + {{2, 2}, 3456}}; + } +}; // struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, YaspGrid<3, +// EquidistantOffsetCoordinates<double, 3>>, anything> + +#if HAVE_DUNE_ALUGRID || HAVE_ALUGRID + +template <class Comm, bool anything> +struct ExpectedResults<ALUGrid<3, 3, cube, conforming, Comm>, + YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + anything> +{ + static int num_coarse_refinements() + { +#if HAVE_ALUGRID + return 0; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static int num_local_refinements() + { +#if HAVE_ALUGRID + return 2; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static std::string id() + { + return "3d_alucubeconforminggrid_yaspgrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { +#if HAVE_ALUGRID + return {32}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return true; + } + + static bool failure_for_higher() + { + return true; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { +#if HAVE_ALUGRID + return {{{0, 0}, 216}, + {{0, 1}, 864}, + {{0, 2}, 3456}, + {{1, 0}, 108}, + {{1, 1}, 864}, + {{1, 2}, 3456}, + {{2, 0}, 108}, + {{2, 1}, 108}, + { {2, 2}, + 3456 }}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return {}; +#endif + } +}; // struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, ALUGrid<3, 3, simplex, nonconforming, +// Comm>, anything> + +template <class Comm, bool anything> +struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + ALUGrid<3, 3, cube, nonconforming, Comm>, + anything> +{ + static int num_coarse_refinements() + { +#if HAVE_ALUGRID + return 0; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static int num_local_refinements() + { +#if HAVE_ALUGRID + return 2; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static std::string id() + { + return "3d_alucubenonconforminggrid_yaspgrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { +#if HAVE_ALUGRID + return {32}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return true; + } + + static bool failure_for_higher() + { + return true; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { +#if HAVE_ALUGRID + return {{{0, 0}, 216}, + {{0, 1}, 864}, + {{0, 2}, 3456}, + {{1, 0}, 108}, + {{1, 1}, 864}, + {{1, 2}, 3456}, + {{2, 0}, 108}, + {{2, 1}, 108}, + { {2, 2}, + 3456 }}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return {}; +#endif + } +}; // ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, ALUGrid<3, 3, cube, nonconforming, Comm>, +// anything> + +template <class Comm, bool anything> +struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + ALUGrid<3, 3, simplex, nonconforming, Comm>, + anything> +{ + static int num_coarse_refinements() + { +#if HAVE_ALUGRID + return 0; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static int num_local_refinements() + { +#if HAVE_ALUGRID + return 2; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static std::string id() + { + return "3d_alusimplexnonconforminggrid_yaspgrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { +#if HAVE_ALUGRID + return {32}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return true; + } + + static bool failure_for_higher() + { + return true; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { +#if HAVE_ALUGRID + return {{{0, 0}, 216}, + {{0, 1}, 864}, + {{0, 2}, 3456}, + {{1, 0}, 108}, + {{1, 1}, 864}, + {{1, 2}, 3456}, + {{2, 0}, 108}, + {{2, 1}, 108}, + { {2, 2}, + 3456 }}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return {}; +#endif + } +}; // struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, ALUGrid<3, 3, simplex, nonconforming, +// Comm>, anything> + +template <class Comm, bool anything> +struct ExpectedResults<ALUGrid<3, 3, cube, conforming, Comm>, ALUGrid<3, 3, simplex, nonconforming, Comm>, anything> +{ + static int num_coarse_refinements() + { +#if HAVE_ALUGRID + return 0; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static int num_local_refinements() + { +#if HAVE_ALUGRID + return 2; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static std::string id() + { + return "3d_alucubeconforminggrid_alusimplexnonconforminggrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { +#if HAVE_ALUGRID + return {32}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return 0; +#endif + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return true; + } + + static bool failure_for_higher() + { + return true; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { +#if HAVE_ALUGRID + return {{{0, 0}, 216}, + {{0, 1}, 864}, + {{0, 2}, 3456}, + {{1, 0}, 108}, + {{1, 1}, 864}, + {{1, 2}, 3456}, + {{2, 0}, 108}, + {{2, 1}, 108}, + { {2, 2}, + 3456 }}; +#else + DUNE_THROW(InvalidStateException, "Please update these for dune-alugrid!"); + return {}; +#endif + } +}; // struct ExpectedResults<ALUGrid<3, 3, cube, conforming, Comm>, ALUGrid<3, 3, simplex, nonconforming, Comm>, +// anything> + +#endif // HAVE_DUNE_ALUGRID || HAVE_ALUGRID +#if HAVE_UG + +template <bool anything> +struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, UGGrid<3>, anything> +{ + static int num_coarse_refinements() + { + return 0; + } + + static int num_local_refinements() + { + return 2; + } + + static std::string id() + { + return "3d_yaspgrid_uggrid"; + } + + static std::set<size_t> num_local_couplings_intersections() + { + // we expect 16 rectangles, each containing two triangles + return {32}; + } + + static bool failure_for_lower_or_equal() + { + return true; + } + + static bool failure_for_equal() + { + return true; + } + + static bool failure_for_higher() + { + return true; + } + + static std::map<std::pair<size_t, size_t>, size_t> results() + { + return {{{0, 0}, 216}, + {{0, 1}, 864}, + {{0, 2}, 3456}, + {{1, 0}, 72}, + {{1, 1}, 864}, + {{1, 2}, 3456}, + {{2, 0}, 72}, + {{2, 1}, 36}, + {{2, 2}, 3456}}; + } +}; // struct ExpectedResults<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, UGGrid<3>, anything> + + +#endif // HAVE_UG + +} // namespace Grid +} // namespace XT +} // namespace Dune + + +using namespace Dune; +using namespace Dune::XT::Grid; + + +// clang-format off +typedef ::testing::Types< std::tuple<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>> +#if HAVE_ALUGRID +// , std::tuple<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, +// ALUGrid<3, 3, cube, conforming>> // <- knwon to fail completely + , std::tuple<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + ALUGrid<3, 3, cube, nonconforming>> +// , std::tuple<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, +// ALUGrid<3, 3, simplex, conforming>> // <- knwon to fail completely + , std::tuple<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, + ALUGrid<3, 3, simplex, nonconforming>> + , std::tuple<ALUGrid<3, 3, cube, conforming>, + YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>> +// , std::tuple<ALUGrid<3, 3, simplex, conforming>, +// ALUGrid<3, 3, simplex, nonconforming>> // <- knwon to fail completely +// , std::tuple<ALUGrid<3, 3, simplex, nonconforming>, +// ALUGrid<3, 3, simplex, nonconforming>> // <- knwon to fail completely + , std::tuple<ALUGrid<3, 3, cube, conforming>, ALUGrid<3, 3, simplex, nonconforming>> +#endif // HAVE_ALUGRID +#if HAVE_UG + , std::tuple<YaspGrid<3, EquidistantOffsetCoordinates<double, 3>>, UGGrid<3>> +#endif + > GridTypes; // clang-format on + +TYPED_TEST_CASE(GluedMultiscaleGridTest, GridTypes); +TYPED_TEST(GluedMultiscaleGridTest, setup_works) +{ + this->setup(); +} +TYPED_TEST(GluedMultiscaleGridTest, visualize_is_callable) +{ + this->visualize_is_callable(); +} +TYPED_TEST(GluedMultiscaleGridTest, couplings_are_of_correct_size) +{ + this->couplings_are_of_correct_size(); +} +TYPED_TEST(GluedMultiscaleGridTest, __STILL_BROKEN__intersections_are_correctly_oriented_for_equal_levels) +{ + this->check_intersection_orientation_for_equal_levels(); +} +TYPED_TEST(GluedMultiscaleGridTest, __STILL_BROKEN__intersections_are_correctly_oriented_for_higher_neighbor_levels) +{ + this->check_intersection_orientation_for_higher_neighbor_levels(); +} +TYPED_TEST(GluedMultiscaleGridTest, + __STILL_BROKEN__intersection_orientation_is_wrong_for_lower_or_equal_neighbor_levels) +{ + this->check_intersection_orientation_for_lower_or_equal_neighbor_levels(); +} + + +#endif // HAVE_DUNE_GRID_GLUE diff --git a/dune/xt/grid/type_traits.hh b/dune/xt/grid/type_traits.hh index b8a2fad25458e27156da2bb1cd6079289830b365..5585330eb57e9cf5965dae7163d32fc475ccb306 100644 --- a/dune/xt/grid/type_traits.hh +++ b/dune/xt/grid/type_traits.hh @@ -68,14 +68,14 @@ struct is_grid<Dune::ALUGrid<dim, dimworld, elType, refineType, Comm>> : public #endif // HAVE_DUNE_ALUGRID -#if HAVE_DUNE_UGGRID +#if HAVE_DUNE_UGGRID || HAVE_UG template <int dim> struct is_grid<Dune::UGGrid<dim>> : public std::true_type { }; -#endif // HAVE_DUNE_UGGRID +#endif // HAVE_DUNE_UGGRID || HAVE_UG #if HAVE_DUNE_SPGRID template <class ct, int dim, template <int> class Ref, class Comm> diff --git a/dune/xt/grid/walker.pbh b/dune/xt/grid/walker.pbh new file mode 100644 index 0000000000000000000000000000000000000000..1df282674c5b5f1823cae838249ead54649b20c9 --- /dev/null +++ b/dune/xt/grid/walker.pbh @@ -0,0 +1,50 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2016 - 2017) + +#ifndef DUNE_XT_GRID_WALKER_PBH +#define DUNE_XT_GRID_WALKER_PBH +#if HAVE_DUNE_PYBINDXI + +#include <dune/pybindxi/pybind11.h> + +#include "type_traits.hh" +#include "walker.hh" + +namespace Dune { +namespace XT { +namespace Grid { + + +template <class G> +typename std::enable_if<is_layer<G>::value, pybind11::class_<Walker<G>>>::type bind_Walker(pybind11::module& m, + const std::string& grid_id) +{ + namespace py = pybind11; + using namespace pybind11::literals; + + typedef Walker<G> C; + + py::class_<Walker<G>> c(m, std::string("Walker__" + grid_id).c_str(), std::string("Walker__" + grid_id).c_str()); + + c.def("clear", &C::clear); + // c.def("append", [](C& self, C& other) { self.append(other); }, "other"_a, py::keep_alive<1,2>()); + c.def("prepare", &C::prepare); + c.def("finalize", &C::finalize); + c.def("walk", [](C& self, const bool use_tbb) { self.walk(use_tbb); }, "use_tbb"_a = false); + + return c; +} // ... bind_BoundaryType(...) + + +} // namespace Grid +} // namespace XT +} // namespace Dune + +#endif // HAVE_DUNE_PYBINDXI +#endif // DUNE_XT_GRID_WALKER_PBH diff --git a/dune/xt/grid/walker/apply-on.pbh b/dune/xt/grid/walker/apply-on.pbh new file mode 100644 index 0000000000000000000000000000000000000000..e4203e08d39da6dc20bd30253d0414aabcf00e32 --- /dev/null +++ b/dune/xt/grid/walker/apply-on.pbh @@ -0,0 +1,135 @@ +// This file is part of the dune-xt-grid project: +// https://github.com/dune-community/dune-xt-grid +// The copyright lies with the authors of this file (see below). +// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +// or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +// with "runtime exception" (http://www.dune-project.org/license.html) +// Authors: +// Felix Schindler (2017) + +#ifndef DUNE_XT_GRID_WALKER_APPLYON_PBH +#define DUNE_XT_GRID_WALKER_APPLYON_PBH +#if HAVE_DUNE_PYBINDXI + +#include <sstream> +#include <type_traits> + +#include <dune/pybindxi/pybind11.h> +#include <dune/pybindxi/operators.h> + +#include <dune/xt/grid/intersection.hh> +#include <dune/xt/grid/gridprovider/provider.hh> +#include <dune/xt/grid/boundaryinfo/interfaces.hh> + +#include "apply-on.hh" + +namespace Dune { +namespace XT { +namespace Grid { +namespace internal { + + +template <class C, class BI, bool with_bi = false> +struct bind_WhichIntersection +{ + pybind11::class_<C> operator()(pybind11::module& m, const std::string& class_id, const std::string& method_id) + { + namespace py = pybind11; + using namespace pybind11::literals; + + py::class_<C, ApplyOn::WhichIntersection<typename C::GridViewType>> c(m, class_id.c_str(), class_id.c_str()); + c.def(py::init<>()); + + m.def(method_id.c_str(), [](const BI& /*boundary_info*/) { return C(); }, "boundary_info"_a); + + return c; + } +}; // struct bind_WhichIntersection + +template <class C, class BI> +struct bind_WhichIntersection<C, BI, true> +{ + pybind11::class_<C> operator()(pybind11::module& m, const std::string& class_id, const std::string& method_id) + { + namespace py = pybind11; + using namespace pybind11::literals; + + py::class_<C, ApplyOn::WhichIntersection<typename C::GridViewType>> c(m, class_id.c_str(), class_id.c_str()); + c.def(py::init<const BI&>()); + + m.def(method_id.c_str(), [](const BI& boundary_info) { return C(boundary_info); }, "boundary_info"_a); + + return c; + } +}; // struct bind_WhichIntersection<..., true> + + +} // namespace internal + + +template <class GV> +void addbind_WhichIntersection(pybind11::module& m, const std::string& grid_id, const std::string& layer) +{ + namespace py = pybind11; + using namespace pybind11::literals; + + py::class_<ApplyOn::WhichIntersection<GV>>(m, + std::string("ApplyOnWhichIntersection__" + grid_id + "_" + layer).c_str(), + std::string("ApplyOnWhichIntersection__" + grid_id + "_" + layer).c_str()); + + internal::bind_WhichIntersection<ApplyOn::AllIntersections<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + false>()( + m, "ApplyOnAllIntersections__" + grid_id + "_" + layer, "make_apply_on_all_intersections__" + layer); + + internal::bind_WhichIntersection<ApplyOn::InnerIntersections<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + false>()( + m, "ApplyOnInnerIntersections__" + grid_id + "_" + layer, "make_apply_on_inner_intersections__" + layer); + + internal::bind_WhichIntersection<ApplyOn::InnerIntersectionsPrimally<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + false>()(m, + "ApplyOnInnerIntersectionsPrimally__" + grid_id + "_" + layer, + "make_apply_on_inner_intersections_primally__" + layer); + + internal::bind_WhichIntersection<ApplyOn::BoundaryIntersections<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + false>()( + m, "ApplyOnBoundaryIntersections__" + grid_id + "_" + layer, "make_apply_on_boundary_intersections__" + layer); + + internal::bind_WhichIntersection<ApplyOn::NonPeriodicBoundaryIntersections<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + false>()(m, + "ApplyOnNonPeriodicBoundaryIntersections__" + grid_id + "_" + layer, + "make_apply_on_nonperiodic_boundary_intersections__" + layer); + + internal::bind_WhichIntersection<ApplyOn::PeriodicIntersections<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + false>()( + m, "ApplyOnPeriodicIntersections__" + grid_id + "_" + layer, "make_apply_on_periodic_intersections__" + layer); + + internal::bind_WhichIntersection<ApplyOn::PeriodicIntersectionsPrimally<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + false>()(m, + "ApplyOnPeriodicIntersectionsPrimally__" + grid_id + "_" + layer, + "make_apply_on_periodic_intersections_primally__" + layer); + + internal::bind_WhichIntersection<ApplyOn::DirichletIntersections<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + true>()( + m, "ApplyOnDirichletIntersections__" + grid_id + "_" + layer, "make_apply_on_dirichlet_intersections__" + layer); + + internal::bind_WhichIntersection<ApplyOn::NeumannIntersections<GV>, + BoundaryInfo<typename Intersection<GV>::type>, + true>()( + m, "ApplyOnNeumannIntersections__" + grid_id + "_" + layer, "make_apply_on_neumann_intersections__" + layer); +} // ... addbind_WhichIntersection(...) + + +} // namespace Grid +} // namespace XT +} // namespace Dune + +#endif // HAVE_DUNE_PYBINDXI +#endif // DUNE_XT_GRID_WALKER_APPLYON_PBH diff --git a/python/dune/__init__.py b/python/dune/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..28ee14c19a1fea9a5c3382c496c210f14055bcd2 --- /dev/null +++ b/python/dune/__init__.py @@ -0,0 +1,10 @@ +# This file is part of the dune-xt-grid project: +# https://github.com/dune-community/dune-xt-grid +# The copyright lies with the authors of this file (see below). +# License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +# or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +# with "runtime exception" (http://www.dune-project.org/license.html) +# Authors: +# Felix Schindler (2016 - 2017) + +__import__('pkg_resources').declare_namespace(__name__) diff --git a/python/dune/xt/__init__.py b/python/dune/xt/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..28ee14c19a1fea9a5c3382c496c210f14055bcd2 --- /dev/null +++ b/python/dune/xt/__init__.py @@ -0,0 +1,10 @@ +# This file is part of the dune-xt-grid project: +# https://github.com/dune-community/dune-xt-grid +# The copyright lies with the authors of this file (see below). +# License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +# or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +# with "runtime exception" (http://www.dune-project.org/license.html) +# Authors: +# Felix Schindler (2016 - 2017) + +__import__('pkg_resources').declare_namespace(__name__) diff --git a/python/dune/xt/grid/__init__.py b/python/dune/xt/grid/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..64ae0bd7eddcc2f66187bb6d9a64c6d9a1e5f4dc --- /dev/null +++ b/python/dune/xt/grid/__init__.py @@ -0,0 +1,41 @@ +# This file is part of the dune-xt-grid project: +# https://github.com/dune-community/dune-xt-grid +# The copyright lies with the authors of this file (see below). +# License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) +# or GPL-2.0+ (http://opensource.org/licenses/gpl-license) +# with "runtime exception" (http://www.dune-project.org/license.html) +# Authors: +# Felix Schindler (2017) + +from importlib import import_module + + +def init_logger(max_info_level=-1, + max_debug_level=-1, + enable_warnings=True, + enable_colors=True, + info_color='blue', + debug_color='darkgray', + warning_color='red'): + from ._grid import init_logger as _init_logger + initializers = [_init_logger] + for module_name in ('xt.common', 'xt.la', 'xt.functions', 'gdt'): + try: + mm = import_module('dune.{}'.format(module_name)) + initializers.append(mm.init_logger) + except ModuleNotFoundError: + pass + for initializer in initializers: + initializer(max_info_level, max_debug_level, enable_warnings, enable_colors, info_color, debug_color, + warning_color) + + +def init_mpi(args=list()): + try: + from dune.gdt import init_mpi as _init_mpi + except ModuleNotFoundError: + from dune.xt.common import init_mpi as _init_mpi + _init_mpi(args) + + +from ._grid import *