3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
geometry/intersectionentityset.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_GEOMETRY_INTERSECTION_ENTITY_SET_HH
26#define DUMUX_GEOMETRY_INTERSECTION_ENTITY_SET_HH
27
28#include <algorithm>
29#include <cassert>
30#include <iostream>
31#include <memory>
32#include <type_traits>
33#include <utility>
34#include <vector>
35
36#include <dune/common/indices.hh>
37#include <dune/common/timer.hh>
38#include <dune/common/iteratorrange.hh>
39#include <dune/common/promotiontraits.hh>
40#include <dune/common/reservedvector.hh>
41#include <dune/geometry/affinegeometry.hh>
42#include <dune/geometry/type.hh>
43
46
47namespace Dumux {
48
53template<class DomainEntitySet, class TargetEntitySet>
55{
58
59 using ctype = typename Dune::PromotionTraits<typename DomainEntitySet::ctype, typename TargetEntitySet::ctype>::PromotedType;
60
61 static constexpr int dimWorld = DomainEntitySet::dimensionworld;
62 static_assert(dimWorld == int(TargetEntitySet::dimensionworld), "Entity sets must have the same world dimension");
63
64 using GlobalPosition = Dune::FieldVector<ctype, dimWorld>;
65
66 static constexpr int dimDomain = DomainEntitySet::Entity::Geometry::mydimension;
67 static constexpr int dimTarget = TargetEntitySet::Entity::Geometry::mydimension;
68 static constexpr bool isMixedDimensional = dimDomain != dimTarget;
69
73 class IntersectionEntity
74 {
75 static constexpr int dimIs = std::min(dimDomain, dimTarget);
76 using Geometry = Dune::AffineGeometry<ctype, dimIs, dimWorld>; // geometries are always simplices
77
78 // we can only have multiple neighbors in the mixeddimensional case and then only for the side with the largest dimension
79 using IndexStorage = std::pair<std::conditional_t<dimDomain <= dimTarget, Dune::ReservedVector<std::size_t, 1>, std::vector<std::size_t>>,
80 std::conditional_t<dimTarget <= dimDomain, Dune::ReservedVector<std::size_t, 1>, std::vector<std::size_t>>>;
81
82 static constexpr auto domainIdx = Dune::index_constant<0>{};
83 static constexpr auto targetIdx = Dune::index_constant<1>{};
84
85 public:
86 IntersectionEntity(const DomainTree& domainTree, const TargetTree& targetTree)
87 : domainTree_(domainTree)
88 , targetTree_(targetTree)
89 {}
90
92 void setCorners(const std::vector<GlobalPosition>& corners)
93 {
94 corners_ = corners;
95 assert(corners.size() == dimIs + 1); // Only simplex intersections are allowed
96 }
97
99 void addNeighbors(std::size_t domain, std::size_t target)
100 {
101 if (numDomainNeighbors() == 0 && numTargetNeighbors() == 0)
102 {
103 std::get<domainIdx>(neighbors_).push_back(domain);
104 std::get<targetIdx>(neighbors_).push_back(target);
105 }
106 else if (dimDomain > dimTarget)
107 std::get<domainIdx>(neighbors_).push_back(domain);
108
109 else if (dimTarget > dimDomain)
110 std::get<targetIdx>(neighbors_).push_back(target);
111
112 else
113 DUNE_THROW(Dune::InvalidStateException, "Cannot add more than one neighbor per side for equidimensional intersection!");
114 }
115
117 Geometry geometry() const
118 { return Geometry(Dune::GeometryTypes::simplex(dimIs), corners_); }
119
121 std::size_t numDomainNeighbors() const
122 { return std::get<domainIdx>(neighbors_).size(); }
123
125 std::size_t numTargetNeighbors() const
126 { return std::get<targetIdx>(neighbors_).size(); }
127
129 typename DomainEntitySet::Entity domainEntity(unsigned int n = 0) const
130 { return domainTree_.entitySet().entity(std::get<domainIdx>(neighbors_)[n]); }
131
133 typename TargetEntitySet::Entity targetEntity(unsigned int n = 0) const
134 { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); }
135
136 private:
137 IndexStorage neighbors_;
138 std::vector<GlobalPosition> corners_;
139
140 const DomainTree& domainTree_;
141 const TargetTree& targetTree_;
142 };
143
144 using Intersections = std::vector<IntersectionEntity>;
145
146public:
148 using Entity = IntersectionEntity;
150 using EntityIterator = typename Intersections::const_iterator;
151
156
160 void build(std::shared_ptr<const DomainEntitySet> domainSet, std::shared_ptr<const TargetEntitySet> targetSet)
161 {
162 domainTree_ = std::make_shared<DomainTree>(domainSet);
163 targetTree_ = std::make_shared<TargetTree>(targetSet);
164 build(*domainTree_, *targetTree_);
165 }
166
170 void build(std::shared_ptr<const DomainTree> domainTree, std::shared_ptr<const TargetTree> targetTree)
171 {
172 // make sure the tree don't get out of scope
173 domainTree_ = domainTree;
174 targetTree_ = targetTree;
175 build(*domainTree_, *targetTree_);
176 }
177
182 void build(const DomainTree& domainTree, const TargetTree& targetTree)
183 {
184 Dune::Timer timer;
185
186 // compute raw intersections
187 const auto rawIntersections = intersectingEntities(domainTree, targetTree);
188
189 // create a map to check whether intersection geometries were already inserted
190 // Note that this can only occur if the grids have different dimensionality.
191 // If this is the case, we keep track of the intersections using the indices of the lower-
192 // dimensional entities which is identical for all geometrically identical intersections.
193 std::vector<std::vector<std::vector<GlobalPosition>>> intersectionMap;
194 std::vector<std::vector<std::size_t>> intersectionIndex;
195 if constexpr (isMixedDimensional)
196 {
197 const auto numLowDimEntities = dimTarget < dimDomain ? targetTree.entitySet().size()
198 : domainTree.entitySet().size();
199 intersectionMap.resize(numLowDimEntities);
200 intersectionIndex.resize(numLowDimEntities);
201 }
202
203 // reserve memory for storing the intersections. In case of grids of
204 // different dimensionality this might be an overestimate. We get rid
205 // of the overhead memory at the end of this function though.
206 intersections_.clear();
207 intersections_.reserve(rawIntersections.size());
208
209 for (const auto& rawIntersection : rawIntersections)
210 {
211 bool add = true;
212
213 // Check if intersection was already inserted.
214 // In this case we only add new neighbor information as the geometry is identical.
215 if constexpr (isMixedDimensional)
216 {
217 const auto lowDimNeighborIdx = getLowDimNeighborIdx_(rawIntersection);
218 for (int i = 0; i < intersectionMap[lowDimNeighborIdx].size(); ++i)
219 {
220 if (rawIntersection.cornersMatch(intersectionMap[lowDimNeighborIdx][i]))
221 {
222 add = false;
223 // only add the pair of neighbors using the insertionIndex
224 auto idx = intersectionIndex[lowDimNeighborIdx][i];
225 intersections_[idx].addNeighbors(rawIntersection.first(), rawIntersection.second());
226 break;
227 }
228 }
229 }
230
231 if(add)
232 {
233 // maybe add to the map
234 if constexpr (isMixedDimensional)
235 {
236 intersectionMap[getLowDimNeighborIdx_(rawIntersection)].push_back(rawIntersection.corners());
237 intersectionIndex[getLowDimNeighborIdx_(rawIntersection)].push_back(intersections_.size());
238 }
239
240 // add new intersection and add the neighbors
241 intersections_.emplace_back(domainTree, targetTree);
242 intersections_.back().setCorners(rawIntersection.corners());
243 intersections_.back().addNeighbors(rawIntersection.first(), rawIntersection.second());
244 }
245 }
246
247 intersections_.shrink_to_fit();
248 std::cout << "Computed " << size() << " intersection entities in " << timer.elapsed() << std::endl;
249 }
250
252 typename Intersections::const_iterator ibegin() const
253 { return intersections_.begin(); }
254
256 typename Intersections::const_iterator iend() const
257 { return intersections_.end(); }
258
260 std::size_t size() const
261 { return intersections_.size(); }
262
268 friend Dune::IteratorRange<EntityIterator> intersections(const IntersectionEntitySet& set)
269 { return {set.ibegin(), set.iend()}; }
270
271private:
272 template<class RawIntersection,
273 bool enable = isMixedDimensional, std::enable_if_t<enable, int> = 0>
274 auto getLowDimNeighborIdx_(const RawIntersection& is)
275 {
276 if constexpr (dimTarget < dimDomain)
277 return is.second();
278 else
279 return is.first();
280 }
281
282 Intersections intersections_;
283
284 std::shared_ptr<const DomainTree> domainTree_;
285 std::shared_ptr<const TargetTree> targetTree_;
286};
287
288} // end namespace Dumux
289
290#endif
std::vector< std::size_t > intersectingEntities(const Dune::FieldVector< ctype, dimworld > &point, const BoundingBoxTree< EntitySet > &tree, bool isCartesianGrid=false)
Compute all intersections between entities and a point.
Definition: geometry/intersectingentities.hh:100
Definition: adapt.hh:29
An axis-aligned bounding box volume tree implementation.
Definition: geometry/boundingboxtree.hh:66
const EntitySet & entitySet() const
the entity set this tree was built with
Definition: geometry/boundingboxtree.hh:133
A class representing the intersection entites two geometric entity sets.
Definition: geometry/intersectionentityset.hh:55
IntersectionEntity Entity
make intersection entity type available
Definition: geometry/intersectionentityset.hh:148
IntersectionEntitySet()=default
Default constructor.
void build(std::shared_ptr< const DomainEntitySet > domainSet, std::shared_ptr< const TargetEntitySet > targetSet)
Build intersections.
Definition: geometry/intersectionentityset.hh:160
std::size_t size() const
the number of intersections
Definition: geometry/intersectionentityset.hh:260
void build(const DomainTree &domainTree, const TargetTree &targetTree)
Build intersections.
Definition: geometry/intersectionentityset.hh:182
void build(std::shared_ptr< const DomainTree > domainTree, std::shared_ptr< const TargetTree > targetTree)
Build intersections.
Definition: geometry/intersectionentityset.hh:170
typename Intersections::const_iterator EntityIterator
make entity iterator type available
Definition: geometry/intersectionentityset.hh:150
friend Dune::IteratorRange< EntityIterator > intersections(const IntersectionEntitySet &set)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: geometry/intersectionentityset.hh:268
Intersections::const_iterator ibegin() const
return begin iterator to intersection container
Definition: geometry/intersectionentityset.hh:252
Intersections::const_iterator iend() const
return end iterator to intersection container
Definition: geometry/intersectionentityset.hh:256
An axis-aligned bounding box volume hierarchy for dune grids.
Algorithms that finds which geometric entites intersect.