My Project
Loading...
Searching...
No Matches
CpGridData.hpp
1//===========================================================================
2//
3// File: CpGridData.hpp
4//
5// Created: Sep 17 21:11:41 2013
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// Bård Skaflestad <bard.skaflestad@sintef.no>
9// Markus Blatt <markus@dr-blatt.de>
10// Antonella Ritorto <antonella.ritorto@opm-op.com>
11//
12// Comment: Major parts of this file originated in dune/grid/CpGrid.hpp
13// and got transfered here during refactoring for the parallelization.
14//
15// $Date$
16//
17// $Revision$
18//
19//===========================================================================
20
21/*
22 Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
23 Copyright 2009, 2010, 2013, 2022-2023 Equinor ASA.
24 Copyright 2013 Dr. Blatt - HPC-Simulation-Software & Services
25
26 This file is part of The Open Porous Media project (OPM).
27
28 OPM is free software: you can redistribute it and/or modify
29 it under the terms of the GNU General Public License as published by
30 the Free Software Foundation, either version 3 of the License, or
31 (at your option) any later version.
32
33 OPM is distributed in the hope that it will be useful,
34 but WITHOUT ANY WARRANTY; without even the implied warranty of
35 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36 GNU General Public License for more details.
37
38 You should have received a copy of the GNU General Public License
39 along with OPM. If not, see <http://www.gnu.org/licenses/>.
40*/
48#ifndef OPM_CPGRIDDATA_HEADER
49#define OPM_CPGRIDDATA_HEADER
50
51
52#include <dune/common/parallel/mpihelper.hh>
53#ifdef HAVE_DUNE_ISTL
54#include <dune/istl/owneroverlapcopy.hh>
55#endif
56
57#include <dune/common/parallel/communication.hh>
58#include <dune/common/parallel/variablesizecommunicator.hh>
59#include <dune/grid/common/gridenums.hh>
60
61#if HAVE_ECL_INPUT
62#include <opm/input/eclipse/EclipseState/Grid/NNC.hpp>
63#endif
64
66
67#include "Entity2IndexDataHandle.hpp"
68#include "CpGridDataTraits.hpp"
69//#include "DataHandleWrappers.hpp"
70//#include "GlobalIdMapping.hpp"
71#include "Geometry.hpp"
72
73#include <array>
74#include <initializer_list>
75#include <set>
76#include <vector>
77
78namespace Opm
79{
80class EclipseState;
81}
82namespace Dune
83{
84class CpGrid;
85
86namespace cpgrid
87{
88
89class IndexSet;
90class IdSet;
91class LevelGlobalIdSet;
92class PartitionTypeIndicator;
93template<int,int> class Geometry;
94template<int> class Entity;
95template<int> class EntityRep;
96}
97}
98
99void refine_and_check(const Dune::cpgrid::Geometry<3, 3>&,
100 const std::array<int, 3>&,
101 bool);
102
103void refinePatch_and_check(const std::array<int,3>&,
104 const std::array<int,3>&,
105 const std::array<int,3>&);
106
107void refinePatch_and_check(Dune::CpGrid&,
108 const std::vector<std::array<int,3>>&,
109 const std::vector<std::array<int,3>>&,
110 const std::vector<std::array<int,3>>&,
111 const std::vector<std::string>&);
112
113void check_global_refine(const Dune::CpGrid&,
114 const Dune::CpGrid&);
115
116void fieldProp_check(const Dune::CpGrid& grid, Opm::EclipseGrid eclGrid, const std::string& deck_string);
117
118void testInactiveCellsLgrs(const std::string&,
119 const std::vector<std::array<int,3>>&,
120 const std::vector<std::array<int,3>>&,
121 const std::vector<std::array<int,3>>&,
122 const std::vector<std::string>&);
123
124namespace Dune
125{
126namespace cpgrid
127{
128namespace mover
129{
130template<class T, int i> struct Mover;
131}
132
138{
139 template<class T, int i> friend struct mover::Mover;
140 friend class GlobalIdSet;
141 friend class HierarchicIterator;
142 friend class Dune::cpgrid::IndexSet;
143 friend class Dune::cpgrid::IdSet;
145
146 friend
147 void ::refine_and_check(const Dune::cpgrid::Geometry<3, 3>&,
148 const std::array<int, 3>&,
149 bool);
150 friend
151 void ::refinePatch_and_check(const std::array<int,3>&,
152 const std::array<int,3>&,
153 const std::array<int,3>&);
154
155 friend
156 void ::refinePatch_and_check(Dune::CpGrid&,
157 const std::vector<std::array<int,3>>&,
158 const std::vector<std::array<int,3>>&,
159 const std::vector<std::array<int,3>>&,
160 const std::vector<std::string>&);
161
162 friend
163 void ::check_global_refine(const Dune::CpGrid&,
164 const Dune::CpGrid&);
165 friend
166 void ::fieldProp_check(const Dune::CpGrid& grid, Opm::EclipseGrid eclGrid, const std::string& deck_string);
167 friend
168 void ::testInactiveCellsLgrs(const std::string&,
169 const std::vector<std::array<int,3>>&,
170 const std::vector<std::array<int,3>>&,
171 const std::vector<std::array<int,3>>&,
172 const std::vector<std::string>&);
173
174private:
175 CpGridData(const CpGridData& g);
176
177public:
178 enum{
179#ifndef MAX_DATA_COMMUNICATED_PER_ENTITY
187#else
192 MAX_DATA_PER_CELL = MAX_DATA_COMMUNICATED_PER_ENTITY
193#endif
194 };
195
196 CpGridData() = delete;
197
201 explicit CpGridData(MPIHelper::MPICommunicator comm, std::vector<std::shared_ptr<CpGridData>>& data);
202
203
204
206 explicit CpGridData(std::vector<std::shared_ptr<CpGridData>>& data);
208 ~CpGridData();
209
210
211
212
214 int size(int codim) const;
215
217 int size (GeometryType type) const
218 {
219 if (type.isCube()) {
220 return size(3 - type.dim());
221 } else {
222 return 0;
223 }
224 }
225
231 void readEclipseFormat(const std::string& filename, bool periodic_extension, bool turn_normals = false);
232
233#if HAVE_ECL_INPUT
242 void processEclipseFormat(const Opm::Deck& deck, bool periodic_extension, bool turn_normals = false, bool clip_z = false,
243 const std::vector<double>& poreVolume = std::vector<double>());
244
258 std::vector<std::size_t>
259 processEclipseFormat(const Opm::EclipseGrid* ecl_grid, Opm::EclipseState* ecl_state,
260 bool periodic_extension, bool turn_normals = false, bool clip_z = false,
261 bool pinchActive = true);
262#endif
263
275 void processEclipseFormat(const grdecl& input_data,
276#if HAVE_ECL_INPUT
277 Opm::EclipseState* ecl_state,
278#endif
279 std::array<std::set<std::pair<int, int>>, 2>& nnc,
280 bool remove_ij_boundary, bool turn_normals, bool pinchActive,
281 double tolerance_unique_points);
282
290 void getIJK(int c, std::array<int,3>& ijk) const
291 {
292 ijk = getIJK(global_cell_[c], logical_cartesian_size_);
293 }
294
295 int cellFace(int cell, int local_index) const
296 {
297 return cell_to_face_[cpgrid::EntityRep<0>(cell, true)][local_index].index();
298 }
299
300 int faceToCellSize(int face) const {
301 Dune::cpgrid::EntityRep<1> faceEntity(face, true);
302 return face_to_cell_[faceEntity].size();
303 }
304
311 const std::vector<int>& globalCell() const
312 {
313 return global_cell_;
314 }
315
322 std::array<int,3> getIJK(int idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const
323 {
324 // idx = k*cells_per_dim_[0]*cells_per_dim_[1] + j*cells_per_dim_[0] + i
325 // with 0<= i < cells_per_dim_[0], 0<= j < cells_per_dim_[1], 0<= k <cells_per_dim_[2].
326 assert(cells_per_dim[0]);
327 assert(cells_per_dim[1]);
328 assert(cells_per_dim[2]);
329
330 std::array<int,3> ijk = {0,0,0};
331 ijk[0] = idx_in_parent_cell % cells_per_dim[0]; idx_in_parent_cell /= cells_per_dim[0];
332 ijk[1] = idx_in_parent_cell % cells_per_dim[1];
333 ijk[2] = idx_in_parent_cell /cells_per_dim[1];
334 return ijk;
335 }
336
342 bool disjointPatches(const std::vector<std::array<int,3>>& startIJK_vec, const std::vector<std::array<int,3>>& endIJK_vec) const;
343
351 std::vector<int>
352 getPatchesCells(const std::vector<std::array<int,3>>& startIJK_vec, const std::vector<std::array<int,3>>& endIJK_vec) const;
353
356 bool hasNNCs(const std::vector<int>& cellIndices) const;
357
364 void validStartEndIJKs(const std::vector<std::array<int,3>>& startIJK_vec, const std::vector<std::array<int,3>>& endIJK_vec) const;
365
367 void checkCuboidShape(const std::vector<int>& cellIdx_vec) const;
368
374 bool patchesShareFace(const std::vector<std::array<int,3>>& startIJK_vec, const std::vector<std::array<int,3>>& endIJK_vec) const;
375
376 int sharedFaceTag(const std::vector<std::array<int,3>>& startIJK_2Patches, const std::vector<std::array<int,3>>& endIJK_2Patches) const;
377
378
391 bool mark(int refCount, const cpgrid::Entity<0>& element);
392
396 int getMark(const cpgrid::Entity<0>& element) const;
397
400 bool preAdapt();
401
403 bool adapt();
404
406 void postAdapt();
407
408private:
424 bool compatibleSubdivisions(const std::vector<std::array<int,3>>& cells_per_dim_vec,
425 const std::vector<std::array<int,3>>& startIJK_vec,
426 const std::vector<std::array<int,3>>& endIJK_vec) const;
427
428 std::array<Dune::FieldVector<double,3>,8> getReferenceRefinedCorners(int idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const;
429
437 std::array<int,3> getPatchDim(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
438
446 std::vector<int> getPatchCorners(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
447
455 std::vector<int> getPatchFaces(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
456
464 std::vector<int> getPatchCells(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
465
473 std::vector<int> getPatchBoundaryCorners(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
474
482 std::array<std::vector<int>,6> getBoundaryPatchFaces(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
483
485 std::array<std::vector<double>,3> getWidthsLengthsHeights(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
486
501 Geometry<3,3> cellifyPatch(const std::array<int,3>& startIJK, const std::array<int,3>& endIJK,
502 const std::vector<int>& patch_cells, DefaultGeometryPolicy& cellifiedPatch_geometry,
503 std::array<int,8>& cellifiedPatch_to_point,
504 std::array<int,8>& allcorners_cellifiedPatch) const;
505
506 // @brief Compute the average of array<double,3>.
507 //
508 // @param [in] vector of array<double,3>
509 // @return array<double,3> (average of the entries of the given vector).
510 std::array<double,3> getAverageArr(const std::vector<std::array<double,3>>& vec) const;
511
512public:
514 int getGridIdx() const {
515 // Not the nicest way of checking if "this" points at the leaf grid view of a mixed grid (with coarse and refined cells).
516 // 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.
517 // 2. Unfortunately, level_ is default initialized by 0. This implies, in particular, that if someone wants to check the value of
518 // "this->level_" when "this" points at the leaf grid view of a grid that has been refined, this value is - unfortunately - equal to 0.
519 // 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
520 // reason we check if child_to_parent_cells_.empty() [true for actual level 0 grid, false for the leaf grid view].
521 // --- TO BE IMPROVED ---
522 if ((level_data_ptr_ ->size() >1) && (level_ == 0) && (!child_to_parent_cells_.empty())) {
523 return level_data_ptr_->size() -1;
524 }
525 return level_;
526 }
528 const std::vector<std::shared_ptr<Dune::cpgrid::CpGridData>>& levelData() const
529 {
530 if (level_data_ptr_->empty()) {
531 OPM_THROW(std::logic_error, "Level data has not been initialized\n");
532 }
533 return *level_data_ptr_;
534 }
535
543 const std::tuple<int,std::vector<int>>& getChildrenLevelAndIndexList(int elemIdx) const {
544 return parent_to_children_cells_[elemIdx];
545 }
546
547 const std::vector<std::tuple<int,std::vector<int>>>& getParentToChildren() const {
548 return parent_to_children_cells_;
549 }
571 std::tuple< const std::shared_ptr<CpGridData>,
572 const std::vector<std::array<int,2>>, // parent_to_refined_corners(~boundary_old_to_new_corners)
573 const std::vector<std::tuple<int,std::vector<int>>>, // parent_to_children_faces (~boundary_old_to_new_faces)
574 const std::tuple<int, std::vector<int>>, // parent_to_children_cells
575 const std::vector<std::array<int,2>>, // child_to_parent_faces
576 const std::vector<std::array<int,2>>> // child_to_parent_cells
577 refineSingleCell(const std::array<int,3>& cells_per_dim, const int& parent_idx) const;
578
591 std::tuple< std::shared_ptr<CpGridData>,
592 const std::vector<std::array<int,2>>, // boundary_old_to_new_corners
593 const std::vector<std::tuple<int,std::vector<int>>>, // boundary_old_to_new_faces
594 const std::vector<std::tuple<int,std::vector<int>>>, // parent_to_children_faces
595 const std::vector<std::tuple<int,std::vector<int>>>, // parent_to_children_cell
596 const std::vector<std::array<int,2>>, // child_to_parent_faces
597 const std::vector<std::array<int,2>>> // child_to_parent_cells
598 refinePatch(const std::array<int,3>& cells_per_dim, const std::array<int,3>& startIJK, const std::array<int,3>& endIJK) const;
599
600 // @breif Compute center of an entity/element/cell in the Eclipse way:
601 // - Average of the 4 corners of the bottom face.
602 // - Average of the 4 corners of the top face.
603 // Return average of the previous computations.
604 // @param [in] int Index of a cell.
605 // @return 'eclipse centroid'
606 std::array<double,3> computeEclCentroid(const int idx) const;
607
608 // @breif Compute center of an entity/element/cell in the Eclipse way:
609 // - Average of the 4 corners of the bottom face.
610 // - Average of the 4 corners of the top face.
611 // Return average of the previous computations.
612 // @param [in] Entity<0> Entity
613 // @return 'eclipse centroid'
614 std::array<double,3> computeEclCentroid(const Entity<0>& elem) const;
615
616 // Make unique boundary ids for all intersections.
617 void computeUniqueBoundaryIds();
618
622 bool uniqueBoundaryIds() const
623 {
624 return use_unique_boundary_ids_;
625 }
626
629 void setUniqueBoundaryIds(bool uids)
630 {
631 use_unique_boundary_ids_ = uids;
632 if (use_unique_boundary_ids_ && unique_boundary_ids_.empty()) {
633 computeUniqueBoundaryIds();
634 }
635 }
636
640 const std::vector<double>& zcornData() const {
641 return zcorn;
642 }
643
644
647 const IndexSet& indexSet() const
648 {
649 return *index_set_;
650 }
651
654 {
655 return *local_id_set_;
656 }
657
660 {
661 return *global_id_set_;
662 }
663
667 const std::array<int, 3>& logicalCartesianSize() const
668 {
669 return logical_cartesian_size_;
670 }
671
675 void distributeGlobalGrid(CpGrid& grid,
676 const CpGridData& view_data,
677 const std::vector<int>& cell_part);
678
684 template<class DataHandle>
685 void communicate(DataHandle& data, InterfaceType iftype, CommunicationDirection dir);
686
687 void computeCellPartitionType();
688
689 void computePointPartitionType();
690
691 void computeCommunicationInterfaces(int noexistingPoints);
692
696 using Communication = CpGridDataTraits::Communication;
697 using CollectiveCommunication = CpGridDataTraits::CollectiveCommunication;
698
701#if HAVE_MPI
703 using Communicator = CpGridDataTraits::Communicator;
704
706 using InterfaceMap = CpGridDataTraits::InterfaceMap;
707
709 using CommunicationType = CpGridDataTraits::CommunicationType;
710
712 using ParallelIndexSet = CpGridDataTraits::ParallelIndexSet;
713
715 using RemoteIndices = CpGridDataTraits::RemoteIndices;
716
720 CommunicationType& cellCommunication()
721 {
722 return cell_comm_;
723 }
724
728 const CommunicationType& cellCommunication() const
729 {
730 return cell_comm_;
731 }
732
733 ParallelIndexSet& cellIndexSet()
734 {
735 return cellCommunication().indexSet();
736 }
737
738 const ParallelIndexSet& cellIndexSet() const
739 {
740 return cellCommunication().indexSet();
741 }
742
743 RemoteIndices& cellRemoteIndices()
744 {
745 return cellCommunication().remoteIndices();
746 }
747
748 const RemoteIndices& cellRemoteIndices() const
749 {
750 return cellCommunication().remoteIndices();
751 }
752#endif
753
755 const std::vector<int>& sortedNumAquiferCells() const
756 {
757 return aquifer_cells_;
758 }
759
760private:
761
763 void populateGlobalCellIndexSet();
764
765#if HAVE_MPI
766
772 template<class DataHandle>
773 void gatherData(DataHandle& data, CpGridData* global_view,
774 CpGridData* distributed_view);
775
776
783 template<int codim, class DataHandle>
784 void gatherCodimData(DataHandle& data, CpGridData* global_data,
785 CpGridData* distributed_data);
786
793 template<class DataHandle>
794 void scatterData(DataHandle& data, const CpGridData* global_data,
795 const CpGridData* distributed_data, const InterfaceMap& cell_inf,
796 const InterfaceMap& point_inf);
797
805 template<int codim, class DataHandle>
806 void scatterCodimData(DataHandle& data, CpGridData* global_data,
807 CpGridData* distributed_data);
808
817 template<int codim, class DataHandle>
818 void communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
819 const Interface& interface);
820
829 template<int codim, class DataHandle>
830 void communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
831 const InterfaceMap& interface);
832
833#endif
834
835 void computeGeometry(const CpGrid& grid,
836 const DefaultGeometryPolicy& globalGeometry,
837 const std::vector<int>& globalAquiferCells,
838 const OrientedEntityTable<0, 1>& globalCell2Faces,
839 DefaultGeometryPolicy& geometry,
840 std::vector<int>& aquiferCells,
841 const OrientedEntityTable<0, 1>& cell2Faces,
842 const std::vector< std::array<int,8> >& cell2Points);
843
844 // Representing the topology
856 Opm::SparseTable<int> face_to_point_;
858 std::vector< std::array<int,8> > cell_to_point_;
865 std::array<int, 3> logical_cartesian_size_{};
872 std::vector<int> global_cell_;
878 typedef FieldVector<double, 3> PointType;
882 cpgrid::EntityVariable<int, 1> unique_boundary_ids_;
884 std::unique_ptr<cpgrid::IndexSet> index_set_;
886 std::shared_ptr<const cpgrid::IdSet> local_id_set_;
888 std::shared_ptr<LevelGlobalIdSet> global_id_set_;
890 std::shared_ptr<PartitionTypeIndicator> partition_type_indicator_;
892 std::vector<int> mark_;
894 int level_{0};
896 std::vector<std::shared_ptr<CpGridData>>* level_data_ptr_;
897 // SUITABLE FOR ALL LEVELS EXCEPT FOR LEAFVIEW
899 std::vector<int> level_to_leaf_cells_; // In entry 'level cell index', we store 'leafview cell index'
// {level LGR, {child0, child1, ...}}
901 std::vector<std::tuple<int,std::vector<int>>> parent_to_children_cells_;
// {# children in x-direction, ... y-, ... z-}
903 std::array<int,3> cells_per_dim_;
904 // SUITABLE ONLY FOR LEAFVIEW
// {level, cell index in that level}
906 std::vector<std::array<int,2>> leaf_to_level_cells_;
908 std::vector<std::array<int,2>> corner_history_;
909 // SUITABLE FOR ALL LEVELS INCLUDING LEAFVIEW
// {level parent cell, parent cell index}
911 std::vector<std::array<int,2>> child_to_parent_cells_;
914 std::vector<int> cell_to_idxInParentCell_;
915
916
917
919 Communication ccobj_;
920
921 // Boundary information (optional).
922 bool use_unique_boundary_ids_;
923
929 std::vector<double> zcorn;
930
932 std::vector<int> aquifer_cells_;
933
934#if HAVE_MPI
935
937 CommunicationType cell_comm_;
938
940 std::tuple<Interface,Interface,Interface,Interface,Interface> cell_interfaces_;
941 /*
942 // code deactivated, because users cannot access face indices and therefore
943 // communication on faces makes no sense!
945 std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
946 face_interfaces_;
947 */
949 std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
950 point_interfaces_;
951
952#endif
953
954 // Return the geometry vector corresponding to the given codim.
955 template <int codim>
956 const EntityVariable<Geometry<3 - codim, 3>, codim>& geomVector() const
957 {
958 return geometry_.geomVector<codim>();
959 }
960
961 friend class Dune::CpGrid;
962 template<int> friend class Entity;
963 template<int> friend class EntityRep;
964 friend class Intersection;
965 friend class PartitionTypeIndicator;
966};
967
968
969
970#if HAVE_MPI
971
972namespace
973{
978template<class T>
979T& getInterface(InterfaceType iftype,
980 std::tuple<T,T,T,T,T>& interfaces)
981{
982 switch(iftype)
983 {
984 case 0:
985 return std::get<0>(interfaces);
986 case 1:
987 return std::get<1>(interfaces);
988 case 2:
989 return std::get<2>(interfaces);
990 case 3:
991 return std::get<3>(interfaces);
992 case 4:
993 return std::get<4>(interfaces);
994 }
995 OPM_THROW(std::runtime_error, "Invalid Interface type was used during communication");
996}
997
998} // end unnamed namespace
999
1000template<int codim, class DataHandle>
1001void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
1002 const Interface& interface)
1003{
1004 this->template communicateCodim<codim>(data, dir, interface.interfaces());
1005}
1006
1007template<int codim, class DataHandle>
1008void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data_wrapper, CommunicationDirection dir,
1009 const InterfaceMap& interface)
1010{
1011 Communicator comm(ccobj_, interface);
1012
1013 if(dir==ForwardCommunication)
1014 comm.forward(data_wrapper);
1015 else
1016 comm.backward(data_wrapper);
1017}
1018#endif
1019
1020template<class DataHandle>
1021void CpGridData::communicate(DataHandle& data, InterfaceType iftype,
1022 CommunicationDirection dir)
1023{
1024#if HAVE_MPI
1025 if(data.contains(3,0))
1026 {
1027 Entity2IndexDataHandle<DataHandle, 0> data_wrapper(*this, data);
1028 communicateCodim<0>(data_wrapper, dir, getInterface(iftype, cell_interfaces_));
1029 }
1030 if(data.contains(3,3))
1031 {
1032 Entity2IndexDataHandle<DataHandle, 3> data_wrapper(*this, data);
1033 communicateCodim<3>(data_wrapper, dir, getInterface(iftype, point_interfaces_));
1034 }
1035#else
1036 // Suppress warnings for unused arguments.
1037 (void) data;
1038 (void) iftype;
1039 (void) dir;
1040#endif
1041}
1042}}
1043
1044#if HAVE_MPI
1045#include <opm/grid/cpgrid/Entity.hpp>
1046#include <opm/grid/cpgrid/Indexsets.hpp>
1047
1048namespace Dune {
1049namespace cpgrid {
1050
1051namespace mover
1052{
1053template<class T>
1054class MoveBuffer
1055{
1056 friend class Dune::cpgrid::CpGridData;
1057public:
1058 void read(T& data)
1059 {
1060 data=buffer_[index_++];
1061 }
1062 void write(const T& data)
1063 {
1064 buffer_[index_++]=data;
1065 }
1066 void reset()
1067 {
1068 index_=0;
1069 }
1070 void resize(std::size_t size)
1071 {
1072 buffer_.resize(size);
1073 index_=0;
1074 }
1075private:
1076 std::vector<T> buffer_;
1077 typename std::vector<T>::size_type index_;
1078};
1079template<class DataHandle,int codim>
1080struct Mover
1081{
1082};
1083
1084template<class DataHandle>
1085struct BaseMover
1086{
1087 explicit BaseMover(DataHandle& data)
1088 : data_(data)
1089 {}
1090 template<class E>
1091 void moveData(const E& from, const E& to)
1092 {
1093 std::size_t size=data_.size(from);
1094 buffer.resize(size);
1095 data_.gather(buffer, from);
1096 buffer.reset();
1097 data_.scatter(buffer, to, size);
1098 }
1099 DataHandle& data_;
1100 MoveBuffer<typename DataHandle::DataType> buffer;
1101};
1102
1103
1104template<class DataHandle>
1105struct Mover<DataHandle,0> : public BaseMover<DataHandle>
1106{
1107 Mover(DataHandle& data, CpGridData* gatherView,
1108 CpGridData* scatterView)
1109 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
1110 {}
1111
1112 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
1113 {
1114 Entity<0> from_entity=Entity<0>(*gatherView_, from_cell_index, true);
1115 Entity<0> to_entity=Entity<0>(*scatterView_, to_cell_index, true);
1116 this->moveData(from_entity, to_entity);
1117 }
1118 CpGridData* gatherView_;
1119 CpGridData* scatterView_;
1120};
1121
1122template<class DataHandle>
1123struct Mover<DataHandle,1> : public BaseMover<DataHandle>
1124{
1125 Mover(DataHandle& data, CpGridData* gatherView,
1126 CpGridData* scatterView)
1127 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
1128 {}
1129
1130 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
1131 {
1132 typedef typename OrientedEntityTable<0,1>::row_type row_type;
1133 EntityRep<0> from_cell=EntityRep<0>(from_cell_index, true);
1134 EntityRep<0> to_cell=EntityRep<0>(to_cell_index, true);
1135 const OrientedEntityTable<0,1>& table = gatherView_->cell_to_face_;
1136 row_type from_faces=table.operator[](from_cell);
1137 row_type to_faces=scatterView_->cell_to_face_[to_cell];
1138
1139 for(int i=0; i<from_faces.size(); ++i)
1140 this->moveData(from_faces[i], to_faces[i]);
1141 }
1142 CpGridData *gatherView_;
1143 CpGridData *scatterView_;
1144};
1145
1146template<class DataHandle>
1147struct Mover<DataHandle,3> : public BaseMover<DataHandle>
1148{
1149 Mover(DataHandle& data, CpGridData* gatherView,
1150 CpGridData* scatterView)
1151 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
1152 {}
1153 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
1154 {
1155 const std::array<int,8>& from_cell_points=
1156 gatherView_->cell_to_point_[from_cell_index];
1157 const std::array<int,8>& to_cell_points=
1158 scatterView_->cell_to_point_[to_cell_index];
1159 for(std::size_t i=0; i<8; ++i)
1160 {
1161 this->moveData(Entity<3>(*gatherView_, from_cell_points[i], true),
1162 Entity<3>(*scatterView_, to_cell_points[i], true));
1163 }
1164 }
1165 CpGridData* gatherView_;
1166 CpGridData* scatterView_;
1167};
1168
1169} // end mover namespace
1170
1171template<class DataHandle>
1172void CpGridData::scatterData(DataHandle& data, const CpGridData* global_data,
1173 const CpGridData* distributed_data, const InterfaceMap& cell_inf,
1174 const InterfaceMap& point_inf)
1175{
1176#if HAVE_MPI
1177 if(data.contains(3,0))
1178 {
1179 Entity2IndexDataHandle<DataHandle, 0> data_wrapper(*global_data, *distributed_data, data);
1180 communicateCodim<0>(data_wrapper, ForwardCommunication, cell_inf);
1181 }
1182 if(data.contains(3,3))
1183 {
1184 Entity2IndexDataHandle<DataHandle, 3> data_wrapper(*global_data, *distributed_data, data);
1185 communicateCodim<3>(data_wrapper, ForwardCommunication, point_inf);
1186 }
1187#endif
1188}
1189
1190template<int codim, class DataHandle>
1191void CpGridData::scatterCodimData(DataHandle& data, CpGridData* global_data,
1192 CpGridData* distributed_data)
1193{
1194 CpGridData *gather_view, *scatter_view;
1195 gather_view=global_data;
1196 scatter_view=distributed_data;
1197
1198 mover::Mover<DataHandle,codim> mover(data, gather_view, scatter_view);
1199
1200
1201 for(auto index=distributed_data->cellIndexSet().begin(),
1202 end = distributed_data->cellIndexSet().end();
1203 index!=end; ++index)
1204 {
1205 std::size_t from=index->global();
1206 std::size_t to=index->local();
1207 mover(from,to);
1208 }
1209}
1210
1211namespace
1212{
1213
1214template<int codim, class T, class F>
1215void visitInterior(CpGridData& distributed_data, T begin, T endit, F& func)
1216{
1217 for(T it=begin; it!=endit; ++it)
1218 {
1219 Entity<codim> entity(distributed_data, it-begin, true);
1220 PartitionType pt = entity.partitionType();
1221 if(pt==Dune::InteriorEntity)
1222 {
1223 func(*it, entity);
1224 }
1225 }
1226}
1227
1228template<class DataHandle>
1229struct GlobalIndexSizeGatherer
1230{
1231 GlobalIndexSizeGatherer(DataHandle& data_,
1232 std::vector<int>& ownedGlobalIndices_,
1233 std::vector<int>& ownedSizes_)
1234 : data(data_), ownedGlobalIndices(ownedGlobalIndices_), ownedSizes(ownedSizes_)
1235 {}
1236
1237 template<class T, class E>
1238 void operator()(T& i, E& entity)
1239 {
1240 ownedGlobalIndices.push_back(i);
1241 ownedSizes.push_back(data.size(entity));
1242 }
1243 DataHandle& data;
1244 std::vector<int>& ownedGlobalIndices;
1245 std::vector<int>& ownedSizes;
1246};
1247
1248template<class DataHandle>
1249struct DataGatherer
1250{
1251 DataGatherer(mover::MoveBuffer<typename DataHandle::DataType>& buffer_,
1252 DataHandle& data_)
1253 : buffer(buffer_), data(data_)
1254 {}
1255
1256 template<class T, class E>
1257 void operator()(T& /* it */, E& entity)
1258 {
1259 data.gather(buffer, entity);
1260 }
1261 mover::MoveBuffer<typename DataHandle::DataType>& buffer;
1262 DataHandle& data;
1263};
1264
1265}
1266
1267template<class DataHandle>
1268void CpGridData::gatherData(DataHandle& data, CpGridData* global_data,
1269 CpGridData* distributed_data)
1270{
1271#if HAVE_MPI
1272 if(data.contains(3,0))
1273 gatherCodimData<0>(data, global_data, distributed_data);
1274 if(data.contains(3,3))
1275 gatherCodimData<3>(data, global_data, distributed_data);
1276#endif
1277}
1278
1279template<int codim, class DataHandle>
1280void CpGridData::gatherCodimData(DataHandle& data, CpGridData* global_data,
1281 CpGridData* distributed_data)
1282{
1283#if HAVE_MPI
1284 // Get the mapping to global index from the global id set
1285 const std::vector<int>& mapping =
1286 distributed_data->global_id_set_->getMapping<codim>();
1287
1288 // Get the global indices and data size for the entities whose data is
1289 // to be sent, i.e. the ones that we own.
1290 std::vector<int> owned_global_indices;
1291 std::vector<int> owned_sizes;
1292 owned_global_indices.reserve(mapping.size());
1293 owned_sizes.reserve(mapping.size());
1294
1295 GlobalIndexSizeGatherer<DataHandle> gisg(data, owned_global_indices, owned_sizes);
1296 visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gisg);
1297
1298 // communicate the number of indices that each processor sends
1299 int no_indices=owned_sizes.size();
1300 // We will take the address of the first elemet for MPI_Allgather below.
1301 // Make sure the containers have such an element.
1302 if ( owned_global_indices.empty() )
1303 owned_global_indices.resize(1);
1304 if ( owned_sizes.empty() )
1305 owned_sizes.resize(1);
1306 std::vector<int> no_indices_to_recv(distributed_data->ccobj_.size());
1307 distributed_data->ccobj_.allgather(&no_indices, 1, &(no_indices_to_recv[0]));
1308 // compute size of the vector capable for receiving all indices
1309 // and allgather the global indices and the sizes.
1310 // calculate displacements
1311 std::vector<int> displ(distributed_data->ccobj_.size()+1, 0);
1312 std::transform(displ.begin(), displ.end()-1, no_indices_to_recv.begin(), displ.begin()+1,
1313 std::plus<int>());
1314 int global_size=displ[displ.size()-1];//+no_indices_to_recv[displ.size()-1];
1315 std::vector<int> global_indices(global_size);
1316 std::vector<int> global_sizes(global_size);
1317 MPI_Allgatherv(&(owned_global_indices[0]), no_indices, MPITraits<int>::getType(),
1318 &(global_indices[0]), &(no_indices_to_recv[0]), &(displ[0]),
1319 MPITraits<int>::getType(),
1320 distributed_data->ccobj_);
1321 MPI_Allgatherv(&(owned_sizes[0]), no_indices, MPITraits<int>::getType(),
1322 &(global_sizes[0]), &(no_indices_to_recv[0]), &(displ[0]),
1323 MPITraits<int>::getType(),
1324 distributed_data->ccobj_);
1325 std::vector<int>().swap(owned_global_indices); // free data for reuse.
1326 // Compute the number of data items to send
1327 std::vector<int> no_data_send(distributed_data->ccobj_.size());
1328 for(typename std::vector<int>::iterator begin=no_data_send.begin(),
1329 i=begin, end=no_data_send.end(); i!=end; ++i)
1330 *i = std::accumulate(global_sizes.begin()+displ[i-begin],
1331 global_sizes.begin()+displ[i-begin+1], std::size_t());
1332 // free at least some memory that can be reused.
1333 std::vector<int>().swap(owned_sizes);
1334 // compute the displacements for receiving with allgatherv
1335 displ[0]=0;
1336 std::transform(displ.begin(), displ.end()-1, no_data_send.begin(), displ.begin()+1,
1337 std::plus<std::size_t>());
1338 // Compute the number of data items we will receive
1339 int no_data_recv = displ[displ.size()-1];//+global_sizes[displ.size()-1];
1340
1341 // Collect the data to send, gather it
1342 mover::MoveBuffer<typename DataHandle::DataType> local_data_buffer, global_data_buffer;
1343 if ( no_data_send[distributed_data->ccobj_.rank()] )
1344 {
1345 local_data_buffer.resize(no_data_send[distributed_data->ccobj_.rank()]);
1346 }
1347 else
1348 {
1349 local_data_buffer.resize(1);
1350 }
1351 global_data_buffer.resize(no_data_recv);
1352
1353 DataGatherer<DataHandle> gatherer(local_data_buffer, data);
1354 visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gatherer);
1355 MPI_Allgatherv(&(local_data_buffer.buffer_[0]), no_data_send[distributed_data->ccobj_.rank()],
1356 MPITraits<typename DataHandle::DataType>::getType(),
1357 &(global_data_buffer.buffer_[0]), &(no_data_send[0]), &(displ[0]),
1358 MPITraits<typename DataHandle::DataType>::getType(),
1359 distributed_data->ccobj_);
1360 Entity2IndexDataHandle<DataHandle, codim> edata(*global_data, data);
1361 int offset=0;
1362 for(int i=0; i< codim; ++i)
1363 offset+=global_data->size(i);
1364
1365 typename std::vector<int>::const_iterator s=global_sizes.begin();
1366 for(typename std::vector<int>::const_iterator i=global_indices.begin(),
1367 end=global_indices.end();
1368 i!=end; ++s, ++i)
1369 {
1370 edata.scatter(global_data_buffer, *i-offset, *s);
1371 }
1372#endif
1373}
1374
1375} // end namespace cpgrid
1376} // end namespace Dune
1377
1378#endif
1379
1380#endif
[ provides Dune::Grid ]
Definition CpGrid.hpp:198
Struct that hods all the data needed to represent a Cpgrid.
Definition CpGridData.hpp:138
void postAdapt()
Clean up refinement/coarsening markers - set every element to the mark 0 which represents 'doing noth...
Definition CpGridData.cpp:2735
const cpgrid::LevelGlobalIdSet & globalIdSet() const
Get the global index set.
Definition CpGridData.hpp:659
@ MAX_DATA_PER_CELL
The maximum data items allowed per cell (DUNE < 2.5.2)
Definition CpGridData.hpp:186
std::vector< int > getPatchesCells(const std::vector< std::array< int, 3 > > &startIJK_vec, const std::vector< std::array< int, 3 > > &endIJK_vec) const
Compute cell indices of selected patches of cells (Cartesian grid required).
Definition CpGridData.cpp:2183
int size(GeometryType type) const
number of leaf entities per geometry type in this process
Definition CpGridData.hpp:217
const std::array< int, 3 > & logicalCartesianSize() const
The logical cartesian size of the grid.
Definition CpGridData.hpp:667
void communicate(DataHandle &data, InterfaceType iftype, CommunicationDirection dir)
communicate objects for all codims on a given level
Definition CpGridData.hpp:1021
bool uniqueBoundaryIds() const
Is the grid currently using unique boundary ids?
Definition CpGridData.hpp:622
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:543
std::array< int, 3 > getIJK(int idx_in_parent_cell, const std::array< int, 3 > &cells_per_dim) const
Extract Cartesian index triplet (i,j,k) given an index between 0 and NXxNYxNZ -1 where NX,...
Definition CpGridData.hpp:322
bool patchesShareFace(const std::vector< std::array< int, 3 > > &startIJK_vec, const std::vector< std::array< int, 3 > > &endIJK_vec) const
Determine if a finite amount of patches (of cells) share a face.
Definition CpGridData.cpp:2020
int size(int codim) const
number of leaf entities per codim in this process
Definition CpGridData.cpp:147
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:311
void readEclipseFormat(const std::string &filename, bool periodic_extension, bool turn_normals=false)
Read the Eclipse grid format ('grdecl').
CpGridDataTraits::Communication Communication
The type of the collective communication.
Definition CpGridData.hpp:696
~CpGridData()
Destructor.
Definition CpGridData.cpp:100
void getIJK(int c, std::array< int, 3 > &ijk) const
Extract Cartesian index triplet (i,j,k) of an active cell.
Definition CpGridData.hpp:290
const IndexSet & indexSet() const
Get the index set.
Definition CpGridData.hpp:647
std::tuple< std::shared_ptr< CpGridData >, const std::vector< std::array< int, 2 > >, const std::vector< std::tuple< int, std::vector< int > > >, const std::vector< std::tuple< int, std::vector< int > > >, const std::vector< std::tuple< int, std::vector< int > > >, const std::vector< std::array< int, 2 > >, const std::vector< std::array< int, 2 > > > refinePatch(const std::array< int, 3 > &cells_per_dim, const std::array< int, 3 > &startIJK, const std::array< int, 3 > &endIJK) const
Refine a (connected block-shaped) patch of cells.
Definition CpGridData.cpp:2475
int getMark(const cpgrid::Entity< 0 > &element) const
Return refinement mark for entity.
Definition CpGridData.cpp:2709
CpGridDataTraits::MPICommunicator MPICommunicator
The type of the mpi communicator.
Definition CpGridData.hpp:694
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)
Read the Eclipse grid format ('grdecl').
Definition processEclipseFormat.cpp:409
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:2344
void distributeGlobalGrid(CpGrid &grid, const CpGridData &view_data, const std::vector< int > &cell_part)
Redistribute a global grid.
Definition CpGridData.cpp:1512
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:528
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:640
void setUniqueBoundaryIds(bool uids)
Set whether we want to have unique boundary ids.
Definition CpGridData.hpp:629
bool preAdapt()
Set mightVanish flags for elements that will be refined in the next adapt() call Need to be called af...
Definition CpGridData.cpp:2714
int getGridIdx() const
Add doc/or remove method and replace it with better approach.
Definition CpGridData.hpp:514
bool hasNNCs(const std::vector< int > &cellIndices) const
Check all cells selected for refinement have no NNCs (no neighbor connections).
Definition CpGridData.cpp:2194
bool mark(int refCount, const cpgrid::Entity< 0 > &element)
Mark entity for refinement or coarsening.
Definition CpGridData.cpp:2692
const cpgrid::IdSet & localIdSet() const
Get the local index set.
Definition CpGridData.hpp:653
bool adapt()
TO DO: Documentation. Triggers the grid refinement process - Currently, returns preAdapt()
Definition CpGridData.cpp:2730
void checkCuboidShape(const std::vector< int > &cellIdx_vec) const
Check that every cell to be refined has cuboid shape.
Definition CpGridData.cpp:1824
bool disjointPatches(const std::vector< std::array< int, 3 > > &startIJK_vec, const std::vector< std::array< int, 3 > > &endIJK_vec) const
Determine if a finite amount of patches (of cells) are disjoint, namely, they do not share any corner...
Definition CpGridData.cpp:1971
const std::vector< int > & sortedNumAquiferCells() const
Get sorted active cell indices of numerical aquifer.
Definition CpGridData.hpp:755
void validStartEndIJKs(const std::vector< std::array< int, 3 > > &startIJK_vec, const std::vector< std::array< int, 3 > > &endIJK_vec) const
Check startIJK and endIJK of each patch of cells to be refined are valid, i.e.
Definition CpGridData.cpp:2218
Definition DefaultGeometryPolicy.hpp:53
const EntityVariable< cpgrid::Geometry< 3 - codim, 3 >, codim > & geomVector() const
Definition DefaultGeometryPolicy.hpp:86
Wrapper that turns a data handle suitable for dune-grid into one based on integers instead of entitie...
Definition Entity2IndexDataHandle.hpp:56
Represents an entity of a given codim, with positive or negative orientation.
Definition EntityRep.hpp:99
A class design to hold a variable with a value for each entity of the given codimension,...
Definition EntityRep.hpp:272
Definition Entity.hpp:71
This class encapsulates geometry for vertices, intersections, and cells.
Definition Geometry.hpp:75
The global id set for Dune.
Definition Indexsets.hpp:450
Only needs to provide interface for doing nothing.
Definition Iterators.hpp:108
Definition Indexsets.hpp:199
Definition Indexsets.hpp:56
Definition Indexsets.hpp:350
Represents the topological relationships between sets of entities, for example cells and faces.
Definition OrientedEntityTable.hpp:139
int size() const
Returns the number of rows in the table.
Definition SparseTable.hpp:121
A class design to hold a variable with a value for each entity of the given codimension,...
Definition EntityRep.hpp:307
A SparseTable stores a table with rows of varying size as efficiently as possible.
Definition SparseTable.hpp:55
The namespace Dune is the main namespace for all Dune code.
Definition CartesianIndexMapper.hpp:10
Holds the implementation of the CpGrid as a pimple.
Definition CellQuadrature.cpp:68
Low-level corner-point processing routines and supporting data structures.
MPIHelper::MPICommunicator MPICommunicator
The type of the collective communication.
Definition CpGridDataTraits.hpp:56
AttributeSet
The type of the set of the attributes.
Definition CpGridDataTraits.hpp:66
Definition CpGridData.hpp:130
Raw corner-point specification of a particular geological model.
Definition preprocess.h:56