diff --git a/opm/grid/CpGrid.hpp b/opm/grid/CpGrid.hpp index 240c915d1..711b072ea 100644 --- a/opm/grid/CpGrid.hpp +++ b/opm/grid/CpGrid.hpp @@ -1744,6 +1744,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/common/GridPartitioning.cpp b/opm/grid/common/GridPartitioning.cpp index dd7998925..015cafff8 100644 --- a/opm/grid/common/GridPartitioning.cpp +++ b/opm/grid/common/GridPartitioning.cpp @@ -420,6 +420,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,26 +448,20 @@ 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 use 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); } } } @@ -571,16 +566,17 @@ 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 && !addCornerCells) { 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); + // correct branch to use in genneral case + 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 958e9cd05..f2fb6bf24 100644 --- a/opm/grid/common/GridPartitioning.hpp +++ b/opm/grid/common/GridPartitioning.hpp @@ -115,7 +115,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 2eadb460b..50d9cb10c 100644 --- a/opm/grid/cpgrid/CpGrid.cpp +++ b/opm/grid/cpgrid/CpGrid.cpp @@ -394,6 +394,7 @@ CpGrid::scatterGrid(EdgeWeightMethod method, comm().barrier(); // first create the overlap + auto vertex_cell_map = this->vertexCell(); auto noImportedOwner = addOverlapLayer(*this, computedCellPart, exportList, @@ -1168,6 +1169,32 @@ int CpGrid::cellFace(int cell, int local_index, int level) const : current_view_data_->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_view_data_->cell_to_face_[cpgrid::EntityRep<0>(cell, true)]; diff --git a/opm/grid/cpgrid/CpGridData.cpp b/opm/grid/cpgrid/CpGridData.cpp index 0c5b29b3c..82bfdf505 100644 --- a/opm/grid/cpgrid/CpGridData.cpp +++ b/opm/grid/cpgrid/CpGridData.cpp @@ -205,6 +205,12 @@ PartitionType getPartitionType(const PartitionTypeIndicator& p, int i, return p.getPartitionType(Entity<3>(grid, i, true)); } +int getIndex(std::vector::const_iterator i) +{ + return *i; +} + + int getIndex(const int* i) { return *i; @@ -924,7 +930,7 @@ struct AttributeDataHandle bool fixedSize() { - return true; + return false; } std::size_t size(std::size_t i) { @@ -933,8 +939,7 @@ 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(); + for(auto f=c2e_[i].begin(), fend=c2e_[i].end(); f!=fend; ++f) { char t=getPartitionType(indicator_, *f, grid_); @@ -945,8 +950,7 @@ 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(); + for(auto f=c2e_[i].begin(), fend=c2e_[i].end(); f!=fend; ++f, --s) { std::pair rank_attr; @@ -1607,6 +1611,8 @@ 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); for(int i=0; igetFacePartitionType(i); PartitionType old_type=PartitionType(partition_type_indicator_->point_indicator_[*p]); - if(old_type==InteriorEntity) - { - if(new_type!=OverlapEntity) - partition_type_indicator_->point_indicator_[*p]=new_type; + 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; + } + } + 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 + 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 + 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==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 @@ -1669,6 +1697,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_, @@ -1680,6 +1709,43 @@ 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) { + 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 } diff --git a/opm/grid/cpgrid/CpGridData.hpp b/opm/grid/cpgrid/CpGridData.hpp index 24d77299c..c2f642622 100644 --- a/opm/grid/cpgrid/CpGridData.hpp +++ b/opm/grid/cpgrid/CpGridData.hpp @@ -759,6 +759,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. @@ -792,7 +795,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). */ diff --git a/opm/grid/cpgrid/DataHandleWrappers.hpp b/opm/grid/cpgrid/DataHandleWrappers.hpp index d75e91f19..5dcf857c3 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 65f788637..80870b187 100644 --- a/opm/grid/cpgrid/PartitionTypeIndicator.cpp +++ b/opm/grid/cpgrid/PartitionTypeIndicator.cpp @@ -58,23 +58,27 @@ 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); - if(cell_part!=OverlapEntity) + 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