3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
discretization/cellcentered/tpfa/fvgridgeometry.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 *****************************************************************************/
26#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_GRID_GEOMETRY_HH
27#define DUMUX_DISCRETIZATION_CCTPFA_FV_GRID_GEOMETRY_HH
28
29#include <algorithm>
30
33
41
42namespace Dumux {
43
50template<class GridView, class MapperTraits = DefaultMapperTraits<GridView>>
52: public MapperTraits
53{
56
57 template<class GridGeometry>
59
60 template<class GridGeometry, bool enableCache>
62
67 static constexpr int maxNumScvfNeighbors = int(GridView::dimension)<int(GridView::dimensionworld) ? 8 : 1<<(GridView::dimension-1);
68};
69
76template<class GridView,
77 bool enableGridGeometryCache = false,
80
87template<class GV, class Traits>
88class CCTpfaFVGridGeometry<GV, true, Traits>
89: public BaseGridGeometry<GV, Traits>
90{
93 using ConnectivityMap = typename Traits::template ConnectivityMap<ThisType>;
94 using GridIndexType = typename IndexTraits<GV>::GridIndex;
95 using Element = typename GV::template Codim<0>::Entity;
96
97 static const int dim = GV::dimension;
98 static const int dimWorld = GV::dimensionworld;
99
100public:
102 using LocalView = typename Traits::template LocalView<ThisType, true>;
104 using SubControlVolume = typename Traits::SubControlVolume;
106 using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
108 using DofMapper = typename Traits::ElementMapper;
109
112
114 static constexpr int maxElementStencilSize = LocalView::maxNumElementScvfs*Traits::maxNumScvfNeighbors + 1;
115
117 using GridView = GV;
118
121 : ParentType(gridView)
122 {
123 // Check if the overlap size is what we expect
125 DUNE_THROW(Dune::InvalidStateException, "The cctpfa discretization method needs at least an overlap of 1 for parallel computations. "
126 << " Set the parameter \"Grid.Overlap\" in the input file.");
127 }
128
131 const DofMapper& dofMapper() const
132 { return this->elementMapper(); }
133
135 std::size_t numScv() const
136 {
137 return scvs_.size();
138 }
139
141 std::size_t numScvf() const
142 {
143 return scvfs_.size();
144 }
145
147 std::size_t numBoundaryScvf() const
148 {
149 return numBoundaryScvf_;
150 }
151
153 std::size_t numDofs() const
154 { return this->gridView().size(0); }
155
157 void update()
158 {
159 ParentType::update();
160
161 // clear containers (necessary after grid refinement)
162 scvs_.clear();
163 scvfs_.clear();
164 scvfIndicesOfScv_.clear();
165 flipScvfIndices_.clear();
166
167 // determine size of containers
168 std::size_t numScvs = numDofs();
169 std::size_t numScvf = 0;
170 for (const auto& element : elements(this->gridView()))
171 numScvf += element.subEntities(1);
172
173 // reserve memory
174 scvs_.resize(numScvs);
175 scvfs_.reserve(numScvf);
176 scvfIndicesOfScv_.resize(numScvs);
177 hasBoundaryScvf_.resize(numScvs, false);
178
179 // Build the scvs and scv faces
180 GridIndexType scvfIdx = 0;
181 numBoundaryScvf_ = 0;
182 for (const auto& element : elements(this->gridView()))
183 {
184 const auto eIdx = this->elementMapper().index(element);
185 scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
186
187 // the element-wise index sets for finite volume geometry
188 std::vector<GridIndexType> scvfsIndexSet;
189 scvfsIndexSet.reserve(element.subEntities(1));
190
191 // for network grids there might be multiple intersection with the same geometryInInside
192 // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
193 using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage;
194 std::vector<ScvfGridIndexStorage> outsideIndices;
195 if (dim < dimWorld)
196 {
198 outsideIndices.resize(element.subEntities(1));
199 std::for_each(outsideIndices.begin(), outsideIndices.end(), [eIdx] (auto& nIndices) { nIndices.push_back(eIdx); });
200
201 // second, insert neighbors
202 for (const auto& intersection : intersections(this->gridView(), element))
203 {
204 if (intersection.neighbor())
205 {
206 const auto nIdx = this->elementMapper().index( intersection.outside() );
207 outsideIndices[intersection.indexInInside()].push_back(nIdx);
208 }
209 }
210 }
211
212 for (const auto& intersection : intersections(this->gridView(), element))
213 {
214 // inner sub control volume faces (includes periodic boundaries)
215 if (intersection.neighbor())
216 {
217 // update the grid geometry if we have periodic boundaries
218 if (intersection.boundary())
219 this->setPeriodic();
220
221 if (dim == dimWorld)
222 {
223 const auto nIdx = this->elementMapper().index(intersection.outside());
224 scvfs_.emplace_back(intersection,
225 intersection.geometry(),
226 scvfIdx,
227 ScvfGridIndexStorage({eIdx, nIdx}),
228 false);
229 scvfsIndexSet.push_back(scvfIdx++);
230 }
231 // this is for network grids
232 // (will be optimized away of dim == dimWorld)
233 else
234 {
235 auto indexInInside = intersection.indexInInside();
236 // check if we already handled this facet
237 if (outsideIndices[indexInInside].empty())
238 continue;
239 else
240 {
241 scvfs_.emplace_back(intersection,
242 intersection.geometry(),
243 scvfIdx,
244 outsideIndices[indexInInside],
245 false);
246 scvfsIndexSet.push_back(scvfIdx++);
247 outsideIndices[indexInInside].clear();
248 }
249 }
250 }
251 // boundary sub control volume faces
252 else if (intersection.boundary())
253 {
254 scvfs_.emplace_back(intersection,
255 intersection.geometry(),
256 scvfIdx,
257 ScvfGridIndexStorage({eIdx, static_cast<GridIndexType>(this->gridView().size(0) + numBoundaryScvf_++)}),
258 true);
259 scvfsIndexSet.push_back(scvfIdx++);
260
261 hasBoundaryScvf_[eIdx] = true;
262 }
263 }
264
265 // Save the scvf indices belonging to this scv to build up fv element geometries fast
266 scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
267 }
268
269 // Make the flip index set for network, surface, and periodic grids
270 if (dim < dimWorld || this->isPeriodic())
271 {
272 flipScvfIndices_.resize(scvfs_.size());
273 for (auto&& scvf : scvfs_)
274 {
275 if (scvf.boundary())
276 continue;
277
278 flipScvfIndices_[scvf.index()].resize(scvf.numOutsideScvs());
279 const auto insideScvIdx = scvf.insideScvIdx();
280 // check which outside scvf has the insideScvIdx index in its outsideScvIndices
281 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
282 flipScvfIndices_[scvf.index()][i] = findFlippedScvfIndex_(insideScvIdx, scvf.outsideScvIdx(i));
283 }
284 }
285
286 // build the connectivity map for an effecient assembly
287 connectivityMap_.update(*this);
288 }
289
291 const SubControlVolume& scv(GridIndexType scvIdx) const
292 {
293 return scvs_[scvIdx];
294 }
295
297 const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const
298 {
299 return scvfs_[scvfIdx];
300 }
301
304 const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
305 {
306 return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]];
307 }
308
310 const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
311 {
312 return scvfIndicesOfScv_[scvIdx];
313 }
314
319 const ConnectivityMap &connectivityMap() const
320 { return connectivityMap_; }
321
323 bool hasBoundaryScvf(GridIndexType eIdx) const
324 { return hasBoundaryScvf_[eIdx]; }
325
326private:
327 // find the scvf that has insideScvIdx in its outsideScvIdx list and outsideScvIdx as its insideScvIdx
328 GridIndexType findFlippedScvfIndex_(GridIndexType insideScvIdx, GridIndexType outsideScvIdx)
329 {
330 // go over all potential scvfs of the outside scv
331 for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
332 {
333 const auto& outsideScvf = this->scvf(outsideScvfIndex);
334 for (unsigned int j = 0; j < outsideScvf.numOutsideScvs(); ++j)
335 if (outsideScvf.outsideScvIdx(j) == insideScvIdx)
336 return outsideScvf.index();
337 }
338
339 DUNE_THROW(Dune::InvalidStateException, "No flipped version of this scvf found!");
340 }
341
343 ConnectivityMap connectivityMap_;
344
346 std::vector<SubControlVolume> scvs_;
347 std::vector<SubControlVolumeFace> scvfs_;
348 std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
349 std::size_t numBoundaryScvf_;
350 std::vector<bool> hasBoundaryScvf_;
351
353 std::vector<std::vector<GridIndexType>> flipScvfIndices_;
354};
355
363template<class GV, class Traits>
364class CCTpfaFVGridGeometry<GV, false, Traits>
365: public BaseGridGeometry<GV, Traits>
366{
369 using ConnectivityMap = typename Traits::template ConnectivityMap<ThisType>;
370
371 using GridIndexType = typename IndexTraits<GV>::GridIndex;
372 using Element = typename GV::template Codim<0>::Entity;
373
374 static const int dim = GV::dimension;
375 static const int dimWorld = GV::dimensionworld;
376
377 using ScvfGridIndexStorage = typename Traits::SubControlVolumeFace::Traits::GridIndexStorage;
378 using NeighborVolVarIndices = typename std::conditional_t< (dim<dimWorld),
379 ScvfGridIndexStorage,
380 Dune::ReservedVector<GridIndexType, 1> >;
381
382public:
384 using LocalView = typename Traits::template LocalView<ThisType, false>;
386 using SubControlVolume = typename Traits::SubControlVolume;
388 using SubControlVolumeFace = typename Traits::SubControlVolumeFace;
390 using DofMapper = typename Traits::ElementMapper;
391
394
396 static constexpr int maxElementStencilSize = LocalView::maxNumElementScvfs*Traits::maxNumScvfNeighbors + 1;
397
399 using GridView = GV;
400
403 : ParentType(gridView)
404 {
405 // Check if the overlap size is what we expect
407 DUNE_THROW(Dune::InvalidStateException, "The cctpfa discretization method needs at least an overlap of 1 for parallel computations. "
408 << " Set the parameter \"Grid.Overlap\" in the input file.");
409 }
410
413 const DofMapper& dofMapper() const
414 { return this->elementMapper(); }
415
417 std::size_t numScv() const
418 {
419 return numScvs_;
420 }
421
423 std::size_t numScvf() const
424 {
425 return numScvf_;
426 }
427
429 std::size_t numBoundaryScvf() const
430 {
431 return numBoundaryScvf_;
432 }
433
435 std::size_t numDofs() const
436 { return this->gridView().size(0); }
437
439 void update()
440 {
441 ParentType::update();
442
443 // clear local data
444 scvfIndicesOfScv_.clear();
445 neighborVolVarIndices_.clear();
446
447 // reserve memory or resize the containers
448 numScvs_ = numDofs();
449 numScvf_ = 0;
450 numBoundaryScvf_ = 0;
451 scvfIndicesOfScv_.resize(numScvs_);
452 neighborVolVarIndices_.resize(numScvs_);
453
454 // Build the SCV and SCV face
455 for (const auto& element : elements(this->gridView()))
456 {
457 const auto eIdx = this->elementMapper().index(element);
458
459 // the element-wise index sets for finite volume geometry
460 auto numLocalFaces = element.subEntities(1);
461 std::vector<GridIndexType> scvfsIndexSet;
462 std::vector<NeighborVolVarIndices> neighborVolVarIndexSet;
463 scvfsIndexSet.reserve(numLocalFaces);
464 neighborVolVarIndexSet.reserve(numLocalFaces);
465
466 // for network grids there might be multiple intersection with the same geometryInInside
467 // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
468 std::vector<NeighborVolVarIndices> outsideIndices;
469 if (dim < dimWorld)
470 {
471 outsideIndices.resize(numLocalFaces);
472 for (const auto& intersection : intersections(this->gridView(), element))
473 {
474 if (intersection.neighbor())
475 {
476 const auto nIdx = this->elementMapper().index(intersection.outside());
477 outsideIndices[intersection.indexInInside()].push_back(nIdx);
478 }
479 }
480 }
481
482 for (const auto& intersection : intersections(this->gridView(), element))
483 {
484 // inner sub control volume faces (includes periodic boundaries)
485 if (intersection.neighbor())
486 {
487 // update the grid geometry if we have periodic boundaries
488 if (intersection.boundary())
489 this->setPeriodic();
490
491 if (dim == dimWorld)
492 {
493 scvfsIndexSet.push_back(numScvf_++);
494 const auto nIdx = this->elementMapper().index(intersection.outside());
495 neighborVolVarIndexSet.emplace_back(NeighborVolVarIndices({nIdx}));
496 }
497 // this is for network grids
498 // (will be optimized away of dim == dimWorld)
499 else
500 {
501 auto indexInInside = intersection.indexInInside();
502 // check if we already handled this facet
503 if (outsideIndices[indexInInside].empty())
504 continue;
505 else
506 {
507 scvfsIndexSet.push_back(numScvf_++);
508 neighborVolVarIndexSet.emplace_back(std::move(outsideIndices[indexInInside]));
509 outsideIndices[indexInInside].clear();
510 }
511 }
512 }
513 // boundary sub control volume faces
514 else if (intersection.boundary())
515 {
516 scvfsIndexSet.push_back(numScvf_++);
517 neighborVolVarIndexSet.emplace_back(NeighborVolVarIndices({static_cast<GridIndexType>(numScvs_ + numBoundaryScvf_++)}));
518 }
519 }
520
521 // store the sets of indices in the data container
522 scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
523 neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
524 }
525
526 // build the connectivity map for an effecient assembly
527 connectivityMap_.update(*this);
528 }
529
530 const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
531 { return scvfIndicesOfScv_[scvIdx]; }
532
534 const std::vector<NeighborVolVarIndices>& neighborVolVarIndices(GridIndexType scvIdx) const
535 { return neighborVolVarIndices_[scvIdx]; }
536
541 const ConnectivityMap &connectivityMap() const
542 { return connectivityMap_; }
543
544private:
545
547 std::size_t numScvs_;
548 std::size_t numScvf_;
549 std::size_t numBoundaryScvf_;
550
552 ConnectivityMap connectivityMap_;
553
555 std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
556 std::vector<std::vector<NeighborVolVarIndices>> neighborVolVarIndices_;
557};
558
559} // end namespace Dumux
560
561#endif
Defines the default element and vertex mapper types.
Defines the index types used for grid and local indices.
Base class for grid geometries.
Check the overlap size for different discretization methods.
The available discretization methods in Dumux.
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Struture to define the index types used for grid and local indices.
Definition: indextraits.hh:38
Base class for all finite volume grid geometries.
Definition: basegridgeometry.hh:49
GV GridView
export the grid view type
Definition: basegridgeometry.hh:64
A simple version of the connectivity map for cellcentered schemes. This implementation works for sche...
Definition: cellcentered/connectivitymap.hh:53
Sub control volumes for cell-centered discretization schemes.
Definition: discretization/cellcentered/subcontrolvolume.hh:62
Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models This builds up th...
Definition: discretization/cellcentered/tpfa/fvelementgeometry.hh:50
The default traits for the tpfa finite volume grid geometry Defines the scv and scvf types and the ma...
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:53
static constexpr int maxNumScvfNeighbors
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:67
The finite volume geometry (scvs and scvfs) for cell-centered TPFA models on a grid view This builds ...
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:79
The finite volume geometry (scvs and scvfs) for cell-centered TPFA models on a grid view This builds ...
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:90
typename Traits::template LocalView< ThisType, true > LocalView
export the type of the fv element geometry (the local view type)
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:102
const ConnectivityMap & connectivityMap() const
Returns the connectivity map of which dofs have derivatives with respect to a given dof.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:319
void update()
update all fvElementGeometries (do this again after grid adaption)
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:157
std::size_t numScv() const
The total number of sub control volumes.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:135
CCTpfaFVGridGeometry(const GridView &gridView)
Constructor.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:120
typename Traits::SubControlVolume SubControlVolume
export the type of sub control volume
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:104
typename Traits::SubControlVolumeFace SubControlVolumeFace
export the type of sub control volume
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:106
std::size_t numDofs() const
The total number of degrees of freedom.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:153
const std::vector< GridIndexType > & scvfIndicesOfScv(GridIndexType scvIdx) const
Get the sub control volume face indices of an scv by global index.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:310
std::size_t numScvf() const
The total number of sub control volume faces.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:141
const SubControlVolumeFace & scvf(GridIndexType scvfIdx) const
Get a sub control volume face with a global scvf index.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:297
bool hasBoundaryScvf(GridIndexType eIdx) const
Returns whether one of the geometry's scvfs lies on a boundary.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:323
std::size_t numBoundaryScvf() const
The total number of boundary sub control volume faces.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:147
const SubControlVolumeFace & flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx=0) const
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:304
typename Traits::ElementMapper DofMapper
export dof mapper type
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:108
const SubControlVolume & scv(GridIndexType scvIdx) const
Get a sub control volume with a global scv index.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:291
const DofMapper & dofMapper() const
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:131
The finite volume geometry (scvs and scvfs) for cell-centered TPFA models on a grid view This builds ...
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:366
typename Traits::template LocalView< ThisType, false > LocalView
export the type of the fv element geometry (the local view type)
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:384
const std::vector< GridIndexType > & scvfIndicesOfScv(GridIndexType scvIdx) const
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:530
typename Traits::ElementMapper DofMapper
export dof mapper type
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:390
const DofMapper & dofMapper() const
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:413
typename Traits::SubControlVolumeFace SubControlVolumeFace
export the type of sub control volume
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:388
typename Traits::SubControlVolume SubControlVolume
export the type of sub control volume
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:386
std::size_t numScvf() const
The total number of sub control volume faces.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:423
const ConnectivityMap & connectivityMap() const
Returns the connectivity map of which dofs have derivatives with respect to a given dof.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:541
void update()
update all fvElementGeometries (do this again after grid adaption)
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:439
CCTpfaFVGridGeometry(const GridView &gridView)
Constructor.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:402
std::size_t numDofs() const
The total number of degrees of freedom.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:435
std::size_t numBoundaryScvf() const
The total number of boundary sub control volume faces.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:429
const std::vector< NeighborVolVarIndices > & neighborVolVarIndices(GridIndexType scvIdx) const
Return the neighbor volVar indices for all scvfs in the scv with index scvIdx.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:534
std::size_t numScv() const
The total number of sub control volumes.
Definition: discretization/cellcentered/tpfa/fvgridgeometry.hh:417
The sub control volume face.
Definition: discretization/cellcentered/tpfa/subcontrolvolumeface.hh:90
Check if the overlap size is valid for a given discretization method.
Definition: checkoverlapsize.hh:40
Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models This builds up th...
Sub control volumes for cell-centered discretization schemes.
The sub control volume face.
Stores the face indices corresponding to the neighbors of an element that contribute to the derivativ...