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

CMakeLists.txt

@@ -82,6 +82,7 @@ macro(opm-grid_tests_hook)
       addLgrsOnDistributedGrid_test
       autoRefine_test
       distribution_test
+      partition_of_unity_test
       level_and_grid_cartesianIndexMappers_test
       logicalCartesianSize_and_refinement_test
       test_communication_utils

CMakeLists_files.cmake

@@ -33,7 +33,6 @@ list(APPEND MAIN_SOURCE_FILES
   opm/grid/cpgrid/CpGridData.cpp
   opm/grid/cpgrid/CpGrid.cpp
   opm/grid/cpgrid/CpGridUtilities.cpp
-  opm/grid/cpgrid/DataHandleWrappers.cpp
   opm/grid/cpgrid/GridHelpers.cpp
   opm/grid/cpgrid/Iterators.cpp
   opm/grid/cpgrid/Indexsets.cpp
@@ -100,6 +99,7 @@ list(APPEND TEST_SOURCE_FILES
   tests/test_sparsetable.cpp
   tests/test_subgridpart.cpp
   tests/cpgrid/distribution_test.cpp
+  tests/cpgrid/partition_of_unity_test.cpp
   tests/cpgrid/entityrep_test.cpp
   tests/cpgrid/entity_test.cpp
   tests/cpgrid/facetag_test.cpp

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,47 @@ 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.
+    // Local: consumed by the handle below and not needed afterwards, so it
+    // costs nothing once the interfaces are built.
+    Opm::SparseTable<int> cell_to_allpoint;
+    const std::size_t nc = cell_to_point_.size();
+    std::vector<int> points;
+    // seen[p] == cell+1 marks p as already taken for this cell, so the
+    // duplicate check is O(1) per point rather than a linear scan.
+    std::vector<std::size_t> seen(noExistingPoints, 0);
+    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 (seen[fv] != cell + 1) {
+                    seen[fv] = cell + 1;
+                    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() )
     {

opm/grid/cpgrid/DataHandleWrappers.hpp

@@ -129,98 +129,6 @@ struct FaceViaCellHandleWrapper
     const C2FTable& c2fGather_, c2f_;
 };
 
-struct PointViaCellWarner
-{
-    static bool printWarn;
-    void warn();
-};
-
-/// \brief A data handle to send data attached to points via cell communication
-///
-/// With it we can use the cell communication to also communicate data attached
-/// to faces.
-/// \warning As we send all points of cell most of the points will be send six times
-/// which will temporarily waste some space and bandwidth.
-///
-/// \tparam Handle The type of the data handle to wrap. It must gather and scatter
-///         only for codim-3 entities.
-template<class Handle>
-struct PointViaCellHandleWrapper : public PointViaCellWarner
-{
-    using DataType = typename Handle::DataType;
-    using C2PTable = std::vector< std::array<int,8> >;
-
-    /// \brief Constructs the data handle
-    ///
-    /// \param handle Handle object wrapped and used for actual gather/scatter
-    /// \param c2pGather Table to determine points when gathering
-    /// \param c2p Table to determine points when scattering
-    PointViaCellHandleWrapper(Handle& handle,
-                             const C2PTable& c2pGather,
-                             const C2PTable& c2p)
-        : handle_(handle), c2pGather_(c2pGather), c2p_(c2p)
-    {}
-    bool fixedSize(int i, int j)
-    {
-        if( ! handle_.fixedSize(i, j))
-        {
-            this->warn();
-        }
-        return handle_.fixedSize(i, j);
-    }
-    template<class T>
-    typename std::enable_if<T::codimension != 0, std::size_t>::type
-    size(const T&)
-    {
-        OPM_THROW(std::logic_error, "This should never throw! We only know sizes for cells");
-        return 1;
-    }
-    std::size_t size(const EntityRep<0>& t)
-    {
-        const auto& points = c2pGather_[t.index()];
-        return std::accumulate(points.begin(), points.end(), std::size_t{0},
-                               [this](const auto acc, const auto& point)
-                               { return acc + handle_.size(EntityRep<3>(point, true)); });
-    }
-    bool contains(std::size_t dim, std::size_t codim)
-    {
-        return dim==3 && codim == 0;
-    }
-    template<class B, class T>
-    typename std::enable_if<T::codimension != 0, void>::type
-    gather(B&, const T&)
-    {
-        OPM_THROW(std::logic_error, "This should never throw! We can only gather cell values and indicate that");
-    }
-    template<class B>
-    void gather(B& buffer, const EntityRep<0>& t)
-    {
-        const auto& points = c2pGather_[t.index()];
-        for (const auto& point : points)
-        {
-            handle_.gather(buffer, EntityRep<3>(point, true));
-        }
-    }
-    template<class B, class T>
-    typename std::enable_if<T::codimension != 0, void>::type
-    scatter(B&, const T&, std::size_t)
-    {
-        OPM_THROW(std::logic_error, "This should never throw! We can only gather cell values and indicate that");
-    }
-    template<class B>
-    void scatter(B& buffer, const EntityRep<0>& t, std::size_t s)
-    {
-        const auto& points = c2p_[t.index()];
-        for (const auto& point : points)
-        {
-            handle_.scatter(buffer, EntityRep<3>(point, true), s/8);
-        }
-    }
-private:
-    Handle& handle_;
-    const C2PTable& c2pGather_, c2p_;
-};
-
 } // end namespace Dune::cpgrid
 
 #endif