Entity.hpp

Go to the documentation of this file.

//===========================================================================
//
// File: Entity.hpp
//
// Created: Fri May 29 20:26:48 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            Brd Skaflestad     <bard.skaflestad@sintef.no>
//            Antonella Ritorto   <antonella.ritorto@opm-op.com>
//
// $Date$
//
// $Revision$
//
//===========================================================================
 
/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010, 2022 Equinor ASA.
 
  This file is part of The Open Porous Media project  (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef OPM_ENTITY_HEADER
#define OPM_ENTITY_HEADER
 
#include <dune/common/version.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/gridenums.hh>
 
#include "PartitionTypeIndicator.hpp"
#include <opm/grid/cpgrid/DefaultGeometryPolicy.hpp>
 
 
namespace Dune
{
namespace cpgrid
{
 
template<int,int> class Geometry;
template<int,PartitionIteratorType> class Iterator;
class IntersectionIterator;
class HierarchicIterator;
class CpGridData;
class LevelGlobalIdSet;
 
template <int codim>
class Entity : public EntityRep<codim>
{
    friend class LevelGlobalIdSet;
    friend class GlobalIdSet;
    friend class HierarchicIterator;
    friend class CpGridData;
 
public:
    enum { codimension = codim };
    enum { dimension = 3 };
    enum { mydimension = dimension - codimension };
    enum { dimensionworld = 3 };
 
    // the official DUNE names
    typedef Entity    EntitySeed;
 
    template <int cd>
    struct Codim
    {
        typedef cpgrid::Entity<cd> Entity;
    };
 
 
    typedef cpgrid::Geometry<3-codim,3> Geometry;
    typedef Geometry LocalGeometry;
 
    typedef cpgrid::IntersectionIterator LeafIntersectionIterator;
    typedef cpgrid::IntersectionIterator LevelIntersectionIterator;
    typedef cpgrid::HierarchicIterator HierarchicIterator;
 
    typedef double ctype;
 
    //             Entity(const CpGridData& grid, int entityrep)
    //              : EntityRep<codim>(entityrep), pgrid_(&grid)
    //          {
    //          }
 
    Entity()
        : EntityRep<codim>(), pgrid_( 0 )
    {
    }
 
    Entity(const CpGridData& grid, EntityRep<codim> entityrep)
        : EntityRep<codim>(entityrep), pgrid_(&grid)
    {
    }
 
    Entity(const CpGridData& grid, int index_arg, bool orientation_arg)
        : EntityRep<codim>(index_arg, orientation_arg), pgrid_(&grid)
    {
    }
 
    Entity(int index_arg, bool orientation_arg)
        : EntityRep<codim>(index_arg, orientation_arg), pgrid_()
    {
    }
 
    bool operator==(const Entity& other) const
    {
        return EntityRep<codim>::operator==(other)  &&  pgrid_ == other.pgrid_;
    }
 
    bool operator!=(const Entity& other) const
    {
        return !operator==(other);
    }
 
    EntitySeed seed() const
    {
        return EntitySeed( impl() );
    }
 
    const Geometry& geometry() const;
 
    int level() const;
 
    bool isLeaf() const;
 
    bool isRegular() const
    {
        return true;
    }
 
    PartitionType partitionType() const;
 
    GeometryType type() const
    {
        return Dune::GeometryTypes::cube(3 - codim);
    }
 
    unsigned int subEntities ( const unsigned int cc ) const;
 
    template <int cc>
    typename Codim<cc>::Entity subEntity(int i) const;
 
    inline LevelIntersectionIterator ilevelbegin() const;
 
    inline LevelIntersectionIterator ilevelend() const;
 
    inline LeafIntersectionIterator ileafbegin() const;
 
    inline LeafIntersectionIterator ileafend() const;
 
 
    HierarchicIterator hbegin(int) const;
 
    HierarchicIterator hend(int) const;
 
    bool isNew() const
    {
        return false;
    }
 
    bool mightVanish() const
    {
        return false;
    }
 
    bool hasFather() const;
 
    Entity<0> father() const;
 
    Dune::cpgrid::Geometry<3,3> geometryInFather() const;
 
    bool hasBoundaryIntersections() const;
 
    // Mimic Dune entity wrapper
    const Entity& impl() const
    {
        return *this;
    }
 
    Entity& impl()
    {
        return *this;
    }
 
    bool isValid () const;
 
    Entity<0> getOrigin() const;
    
    Entity<0> getLevelElem() const;
 
    int getLevelCartesianIdx() const;
 
protected:
    const CpGridData* pgrid_;
private:
    // static    to not need any extra storage per Enitity. One object used for all instances
    // constexpr to allow for in-class instantiation
    static constexpr std::array<int,8> in_father_reference_elem_corner_indices_ = {0,1,2,3,4,5,6,7};
};
 
} // namespace cpgrid
} // namespace Dune
 
// now we include the Iterators.hh We need to do this here because for hbegin/hend the compiler
// needs to know the size of hierarchicIterator
#include "Iterators.hpp"
#include "Intersection.hpp"
namespace Dune
{
namespace cpgrid
{
template<int codim>
typename Entity<codim>::LevelIntersectionIterator Entity<codim>::ilevelbegin() const
{
    static_assert(codim == 0, "");
    return LevelIntersectionIterator(*pgrid_, *this, false);
}
 
template<int codim>
typename Entity<codim>::LevelIntersectionIterator Entity<codim>::ilevelend() const
{
    static_assert(codim == 0, "");
    return LevelIntersectionIterator(*pgrid_, *this, true);
}
 
template<int codim>
typename Entity<codim>::LeafIntersectionIterator Entity<codim>::ileafbegin() const
{
    static_assert(codim == 0, "");
    return LeafIntersectionIterator(*pgrid_, *this, false);
}
 
template<int codim>
typename Entity<codim>::LeafIntersectionIterator Entity<codim>::ileafend() const
{
    static_assert(codim == 0, "");
    return LeafIntersectionIterator(*pgrid_, *this, true);
}
 
 
template<int codim>
HierarchicIterator Entity<codim>::hbegin(int maxLevel) const
{
    // Creates iterator with first child as target if there is one. Otherwise empty stack and target.
    return HierarchicIterator(*this, maxLevel);
}
 
template<int codim>
HierarchicIterator Entity<codim>::hend(int maxLevel) const
{
    // Creates iterator with empty stack and target.
    return HierarchicIterator(maxLevel);
}
 
template <int codim>
PartitionType Entity<codim>::partitionType() const
{
    return pgrid_->partition_type_indicator_->getPartitionType(*this);
}
} // namespace cpgrid
} // namespace Dune
 
 
#include <opm/grid/cpgrid/CpGridData.hpp>
 
namespace Dune {
namespace cpgrid {
 
template<int codim>
unsigned int Entity<codim>::subEntities ( const unsigned int cc ) const
{
    if (cc == codim) {
        return 1;
    }
    else if ( codim == 0 ){ // Cell/element/Entity<0>
        if ( cc == 3 ) { // Get number of corners of the element.
            return 8;
        }
    }
    return 0;
}
 
template <int codim>
const typename Entity<codim>::Geometry& Entity<codim>::geometry() const
{
    return pgrid_->geomVector<codim>()[*this];
}
 
template <int codim>
template <int cc>
typename Entity<codim>::template Codim<cc>::Entity Entity<codim>::subEntity(int i) const
{
    static_assert(codim == 0, "");
    if (cc == 0) { // Cell/element/Entity<0>
        assert(i == 0);
        typename Codim<cc>::Entity se(*pgrid_, this->index(), this->orientation());
        return se;
    } else if (cc == 3) { // Corner/Entity<3>
        assert(i >= 0 && i < 8);
        int corner_index = pgrid_->cell_to_point_[this->index()][i];
        typename Codim<cc>::Entity se(*pgrid_, corner_index, true);
        return se;
    }
    else {
        OPM_THROW(std::runtime_error,
                  "No subentity exists of codimension " + std::to_string(cc));
    }
}
 
template <int codim>
bool Entity<codim>::hasBoundaryIntersections() const
{
    // Copied implementation from EntityDefaultImplementation,
    // except for not checking LevelIntersectionIterators.
    typedef LeafIntersectionIterator Iter;
    Iter end = ileafend();
    for (Iter it = ileafbegin(); it != end; ++it) {
        if (it->boundary()) return true;
    }
    return false;
}
 
template <int codim>
bool Entity<codim>::isValid() const
{
    return pgrid_ ?  EntityRep<codim>::index() < pgrid_->size(codim) : false;
}
 
 
// level() It simply returns the level of the entity in the grid hierarchy.
template <int codim>
int Entity<codim>::level() const
{
    // if distributed_data_ is not empty, level_data_ptr_ has size 1.
    if ((*(pgrid_ -> level_data_ptr_)).size() == 1){
        return 0; // when there is no refinenment, level_ is not automatically instantiated
    }
    if (pgrid_ == (*(pgrid_->level_data_ptr_)).back().get()) { // entity on the leafview -> get the level where it was born:
        return pgrid_ -> leaf_to_level_cells_[this-> index()][0]; // leaf_to_level_cells_ leafIdx -> {level/LGR, cell index in LGR}
    }
    else {
        return pgrid_-> level_;
    }
}
 
 
// isLeaf()
// - if distributed_data_ is empty: an element is a leaf <-> hbegin and hend return the same iterator. Then,
//    *cells from level 0 (coarse grid) that are not parents, are Leaf.
//    *cells from any LGR are Leaf, since they do not children (nested refinement not supported).
// - if distrubuted_data_ is NOT empty: there may be children on a different process. Therefore,
// isLeaf() returns true <-> the element is a leaf entity of the global refinement hierarchy. Equivalently,
// it can be checked whether parent_to_children_cells_ is empty.
 
template<int codim>
bool Entity<codim>::isLeaf() const
{
    if ((pgrid_ -> parent_to_children_cells_).empty()){ // LGR cells
        return true;
    }
    else {
        return (std::get<0>((pgrid_ -> parent_to_children_cells_)[this-> index()]) == -1);  // Cells from GLOBAL, not involved in any LGR
    }
}
 
template<int codim>
bool Entity<codim>::hasFather() const
{
    if ((pgrid_ -> child_to_parent_cells_.empty()) || (pgrid_ -> child_to_parent_cells_[this->index()][0] == -1)){
        return false;
    }
    else{
        return true;
    }
}
 
template<int codim>
Entity<0> Entity<codim>::father() const
{
    if (this->hasFather()){
        const int& coarse_level = pgrid_ -> child_to_parent_cells_[this->index()][0]; // currently, always 0
        const int& parent_index = pgrid_ -> child_to_parent_cells_[this->index()][1];
        const auto& coarse_grid = (*(pgrid_ -> level_data_ptr_))[coarse_level].get();
        return Entity<0>( *coarse_grid, parent_index, true);
    }
    else{
        OPM_THROW(std::logic_error, "Entity has no father.");
    }
}
 
template<int codim>
Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() const
{
    if (!(this->hasFather())){
        OPM_THROW(std::logic_error, "Entity has no father.");
    }
    else{
        // Get IJK index of the entity.
        std::array<int,3> eIJK;
        // Get the amount of children cell in each direction of the parent cell of the entity (same for all parents of each LGR)
        std::array<int,3> cells_per_dim;
        // Get "child0" IJK in the LGR
        std::array<int,3> child0_IJK;
        const auto& level0_grid =  (*(pgrid_->level_data_ptr_))[0];
        const auto& child0_Idx = std::get<1>((*level0_grid).parent_to_children_cells_[this->father().index()])[0];
        // 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.
        if (pgrid_ == (*(pgrid_->level_data_ptr_)).back().get()) // checking if pgrid_ is the LeafView
        {
            const auto& lgr_grid = (*(pgrid_->level_data_ptr_))[this -> level()];
            cells_per_dim = (*lgr_grid).cells_per_dim_;
 
            const auto& entity_lgrIdx = pgrid_ -> leaf_to_level_cells_[this->index()][1]; // leaf_to_level_cells_[cell] = {level, index}
            (*lgr_grid).getIJK(entity_lgrIdx, eIJK);
            (*lgr_grid).getIJK(child0_Idx, child0_IJK);
        }
        else // Getting grid dimension and IJK entity index when pgrid_ is an LGR
        {
            pgrid_ -> getIJK(this->index(), eIJK);
            cells_per_dim = pgrid_ -> cells_per_dim_;
            pgrid_ -> getIJK(child0_Idx, child0_IJK);
        }
        // Transform the local coordinates that comes from the refinemnet in such a way that the
        // reference element of each parent cell is the unit cube. Here, (eIJK[*]-"shift")/cells_per_dim[*]
        // Get the local coordinates of the entity (in the reference unit cube).
        FieldVector<double, 3> corners_in_father_reference_elem_temp[8] = {
            // corner '0'
            { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
            // corner '1'
            { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
            // corner '2'
            { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
            // corner '3'
            { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2])/cells_per_dim[2] },
            // corner '4'
            { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] },
            // corner '5'
            { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1])/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] },
            // corner '6'
            { double(eIJK[0]-child0_IJK[0])/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] },
            // corner '7'
            { double(eIJK[0]-child0_IJK[0]+1)/cells_per_dim[0], double(eIJK[1]-child0_IJK[1]+1)/cells_per_dim[1],
              double(eIJK[2]-child0_IJK[2]+1)/cells_per_dim[2] }};
        auto in_father_reference_elem_corners = std::make_shared<EntityVariable<cpgrid::Geometry<0, 3>, 3>>();
        EntityVariableBase<cpgrid::Geometry<0, 3>>& mutable_in_father_reference_elem_corners = *in_father_reference_elem_corners;
        // Assign the corners. Make use of the fact that pointers behave like iterators.
        mutable_in_father_reference_elem_corners.assign(corners_in_father_reference_elem_temp,
                                                        corners_in_father_reference_elem_temp + 8);
        // Compute the center of the 'local-entity'.
        Dune::FieldVector<double, 3> center_in_father_reference_elem = {0., 0.,0.};
        for (int corn = 0; corn < 8; ++corn) {
            for (int c = 0; c < 3; ++c)
            {
                center_in_father_reference_elem[c] += corners_in_father_reference_elem_temp[corn][c]/8.;
            }
        }
        // Compute the volume of the 'local-entity'.
        double volume_in_father_reference_elem = double(1)/(cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2]);
        // Construct (and return) the Geometry<3,3> of 'child-cell in the reference element of its father (unit cube)'.
        return Dune::cpgrid::Geometry<3,3>(center_in_father_reference_elem, volume_in_father_reference_elem,
                                           in_father_reference_elem_corners, in_father_reference_elem_corner_indices_.data());
    }
}
 
 
template<int codim>
Dune::cpgrid::Entity<0> Dune::cpgrid::Entity<codim>::getOrigin() const
{
    if (hasFather())
    {
        return this->father(); // currently, always a level 0 entity.
    }
    if (!(pgrid_ -> leaf_to_level_cells_.empty())) // entity on the LeafGridView
    {
        const int& entityIdxInLevel0 = pgrid_->leaf_to_level_cells_[this->index()][1]; // leaf_to_level_cells_ [leaf idx] = {0, cell idx}
        const auto& coarse_grid = (*(pgrid_ -> level_data_ptr_))[0].get();
        return Dune::cpgrid::Entity<0>( *coarse_grid, entityIdxInLevel0, true);
    }
    else
    {
        return *this;
    }
}
 
template<int codim>
Dune::cpgrid::Entity<0> Dune::cpgrid::Entity<codim>::getLevelElem() const
{
    // Check that the element belongs to the leaf grid view
    // This is needed to get the index of the element in the level it was born.
    // For cells born in level 0, it is enough to do elem.getOrigin().index().
    // For cells born in LGRs/levels>0, the level index is stored in leaf_to_level_cells_
    // leaf_to_level_cells_ [leaf idx] = {level, cell idx}
    if (!(pgrid_ -> leaf_to_level_cells_.empty())) // entity on the LeafGridView
    {
        const auto entityLevel = this -> level();
        const int& entityLevelIdx = pgrid_->leaf_to_level_cells_[this->index()][1];
        const auto& level_grid = (*(pgrid_ -> level_data_ptr_))[entityLevel].get();
        return Dune::cpgrid::Entity<0>( *level_grid, entityLevelIdx, true);
    }
    else {
        throw std::invalid_argument("The entity provided does not belong to the leaf grid view. ");
    }
}
 
template<int codim>
int Dune::cpgrid::Entity<codim>::getLevelCartesianIdx() const
{
    const auto entityLevel = this -> level();
    const auto level = (*(pgrid_ -> level_data_ptr_))[entityLevel].get();
    const auto& elemInLevel = this->getLevelElem(); // throws when the entity does not belong to the leaf grid view.
    return level -> global_cell_[elemInLevel.index()];
}
 
} // namespace cpgrid
} // namespace Dune
 
 
#endif // OPM_ENTITY_HEADER