Give correct labeling of nodes in parallel

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<int> >, 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;

opm/grid/common/GridPartitioning.cpp

@@ -420,6 +420,7 @@ void addOverlapLayer(const CpGrid& grid,
                      const std::vector<int>& cell_part,
                      std::vector<std::tuple<int,int,char>>& exportList,
                      bool addCornerCells,
+                     const std::array<std::vector< std::set<int> >,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<int> 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<std::vector<std::set<int>>,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

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
 {

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<std::vector<std::set<int>>, 2>
+CpGrid::vertexCell() const
+{
+    //int level = -1;
+    int nc = numCells();//(level)
+    int nv = this->numVertices();
+    std::vector<std::set<int>> vertex_cell(nv);
+    std::vector<std::set<int>> 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<std::vector<std::set<int>>,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)];

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<int>::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<class B>
     void gather(B& buffer, std::size_t i)
     {
-        typedef typename GetRowType<T>::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<class B>
     void scatter(B& buffer, std::size_t i, std::size_t s)
     {
-        typedef typename GetRowType<T>::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<int,char> 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<Dune::PartitionType,4> type_order({InteriorEntity,OverlapEntity,BorderEntity,FrontEntity});
     partition_type_indicator_->point_indicator_.resize(geometry_.geomVector<3>().size(),
                                                        InteriorEntity);
     for(int i=0; i<face_to_point_.size(); ++i)
@@ -1616,15 +1622,37 @@ void CpGridData::computePointPartitionType()
         {
             PartitionType new_type=partition_type_indicator_->getFacePartitionType(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<std::map<int,char> >().swap(face_attributes);
     */
+    /*
     std::vector<std::map<int,char> > point_attributes(noExistingPoints);
     AttributeDataHandle<std::vector<std::array<int,8> > >
         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<int> 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<int> 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<std::map<int,char> > allpoint_attributes(noExistingPoints);
+    AttributeDataHandle<Opm::SparseTable<int> >
+        allpoint_handle(ccobj_.rank(), *partition_type_indicator_,
+                     allpoint_attributes, cell_to_allpoint_, *this);
+    if( static_cast<const Dune::Interface&>(std::get<All_All_Interface>(cell_interfaces_))
+        .interfaces().size() )
+    {
+        comm.forward(allpoint_handle);
+    }
+    createInterfaces(allpoint_attributes, partition_type_indicator_->point_indicator_.begin(),
+                     point_interfaces_);
 #endif
 }
 

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<int>             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<int>             cell_to_allpoint_;
     /** @brief Vector that contains an arrays of the points of each cell*/
     std::vector< std::array<int,8> >       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<LevelGlobalIdSet> global_id_set_;
     /** @brief The indicator of the partition type of the entities */
+    public:
     std::shared_ptr<PartitionTypeIndicator> partition_type_indicator_;
+    private:
     /** Mark elements to be refined **/
     std::vector<int> mark_;
     /** Level of the current CpGridData (0 when it's "GLOBAL", 1,2,.. for LGRs). */

opm/grid/cpgrid/DataHandleWrappers.hpp

@@ -148,7 +148,7 @@ template<class Handle>
 struct PointViaCellHandleWrapper : public PointViaCellWarner
 {
     using DataType = typename Handle::DataType;
-    using C2PTable = std::vector< std::array<int,8> >;
+    using C2PTable = Opm::SparseTable<int>;//std::vector< std::array<int,8> >;
 
     /// \brief Constructs the data handle
     ///

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<cell_to_face.size(); ++subindex, ++intersection)
+                for(int subindex=0; subindex<cell_to_face.size(); ++subindex, ++intersection) // Correct ?
                     if(cell_to_face[subindex].index()==i)
                         break;
                 assert(intersection!=cell0.ilevelend());
-                if(intersection.boundary())
-                    return FrontEntity;
-                else
+                if(intersection.boundary()){
+                    return OverlapEntity;
+                }else{
                     return cell_part;
+                }
             }
         }
         else