Loading...
Searching...
No Matches
CpGridData.hpp
//===========================================================================
//
// File: CpGridData.hpp
//
// Created: Sep 17 21:11:41 2013
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <bard.skaflestad@sintef.no>
// Markus Blatt <markus@dr-blatt.de>
// Antonella Ritorto <antonella.ritorto@opm-op.com>
//
// Comment: Major parts of this file originated in dune/grid/CpGrid.hpp
// and got transfered here during refactoring for the parallelization.
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2013, 2022-2023 Equinor ASA.
Copyright 2013 Dr. Blatt - HPC-Simulation-Software & Services
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_CPGRIDDATA_HEADER
#define OPM_CPGRIDDATA_HEADER
#include <dune/common/parallel/mpihelper.hh>
#ifdef HAVE_DUNE_ISTL
#include <dune/istl/owneroverlapcopy.hh>
#endif
#include <dune/common/parallel/communication.hh>
#include <dune/common/parallel/variablesizecommunicator.hh>
#include <dune/grid/common/gridenums.hh>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/input/eclipse/EclipseState/Grid/NNC.hpp>
#endif
#include <opm/grid/cpgpreprocess/preprocess.h>
#include "Entity2IndexDataHandle.hpp"
#include "CpGridDataTraits.hpp"
//#include "DataHandleWrappers.hpp"
//#include "GlobalIdMapping.hpp"
#include "Geometry.hpp"
#include <array>
#include <initializer_list>
#include <set>
#include <vector>
namespace Opm
{
class EclipseState;
}
namespace Dune
{
class CpGrid;
namespace cpgrid
{
class IndexSet;
class IdSet;
class LevelGlobalIdSet;
class PartitionTypeIndicator;
template<int,int> class Geometry;
template<int> class Entity;
template<int> class EntityRep;
}
}
void refine_and_check(const Dune::cpgrid::Geometry<3, 3>&,
const std::array<int, 3>&,
bool);
namespace Dune
{
namespace cpgrid
{
namespace mover
{
template<class T, int i> struct Mover;
}
class CpGridData
{
template<class T, int i> friend struct mover::Mover;
friend class GlobalIdSet;
friend class HierarchicIterator;
friend class Dune::cpgrid::IndexSet;
friend class Dune::cpgrid::IdSet;
friend class Dune::cpgrid::LevelGlobalIdSet;
friend
void ::refine_and_check(const Dune::cpgrid::Geometry<3, 3>&,
const std::array<int, 3>&,
bool);
private:
CpGridData(const CpGridData& g);
public:
enum{
#ifndef MAX_DATA_COMMUNICATED_PER_ENTITY
MAX_DATA_PER_CELL = 2
#else
MAX_DATA_PER_CELL = MAX_DATA_COMMUNICATED_PER_ENTITY
#endif
};
CpGridData() = delete;
explicit CpGridData(MPIHelper::MPICommunicator comm, std::vector<std::shared_ptr<CpGridData>>& data);
explicit CpGridData(std::vector<std::shared_ptr<CpGridData>>& data);
~CpGridData();
int size(int codim) const;
int size (GeometryType type) const
{
if (type.isCube()) {
return size(3 - type.dim());
} else {
return 0;
}
}
void readEclipseFormat(const std::string& filename,
bool periodic_extension,
bool turn_normals = false,
bool edge_conformal = false);
#if HAVE_ECL_INPUT
void processEclipseFormat(const Opm::Deck& deck,
bool periodic_extension,
bool turn_normals = false,
bool clip_z = false,
const std::vector<double>& poreVolume = std::vector<double>{},
bool edge_conformal = false);
std::vector<std::size_t>
processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
Opm::EclipseState* ecl_state,
bool periodic_extension,
bool turn_normals = false,
bool clip_z = false,
bool pinchActive = true,
bool edge_conformal = false);
#endif
void processEclipseFormat(const grdecl& input_data,
#if HAVE_ECL_INPUT
Opm::EclipseState* ecl_state,
#endif
std::array<std::set<std::pair<int, int>>, 2>& nnc,
bool remove_ij_boundary,
bool turn_normals,
bool pinchActive,
double tolerance_unique_points,
bool edge_conformal);
void getIJK(int c, std::array<int,3>& ijk) const;
int cellFace(int cell, int local_index) const
{
return cell_to_face_[cpgrid::EntityRep<0>(cell, true)][local_index].index();
}
auto cellToFace(int cellIdx) const
{
return cell_to_face_[cpgrid::EntityRep<0>(cellIdx, true)];
}
const auto& cellToPoint() const
{
return cell_to_point_;
}
const auto& cellToPoint(int cellIdx) const
{
return cell_to_point_[cellIdx];
}
int faceToCellSize(int face) const {
Dune::cpgrid::EntityRep<1> faceRep(face, true);
return face_to_cell_[faceRep].size();
}
auto faceTag(int faceIdx) const
{
Dune::cpgrid::EntityRep<1> faceRep(faceIdx, true);
return face_tag_[faceRep];
}
auto faceNormals(int faceIdx) const
{
Dune::cpgrid::EntityRep<1> faceRep(faceIdx, true);
return face_normals_[faceRep];
}
auto faceToPoint(int faceIdx) const
{
return face_to_point_[faceIdx];
}
int numFaces() const
{
return face_to_cell_.size();
}
auto cornerHistorySize() const
{
return corner_history_.size();
}
const auto& getCornerHistory(int cornerIdx) const
{
if(cornerHistorySize()) {
return corner_history_[cornerIdx];
}
else {
OPM_THROW(std::logic_error, "Vertex has no history record.\n");
}
}
const std::vector<int>& globalCell() const
{
return global_cell_;
}
bool hasNNCs(const std::vector<int>& cellIndices) const;
bool mark(int refCount, const cpgrid::Entity<0>& element);
int getMark(const cpgrid::Entity<0>& element) const;
bool preAdapt();
bool adapt();
void postAdapt();
private:
std::array<Dune::FieldVector<double,3>,8> getReferenceRefinedCorners(int idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const;
public:
int getGridIdx() const {
// Not the nicest way of checking if "this" points at the leaf grid view of a mixed grid (with coarse and refined cells).
// 1. When the grid has been refined at least onece, level_data_ptr_ ->size() >1. Therefore, there is a chance of "this" pointing at the leaf grid view.
// 2. Unfortunately, level_ is default initialized by 0. This implies, in particular, that if someone wants to check the value of
// "this->level_" when "this" points at the leaf grid view of a grid that has been refined, this value is - unfortunately - equal to 0.
// 3. Due to 2. we need an extra bool value to distinguish between the actual level 0 grid and such a leaf grid view (with incorrect level_ == 0). For this
// reason we check if child_to_parent_cells_.empty() [true for actual level 0 grid, false for the leaf grid view].
// --- TO BE IMPROVED ---
if ((level_data_ptr_ ->size() >1) && (level_ == 0) && (!child_to_parent_cells_.empty())) {
return level_data_ptr_->size() -1;
}
return level_;
}
const std::vector<std::shared_ptr<Dune::cpgrid::CpGridData>>& levelData() const
{
if (level_data_ptr_->empty()) {
OPM_THROW(std::logic_error, "Level data has not been initialized\n");
}
return *level_data_ptr_;
}
const std::tuple<int,std::vector<int>>& getChildrenLevelAndIndexList(int elemIdx) const {
return parent_to_children_cells_[elemIdx];
}
const std::vector<std::tuple<int,std::vector<int>>>& getParentToChildren() const {
return parent_to_children_cells_;
}
const cpgrid::DefaultGeometryPolicy getGeometry() const
{
return geometry_;
}
int getLeafIdxFromLevelIdx(int level_cell_idx) const
{
if (level_to_leaf_cells_.empty()) {
OPM_THROW(std::logic_error, "Grid has no LGRs. No mapping to the leaf.\n");
}
return level_to_leaf_cells_[level_cell_idx];
}
std::tuple< const std::shared_ptr<CpGridData>,
const std::vector<std::array<int,2>>, // parent_to_refined_corners(~boundary_old_to_new_corners)
const std::vector<std::tuple<int,std::vector<int>>>, // parent_to_children_faces (~boundary_old_to_new_faces)
const std::tuple<int, std::vector<int>>, // parent_to_children_cells
const std::vector<std::array<int,2>>, // child_to_parent_faces
const std::vector<std::array<int,2>>> // child_to_parent_cells
refineSingleCell(const std::array<int,3>& cells_per_dim, const int& parent_idx) const;
// @breif Compute center of an entity/element/cell in the Eclipse way:
// - Average of the 4 corners of the bottom face.
// - Average of the 4 corners of the top face.
// Return average of the previous computations.
// @param [in] int Index of a cell.
// @return 'eclipse centroid'
std::array<double,3> computeEclCentroid(const int idx) const;
// @breif Compute center of an entity/element/cell in the Eclipse way:
// - Average of the 4 corners of the bottom face.
// - Average of the 4 corners of the top face.
// Return average of the previous computations.
// @param [in] Entity<0> Entity
// @return 'eclipse centroid'
std::array<double,3> computeEclCentroid(const Entity<0>& elem) const;
// Make unique boundary ids for all intersections.
void computeUniqueBoundaryIds();
bool uniqueBoundaryIds() const
{
return use_unique_boundary_ids_;
}
void setUniqueBoundaryIds(bool uids)
{
use_unique_boundary_ids_ = uids;
if (use_unique_boundary_ids_ && unique_boundary_ids_.empty()) {
computeUniqueBoundaryIds();
}
}
const std::vector<double>& zcornData() const {
return zcorn;
}
const IndexSet& indexSet() const
{
return *index_set_;
}
const cpgrid::IdSet& localIdSet() const
{
return *local_id_set_;
}
const cpgrid::LevelGlobalIdSet& globalIdSet() const
{
return *global_id_set_;
}
const std::array<int, 3>& logicalCartesianSize() const
{
return logical_cartesian_size_;
}
void distributeGlobalGrid(CpGrid& grid,
const CpGridData& view_data,
const std::vector<int>& cell_part);
template<class DataHandle>
void communicate(DataHandle& data, InterfaceType iftype, CommunicationDirection dir);
void computeCellPartitionType();
void computePointPartitionType();
void computeCommunicationInterfaces(int noexistingPoints);
using MPICommunicator = CpGridDataTraits::MPICommunicator ;
using Communication = CpGridDataTraits::Communication;
using CollectiveCommunication = CpGridDataTraits::CollectiveCommunication;
using AttributeSet = CpGridDataTraits::AttributeSet;
#if HAVE_MPI
using Communicator = CpGridDataTraits::Communicator;
using InterfaceMap = CpGridDataTraits::InterfaceMap;
using CommunicationType = CpGridDataTraits::CommunicationType;
using ParallelIndexSet = CpGridDataTraits::ParallelIndexSet;
using RemoteIndices = CpGridDataTraits::RemoteIndices;
CommunicationType& cellCommunication()
{
return cell_comm_;
}
const CommunicationType& cellCommunication() const
{
return cell_comm_;
}
ParallelIndexSet& cellIndexSet()
{
return cellCommunication().indexSet();
}
const ParallelIndexSet& cellIndexSet() const
{
return cellCommunication().indexSet();
}
RemoteIndices& cellRemoteIndices()
{
return cellCommunication().remoteIndices();
}
const RemoteIndices& cellRemoteIndices() const
{
return cellCommunication().remoteIndices();
}
#endif
const std::vector<int>& sortedNumAquiferCells() const
{
return aquifer_cells_;
}
private:
void populateGlobalCellIndexSet();
#if HAVE_MPI
template<class DataHandle>
void gatherData(DataHandle& data, CpGridData* global_view,
CpGridData* distributed_view);
template<int codim, class DataHandle>
void gatherCodimData(DataHandle& data, CpGridData* global_data,
CpGridData* distributed_data);
template<class DataHandle>
void scatterData(DataHandle& data, const CpGridData* global_data,
const CpGridData* distributed_data, const InterfaceMap& cell_inf,
const InterfaceMap& point_inf);
template<int codim, class DataHandle>
void scatterCodimData(DataHandle& data, CpGridData* global_data,
CpGridData* distributed_data);
template<int codim, class DataHandle>
void communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
const Interface& interface);
template<int codim, class DataHandle>
void communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
const InterfaceMap& interface);
#endif
void computeGeometry(const CpGrid& grid,
const DefaultGeometryPolicy& globalGeometry,
const std::vector<int>& globalAquiferCells,
const OrientedEntityTable<0, 1>& globalCell2Faces,
DefaultGeometryPolicy& geometry,
std::vector<int>& aquiferCells,
const OrientedEntityTable<0, 1>& cell2Faces,
const std::vector< std::array<int,8> >& cell2Points);
// Representing the topology
cpgrid::OrientedEntityTable<0, 1> cell_to_face_;
cpgrid::OrientedEntityTable<1, 0> face_to_cell_;
Opm::SparseTable<int> face_to_point_;
std::vector< std::array<int,8> > cell_to_point_;
std::array<int, 3> logical_cartesian_size_{};
std::vector<int> global_cell_;
cpgrid::EntityVariable<enum face_tag, 1> face_tag_;
cpgrid::DefaultGeometryPolicy geometry_;
typedef FieldVector<double, 3> PointType;
cpgrid::SignedEntityVariable<PointType, 1> face_normals_;
cpgrid::EntityVariable<int, 1> unique_boundary_ids_;
std::unique_ptr<cpgrid::IndexSet> index_set_;
std::shared_ptr<const cpgrid::IdSet> local_id_set_;
std::shared_ptr<LevelGlobalIdSet> global_id_set_;
std::shared_ptr<PartitionTypeIndicator> partition_type_indicator_;
std::vector<int> mark_;
int level_{0};
std::vector<std::shared_ptr<CpGridData>>* level_data_ptr_;
// SUITABLE FOR ALL LEVELS EXCEPT FOR LEAFVIEW
std::vector<int> level_to_leaf_cells_; // In entry 'level cell index', we store 'leafview cell index' // {level LGR, {child0, child1, ...}}
std::vector<std::tuple<int,std::vector<int>>> parent_to_children_cells_; // {# children in x-direction, ... y-, ... z-}
std::array<int,3> cells_per_dim_;
// SUITABLE ONLY FOR LEAFVIEW // {level, cell index in that level}
std::vector<std::array<int,2>> leaf_to_level_cells_;
std::vector<std::array<int,2>> corner_history_;
// SUITABLE FOR ALL LEVELS INCLUDING LEAFVIEW // {level parent cell, parent cell index}
std::vector<std::array<int,2>> child_to_parent_cells_;
std::vector<int> cell_to_idxInParentCell_;
int refinement_max_level_{0};
Communication ccobj_;
// Boundary information (optional).
bool use_unique_boundary_ids_;
std::vector<double> zcorn;
std::vector<int> aquifer_cells_;
#if HAVE_MPI
CommunicationType cell_comm_;
std::tuple<Interface,Interface,Interface,Interface,Interface> cell_interfaces_;
/*
// code deactivated, because users cannot access face indices and therefore
// communication on faces makes no sense!
std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
face_interfaces_;
*/
std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
point_interfaces_;
#endif
// Return the geometry vector corresponding to the given codim.
template <int codim>
const EntityVariable<Geometry<3 - codim, 3>, codim>& geomVector() const
{
return geometry_.geomVector<codim>();
}
friend class Dune::CpGrid;
template<int> friend class Entity;
template<int> friend class EntityRep;
friend class Intersection;
friend class PartitionTypeIndicator;
};
#if HAVE_MPI
namespace
{
template<class T>
T& getInterface(InterfaceType iftype,
std::tuple<T,T,T,T,T>& interfaces)
{
switch(iftype)
{
case 0:
return std::get<0>(interfaces);
case 1:
return std::get<1>(interfaces);
case 2:
return std::get<2>(interfaces);
case 3:
return std::get<3>(interfaces);
case 4:
return std::get<4>(interfaces);
}
OPM_THROW(std::runtime_error, "Invalid Interface type was used during communication");
}
} // end unnamed namespace
template<int codim, class DataHandle>
void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
const Interface& interface)
{
this->template communicateCodim<codim>(data, dir, interface.interfaces());
}
template<int codim, class DataHandle>
void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data_wrapper, CommunicationDirection dir,
const InterfaceMap& interface)
{
Communicator comm(ccobj_, interface);
if(dir==ForwardCommunication)
comm.forward(data_wrapper);
else
comm.backward(data_wrapper);
}
#endif
template<class DataHandle>
void CpGridData::communicate(DataHandle& data, InterfaceType iftype,
CommunicationDirection dir)
{
#if HAVE_MPI
if(data.contains(3,0))
{
Entity2IndexDataHandle<DataHandle, 0> data_wrapper(*this, data);
communicateCodim<0>(data_wrapper, dir, getInterface(iftype, cell_interfaces_));
}
if(data.contains(3,3))
{
Entity2IndexDataHandle<DataHandle, 3> data_wrapper(*this, data);
communicateCodim<3>(data_wrapper, dir, getInterface(iftype, point_interfaces_));
}
#else
// Suppress warnings for unused arguments.
(void) data;
(void) iftype;
(void) dir;
#endif
}
}}
#if HAVE_MPI
#include <opm/grid/cpgrid/Entity.hpp>
#include <opm/grid/cpgrid/Indexsets.hpp>
namespace Dune {
namespace cpgrid {
namespace mover
{
template<class T>
class MoveBuffer
{
friend class Dune::cpgrid::CpGridData;
public:
void read(T& data)
{
data=buffer_[index_++];
}
void write(const T& data)
{
buffer_[index_++]=data;
}
void reset()
{
index_=0;
}
void resize(std::size_t size)
{
buffer_.resize(size);
index_=0;
}
private:
std::vector<T> buffer_;
typename std::vector<T>::size_type index_;
};
template<class DataHandle,int codim>
struct Mover
{
};
template<class DataHandle>
struct BaseMover
{
explicit BaseMover(DataHandle& data)
: data_(data)
{}
template<class E>
void moveData(const E& from, const E& to)
{
std::size_t size=data_.size(from);
buffer.resize(size);
data_.gather(buffer, from);
buffer.reset();
data_.scatter(buffer, to, size);
}
DataHandle& data_;
MoveBuffer<typename DataHandle::DataType> buffer;
};
template<class DataHandle>
struct Mover<DataHandle,0> : public BaseMover<DataHandle>
{
Mover(DataHandle& data, CpGridData* gatherView,
CpGridData* scatterView)
: BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
{}
void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
{
Entity<0> from_entity=Entity<0>(*gatherView_, from_cell_index, true);
Entity<0> to_entity=Entity<0>(*scatterView_, to_cell_index, true);
this->moveData(from_entity, to_entity);
}
CpGridData* gatherView_;
CpGridData* scatterView_;
};
template<class DataHandle>
struct Mover<DataHandle,1> : public BaseMover<DataHandle>
{
Mover(DataHandle& data, CpGridData* gatherView,
CpGridData* scatterView)
: BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
{}
void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
{
typedef typename OrientedEntityTable<0,1>::row_type row_type;
EntityRep<0> from_cell=EntityRep<0>(from_cell_index, true);
EntityRep<0> to_cell=EntityRep<0>(to_cell_index, true);
const OrientedEntityTable<0,1>& table = gatherView_->cell_to_face_;
row_type from_faces=table.operator[](from_cell);
row_type to_faces=scatterView_->cell_to_face_[to_cell];
for(int i=0; i<from_faces.size(); ++i)
this->moveData(from_faces[i], to_faces[i]);
}
CpGridData *gatherView_;
CpGridData *scatterView_;
};
template<class DataHandle>
struct Mover<DataHandle,3> : public BaseMover<DataHandle>
{
Mover(DataHandle& data, CpGridData* gatherView,
CpGridData* scatterView)
: BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
{}
void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
{
const std::array<int,8>& from_cell_points=
gatherView_->cell_to_point_[from_cell_index];
const std::array<int,8>& to_cell_points=
scatterView_->cell_to_point_[to_cell_index];
for(std::size_t i=0; i<8; ++i)
{
this->moveData(Entity<3>(*gatherView_, from_cell_points[i], true),
Entity<3>(*scatterView_, to_cell_points[i], true));
}
}
CpGridData* gatherView_;
CpGridData* scatterView_;
};
} // end mover namespace
template<class DataHandle>
void CpGridData::scatterData(DataHandle& data, const CpGridData* global_data,
const CpGridData* distributed_data, const InterfaceMap& cell_inf,
const InterfaceMap& point_inf)
{
#if HAVE_MPI
if(data.contains(3,0))
{
Entity2IndexDataHandle<DataHandle, 0> data_wrapper(*global_data, *distributed_data, data);
communicateCodim<0>(data_wrapper, ForwardCommunication, cell_inf);
}
if(data.contains(3,3))
{
Entity2IndexDataHandle<DataHandle, 3> data_wrapper(*global_data, *distributed_data, data);
communicateCodim<3>(data_wrapper, ForwardCommunication, point_inf);
}
#endif
}
template<int codim, class DataHandle>
void CpGridData::scatterCodimData(DataHandle& data, CpGridData* global_data,
CpGridData* distributed_data)
{
CpGridData *gather_view, *scatter_view;
gather_view=global_data;
scatter_view=distributed_data;
mover::Mover<DataHandle,codim> mover(data, gather_view, scatter_view);
for(auto index=distributed_data->cellIndexSet().begin(),
end = distributed_data->cellIndexSet().end();
index!=end; ++index)
{
std::size_t from=index->global();
std::size_t to=index->local();
mover(from,to);
}
}
namespace
{
template<int codim, class T, class F>
void visitInterior(CpGridData& distributed_data, T begin, T endit, F& func)
{
for(T it=begin; it!=endit; ++it)
{
Entity<codim> entity(distributed_data, it-begin, true);
PartitionType pt = entity.partitionType();
if(pt==Dune::InteriorEntity)
{
func(*it, entity);
}
}
}
template<class DataHandle>
struct GlobalIndexSizeGatherer
{
GlobalIndexSizeGatherer(DataHandle& data_,
std::vector<int>& ownedGlobalIndices_,
std::vector<int>& ownedSizes_)
: data(data_), ownedGlobalIndices(ownedGlobalIndices_), ownedSizes(ownedSizes_)
{}
template<class T, class E>
void operator()(T& i, E& entity)
{
ownedGlobalIndices.push_back(i);
ownedSizes.push_back(data.size(entity));
}
DataHandle& data;
std::vector<int>& ownedGlobalIndices;
std::vector<int>& ownedSizes;
};
template<class DataHandle>
struct DataGatherer
{
DataGatherer(mover::MoveBuffer<typename DataHandle::DataType>& buffer_,
DataHandle& data_)
: buffer(buffer_), data(data_)
{}
template<class T, class E>
void operator()(T& /* it */, E& entity)
{
data.gather(buffer, entity);
}
mover::MoveBuffer<typename DataHandle::DataType>& buffer;
DataHandle& data;
};
}
template<class DataHandle>
void CpGridData::gatherData(DataHandle& data, CpGridData* global_data,
CpGridData* distributed_data)
{
#if HAVE_MPI
if(data.contains(3,0))
gatherCodimData<0>(data, global_data, distributed_data);
if(data.contains(3,3))
gatherCodimData<3>(data, global_data, distributed_data);
#endif
}
template<int codim, class DataHandle>
void CpGridData::gatherCodimData(DataHandle& data, CpGridData* global_data,
CpGridData* distributed_data)
{
#if HAVE_MPI
// Get the mapping to global index from the global id set
const std::vector<int>& mapping =
distributed_data->global_id_set_->getMapping<codim>();
// Get the global indices and data size for the entities whose data is
// to be sent, i.e. the ones that we own.
std::vector<int> owned_global_indices;
std::vector<int> owned_sizes;
owned_global_indices.reserve(mapping.size());
owned_sizes.reserve(mapping.size());
GlobalIndexSizeGatherer<DataHandle> gisg(data, owned_global_indices, owned_sizes);
visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gisg);
// communicate the number of indices that each processor sends
int no_indices=owned_sizes.size();
// We will take the address of the first elemet for MPI_Allgather below.
// Make sure the containers have such an element.
if ( owned_global_indices.empty() )
owned_global_indices.resize(1);
if ( owned_sizes.empty() )
owned_sizes.resize(1);
std::vector<int> no_indices_to_recv(distributed_data->ccobj_.size());
distributed_data->ccobj_.allgather(&no_indices, 1, &(no_indices_to_recv[0]));
// compute size of the vector capable for receiving all indices
// and allgather the global indices and the sizes.
// calculate displacements
std::vector<int> displ(distributed_data->ccobj_.size()+1, 0);
std::transform(displ.begin(), displ.end()-1, no_indices_to_recv.begin(), displ.begin()+1,
std::plus<int>());
int global_size=displ[displ.size()-1];//+no_indices_to_recv[displ.size()-1];
std::vector<int> global_indices(global_size);
std::vector<int> global_sizes(global_size);
MPI_Allgatherv(&(owned_global_indices[0]), no_indices, MPITraits<int>::getType(),
&(global_indices[0]), &(no_indices_to_recv[0]), &(displ[0]),
MPITraits<int>::getType(),
distributed_data->ccobj_);
MPI_Allgatherv(&(owned_sizes[0]), no_indices, MPITraits<int>::getType(),
&(global_sizes[0]), &(no_indices_to_recv[0]), &(displ[0]),
MPITraits<int>::getType(),
distributed_data->ccobj_);
std::vector<int>().swap(owned_global_indices); // free data for reuse.
// Compute the number of data items to send
std::vector<int> no_data_send(distributed_data->ccobj_.size());
for(typename std::vector<int>::iterator begin=no_data_send.begin(),
i=begin, end=no_data_send.end(); i!=end; ++i)
*i = std::accumulate(global_sizes.begin()+displ[i-begin],
global_sizes.begin()+displ[i-begin+1], std::size_t());
// free at least some memory that can be reused.
std::vector<int>().swap(owned_sizes);
// compute the displacements for receiving with allgatherv
displ[0]=0;
std::transform(displ.begin(), displ.end()-1, no_data_send.begin(), displ.begin()+1,
std::plus<std::size_t>());
// Compute the number of data items we will receive
int no_data_recv = displ[displ.size()-1];//+global_sizes[displ.size()-1];
// Collect the data to send, gather it
mover::MoveBuffer<typename DataHandle::DataType> local_data_buffer, global_data_buffer;
if ( no_data_send[distributed_data->ccobj_.rank()] )
{
local_data_buffer.resize(no_data_send[distributed_data->ccobj_.rank()]);
}
else
{
local_data_buffer.resize(1);
}
global_data_buffer.resize(no_data_recv);
DataGatherer<DataHandle> gatherer(local_data_buffer, data);
visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gatherer);
MPI_Allgatherv(&(local_data_buffer.buffer_[0]), no_data_send[distributed_data->ccobj_.rank()],
MPITraits<typename DataHandle::DataType>::getType(),
&(global_data_buffer.buffer_[0]), &(no_data_send[0]), &(displ[0]),
MPITraits<typename DataHandle::DataType>::getType(),
distributed_data->ccobj_);
Entity2IndexDataHandle<DataHandle, codim> edata(*global_data, data);
int offset=0;
for(int i=0; i< codim; ++i)
offset+=global_data->size(i);
typename std::vector<int>::const_iterator s=global_sizes.begin();
for(typename std::vector<int>::const_iterator i=global_indices.begin(),
end=global_indices.end();
i!=end; ++s, ++i)
{
edata.scatter(global_data_buffer, *i-offset, *s);
}
#endif
}
} // end namespace cpgrid
} // end namespace Dune
#endif
#endif
Dune::CpGrid
[ provides Dune::Grid ]
Definition CpGrid.hpp:203
Dune::cpgrid::CpGridData
Struct that hods all the data needed to represent a Cpgrid.
Definition CpGridData.hpp:118
Dune::cpgrid::CpGridData::postAdapt
void postAdapt()
Clean up refinement/coarsening markers - set every element to the mark 0 which represents 'doing noth...
Definition CpGridData.cpp:1950
Dune::cpgrid::CpGridData::globalIdSet
const cpgrid::LevelGlobalIdSet & globalIdSet() const
Get the global index set.
Definition CpGridData.hpp:581
Dune::cpgrid::CpGridData::size
int size(GeometryType type) const
number of leaf entities per geometry type in this process
Definition CpGridData.hpp:175
Dune::cpgrid::CpGridData::logicalCartesianSize
const std::array< int, 3 > & logicalCartesianSize() const
The logical cartesian size of the grid.
Definition CpGridData.hpp:589
Dune::cpgrid::CpGridData::communicate
void communicate(DataHandle &data, InterfaceType iftype, CommunicationDirection dir)
communicate objects for all codims on a given level
Definition CpGridData.hpp:943
Dune::cpgrid::CpGridData::uniqueBoundaryIds
bool uniqueBoundaryIds() const
Is the grid currently using unique boundary ids?
Definition CpGridData.hpp:544
Dune::cpgrid::CpGridData::getChildrenLevelAndIndexList
const std::tuple< int, std::vector< int > > & getChildrenLevelAndIndexList(int elemIdx) const
Retrieves the level and child indices of a given parent cell.
Definition CpGridData.hpp:472
Dune::cpgrid::CpGridData::size
int size(int codim) const
number of leaf entities per codim in this process
Definition CpGridData.cpp:149
Dune::cpgrid::CpGridData::readEclipseFormat
void readEclipseFormat(const std::string &filename, bool periodic_extension, bool turn_normals=false, bool edge_conformal=false)
Read the Eclipse grid format ('grdecl').
Dune::cpgrid::CpGridData::globalCell
const std::vector< int > & globalCell() const
Return global_cell_ of any level grid, or the leaf grid view (in presence of refinement).
Definition CpGridData.hpp:393
Dune::cpgrid::CpGridData::Communication
CpGridDataTraits::Communication Communication
The type of the collective communication.
Definition CpGridData.hpp:618
Dune::cpgrid::CpGridData::~CpGridData
~CpGridData()
Destructor.
Definition CpGridData.cpp:102
Dune::cpgrid::CpGridData::getIJK
void getIJK(int c, std::array< int, 3 > &ijk) const
Extract Cartesian index triplet (i,j,k) of an active cell.
Definition CpGridData.cpp:1713
Dune::cpgrid::CpGridData::indexSet
const IndexSet & indexSet() const
Get the index set.
Definition CpGridData.hpp:569
Dune::cpgrid::CpGridData::getMark
int getMark(const cpgrid::Entity< 0 > &element) const
Return refinement mark for entity.
Definition CpGridData.cpp:1901
Dune::cpgrid::CpGridData::MPICommunicator
CpGridDataTraits::MPICommunicator MPICommunicator
The type of the mpi communicator.
Definition CpGridData.hpp:616
Dune::cpgrid::CpGridData::refineSingleCell
std::tuple< const std::shared_ptr< CpGridData >, const std::vector< std::array< int, 2 > >, const std::vector< std::tuple< int, std::vector< int > > >, const std::tuple< int, std::vector< int > >, const std::vector< std::array< int, 2 > >, const std::vector< std::array< int, 2 > > > refineSingleCell(const std::array< int, 3 > &cells_per_dim, const int &parent_idx) const
Refine a single cell and return a shared pointer of CpGridData type.
Definition CpGridData.cpp:1763
Dune::cpgrid::CpGridData::distributeGlobalGrid
void distributeGlobalGrid(CpGrid &grid, const CpGridData &view_data, const std::vector< int > &cell_part)
Redistribute a global grid.
Definition CpGridData.cpp:1489
Dune::cpgrid::CpGridData::levelData
const std::vector< std::shared_ptr< Dune::cpgrid::CpGridData > > & levelData() const
Add doc/or remove method and replace it with better approach.
Definition CpGridData.hpp:457
Dune::cpgrid::CpGridData::zcornData
const std::vector< double > & zcornData() const
Return the internalized zcorn copy from the grid processing, if no cells were adjusted during the min...
Definition CpGridData.hpp:562
Dune::cpgrid::CpGridData::setUniqueBoundaryIds
void setUniqueBoundaryIds(bool uids)
Set whether we want to have unique boundary ids.
Definition CpGridData.hpp:551
Dune::cpgrid::CpGridData::preAdapt
bool preAdapt()
Set mightVanish flags for elements that will be refined in the next adapt() call Need to be called af...
Definition CpGridData.cpp:1906
Dune::cpgrid::CpGridData::MAX_DATA_PER_CELL
@ MAX_DATA_PER_CELL
The maximum data items allowed per cell (DUNE < 2.5.2).
Definition CpGridData.hpp:143
Dune::cpgrid::CpGridData::getGridIdx
int getGridIdx() const
Add doc/or remove method and replace it with better approach.
Definition CpGridData.hpp:443
Dune::cpgrid::CpGridData::processEclipseFormat
void processEclipseFormat(const grdecl &input_data, std::array< std::set< std::pair< int, int > >, 2 > &nnc, bool remove_ij_boundary, bool turn_normals, bool pinchActive, double tolerance_unique_points, bool edge_conformal)
Read the Eclipse grid format ('grdecl').
Definition processEclipseFormat.cpp:438
Dune::cpgrid::CpGridData::hasNNCs
bool hasNNCs(const std::vector< int > &cellIndices) const
Check all cells selected for refinement have no NNCs (no neighbor connections).
Definition CpGridData.cpp:1733
Dune::cpgrid::CpGridData::mark
bool mark(int refCount, const cpgrid::Entity< 0 > &element)
Mark entity for refinement or coarsening.
Definition CpGridData.cpp:1884
Dune::cpgrid::CpGridData::localIdSet
const cpgrid::IdSet & localIdSet() const
Get the local index set.
Definition CpGridData.hpp:575
Dune::cpgrid::CpGridData::adapt
bool adapt()
TO DO: Documentation. Triggers the grid refinement process - Currently, returns preAdapt().
Definition CpGridData.cpp:1945
Dune::cpgrid::CpGridData::AttributeSet
CpGridDataTraits::AttributeSet AttributeSet
The type of the set of the attributes.
Definition CpGridData.hpp:622
Dune::cpgrid::CpGridData::sortedNumAquiferCells
const std::vector< int > & sortedNumAquiferCells() const
Get sorted active cell indices of numerical aquifer.
Definition CpGridData.hpp:677
Dune::cpgrid::DefaultGeometryPolicy
Definition DefaultGeometryPolicy.hpp:53
Dune::cpgrid::Entity2IndexDataHandle
Wrapper that turns a data handle suitable for dune-grid into one based on integers instead of entitie...
Definition Entity2IndexDataHandle.hpp:56
Dune::cpgrid::EntityRep
Represents an entity of a given codim, with positive or negative orientation.
Definition EntityRep.hpp:99
Dune::cpgrid::EntityVariable
A class design to hold a variable with a value for each entity of the given codimension,...
Definition EntityRep.hpp:272
Dune::cpgrid::Entity
Definition Entity.hpp:72
Dune::cpgrid::Geometry
This class encapsulates geometry for vertices, intersections, and cells.
Definition Geometry.hpp:76
Dune::cpgrid::IdSet
Definition Indexsets.hpp:201
Dune::cpgrid::IndexSet
Definition Indexsets.hpp:56
Dune::cpgrid::LevelGlobalIdSet
Definition Indexsets.hpp:371
Dune::cpgrid::OrientedEntityTable
Represents the topological relationships between sets of entities, for example cells and faces.
Definition OrientedEntityTable.hpp:139
Dune::cpgrid::PartitionTypeIndicator
Definition PartitionTypeIndicator.hpp:50
Dune::cpgrid::SignedEntityVariable
A class design to hold a variable with a value for each entity of the given codimension,...
Definition EntityRep.hpp:307
Opm::SparseTable
A SparseTable stores a table with rows of varying size as efficiently as possible.
Definition SparseTable.hpp:55
Dune::cpgrid
Copyright 2019 Equinor AS.
Definition GridPartitioning.cpp:673
Dune
The namespace Dune is the main namespace for all Dune code.
Definition CartesianIndexMapper.hpp:10
Opm
Holds the implementation of the CpGrid as a pimple.
Definition CellQuadrature.cpp:68
preprocess.h
Low-level corner-point processing routines and supporting data structures.
Dune::cpgrid::CpGridDataTraits::MPICommunicator
MPIHelper::MPICommunicator MPICommunicator
The type of the collective communication.
Definition CpGridDataTraits.hpp:56
Dune::cpgrid::CpGridDataTraits::AttributeSet
AttributeSet
The type of the set of the attributes.
Definition CpGridDataTraits.hpp:66
Dune::cpgrid::mover::Mover
Definition CpGridData.hpp:110
grdecl
Raw corner-point specification of a particular geological model.
Definition preprocess.h:56