Line | Branch | Exec | Source |
---|---|---|---|
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 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
6 | // | ||
7 | /*! | ||
8 | * \file | ||
9 | * \ingroup BoxDFMModel | ||
10 | * \brief Base class for the local finite volume geometry for the box discrete | ||
11 | * fracture model. | ||
12 | * | ||
13 | * This builds up the sub control volumes and sub control | ||
14 | * volume faces for an element. | ||
15 | */ | ||
16 | |||
17 | #ifndef DUMUX_POROUSMEDIUMFLOW_BOXDFM_FV_ELEMENT_GEOMETRY_HH | ||
18 | #define DUMUX_POROUSMEDIUMFLOW_BOXDFM_FV_ELEMENT_GEOMETRY_HH | ||
19 | |||
20 | #include <optional> | ||
21 | #include <utility> | ||
22 | |||
23 | #include <dune/geometry/type.hh> | ||
24 | #include <dune/localfunctions/lagrange/pqkfactory.hh> | ||
25 | |||
26 | #include <dumux/discretization/scvandscvfiterators.hh> | ||
27 | #include "geometryhelper.hh" | ||
28 | |||
29 | namespace Dumux { | ||
30 | |||
31 | /*! | ||
32 | * \ingroup BoxDFMModel | ||
33 | * \brief Base class for the finite volume geometry vector for box discrete fracture model. | ||
34 | * | ||
35 | * This builds up the sub control volumes and sub control volume faces for each element. | ||
36 | * | ||
37 | * \tparam GG the finite volume grid geometry type | ||
38 | * \tparam enableGridGeometryCache if the grid geometry is cached or not | ||
39 | */ | ||
40 | template<class GG, bool enableGridGeometryCache> | ||
41 | class BoxDfmFVElementGeometry; | ||
42 | |||
43 | //! Specialization in case the FVElementGeometries are stored | ||
44 | template<class GG> | ||
45 | class BoxDfmFVElementGeometry<GG, true> | ||
46 | { | ||
47 | using GridView = typename GG::GridView; | ||
48 | static constexpr int dim = GridView::dimension; | ||
49 | static constexpr int dimWorld = GridView::dimensionworld; | ||
50 | using GridIndexType = typename GridView::IndexSet::IndexType; | ||
51 | using CoordScalar = typename GridView::ctype; | ||
52 | using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType; | ||
53 | public: | ||
54 | //! export type of the element | ||
55 | using Element = typename GridView::template Codim<0>::Entity; | ||
56 | //! Export type of subcontrol volume | ||
57 | using SubControlVolume = typename GG::SubControlVolume; | ||
58 | //! Export type of subcontrol volume face | ||
59 | using SubControlVolumeFace = typename GG::SubControlVolumeFace; | ||
60 | //! Export type of finite volume grid geometry | ||
61 | using GridGeometry = GG; | ||
62 | |||
63 | //! The maximum number of scvs per element (2^dim for cubes) | ||
64 | //! multiplied by 3 for the maximum number of fracture scvs per vertex | ||
65 | static constexpr std::size_t maxNumElementScvs = (1<<dim)*3; | ||
66 | |||
67 | //! Constructor | ||
68 | BoxDfmFVElementGeometry(const GridGeometry& gridGeometry) | ||
69 | : gridGeometryPtr_(&gridGeometry) {} | ||
70 | |||
71 | //! Get a sub control volume with a local scv index | ||
72 | const SubControlVolume& scv(std::size_t scvIdx) const | ||
73 | { return gridGeometry().scvs(eIdx_)[scvIdx]; } | ||
74 | |||
75 | //! Get a sub control volume face with a local scvf index | ||
76 | const SubControlVolumeFace& scvf(std::size_t scvfIdx) const | ||
77 | { return gridGeometry().scvfs(eIdx_)[scvfIdx]; } | ||
78 | |||
79 | /*! | ||
80 | * \brief Iterator range for sub control volumes. | ||
81 | * | ||
82 | * Iterates over all scvs of the bound element. | ||
83 | * This is a free function found by means of ADL. | ||
84 | * To iterate over all sub control volumes of this FVElementGeometry use | ||
85 | * for (auto&& scv : scvs(fvGeometry)). | ||
86 | */ | ||
87 | friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator> | ||
88 | scvs(const BoxDfmFVElementGeometry& fvGeometry) | ||
89 | { | ||
90 | const auto& g = fvGeometry.gridGeometry(); | ||
91 | using Iter = typename std::vector<SubControlVolume>::const_iterator; | ||
92 | return Dune::IteratorRange<Iter>(g.scvs(fvGeometry.eIdx_).begin(), g.scvs(fvGeometry.eIdx_).end()); | ||
93 | } | ||
94 | |||
95 | /*! | ||
96 | * \brief Iterator range for sub control volumes faces. | ||
97 | * | ||
98 | * Iterates over all scvfs of the bound element. | ||
99 | * This is a free function found by means of ADL. | ||
100 | * To iterate over all sub control volume faces of this FVElementGeometry use | ||
101 | * for (auto&& scvf : scvfs(fvGeometry)). | ||
102 | */ | ||
103 | friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator> | ||
104 | scvfs(const BoxDfmFVElementGeometry& fvGeometry) | ||
105 | { | ||
106 | const auto& g = fvGeometry.gridGeometry(); | ||
107 | using Iter = typename std::vector<SubControlVolumeFace>::const_iterator; | ||
108 | return Dune::IteratorRange<Iter>(g.scvfs(fvGeometry.eIdx_).begin(), g.scvfs(fvGeometry.eIdx_).end()); | ||
109 | } | ||
110 | |||
111 | //! Get a local finite element basis | ||
112 | const FeLocalBasis& feLocalBasis() const | ||
113 | { return gridGeometry().feCache().get(element_->type()).localBasis(); } | ||
114 | |||
115 | //! The total number of sub control volumes | ||
116 | std::size_t numScv() const | ||
117 | { return gridGeometry().scvs(eIdx_).size(); } | ||
118 | |||
119 | //! The total number of sub control volume faces | ||
120 | std::size_t numScvf() const | ||
121 | { return gridGeometry().scvfs(eIdx_).size(); } | ||
122 | |||
123 | /*! | ||
124 | * \brief bind the local view (r-value overload) | ||
125 | * This overload is called when an instance of this class is a temporary in the usage context | ||
126 | * This allows a usage like this: `const auto view = localView(...).bind(element);` | ||
127 | */ | ||
128 | BoxDfmFVElementGeometry bind(const Element& element) && | ||
129 | { | ||
130 | this->bindElement(element); | ||
131 | return std::move(*this); | ||
132 | } | ||
133 | |||
134 | //! This function is for compatibility reasons with cc methods | ||
135 | //! The box stencil is always element-local so bind and bindElement are identical. | ||
136 | void bind(const Element& element) & | ||
137 | { this->bindElement(element); } | ||
138 | |||
139 | /*! | ||
140 | * \brief bind the local view (r-value overload) | ||
141 | * This overload is called when an instance of this class is a temporary in the usage context | ||
142 | * This allows a usage like this: `const auto view = localView(...).bindElement(element);` | ||
143 | */ | ||
144 | BoxDfmFVElementGeometry bindElement(const Element& element) && | ||
145 | { | ||
146 | this->bindElement(element); | ||
147 | return std::move(*this); | ||
148 | } | ||
149 | |||
150 | /*! | ||
151 | * \brief Binding of an element, has to be called before using the fvgeometries | ||
152 | * | ||
153 | * Prepares all the volume variables within the element. | ||
154 | * For compatibility reasons with the FVGeometry cache being disabled. | ||
155 | */ | ||
156 | void bindElement(const Element& element) & | ||
157 | { | ||
158 | element_ = element; | ||
159 | eIdx_ = gridGeometry().elementMapper().index(element); | ||
160 | } | ||
161 | |||
162 | //! The global finite volume geometry we are a restriction of | ||
163 | const GridGeometry& gridGeometry() const | ||
164 | { return *gridGeometryPtr_; } | ||
165 | |||
166 | //! Returns true if bind/bindElement has already been called | ||
167 | bool isBound() const | ||
168 | { return static_cast<bool>(element_); } | ||
169 | |||
170 | //! The bound element | ||
171 | const Element& element() const | ||
172 | { return *element_; } | ||
173 | |||
174 | // suppress warnings due to current implementation | ||
175 | // these interfaces should be used! | ||
176 | #pragma GCC diagnostic push | ||
177 | #pragma GCC diagnostic ignored "-Wdeprecated-declarations" | ||
178 | |||
179 | //! Create the geometry of a given sub control volume | ||
180 | typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const | ||
181 | { return scv.geometry(); } | ||
182 | |||
183 | //! Create the geometry of a given sub control volume face | ||
184 | typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const | ||
185 | { return scvf.geometry(); } | ||
186 | |||
187 | #pragma GCC diagnostic pop | ||
188 | |||
189 | private: | ||
190 | const GridGeometry* gridGeometryPtr_; | ||
191 | |||
192 | std::optional<Element> element_; | ||
193 | GridIndexType eIdx_; | ||
194 | }; | ||
195 | |||
196 | //! Specialization in case the FVElementGeometries are not stored | ||
197 | template<class GG> | ||
198 | class BoxDfmFVElementGeometry<GG, false> | ||
199 | { | ||
200 | using GridView = typename GG::GridView; | ||
201 | static constexpr int dim = GridView::dimension; | ||
202 | static constexpr int dimWorld = GridView::dimensionworld; | ||
203 | |||
204 | using GridIndexType = typename GridView::IndexSet::IndexType; | ||
205 | |||
206 | using CoordScalar = typename GridView::ctype; | ||
207 | using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType; | ||
208 | using GeometryHelper = typename GG::GeometryHelper; | ||
209 | public: | ||
210 | //! export type of the element | ||
211 | using Element = typename GridView::template Codim<0>::Entity; | ||
212 | //! Export type of subcontrol volume | ||
213 | using SubControlVolume = typename GG::SubControlVolume; | ||
214 | //! Export type of subcontrol volume face | ||
215 | using SubControlVolumeFace = typename GG::SubControlVolumeFace; | ||
216 | //! Export type of finite volume grid geometry | ||
217 | using GridGeometry = GG; | ||
218 | //! The maximum number of scvs per element (2^dim for cubes) | ||
219 | //! multiplied by 3 for the maximum number of fracture scvs per vertex | ||
220 | static constexpr std::size_t maxNumElementScvs = (1<<dim)*3; | ||
221 | |||
222 | //! Constructor | ||
223 | BoxDfmFVElementGeometry(const GridGeometry& gridGeometry) | ||
224 |
12/48✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 25 taken 2458684 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 2458684 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 2458684 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 2458684 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 6 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 6 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 6 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 6 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 6 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 6 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 6 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 6 times.
✗ Branch 59 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
|
11143040 | : gridGeometryPtr_(&gridGeometry) {} |
225 | |||
226 | //! Get a sub control volume with a local scv index | ||
227 | const SubControlVolume& scv(std::size_t scvIdx) const | ||
228 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 132628712 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 132628712 times.
✓ Branch 4 taken 132628712 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 132628712 times.
✗ Branch 7 not taken.
|
544414032 | { return scvs_[scvIdx]; } |
229 | |||
230 | //! Get a sub control volume face with a local scvf index | ||
231 | const SubControlVolumeFace& scvf(std::size_t scvfIdx) const | ||
232 | { return scvfs_[scvfIdx]; } | ||
233 | |||
234 | /*! | ||
235 | * \brief Iterator range for sub control volumes. | ||
236 | * | ||
237 | * Iterates over all scvs of the bound element. | ||
238 | * This is a free function found by means of ADL. | ||
239 | * To iterate over all sub control volumes of this FVElementGeometry use | ||
240 | * for (auto&& scv : scvs(fvGeometry)). | ||
241 | */ | ||
242 | friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator> | ||
243 | scvs(const BoxDfmFVElementGeometry& fvGeometry) | ||
244 | { | ||
245 | using Iter = typename std::vector<SubControlVolume>::const_iterator; | ||
246 | 627091314 | return Dune::IteratorRange<Iter>(fvGeometry.scvs_.begin(), fvGeometry.scvs_.end()); | |
247 | } | ||
248 | |||
249 | /*! | ||
250 | * \brief Iterator range for sub control volumes faces. | ||
251 | * | ||
252 | * Iterates over all scvfs of the bound element. | ||
253 | * This is a free function found by means of ADL. | ||
254 | * To iterate over all sub control volume faces of this FVElementGeometry use | ||
255 | * for (auto&& scvf : scvfs(fvGeometry)). | ||
256 | */ | ||
257 | friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator> | ||
258 | scvfs(const BoxDfmFVElementGeometry& fvGeometry) | ||
259 | { | ||
260 | using Iter = typename std::vector<SubControlVolumeFace>::const_iterator; | ||
261 | 63739752 | return Dune::IteratorRange<Iter>(fvGeometry.scvfs_.begin(), fvGeometry.scvfs_.end()); | |
262 | } | ||
263 | |||
264 | //! Get a local finite element basis | ||
265 | 4376172 | const FeLocalBasis& feLocalBasis() const | |
266 |
4/8✓ Branch 1 taken 4376172 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4376172 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4376172 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4376172 times.
✗ Branch 11 not taken.
|
30633204 | { return gridGeometry().feCache().get(element_->type()).localBasis(); } |
267 | |||
268 | //! The total number of sub control volumes | ||
269 | std::size_t numScv() const | ||
270 |
2/4✓ Branch 3 taken 8752344 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8752344 times.
✗ Branch 7 not taken.
|
65242288 | { return scvs_.size(); } |
271 | |||
272 | //! The total number of sub control volume faces | ||
273 | std::size_t numScvf() const | ||
274 | 2458684 | { return scvfs_.size(); } | |
275 | |||
276 | /*! | ||
277 | * \brief bind the local view (r-value overload) | ||
278 | * This overload is called when an instance of this class is a temporary in the usage context | ||
279 | * This allows a usage like this: `const auto view = localView(...).bind(element);` | ||
280 | */ | ||
281 | BoxDfmFVElementGeometry bind(const Element& element) && | ||
282 | { | ||
283 | this->bindElement(element); | ||
284 | return std::move(*this); | ||
285 | } | ||
286 | |||
287 | /*! | ||
288 | * \brief Binding of an element, has to be called before using the fvgeometries | ||
289 | * Prepares all the volume variables within the element. | ||
290 | * \note For the box scheme, bind() and bindElement() are identical, but the | ||
291 | * distinction is here for the sake of compatibility with cc schemes. | ||
292 | */ | ||
293 | void bind(const Element& element) & | ||
294 |
1/10✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 8 taken 8852 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
2467536 | { this->bindElement(element); } |
295 | |||
296 | /*! | ||
297 | * \brief bind the local view (r-value overload) | ||
298 | * This overload is called when an instance of this class is a temporary in the usage context | ||
299 | * This allows a usage like this: `const auto view = localView(...).bindElement(element);` | ||
300 | */ | ||
301 | BoxDfmFVElementGeometry bindElement(const Element& element) && | ||
302 | { | ||
303 | this->bindElement(element); | ||
304 | return std::move(*this); | ||
305 | } | ||
306 | |||
307 | /*! | ||
308 | * \brief Binding of an element, has to be called before using the fvgeometries | ||
309 | * Prepares all the volume variables within the element. | ||
310 | */ | ||
311 | 2803452 | void bindElement(const Element& element) & | |
312 | { | ||
313 |
2/2✓ Branch 0 taken 8846 times.
✓ Branch 1 taken 1392880 times.
|
2803452 | element_ = element; |
314 | 5606904 | eIdx_ = gridGeometry().elementMapper().index(element); | |
315 | 2803452 | makeElementGeometries_(); | |
316 | 2803452 | } | |
317 | |||
318 | //! The global finite volume geometry we are a restriction of | ||
319 | ✗ | const GridGeometry& gridGeometry() const | |
320 | ✗ | { return *gridGeometryPtr_; } | |
321 | |||
322 | //! Returns true if bind/bindElement has already been called | ||
323 | bool isBound() const | ||
324 | { return static_cast<bool>(element_); } | ||
325 | |||
326 | //! The bound element | ||
327 | const Element& element() const | ||
328 | 32658432 | { return *element_; } | |
329 | |||
330 | // suppress warnings due to current implementation | ||
331 | // these interfaces should be used! | ||
332 | #pragma GCC diagnostic push | ||
333 | #pragma GCC diagnostic ignored "-Wdeprecated-declarations" | ||
334 | |||
335 | //! Create the geometry of a given sub control volume | ||
336 | typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const | ||
337 | { return scv.geometry(); } | ||
338 | |||
339 | //! Create the geometry of a given sub control volume face | ||
340 | typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const | ||
341 | { return scvf.geometry(); } | ||
342 | |||
343 | #pragma GCC diagnostic pop | ||
344 | |||
345 | private: | ||
346 | |||
347 | 2803452 | void makeElementGeometries_() | |
348 | { | ||
349 | // get the element geometry | ||
350 | 2803452 | const auto& element = *element_; | |
351 | 2803452 | const auto elementGeometry = element.geometry(); | |
352 |
1/2✓ Branch 1 taken 1401726 times.
✗ Branch 2 not taken.
|
2803452 | const auto refElement = referenceElement(element); |
353 | |||
354 | // get the sub control volume geometries of this element | ||
355 |
1/2✓ Branch 1 taken 1401726 times.
✗ Branch 2 not taken.
|
2803452 | GeometryHelper geometryHelper(elementGeometry); |
356 | |||
357 | // construct the sub control volumes | ||
358 |
2/4✓ Branch 1 taken 1401726 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1401726 times.
✗ Branch 5 not taken.
|
5606904 | scvs_.resize(elementGeometry.corners()); |
359 | using LocalIndexType = typename SubControlVolumeFace::Traits::LocalIndexType; | ||
360 |
4/4✓ Branch 0 taken 8761392 times.
✓ Branch 1 taken 2803452 times.
✓ Branch 2 taken 8761392 times.
✓ Branch 3 taken 2803452 times.
|
23129688 | for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx) |
361 | { | ||
362 | // get associated dof index | ||
363 |
2/4✓ Branch 1 taken 4380696 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4380696 times.
✗ Branch 5 not taken.
|
17522784 | const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim); |
364 | |||
365 | // add scv to the local container | ||
366 |
1/2✓ Branch 1 taken 4380696 times.
✗ Branch 2 not taken.
|
8761392 | scvs_[scvLocalIdx] = SubControlVolume(geometryHelper, |
367 | scvLocalIdx, | ||
368 | eIdx_, | ||
369 | dofIdxGlobal); | ||
370 | } | ||
371 | |||
372 | // construct the sub control volume faces | ||
373 |
1/2✓ Branch 1 taken 1401726 times.
✗ Branch 2 not taken.
|
2803452 | const auto numInnerScvf = (dim==1) ? 1 : element.subEntities(dim-1); |
374 |
1/2✓ Branch 1 taken 1401726 times.
✗ Branch 2 not taken.
|
2803452 | scvfs_.resize(numInnerScvf); |
375 | |||
376 | 2803452 | unsigned int scvfLocalIdx = 0; | |
377 |
2/2✓ Branch 0 taken 9309240 times.
✓ Branch 1 taken 2803452 times.
|
12112692 | for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx) |
378 | { | ||
379 | // find the local scv indices this scvf is connected to | ||
380 |
1/4✓ Branch 3 taken 9309240 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
37236960 | std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)), |
381 |
1/2✓ Branch 2 taken 4654620 times.
✗ Branch 3 not taken.
|
9309240 | static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))}); |
382 | |||
383 |
2/4✗ Branch 2 not taken.
✓ Branch 3 taken 9309240 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 9309240 times.
|
18618480 | scvfs_[scvfLocalIdx] = SubControlVolumeFace(geometryHelper, |
384 | element, | ||
385 | elementGeometry, | ||
386 | scvfLocalIdx, | ||
387 |
1/2✓ Branch 1 taken 9309240 times.
✗ Branch 2 not taken.
|
9309240 | std::move(localScvIndices)); |
388 | } | ||
389 | |||
390 | // construct the ... | ||
391 | // ... sub-control volume faces on the domain boundary | ||
392 | // ... sub-control volumes on fracture facets | ||
393 | // ... sub-control volume faces on fracture facets | ||
394 | // NOTE We do not construct fracture scvfs on boundaries here! | ||
395 | // That means specifying Neumann fluxes on only the fractures is not possible | ||
396 | // However, it is difficult to interpret this here as a fracture ending on the boundary | ||
397 | // could also be connected to a facet which is both boundary and fracture at the same time! | ||
398 | // In that case, the fracture boundary scvf wouldn't make sense. In order to do it properly | ||
399 | // we would have to find only those fractures that are at the boundary and aren't connected | ||
400 | // to a fracture which is a boundary. | ||
401 |
1/2✓ Branch 1 taken 1401726 times.
✗ Branch 2 not taken.
|
2803452 | LocalIndexType scvLocalIdx = element.subEntities(dim); |
402 |
7/14✓ Branch 1 taken 1401726 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1401726 times.
✓ Branch 4 taken 1401726 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1401726 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 11564844 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 13 taken 4380696 times.
✓ Branch 14 taken 4380696 times.
✗ Branch 15 not taken.
|
27334866 | for (const auto& intersection : intersections(gridGeometry().gridView(), element)) |
403 | { | ||
404 | // first, obtain all vertex indices on this intersection | ||
405 |
2/4✓ Branch 1 taken 8761392 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8761392 times.
✗ Branch 5 not taken.
|
17522784 | const auto& isGeometry = intersection.geometry(); |
406 |
2/3✓ Branch 0 taken 154224 times.
✓ Branch 1 taken 8607168 times.
✗ Branch 2 not taken.
|
8761392 | const auto numCorners = isGeometry.corners(); |
407 |
2/3✓ Branch 0 taken 154224 times.
✓ Branch 1 taken 8607168 times.
✗ Branch 2 not taken.
|
8761392 | const auto idxInInside = intersection.indexInInside(); |
408 | |||
409 |
4/10✓ Branch 1 taken 8761392 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8761392 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 8761392 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4380696 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
35045568 | std::vector<GridIndexType> isVertexIndices(numCorners); |
410 |
2/2✓ Branch 0 taken 18618480 times.
✓ Branch 1 taken 8761392 times.
|
27379872 | for (unsigned int vIdxLocal = 0; vIdxLocal < numCorners; ++vIdxLocal) |
411 |
1/2✓ Branch 3 taken 18618480 times.
✗ Branch 4 not taken.
|
37236960 | isVertexIndices[vIdxLocal] = gridGeometry().vertexMapper().subIndex(element, |
412 | refElement.subEntity(idxInInside, 1, vIdxLocal, dim), | ||
413 | dim); | ||
414 | |||
415 |
3/3✓ Branch 0 taken 348592 times.
✓ Branch 1 taken 4134950 times.
✓ Branch 2 taken 4277850 times.
|
8761392 | if (intersection.boundary()) |
416 | { | ||
417 |
2/2✓ Branch 0 taken 513444 times.
✓ Branch 1 taken 205692 times.
|
719136 | for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < numCorners; ++isScvfLocalIdx) |
418 | { | ||
419 | // find the scv this scvf is connected to | ||
420 | 513444 | const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, isScvfLocalIdx, dim)); | |
421 |
2/6✓ Branch 1 taken 513444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 513444 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1540332 | std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx}; |
422 | |||
423 | 1026888 | scvfs_.emplace_back(geometryHelper, | |
424 | intersection, | ||
425 | isGeometry, | ||
426 | isScvfLocalIdx, | ||
427 | scvfLocalIdx, | ||
428 |
1/2✓ Branch 1 taken 513444 times.
✗ Branch 2 not taken.
|
513444 | std::move(localScvIndices)); |
429 | |||
430 | // increment local counter | ||
431 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 513444 times.
|
513444 | scvfLocalIdx++; |
432 | } | ||
433 | } | ||
434 | |||
435 | // maybe add fracture scvs & scvfs | ||
436 |
3/4✓ Branch 1 taken 8761392 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 286404 times.
✓ Branch 4 taken 8474988 times.
|
8761392 | if (this->gridGeometry().isOnFracture(element, intersection)) |
437 | { | ||
438 | // add fracture scv for each vertex of intersection | ||
439 |
1/2✓ Branch 1 taken 286404 times.
✗ Branch 2 not taken.
|
286404 | const auto curNumScvs = scvs_.size(); |
440 |
1/2✓ Branch 1 taken 286404 times.
✗ Branch 2 not taken.
|
286404 | scvs_.reserve(curNumScvs+numCorners); |
441 |
2/2✓ Branch 0 taken 600528 times.
✓ Branch 1 taken 286404 times.
|
886932 | for (unsigned int vIdxLocal = 0; vIdxLocal < numCorners; ++vIdxLocal) |
442 | 1801584 | scvs_.emplace_back(geometryHelper, | |
443 | intersection, | ||
444 | isGeometry, | ||
445 | vIdxLocal, | ||
446 |
1/4✓ Branch 1 taken 600528 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
600528 | static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, vIdxLocal, dim)), |
447 | scvLocalIdx++, | ||
448 | idxInInside, | ||
449 | 600528 | eIdx_, | |
450 | 1201056 | isVertexIndices[vIdxLocal]); | |
451 | |||
452 | // add fracture scvf for each edge of the intersection in 3d | ||
453 | if (dim == 3) | ||
454 | { | ||
455 |
1/2✓ Branch 1 taken 27720 times.
✗ Branch 2 not taken.
|
27720 | const auto& faceRefElement = referenceElement(isGeometry); |
456 |
4/4✓ Branch 0 taken 83160 times.
✓ Branch 1 taken 27720 times.
✓ Branch 2 taken 83160 times.
✓ Branch 3 taken 27720 times.
|
194040 | for (unsigned int edgeIdx = 0; edgeIdx < faceRefElement.size(1); ++edgeIdx) |
457 | { | ||
458 | // inside/outside scv indices in face local node numbering | ||
459 | 332640 | std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(faceRefElement.subEntity(edgeIdx, 1, 0, dim-1)), | |
460 |
1/2✓ Branch 2 taken 83160 times.
✗ Branch 3 not taken.
|
83160 | static_cast<LocalIndexType>(faceRefElement.subEntity(edgeIdx, 1, 1, dim-1))}); |
461 | |||
462 | // add offset to get the right scv indices | ||
463 | 249480 | std::for_each( localScvIndices.begin(), | |
464 | localScvIndices.end(), | ||
465 | 166320 | [curNumScvs] (auto& elemLocalIdx) { elemLocalIdx += curNumScvs; } ); | |
466 | |||
467 | // add scvf | ||
468 |
4/6✓ Branch 0 taken 7182 times.
✓ Branch 1 taken 34398 times.
✓ Branch 2 taken 41580 times.
✓ Branch 3 taken 41580 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
166320 | scvfs_.emplace_back(geometryHelper, |
469 | intersection, | ||
470 | isGeometry, | ||
471 | edgeIdx, | ||
472 | scvfLocalIdx++, | ||
473 |
1/2✓ Branch 1 taken 83160 times.
✗ Branch 2 not taken.
|
83160 | std::move(localScvIndices), |
474 | intersection.boundary()); | ||
475 | } | ||
476 | } | ||
477 | |||
478 | // dim == 2, intersection is an edge, make 1 scvf | ||
479 | else | ||
480 | { | ||
481 | // inside/outside scv indices in face local node numbering | ||
482 |
2/6✓ Branch 1 taken 258684 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 258684 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
776052 | std::vector<LocalIndexType> localScvIndices({0, 1}); |
483 | |||
484 | // add offset such that the fracture scvs above are addressed | ||
485 | 776052 | std::for_each( localScvIndices.begin(), | |
486 | localScvIndices.end(), | ||
487 | 517368 | [curNumScvs] (auto& elemLocalIdx) { elemLocalIdx += curNumScvs; } ); | |
488 | |||
489 | // add scvf | ||
490 |
2/2✓ Branch 0 taken 2146 times.
✓ Branch 1 taken 127196 times.
|
258684 | scvfs_.emplace_back(geometryHelper, |
491 | intersection, | ||
492 | isGeometry, | ||
493 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 258684 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
258684 | /*idxOnIntersection*/0, |
494 | scvfLocalIdx++, | ||
495 |
1/2✓ Branch 1 taken 258684 times.
✗ Branch 2 not taken.
|
258684 | std::move(localScvIndices), |
496 | intersection.boundary()); | ||
497 | } | ||
498 | } | ||
499 | } | ||
500 | 2803452 | } | |
501 | |||
502 | //! The bound element | ||
503 | std::optional<Element> element_; | ||
504 | GridIndexType eIdx_; | ||
505 | |||
506 | //! The global geometry this is a restriction of | ||
507 | const GridGeometry* gridGeometryPtr_; | ||
508 | |||
509 | //! The element-local scvs & scvfs | ||
510 | std::vector<SubControlVolume> scvs_; | ||
511 | std::vector<SubControlVolumeFace> scvfs_; | ||
512 | }; | ||
513 | |||
514 | } // end namespace Dumux | ||
515 | |||
516 | #endif | ||
517 |