Store a variant in the geometry

To avoid copying a shared pointer, which is rather costly, the geometry now stores a variant holding either the ownership or a non-owning view to the entity variable. In the case of entity variables owned by the grid, a non-owning storage is fine as the geometry should not outlive the grid. On grid refinement the grid points are only owned by the geometry, so a move into the geometry is preferred.

Author: SoilRos

opm/grid/cpgrid/CpGrid.cpp

@@ -3845,7 +3845,7 @@ void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geom
 
         // Create a pointer to the first element of "adapted_cell_to_point" (required as the fourth argement to construct a Geometry<3,3> type object).
         int* indices_storage_ptr = adapted_cell_to_point[cell].data();
-        adapted_cells[cell] = cpgrid::Geometry<3,3>(cellGeom.center(), cellGeom.volume(), allCorners, indices_storage_ptr);
+        adapted_cells[cell] = cpgrid::Geometry<3,3>(cellGeom.center(), cellGeom.volume(), allCorners.get(), indices_storage_ptr);
     } // adapted_cells
 
     // Adapted/Leaf-grid-view face to cell.
@@ -3966,7 +3966,7 @@ void CpGrid::populateRefinedCells(std::vector<Dune::cpgrid::EntityVariableBase<c
 
             // Create a pointer to the first element of "refined_cell_to_point" (required as the fourth argement to construct a Geometry<3,3> type object).
             int* indices_storage_ptr = refined_cell_to_point_vec[shiftedLevel][cell].data();
-            refined_cells_vec[shiftedLevel][cell] = cpgrid::Geometry<3,3>(elemLgrGeom.center(), elemLgrGeom.volume(), allLevelCorners, indices_storage_ptr);
+            refined_cells_vec[shiftedLevel][cell] = cpgrid::Geometry<3,3>(elemLgrGeom.center(), elemLgrGeom.volume(), allLevelCorners.get(), indices_storage_ptr);
         } // refined_cells
         // Refined face to cell.
         refined_cell_to_face_vec[shiftedLevel].makeInverseRelation(refined_face_to_cell_vec[shiftedLevel]);

opm/grid/cpgrid/CpGridData.cpp

@@ -535,7 +535,7 @@ struct CellGeometryHandle
             buffer.read(pos[i]);
 
         buffer.read(vol);
-        scatterCont_[t] = Geom(pos, vol, pointGeom_, cell2Points_[t.index()].data());
+        scatterCont_[t] = Geom(pos, vol, pointGeom_.get(), cell2Points_[t.index()].data());
         double isAquifer;
         buffer.read(isAquifer);
         if (isAquifer == 1.0)
@@ -1577,8 +1577,8 @@ void CpGridData::distributeGlobalGrid(CpGrid& grid,
 
     // Compute partition type for points
     computePointPartitionType();
-   
-    computeCommunicationInterfaces(noExistingPoints);   
+
+    computeCommunicationInterfaces(noExistingPoints);
 #else // #if HAVE_MPI
     static_cast<void>(grid);
     static_cast<void>(view_data);
@@ -2065,9 +2065,9 @@ int CpGridData::sharedFaceTag(const std::vector<std::array<int,3>>& startIJK_2Pa
     assert(endIJK_2Patches.size() == 2);
 
     int faceTag = -1; // 0 represents I_FACE, 1 J_FACE, and 2 K_FACE. Use -1 for no sharing face case.
-     
+
     if (patchesShareFace(startIJK_2Patches, endIJK_2Patches)) {
-        
+
         const auto& detectSharing = [](const std::vector<int>& faceIdxs, const std::vector<int>& otherFaceIdxs){
             bool faceIsShared = false;
             for (const auto& face : faceIdxs) {
@@ -2080,7 +2080,7 @@ int CpGridData::sharedFaceTag(const std::vector<std::array<int,3>>& startIJK_2Pa
             }
             return faceIsShared; // should be false here
         };
-     
+
         const auto& [iFalse, iTrue, jFalse, jTrue, kFalse, kTrue] = this->getBoundaryPatchFaces(startIJK_2Patches[0], endIJK_2Patches[0]);
         const auto& [iFalseOther, iTrueOther, jFalseOther, jTrueOther, kFalseOther, kTrueOther] =
             this->getBoundaryPatchFaces(startIJK_2Patches[1], endIJK_2Patches[1]);
@@ -2306,7 +2306,7 @@ Geometry<3,3> CpGridData::cellifyPatch(const std::array<int,3>& startIJK, const
         const int* cellifiedPatch_indices_storage_ptr = &allcorners_cellifiedPatch[0];
         // Construct (and return) the Geometry<3,3> of the 'cellified patch'.
         return Geometry<3,3>(cellifiedPatch_center, cellifiedPatch_volume,
-                             cellifiedPatch_geometry.geomVector(std::integral_constant<int,3>()),
+                             cellifiedPatch_geometry.geomVector(std::integral_constant<int,3>()).get(),
                              cellifiedPatch_indices_storage_ptr);
     }
 }
@@ -2699,7 +2699,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/Entity.hpp

@@ -539,10 +539,9 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() cons
         const auto& auxArr = pgrid_ -> getReferenceRefinedCorners(idx_in_parent_cell, cells_per_dim);
         FieldVector<double, 3> corners_in_father_reference_elem_temp[8] =
             { auxArr[0], auxArr[1], auxArr[2], auxArr[3], auxArr[4], auxArr[5], auxArr[6], auxArr[7]};
-        auto in_father_reference_elem_corners = std::make_shared<EntityVariable<cpgrid::Geometry<0, 3>, 3>>();
-        EntityVariableBase<cpgrid::Geometry<0, 3>>& mutable_in_father_reference_elem_corners = *in_father_reference_elem_corners;
+        EntityVariable<cpgrid::Geometry<0, 3>, 3> in_father_reference_elem_corners;
         // Assign the corners. Make use of the fact that pointers behave like iterators.
-        mutable_in_father_reference_elem_corners.assign(corners_in_father_reference_elem_temp,
+        in_father_reference_elem_corners.assign(corners_in_father_reference_elem_temp,
                                                         corners_in_father_reference_elem_temp + 8);
         // Compute the center of the 'local-entity'.
         Dune::FieldVector<double, 3> center_in_father_reference_elem = {0., 0.,0.};
@@ -556,7 +555,7 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() cons
         double volume_in_father_reference_elem = double(1)/(cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2]);
         // Construct (and return) the Geometry<3,3> of 'child-cell in the reference element of its father (unit cube)'.
         return Dune::cpgrid::Geometry<3,3>(center_in_father_reference_elem, volume_in_father_reference_elem,
-                                           in_father_reference_elem_corners, in_father_reference_elem_corner_indices.data());
+                                           std::move(in_father_reference_elem_corners), in_father_reference_elem_corner_indices.data());
     }
     else {
         OPM_THROW(std::logic_error, "Entity has no father.");

opm/grid/cpgrid/Geometry.hpp

@@ -38,11 +38,13 @@
 #define OPM_GEOMETRY_HEADER
 
 #include <cmath>
+#include <variant>
 
 // Warning suppression for Dune includes.
 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
 
 #include <dune/common/version.hh>
+#include <dune/common/overloadset.hh>
 #include <dune/geometry/referenceelements.hh>
 #include <dune/grid/common/geometry.hh>
 
@@ -405,6 +407,10 @@ namespace Dune
 
             /// @brief Construct from center, volume (1- and 0-moments) and
             ///        corners.
+            /// @warning This constructor does not own the corners or indices,
+            ///         thus, the pointers must remain valid for the lifetime of
+            ///         the Geometry object.
+            ///
             /// @param pos the centroid of the entity
             /// @param vol the volume(area) of the entity
             /// @param allcorners pointer of all corner positions in the grid
@@ -413,17 +419,34 @@ namespace Dune
             ///                       by (kji), i.e. i running fastest.
             Geometry(const GlobalCoordinate& pos,
                      ctype vol,
-                     std::shared_ptr<const EntityVariable<cpgrid::Geometry<0, 3>, 3>> allcorners_ptr,
+                     EntityVariable<cpgrid::Geometry<0, 3>, 3> const * allcorners_view,
                      const int* corner_indices)
                 : pos_(pos), vol_(vol),
-                  allcorners_(allcorners_ptr), cor_idx_(corner_indices)
+                  allcorners_(allcorners_view), cor_idx_(corner_indices)
             {
-                assert(allcorners_ && corner_indices);
+                assert(allcorners_view && corner_indices);
             }
 
+            /// @brief Construct from center, volume (1- and 0-moments) and
+            ///        corners.
+            ///
+            /// @param pos the centroid of the entity
+            /// @param vol the volume(area) of the entity
+            /// @param allcorners pointer of all corner positions in the grid
+            /// @param corner_indices array of 8 indices into allcorners. The
+            ///                       indices must be given in lexicographical order
+            ///                       by (kji), i.e. i running fastest.
+            Geometry(const GlobalCoordinate& pos,
+                     ctype vol,
+                     const EntityVariable<cpgrid::Geometry<0, 3>, 3>& allcorners_storage,
+                     const int* corner_indices)
+                : pos_(pos), vol_(vol),
+                  allcorners_(allcorners_storage), cor_idx_(corner_indices)
+            {}
+
             /// Default constructor, giving a non-valid geometry.
             Geometry()
-                : pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
+                : pos_(0.0), vol_(0.0), allcorners_(nullptr), cor_idx_(nullptr)
             {
             }
 
@@ -511,8 +534,16 @@ namespace Dune
             /// @brief Get the cor-th of 8 corners of the hexahedral base cell.
             GlobalCoordinate corner(int cor) const
             {
-                assert(allcorners_ && cor_idx_);
-                return (allcorners_->data())[cor_idx_[cor]].center();
+                return std::visit(
+                    Dune::overload(
+                        [](const EntityVariable<Geometry<0, 3>, 3>& corners) -> const EntityVariable<Geometry<0, 3>, 3>& {
+                            return corners;
+                        },
+                        [](EntityVariable<Geometry<0, 3>, 3> const * corners) -> const EntityVariable<Geometry<0, 3>, 3>& {
+                            return *corners;
+                        }
+                    ), allcorners_
+                ).data()[cor_idx_[cor]].center();
             }
 
             /// Cell volume.
@@ -1022,7 +1053,7 @@ namespace Dune
                             refined_cells[refined_cell_idx] =
                                 Geometry<3,cdim>(refined_cell_center,
                                                  refined_cell_volume,
-                                                 all_geom.geomVector(std::integral_constant<int,3>()),
+                                                 all_geom.geomVector(std::integral_constant<int,3>()).get(),
                                                  indices_storage_ptr);
                         } // end i-for-loop
                     }  // end j-for-loop
@@ -1043,8 +1074,10 @@ namespace Dune
         private:
             GlobalCoordinate pos_;
             double vol_;
-            std::shared_ptr<const EntityVariable<Geometry<0, 3>,3>> allcorners_; // For dimension 3 only
-            const int* cor_idx_; // For dimension 3 only
+
+            // store either a (view) pointer or the variable itself. For dimension 3 only
+            std::variant<EntityVariable<Geometry<0, 3>,3> const *, EntityVariable<Geometry<0, 3>,3>> allcorners_;
+            int const * cor_idx_; // For dimension 3 only
 
             /// @brief
             ///   Auxiliary function to get refined_face information: tag, index, face_to_point_, face_to_cell, face centroid,

opm/grid/cpgrid/processEclipseFormat.cpp

@@ -1223,7 +1223,7 @@ namespace cpgrid
                                                       double vol,
                                                       const std::array<int,8>& corner_indices)
             {
-                return cpgrid::Geometry<3, 3>(pos, vol, allcorners_, &corner_indices[0]);
+                return cpgrid::Geometry<3, 3>(pos, vol, allcorners_.get(), &corner_indices[0]);
             }
         };
 

tests/cpgrid/geometry_test.cpp

@@ -152,7 +152,7 @@ BOOST_AUTO_TEST_CASE(cellgeom)
     }
 
     int cor_idx[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
-    Geometry g(c, v, pg, cor_idx);
+    Geometry g(c, v, pg.get(), cor_idx);
 
     // Verification of properties.
     BOOST_CHECK(g.type().isCube());
@@ -206,7 +206,7 @@ BOOST_AUTO_TEST_CASE(cellgeom)
     {
         (*pg1).push_back(cpgrid::Geometry<0, 3>(crn));
     }
-    g = Geometry(c, v, pg1, cor_idx);
+    g = Geometry(c, v, pg1.get(), cor_idx);
 
     // Verification of properties.
     BOOST_CHECK(g.type().isCube());
@@ -525,7 +525,7 @@ BOOST_AUTO_TEST_CASE(refine_simple_cube)
     }
 
     int cor_idx[8] = {0, 1, 2, 3, 4, 5, 6, 7};
-    Geometry g(c, v, pg, cor_idx);
+    Geometry g(c, v, pg.get(), cor_idx);
 
     refine_and_check(g, {1, 1, 1}, true);
     refine_and_check(g, {2, 3, 4}, true);
@@ -567,7 +567,7 @@ BOOST_AUTO_TEST_CASE(refine_distorted_cube)
     }
 
     int cor_idx[8] = {0, 1, 2, 3, 4, 5, 6, 7};
-    Geometry g(center, v, pg, cor_idx);
+    Geometry g(center, v, pg.get(), cor_idx);
     refine_and_check(g, {1, 1, 1});
     refine_and_check(g, {2, 3, 4});