My Project
Loading...
Searching...
No Matches
Entity.hpp
1//===========================================================================
2//
3// File: Entity.hpp
4//
5// Created: Fri May 29 20:26:48 2009
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// Brd Skaflestad <bard.skaflestad@sintef.no>
9// Antonella Ritorto <antonella.ritorto@opm-op.com>
10//
11// $Date$
12//
13// $Revision$
14//
15//===========================================================================
16
17/*
18 Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
19 Copyright 2009, 2010, 2022 Equinor ASA.
20
21 This file is part of The Open Porous Media project (OPM).
22
23 OPM is free software: you can redistribute it and/or modify
24 it under the terms of the GNU General Public License as published by
25 the Free Software Foundation, either version 3 of the License, or
26 (at your option) any later version.
27
28 OPM is distributed in the hope that it will be useful,
29 but WITHOUT ANY WARRANTY; without even the implied warranty of
30 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 GNU General Public License for more details.
32
33 You should have received a copy of the GNU General Public License
34 along with OPM. If not, see <http://www.gnu.org/licenses/>.
35*/
36
37#ifndef OPM_ENTITY_HEADER
38#define OPM_ENTITY_HEADER
39
40#include <dune/common/version.hh>
41#include <dune/geometry/type.hh>
42#include <dune/grid/common/gridenums.hh>
43
44#include "PartitionTypeIndicator.hpp"
45#include <opm/grid/cpgrid/DefaultGeometryPolicy.hpp>
46
47
48namespace Dune
49{
50namespace cpgrid
51{
52
53template<int,int> class Geometry;
54template<int,PartitionIteratorType> class Iterator;
55class IntersectionIterator;
56class HierarchicIterator;
57class CpGridData;
58class LevelGlobalIdSet;
59
63template <int codim>
64class Entity : public EntityRep<codim>
65{
66 friend class LevelGlobalIdSet;
67 friend class GlobalIdSet;
68 friend class HierarchicIterator;
69 friend class CpGridData;
70
71public:
74 enum { codimension = codim };
75 enum { dimension = 3 };
76 enum { mydimension = dimension - codimension };
77 enum { dimensionworld = 3 };
78
79 // the official DUNE names
80 typedef Entity EntitySeed;
81
85 template <int cd>
86 struct Codim
87 {
89 };
90
91
92 typedef cpgrid::Geometry<3-codim,3> Geometry;
93 typedef Geometry LocalGeometry;
94
98
99 typedef double ctype;
100
105 // Entity(const CpGridData& grid, int entityrep)
106 // : EntityRep<codim>(entityrep), pgrid_(&grid)
107 // {
108 // }
109
112 : EntityRep<codim>(), pgrid_( 0 )
113 {
114 }
115
117 Entity(const CpGridData& grid, EntityRep<codim> entityrep)
118 : EntityRep<codim>(entityrep), pgrid_(&grid)
119 {
120 }
121
123 Entity(const CpGridData& grid, int index_arg, bool orientation_arg)
124 : EntityRep<codim>(index_arg, orientation_arg), pgrid_(&grid)
125 {
126 }
127
129 Entity(int index_arg, bool orientation_arg)
130 : EntityRep<codim>(index_arg, orientation_arg), pgrid_()
131 {
132 }
133
135 bool operator==(const Entity& other) const
136 {
137 return EntityRep<codim>::operator==(other) && pgrid_ == other.pgrid_;
138 }
139
141 bool operator!=(const Entity& other) const
142 {
143 return !operator==(other);
144 }
145
150 {
151 return EntitySeed( impl() );
152 }
153
155 const Geometry& geometry() const;
156
158 int level() const;
159
165 bool isLeaf() const;
166
168 bool isRegular() const
169 {
170 return true;
171 }
172
175 PartitionType partitionType() const;
176
179 GeometryType type() const
180 {
181 return Dune::GeometryTypes::cube(3 - codim);
182 }
183
185 unsigned int subEntities ( const unsigned int cc ) const;
186
189 template <int cc>
190 typename Codim<cc>::Entity subEntity(int i) const;
191
194
197
200
203
204
207
210
212 bool isNew() const
213 {
214 return false;
215 }
216
218 bool mightVanish() const
219 {
220 return false;
221 }
222
228 bool hasFather() const;
229
234
241
247
248 // Mimic Dune entity wrapper
250 const Entity& impl() const
251 {
252 return *this;
253 }
254
255 Entity& impl()
256 {
257 return *this;
258 }
259
262 bool isValid () const;
263
268
269protected:
270 const CpGridData* pgrid_;
271private:
272 // static to not need any extra storage per Enitity. One object used for all instances
273 // constexpr to allow for in-class instantiation
275 static constexpr std::array<int,8> in_father_reference_elem_corner_indices_ = {0,1,2,3,4,5,6,7};
276};
277
278} // namespace cpgrid
279} // namespace Dune
280
281// now we include the Iterators.hh We need to do this here because for hbegin/hend the compiler
282// needs to know the size of hierarchicIterator
283#include "Iterators.hpp"
284#include "Intersection.hpp"
285namespace Dune
286{
287namespace cpgrid
288{
289template<int codim>
291{
292 static_assert(codim == 0, "");
293 return LevelIntersectionIterator(*pgrid_, *this, false);
294}
295
296template<int codim>
298{
299 static_assert(codim == 0, "");
300 return LevelIntersectionIterator(*pgrid_, *this, true);
301}
302
303template<int codim>
305{
306 static_assert(codim == 0, "");
307 return LeafIntersectionIterator(*pgrid_, *this, false);
308}
309
310template<int codim>
312{
313 static_assert(codim == 0, "");
314 return LeafIntersectionIterator(*pgrid_, *this, true);
315}
316
317
318template<int codim>
320{
321 // Creates iterator with first child as target if there is one. Otherwise empty stack and target.
322 return HierarchicIterator(*this, maxLevel);
323}
324
326template<int codim>
328{
329 // Creates iterator with empty stack and target.
330 return HierarchicIterator(maxLevel);
331}
332
333template <int codim>
334PartitionType Entity<codim>::partitionType() const
335{
336 return pgrid_->partition_type_indicator_->getPartitionType(*this);
337}
338} // namespace cpgrid
339} // namespace Dune
340
341
342#include <opm/grid/cpgrid/CpGridData.hpp>
343
344namespace Dune {
345namespace cpgrid {
346
347template<int codim>
348unsigned int Entity<codim>::subEntities ( const unsigned int cc ) const
349{
350 if (cc == codim) {
351 return 1;
352 }
353 else if ( codim == 0 ){ // Cell/element/Entity<0>
354 if ( cc == 3 ) { // Get number of corners of the element.
355 return 8;
356 }
357 }
358 return 0;
359}
360
361template <int codim>
363{
364 return pgrid_->geomVector<codim>()[*this];
365}
366
367template <int codim>
368template <int cc>
369typename Entity<codim>::template Codim<cc>::Entity Entity<codim>::subEntity(int i) const
370{
371 static_assert(codim == 0, "");
372 if (cc == 0) { // Cell/element/Entity<0>
373 assert(i == 0);
374 typename Codim<cc>::Entity se(*pgrid_, this->index(), this->orientation());
375 return se;
376 } else if (cc == 3) { // Corner/Entity<3>
377 assert(i >= 0 && i < 8);
378 int corner_index = pgrid_->cell_to_point_[this->index()][i];
379 typename Codim<cc>::Entity se(*pgrid_, corner_index, true);
380 return se;
381 }
382 else {
383 OPM_THROW(std::runtime_error,
384 "No subentity exists of codimension " + std::to_string(cc));
385 }
386}
387
388template <int codim>
390{
391 // Copied implementation from EntityDefaultImplementation,
392 // except for not checking LevelIntersectionIterators.
393 typedef LeafIntersectionIterator Iter;
394 Iter end = ileafend();
395 for (Iter it = ileafbegin(); it != end; ++it) {
396 if (it->boundary()) return true;
397 }
398 return false;
399}
400
401template <int codim>
403{
404 return pgrid_ ? EntityRep<codim>::index() < pgrid_->size(codim) : false;
405}
406
407
408// level() It simply returns the level of the entity in the grid hierarchy.
409template <int codim>
411{
412 // if distributed_data_ is not empty, level_data_ptr_ has size 1.
413 if ((*(pgrid_ -> level_data_ptr_)).size() == 1){
414 return 0; // when there is no refinenment, level_ is not automatically instantiated
415 }
416 if (pgrid_ == (*(pgrid_->level_data_ptr_)).back().get()) { // entity on the leafview -> get the level where it was born:
417 return pgrid_ -> leaf_to_level_cells_[this-> index()][0]; // leaf_to_level_cells_ leafIdx -> {level/LGR, cell index in LGR}
418 }
419 else {
420 return pgrid_-> level_;
421 }
422}
423
424
425// isLeaf()
426// - if distributed_data_ is empty: an element is a leaf <-> hbegin and hend return the same iterator. Then,
427// *cells from level 0 (coarse grid) that are not parents, are Leaf.
428// *cells from any LGR are Leaf, since they do not children (nested refinement not supported).
429// - if distrubuted_data_ is NOT empty: there may be children on a different process. Therefore,
430// isLeaf() returns true <-> the element is a leaf entity of the global refinement hierarchy. Equivalently,
431// it can be checked whether parent_to_children_cells_ is empty.
432
433template<int codim>
435{
436 if ((pgrid_ -> parent_to_children_cells_).empty()){ // LGR cells
437 return true;
438 }
439 else {
440 return (std::get<0>((pgrid_ -> parent_to_children_cells_)[this-> index()]) == -1); // Cells from GLOBAL, not involved in any LGR
441 }
442}
443
444template<int codim>
446{
447 if ((pgrid_ -> child_to_parent_cells_.empty()) || (pgrid_ -> child_to_parent_cells_[this->index()][0] == -1)){
448 return false;
449 }
450 else{
451 return true;
452 }
453}
454
455template<int codim>
457{
458 if (this->hasFather()){
459 const int& coarse_level = pgrid_ -> child_to_parent_cells_[this->index()][0]; // currently, always 0
460 const int& parent_index = pgrid_ -> child_to_parent_cells_[this->index()][1];
461 const auto& coarse_grid = (*(pgrid_ -> level_data_ptr_))[coarse_level].get();
462 return Entity<0>( *coarse_grid, parent_index, true);
463 }
464 else{
465 OPM_THROW(std::logic_error, "Entity has no father.");
466 }
467}
468
469template<int codim>
471{
472 if (!(this->hasFather())){
473 OPM_THROW(std::logic_error, "Entity has no father.");
474 }
475 else{
476 // Get IJK index of the entity.
477 std::array<int,3> eIJK;
478 // Get the amount of children cell in each direction of the parent cell of the entity (same for all parents of each LGR)
479 std::array<int,3> cells_per_dim;
480 // Get "child0" IJK in the LGR
481 std::array<int,3> child0_IJK;
482 const auto& level0_grid = (*(pgrid_->level_data_ptr_))[0];
483 const auto& child0_Idx = std::get<1>((*level0_grid).parent_to_children_cells_[this->father().index()])[0];
484 // If pgrid_ is the leafview, go to the LGR where the entity was born to get its IJK index in the LGR and the LGR dimension.
485 if (pgrid_ == (*(pgrid_->level_data_ptr_)).back().get()) // checking if pgrid_ is the LeafView
486 {
487 const auto& lgr_grid = (*(pgrid_->level_data_ptr_))[this -> level()];
488 cells_per_dim = (*lgr_grid).cells_per_dim_;
489
490 const auto& entity_lgrIdx = pgrid_ -> leaf_to_level_cells_[this->index()][1]; // leaf_to_level_cells_[cell] = {level, index}
491 (*lgr_grid).getIJK(entity_lgrIdx, eIJK);
492 (*lgr_grid).getIJK(child0_Idx, child0_IJK);
493 }
494 else // Getting grid dimension and IJK entity index when pgrid_ is an LGR
495 {
496 pgrid_ -> getIJK(this->index(), eIJK);
497 cells_per_dim = pgrid_ -> cells_per_dim_;
498 pgrid_ -> getIJK(child0_Idx, child0_IJK);
499 }
500 // Transform the local coordinates that comes from the refinemnet in such a way that the
501 // reference element of each parent cell is the unit cube. Here, (eIJK[*]-"shift")/cells_per_dim[*]
502 // Get the local coordinates of the entity (in the reference unit cube).
503 FieldVector<double, 3> corners_in_father_reference_elem_temp[8] = {
504 // corner '0'
505 { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
506 double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
507 // corner '1'
508 { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
509 double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
510 // corner '2'
511 { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
512 double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
513 // corner '3'
514 { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
515 double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
516 // corner '4'
517 { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
518 double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] },
519 // corner '5'
520 { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
521 double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] },
522 // corner '6'
523 { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
524 double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] },
525 // corner '7'
526 { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
527 double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] }};
528 auto in_father_reference_elem_corners = std::make_shared<EntityVariable<cpgrid::Geometry<0, 3>, 3>>();
529 EntityVariableBase<cpgrid::Geometry<0, 3>>& mutable_in_father_reference_elem_corners = *in_father_reference_elem_corners;
530 // Assign the corners. Make use of the fact that pointers behave like iterators.
531 mutable_in_father_reference_elem_corners.assign(corners_in_father_reference_elem_temp,
532 corners_in_father_reference_elem_temp + 8);
533 // Compute the center of the 'local-entity'.
534 Dune::FieldVector<double, 3> center_in_father_reference_elem = {0., 0.,0.};
535 for (int corn = 0; corn < 8; ++corn) {
536 for (int c = 0; c < 3; ++c)
537 {
538 center_in_father_reference_elem[c] += corners_in_father_reference_elem_temp[corn][c]/8.;
539 }
540 }
541 // Compute the volume of the 'local-entity'.
542 double volume_in_father_reference_elem = double(1)/(cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2]);
543 // Construct (and return) the Geometry<3,3> of 'child-cell in the reference element of its father (unit cube)'.
544 return Dune::cpgrid::Geometry<3,3>(center_in_father_reference_elem, volume_in_father_reference_elem,
545 in_father_reference_elem_corners, in_father_reference_elem_corner_indices_.data());
546 }
547}
548
549
550template<int codim>
552{
553 if (hasFather())
554 {
555 return this->father(); // currently, always a level 0 entity.
556 }
557 if (!(pgrid_ -> leaf_to_level_cells_.empty())) // entity on the LeafGridView
558 {
559 const int& entityIdxInLevel0 = pgrid_->leaf_to_level_cells_[this->index()][1]; // leaf_to_level_cells_ [leaf idx] = {0, cell idx}
560 const auto& coarse_grid = (*(pgrid_ -> level_data_ptr_))[0].get();
561 return Dune::cpgrid::Entity<0>( *coarse_grid, entityIdxInLevel0, true);
562 }
563 else
564 {
565 return *this;
566 }
567}
568
569} // namespace cpgrid
570} // namespace Dune
571
572
573#endif // OPM_ENTITY_HEADER
Struct that hods all the data needed to represent a Cpgrid.
Definition CpGridData.hpp:135
Represents an entity of a given codim, with positive or negative orientation.
Definition EntityRep.hpp:99
bool operator==(const EntityRep &other) const
Equality operator.
Definition EntityRep.hpp:179
int index() const
The (positive) index of an entity.
Definition EntityRep.hpp:126
Base class for EntityVariable and SignedEntityVariable.
Definition EntityRep.hpp:219
Definition Entity.hpp:65
bool mightVanish() const
Returns true, if entity might disappear during the next call to adapt(). Dummy.
Definition Entity.hpp:218
Entity< 0 > father() const
ONLY FOR CELLS (Entity<0>).
Definition Entity.hpp:456
unsigned int subEntities(const unsigned int cc) const
Return the number of all subentities of the entity of a given codimension cc.
Definition Entity.hpp:348
LeafIntersectionIterator ileafend() const
End leaf-iterator for the cell-cell intersections of this entity.
Definition Entity.hpp:311
EntitySeed seed() const
Return an entity seed (light-weight entity).
Definition Entity.hpp:149
HierarchicIterator hend(int) const
Iterator end over the children/beyond last child iterator.
Definition Entity.hpp:327
const Geometry & geometry() const
Return the geometry of the entity (does not depend on its orientation).
Definition Entity.hpp:362
Entity(const CpGridData &grid, int index_arg, bool orientation_arg)
Constructor taking a grid, entity index, and orientation.
Definition Entity.hpp:123
bool hasFather() const
ONLY FOR CELLS (Entity<0>) Check if the entity comes from an LGR, i.e., it has been created via refin...
Definition Entity.hpp:445
bool isValid() const
isValid method for EntitySeed
Definition Entity.hpp:402
Entity< 0 > getOrigin() const
getOrigin() Returns parent entity in level 0, if the entity was born in any LGR.
Definition Entity.hpp:551
bool operator==(const Entity &other) const
Equality.
Definition Entity.hpp:135
bool isRegular() const
Refinement is not defined for CpGrid.
Definition Entity.hpp:168
Entity()
Constructor taking a grid and an integer entity representation.
Definition Entity.hpp:111
HierarchicIterator hbegin(int) const
Iterator begin over the children. [If requested, also over descendants more than one generation away....
Definition Entity.hpp:319
const Entity & impl() const
Access the actual implementation class behind Entity interface class.
Definition Entity.hpp:250
PartitionType partitionType() const
For now, the grid is serial and the only partitionType() is InteriorEntity.
Definition Entity.hpp:334
Dune::cpgrid::Geometry< 3, 3 > geometryInFather() const
Return LocalGeometry representing the embedding of the entity into its father (when hasFather() is tr...
Definition Entity.hpp:470
Entity(int index_arg, bool orientation_arg)
Constructor taking a entity index, and orientation.
Definition Entity.hpp:129
LeafIntersectionIterator ileafbegin() const
Start leaf-iterator for the cell-cell intersections of this entity.
Definition Entity.hpp:304
Entity(const CpGridData &grid, EntityRep< codim > entityrep)
Constructor taking a grid and an entity representation.
Definition Entity.hpp:117
LevelIntersectionIterator ilevelend() const
End level-iterator for the cell-cell intersections of this entity.
Definition Entity.hpp:297
bool isNew() const
Returns true, if the entity has been created during the last call to adapt(). Dummy.
Definition Entity.hpp:212
LevelIntersectionIterator ilevelbegin() const
Start level-iterator for the cell-cell intersections of this entity.
Definition Entity.hpp:290
GeometryType type() const
Return marker object (GeometryType object) representing the reference element of the entity.
Definition Entity.hpp:179
int level() const
Return the level of the entity in the grid hierarchy. Level = 0 represents the coarsest grid.
Definition Entity.hpp:410
bool isLeaf() const
Check if the entity is in the leafview.
Definition Entity.hpp:434
bool operator!=(const Entity &other) const
Inequality.
Definition Entity.hpp:141
bool hasBoundaryIntersections() const
Returns true if any of my intersections are on the boundary.
Definition Entity.hpp:389
Codim< cc >::Entity subEntity(int i) const
Obtain subentity.
This class encapsulates geometry for vertices, intersections, and cells.
Definition Geometry.hpp:75
The global id set for Dune.
Definition Indexsets.hpp:304
Only needs to provide interface for doing nothing.
Definition Iterators.hpp:108
Definition Intersection.hpp:278
Definition Indexsets.hpp:248
Copyright 2019 Equinor AS.
Definition CartesianIndexMapper.hpp:10
Definition Entity.hpp:87