Entity.hpp

//===========================================================================
//
// 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/geometry/referenceelements.hh>
#include <dune/grid/common/gridenums.hh>
 
#include "PartitionTypeIndicator.hpp"
#include <opm/grid/cpgrid/DefaultGeometryPolicy.hpp>
 
// To be able to test local and global ids of vertices
void refinePatch_and_check(Dune::CpGrid&,
                           const std::vector<std::array<int,3>>&,
                           const std::vector<std::array<int,3>>&,
                           const std::vector<std::array<int,3>>&,
                           const std::vector<std::string>&);
 
namespace Dune
{
namespace cpgrid
{
 
template <int codim>
class Entity;
 
namespace Impl {
// Forward declaration of traits class for codimensions.
template <int codim>
struct CodimTraits;
 
// Specialization of traits class for codimension 0 (cells).
template<>
struct CodimTraits<0>
{
    using Entity = ::Dune::cpgrid::Entity<0>;
};
 
#if DUNE_VERSION_LT(DUNE_GRID, 2, 11)
// Specialization of traits class for codimension 1 (faces).
// This is only needed for DUNE < 2.11, as the face entity type was needed for the VTK writer, which was fixed in DUNE 2.11.
// See: https://gitlab.dune-project.org/core/dune-grid/-/merge_requests/793
template<>
struct CodimTraits<1>
{
    using Entity = ::Dune::cpgrid::Entity<1>;
};
#endif
 
// Specialization of traits class for codimension 3 (vertices).
template<>
struct CodimTraits<3>
{
    using Entity = ::Dune::cpgrid::Entity<3>;
};
 
}
 
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;
    friend void ::refinePatch_and_check(Dune::CpGrid&,
                                        const std::vector<std::array<int,3>>&,
                                        const std::vector<std::array<int,3>>&,
                                        const std::vector<std::array<int,3>>&,
                                        const std::vector<std::string>&);
 
public:
    static constexpr int codimension = codim;
    static constexpr int dimension = 3;
    static constexpr int mydimension = dimension - codimension;
    static constexpr int dimensionworld = 3;
 
    // the official DUNE names
    typedef Entity    EntitySeed;

    template <int cd>
    using Codim = typename Impl::CodimTraits<cd>;
 
    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;

    bool mightVanish() const;

    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;
 
    int getIdxInParentCell() const;
 
protected:
    const CpGridData* pgrid_;
};
 
// Reference element for a given entity
template <int codim>
auto referenceElement(const Entity<codim>&)
{
    return Dune::referenceElement<double,3-codim>(Dune::GeometryTypes::cube(3-codim));
}
 
 
} // 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;
        }
    }
    DUNE_THROW(NotImplemented, "subEntities not implemented for this codimension");
}
 
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
{
    // Parallel and LGRs:
    // If the grid has been distributed and - after that - LGRs have been added, then level_data_ptr_
    // points at distributed_data_ which has size > 1.
    //
    // Serial and LGRs:
    // If the grid is not distributed and there LGRs have been added, level_data_ptr_ points at data_ which
    // has size > 1.
    //
    // - leaf_to_level_cells_ is non-empty only on the leaf grid view of a grid that has been refined.
    // - level_ is set equal to zero when instantiating a pgrid, and rewriten when it corresponds to a refined level grid.
    return pgrid_->leaf_to_level_cells_.empty()? pgrid_->level_ : pgrid_->leaf_to_level_cells_[this-> index()][0];
}
 
// 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()){ // Grid cells without nested refinement
        return true;
    }
    else {
        return (std::get<0>((pgrid_ -> parent_to_children_cells_)[this-> index()]) == -1);  // Not involved in any LGR
    }
}
 
template<int codim>
bool Entity<codim>::isNew() const
{
    int data_count = pgrid_->level_data_ptr_->size();
    if (data_count == 1) { // CpGrid has only level zero grid.
        return false; // no element is new
    }
    else {
        assert(data_count >= 2); // Level zero and leaf grids, plus at least one refined level grid
        // In nested refinement, an element created during the last refinement step
        // may itself be refined further within the same step at a higher level.
        // To handle this case, we also check whether the element is a leaf
        // (i.e., it has no children).
        return isLeaf() && (pgrid_->refinement_max_level_ == data_count - 2); // minus "level zero" and "leaf"
    }
}
 
template<int codim>
bool Entity<codim>::mightVanish() const
{
    return false;
}
 
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& coarser_level = pgrid_ -> child_to_parent_cells_[this->index()][0];
        const int& parent_cell_index = pgrid_ -> child_to_parent_cells_[this->index()][1];
        return Entity<0>( *((*(pgrid_ -> level_data_ptr_))[coarser_level].get()), parent_cell_index, true);
    }
    else{
        OPM_THROW(std::logic_error, "Entity has no father.");
    }
}
 
template<int codim>
int Dune::cpgrid::Entity<codim>::getIdxInParentCell() const
{
    return pgrid_ -> cell_to_idxInParentCell_[this->index()];
}
 
 
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.");
    }
 
    // Indices of corners in entity's geometry in father reference element.
    static constexpr std::array<int,8> in_father_reference_elem_corner_indices = {0,1,2,3,4,5,6,7};
    // 'static': The returned object Geometry<3,3> stores a pointer to in_father_reference_elem_corner_indices. Therefore,
    // this variable is declared static to prolongate its lifetime beyond this function (static storage duration).
 
    auto idx_in_parent_cell = pgrid_ -> cell_to_idxInParentCell_[this->index()];
    if (idx_in_parent_cell !=-1) {
        const auto& cells_per_dim =  (*(pgrid_ -> level_data_ptr_))[this->level()] -> cells_per_dim_;
        const auto& auxArr = pgrid_ -> getReferenceRefinedCorners(idx_in_parent_cell, cells_per_dim);
        FieldVector<double, 3> corners_in_father_reference_elem_temp[8] =
            { auxArr[0], auxArr[1], auxArr[2], auxArr[3], auxArr[4], auxArr[5], auxArr[6], auxArr[7]};
        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());
    }
    else {
        OPM_THROW(std::logic_error, "Entity has no father.");
    }
}
 
template<int codim>
Dune::cpgrid::Entity<0> Dune::cpgrid::Entity<codim>::getOrigin() const
{
    if (hasFather()) { // Entit is a refined cell belonging to
        // a refined level grid (level>0) or to the leaf grid.
        auto ancestor = this->father();
        while (ancestor.hasFather()){
            ancestor = ancestor.father();
        }
        assert(ancestor.level() == 0);
        return ancestor;
    }
    else if (!(pgrid_ -> leaf_to_level_cells_.empty())) { // Entity is a coarse cell
        // (born in level zero grid, never involved in refinement) belonging to the leaf grid.
        // leaf_to_level_cells_ [leaf idx] = { level where entity was born, cell idx in that level}
        const int& level = pgrid_->leaf_to_level_cells_[this->index()][0];
        assert(level == 0);
        const int& levelElemIdx = pgrid_->leaf_to_level_cells_[this->index()][1];
        return Dune::cpgrid::Entity<0>( *((*(pgrid_ -> level_data_ptr_))[level].get()), levelElemIdx, true);
    }
    else { // Entity belongs to level zero grid.
        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.
    // leaf_to_level_cells_ [leaf idx] = {level where the entity was born, equivalent cell idx in that level}
    if (!(pgrid_ -> leaf_to_level_cells_.empty())) // entity on the LeafGridView
    {
        const int& entityLevelIdx = pgrid_->leaf_to_level_cells_[this->index()][1];
        return Dune::cpgrid::Entity<0>( *((*(pgrid_ -> level_data_ptr_))[this->level()].get()), entityLevelIdx, true);
    }
    else {
        return *this;
    }
}
 
template<int codim>
int Dune::cpgrid::Entity<codim>::getLevelCartesianIdx() const
{
    const auto& level_data = (*(pgrid_ -> level_data_ptr_))[level()].get();
    return level_data -> global_cell_[getLevelElem().index()];
}
 
} // namespace cpgrid
} // namespace Dune
 
 
#endif // OPM_ENTITY_HEADER