diff --git a/opm/grid/CpGrid.hpp b/opm/grid/CpGrid.hpp index 732d0259a..1e7463641 100644 --- a/opm/grid/CpGrid.hpp +++ b/opm/grid/CpGrid.hpp @@ -1221,6 +1221,9 @@ namespace Dune /// of the vertex. int faceVertex(int face, int local_index) const; + + /// \brief Get maps from vertexToCell out[0] and cellToVertex out[1]; + std::array< std::vector< std::set >, 2 > vertexCell() const; /// \brief Get vertical position of cell center ("zcorn" average). /// \brief cell_index The index of the specific cell. double cellCenterDepth(int cell_index) const; diff --git a/opm/grid/MinpvProcessor.cpp b/opm/grid/MinpvProcessor.cpp index 87702181b..9e0588be9 100644 --- a/opm/grid/MinpvProcessor.cpp +++ b/opm/grid/MinpvProcessor.cpp @@ -125,7 +125,7 @@ MinpvProcessor::process(const std::vector& thickness, bool c_active = actnum.empty() || actnum[c]; bool c_thin = (thickness[c] <= z_tolerance); bool c_thin_inactive = !c_active && c_thin; - bool c_low_pv_active = pv[c] < minpvv[c] && c_active; + bool c_low_pv_active = (pv[c] < minpvv[c] && c_active) || (c_thin && c_active); if (c_low_pv_active || c_thin_inactive) { std::array cz = getCellZcorn(ii, jj, kk, zcorn); @@ -168,7 +168,7 @@ MinpvProcessor::process(const std::vector& thickness, bool active = actnum.empty() || actnum[c_below]; bool thin = (thickness[c_below] <= z_tolerance); bool thin_inactive = !active && thin; - bool low_pv_active = pv[c_below] < minpvv[c_below] && active; + bool low_pv_active = ((pv[c_below] < minpvv[c_below]) && active) || (thin && active); while ( (thin_inactive || low_pv_active) && kk_iter < dims_[2] ) @@ -225,13 +225,26 @@ MinpvProcessor::process(const std::vector& thickness, active = actnum.empty() || actnum[c_below]; thin = (thickness[c_below] <= z_tolerance); thin_inactive = (!actnum.empty() && !actnum[c_below]) && thin; - low_pv_active = pv[c_below] < minpvv[c_below] && active; + low_pv_active = ((pv[c_below] < minpvv[c_below]) && active) || (thin && active); } // create nnc if false or merge the cells if true - if (mergeMinPVCells && c_low_pv_active) { + if (mergeMinPVCells){// && c_low_pv_active) { + // try to make a topological connected grid + // in Flow this is currently called only for edge_conformal grids + // however zcorn inbetween is not modified to make zcorn sorted // Set lower k coordinates of cell below to upper cells's coordinates. // i.e fill the void using the cell below + if (kk==0 || kk_iter == dims_[2]) { + kk = kk_iter; + continue; + } + //bottom cell not active, hence no nnc is created + if (!actnum.empty() && !actnum[c_below]) { + kk = kk_iter; + continue; + } + std::array cz_below = getCellZcorn(ii, jj, kk_iter, zcorn); for (int count = 0; count < 4; ++count) { cz_below[count] = cz[count]; @@ -260,7 +273,7 @@ MinpvProcessor::process(const std::vector& thickness, auto above_active = actnum.empty() || actnum[c_above]; auto above_inactive = !actnum.empty() && !actnum[c_above]; auto above_thin = thickness[c_above] < z_tolerance; - auto above_small_pv = pv[c_above] < minpvv[c_above]; + auto above_small_pv = (pv[c_above] < minpvv[c_above]) || above_thin; if ((above_inactive && above_thin) || (above_active && above_small_pv && (!pinchNOGAP || above_thin) ) ) { @@ -268,7 +281,7 @@ MinpvProcessor::process(const std::vector& thickness, c_above = ii + dims_[0] * (jj + dims_[1] * (k_above)); above_active = actnum.empty() || actnum[c_above]; above_inactive = !actnum.empty() && !actnum[c_above]; - auto above_significant_pv = pv[c_above] > minpvv[c_above]; + auto above_significant_pv = !((pv[c_above] < minpvv[c_above]) || (thickness[c_above] < z_tolerance)); auto above_broad = thickness[c_above] > z_tolerance; // \todo if condition seems wrong and should be the negation of above? @@ -293,9 +306,12 @@ MinpvProcessor::process(const std::vector& thickness, option4ALLZero = option4ALLZero || (!permz.empty() && permz[c_above] == 0.0) || multz(c_above) == 0.0; nnc_allowed = nnc_allowed && (computeGap(cz_above, cz_below) < max_gap) && (!pinchOption4ALL || !option4ALLZero) ; + //bool + above_small_pv = (pv[c_above] < minpvv[c_above]) || (thickness[c_above] < z_tolerance); + bool below_small_pv = (pv[c_below] < minpvv[c_below]) || (thickness[c_below] < z_tolerance); if ( nnc_allowed && (actnum.empty() || (actnum[c_above] && actnum[c_below])) && - pv[c_above] > minpvv[c_above] && pv[c_below] > minpvv[c_below]) { + !(above_small_pv) && !(below_small_pv) ){ result.add_nnc(c_above, c_below); } kk = kk_iter; @@ -303,7 +319,7 @@ MinpvProcessor::process(const std::vector& thickness, } else { - if (kk < dims_[2] - 1 && (actnum.empty() || actnum[c]) && pv[c] > minpvv[c] && + if (kk < dims_[2] - 1 && (actnum.empty() || actnum[c]) && !((pv[c] < minpvv[c]) || (thickness[c] < z_tolerance)) && multz(c) != 0.0) { // Check whether there is a gap to the neighbor below whose thickness is less @@ -311,7 +327,9 @@ MinpvProcessor::process(const std::vector& thickness, int kk_below = kk + 1; int c_below = ii + dims_[0] * (jj + dims_[1] * kk_below); - if ((actnum.empty() || actnum[c_below]) && pv[c_below] > minpvv[c_below]) + if ( (actnum.empty() || actnum[c_below]) + && + !((pv[c_below] < minpvv[c_below]) || (thickness[c_below] < z_tolerance) ) ) { // Check MAX_GAP threshold std::array cz = getCellZcorn(ii, jj, kk, zcorn); diff --git a/opm/grid/common/GridPartitioning.cpp b/opm/grid/common/GridPartitioning.cpp index e789ca77c..8409d92e2 100644 --- a/opm/grid/common/GridPartitioning.cpp +++ b/opm/grid/common/GridPartitioning.cpp @@ -260,6 +260,7 @@ void addOverlapCornerCell(const CpGrid& grid, std::vector >& cell_overlap, int level) { + assert(false); bool validLevel = (level>-1) && (level <= grid.maxLevel()); const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet(); int my_index = ix.index(from); @@ -291,6 +292,7 @@ void addOverlapCornerCell(const CpGrid& grid, std::vector>& exportList, int level) { + assert(false); // Add corner cells to the overlap layer. Example of a subdomain of a 4x4 grid // with and without corner cells in the overlap is given below. Note that the // corner cell is not needed for cell centered finite volume schemes. @@ -334,6 +336,7 @@ void addOverlapLayer(const CpGrid& grid, int recursion_deps, int level) { + assert(false); bool validLevel = (level>-1) && (level <= grid.maxLevel()); const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet(); @@ -390,6 +393,7 @@ void addOverlapLayer(const CpGrid& grid, bool all, int level) { + assert(false); cell_overlap.resize(cell_part.size()); bool validLevel = (level>-1) && (level <= grid.maxLevel()); const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet(); @@ -420,6 +424,7 @@ void addOverlapLayer(const CpGrid& grid, const std::vector& cell_part, std::vector>& exportList, bool addCornerCells, + const std::array >,2>& vertex_cell_maps, int recursion_deps, int level) { @@ -447,29 +452,43 @@ void addOverlapLayer(const CpGrid& grid, cell_part, exportList, addCornerCells, + vertex_cell_maps, recursion_deps-1, level); + } else if (addCornerCells) { - // Add cells to the overlap that just share a corner with e. - auto iit2 = validLevel? iit->outside().ilevelbegin() : iit->outside().ileafbegin(); - const auto& endIit2 = validLevel? iit->outside().ilevelend() : iit->outside().ileafend(); - - for (; iit2 != endIit2; ++iit2) { - if ( iit2->neighbor() ) - { - int nb_index2 = ix.index(iit2->outside()); - if( cell_part[nb_index2]!=owner ) { - addOverlapCornerCell(grid, - owner, - e, - iit2->outside(), - cell_part, - exportList, - level); + auto vertices = vertex_cell_maps[1][index]; //could probably used the subindices but do not trust them here + for(auto vertex: vertices){ + std::set cells = vertex_cell_maps[0][vertex]; + for(auto& cell: cells){ + if(cell_part[cell] != owner){ + // cell is the neighbor to owner + exportList.emplace_back(cell, owner, AttributeSet::copy); + // do this since it is allways done probably an ERROR for + exportList.emplace_back(index, cell_part[cell], AttributeSet::copy); } } } + // Add cells to the overlap that just share a corner with e. + // auto iit2 = validLevel? iit->outside().ilevelbegin() : iit->outside().ileafbegin(); + // const auto& endIit2 = validLevel? iit->outside().ilevelend() : iit->outside().ileafend(); + + // for (; iit2 != endIit2; ++iit2) { + // if ( iit2->neighbor() ) + // { + // int nb_index2 = ix.index(iit2->outside()); + // if( cell_part[nb_index2]!=owner ) { + // addOverlapCornerCell(grid, + // owner, + // e, + // iit2->outside(), + // cell_part, + // exportList, + // level); + // } + // } + // } } } } @@ -487,6 +506,7 @@ void addOverlapLayerNoZeroTrans(const CpGrid& grid, const double* trans, int level) { + assert(false); using AttributeSet = Dune::OwnerOverlapCopyAttributeSet::AttributeSet; bool validLevel = (level>-1) && (level <= grid.maxLevel()); const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet(); @@ -571,16 +591,16 @@ int addOverlapLayer([[maybe_unused]] const CpGrid& grid, auto it = validLevel? grid.template lbegin<0>(level) : grid.template leafbegin<0>(); const auto& endIt = validLevel? grid.template lend<0>(level) : grid.template leafend<0>(); - + const std::array>,2> vertex_cell_maps = grid.vertexCell(); for (; it != endIt; ++it) { int index = ix.index(*it); auto owner = cell_part[index]; exportProcs.insert(std::make_pair(owner, 0)); - if ( trans ) { + if ( trans && false) { addOverlapLayerNoZeroTrans(grid, index, *it, owner, cell_part, exportList, addCornerCells, layers-1, trans, level); } else { - addOverlapLayer(grid, index, *it, owner, cell_part, exportList, addCornerCells, layers-1, level); + addOverlapLayer(grid, index, *it, owner, cell_part, exportList, addCornerCells, vertex_cell_maps, layers-1, level); } } // remove multiple entries diff --git a/opm/grid/common/GridPartitioning.hpp b/opm/grid/common/GridPartitioning.hpp index 84329bbe2..5a7d91746 100644 --- a/opm/grid/common/GridPartitioning.hpp +++ b/opm/grid/common/GridPartitioning.hpp @@ -118,7 +118,8 @@ namespace Dune bool addCornerCells, const double* trans, int layers = 1, - int level = -1); + int level = -1 + ); namespace cpgrid { diff --git a/opm/grid/cpgrid/CpGrid.cpp b/opm/grid/cpgrid/CpGrid.cpp index 38568e211..9848a8925 100644 --- a/opm/grid/cpgrid/CpGrid.cpp +++ b/opm/grid/cpgrid/CpGrid.cpp @@ -390,6 +390,7 @@ CpGrid::scatterGrid(EdgeWeightMethod method, comm().barrier(); // first create the overlap + auto vertex_cell_map = this->vertexCell(); auto noImportedOwner = addOverlapLayer(*this, computedCellPart, exportList, @@ -1176,6 +1177,30 @@ int CpGrid::cellFace(int cell, int local_index, int level) const : current_data_->back()->cell_to_face_[cpgrid::EntityRep<0>(cell, true)][local_index].index(); } +std::array>,2> CpGrid::vertexCell() const{ + //int level = -1; + int nc = numCells();//(level) + int nv = this->numVertices(); + std::vector> vertex_cell(nv); + std::vector> cell_vertex(nc); + for(int cell=0; cell< nc; ++cell){ //cells + int nlf = numCellFaces(cell);//level); + for(int lf=0; lf < nlf; ++lf){//local faces of cell + int face = this->cellFace(cell,lf);//level); + int nlv = numFaceVertices(face);//numFaceVertices(face,level); + for(int lv=0; lv < nlv; ++lv){//local nodes of face + int vertex = faceVertex(face,lv); + vertex_cell[vertex].insert(cell); + cell_vertex[cell].insert(vertex); + } + } + } + std::array>,2> out; + out[0] = vertex_cell; + out[1] = cell_vertex; + return out; +} + const cpgrid::OrientedEntityTable<0,1>::row_type CpGrid::cellFaceRow(int cell) const { return current_data_->back()->cell_to_face_[cpgrid::EntityRep<0>(cell, true)]; diff --git a/opm/grid/cpgrid/CpGridData.cpp b/opm/grid/cpgrid/CpGridData.cpp index c4fdabd06..68bd8f5a5 100644 --- a/opm/grid/cpgrid/CpGridData.cpp +++ b/opm/grid/cpgrid/CpGridData.cpp @@ -206,15 +206,35 @@ PartitionType getPartitionType(const PartitionTypeIndicator& p, int i, return p.getPartitionType(Entity<3>(grid, i, true)); } +// template +// int getIndex(const T* i) +// { +// return *i; +// } + + +int getIndex(std::vector::const_iterator i) +{ + return *i; +} + + int getIndex(const int* i) { - return *i; + return *i; } +// template +// int getIndex(const EntityRep* i) +// { +// return i->index(); +// } + + template int getIndex(T i) { - return i->index(); + return i->index(); } template @@ -934,7 +954,8 @@ struct AttributeDataHandle bool fixedSize() { - return true; + return false; + //return true; } std::size_t size(std::size_t i) { @@ -943,8 +964,8 @@ struct AttributeDataHandle template void gather(B& buffer, std::size_t i) { - typedef typename GetRowType::type::const_iterator RowIter; - for(RowIter f=c2e_[i].begin(), fend=c2e_[i].end(); + //typedef typename GetRowType::type::const_iterator RowIter; + for(auto f=c2e_[i].begin(), fend=c2e_[i].end(); f!=fend; ++f) { char t=getPartitionType(indicator_, *f, grid_); @@ -955,8 +976,8 @@ struct AttributeDataHandle template void scatter(B& buffer, std::size_t i, std::size_t s) { - typedef typename GetRowType::type::const_iterator RowIter; - for(RowIter f=c2e_[i].begin(), fend=c2e_[i].end(); + //typedef typename GetRowType::type::const_iterator RowIter; + for(auto f=c2e_[i].begin(), fend=c2e_[i].end(); f!=fend; ++f, --s) { std::pair rank_attr; @@ -1588,8 +1609,8 @@ void CpGridData::distributeGlobalGrid(CpGrid& grid, // Compute partition type for points computePointPartitionType(); - - computeCommunicationInterfaces(noExistingPoints); + + computeCommunicationInterfaces(noExistingPoints); #else // #if HAVE_MPI static_cast(grid); static_cast(view_data); @@ -1618,41 +1639,61 @@ void CpGridData::computePointPartitionType() // type border, then the type of the point is overwritten with border. In the other cases // we set the type of the point to the one of the face as long as the type of the point is // not border. + // ghost should never happen; + std::array type_order({InteriorEntity,OverlapEntity,BorderEntity,FrontEntity}); partition_type_indicator_->point_indicator_.resize(geometry_.geomVector<3>().size(), InteriorEntity); - // marking border entities first for(int i=0; igetFacePartitionType(i); - if (face_type == BorderEntity) + for(auto p=face_to_point_[i].begin(), + pend=face_to_point_[i].end(); p!=pend; ++p) { - // all vertices are border - for(auto p=face_to_point_[i].begin(), - pend=face_to_point_[i].end(); p!=pend; ++p) - { - partition_type_indicator_->point_indicator_[*p]=face_type; + PartitionType new_type=partition_type_indicator_->getFacePartitionType(i); + PartitionType old_type=PartitionType(partition_type_indicator_->point_indicator_[*p]); + int new_order=-1; + int old_order=-1; + for(int j=0; j < 4; ++j){ + if(new_type == type_order[j]){ + new_order = j; + } + if(old_type == type_order[j]){ + old_order = j; + } } - } - } - - for(int i=0; igetFacePartitionType(i); - - if (new_type != InteriorEntity && new_type != BorderEntity) - { - for(auto p=face_to_point_[i].begin(), - pend=face_to_point_[i].end(); p!=pend; ++p) - { - PartitionType old_type=PartitionType(partition_type_indicator_->point_indicator_[*p]); - - if ( old_type != BorderEntity && new_type == FrontEntity) - { - partition_type_indicator_->point_indicator_[*p] = FrontEntity; - } else if(new_type == OverlapEntity && old_type != FrontEntity && old_type != BorderEntity) { - partition_type_indicator_->point_indicator_[*p] = OverlapEntity; + assert(new_order>=0); + assert(old_order>=0); + if(old_order== 3){//font + // something which is in contact with overlap or front can not be in contact with interior face + //assert(new_order !=0 ); + if(new_order == 0){ + std::cout << "Front and Interior face at same node" << std::endl; + } + if(new_order == 2){ + std::cout << "Front and Border face at same node" << std::endl; } + } + if(old_order == 2){// border + //assert(new_order != 3);// a border should not be should not be in contact with front. + if(new_order ==3){ + std::cout << "Border and Front face at same node" << std::endl; + } + } + + + if(new_order> old_order){ + partition_type_indicator_->point_indicator_[*p]= type_order[new_order]; + } + + // if(old_type==InteriorEntity) + // { + // if(new_type!=OverlapEntity) + // partition_type_indicator_->point_indicator_[*p]=new_type; + // } + // if(old_type==OverlapEntity) + // partition_type_indicator_->point_indicator_[*p]=new_type; + // if(old_type==FrontEntity && new_type==BorderEntity) + // partition_type_indicator_->point_indicator_[*p]=new_type; } } #endif @@ -1697,6 +1738,7 @@ void CpGridData::computeCommunicationInterfaces([[maybe_unused]] int noExistingP face_interfaces_); std::vector >().swap(face_attributes); */ + /* std::vector > point_attributes(noExistingPoints); AttributeDataHandle > > point_handle(ccobj_.rank(), *partition_type_indicator_, @@ -1708,6 +1750,51 @@ void CpGridData::computeCommunicationInterfaces([[maybe_unused]] int noExistingP } createInterfaces(point_attributes, partition_type_indicator_->point_indicator_.begin(), point_interfaces_); + */ + //point inter_face all + size_t nc = cell_to_point_.size(); + std::vector points; + for (size_t cell = 0; cell < nc; ++cell){ + if (false) { + points.clear(); + for (int i = 0; i < 8; ++i) { + points.push_back(cell_to_point_[cell][i]); + } + cell_to_allpoint_.appendRow(points.begin(), points.end()); + } else { + const auto& faces = cell_to_face_[cpgrid::EntityRep<0>(cell, true)]; + int nf = faces.size(); + // new code + points.clear(); + for (int f = 0; f < nf; ++f) { + auto face = faces[f].index(); + const auto& fnodes = face_to_point_[face]; + //auto faceSize = fnodes.size(); + // std::set nodes; + for(const auto& fv: fnodes){ + //for (int v = 0; v < faceSize; ++v) { + // int fv = fnodes[v]; + points.push_back(fv); + } + } + std::sort(points.begin(), points.end()); + auto unique_end = std::unique(points.begin(), points.end()); + points.erase(unique_end, points.end()); + cell_to_allpoint_.appendRow(points.begin(), points.end()); + } + } + + std::vector > allpoint_attributes(noExistingPoints); + AttributeDataHandle > + allpoint_handle(ccobj_.rank(), *partition_type_indicator_, + allpoint_attributes, cell_to_allpoint_, *this); + if( static_cast(std::get(cell_interfaces_)) + .interfaces().size() ) + { + comm.forward(allpoint_handle); + } + createInterfaces(allpoint_attributes, partition_type_indicator_->point_indicator_.begin(), + point_interfaces_); #endif } @@ -1781,6 +1868,43 @@ bool CpGridData::hasNNCs(const std::vector& cellIndices) const return hasNNC; } +bool CpGridData::compatibleSubdivisions(const std::vector>& cells_per_dim_vec, + const std::vector>& startIJK_vec, + const std::vector>& endIJK_vec) const +{ + bool compatibleSubdivisions = true; + if (startIJK_vec.size() > 1) { + bool notAllowedYet = false; + for (std::size_t level = 0; level < startIJK_vec.size(); ++level) { + for (std::size_t otherLevel = level+1; otherLevel < startIJK_vec.size(); ++otherLevel) { + const auto& sharedFaceTag = Opm::Lgr::sharedFaceTag({startIJK_vec[level], startIJK_vec[otherLevel]}, + {endIJK_vec[level],endIJK_vec[otherLevel]}, + this->logicalCartesianSize()); + if(sharedFaceTag == -1){ + break; // Go to the next "other patch" + } + if (sharedFaceTag == 0 ) { + notAllowedYet = notAllowedYet || + ((cells_per_dim_vec[level][1] != cells_per_dim_vec[otherLevel][1]) || (cells_per_dim_vec[level][2] != cells_per_dim_vec[otherLevel][2])); + } + if (sharedFaceTag == 1) { + notAllowedYet = notAllowedYet || + ((cells_per_dim_vec[level][0] != cells_per_dim_vec[otherLevel][0]) || (cells_per_dim_vec[level][2] != cells_per_dim_vec[otherLevel][2])); + } + if (sharedFaceTag == 2) { + notAllowedYet = notAllowedYet || + ((cells_per_dim_vec[level][0] != cells_per_dim_vec[otherLevel][0]) || (cells_per_dim_vec[level][1] != cells_per_dim_vec[otherLevel][1])); + } + if (notAllowedYet){ + compatibleSubdivisions = false; + break; + } + } // end-otherLevel-for-loop + } // end-level-for-loop + }// end-if-patchesShareFace + return compatibleSubdivisions; +} + std::tuple< const std::shared_ptr, const std::vector>, // parent_to_refined_corners(~boundary_old_to_new_corners) const std::vector>> > // parent_to_children_faces (~boundary_old_to_new_faces) @@ -1925,7 +2049,7 @@ bool CpGridData::preAdapt() if (local_empty) mark_.resize(size(0)); } - + // Detect the maximum mark across processes, and rewrite // the local entry in mark_, i.e., // mark_[ element.index() ] = max{ local marks in processes where this element belongs to}. diff --git a/opm/grid/cpgrid/CpGridData.hpp b/opm/grid/cpgrid/CpGridData.hpp index 64bdbc635..3f2d29f87 100644 --- a/opm/grid/cpgrid/CpGridData.hpp +++ b/opm/grid/cpgrid/CpGridData.hpp @@ -285,8 +285,6 @@ class CpGridData /// /// \param[in,out] nnc Non-neighboring connections. /// - /// \param[in] remove_ij_boundary True to remove outer cell layer - /// /// \param[in] turn_normals Whether or not to turn all normals. /// This is intended for handling inputs with wrong orientations. /// @@ -336,7 +334,7 @@ class CpGridData { return cell_to_point_; } - + const auto& cellToPoint(int cellIdx) const { return cell_to_point_[cellIdx]; @@ -383,7 +381,7 @@ class CpGridData OPM_THROW(std::logic_error, "Vertex has no history record.\n"); } } - + /// Return global_cell_ of any level grid, or the leaf grid view (in presence of refinement). /// global_cell_ has size number of cells present on a process and maps to the underlying Cartesian Grid. /// @@ -437,6 +435,25 @@ class CpGridData void postAdapt(); private: + /// @brief Check compatibility of number of subdivisions of neighboring LGRs. + /// + /// Check shared faces on boundaries of LGRs. Not optimal since the code below does not take into account + /// active/inactive cells, instead, relies on "ijk-computations". + /// + /// @param [in] cells_per_dim_vec Vector of expected subdivisions per cell, per direction, in each LGR. + /// @param [in] startIJK_vec Vector of Cartesian triplet indices where each patch starts. + /// @param [in] endIJK_vec Vector of Cartesian triplet indices where each patch ends. + /// Last cell part of the lgr will be {endIJK_vec[patch][0]-1, ..., endIJK_vec[patch][2]-1}. + /// @return True if all block of cells either do not share faces on their boundaries, or they may share faces with compatible + /// subdivisions. Example: block1 and block2 share an I_FACE, then number of subdivisions NY NZ should coincide, i.e. + /// if block1, block2 cells_per_dim values are {NX1, NY1, NZ1}, {NX2, NY2, NZ2}, respectively, then NY1 == NY2 and + /// NZ1 == NZ2. + /// False if at least two blocks share a face and their subdivions are not compatible. In the example above, + /// if NY1 != NY2 or NZ1 != NZ2. + bool compatibleSubdivisions(const std::vector>& cells_per_dim_vec, + const std::vector>& startIJK_vec, + const std::vector>& endIJK_vec) const; + std::array,8> getReferenceRefinedCorners(int idx_in_parent_cell, const std::array& cells_per_dim) const; public: @@ -490,7 +507,7 @@ class CpGridData } return level_to_leaf_cells_[level_cell_idx]; } - + /// @brief Refine a single cell and return a shared pointer of CpGridData type. /// /// refineSingleCell() takes a cell and refines it in a chosen amount of cells (per direction); creating the @@ -773,6 +790,9 @@ class CpGridData cpgrid::OrientedEntityTable<1, 0> face_to_cell_; /** @brief Container for the lookup of the points for each face. */ Opm::SparseTable face_to_point_; + /** @brief Container for the lookup of the ALL points i.e. not only the cannonical + * for each cell used for communication. */ + Opm::SparseTable cell_to_allpoint_; /** @brief Vector that contains an arrays of the points of each cell*/ std::vector< std::array > cell_to_point_; /** @brief The size of the underlying logical cartesian grid. @@ -806,7 +826,9 @@ class CpGridData /** @brief The global id set (used also as local id set). */ std::shared_ptr global_id_set_; /** @brief The indicator of the partition type of the entities */ + public: std::shared_ptr partition_type_indicator_; + private: /** Mark elements to be refined **/ std::vector mark_; /** Level of the current CpGridData (0 when it's "GLOBAL", 1,2,.. for LGRs). */ @@ -834,6 +856,7 @@ class CpGridData /** To keep track of refinement processes */ int refinement_max_level_{0}; + /// \brief Object for collective communication operations. Communication ccobj_; diff --git a/opm/grid/cpgrid/DataHandleWrappers.hpp b/opm/grid/cpgrid/DataHandleWrappers.hpp index e89e08744..0f725e1e4 100644 --- a/opm/grid/cpgrid/DataHandleWrappers.hpp +++ b/opm/grid/cpgrid/DataHandleWrappers.hpp @@ -148,7 +148,7 @@ template struct PointViaCellHandleWrapper : public PointViaCellWarner { using DataType = typename Handle::DataType; - using C2PTable = std::vector< std::array >; + using C2PTable = Opm::SparseTable;//std::vector< std::array >; /// \brief Constructs the data handle /// diff --git a/opm/grid/cpgrid/PartitionTypeIndicator.cpp b/opm/grid/cpgrid/PartitionTypeIndicator.cpp index 0879f9cc6..2eab07054 100644 --- a/opm/grid/cpgrid/PartitionTypeIndicator.cpp +++ b/opm/grid/cpgrid/PartitionTypeIndicator.cpp @@ -60,20 +60,46 @@ PartitionType PartitionTypeIndicator::getFacePartitionType(int i) const int cell_index = cells_of_face[0].index(); Entity<0> cell0(*grid_data_, cell_index, true); PartitionType cell_part = getPartitionType(cell0); - return cell_part; - } else { - // One of the cells might have index std::numeric_limits::max(). - // That means that the neighbor is not on this process but on - // another one. In this case the intersection is Front - auto idx0 = cells_of_face[0].index(); - auto idx1 = cells_of_face[1].index(); - if (idx0 == std::numeric_limits::max() || idx1 == std::numeric_limits::max()) { - // One neighbor is on another process => FrontEntity - return FrontEntity; + if(cell_part!=OverlapEntity){ + assert(cell_part == InteriorEntity); + return cell_part; + } + else + { + assert(cell_part == OverlapEntity); + // If the cell is in the overlap and the face is on the boundary, + // then the partition type has to Front! Here we check whether + // we are at the boundary. + OrientedEntityTable<0,1>::row_type cell_to_face=grid_data_->cell_to_face_[cell0]; + Entity<0>::LeafIntersectionIterator intersection=cell0.ilevelbegin(); + for(int subindex=0; subindex cell0(*grid_data_, cells_of_face[0].index(), true); + Entity<0> cell1(*grid_data_, cells_of_face[1].index(), true); + if(cells_of_face[0].index()==std::numeric_limits::max()) + { + assert(cells_of_face[1].index()!=std::numeric_limits::max()); + // At the boder of the processor's but not the global domain + return getProcessorBoundaryPartitionType(getPartitionType(cell1)); + } + if(cells_of_face[1].index()==std::numeric_limits::max()) + { + assert(cells_of_face[0].index()!=std::numeric_limits::max()); + // At the boder of the processor's but not the global domain + return getProcessorBoundaryPartitionType(getPartitionType(cell0)); } - Entity<0> cell0(*grid_data_, idx0, true); - Entity<0> cell1(*grid_data_, idx1, true); if(getPartitionType(cell0) == getPartitionType(cell1)) return getPartitionType(cell0); else diff --git a/opm/grid/cpgrid/PartitionTypeIndicator.hpp b/opm/grid/cpgrid/PartitionTypeIndicator.hpp index 2a0b0b262..d8327ab48 100644 --- a/opm/grid/cpgrid/PartitionTypeIndicator.hpp +++ b/opm/grid/cpgrid/PartitionTypeIndicator.hpp @@ -69,7 +69,8 @@ class PartitionTypeIndicator /// \return The partition type of the point. PartitionType getPartitionType(const EntityRep<3>& point_entity) const; -private: +//private: +public: /// Get the partition type of a face by its index /// \param i The index of the face. /// \return The partition type of the face associated with this index. diff --git a/opm/grid/cpgrid/processEclipseFormat.cpp b/opm/grid/cpgrid/processEclipseFormat.cpp index 6db619602..224406d84 100644 --- a/opm/grid/cpgrid/processEclipseFormat.cpp +++ b/opm/grid/cpgrid/processEclipseFormat.cpp @@ -205,8 +205,10 @@ namespace cpgrid return transMult.getMultiplier(cartindex, ::Opm::FaceDir::ZPlus) * transMult.getMultiplier(cartindex, ::Opm::FaceDir::ZMinus); }; + bool mergeMinPVCells = edge_conformal;// try to make geometrically connected grids + minpv_result = mp.process(thickness, z_tolerance, ecl_grid.getPinchMaxEmptyGap(), - poreVolume, ecl_grid.getMinpvVector(), actnumData, false, + poreVolume, ecl_grid.getMinpvVector(), actnumData, edge_conformal, zcornData.data(), nogap, pinchOptionALL, permZ, multZ, tolerance_unique_points); if (!minpv_result.nnc.empty()) { @@ -224,6 +226,9 @@ namespace cpgrid // Add PINCH NNCs. std::vector pinchedNNCs; + if(edge_conformal){ + assert(minpv_result.nnc.size() == 0); + } for (const auto& [cell1, cell2] : minpv_result.nnc) { nnc_cells[PinchNNC].insert({cell1, cell2}); @@ -416,12 +421,19 @@ namespace cpgrid } else { // Make the grid. + auto pinchActive_copy = pinchActive; + if(edge_conformal){ + // In edge conformal grids we need to set the pinchActive flag to true + pinchActive_copy = true; + assert(nnc_cells[PinchNNC].empty()); + assert(nnc_cells[ExplicitNNC].empty()); + } this->processEclipseFormat(g, ecl_state, nnc_cells, false, turn_normals, - pinchActive, + pinchActive_copy, tolerance_unique_points, edge_conformal); }