My Project
Loading...
Searching...
No Matches
Indexsets.hpp
1//===========================================================================
2//
3// File: Indexsets.hpp
4//
5// Created: Fri May 29 23:30:01 2009
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// Bård Skaflestad <bard.skaflestad@sintef.no>
9//
10// $Date$
11//
12// $Revision$
13//
14//===========================================================================
15
16/*
17Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18Copyright 2009, 2010, 2022 Equinor ASA.
19
20This file is part of The Open Porous Media project (OPM).
21
22OPM is free software: you can redistribute it and/or modify
23it under the terms of the GNU General Public License as published by
24the Free Software Foundation, either version 3 of the License, or
25(at your option) any later version.
26
27OPM is distributed in the hope that it will be useful,
28but WITHOUT ANY WARRANTY; without even the implied warranty of
29MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30GNU General Public License for more details.
31
32You should have received a copy of the GNU General Public License
33along with OPM. If not, see <http://www.gnu.org/licenses/>.
34*/
35
36#ifndef OPM_INDEXSETS_HEADER
37#define OPM_INDEXSETS_HEADER
38
39#include <dune/geometry/type.hh>
40#include <opm/common/ErrorMacros.hpp>
41#include "GlobalIdMapping.hpp"
42#include "Intersection.hpp"
43
44#include <unordered_map>
45namespace Dune
46{
47 namespace cpgrid
48 {
53 {
54 public:
57 typedef int IndexType;
58
60 template <int cc>
61 struct Codim
62 {
64 };
65
68 typedef std::vector<GeometryType> Types;
69
73 IndexSet() : IndexSet(0,0){}
74
75 IndexSet(std::size_t numCells, std::size_t numPoints)
76 {
77 geom_types_[0].emplace_back(Dune::GeometryTypes::cube(3));
78 geom_types_[3].emplace_back(Dune::GeometryTypes::cube(0));
79 size_codim_map_[0] = numCells;
80 size_codim_map_[3] = numPoints;
81 }
82
85 {}
86
91 const Types& geomTypes(int codim) const
92 {
93 return geom_types_[codim];
94 }
95
100 const Types& types(int codim) const
101 {
102 return geom_types_[codim];
103 }
104
109 int size(GeometryType type) const
110 {
111 if (type.isCube()) {
112 return size(3 - type.dim()); // return grid_.size(type);
113 } else {
114 return 0;
115 }
116 }
117
118
123 int size(int codim) const
124 {
125 return size_codim_map_[codim]; //grid_.size(codim)
126 }
127
128
134 template<int cd>
136 {
137 return e.index();
138 }
139
145 template<class EntityType>
146 IndexType index(const EntityType& e) const
147 {
148 return e.index();
149 }
150
156 template <int cc>
157 IndexType subIndex(const cpgrid::Entity<0>& e, int i) const
158 {
159 return index(e.template subEntity<cc>(i));
160 }
161
167 IndexType subIndex(const cpgrid::Entity<0>& e, int i, unsigned int cc) const;
168
169
170 template<int codim>
171 IndexType subIndex(const cpgrid::Entity<codim>& /* e */, int /* i */, unsigned int /* cc */) const
172 {
173 DUNE_THROW(NotImplemented, "subIndex not implemented for codim"
174 << codim << "entities.");
175 }
181 template <class EntityType>
182 bool contains(const EntityType& e) const
183 {
184 // return index(e) >= 0 && index(e) < grid_.size(EntityType::codimension); //EntityType::codimension == 0;
185 return index(e) >= 0 && index(e) < this->size(EntityType::codimension);
186 }
187
188 private:
189 // const CpGridData& grid_;
190 Types geom_types_[4];
191 std::array<int,4> size_codim_map_{0,0,0,0};
192 };
193
194
195 class IdSet
196 {
197 friend class ReversePointGlobalIdSet;
198 friend class Dune::cpgrid::CpGridData;
199 //friend class Dune::cpgrid::LevelGlobalIdSet; Not needed due to repeated code in LevelGlobalIdSet (computeId_cell and computeId_point)
200 public:
201 typedef int IdType;
202
203 IdSet(const CpGridData& grid)
204 : grid_(grid)
205 {
206 }
207
208 IdType id(const cpgrid::Entity<0>& e) const
209 {
210 return computeId_cell(e);
211 }
212
213 IdType id(const cpgrid::Entity<3>& e) const
214 {
215 return computeId_point(e);
216 }
217
218 template<class EntityType>
219 IdType id(const EntityType& e) const
220 {
221 return computeId(e);
222 }
223
225 IdType id( const cpgrid::Intersection& intersection ) const
226 {
227 return intersection.id();
228 }
229
230 template<int cc>
231 IdType subId(const cpgrid::Entity<0>& e, int i) const
232 {
233 return id(e.template subEntity<cc>(i));
234 }
235
236 IdType subId(const cpgrid::Entity<0>& e, int i, int cc) const;
237
238 private:
239
240 template<class EntityType>
241 IdType computeId(const EntityType& e) const
242 {
243 IdType myId = 0;
244 for( int c=0; c<EntityType::codimension; ++c )
245 myId += grid_.indexSet().size( c );
246 return myId + e.index();
247 }
248
249 const CpGridData& grid_;
250
251 IdType computeId_cell(const cpgrid::Entity<0>& e) const
252 {
253 IdType myId = 0;
254 // Case: Leaf grid view is a mixed of coarse and fined cells.
255 if (grid_.levelData().size() > 1) {
256 const auto& gridIdx = grid_.getGridIdx();
257 // Level zero grid
258 if ( gridIdx == 0 ) {
259 return myId + e.index();
260 }
261 // Level 1, 2, ...., maxLevel refined grids
262 if ( (gridIdx>0) && (gridIdx < static_cast<int>(grid_.levelData().size() -1)) ) {
263 if ((e.level() != gridIdx)) { // cells equiv to pre-existing cells
264 return grid_.levelData()[e.level()]->localIdSet().id(e.getEquivLevelElem());
265 }
266 else {
267 // Count (and add to myId) all the entities of all the codimensions (for CpGrid, only 0 and 3)
268 // from all the "previous" level grids.
269 for (int lowerLevel = 0; lowerLevel< gridIdx; ++lowerLevel) {
270 for( int c=0; c<4; ++c ) {
271 myId += grid_.levelData()[lowerLevel]->indexSet().size( c );
272 }
273 }
274 return myId + e.index();
275 }
276 }
277 else { // Leaf grid view (grid view with mixed coarse and refined cells).
278 assert( grid_.getGridIdx() == (static_cast<int>(grid_.levelData().size()) -1) );
279 // In this case, we search for the ids defined in previous levels
280 // (since each entities must keep its id along the entire hiearchy)
281 std::array<int,2> level_levelIdx = {0,0};
282 level_levelIdx = grid_.leaf_to_level_cells_[e.index()];
283 const auto& levelEntity = cpgrid::Entity<0>(*(grid_.levelData()[level_levelIdx[0]]), level_levelIdx[1], true);
284 return grid_.levelData()[level_levelIdx[0]]->local_id_set_ ->id(levelEntity);
285 }
286 } // end-if-data_.size()>1
287 else { // Case: No LGRs / No refined level grids. Only level 0 grid (GLOBAL grid).
288 return myId + e.index();
289 }
290 }
291
292 IdType computeId_point(const cpgrid::Entity<3>& e) const
293 {
294 IdType myId = 0;
295 // Case: Leaf grid view is a mixed of coarse and fined cells.
296 if (grid_.levelData().size() > 1) {
297 const auto& gridIdx = grid_.getGridIdx();
298 // Level zero grid
299 if ( gridIdx == 0 ) {
300 // Count all the entities of (all the levels) level 0 of all codimensions lower than 3 (for CpGrid, only codim = 0 cells).
301 for( int c=0; c<3; ++c ) {
302 myId += grid_.indexSet().size( c );
303 }
304 return myId + e.index();
305 }
306 // Level 1, 2, ...., maxLevel refined grids.
307 if ( (gridIdx>0) && (gridIdx < static_cast<int>(grid_.levelData().size() -1)) ) {
308 const auto& level_levelIdx = grid_.corner_history_[e.index()];
309 if(level_levelIdx[0] != -1) { // corner equiv to a pre-exisiting level corner
310 const auto& levelEntity = cpgrid::Entity<3>(*(grid_.levelData()[level_levelIdx[0]]), level_levelIdx[1], true);
311 return grid_.levelData()[level_levelIdx[0]]->localIdSet().id(levelEntity);
312 }
313 else {
314 // Count (and add to myId) all the entities of all the codimensions (for CpGrid, only 0 and 3)
315 // from all the "previous" level grids.
316 for (int lowerLevel = 0; lowerLevel< gridIdx; ++lowerLevel) {
317 for( int c=0; c<4; ++c ) {
318 myId += grid_.levelData()[lowerLevel]->indexSet().size( c );
319 }
320 }
321 // Count (and add to myId) all the entities of the refined level grid of codim < 3.
322 for( int c=0; c<3; ++c ) {
323 myId += grid_.indexSet().size( c );
324 }
325 return myId + e.index();
326 }
327 }
328 else { // Leaf grid view (grid view with mixed coarse and refined cells).
329 assert( grid_.getGridIdx() == (static_cast<int>(grid_.levelData().size()) -1) );
330 // In this case, we search for the ids defined in previous levels
331 // (since each entities must keep its id along the entire hiearchy)
332 std::array<int,2> level_levelIdx = {0,0};
333 level_levelIdx = grid_.corner_history_[e.index()];
334 const auto& levelEntity = cpgrid::Entity<3>(*(grid_.levelData()[level_levelIdx[0]]), level_levelIdx[1], true);
335 return grid_.levelData()[level_levelIdx[0]]->local_id_set_ ->id(levelEntity);
336 }
337 } // end-if-data_.size()>1
338 else { // Case: No LGRs / No refined level grids. Only level 0 grid (GLOBAL grid).
339 for( int c=0; c<3; ++c ) {
340 myId += grid_.indexSet().size( c );
341 }
342 return myId + e.index();
343 }
344 }
345 };
346
347
349 {
350 friend class CpGridData;
351 friend class ReversePointGlobalIdSet;
352 public:
353 typedef int IdType;
354
355 void swap(std::vector<int>& cellMapping,
356 std::vector<int>& faceMapping,
357 std::vector<int>& pointMapping)
358 {
359 idSet_=nullptr;
360 GlobalIdMapping::swap(cellMapping,
361 faceMapping,
362 pointMapping);
363 }
364 LevelGlobalIdSet(std::shared_ptr<const IdSet> ids, const CpGridData* view)
365 : idSet_(ids), view_(view)
366 {}
368 : idSet_(), view_()
369 {}
370 template<int codim>
371 IdType id(const Entity<codim>& e) const
372 {
373 assert(view_ == e.pgrid_);
374 // We need to ask the local id set with the full entity
375 // as it needs to be able to determine the level and other
376 // things that are not available in EntityRep.
377 if(idSet_)
378 return idSet_->id(e);
379 else
380 // This a parallel grid and we need to use the mapping
381 // build from the ids of the sequential grid
382 return this->template getMapping<codim>()[e.index()];
383 }
384 template<int codim>
385 IdType id(const EntityRep<codim>& e) const
386 {
387 if(idSet_)
388 return idSet_->id(e);
389 else
390 return this->template getMapping<codim>()[e.index()];
391 }
392
393 template<int cc>
394 IdType subId(const cpgrid::Entity<0>& e, int i) const
395 {
396 assert(view_ == e.pgrid_);
397 return id(e.template subEntity<cc>(i));
398 }
399
400 IdType subId(const cpgrid::Entity<0>& e, int i, int cc) const;
401
402 template<int codim>
403 IdType getMaxCodimGlobalId()
404 {
405 auto max_elem_codim = std::max_element(this->template getMapping<codim>().begin(),
406 this->template getMapping<codim>().end());
407 return *max_elem_codim;
408 }
409
410 IdType getMaxGlobalId()
411 {
412 /*// Ignore faces
413 auto max_globalId = std::max(getMaxCodimGlobalId<0>(), getMaxCodimGlobalId<1>());
414 max_globalId = std::max(max_globalId, getMaxCodimGlobalId<3>());
415 return max_globalId;*/
416 return std::max(getMaxCodimGlobalId<0>(), getMaxCodimGlobalId<3>());
417 }
418
419 private:
420 std::shared_ptr<const IdSet> idSet_;
421 const CpGridData* view_;
422 };
423
431 {
432 public:
434 using IdType = typename LevelGlobalIdSet::IdType;
435
436 GlobalIdSet(const CpGridData& view);
437
438 template<int codim>
439 IdType id(const Entity<codim>& e) const
440 {
441 return levelIdSet(e.pgrid_).id(e);
442 }
443
444 template<int cc>
445 IdType subId(const cpgrid::Entity<0>& e, int i) const
446 {
447 return levelIdSet(e.pgrid_).template subId<cc>(e, i);
448 }
449
450 IdType subId(const cpgrid::Entity<0>& e, int i, int cc) const;
451
452 void insertIdSet(const CpGridData& view);
453 private:
455 const LevelGlobalIdSet& levelIdSet(const CpGridData* const data) const
456 {
457 auto candidate = idSets_.find(data);
458 assert(candidate != idSets_.end());
459 return *candidate->second;
460 }
462 std::map<const CpGridData* const, std::shared_ptr<const LevelGlobalIdSet>> idSets_;
463 };
464
466 {
467 public:
469 {
470 if(idSet.idSet_)
471 {
472 grid_ = &(idSet.idSet_->grid_);
473 }
474 else
475 {
476 mapping_.reset(new std::unordered_map<int,int>);
477 int localId = 0;
478 for (const auto& globalId: idSet.template getMapping<3>())
479 (*mapping_)[globalId] = localId++;
480 }
481 }
482 int operator[](int i) const
483 {
484 if (mapping_)
485 {
486 return(*mapping_)[i];
487 }
488 else if (grid_)
489 {
490 return i - grid_->size(0) - grid_->size(1) - grid_->size(2);
491 }
492
493 OPM_THROW(std::runtime_error, "No grid or mapping. Should not be here!");
494 }
495 void release()
496 {
497 mapping_.reset(nullptr);
498 }
499 private:
500 std::unique_ptr<std::unordered_map<int,int> > mapping_;
501 const CpGridData* grid_ = nullptr;
502 };
503
504 } // namespace cpgrid
505} // namespace Dune
506
507#endif // OPM_INDEXSETS_HEADER
Struct that hods all the data needed to represent a Cpgrid.
Definition CpGridData.hpp:147
const IndexSet & indexSet() const
Get the index set.
Definition CpGridData.hpp:607
Represents an entity of a given codim, with positive or negative orientation.
Definition EntityRep.hpp:99
int index() const
The (positive) index of an entity.
Definition EntityRep.hpp:126
Definition Entity.hpp:71
Class managing the mappings of local indices to global ids.
Definition GlobalIdMapping.hpp:31
void swap(std::vector< int > &cellMapping, std::vector< int > &faceMapping, std::vector< int > &pointMapping)
Swap data for initialization.
Definition GlobalIdMapping.hpp:38
The global id set for Dune.
Definition Indexsets.hpp:431
typename LevelGlobalIdSet::IdType IdType
The type of the id.
Definition Indexsets.hpp:434
Definition Indexsets.hpp:196
IdType id(const cpgrid::Intersection &intersection) const
return id of intersection (here face number)
Definition Indexsets.hpp:225
Definition Indexsets.hpp:53
bool contains(const EntityType &e) const
Definition Indexsets.hpp:182
int size(int codim) const
Definition Indexsets.hpp:123
IndexType subIndex(const cpgrid::Entity< 0 > &e, int i) const
Definition Indexsets.hpp:157
const Types & types(int codim) const
Definition Indexsets.hpp:100
~IndexSet()
Destructor.
Definition Indexsets.hpp:84
std::vector< GeometryType > Types
Definition Indexsets.hpp:68
IndexType index(const EntityType &e) const
Definition Indexsets.hpp:146
IndexType index(const cpgrid::Entity< cd > &e) const
Definition Indexsets.hpp:135
int size(GeometryType type) const
Definition Indexsets.hpp:109
IndexSet()
Definition Indexsets.hpp:73
const Types & geomTypes(int codim) const
Definition Indexsets.hpp:91
int IndexType
Definition Indexsets.hpp:57
Definition Intersection.hpp:66
Definition Indexsets.hpp:349
Definition Indexsets.hpp:466
Copyright 2019 Equinor AS.
Definition CartesianIndexMapper.hpp:10
Export the type of the entity used as parameter in the index(...) method.
Definition Indexsets.hpp:62