47 #ifndef OPM_CPGRIDDATA_HEADER
48 #define OPM_CPGRIDDATA_HEADER
51 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
53 #include <dune/common/version.hh>
54 #include <dune/common/parallel/mpihelper.hh>
56 #include <dune/istl/owneroverlapcopy.hh>
59 #include <dune/common/parallel/collectivecommunication.hh>
60 #include <dune/common/parallel/indexset.hh>
61 #include <dune/common/parallel/interface.hh>
62 #include <dune/common/parallel/plocalindex.hh>
63 #include <dune/common/parallel/variablesizecommunicator.hh>
64 #include <dune/grid/common/gridenums.hh>
66 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
72 #include "OrientedEntityTable.hpp"
73 #include "DefaultGeometryPolicy.hpp"
76 #include <opm/grid/utility/OpmParserIncludes.hpp>
78 #include "Entity2IndexDataHandle.hpp"
79 #include "GlobalIdMapping.hpp"
91 class PartitionTypeIndicator;
98 template<
class T,
int i>
struct Mover;
127 int size(
int codim)
const;
130 int size (GeometryType type)
const
133 return size(3 - type.dim());
153 void readEclipseFormat(
const std::string& filename,
bool periodic_extension,
bool turn_normals =
false);
156 void processEclipseFormat(
const Opm::Deck& deck,
bool periodic_extension,
bool turn_normals =
false,
bool clip_z =
false,
165 const std::vector<double>& poreVolume = std::vector<double>());
175 void processEclipseFormat(
const Opm::EclipseGrid& ecl_grid,
bool periodic_extension,
bool turn_normals =
false,
bool clip_z =
false,
176 const std::vector<double>& poreVolume = std::vector<double>());
194 void getIJK(
int c, std::array<int,3>& ijk)
const
196 int gc = global_cell_[c];
197 ijk[0] = gc % logical_cartesian_size_[0]; gc /= logical_cartesian_size_[0];
198 ijk[1] = gc % logical_cartesian_size_[1];
199 ijk[2] = gc / logical_cartesian_size_[1];
203 void computeUniqueBoundaryIds();
210 return use_unique_boundary_ids_;
217 use_unique_boundary_ids_ = uids;
218 if (use_unique_boundary_ids_ && unique_boundary_ids_.empty()) {
219 computeUniqueBoundaryIds();
242 return logical_cartesian_size_;
250 const std::vector<int>& cell_part,
258 template<
class DataHandle>
259 void communicate(DataHandle& data, InterfaceType iftype, CommunicationDirection dir);
270 template<
class DataHandle>
271 void gatherData(DataHandle& data,
CpGridData* global_view,
281 template<
int codim,
class DataHandle>
282 void gatherCodimData(DataHandle& data,
CpGridData* global_data,
291 template<
class DataHandle>
292 void scatterData(DataHandle& data,
CpGridData* global_data,
302 template<
int codim,
class DataHandle>
303 void scatterCodimData(DataHandle& data,
CpGridData* global_data,
307 typedef VariableSizeCommunicator<>::InterfaceMap InterfaceMap;
317 template<
int codim,
class DataHandle>
318 void communicateCodim(DataHandle& data, CommunicationDirection dir,
319 const Interface& interface);
329 template<
int codim,
class DataHandle>
330 void communicateCodim(DataHandle& data, CommunicationDirection dir,
331 const InterfaceMap& interface);
347 std::vector< array<int,8> > cell_to_point_;
354 std::array<int, 3> logical_cartesian_size_;
361 std::vector<int> global_cell_;
367 typedef FieldVector<double, 3> PointType;
371 std::vector<PointType> allcorners_;
384 typedef MPIHelper::MPICommunicator MPICommunicator;
385 typedef Dune::CollectiveCommunication<MPICommunicator> CollectiveCommunication;
387 CollectiveCommunication ccobj_;
390 bool use_unique_boundary_ids_;
397 std::vector<double> zcorn;
399 #ifdef HAVE_DUNE_ISTL
400 typedef Dune::OwnerOverlapCopyAttributeSet::AttributeSet AttributeSet;
402 enum AttributeSet{owner, overlap, copy};
409 typedef Dune::ParallelIndexSet<int,ParallelLocalIndex<AttributeSet>,512> ParallelIndexSet;
412 ParallelIndexSet cell_indexset_;
415 typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
418 RemoteIndices cell_remote_indices_;
421 std::tuple<Interface,Interface,Interface,Interface,Interface> cell_interfaces_;
430 std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
433 InterfaceMap cell_gather_scatter_interface;
435 InterfaceMap point_gather_scatter_interface;
441 const EntityVariable<Geometry<3 - codim, 3>, codim>& geomVector()
const
447 template<
int>
friend class Entity;
448 template<
int>
friend class EntityRep;
449 template<
int>
friend class EntityPointer;
450 friend class Intersection;
451 friend class PartitionTypeIndicator;
463 T& getInterface(InterfaceType iftype,
464 std::tuple<T,T,T,T,T>& interfaces)
469 return std::get<0>(interfaces);
471 return std::get<1>(interfaces);
473 return std::get<2>(interfaces);
475 return std::get<3>(interfaces);
477 return std::get<4>(interfaces);
479 OPM_THROW(std::runtime_error,
"Invalid Interface type was used during communication");
484 template<
int codim,
class DataHandle>
485 void CpGridData::communicateCodim(DataHandle& data, CommunicationDirection dir,
486 const Interface& interface)
488 this->
template communicateCodim<codim>(data, dir, interface.interfaces());
491 template<
int codim,
class DataHandle>
492 void CpGridData::communicateCodim(DataHandle& data, CommunicationDirection dir,
493 const InterfaceMap& interface)
495 if ( interface.empty() )
502 Entity2IndexDataHandle<DataHandle, codim> data_wrapper(*
this, data);
503 VariableSizeCommunicator<> comm(ccobj_, interface);
504 if(dir==ForwardCommunication)
505 comm.forward(data_wrapper);
507 comm.backward(data_wrapper);
511 template<
class DataHandle>
513 CommunicationDirection dir)
516 if(data.contains(3,0))
517 communicateCodim<0>(data, dir, getInterface(iftype, cell_interfaces_));
518 if(data.contains(3,3))
519 communicateCodim<3>(data, dir, getInterface(iftype, point_interfaces_));
530 #include <dune/grid/cpgrid/Entity.hpp>
531 #include <dune/grid/cpgrid/Indexsets.hpp>
545 data=buffer_[index_++];
547 void write(
const T& data)
549 buffer_[index_++]=data;
555 void resize(std::size_t size)
557 buffer_.resize(size);
561 std::vector<T> buffer_;
562 typename std::vector<T>::size_type index_;
564 template<
class DataHandle,
int codim>
569 template<
class DataHandle>
572 BaseMover(DataHandle& data)
576 void moveData(
const E& from,
const E& to)
578 std::size_t size=data_.size(from);
580 data_.gather(buffer, from);
582 data_.scatter(buffer, to, size);
585 MoveBuffer<typename DataHandle::DataType> buffer;
589 template<
class DataHandle>
590 struct Mover<DataHandle,0> :
public BaseMover<DataHandle>
592 Mover<DataHandle,0>(DataHandle& data, CpGridData* gatherView,
593 CpGridData* scatterView)
594 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
597 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
599 Entity<0> from_entity=Entity<0>(*gatherView_, from_cell_index,
true);
600 Entity<0> to_entity=Entity<0>(*scatterView_, to_cell_index,
true);
601 this->moveData(from_entity, to_entity);
603 CpGridData* gatherView_;
604 CpGridData* scatterView_;
607 template<
class DataHandle>
608 struct Mover<DataHandle,1> :
public BaseMover<DataHandle>
610 Mover<DataHandle,1>(DataHandle& data, CpGridData* gatherView,
611 CpGridData* scatterView)
612 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
615 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
617 typedef typename OrientedEntityTable<0,1>::row_type row_type;
618 EntityRep<0> from_cell=EntityRep<0>(from_cell_index,
true);
619 EntityRep<0> to_cell=EntityRep<0>(to_cell_index,
true);
620 OrientedEntityTable<0,1>& table = gatherView_->cell_to_face_;
621 row_type from_faces=table.operator[](from_cell);
622 row_type to_faces=scatterView_->cell_to_face_[to_cell];
624 for(
int i=0; i<from_faces.size(); ++i)
625 this->moveData(from_faces[i], to_faces[i]);
627 CpGridData *gatherView_;
628 CpGridData *scatterView_;
631 template<
class DataHandle>
632 struct Mover<DataHandle,3> :
public BaseMover<DataHandle>
634 Mover<DataHandle,3>(DataHandle& data, CpGridData* gatherView,
635 CpGridData* scatterView)
636 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
638 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
640 const std::array<int,8>& from_cell_points=
641 gatherView_->cell_to_point_[from_cell_index];
642 const std::array<int,8>& to_cell_points=
643 scatterView_->cell_to_point_[to_cell_index];
644 for(std::size_t i=0; i<8; ++i)
646 this->moveData(Entity<3>(*gatherView_, from_cell_points[i],
true),
647 Entity<3>(*scatterView_, to_cell_points[i],
true));
650 CpGridData* gatherView_;
651 CpGridData* scatterView_;
656 template<
class DataHandle>
657 void CpGridData::scatterData(DataHandle& data, CpGridData* global_data,
658 CpGridData* distributed_data)
661 if(data.contains(3,0))
662 scatterCodimData<0>(data, global_data, distributed_data);
663 if(data.contains(3,3))
664 scatterCodimData<3>(data, global_data, distributed_data);
668 template<
int codim,
class DataHandle>
669 void CpGridData::scatterCodimData(DataHandle& data, CpGridData* global_data,
670 CpGridData* distributed_data)
672 CpGridData *gather_view, *scatter_view;
673 gather_view=global_data;
674 scatter_view=distributed_data;
676 mover::Mover<DataHandle,codim> mover(data, gather_view, scatter_view);
678 typedef typename ParallelIndexSet::const_iterator Iter;
679 for(Iter index=distributed_data->cell_indexset_.begin(),
680 end = distributed_data->cell_indexset_.end();
683 std::size_t from=index->global();
684 std::size_t to=index->local();
692 template<
int codim,
class T,
class F>
693 void visitInterior(CpGridData& distributed_data, T begin, T endit, F& func)
695 for(T it=begin; it!=endit; ++it)
697 Entity<codim> entity(distributed_data, it-begin,
true);
698 PartitionType pt = entity.partitionType();
699 if(pt==Dune::InteriorEntity)
706 template<
class DataHandle>
707 struct GlobalIndexSizeGatherer
709 GlobalIndexSizeGatherer(DataHandle& data_,
710 std::vector<int>& ownedGlobalIndices_,
711 std::vector<int>& ownedSizes_)
712 : data(data_), ownedGlobalIndices(ownedGlobalIndices_), ownedSizes(ownedSizes_)
715 template<
class T,
class E>
716 void operator()(T& i, E& entity)
718 ownedGlobalIndices.push_back(i);
719 ownedSizes.push_back(data.size(entity));
722 std::vector<int>& ownedGlobalIndices;
723 std::vector<int>& ownedSizes;
726 template<
class DataHandle>
729 DataGatherer(mover::MoveBuffer<typename DataHandle::DataType>& buffer_,
731 : buffer(buffer_), data(data_)
734 template<
class T,
class E>
735 void operator()(T& , E& entity)
737 data.gather(buffer, entity);
739 mover::MoveBuffer<typename DataHandle::DataType>& buffer;
745 template<
class DataHandle>
746 void CpGridData::gatherData(DataHandle& data, CpGridData* global_data,
747 CpGridData* distributed_data)
750 if(data.contains(3,0))
751 gatherCodimData<0>(data, global_data, distributed_data);
752 if(data.contains(3,3))
753 gatherCodimData<3>(data, global_data, distributed_data);
757 template<
int codim,
class DataHandle>
758 void CpGridData::gatherCodimData(DataHandle& data, CpGridData* global_data,
759 CpGridData* distributed_data)
763 const std::vector<int>& mapping =
764 distributed_data->global_id_set_->getMapping<codim>();
768 std::vector<int> owned_global_indices;
769 std::vector<int> owned_sizes;
770 owned_global_indices.reserve(mapping.size());
771 owned_sizes.reserve(mapping.size());
773 GlobalIndexSizeGatherer<DataHandle> gisg(data, owned_global_indices, owned_sizes);
774 visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gisg);
777 int no_indices=owned_sizes.size();
780 if ( owned_global_indices.empty() )
781 owned_global_indices.resize(1);
782 if ( owned_sizes.empty() )
783 owned_sizes.resize(1);
784 std::vector<int> no_indices_to_recv(distributed_data->ccobj_.size());
785 distributed_data->ccobj_.allgather(&no_indices, 1, &(no_indices_to_recv[0]));
789 std::vector<int> displ(distributed_data->ccobj_.size()+1, 0);
790 std::transform(displ.begin(), displ.end()-1, no_indices_to_recv.begin(), displ.begin()+1,
792 int global_size=displ[displ.size()-1];
793 std::vector<int> global_indices(global_size);
794 std::vector<int> global_sizes(global_size);
795 MPI_Allgatherv(&(owned_global_indices[0]), no_indices, MPITraits<int>::getType(),
796 &(global_indices[0]), &(no_indices_to_recv[0]), &(displ[0]),
797 MPITraits<int>::getType(),
798 distributed_data->ccobj_);
799 MPI_Allgatherv(&(owned_sizes[0]), no_indices, MPITraits<int>::getType(),
800 &(global_sizes[0]), &(no_indices_to_recv[0]), &(displ[0]),
801 MPITraits<int>::getType(),
802 distributed_data->ccobj_);
803 std::vector<int>().swap(owned_global_indices);
805 std::vector<int> no_data_send(distributed_data->ccobj_.size());
806 for(
typename std::vector<int>::iterator begin=no_data_send.begin(),
807 i=begin, end=no_data_send.end(); i!=end; ++i)
808 *i = std::accumulate(global_sizes.begin()+displ[i-begin],
809 global_sizes.begin()+displ[i-begin+1], std::size_t());
811 std::vector<int>().swap(owned_sizes);
814 std::transform(displ.begin(), displ.end()-1, no_data_send.begin(), displ.begin()+1,
815 std::plus<std::size_t>());
817 int no_data_recv = displ[displ.size()-1];
820 mover::MoveBuffer<typename DataHandle::DataType> local_data_buffer, global_data_buffer;
821 if ( no_data_send[distributed_data->ccobj_.rank()] )
823 local_data_buffer.resize(no_data_send[distributed_data->ccobj_.rank()]);
827 local_data_buffer.resize(1);
829 global_data_buffer.resize(no_data_recv);
831 DataGatherer<DataHandle> gatherer(local_data_buffer, data);
832 visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gatherer);
833 MPI_Allgatherv(&(local_data_buffer.buffer_[0]), no_data_send[distributed_data->ccobj_.rank()],
834 MPITraits<typename DataHandle::DataType>::getType(),
835 &(global_data_buffer.buffer_[0]), &(no_data_send[0]), &(displ[0]),
836 MPITraits<typename DataHandle::DataType>::getType(),
837 distributed_data->ccobj_);
838 Entity2IndexDataHandle<DataHandle, codim> edata(*global_data, data);
840 for(
int i=0; i< codim; ++i)
841 offset+=global_data->size(i);
843 typename std::vector<int>::const_iterator s=global_sizes.begin();
844 for(
typename std::vector<int>::const_iterator i=global_indices.begin(),
845 end=global_indices.end();
848 edata.scatter(global_data_buffer, *i-offset, *s);
const std::array< int, 3 > & logicalCartesianSize() const
The logical cartesian size of the grid.
Definition: CpGridData.hpp:240
Definition: CpGridData.hpp:93
int size(GeometryType type) const
number of leaf entities per geometry type in this process
Definition: CpGridData.hpp:130
Definition: DefaultGeometryPolicy.hpp:54
Raw corner-point specification of a particular geological model.
Definition: preprocess.h:56
void readSintefLegacyFormat(const std::string &grid_prefix)
Read the Sintef legacy grid format ('topogeom').
Definition: readSintefLegacyFormat.cpp:66
Definition: Indexsets.hpp:192
bool uniqueBoundaryIds() const
Is the grid currently using unique boundary ids?
Definition: CpGridData.hpp:208
void readEclipseFormat(const std::string &filename, bool periodic_extension, bool turn_normals=false)
Read the Eclipse grid format ('grdecl').
[ provides Dune::Grid ]
Definition: CpGrid.hpp:213
This class encapsulates geometry for both vertices, intersections and cells.
Definition: CpGridData.hpp:92
Struct that hods all the data needed to represent a Cpgrid.
Definition: CpGridData.hpp:105
void setUniqueBoundaryIds(bool uids)
Set whether we want to have unique boundary ids.
Definition: CpGridData.hpp:215
void communicate(DataHandle &data, InterfaceType iftype, CommunicationDirection dir)
communicate objects for all codims on a given level
Definition: CpGridData.hpp:512
CpGridData()
Constructor.
Definition: CpGridData.cpp:38
const EntityVariable< cpgrid::Geometry< 3-codim, 3 >, codim > & geomVector() const
Definition: DefaultGeometryPolicy.hpp:79
void writeSintefLegacyFormat(const std::string &grid_prefix) const
Write the Sintef legacy grid format ('topogeom').
Definition: writeSintefLegacyFormat.cpp:72
Definition: PartitionTypeIndicator.hpp:49
void processEclipseFormat(const grdecl &input_data, double z_tolerance, bool remove_ij_boundary, bool turn_normals=false)
Read the Eclipse grid format ('grdecl').
Definition: processEclipseFormat.cpp:208
int size(int codim) const
number of leaf entities per codim in this process
Definition: CpGridData.cpp:132
Definition: Indexsets.hpp:250
const IndexSet & indexSet() const
Get the index set.
Definition: CpGridData.hpp:233
void getIJK(int c, std::array< int, 3 > &ijk) const
Extract Cartesian index triplet (i,j,k) of an active cell.
Definition: CpGridData.hpp:194
Low-level corner-point processing routines and supporting data structures.
Represents an entity of a given codim, with positive or negative orientation.
Definition: CpGridData.hpp:94
Definition: Indexsets.hpp:52
Definition: CpGridData.hpp:98
void distributeGlobalGrid(const CpGrid &grid, const CpGridData &view_data, const std::vector< int > &cell_part, int overlap_layers)
Redistribute a global grid.
Definition: CpGridData.cpp:518
~CpGridData()
Destructor.
Definition: CpGridData.cpp:94
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:226