CpGrid: build point communication interfaces from all points of a cell

The codim-3 (point) communication interface and the point-via-cell data
handle were built from cell_to_point_, i.e. only the eight canonical corners
of each cell. On corner-point grids with hanging nodes a face can reference
nodes that are not canonical corners of the neighbouring cell; such nodes
never entered point_interfaces_, so they were invisible to codim-3
communication and their ownership failed to be a partition of unity across
ranks (observed as 'Owner is not a partition of unity' in vertex-based
mechanics assembly).

Introduce cell_to_allpoint_ (SparseTable: every node of every face of each
cell, sorted and deduplicated) and use it for the point interface and for
PointViaCellHandleWrapper. AttributeDataHandle::fixedSize() now returns
false since rows vary in length, and row iteration uses the SparseTable
iterators.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: hnil

opm/grid/cpgrid/CpGridData.cpp

@@ -206,6 +206,12 @@ PartitionType getPartitionType(const PartitionTypeIndicator& p, int i,
     return p.getPartitionType(Entity<3>(grid, i, true));
 }
 
+// Row entries of Opm::SparseTable<int> are visited through this iterator type.
+int getIndex(std::vector<int>::const_iterator i)
+{
+    return *i;
+}
+
 int getIndex(const int* i)
 {
     return *i;
@@ -934,7 +940,8 @@ struct AttributeDataHandle
 
     bool fixedSize()
     {
-        return true;
+        // SparseTable rows have varying length (all points of a cell).
+        return false;
     }
     std::size_t size(std::size_t i)
     {
@@ -943,8 +950,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_);
@@ -955,8 +961,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;
@@ -1697,10 +1702,40 @@ void CpGridData::computeCommunicationInterfaces([[maybe_unused]] int noExistingP
     face_interfaces_);
     std::vector<std::map<int,char> >().swap(face_attributes);
     */
+    // Build the point (codim 3) communication interface from ALL points of
+    // each cell -- every node of every face -- not only the eight canonical
+    // corners stored in cell_to_point_.  On corner-point grids with hanging
+    // nodes a face can reference nodes that are not canonical corners of the
+    // neighbouring cell; interfaces built from the canonical corners only
+    // miss those nodes, leaving them without communication (node ownership
+    // then fails to be a partition of unity).
+    // The gather/scatter of AttributeDataHandle pairs entries of the
+    // sender's and the receiver's row for the same cell BY POSITION, so the
+    // rows must list the points in an order that is identical on every rank
+    // sharing the cell.  Local point indices differ between ranks; the
+    // face/point traversal order is copied from the global grid during
+    // distribution and is rank-independent.  Deduplicate while preserving
+    // the first-occurrence traversal order -- do NOT sort by local index.
+    const std::size_t nc = cell_to_point_.size();
+    std::vector<int> points;
+    for (std::size_t cell = 0; cell < nc; ++cell) {
+        points.clear();
+        const auto& faces = cell_to_face_[cpgrid::EntityRep<0>(cell, true)];
+        const int nf = faces.size();
+        for (int f = 0; f < nf; ++f) {
+            for (const auto& fv : face_to_point_[faces[f].index()]) {
+                if (std::find(points.begin(), points.end(), fv) == points.end()) {
+                    points.push_back(fv);
+                }
+            }
+        }
+        cell_to_allpoint_.appendRow(points.begin(), points.end());
+    }
+
     std::vector<std::map<int,char> > point_attributes(noExistingPoints);
-    AttributeDataHandle<std::vector<std::array<int,8> > >
+    AttributeDataHandle<Opm::SparseTable<int> >
         point_handle(ccobj_.rank(), *partition_type_indicator_,
-                     point_attributes, cell_to_point_, *this);
+                     point_attributes, cell_to_allpoint_, *this);
     if( static_cast<const Dune::Interface&>(std::get<All_All_Interface>(cell_interfaces_))
         .interfaces().size() )
     {
@@ -1928,7 +1963,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}.

opm/grid/cpgrid/CpGridData.hpp

@@ -775,6 +775,11 @@ 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 All points of each cell, i.e. every point of every face of the
+     *   cell, not only the eight canonical corners.  On corner-point grids
+     *   with hanging nodes the extra points are essential for building
+     *   complete codim-3 communication interfaces. */
+    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.

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>;
 
     /// \brief Constructs the data handle
     ///